Apalab -- A Matlab toolbox for Approximate Polynomial Algebra
Copyright � 2007 by Zhonggang Zeng.  All rights reserved.

(Updated on June 25, 2007)

Apalab  is a Matlab toolbox for numerical and symbolic computation in Approximate Polynomial Algebra. Its Maple counterpart is named  ApaTools.  These packages are integral components of a research project in development by Zhonggang Zeng, Department of Mathematics, Northeastern Illinois University, Chicago, IL 60625, supported by National Science Foundation under grants DMS-0715127 and DMS-0412003.  Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation.

website: http://www.neiu.edu/~zzeng                   email: zzeng@neiu.edu

Preprint:   ApaTools:  A software toolbox for approximate polynomial algebra,  Zhonggang Zeng, 2007  

Disclaimer:   Apalab is an on-going research project, not yet a finished product.  Bugs and errors almost certainly exist.  The author is not responsible for any damage  or loss from using this software. Users should verify their results through other means.
 

Table of contents

• Software tools currently available in the box
• Download and installation
• Representation of polynomials in Apalab
• Modules with brief explanations and examples

Software tools currently available in the box:
          (click each module name for a brief explanation and examples)

Univariate polynomial:         uvGCD     ExtendedGCD     uvFactor
Multivariate polynomial:       mvGCD      SquarefreeFactor     mvFactorRefine 
Matrix computation:             ApproxiJordanForm      NulVector        ApproxiRank  
                                          ApproxiRankRowUpdate   ApproxiRankRowDowndate  ApproxiRankColumnUpdate   ApproxiRankColumnDowndate     
Utilities:            
      Polynomial manipulation:       Polynomial2CoefficientMatrix
                                                  mvPolynDegree            mvPolynTimesMonom      mvPolynAdd        mvPolynMultiply      mvPolynEvaluate
                                                  mvPolynDerivative    mvPolynClear     mvPolynPower
      Matrix computation:               HouseholderVector     HouseholderTransform                    HouseholderTransformRight      GivensMatrix
                                                IntegerInverse                MatrixWithGivenJordanForm      ScaledLeastSquares

Download and installation:

1. Download the package zip file according to your Matlab version:
   (Support for older versions of Matlab are possible upon request)
        -- For Matlab 2008:    ApalabForMatlab2008.zip
        -- For Matlab 2009:    ApalabForMatlab2009.zip
        -- For Matlab 2010:    ApalabForMatlab2010.zip
        -- For Matlab 2011:    ApalabForMatlab2011.zip
2. Unzip  the zip file in a directory,  say  "c:\polynomials",  so that a director  "c:\polynomials\ApalabForMatlab61" or  "c:\polynomials\ApalabForMatlab65"  etc will be created and contains all the necessary files.
3. Start Matlab.
4. Set path to the directory, e.g.    >> path(path,'c:/polynomials/ApalabForMatlab61')  
5. Now all the tool modules will be accessible. 
6. To access the help file of each module, say uvGCD,   try    >> help uvGCD    

Representation of polynomials in Apalab

Coefficient vectors of univariate polynomials:  A univariate polynomial is represented by its coefficient vector.  For example,  f(x) = x³ + 2x² + 3x + 4  is represented as
               f  =  [1, 2, 3, 4]   or  its transpose
Polynomial multiplication  h = f*g,  is carried out in Matlab by 
             >>  h = conv(f,g)
with  f, g, h  being coefficient vectors.

Coefficient matrix:  An  m-multivariate polynomial  f  of  n  terms is represented by a coefficient matrix  F  of size  (m+1) � n  according to the following rules:
     1.  Every column of  F  represents one term (monomial)
     2.  F(i,j)    for i<m+1   is the exponent of the i-th variable on j-th term
     3.  F(m+1,j)  is the coefficient of the j-th term
For example, let
          p(x, y, z)  =  8  +  3x³ y  +  5 x² z⁵   -  2 y³   + 6 y³ z²  - 7 x y z
Then it is represented as the following coefficient matrix in Matlab
 
>> P = [0 3 2 0 0 1; 0 1 0 3 3 1; 0 0 5 0 2 1; 8 3 5 -2 6 -7]

P =

      0    3    2    0    0     1
      0    1    0    3    3     1
      0    0    5    0    2     1
      8    3    5   -2    6    -7

Coefficient matrices can be constructed using the Symbolic Math Toolbox and function Polynomial2CoefficientMatrix, as shown in the following example

>> syms x y z  % declare symbolic variables for the polynomial
>> p = 8 + 3*x^3*y + 5*x^2*z^5 - 2*y^3 + 6*y^3*z^2 - 7*x*y*z;
>> pretty(p)   % showcase the polynomial

                       3        2  5      3      3  2
                8 + 3 x  y + 5 x  z  - 2 y  + 6 y  z  - 7 x y z

>> F = Polynomial2CoefficientMatrix(p,[x,y,z])

 F =

     0     3     2     0     0     1
     0     1     0     3     3     1
     0     0     5     0     2     1
     8     3     5    -2     6    -7

Maple  module  Polynomial2CoefficientMatrix  is part of package  ApaTools  for converting a polynomial into a coefficient matrix. Its simplified version is provided in a Maple classic worksheet  Polyn2CoeffMat.mws.  Right-click the link, save as a .mws file, then open it in Maple. Here is a screen shot of the worksheet.

Basic tools for manipulating multivariate polynomials in coefficient matrices:  One can add/subtract, multiply two multivariate polynomials, evaluate a polynomial, etc. using utility tools  mvPolynDegree,  mvPolynTimesMonom,   mvPolynAdd,  mvPolynMultiply,   mvPolynEvaluate,  mvPolynDerivative,  mvPolynClear,   mvPolynPower 
                                    

Total degree  and  tuple degree:   The total degree of a monomial is the sum of the exponents on all variables in the monomial.  The total degree of a polynomial is the maximum of total degrees of all terms in the polynomial.

For an  m-variable polynomial  f( x₁ , ..., x_m ),  its m-tuple degree is a vector  [d₁ , d₂ , ..., d_m ]  where  d_j  is the degree of  f  in  x_j.   We call this type of degree a tuple degree, or  an  m-tuple degree  if  m  is to be specified.

