SVDC

The SVDC procedure computes the Singular Value Decomposition (SVD) of a square ( n  x  n ) or non-square ( n  x  m ) array as the product of orthogonal and diagonal arrays. SVD is a very powerful tool for the solution of linear systems, and is often used when a solution cannot be determined by other numerical algorithms.

The SVD of an ( m x n ) non-square array A is computed as the product of an ( m  x  n ) column orthogonal array U , an ( n x n ) diagonal array SV , composed of the singular values, and the transpose of an ( n x n ) orthogonal array V: A = U  SV  V T

SVDC is based on the routine ``` svdcmp``` described in section 2.6 of Numerical Recipes in C: The Art of Scientific Computing (Second Edition), published by Cambridge University Press, and is used by permission.

SVDC, A, W, U, V

Arguments

A

The square ( n x n ) or non-square ( n x m ) single- or double-precision floating-point array to decompose.

W

On output, W is an n -element output vector containing the "singular values."

U

On output, U is an n -column, m -row orthogonal array used in the decomposition of A .

V

On output, V is an n -column, n -row orthogonal array used in the decomposition of A .

Keywords

COLUMN

Set this keyword if the input array A is in column-major format (composed of column vectors) rather than in row-major format (composed of row vectors).

DOUBLE

Set this keyword to force the computation to be done in double-precision arithmetic.

Example

To find the singular values of an array A:

A = [[1.0, 2.0, -1.0, 2.5], \$ ; Define the array A.

[1.5, 3.3, -0.5, 2.0], \$

[3.1, 0.7,  2.2, 0.0], \$

[0.0, 0.3, -2.0, 5.3], \$

[2.1, 1.0,  4.3, 2.2], \$

[0.0, 5.5,  3.8, 0.2]]

SVDC, A, W, U, V ; Compute the Singular Value Decomposition.

PRINT, W ; Print the singular values.

IDL prints:

8.81973 2.65502 4.30598 6.84484

To verify the decomposition, use the relationship A = U ## SV ## TRANSPOSE(V), where SV is a diagonal array created from the output vector W.

sv = FLTARR(4, 4)

FOR K = 0, 3 DO sv(K,K) = W[K]

result = U ## sv ## TRANSPOSE(V)

PRINT, result

IDL prints:

1.00000 2.00000 -1.00000 2.50000

1.50000 3.30000 -0.500001 2.00000

3.10000 0.700000 2.20000 0.00000

2.23517e-08 0.300000 -2.00000 5.30000

2.10000 0.999999 4.30000 2.20000

-3.91155e-07 5.50000 3.80000 0.200000

This is the input array, to within machine precision.