Recent Advances in Linear Algebra Functionality

Zbigniew Leyk

Wolfram Research, Inc.

Abstract

This talk will discuss  sparse matrix algorithms implemented recently in Mathematica. We can always convert a sparse matrix to a corresponding dense matrix and use dense-matrix algorithms. But for large sparse matrices this approach is very inefficient because most computations are performed on elements which are zero. Therefore, sparse matrices require special algorithms to solve linear equations,  find eigenvalues or compute a least squares solution. Such algorithms take advantage of the sparse structure of the matrix. Numerical examples are presented to show that the sparse-matrix algorithms are much faster than their dense-matrix counterparts.

Sparse Matrices

Sparse matrices were introduced in Mathematica Version 5. These are matrices with a small number of  nonzero elements. It is possible to apply general linear algebra methods to sparse matrices, but it is wasteful because most arithmetic operations is performed on zeros. Therefore, it is important to have special methods which are designed to work with sparse matrices.

As an example of a sparse matrix consider this tridiagonal matrix with the number 3 on the main diagonal, 2 on the first superdiagonal, and 1 on the first subdiagonal.

A  =  
      SparseArray[{{i_, i_} -> 
            3. ,  {i_, j_} /; i == j + 1 -> 
            2. ,  {i_, j_} /; i == j - 1 -> 1. },  {100, 
          100}]

The matrix has only 298 nonzero elements (which are stored) and 9702 zero elements (which are not stored).

Singular Value Decomposition

We start with presenting a method which is used for computing singular values and singular vectors of a sparse matrix. Any m×n matrix A can be factored into the form

A = U . S . (V^T)

where U and V are orthonormal matrices of sizes m×m and n×n, respectively, and the diagonal matrix S contains singular values. There is a close relation between singular values of A and eigenvalues of (A^T) . A. From the above relation we get

(A^T)  .  
          A = (V . (S^T) . (U^T)) . U . S . (V^T) = 
          V . (S^T) . S . (V^T)

which means that V is the orthonormal matrix of eigenvectors of (A^T) . A and (S^T) . S is the diagonal matrix containing n (non-negative) eigenvalues.

In many applications one is interested in computing a few, say k, largest or smallest singular values and corresponding vectors. For these few singular values and vectors  the following relations are still valid ((Subscript[u, i]) is the (i^th) column of U, (Subscript[v, i]) is the (i^th) column of U, i=1,…,k):

A  .  (Subscript[v, i]) = (Subscript[u, i])  .  (Subscript[s, i])
   (A^T)  .  (Subscript[u, i]) = (Subscript[v, i])  .  (Subscript[s, i])

Arnoldi Method

Mathematica Version 5 uses the numerical package ARPACK to compute a few eigenvalues and/or eigenvectors of a sparse matrix. ARPACK is a collection of  subroutines designed to solve large scale eigenvalue problems. It is most appropriate for large sparse or structured matrices. ARPACK is based on the Arnoldi method. This package can also be used for computing singular values and vectors of a sparse matrix if we transformation the matrix to a suitable form. These are the two mostly used transformations.

•  One approach uses the fact that singular values are equal to the square root of eigenvalues of the matrix (A^T) . A. With this approach there is a problem with computing left singular vectors (U) — especially those corresponding to zero singular values. 
•  Another approach is to use the (m + n)×(m + n) augmented matrix of the form:
   (({{0, A}, {(A^T), 0}}))
   With this approach we directly compute both the left and right singular vectors as well as singular values as eigenvectors and eigenvalues of the augmented matrix. The only problem we can have is the fact that we now have two sets of singular values and, possibly, additional zero singular values (=m - n, if m ≥ n).

Computing Singular Values and Vectors

To compute singular values/vectors of a sparse matrix in Mathematica Version 5.1 we use the functions SingularValueDecomposition and SingularValueList with three arguments: matrix argument, the number of required singular values/vectors and the Arnoldi method. If the second argument is positive, we want to compute a few largest singular values. If it is negative, we want to compute a few smallest singular values.

SingularValueList[A, 3,  Method -> Arnoldi] // 
      InputForm

The result is almost the same as the one generated by the corresponding dense matrix.