Utility module   mvPolynDegree  calculates either the total or tuple degree of a given polynomial in coefficient matrix.

Modules with brief explanations and examples

• NulVector
  The module for computing the smallest singular value of an upper-triangular matrix along with the corresponding left and right singular vectors. In many cases, it is desirable to determine if a matrix is approximately rank-deficient without computing a full-blown singular value decomposition. Module NulVector  is designed for this purpose. It also serves as a submodule for  ApproxiRank  that computes the approximate rank and approximate kernel.
  
   syntax: >> LHS = RHS
  
    Input parameters:
                R -- m�n upper triangular matrix (m ≥ n)
            theta -- (optional) the rank threshold
            scale -- (optional) ||R||the rank threshold
             m, n -- (optional) size of R
  
   Output parameters:
                s -- an approximation to the smallest singular value
                x -- the unit approximate null vector
                y -- the unit left approximate null vector
  
  Example:
  
  >> A  % a seemingly full-rank matrix
  
  A =
      1  10   0   0   0   0   0   0   0   0
      0   1  10   0   0   0   0   0   0   0
      0   0   1  10   0   0   0   0   0   0
      0   0   0   1  10   0   0   0   0   0
      0   0   0   0   1  10   0   0   0   0
      0   0   0   0   0   1  10   0   0   0
      0   0   0   0   0   0   1  10   0   0
      0   0   0   0   0   0   0   1  10   0
      0   0   0   0   0   0   0   0   1  10
      0   0   0   0   0   0   0   0   0   1
  
  >> [s,x,y] = NulVector(A,1e-10);
  >> s        % the smallest singular value is tiny:
  
  s =
     9.900000000000001e-010
  
  >> [x,y]  % the right and left singular vectors
  
  ans =
      -0.99498743710662     -0.00000000098504
       0.09949874371066      0.00000000994888
      -0.00994987437107     -0.00000009949864
       0.00099498743711      0.00000099498743
      -0.00009949874371     -0.00000994987437
       0.00000994987437      0.00009949874371
      -0.00000099498743     -0.00099498743711
       0.00000009949864      0.00994987437107
      -0.00000000994888     -0.09949874371066
       0.00000000098504      0.99498743710662
  
  >> [norm(A*x), norm(y'*A)]  % verify them as approxi-nulvectors
  
  ans =
    1.0e-009 *
      0.99000000000001 0.99000000000001
  
   
• ApproxiRank
  The module  ApproxiRank  computes the approximate rank rank_θ (A) and the approximate kernel  K_θ (A)  of a given matrix  A within a given rank threshold  θ.   The definitions of the approximate rank and approximate kernel are
              rank_θ (A)  = the smallest rank of all the matrix  B   such that  ||B-A|| < θ.
              K_θ (A)   =  the kernel of the matrix  B  that is the nearest to  A  among all matrices whose exact rank equals to rank_θ (A).
  Module  ApproxiRank  is efficient when the approximate nullity of  A  is low.
  
  Reference:   A rank-revealing method with updating, downdating and applications, T.Y. Li  and  Z. Zeng,  SIAM Journal on Matrix Analysis and Applications, 26(2005),  pp918-946.
  
  Computing the approximate rank and null space of a matrix
  
  Syntax:
      fast execution without complete QR
          >>[ark,s,X] = ApproxiRank(A) % using default tol
          >>[ark,s,X] = ApproxiRank(A,tol)
  complete output, required for updating/downdating
          >>[ark,s,X,R,Q,scale] = ApproxiRank(A,tol)
  
  <Input Parameters>
          A -- mxn matrix (m>=n)
        tol -- (optional) the rank decision threshold
  
  <Output parameters>
        ark -- the approximate rank
          s -- the approximate singular values of A below the threshold
          X -- matrix whose columns form an orthonormal bases for the approximate null space
       R, Q -- (optional) the QR decomposition used for up-/downdating
      scale -- (optional) scaling factor
  
   
• ApproxiRankRowUpdate
  Row updating module for ApproxiRank.  After computing the approximate rank and kernel of a matrix  A,  if a new row is inserted into  A to form a new matrix  B,   ApproxiRankRowUpdate  calculates the approximate rank/kernel of  B  by updating the approximate rank/kernel of  A. 
  
  Reference:   A rank-revealing method with updating, downdating and applications, T.Y. Li  and  Z. Zeng,  SIAM Journal on Matrix Analysis and Applications, 26(2005),  pp918-946.
  
  Syntax
      >>[ark,Y] = ApproxiRankRowUpdate(a,k,X,Q,R,tol)
      >>[ark,Y,R,Q] = ApproxiRankRowUpdate(a,k,X,Q,R,tol)
  
  Input:
      a -- the new row to be inserted
      k -- the row index for the inserted column
      X -- the orthonormal (column) basis the null space of A before a is inserted
    Q,R -- the output of ApproxiRank
    tol -- the rank-decision threshold
  scale -- the output of ApproxiRank
  
  Output
    ark -- the approximate rank of updated matrix A
      Y -- the updated orthonomal basis for the null space
    R,Q -- same as the output of ApproxiRank for updated A
  
   
• ApproxiRankRowDowndate
  Row downdating module for ApproxiRank.  After computing the approximate rank and kernel of a matrix  A,  if a new row is deleted from  A to form a new matrix  B,   ApproxiRankRowDowndate  calculates the approximate rank/kernel of  B  by updating the approximate rank/kernel of  A. 
  
  Reference:   A rank-revealing method with updating, downdating and applications, T.Y. Li  and  Z. Zeng,  SIAM Journal on Matrix Analysis and Applications, 26(2005),  pp918-946.
  
  Syntax
      >>[ark,Y] = ApproxiRankRowDowndate(k,X,Q,R,tol)
      >>[ark,Y,R,Q] = ApproxiRankRowDowndate(k,X,Q,R,tol)
  
  Input:
       k -- the row index for the row to be deleted
       X -- the orthonormal (column) basis the null space of A before a is deleted
     Q,R -- the output of ApproxiRank
     tol -- the rank-decision threshold
   scale -- the output of ApproxiRank
  
  Output
    ark -- the approximate rank of updated matrix A
      Y -- the updated orthonomal basis for the null space
    R,Q -- same as the output of ApproxiRank for updated A
  
   
• ApproxiRankColumnUpdate
  Columnwise updating module for ApproxiRank.  After computing the approximate rank and kernel of a matrix  A,  if a new column is inserted into  A to form a new matrix  B,   ApproxiRankColumnUpdate  calculates the approximate rank/kernel of  B  by updating the approximate rank/kernel of  A. 
  
  Reference:   A rank-revealing method with updating, downdating and applications, T.Y. Li  and  Z. Zeng,  SIAM Journal on Matrix Analysis and Applications, 26(2005),  pp918-946.
  
  Syntax
      >>[ark,Y] = ApproxiRankColumnUpdate(a,X,Q,R,tol,scale)
      >>[ark,Y] = ApproxiRankColumnUpdate(a,X,Q,R,tol,scale,k)
      >>[ark,Y,R,Q] = ApproxiRankColumnUpdate(a,X,Q,R,tol,scale)
      >>[ark,Y,R,Q] = ApproxiRankColumnUpdate(a,X,Q,R,tol,scale,k)
  
  Input:
       k -- the row index for the column to be inserted
       X -- the orthonormal (column) basis the null space of A before a is inserted
     Q,R -- the output of ApproxiRank
     tol -- the rank-decision threshold
   scale -- the output of ApproxiRank
  
  Output
    ark -- the approximate rank of updated matrix A
      Y -- the updated orthonormal basis for the null space
    R,Q -- same as the output of ApproxiRank for updated A
  
   
• ApproxiRankColumnDowndate
  Columnwise downdating module for ApproxiRank.  After computing the approximate rank and kernel of a matrix  A,  if a new column  is deleted from  A  to form a new matrix  B,   ApproxiRankColumnDowndate  calculates the approximate rank/kernel of  B  by updating the approximate rank/kernel of  A. 
  
  Reference:   A rank-revealing method with updating, downdating and applications, T.Y. Li  and  Z. Zeng,  SIAM Journal on Matrix Analysis and Applications, 26(2005),  pp918-946.
  
  Syntax
      >>[ark,Y] = ApproxiRankColumnDownDate(k,X,Q,R,tol)
      >>[ark,Y,R,Q] = ApproxiRankColumnDowndate(k,X,Q,R,tol)
  
  Input:
       k -- the row index for the column to be deleted
       X -- the orthonormal (column) basis the null space of A before a is deleted
     Q,R -- the output of ApproxiRank
     tol -- the rank-decision threshold
   scale -- the output of ApproxiRank
  
  Output
    ark -- the approximate rank of updated matrix A
      Y -- the updated orthonormal basis for the null space
    R,Q -- same as the output of ApproxiRank for updated A
  
   
• ApproxiJordanForm
  Computing an Approximate Jordan decomposition
            X*J*inv(X)
  of a given matrix A within a distance threshold theta, formulated under a 'three-strikes principle'. For details, see
  
      Z. Zeng and T.Y. Li, A numerical algorithm for computing the Jordan
            Canonical Form, preprint, 2007
  
  In a nutshell, let A be a matrix whose entries are given approximately
  with error magnitude being small. The exact JCF of A  will be degraded.
  However, it is possible to recover the underlying Jordan structure and
  approximate X and J by ApproxiJoranForm
  
  Syntax: There are several choices for either LHS or RHS
      >>    LHS       =      RHS
  ------------------- | --------------------------------------
       [J,X]          |   ApproxiJordanForm(A)
       [J,X,e,s,t]    |   ApproxiJordanForm(A,theta)
                      |   ApproxiJordanForm(A,theta,tau,gap)
  
  Input:   A -- matrix whose AJCF is to be computed
       theta -- (optional) distance threshold
         tau -- (optional) deflation threshold, simple eigenvalues
                   whose geometric condition numbers above tau will be deflated
         gap -- (optional) singular value gap used in rank revealing
  
  Output: J -- The Approximate Jordan Canonical Form
          X -- The principle vector matrix
          e -- the list of distinct eigenvalues
          s -- The Jordan block sizes of the eigenvalues
          t -- the residuals (backward errors) and the staircase condition
                   numbers of the eigenvalues
  
  Example:  One can construct a test matrix with a given Jordan Form by
  
  >>  A = MatrixWithGivenJordanForm([2 2 2 3 3],[3 2 1 2 2])
  A =
      0   5   0  -3   0   0  -1   0   0  -2
    -10   1  10   0   9   0   1  -4  -3   0
     -2   5   2  -3   0   0  -1   0   0  -2
    -17   0  16   1  15   0   1  -7  -5  -1
      6  -1  -5   1  -3   0   0   3   2   1
      0   0   0   0   0   2   0   0   0   0
     -7  -2   7   0   6   0   4  -2  -2   0
      6   0  -6   0  -6   0   0   6   2   0
      9  -3  -6   3  -6   0   0   3   5   3
     10  -7  -6   5  -6   0   1   3   2   6
  
  Its AJCF can be computed as follows:
  
  >> [J,X,e,s,t] = ApproxiJordanForm(A);
  >> J
  
  J =
     3.0000  1.0835  0       0       0       0       0       0       0       0
     0       3.0000  0       0       0       0       0       0       0       0
     0       0       3.0000  0.9464  0       0       0       0       0       0
     0       0       0       3.0000  0       0       0       0       0       0
     0       0       0       0       2.0000  0.6457  0       0       0       0
     0       0       0       0       0       2.0000 -3.5426  0       0       0
     0       0       0       0       0       0       2.0000  0       0       0
     0       0       0       0       0       0       0       2.0000 -1.1010  0
     0       0       0       0       0       0       0       0       2.0000  0
     0       0       0       0       0       0       0       0       0       2.0000
  
  The distinct eigenvalues:
  
  >> e
  
  e =
      3.0000
      2.0000
  
  The Jordan block sizes
  
  >> s
  
  s =
      2  2  0
      3  2  1
  
  The backward error
  
  >> t(:,1)
  
  ans =
  
  1.0e-016 *
  
      0.4897
      0.6711
  
  The staircase condition numbers
  
  >> t(:,2)
  
  ans =
  
  33.5184
  78.4598
   
• uvGCD
  Computing the approximate greatest common divisor (AGCD) of a univariate polynomial pair (f,g) within a distance threshold q, denoted by  u=GCD_q(f,g). 
  
  Computing exact GCD is an ill-posed problem whose solution will generically degrade drastically by a tiny perturbation on the given polynomial coefficients.  The AGCD is, on the other hand, continuously dependent on the problem data, and  uvGCD  can accurately computes the underlying GCD even if the given polynomial pair is inexact. 
  
  The choice of parameter q:  Let (f,g) be the given polynomial pair that approximates an underlying pair (p,q) whose GCD is to be computed.  Then the threshold q needs to satisfy     ||(f,g) - (p,q)||  < q < t   where t  is the distance from  (f,g)  to the nearest polynomial pair (r,s) whose GCD has a higher degree.  In other words, the lower bound of q  is the magnitude of data error and the upper bound is generally an unknown number.  The rule of thumb is that to choose  q  slightly larger than the magnitude of data error.  For example, if the coefficients of (f,g) are accurate to around machine precision (2.2e-16),  then it is likely safe to set  q  = 1e-10.

  See the reference:  Zhonggang Zeng,  The approximate GCD of inexact polynomials, I: a univariate algorithm, Preprint (2004)

  Example:

  >> f
  
  f =
      3 8 14 20 11 4
  
  >> g
  
  g =
      5 14 26 40 30 20 11 4

  >> [u,v,w,res,cond] = uvGCD(f,g,1.0e-10)
  
  u =
      8.30662386291807
      16.61324772583615
      24.91987158875423
      33.22649545167231
  
  v =
      0.36115755925731
      0.24077170617154
      0.12038585308577
  
  w =
      0.60192926542885
      0.48154341234308
      0.36115755925731
      0.24077170617154
      0.12038585308577
  
  res =
      1.614869854000228e-016
  
  cond =
      45.53195717395850
   

  Here,   the syntax:   >> [u, v, w,res, cond] = uvGCD(f,g,theta)

  Input:

   f, g  -- coefficient vectors of polynomial f and g
  theta  -- tolerance within which maximum degree GCD is sought

  Output:

  u    --  the approximate GCD
  v    --  the cofactor for f
  w    --  the cofactor for g
  res  -- the backward error, or residual
  cond -- the condition number of the approximate GCD

  
  More specifically,  for given  (f,g) and q,  the module uvGCD calculates a triplet (u,v,w) of polynomials satisfying the following "three-strikes" principles:
       (a)  Backward nearness:   ||(f,g) - (uv,uw)|| < q
       (b)  Maximum co-dimension:  Polynomial pairs (uv,uw) belongs to the manifold  P_k  = {(p,q) | GCD(p,q) is of degree k}  that has the highest co-dimension amoung all such manifolds intersecting the q-neighborhood of (f,g).
       (c)  Minimum distance:  The pair (uv,uw) is the nearest to (f,g) amoung all polynomial pairs in P_k.
  The polynomial   u=GCD(uv,uw)  is called an AGCD of  (f,g).
  
   

• ExtendedGCD
  ExtendedGCD computes the extended GCD in approximate sense.  For a given polynomial pair (f,g), it computes a polynomial pair (p,q)  such that   p f + qg = GCD(f,g)  in approximate sense.
  
  Calling syntax:       >> [p, q] = ExtendedGCD(v,w);
  
  Input:  v, w -- approximate cofactors, as in the output of
                      >>[u,v,w] = uvGCD(f,g,tol)
  
  Output  p, q -- polynomials such that p*f+q*g = GCD(f,g) approximately
  
  Example:
  
  >> f = [3 8 14 20 11 4];          % set polynomials f and g
  >> g = [5 14 26 40 30 20 11 4];
  >> [u,v,w,res,cond] = uvGCD(f,g); % calculate the GCD triplet (u,v,w)
  >> [p,q] = ExtendedGCD(v,w)       % calculate the extended GCD using v and w
  
  p =
      -4.75485272428887
      -0.98068837438458
      -0.89153488580416
      -1.18871318107222
  
  q =
      2.85291163457332
      0.20802480668764
  
  >> h = conv(p,f)+conv(q,g);    % compute p*f+q*g
  >> h'
  
  ans =
      0.00000000000000
      0.00000000000001
      0
      0.00000000000001
      0
     -0.98068837438458
     -1.96137674876916
     -2.94206512315374
     -3.92275349753832
  
  >> u                           % compare with u=GCD(f,g)
  
  u =
     -0.98068837438458
     -1.96137674876916
     -2.94206512315374
     -3.92275349753832
   
• uvFactor
  The module  uvFactor  computes an approximate factorization     (a₁ x  -  b₁ )^m1   (a₂ x  -  b₂ )^m2   ... (a_k x  -  b_k )^mk     of a given polynomial  f(x)  within a given tolerance  θ  for data perturbation. 
  
  When the coefficients of  f(x)  are approximate (inexact) and the multiplicities  m_j's are nontrivial, finding such a factorization has been a long standing challenge in numerical computation since the factorization problem in exact sense is ill-posed.  Namely, a tiny perturbation on the coefficients of the given polynomial  f  almost always degrades the accuracy of the factors and multiplicities drastically.   The module  uvFactor  in  Apalab,  however, is designed to calculate the factors and multiplicities accurately even if the input polynomial is perturbed.
  
  Syntax:  >> [F,res] = uvFactor(f, theta )
           >> [F,res] = uvFactor(f, theta, 1 )
  
  Input:     p -- coefficient vector of the polynomial to be factored
           tol -- (optional) backward error tolerance
     showroots -- (optional, 0 or 1) showing roots or not
  
  Output     F -- kx3 matrix containing factors with
                    each row [a_i , b_i ,m_i ] representing a factor (a_ix + b_i)^mi
  
  Example:   We construct polynomial    p(x)  =  ( x-1 )⁴⁰  ( x-2 )³⁰ ( x-3 )²⁰  ( x-4 )¹⁰    with coefficients round to 16 digits.
  
  >> p = poly([ones(1,40),2*ones(1,30),3*ones(1,20),4*ones(1,10)]);
  >> [F,res] = uvFactor(p,1e-9)
  
  F =
      0.51445747892761   -2.05782991571042   10.00000000000000
      0.67077048679987   -2.01231146039961   20.00000000000000
      0.94861271967198   -1.89722543934395   30.00000000000000
      1.49988840579248   -1.49988840579248   40.00000000000000
  
  
  res =
      8.584956337989512e-016
  
  The roots and multiplicities can be shown by
  
  >> [-F(:,2)./F(:,1),F(:,3)]
  
  ans =
  
      4.00000000000002    10.00000000000000
      2.99999999999998    20.00000000000000
      2.00000000000001    30.00000000000000
      1.00000000000000    40.00000000000000
  
  Or, use the root-showing flag parameter
  
  >> [F,res] = uvFactor(p,1e-9,1);
  
  THE CONDITION NUMBER: 6.70133
  THE BACKWARD ERROR: 3.51e-015
  THE ESTIMATED FORWARD ROOT ERROR: 4.70e-014
  
  FACTORS
  
      ( x - 3.999999999999993 )^10
      ( x - 3.000000000000007 )^20
      ( x - 1.999999999999998 )^30
      ( x - 1.000000000000000 )^40
   
• mvGCD
  Computing the approximate greatest common divisor (AGCD) of a multivariate polynomial pair (f,g) within a distance threshold q, denoted by  u=GCD_q(f,g). 
  
  Computing exact GCD is an ill-posed problem whose solution will generically degrade drastically by a tiny perturbation on the given polynomial coefficients.  The AGCD is, on the other hand, continuously dependent on the problem data, and  mvGCD  can accurately computes the underlying GCD even if the given polynomial pair is inexact. 
  
  The choice of parameter q:  Let (f,g) be the given polynomial pair that approximates an underlying pair (p,q) whose GCD is to be computed.  Then the threshold q needs to satisfy     ||(f,g) - (p,q)||  < q < t   where t  is the distance from  (f,g)  to the nearest polynomial pair (r,s) whose GCD has a higher degree.  In other words, the lower bound of q  is the magnitude of data error and the upper bound is generally an unknown number.  The rule of thumb is that to choose  q  slightly larger than the magnitude of data error.  For example, if the coefficients of (f,g) are accurate to around machine precision (2.2e-16),  then it is likely safe to set  q  = 1e-10.

  Reference:  Z. Zeng and B.H. Dayton,  The approximate GCD of inexact polynomials, II: a multivariate algorithm, Proc. of ISSAC 04, pp320-327, 2004.

  Syntax:    >> [u,v,w,res,cond] = mvGCD(p, q, theta, ctrl)

  where
        Input:   p, q�--- coefficient matrices of p and q
  ������������  theta ---�(optional) the residual tolerance
                 ctrl�---�(optional) control parameter,
                            if ctrl = 1 then refinement will be performed.
                            otherwise the code stops without refinement.
  
        Output:     u --- the AGCD of p and q
  �������������� v, w --- the co-factors
                  res --- the residual
                c ond --- the AGCD condition number (if ctrl = 1 ).
  
   

  Example:

  f =
  ���� 0 �1 �2 �3� 4� 0 �1 ��2 ��0 �1 �2 �0� 0
  ���� 0 �0 �0 �0 �0� 1 �1 ��1 ��2 �2 �2 �3 �4
  ���� 4� 4 -3 -2� 1 12 �6 �-6 �13 �2 -2 �6 �1

  g =
  ���� 0 ��1�� 2� 3� 4� �0� �1� �2 �3� 0 ��1 �2� 0� 1 �0
  ���� 0 ��0 ��0� 0� 0 ��1� �1 ��1� 1� 2 ��2 �2 �3 �3 �4
  ���� 4 �12 �13� 6 �1 �12 �26 �18 �4 13 �18 �6 �6 �4 �1

  >> [u,v,w,res,cond] = mvGCD(f,g,1e-10,1)

  u =
  ����������� 0� 1.0000e+000� 2.0000e+000����������� 0� 1.0000e+000����������� 0
  ����������� 0����������� 0����������� 0� 1.0000e+000� 1.0000e+000� 2.0000e+000
  � 1.2060e+000� 2.4120e+000� 1.2060e+000� 2.4120e+000� 2.4120e+000� 1.2060e+000

  v =
  ����������� 0� 1.0000e+000� 2.0000e+000����������� 0� 1.0000e+000����������� 0
  ����������� 0����������� 0����������� 0� 1.0000e+000� 1.0000e+000� 2.0000e+000
  � 3.3167e+000 -3.3167e+000� 8.2917e-001� 3.3167e+000 -1.6583e+000� 8.2917e-001

  w =
  ����������� 0� 1.0000e+000� 2.0000e+000����������� 0� 1.0000e+000����������� 0
  ����������� 0����������� 0����������� 0� 1.0000e+000� 1.0000e+000� 2.0000e+000
  � 3.3167e+000� 3.3167e+000� 8.2917e-001� 3.3167e+000� 1.6583e+000� 8.2917e-001

  res =
  � 5.7499e-016

  cond =
  � 9.9903e+000

   

• SquarefreeFactor
  For a given multivariate polynomial  p,  the module  SquarefreeFactor  computes a squarefree factorization
                 p₀  (p₁)^m1  (p₂)^m2  ... (p_k)^mk
  of   p  along with multiplicities  m₁ , m₂ , ... m_k  accurately even if the polynomial is perturbed.   Here  p₀  is a constant factor and  p_j 's are nontrivial polynomial factors for  j > 0.
  
   Syntax:   >> [F,res,cond] = SquarefreeFactor(p,tol);
  
   Input   p -- a multivariate polynomial (as a coeff. matrix) to be factored
         tol -- (optional) backward error tolerance
  
   Output  F -- cell of size kx2. For each j = 1, 2, ..., k, F{j,1}, F{j,2}
                   are a squarefree factor and its multiplicity respectively
         res -- residual
        cond -- condition number of the factorization
  
  Example:
  
  We construct a factorable polynomial
  
  >> p = [2 1 0 0; 0 1 2 0; 1 1 -1 1];
  >> q = [3 0 3 0; 2 2 0 0; 1 1 -2 -2];
  >> r = [3 0 1 0; 0 3 2 0; 1 -1 -3 2];
  >> u = mvPolynPower(p,4);
  >> v = mvPolynPower(q,3);
  >> w = mvPolynMultiply(u,v);
  >> f = mvPolynMultiply(w,r);
  >> [F,res,cond] = SquarefreeFactor(f,1e-10)
  
  F =
  
      [3x1 double] [1]
      [3x4 double] [1]
      [3x4 double] [3]
      [3x4 double] [4]
  
  
  res =
  
      9.094947017729282e-013
  
  
  cond =
  
      0.00417001580803
  
  Factors can be retrieved as follows:
  
  >> F{1,1}  % the constant factor
  
  ans =
  
      0
      0
      2.32164695714406
  
  >> F{2,1}  % factor of multiplicity F{2,2} = 1
  
  ans =
  
      0                  3.00000000000000  1.00000000000000  0
      0                  0                 2.00000000000000  3.00000000000000
     -1.19889333343897  -0.59944666671943  1.79834000015840  0.59944666671949
  
  >> F{3,1}  % factor of multiplicity F{3,2} = 3
  
  ans =
  
      0                 3.00000000000000   0                  3.00000000000000
      0                 0                  2.00000000000000   2.00000000000000
      1.46833846147498  1.46833846147490  -0.73416923073738  -0.73416923073740
  
  >> F{3,1}
  
  ans =
  
      0                 2.00000000000000   1.00000000000000   0
      0                 0                  1.00000000000000   2.00000000000000
      1.16082347857208  1.16082347857205   1.16082347857192  -1.16082347857206
  
  
  
  Even if one make a substantial perturbation, a squarefree factorization can still be calculated up to the data accuracy:
  
  >> g(3,:) = f(3,:)+1e-6*rand(1,size(f,2));   % perturbed polynomial
  >> [F,res,cond] = SquarefreeFactor(g,1e-4,1)
  
  F =
  
      [3x1 double] [1]
      [3x4 double] [1]
      [3x4 double] [3]
      [3x4 double] [4]
  
  
  res =
  
      1.216959617522662e-006
  
  
  cond =
  
      0.00417001629401
  
  >> F{2,1}
  
  ans =
  
      0                  3.00000000000000  1.00000000000000  0
      0                  0                 2.00000000000000  3.00000000000000
     -1.19889332029971  -0.59944665012542  1.79834000569698  0.59944664442427
  
  
   
• mvFactorRefine
  For a given initial factorization (stored in F)
            p₀  (p₁)^m1  (p₂)^m2  ... (p_k)^mk
  of multivariate polynomial p, the module mvFactorRefine attempts to improve the accuracy by applying the Gauss-Newton iteration module GaussNewton.
  
  This module is used as a submodule for the module SquarefreeFactor. 
  
  
  Syntax: >>[G,res,cond] = mvFactorRefine(F, p);
  
     Input:  F -- a kx2 cell containing the factors and multiplicities of the
                        initial factorization of p. For each j = 1, 2, ..., k,
                        F{j,1} is the j-th factor with multiplicity F{j,2}.
             p -- multivariate polynomial in coeff. matrix
  
     Output:
             G -- a cell containing the refined factors in the same format as F
           res -- the residual, i.e., backward error of the factorization
          cond -- the condition number of the factorization
  
   
• Polynomial2CoefficientMatrix
     convert a polynomial to coefficient matrix
     (requires Symbolic Math Toolbox)
     A coefficient matrix F for an m-variate polynomial p of n terms satisfies:
        (1) (m+1) x n
        (2) each term of p is represented by a column of F
        (3) F(m+1,k) is the coefficient of the k-th term of p
        (4) F(j,k) for j=1:m is the exponent for the j-th variable
                in k-th term  

     Syntax:   >> F = Polynomial2CoefficientMatrix(p,u)
               >> F = Polynomial2CoefficientMatrix(p,[x,y,z])

     Input:
                 p -- polynomial
      u or [x,y,z] -- vector of variables

     Oupput:     F -- coefficient matrix

     Example:

    >> syms x y z  % >> p = 3 + x^3*y^2 - 4*y*z^5 + 8*x^2*y^4*z^4
    >> F = Polynomial2CoefficientMatrix(p,[x,y,z])

    F =

        0     3     0     2
        0     2     1     4

        0     0     5     4

        3     1    -4     8

• mvPolynDegree
  Computing the total or tuple degree
   
   Syntax: >> d = mvPolynDegree(p,'total')
           >> d = mvPolynDegree(p,'tuple')
  
  Example:
  
  >> f
  
  f =
     0   3   6   0   3   0   0   3   0   0
     0   0   0   3   3   6   0   0   3   0
     0   0   0   0   0   0   3   3   3   6
    -2  -1   1  -1   2   1  -1   2   2   1
  
  >> d = mvPolynDegree(f,'total')
  
  d =
      6
  
  >> d = mvPolynDegree(f,'tuple')
  
  d =
      6   6   6
  
   
• mvPolynTimesMonom:
  Compute a multivariate polynomial p times a monomial t
  
  Syntax: >> f = mvPolynTimesMonom(p,t)
  
  Example:
  
  >> f
  
  f =
     0   3   6   0   3   0   0   3   0   0
     0   0   0   3   3   6   0   0   3   0
     0   0   0   0   0   0   3   3   3   6
    -2  -1   1  -1   2   1  -1   2   2   1
  
  >> t = [2,3,1,5]'
  
  t =
     2
     3
     1
     5
  
  >> mvPolynTimesMonom(f,t)
  
  ans =
     2   5   8   2   5   2   2   5   2   2
     3   3   3   6   6   9   3   3   6   3
     1   1   1   1   1   1   4   4   4   7
   -10  -5   5  -5  10   5  -5  10  10   5
   
• mvPolynAdd:
  Compute the addition f+g of two multivariate polynomials f and g,
  or the linear combination s*f + t*g if scalars s and t are also given
  
  Syntax: >> d = mvPolynAdd(f,g)
          >> d = mvPolynAdd(f,g,s,t)
  
  Input: f, g -- multivariate polynomials in coeff. matrices
         s, t -- (optional) complex numbers
  
  Output: multivariate polynomial f+g, or s*f + t*g
  
  Example:
  
  >> f
  
  f =
     0   2  4   0  2  0
     0   0  0   2  2  4
    -2  -1  1  -1  2  1
  
  >> p
  
  p =
     3   5   7   3   5   3
     1   1   1   3   3   5
     8  12   4  12   8   4
  
  >> r = mvPolynAdd(f,p,1,-1) % compute r=f-p
  
  r =
      0   2   4   3   5   7   0   2   3   5   0   3
      0   0   0   1   1   1   2   2   3   3   4   5
     -2  -1   1  -8 -12  -4  -1   2 -12  -8   1  -4
   
• mvPolynMultiply:
  Compute the product f*g of two multivariate polynomials f and g
   
  Syntax: >> h = mvPolynMultiply(f,g)
  Input: f, g -- multivariate polynomials in coeff. matrices
  
  Output: multivariate polynomial f*g
  
  Example:
  
  >> f
  
  f =
     0   2   4   0
     0   0   0   2
    -4  -2   2  -2
  
  >> g
  
  g =
     3   5   7
     1   1   1
     8  12   4
  
  >> mvPolynMultiply(f,g)
  
  ans =
       3    5    7    9   11    3    5    7
       1    1    1    1    1    3    3    3
     -32  -64  -24   16    8  -16  -24   -8
   
• mvPolynPower
  Compute the power f^k of the multivariate polynomials f
  
  Syntax:   >> h = mvPolynPower(f,k)
  
  Input:  f -- multivariate polynomials in coeff. matrices
          k -- integer exponent of the power
  
  Output: multivariate polynomial f^k
  
  Example:
  
  >> f
  
  f =
      0   2   4   0
      0   0   0   2
     -4  -2   2  -2
  
  >> mvPolynPower(f,3)
  
  ans =
      0    2    4    6    8   10   12    0    2    4    6    8    0    2    4    0
      0    0    0    0    0    0    0    2    2    2    2    2    4    4    4    6
    -64  -96   48   88  -24  -24    8  -96  -96   72   48  -24  -48  -24   24   -8
  
  
   
• mvPolynEvaluate:
   Evaluate multivariate polynomial f at x
  
      Syntax:   >> y = mvPolynEvaluate(f,x)
  
      Input:   f -- multivariate polynomial (as a coeff. matrix)
                   x -- vector value for the variables
  
      utput: the value of f(x)
  
  Example:
  
  >> f
  
  f =
     0    2    4    0
     0    0    0    2
    -4   -2    2   -2
  
  >> mvPolynEvaluate(f,[1;1])
  
  ans =
  
     -6
  
   
• mvPolynDerivative
  Compute the partial derivative of the multivariate polynomial p
  with respect to the j-th variable
  
  Syntax: >> q = mvPolynDerivative(p,j)
  
  Input:
  p -- polynomial in coefficient matrix
  j -- index of the variable w.r.t. which the derivative is to be taken
  
  Output: multivariate polynomial as the partial derivative
  
  Example: 
  
  >> p = [2 1 0 0; 0 1 2 0; 1 1 -1 1]
  
  p =
     2   1   0   0
     0   1   2   0
     1   1  -1   1
  
  >> mvPolynDerivative(p,1)
  
  ans =
      1   0
      0   1
      2   1
  
  >> mvPolynDerivative(p,2)
  
  ans =
      1   0
      0   1
      1  -2
  
   
• mvPolynClear
  Clear tiny/zero coefficients of a multivariate polynomial f
  
  Syntax:
       >> g = mvPolynClear(f)         % clear zero coefficients of f
       >> g = mvPolynClear(f,epsilon) % clear coefficients whose magnitudes
  are less than or equal to epsilon
  
  Input:    f -- multivariate polynomial (as a coeff. matrix)
      epsilon -- magnitude threshold of coefficients to be cleared
  
  Output: multivariate polynomial in coefficient matrix
  
  Example:
  
  >> h
  
  h =
      3    5    7    9   11    3    5    7
      1    1    1    1    1    3    3    3
    -32    0  -24   16    0  -16    0   -8
  
  >> mvPolynClear(h)
  
  ans =
      3    7    9    3    7
      1    1    1    3    3
    -32  -24   16  -16   -8
  
  >> g
  
  g =
    3.0000e+000  5.0000e+000  7.0000e+000  3.0000e+000  5.0000e+000  3.0000e+000
    1.0000e+000  1.0000e+000  1.0000e+000  3.0000e+000  3.0000e+000  5.0000e+000
    8.0000e+000  9.5013e-013  4.0000e+000  1.2000e+001  8.0000e+000  2.3114e-013
  
  >> mvPolynClear(g,1e-10)
  
  ans =
      3   7   3   5
      1   1   3   3
      8   4  12   8
  
  
   
• HouseholderVector:
  computing the Householder vector u such that
           H = I - 2*u*u'
  makes Hv = [*;0;...;0]
  
  syntax >> u = HouseholderVector(v)
  
  Example: 
  
  >> v = rand(5,1)
  
  v =
     0.8381
     0.0196
     0.6813
     0.3795
     0.8318
  
  >> u = HouseholderVector(v)
  
  u =
     0.8922
     0.0078
     0.2698
     0.1503
     0.3294
  
  >> (eye(5)-2*u*u')*v
  
  ans =
    -1.4152
    -0.0000
    -0.0000
    -0.0000
    -0.0000
   
• HouseholderTransform:
  computing the Householder transformation w = H*v with
       H = I - 2*u*u'
  
  syntax >> w = HouseholderTransform(u,v)
  
  Example:
  
  >> v
  
  v =
     0.8381
     0.0196
     0.6813
     0.3795
     0.8318
  
  >> u = HouseholderVector(v)
  
  u =
     0.8922
     0.0078
     0.2698
     0.1503
     0.3294
  
  >> w = HouseholderTransform(u,v)
  
  w =
    -1.4152
    -0.0000
    -0.0000
    -0.0000
    -0.0000
  
   
• HouseholderTransformRight:
  computing the Householder transformation from right side w = v*H  with
      H = I - 2*u*u'
  
  syntax >> w = HouseholderTransformRight(u,v)
  
  Example:
  
  >> a = v'
  
  a =
     0.8381   0.0196   0.6813   0.3795   0.8318
  
  >> u = HouseholderVector(v);
  >> w = HouseholderTransformRight(u,a)
  
  w =
    -1.4152  -0.0000  -0.0000  -0.0000  -0.0000
  
   
• GivensMatrix:
  Generating a Givens' (2x2) matrix T such that
              T*v = [d; 0]
  
  Input: v -- a vector of dimension 2
  
  Output: T -- Givens' matrix
  d -- the 2-norm of v
  
  syntax >> [T,d] = GivensMatrix(v)
  
  Example:
  
  >> v
  
  v =
     0.7907
    -1.5935
  
  >> [T,d] = GivensMatrix(v)
  
  T =
     0.4445  -0.8958
     0.8958   0.4445
  
  
  d =
     1.7789
  
  >> T*v
  
  ans =
     1.7789
     0.0000
  
   
• IntegerInverse:
  generate a pair of nxn integer matrices that are inverse to each other
  
  syntax >> [X,Y] = IntegerInverse(n)
         >> [X,Y] = IntegerInverse(n,m) % entries in interval [-m,m]
  
  This module is mainly used to generate test cases
  
  Example:
  
  >> [X,Y] = IntegerInverse(5)
  
  X =
     1   0  -1   0   0
     0   1   0  -1  -1
    -1   0   2   0   0
     0  -1   0   2   1
     0  -1   0   1   2
  
  
  Y =
     2   0   1   0   0
     0   3   0   1   1
     1   0   1   0   0
     0   1   0   1   0
     0   1   0   0   1
  
  >> X*Y
  
  ans =
     1   0   0   0   0
     0   1   0   0   0
     0   0   1   0   0
     0   0   0   1   0
     0   0   0   0   1
   
• MatrixWithGivenJordanForm
  Construct a matrix with given Jordan Canonical Form
  
  Syntax: There are several choices for either LHS or RHS
     >> LHS           =           RHS
  ------------------- | --------------------------------------
        A             |   MatrixWithGivenJordanForm(s,j)
       [A,X]          |   MatrixWithGivenJordanForm(s,j,k)
       [A,X,Y]        |
  
  Input:
       s -- eigenvalues in each Jordan block
       j -- size of each Jordan block
       k -- (optional) number of (simple) eigenvalues of a kxk random matrix
  
  Output:
       A -- a matrix whose eigenvalues and Jordan structure are defined by
              input s, j and k
     X,Y -- matrices such that Y*A*X is the Jordan Canonical Form of A
  
  Example:  To generate a matrix with
         eigenvalues | Jordan block sizes
              2      |     4, 3, 2
              3      |     2, 2
  
  >> [A,X,Y] = MatrixWithGivenJordanForm([2 2 2 3 3],[3 2 1 2 2])
  
  A =
     -2   8  -4  -3   1  -1  -1  -1   1   1
      4  -4   5   1   0   0   1   2   1  -2
      8 -15  10   6  -3   2   2   2  -2  -2
     -1   5  -1  -1   4  -2  -1  -1   2   0
      0   0   0   0   2   0   0   0   0   0
      0   2   1  -2   1   1   0   1   2  -1
      2  -6   4   1  -2   1   4   4   1  -2
      0   0   0   0   0   0   0   3   0   0
    -10  21 -11  -7   5  -3  -3  -5   4   4
     -5   8  -8   2  -2   1  -1  -4  -4   7
  
  
  X =
      1   0  -1   0   0   1   0   0   1   0
      0   1   1   0   0   1   0   0  -1  -1
     -1   1   3  -1   0   0   0   0  -3   0
      0   0  -1   2   0   0  -1   0   2  -1
      0   0   0   0   1   1   0   0   0   0
      1   1   0   0   1   4   0   0   0  -1
      0   0   0  -1   0   0   2   0  -1   0
      0   0   0   0   0   0   0   1   0   1
      1  -1  -3   2   0   0  -1   0   5   1
      0  -1   0  -1   0  -1   0   1   1   5
  
  
  Y =
      6  -3   3   3   1  -1   1   0  -1   0
     -3   7  -4  -2   1  -1  -1  -1   0   1
      3  -4   4   1   0   0   1   1   1  -1
      3  -2   1   4   0   0   1  -1  -2   1
      1   1   0   0   2  -1   0   0   0   0
     -1  -1   0   0  -1   1   0   0   0   0
      1  -1   1   1   0   0   1   0   0   0
      0  -1   1  -1   0   0   0   2   1  -1
     -1   0   1  -2   0   0   0   1   2  -1
      0   1  -1   1   0   0   0  -1  -1   1
  
  The results can be verified:
  
  >> Y*A*X
  
  ans =
  
      2   1   0   0   0   0   0   0   0   0
      0   2   1   0   0   0   0   0   0   0
      0   0   2   0   0   0   0   0   0   0
      0   0   0   2   1   0   0   0   0   0
      0   0   0   0   2   0   0   0   0   0
      0   0   0   0   0   2   0   0   0   0
      0   0   0   0   0   0   3   1   0   0
      0   0   0   0   0   0   0   3   0   0
      0   0   0   0   0   0   0   0   3   1
      0   0   0   0   0   0   0   0   0   3
   
• ScaledLeastSquares:
  Scaled least squares solver for A*x = b
  
  Input: A -- matrix
         v -- right-side vector
         w -- (optional) weight vector
  
  Output: (return) -- the least squares solution
  
  syntax >> x = ScaledLeastSquares(A,b)
         >> x = ScaledLeastSquares(A,b,w)