Computes the singular value decomposition (SVD) of a rectangular matrix, A. This gives the decomposition
where V is an n×n orthogonal matrix, U is an m×m orthogonal matrix, and S is a real, rectangular diagonal matrix.
Required Arguments
A — Array of size m×n containing the matrix. (Input [/Output]) If the packaged option lin_svd_overwrite_input is used, this array is not saved on output.
S — Array of size min(m, n) containing the real singular values. These nonnegative values are in non-increasing order. (Output)
U — Array of size m×m containing the singular vectors, U. (Output)
V — Array of size n×n containing the singular vectors, V. (Output)
Optional Arguments
MROWS= m (Input) Uses array A(1:m, 1:n) for the input matrix. Default: m = size(A, 1)
NCOLS= n (Input) Uses array A(1:m, 1:n) for the input matrix. Default: n = size(A, 2)
RANK= k (Output) Number of singular values that exceed the value Small. RANK will satisfy k <= min(m, n).
iopt= iopt(:) (Input) Derived type array with the same precision as the input matrix. Used for passing optional data to the routine. The options are as follows:
Packaged Options forLIN_SVD
Option Prefix = ?
Option Name
Option Value
S_, d_, c_, z_
lin_svd_set_small
1
S_, d_, c_, z_
lin_svd_overwrite_input
2
S_, d_, c_, z_
lin_svd_scan_for_NaN
3
S_, d_, c_, z_
lin_svd_use_qr
4
S_, d_, c_, z_
lin_svd_skip_orth
5
S_, d_, c_, z_
lin_svd_use_gauss_elim
6
S_, d_, c_, z_
lin_svd_set_perf_ratio
7
iopt(IO)= ?_options(?_lin_svd_set_small, Small) If a singular value is smaller than Small, it is defined as zero for the purpose of computing the rank of A. Default: the smallest number that can be reciprocated safely
iopt(IO)= ?_options(?_lin_svd_overwrite_input, ?_dummy) Does not save the input array A(:, :).
iopt(IO)= ?_options(?_lin_svd_scan_for_NaN, ?_dummy) Examines each input array entry to find the first value such that
isNaN(a(i,j)) == .true.
See the isNaN() function, Chapter 10. Default: The array is not scanned for NaNs.
iopt(IO)= ?_options(?_lin_svd_use_qr, ?_dummy) Uses a rational QR algorithm to compute eigenvalues. Accumulate the singular vectors using this algorithm. Default: singular vectors computed using inverse iteration
iopt(IO)= ?_options(?_lin_svd_skip_Orth, ?_dummy) If the eigenvalues are computed using inverse iteration, skips the final orthogonalization of the vectors. This method results in a more efficient computation. However, the singular vectors, while a complete set, may not be orthogonal. Default: singular vectors are orthogonalized if obtained using inverse iteration
iopt(IO)= ?_options(?_lin_svd_use_gauss_elim, ?_dummy) If the eigenvalues are computed using inverse iteration, uses standard elimination with partial pivoting to solve the inverse iteration problems. Default: singular vectors computed using cyclic reduction
iopt(IO)= ?_options(?_lin_svd_set_perf_ratio, perf_ratio) Uses residuals for approximate normalized singular vectors if they have a performance index no larger than perf_ratio. Otherwise an alternate approach is taken and the singular vectors are computed again: Standard elimination is used instead of cyclic reduction, or the standard QR algorithm is used as a backup procedure to inverse iteration. Larger values of perf_ratio are less likely to cause these exceptions. Default: perf_ratio = 4
FORTRAN 90 Interface
Generic: CALL LIN_SVD (A, S, U, V[, …])
Specific: The specific interface names are S_LIN_SVD, D_LIN_SVD, C_LIN_SVD, and Z_LIN_SVD.
Description
Routine lin_svd is an implementation of the QR algorithm for computing the SVD of rectangular matrices. An orthogonal reduction of the input matrix to upper bidiagonal form is performed. Then, the SVD of a real bidiagonal matrix is calculated. The orthogonal decomposition AV = US results from products of intermediate matrix factors. See Golub and Van Loan (1989, Chapter 8) for details.
Fatal, Terminal, and Warning Error Messages
See the messages.gls file for error messages for LIN_SVD. These error messages are numbered 1001–1010; 1021–1030; 1041–1050; 1061–1070.
Examples
Example 1: Computing the SVD
The SVD of a square, random matrix A is computed. The residuals R = AV – US are small with respect to working precision. Also, see operator_ex21, supplied with the product examples.
use lin_svd_int
use rand_gen_int
implicit none
! This is Example 1 for LIN_SVD.
integer, parameter :: n=32
real(kind(1d0)), parameter :: one=1d0
real(kind(1d0)) err
real(kind(1d0)), dimension(n,n) :: A, U, V, S(n), y(n*n)
! Generate a random n by n matrix.
call rand_gen(y)
A = reshape(y,(/n,n/))
! Compute the singular value decomposition.
call lin_svd(A, S, U, V)
! Check for small residuals of the expression A*V - U*S.
Example 2: Linear Least Squares with a Quadratic Constraint
An m×n matrix equation Ax≅b, m > n, is approximated in a least-squares sense. The matrix b is size m×k. Each of the k solution vectors of the matrix x is constrained to have Euclidean length of value αj > 0. The value of αi is chosen so that the constrained solution is 0.25 the length of the nonregularized or standard least-squares equation. See Golub and Van Loan (1989, Chapter 12) for more details. In the Example 2 code, Newton’s method is used to solve for each regularizing parameter of the k systems. The solution is then computed and its length is checked. Also, see operator_ex22, supplied with the product examples.
Example 3: Generalized Singular Value Decomposition
The n×n matrices A and B are expanded in a Generalized Singular Value Decomposition (GSVD). Two n×n orthogonal matrices, U and V, and a nonsingular matrix X are computed such that
AX = Udiag(c1, .…, cn)
and
BX = Vdiag(s1, …, sn)
The values si and ci are normalized so that
The ci are nonincreasing, and the si are nondecreasing. See Golub and Van Loan (1989, Chapter 8) for more details. Our method is based on computing three SVDs as opposed to the QR decomposition and two SVDs outlined in Golub and Van Loan. As a bonus, an SVD of the matrix X is obtained, and you can use this information to answer further questions about its conditioning. This form of the decomposition assumes that the matrix
has all its singular values strictly positive. For alternate problems, where some singular values of D are zero, the GSVD becomes
UTA = diag(c1, .…, cn)W
and
VTB = diag(s1, …, sn)W
The matrix W has the same singular values as the matrix D. Also, see operator_ex23, supplied with the product examples.
Example 4: Ridge Regression as Cross-Validation with Weighting
This example illustrates a particular choice for the ridge regression problem: The least-squares problem Ax≅b is modified by the addition of a regularizing term to become
The solution to this problem, with row k deleted, is denoted by xk(λ). Using nonnegative weights (w1, …, wm), the cross-validation squared error C(λ) is given by:
With the SVD A = USVT and product g = UTb, this quantity can be written as
This expression is minimized. See Golub and Van Loan (1989, Chapter 12) for more details. In the Example 4 code, mC(λ), at p = 10 grid points are evaluated using a log-scale with respect to λ, . Array operations and intrinsics are used to evaluate the function and then to choose an approximate minimum. Following the computation of the optimum λ, the regularized solutions are computed. Also, see operator_ex24, supplied with the product examples.