SingularValueList[Normal[A], 3,  Method -> Direct] // 
      InputForm

Computing Singular Values and Vectors (cont)

Unfortunately, the example presented before is not what you will normally see when working with the Arnoldi method. The Arnoldi method is an iterative method and requires many iterations to obtain a result within a given tolerance.

•  The default tolerance = Sqrt[$MachineEpsilon] = 1.49*^-8. 
•  It is often necessary to increase the number of iterative steps (default limit is 1000). 
•  Sometimes it is necessary to decrease/increase the tolerance.

This changes can be done using suboptions in the option Method.

Computing Singular Values and Vectors (cont)

We often need to increase the number of iterative steps and/or decrease tolerance in order to obtain an acceptable result.

A  =  
        SparseArray[{{i_, i_} -> 
              3. ,  {i_, j_} /; i == j + 1 -> 
              2. ,  {i_, j_} /; i == j - 1 -> 1. },  {500, 
            500}];
(ss  = 
          SingularValueList[A, 3, 
            Method -> {Arnoldi,  Tolerance -> 10^-6, 
                MaxIterations -> 10000}]) // InputForm

This is the result generated by the corresponding dense matrix.

(sd  =  SingularValueList[Normal[A], 3]) // 
      InputForm

And this is the relative error between the computed singular values.

Max[Abs[ss  -  sd]]/Max[sd]

Computing Singular Values and Vectors (cont)

Another useful suboption of the option Method is Shift. This suboption is especially helpful when we want to compute a few smallest singular values.

Increasing the number of iterative steps usually does not work well for the smallest singular values.

SingularValueList[A, -3, 
      Method -> {Arnoldi,  MaxIterations -> 10000}]

But when we use a shift, then we can expect to get a better result (note that for Shift we need to write the second argument as 3, not -3).

(ss = 
          SingularValueList[A, 3, 
            Method -> {Arnoldi, Shift -> 0.01}]) // 
      InputForm

Compare this result with the one generated by the corresponding dense matrix.

(sd  =  SingularValueList[Normal[A],  -3]) // 
      InputForm

Max[Abs[ss  -  sd]]/Max[sd]

Least Squares

