The Singular Value Decomposition

The singular value decomposition of a matrix is represented by the by the SingularValueDecompositionT 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();

Dim aSVD = Matrix.Create(4, 4, New Double() {
    1.8, 5.25, 1.58, -1.11,
    2.88, -2.95, -2.69, -0.66,
    2.05, -0.95, -2.9, -0.59,
    -0.89, -3.8, -1.04, 0.8}, MatrixElementOrder.ColumnMajor)
Dim svd = aSVD.GetSingularValueDecomposition()
'

let aSVD = Matrix.Create(4, 4, 
            [|
                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)
let 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:

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();

Dim 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)
Dim x = svd.Solve(b)
Dim invA = svd.GetInverse()
Dim detA = svd.GetDeterminant()
'

let b = Matrix.Create(4, 2, 
         [|
            9.52;24.35;0.77;-6.22;
            18.47;2.25;-13.28;-6.21 
         |], MatrixElementOrder.ColumnMajor)
let x = svd.Solve(b)
let invA = svd.GetInverse()
let 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;

Dim S = svd.SingularValues
Dim U = svd.LeftSingularVectors
Dim V = svd.RightSingularVectors
'

let S = svd.SingularValues
let U = svd.LeftSingularVectors
let V = svd.RightSingularVectors
// //

The implementation of the singular value decomposition is based on the LAPACK routine DGESDD.

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 SparseMatrixT 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.