Fast Linear Algebra In Mathematica

Zbigniew Leyk (zleyk@wolfram.com)

Plan of the talk

LAPACK

•  LAPACK is a (dense) linear algebra package of Fortran 77 routines for solving the most commonly occurring problems in numerical linear algebra.  
•  It was written from scratch in late 1980s as a replacement for two other linear algebra packages: LINPACK and EISPACK. 
•  It was designed for achieving greater efficiency and robustness on modern computers.
•  Current version is 3.0.
•  C version is also available (it is called CLAPACK). This version is used in Mathematica.

Advantages of LAPACK:

•  implemented using BLAS (basic linear algebra subroutines) which can be highly optimized. Such highly-optimized BLAS versions exist for many platforms (among others for Windows and Linux Intel 32-bit machines)
•  some LAPACK routines can also be highly optimized
•  makes use of modern linear algebra algorithms
•  improves on the accuracy of the standard algorithms comparing with EISPACK. Among others, it implemented high accuracy algorithms for finding singular values and eigenvalues.
•  contains routines for generalized functions (among others, generalized eigenproblem and generalized singular value decomposition)
•  designed to be used for serious computations, especially as a source of building blocks for larger applications
•  portable over a large class of modern computers

LAPACK and Machine-Precision Linear Algebra Functions

•  We use unmodified BLAS and LAPACK packages in Mathematica 5. They replaced LINPACK and EISPACK packages previously used in the kernel up to and including the version 4.2.
•  We use highly optimized BLAS libraries for platforms for which they exist (for instance, Intel's MKL library on 32-bit Intel machines).
•  If a highly optimized BLAS does not exist for a given platform, we generate ATLAS BLAS for this platform (ATLAS is a self-optimizing BLAS generator that produces BLAS automatically tuned for a given machine). For all supported platforms, ATLAS achieves performance on par with machine-specific tuned libraries.
•  If a highly optimized library contains optimized LAPACK routines, we also use them in place of unoptimized ones.

This is the list of Mathematica 5 linear algebra functions which use machine-precision BLAS/LAPACK routines (in parentheses).
•  CholeskyDecomposition (potrf)
•  Dot (BLAS 1,2,3: dot, gemv, gemm)
•  Det (getrf)
•  Eigensystem, Eigenvalues, Eigenvectors (geev)
•  Developer`HessenbergDecomposition (gehrd, dorghr/zunghr, gebal)
•  Inverse (getrf, getri, gecon)
•  LinearSolve (gesv, gelss, trcon, potv)
•  LUDecomposition (getrf, gecon)
•  LUBackSubstitution (getrs)
•  NullSpace (gesdd)
•  PseudoInverse (gesdd)
•  QRDecomposition (geqp3, geqrf, dorgqr, zungqr)
•  RowReduce (BLAS 1, 2: swap, scal, dger/zgeru)
•  SchurDecomposition (gees, gebal)
•  SingularValues, SingularValueDecomposition, SingularValueList (gesdd)

Are They Really Faster?

These are timings of some linear algebra functions available in Mathematica 5 and the corresponding functions in Mathematica 4.2.  Timings provided here were obtained on a Pentium 3 (Coppermine) machine with 500 MHz CPU. On this machine Mathematica 5 uses Intel's MKL library (5.2sp1). The matrices a and b, which are used in tests, are random matrices of size 1000×1000. The vector v has length 1000.
The functions using BLAS3  are roughly 7-8 times faster in the version 5, and the functions using BLAS2 are about 2-3 times faster. SingularValues uses a new faster algorithm (divide-and-conquer, gelsd) that is why it is almost 6 times faster.

({{                 Function,    V4  .2, V5}, {Dot[a,  b] (BLAS3), 48.04  s, 5.21  s}, {Dot[a, v] (BLAS2), 0.049  s, 0.015  s}, {Det[a] (BLAS3), 15.9  s, 1.96  s}, {LinearSolve[a, v] (BLAS3), 16.1  s, 2.15  s}, {Inverse[a] (BLAS3), 47.2  s, 5.99  s}, {Eigenvalues[a] (BLAS2), 132.9  s, 52.5  s}, {Eigensystem[a] (BLAS2), 309.2  s, 140.6  s}, {SingularValues[a] (BLAS2), 201.4  s, 36.3  s}, {QRDecomposition[a] (BLAS3), 53.3  s, 8.5  s}})

New Linear Algebra Functions

A few new linear algebra functions are added in Mathematica 5.
•  generalized functions
•  MatrixRank[a] returns the rank of a matrix a
•  Norm[v, p] returns the p-norm (1 ⩽ p ⩽ ∞) of a vector v
•  Norm[a, p] returns the p-norm (1 ⩽ p ⩽ ∞) of a matrix a
•  SingularValueList[a] returns a list of nonzero singular values of a matrix a
•  SingularValueDecomposition[a] returns the full singular value decomposition of a matrix a in the form {u, s, v} where u and v are orthogonal/unitary matrices, and s is a (rectangular) diagonal matrix containing on its diagonal all singular values (including zero singular values)
The last two functions are a replacement for SingularValues.

Examples

a = (({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}));

MatrixRank[a]

{Norm[N[a], 1],  Norm[N[a],  1.5],  Norm[N[a], 2. ],  
      Norm[N[a],  2.5], Norm[N[a], ∞]}

SingularValueList[N[a]]

Examples

{u, s, v}  =  SingularValueDecomposition[N[a]];  
    Map[MatrixForm,  {u, s, v}]

{u1, s1, v1}  =  SingularValues[N[a]];  
    Map[MatrixForm,  {u1, s1, v1}]

Tr[s,  List]

Generalized Functions

There are three generalized functions added in Mathematica 5.
•  generalized eigenproblem
•  generalized Schur decomposition (QZ decomposition)
•  generalized singular value decomposition
In the version 5 they work only with machine-precision matrices. Support for arbitrary-precision matrices will be added in the future.

Generalized Eigenproblem

Given two n×n matrices A and B, the generalized (nonsymmetric) eigenproblem is defined as finding eigenvalues and corresponding nonzero eigenvectors such that

The numbers (Subscript[λ, k]) are called the generalized eigenvalues and the vectors (Subscript[x, k]) are the generalized eigenvalues. A generalized eigenvalue λ is a number satisfying the equation Det[A  - λB] == 0. The generalized eigenvalues can have an infinite value.

Mathematica 5 uses the following syntax for the generalized eigenproblem (the pair {a,b} is called a matrix pencil).

a = N[(({{1, 3, 5}, {-1, 2, 4}, {2, 3, -2}}))]; b = N[(({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}))];
Eigenvalues[{a, b}]

Note that the matrix b is singular (Det[b] == -6.661338147750939`*^-16) so you cannot use here the standard function Eigenvalues.

