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:

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:

Table 3-67: LAPACK Routine Basis for LA_SVD

Table 3-67: LAPACK Routine Basis for LA_SVD
Output Type	LAPACK Routine
Output Type	QR Iteration	Divide-and-conquer
Float	sgesvd	sgesdd
Double	dgesvd	dgesdd
Complex	cgesvd	cgesdd
Double complex	zgesvd	zgesdd

Syntax

Arguments

Array

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

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.

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:

STATUS = 0: The computation was successful.

STATUS > 0: The computation did not converge. The STATUS value specifies how many superdiagonals did not converge to zero.

Note
If STATUS is not specified, any error messages will output to the screen.

Examples

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

LA_SVD error:  1.0103030e-014  
LA_SVD variance:       82.802816