LA_SVD

The LA_SVD procedure computes the singular value decomposition (SVD) of an n-columns by m-row array as the product of orthogonal and diagonal arrays:

A is real: A = U S V^T

A is complex: A = U S V^H

where U is an orthogonal array containing the left singular vectors, S is a diagonal array containing the singular values, and V is an orthogonal array containing the right singular vectors. The superscript T represents the transpose while the superscript H represents the Hermitian, or transpose complex conjugate.

If n < m then U has dimensions (n x m), S has dimensions (n x n), and V^H has dimensions (n x n). If n ³ m then U has dimensions (m x m), S has dimensions (m x m), and V^H has dimensions (n x m). The following diagram shows the array dimensions:

LA_SVD is based on the following LAPACK routines:

For details see Anderson et al., LAPACK Users' Guide, 3rd ed., SIAM, 1999.

Syntax

LA_SVD, Array, W, U, V [, /DOUBLE] [, /DIVIDE_CONQUER] [, STATUS=variable]

Arguments

Array

The real or complex array to decompose.

W

On output, W is a vector with MIN(m, n) elements containing the singular values.

U

On output, U is an orthogonal array with MIN(m, n) columns and m rows used in the decomposition of Array. If Array is complex then U will be complex, otherwise U will be real.

V

On output, V is an orthogonal array with MIN(m, n) columns and n rows used in the decomposition of Array. If Array is complex then V will be complex, otherwise V will be real.

Keywords

DIVIDE_CONQUER

If this keyword is set, then the divide-and-conquer method is used to compute the singular vectors, otherwise, QR iteration is used. The divide-and-conquer method is faster at computing singular vectors of large matrices, but uses more memory and may produce less accurate singular values.

DOUBLE

Set this keyword to use double-precision for computations and to return a double-precision (real or complex) result. Set DOUBLE = 0 to use single-precision for computations and to return a single-precision (real or complex) result. The default is /DOUBLE if Array is double precision, otherwise the default is DOUBLE = 0.

STATUS

Set this keyword to a named variable that will contain the status of the computation. Possible values are:

Examples

Construct a sample input array A, consisting of smoothed random values:

PRO ExLA_SVD 
; Create a smoothed random array: 
n = 100 
m = 200 
seed = 12321 
a = SMOOTH(RANDOMN(seed, n, m, /DOUBLE), 5) 
 
; Compute the SVD and check reconstruction error: 
LA_SVD, a, w, u, v 
arecon = u ## DIAG_MATRIX(w) ## TRANSPOSE(v) 
PRINT, 'LA_SVD error:', MAX(ABS(arecon - a)) 
 
; Keep only the 15 largest singular values 
wfiltered = w 
wfiltered[15:*] = 0.0 
; Reconstruct the array: 
afiltered = u ## DIAG_MATRIX(wfiltered) ## TRANSPOSE(v) 
percentVar = 100*(w^2)/TOTAL(w^2) 
PRINT, 'LA_SVD Variance:', TOTAL(percentVar[0:14]) 
END 

When this program is compiled and run, IDL prints:

LA_SVD error:  1.0103030e-014 
LA_SVD variance:       82.802816 

Version History

Introduced 5.6

See Also

LA_CHOLDC, LA_LUDC, SVDC