The numerical algorithms

The algorithms underlying the basic MATLAB functions are described in the LINPACK and EISPACK guides [1-3]. The following list gives the subroutines used by these functions.

   INV(A)          - CGECO,CGEDI
   DET(A)          - CGECO,CGEDI
   LU(A)           - CGEFA
   RCOND(A)        - CGECO
   CHOL(A)         - CPOFA
   SVD(A)          - CSVDC
   COND(A)         - CSVDC
   NORM(A,2)       - CSVDC
   PINV(A,eps)     - CSVDC
   RANK(A,eps)     - CSVDC
   QR(A)           - CQRDC,CQRSL
   ORTH(A)         - CQRDC,CSQSL
   \verb!A\B! and \verb!B/A!     - CGECO,CGESL if A is square.
                   - CQRDC,CQRSL if A is not square.
   EIG(A)          - HTRIDI,IMTQL2,HTRIBK if A is Hermitian.
                   - CORTH,COMQR2         if A is not Hermitian.
   SCHUR(A)        - same as EIG.
   HESS(A)         - same as EIG.

Minor modifications were made to all these subroutines. The LINPACK routines were changed to replace the Fortran complex arithmetic with explicit references to real and imaginary parts. Since most of the floating point arithmetic is concentrated in a few low-level subroutines which perform vector operations (the Basic Linear Algebra Subprograms), this was not an extensive change. It also facilitated implementation of the FLOP and CHOP features which count and optionally truncate each floating point operation.

The EISPACK subroutine COMQR2 was modified to allow access to the Schur triangular form, ordinarily just an intermediate result. IMTQL2 was modified to make computation of the eigenvectors optional. Both subroutines were modified to eliminate the machine-dependent accuracy parameter and all the EISPACK subroutines were changed to include FLOP and CHOP.

The algorithms employed for the POLY and ROOTS functions illustrate an interesting aspect of the modern approach to eigenvalue computation. POLY(A) generates the characteristic polynomial of A and ROOTS(POLY(A)) finds the roots of that polynomial, which are, of course, the eigenvalues of A . But both POLY and ROOTS use EISPACK eigenvalues subroutines, which are based on similarity transformations. So the classical approach which characterizes eigenvalues as roots of the characteristic polynomial is actually reversed.

If A is an n by n matrix, POLY(A) produces the coefficients C(1) through C(n+1), with C(1) = 1, in

      DET(z*EYE-A) = C(1)*z**n + ... + C(n)*z + C(n+1) .

The algorithm can be expressed compactly using MATLAB:

      Z = EIG(A);
      C = 0*ONES(n+1,1);  C(1) = 1;
      for j = 1:n, C(2:j+1) = C(2:j+1) - Z(j)*C(1:j);
      C

This recursion is easily derived by expanding the product

      (z - z(1))*(z - z(2))* ... * (z-z(n)) .

It is possible to prove that POLY(A) produces the coefficients in the characteristic polynomial of a matrix within roundoff error of A . This is true even if the eigenvalues of A are badly conditioned. The traditional algorithms for obtaining the characteristic polynomial which do not use the eigenvalues do not have such satisfactory numerical properties.

If C is a vector with n+1 components, ROOTS(C) finds the roots of the polynomial of degree n ,

       p(z) = C(1)*z**n + ... + C(n)*z + C(n+1) .

The algorithm simply involves computing the eigenvalues of the companion matrix:

      A = 0*ONES(n,n)
      for j = 1:n, A(1,j) = -C(j+1)/C(1);
      for i = 2:n, A(i,i-1) = 1;
      EIG(A)

It is possible to prove that the results produced are the exact eigenvalues of a matrix within roundoff error of the companion matrix A, but this does not mean that they are the exact roots of a polynomial with coefficients within roundoff error of those in C . There are more accurate, more efficient methods for finding polynomial roots, but this approach has the crucial advantage that it does not require very much additional code.

The elementary functions EXP, LOG, SQRT, SIN, COS and ATAN are applied to square matrices by diagonalizing the matrix, applying the functions to the individual eigenvalues and then transforming back. For example, EXP(A) is computed by

      <X,D> = EIG(A);
      for j = 1:n, D(j,j) = EXP(D(j,j));
      X*D/X

This is essentially method number 14 out of the 19 'dubious' possibilities described in [8]. It is dubious because it doesn't always work. The matrix of eigenvectors X can be arbitrarily badly conditioned and all accuracy lost in the computation of X*D/X. A warning message is printed if RCOND(X) is very small, but this only catches the extreme cases. An example of a case not detected is when A has a double eigenvalue, but theoretically only one linearly independent eigenvector associated with it. The computed eigenvalues will be separated by something on the order of the square root of the roundoff level. This separation will be reflected in RCOND(X) which will probably not be small enough to trigger the error message. The computed EXP(A) will be accurate to only half precision. Better methods are known for computing EXP(A), but they do not easily extend to the other five functions and would require a considerable amount of additional code.

The expression A**p is evaluated by repeated multiplication if p is an integer greater than 1. Otherwise it is evaluated by

      <X,D> = EIG(A);
      for j = 1:n, D(j,j) = EXP(p*LOG(D(j,j)))
      X*D/X

This suffers from the same potential loss of accuracy if X is badly conditioned. It was partly for this reason that the case p = 1 is included in the general case. Comparison of A**1 with A gives some idea of the loss of accuracy for other values of p and for the elementary functions.

RREF, the reduced row echelon form, is of some interest in theoretical linear algebra, although it has little computational value. It is included in MATLAB for pedagogical reasons. The algorithm is essentially Gauss-Jordan elimination with detection of negligible columns applied to rectangular matrices.

There are three separate places in MATLAB where the rank of a matrix is implicitly computed: in RREF(A), in A\B for non- square A, and in the pseudoinverse PINV(A). Three different algorithms with three different criteria for negligibility are used and so it is possible that three different values could be produced for the same matrix. With RREF(A), the rank of A is the number of nonzero rows. The elimination algorithm used for RREF is the fastest of the three rank-determining algorithms, but it is the least sophisticated numerically and the least reliable. With A\B, the algorithm is essentially that used by example subroutine SQRST in chapter 9 of the LINPACK guide. With PINV(A), the algorithm is based on the singular value decomposition and is described in chapter 11 of the LINPACK guide. The SVD algorithm is the most time-consuming, but the most reliable and is therefore also used for RANK(A).

The uniformly distributed random numbers in RAND are obtained from the machine-independent random number generator URAND described in [9]. It is possible to switch to normally distributed random numbers, which are obtained using a transformation also described in [9].

The computation of $\sqrt{a^2 + b^2}$ is required in many matrix algorithms, particularly those involving complex arithmetic. A new approach to carrying out this operation is described by Moler and Morrison [10]. It is a cubically convergent algorithm which starts with a and b , rather than with their squares, and thereby avoids destructive arithmetic underflows and overflows. In MATLAB, the algorithm is used for complex modulus, Euclidean vector norm, plane rotations, and the shift calculation in the eigenvalue and singular value iterations.