Similar syntax is used for the generalized singular value decomposition and generalized Schur decomposition (QZ decomposition).

SingularValueList[{a, b}]

{s, p, q, z}  =  SchurDecomposition[{a, b}];

Arbitrary-Precision Linear Algebra Functions

The changes in the version 5 are not limited to machine-precision linear algebra functions. Almost all arbitrary-precision linear algebra functions were changed. The linear algebra functions in the versions up to 4.2 were based on LINPACK and EISPACK. In the version 5 they are based on LAPACK. Only these three arbitrary-precision functions are still based on LINPACK/EISPACK: SingularValues (and their new replacements SingularValueList and SingularValueDecomposition), NullSpace and PseudoInverse. They will be replaced with the LAPACK versions in the near future. Also the generalized functions do not have arbitrary-precision counterparts.

The advantage of switching to LAPACK is evident: the linear algebra functions are more accurate and faster than the corresponding functions based on LINPACK/EISPACK.

How Fast Are They?

These timings were obtained on a Pentium 3 computer with 500 MHz CPU. The matrices a and b are real-valued random matrices of size 100×100 and precision 20. The vector v is a real-valued random vector of length 100 and precision 20.

({{         Functions, V4  .2, V5}, {Dot[a, b], 14.4  s, 7.87  s}, {Dot[a, v], 0.16  s, 0.078  s}, {Det[a], 6.12  s, 2.42  s}, {LinearSolve[a, v], 6.67  s, 2.60  s}, {Inverse[a], 20.1  s, 4.21  s}, {Eigenvalues[a], 89.4  s, 36.4  s}, {Eigensystem[a], 212.7  s, 89.0  s}, {QRDecomposition[a], 26.1  s, 8.32  s}})

Top Level BLAS and LAPACK Functions

Most BLAS functions (21) and some LAPACK functions (45) are available as top level Mathematica functions.  Although they are treated (as for now) as internal functions and are not documented, they are already used in Mathematica programs (mostly in the MatrixManipulation package).  The main difference between these top level BLAS/LAPACK functions and other Mathematica functions is in syntax. Their syntax is similar to the Fortran BLAS/LAPACK routines. They are placed in LinearAlgebra`BLAS and LinearAlgebra`LAPACK context paths.

Example 1

This is a BLAS AXPY function which do the following operations on vectors y and x (α is a scalar): y = y + αx. The vectors y and x must be initialized before calling AXPY, and y is overwritten on exit.

x  =  {1, 2};  y  =  {-2, 3};
LinearAlgebra`BLAS`AXPY[2, x, y];  y

Example 2

This is Gaussian elimination.

A = AF = (({{1`20, 2, 3}, {4, 5, 6}, {7, 8, 10}})); x = b = (({{1}, {-1}, {2}}));

First we factorize the matrix AF.

LinearAlgebra`LAPACK`GETRF[AF,  pivot,  singular];  
    singular

Then we solve the equation (the solution is returned in x).

LinearAlgebra`LAPACK`GETRS["N",  AF,  pivot,  x];  
    MatrixForm[x]

And we can check that the solution is correct:

A . x == b

Example 3

Sometimes it is useful to apply the top level BLAS/LAPACK functions in place of regular Mathematica functions. This is an example showing that when we use the LAPACK function TRTRI which computes the inverse of an upper or lower triangular matrix we can get a significant speed-up.

First we create a triangular matrix (this is not necessary for TRTRI but it is necessary for Inverse).

a  =  ainv1  = 
          Table[If[i  <=  j, 1.  + 0.1 Random[], 0. ],  {i, 1000},  {j, 
              1000}];

Timing[ainv  =  Inverse[a];]

Timing[LinearAlgebra`LAPACK`TRTRI["U", "N", ainv1,  singular]; 
      singular]

We can check that the answer is correct.

Max[Abs[ainv1 . a - IdentityMatrix[1000]]]