## Singular Value Decomposition.

This routine decomposes an*m×n*matrix

**A**, with

*m*≥

*n*, into a product of the three matrices

**U**,

**D**, and

**V**, i.e.

^{T}**A**=

**U D V**, where

^{T}**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.

**V**denotes the transpose of

^{T}**V**. If

*m*<

*n*, then the procedure may be used for the matrix

**A**. The singular values of

^{T}**A**are the diagonal elements of the diagonal matrix

**D**and correspond to the positive square roots of the eigenvalues of the matrix

**A**.

^{T}AThis 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 V**is the singular value decomposition of^{T}**A**as calculated by the routine Singular_Value_Decomposition. Given**U D V**=^{T}x**B**, then**x**=**V D**, where^{#}U^{T}B**D**is the pseudo-inverse of^{#}**D**, i.e. if*D*_{i}> 0 then*D*^{#}_{i}= 1 /*D*_{i}and if*D*_{i}= 0, then*D*^{#}_{i}= 0. Since the singular values are subject to round-off error, a user-specified tolerance is provided so that if*D*_{i}<*tolerance*, then*D*_{i}is treated as if it is 0. The default tolerance is*D*_{0}* 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 V**. where^{T}**U**,**D**,**V**constitute the singular value decomposition of**A**. Let**Astar**be the pseudo-inverse then**Astar**=**V D**, where^{#}U^{T}**D**is the pseudo-inverse of^{#}**D**, i.e. if*D*_{i}> 0 then*D*^{#}_{i}= 1 /*D*_{i}and if*D*_{i}= 0, then*D*^{#}_{i}= 0. Since the singular values are subject to round-off error, a user-specified tolerance is provided so that if*D*_{i}<*tolerance*, then*D*_{i}is treated as if it is 0. The default tolerance is*D*_{0}* 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( ).