Systems of Linear Equations



Singular Value Decomposition

Singular Value Decomposition.

This routine decomposes an m×n matrix A, with mn, into a product of the three matrices U, D, and VT, i.e. A = U D VT, where U is an m×n matrix whose columns are mutually orthogonal, D is an n×n diagonal matrix, and V is an n×n orthogonal matrix. VT denotes the transpose of V. If m < n, then the procedure may be used for the matrix AT. The singular values of A are the diagonal elements of the diagonal matrix D and correspond to the positive square roots of the eigenvalues of the matrix ATA.

This procedure programmed here is based on the method of Golub and Reinsch as given on pages 134 - 151 of the "Handbook for Automatic Computation vol II - Linear Algebra" edited by Wilkinson and Reinsch and published by Springer-Verlag, 1971.

Function List

  • int Singular_Value_Decomposition( double* A, int nrows, int ncols, double* U, double* singular_values, double* V, double* dummy_array)

    The routine Singular_Value_Decomposition decomposes the nrows×ncols matrix A into the product of an nrows×ncols matrix U whose columns are mutually orthogonal, a diagonal matrix with diagonal elements given by the array singular_values, and an orthogonal ncols×ncols matrix V. The singular values are sorted from largest to smallest. The argument dummy_array is working storage and should be dimensioned ncols in the calling routine. The function returns a 0 if successful and a -1 if the procedure fails to converge within a set number of iterations.

  • void Singular_Value_Decomposition_Solve( double* U, double* D, double* V, double tolerance, int nrows, int ncols, double *B, double* x )

    The routine Singular_Value_Decomposition_Solve solves the system of linear equations A x = B where A = U D VT is the singular value decomposition of A as calculated by the routine Singular_Value_Decomposition. Given U D VTx = B, then
    x = V D# UTB, where D# is the pseudo-inverse of D, i.e. if Di > 0 then D#i = 1 / Di and if Di = 0, then D#i = 0. Since the singular values are subject to round-off error, a user-specified tolerance is provided so that if Di < tolerance, then Di is treated as if it is 0. The default tolerance is D0 * DBL_EPSILON * ncols, if the user-specified tolerance is less than the default tolerance, the default tolerance is used.

  • void Singular_Value_Decomposition_Inverse( double* U, double* D, double* V, double tolerance, int nrows, int ncols, double *Astar )

    The routine Singular_Value_Decomposition_Inverse calculates the pseudo-inverse of the matrix A = U D VT. where U, D, V constitute the singular value decomposition of A. Let Astar be the pseudo-inverse then Astar = V D# UT, where D# is the pseudo-inverse of D, i.e. if Di > 0 then D#i = 1 / Di and if Di = 0, then D#i = 0. Since the singular values are subject to round-off error, a user-specified tolerance is provided so that if Di < tolerance, then Di is treated as if it is 0. The default tolerance is D0 * DBL_EPSILON * ncols, if the user-specified tolerance is less than the default tolerance, the default tolerance is used.

C Source Code

  • The file, singular_value_decomposition.c, contains the routines Singular_Value_Decomposition( ), Singular_Value_Decomposition_Solve( ), and Singular_Value_Decomposition_Inverse( ).