Wolfram Research, Inc.

3.7.10 Advanced Topic: Matrix Decompositions

Singular value decomposition and pseudoinverse.

Singular value decomposition is an important element of many numerical matrix algorithms. The basic idea is to write any matrix in the form , where is a diagonal matrix, and are row orthonormal matrices, and is the Hermitian transpose of u.

The function SingularValues[m] returns a list containing the matrix , the list of diagonal elements of , and the matrix .

The diagonal elements of are known as the singular values of the matrix . One interpretation of the singular values is as follows. If you take a unit sphere in -dimensional space, and multiply each vector in it by an matrix , you will get an ellipsoid in -dimensional space. The singular values give the lengths of the principal axes of the ellipsoid. If the matrix is singular in some way, this will be reflected in the shape of the ellipsoid. In fact, the ratio of the largest singular value of a matrix to the smallest one gives a condition number of the matrix, which determines, for example, the accuracy of numerical matrix inverses.

Very small singular values are usually numerically meaningless. SingularValues removes any singular values that are smaller than a certain tolerance multiplied by the largest singular value. The option Tolerance specifies the tolerance to use.

Here is the singular value decomposition for a non-singular matrix.

In[1]:= SingularValues[N[{{1, 3}, {-4, 3}}]]

Out[1]=

Here is a convenient way to pick out the pieces of the result.

In[2]:= {u, md, v} = %

Out[2]=

This gives something close to the original matrix.

In[3]:= Transpose[u] . DiagonalMatrix[md] . v

Out[3]=

Here is the singular value decomposition of a singular matrix. Only one singular value is given.

In[4]:= SingularValues[N[{{1, 2}, {1, 2}}]]

Out[4]=

The standard definition for the inverse of a matrix fails if the matrix is not square. Using singular value decomposition, however, it is possible to define a pseudoinverse even for non-square matrices, or for singular square ones. The pseudoinverse is defined in terms of the objects , and as . The pseudoinverse has the property that the sum of the squares of all the entries in , where is an identity matrix, is minimized. The pseudoinverse found in this way is important, for example, in carrying out fits to numerical data. The pseudoinverse is sometimes known as the generalized inverse, or the Moore-Penrose inverse.

Other matrix decompositions.

Singular value decomposition writes any matrix as a product of a diagonal matrix with row and column orthonormal matrices. In some algorithms, it is also important to be able to decompose matrices as products involving triangular matrices.

QR decomposition writes any matrix as a product , where is an orthonormal matrix, denotes Hermitian transpose, and is a triangular matrix, in which all entries below the leading diagonal are zero. The function QRDecomposition[m] returns a list containing the matrices and . QR decomposition is often used in solving least-squares fitting problems, and is typically faster than singular value decomposition.

This computes the QR decomposition of a matrix, then extracts the matrices and .

In[5]:= {q, r} = QRDecomposition[ N[{{1, 5}, {-7, 3}}] ]

Out[5]=

Since is orthogonal, this gives the identity matrix.

In[6]:= Transpose[q] . q // Chop

Out[6]=

This gives the original matrix.

In[7]:= Transpose[q] . r

Out[7]=

Schur decomposition writes any matrix in the form , where is an orthogonal matrix, and is block-upper triangular. The function SchurDecomposition[m] returns a list containing the matrices and . Schur decomposition is often used in evaluating functions of matrices.

The stability of numerical matrix algorithms is sometimes improved by "pivoting" or "balancing", in which the rows and columns of the matrix are permuted and perhaps rescaled before QR or Schur decomposition is done. Setting the option Pivoting -> True tells Mathematica to perform such pivoting or balancing, and to give the matrix required as the last element of the list returned.

Jordan decomposition writes any square matrix in the form . The matrix is usually just a diagonal matrix of eigenvalues, but in general it can also have ones directly above the diagonal. The function JordanDecomposition[m] returns a list containing the matrices and . Jordan decomposition is often used in evaluating functions of exact matrices.

LU decomposition is convenient for preprocessing matrices that will appear repeatedly in collections of linear equations, as discussed in Section 3.7.8.