Another new function which will work with sparse matrices in Mathematica Version 5.1 is LinearAlgebra`LeastSquares. It is used for solving the least squares problems: Find a vector x such that A . x == b where the matrix A is of size m×n and b is of size m. When m > n then the system is overdetermined and usually it has no exact solution. Therefore, we try to minimize 2-norm ‖b - 
          A . x (Subscript[‖, 2]), that is, we try to find x which makes this norm minimal.

If A has full rank, then there is a unique least squares solution (Subscript[x, 0]) which solves the following system (normal equations):

(A^T)  .  A  .  (Subscript[x, 0]) = (A^T)  .  b

For machine precision and arbitrary precision the normal equations are often ill-conditioned. Therefore, we cannot use this form to compute a least squares solution. We need to use other methods to solve it — most of them is based on QR or SVD decompositions. 
Moreover, if A is sparse, then the product (A^T) . A can be less sparse or even dense.

Augmented System

This is an explanation how LeastSquares work. For (machine and arbitrary precision) sparse matrices, to get the least squares solution the augmented linear system of the following form is built ((Subscript[I, m]) = identity matrix):

(({{(Subscript[I, m]), A}, {(A^T), 0}})) . (({{r}, {x}}))=(({{b}, {0}}))

where r = b - A . x is the residual. Then this system is solved using LinearSolve.

Machine-Precision Sparse Matrices

Consider the following tridiagonal matrix of machine precision.

m = 2000;  n = 1000;

    A  = 
        SparseArray[{{i_, j_} :>  
              Random[] /; Abs[i - j] ≤  1},  {m, 
            n}];

    b = Table[Random[], {m}];

In version 5, we could use PseudoInverse to get a least squares solution. Now, in version 5.1, we can also use LeastSquares. Here we compare their timings.

{Timing[x1 = PseudoInverse[A] . b; ], 
        Timing[x = LinearAlgebra`LeastSquares[A, b]]} // 
      Map[First,  #] &

As we see LeastSquares is hundreds times faster than PseudoInverse. This is because LeastSquares can work with sparse matrices whereas PseudoInverse cannot. Moreover, the results generated by both the functions are close.

{Max[Abs[Transpose[A] . (b - A . x)]], 
Max[
        Abs[Transpose[A] . (b - A . x1)]], 
Max[
        Abs[x - x1]]}

Arbitrary-Precision Sparse Matrices

We can use LeastSquares to get least squares solutions for arbitrary-precision matrices. This is an example.

m = 100;  n = 50;
A = 
        SparseArray[
          Table[Random[Real,  {-100, 100}, 
              20], {m}, {n}]];
b = Table[Random[Real,  {-100, 100}, 20], {m}];

Compare timings of LeastSquares and PseudoInverse for this matrix.

{Timing[x = 
              LinearAlgebra`LeastSquares[A, 
                b];], 
Timing[x1 = 
              PseudoInverse[A] . b;]} // Map[First, #] &

This time an explanation why LeastSquares is faster than PseudoInverse is different.  In the above example LeastSquares uses the iterative refinement method (see the next slide), whereas PseudoInverse cannot use it.  The results computed by the two functions are close.

{Max[Abs[Transpose[A] . (b - A . x)]], 
Max[
        Abs[Transpose[A] . (b - A . x1)]], 
Max[
        Abs[x - x1]]}

Iterative Refinement Method

For some linear systems of equations of arbitrary precision we can use the iterative refinement (or iterative improvement) method. In this method most computations is performed in machine precision. Only two operations are carried out in arbitrary precision: 
•  the matrix by vector product and 
•  adding two vectors

In total, there are O((n^2)) arithmetic operations in arbitrary precision for iterative refinement. We do not count arithmetic operations in machine precision because they are much faster than the corresponding operations in arbitrary precision.
Because solving an arbitrary-precision linear system of equations requires O((n^3)) arithmetic operations, iterative refinement will always be faster for larger matrices than any decomposition method.

The only problem with the iterative refinement method is that it works only for linear systems with small or moderate condition numbers. Ill-conditioned linear systems of equations still require Gaussian elimination done in arbitrary precision.

Example

Consider the following symmetric tridiagonal system of precision 120.

n = 1000;
mat = 
        SparseArray[{{i_, i_} -> -2`120, {i_, j_} /; 
                Abs[i - j] == 1 -> 1`120},  {n, 
            n}];
b  =  Table[1, {n}];
Timing[x = 
            LinearSolve[mat, b,  Method -> IterativeRefinement];] // 
      First

We check that this method computes the correct solution.

Max[Abs[mat . x - b]]

If we use the direct method, we get the same result but it will take much more time.

Timing[x = 
          LinearSolve[mat, b, 
            Method -> Direct];]
Max[Abs[mat . x - b]]

Therefore, for this matrix, the iterative refinement method is automatically chosen by LinearSolve[A,b].

Example (cont)

The iterative refinement method cannot solve linear systems with large condition numbers.

<< LinearAlgebra`MatrixManipulation`
n = 20;  prec = 30;  SeedRandom[322];
mat = N[HilbertMatrix[n],  prec];
b  =  
        Table[Random[Real,  {-1, 1},  
            prec],  {n}];
(x = LinearSolve[mat, b, Method -> IterativeRefinement]) // 
      Head

We check that the condition number of the matrix is large.

{LinearAlgebra`MatrixConditionNumber[mat] // 
        N, 
LinearAlgebra`MatrixConditionNumber[
        N[mat]]}

Example (cont)

Also, for some small matrices, the iterative refinement method can be slower than the direct method as it can be seen on this example.

n = 40;  prec = 200;  SeedRandom[123];
mat = 
        Table[Random[Real,  {-1, 1},  
            prec],  {n}, {n}];
mat[[n]] = 
        mat[[1]] + 
          Table[Random[Real, {-1, 1}, 
                prec] 10^-13, {n}];
b  =  
        Table[Random[Real, {-1, 1}, 
            prec], {n}];
{Timing[x1 = 
              LinearSolve[mat, b, 
                Method -> 
                  IterativeRefinement];], 
Timing[x2 = 
              LinearSolve[mat, b, Method -> Direct];]} // 
      Map[First, #] &

LinearAlgebra`MatrixConditionNumber[mat] // N