The Singular Value Decomposition
The Singular Value Decomposition or SVD of a matrix expresses the matrix as the product of an orthogonal (or unitary) matrix, a diagonal matrix, and an orthogonal (or unitary) matrix. The matrix can be of any type and of any shape. The SVD decomposition is usually written as
A = UΣVT,
where U and V are square, orthogonal matrices, and Σ is a diagonal matrix with the same shape as A, or
A = UΣVH,
in the complex case. The singular values of a complex matrix are real.
The singular value decomposition is the most stable of all decompositions. Aside from its own, unique, uses, it is used in situations where the highest possible accuracy is required. The numerical rank is determined using the singular value decomposition, as is the exact condition number, which is the ratio of the largest to the smallest singular value. On the downside, the method is slower than all other decompositions.
Working with Singular Value Decompositions
The singular value decomposition of a matrix is represented by the by the SingularValueDecomposition<T> class. It has no constructors. Instead, it is created by calling the GetSingularValueDecomposition method on the matrix. This method has four overloads. The first overload takes no arguments: The second overload takes a Boolean value that specifies whether the contents of the matrix may be overwritten during the SVD calculation. The default is false.
var aSVD = Matrix.Create(4, 4, new double[] {
1.80, 5.25, 1.58, -1.11,
2.88,-2.95, -2.69, -0.66,
2.05,-0.95, -2.90, -0.59,
-0.89,-3.80,-1.04, 0.80 }, MatrixElementOrder.ColumnMajor);
var svd = aSVD.GetSingularValueDecomposition();The two remaining overloads let you specify which calculations are to be performed. Because the calculation of the left and right singular vectors is relatively expensive and not always needed, it is possible to specify that only the singular values should be computed. This is done in one of two ways. The first is by passing a value of type SingularValueDecompositionFactors to the method. Alternatively, you can set the RequestedFactors property to a value of the same type before the decomposition is computed. The possible values for the SingularValueDecompositionFactors enumeration are listed below:
Value | Description |
|---|---|
SingularValues | The singular values are required. This is always included. |
LeftSingularVectors | The left singular vectors are required. |
RightSingularVectors | The right singular vectors are required. |
Thin | The singular values as well as only the corresponding left and right singular vectors are requested. |
All | The singular values as well as the left and right singular vectors are required. |
A common option is Thin, which is often used for matrices that have more rows than columns. In this case, only a few of the left singular vectors correspond to singular values. The remaining singular vectors correspond to zero rows in the matrix of singular values, and are rarely used.
The Decompose method performs the actual decomposition. This method copies the matrix if necessary. It then calls the appropriate routine to perform the actual decomposition. This method is called by other methods as needed. You will rarely need to call it explicitly.
Once the decomposition is computed, a number of operations can be performed in much less time. You can repeatedly solve a system of linear equations or a least squares problem with different right-hand sides. You can also calculate the determinant (only the magnitude) and the inverse of the base matrix:
var b = Matrix.Create(4, 2, new double[] {
9.52,24.35,0.77,-6.22,
18.47,2.25,-13.28,-6.21 }, MatrixElementOrder.ColumnMajor);
var x = svd.Solve(b);
var invA = svd.GetInverse();
var detA = svd.GetDeterminant();A number of methods are unique to the singular value decomposition. The GetPseudoInverse method returns a matrix that is the Moore-Penrose pseudo-inverse of the matrix. A pseudo-inverse is a generalization of the inverse of a matrix to matrices that are not square or singular. The Moore-Penrose pseudo-inverse is the most common, and is defined as the unique matrix A+ that satisfies the following four conditions:
A+AA+ = A+
(A+A)T = A+A
The SingularValues property returns a vector containing the singular values in decreasing order. The SingularValueMatrix returns the matrix Σ in the decomposition: a diagonal matrix of the same shape as the decomposed matrix with the singular values on the diagonal.
The LeftSingularVectors and RightSingularVectors return the matrices U and V, respectively. The singular vectors are stored in the columns of these matrices. When a thin SVD is requested, the meaning is slightly different. For a matrix with more rows than columns, the left singular matrix contains only n columns, where n is the number of columns in the decomposed matrix. For a matrix with more columns than rows, the right singular matrix (which appears transposed in the decomposition) has only m columns, where m is the number of rows in the decomposed matrix. The columns of these reduced matrices are still orthogonal.
var S = svd.SingularValues;
var U = svd.LeftSingularVectors;
var V = svd.RightSingularVectors;The implementation of the singular value decomposition is based on the LAPACK routine DGESDD.
The Singular Value Decomposition of a Sparse Matrix
In general, the singular value decomposition of a sparse matrix is not sparse. Computing the full SVD of a large, sparse matrix is therefore often all but impossible. Most calculations only require the largest singular values and their associated singular vectors. These can be computed efficiently for sparse matrices. The SparseMatrix<T> class has an extra overload of the GetSingularValueDecomposition class that returns the partial singular value decomposition of a sparse matrix. It uses an implicitly rescaled Lanczos method to compute the largest singular values based on the PROPACK package. Only real double-precision matrices are supported at this time.
The method has one additional required argument which specifies how many singular values should be returned. The method has many optional arguments that let you fine tune the calculation.