Singular Value Decomposition

Let A be a general real m-by-n matrix. The singular value decomposition (SVD) of A is the factorization $A=U \Sigma V^T$, where U and V are orthogonal, and $\Sigma = {\mbox {\rm diag}}( \sigma_1 , \ldots , \sigma_r )$, $r = \min (m,n)$, with $\sigma_1 \geq \cdots \geq \sigma_r \geq 0$. If A is complex, then its SVD is $A=U \Sigma V^H$ where U and V are unitary, and $\Sigma$ is as before with real diagonal elements. The $\sigma _ i$ are called the singular values, the first r columns of V the right singular vectors and the first r columns of U the left singular vectors.

The routines described in this section, and listed in Table 2.12, are used to compute this decomposition. The computation proceeds in the following stages:

1.

The matrix A is reduced to bidiagonal form: A=U₁ B V₁^T if A is real (A=U₁ B V₁^H if A is complex), where U₁ and V₁ are orthogonal (unitary if A is complex), and B is real and upper-bidiagonal when $m \geq n$ and lower bidiagonal when m < n, so that B is nonzero only on the main diagonal and either on the first superdiagonal (if $m \geq n$) or the first subdiagonal (if m<n).

2.

The SVD of the bidiagonal matrix B is computed: $B=U_2 \Sigma V_2^T$, where U₂ and V₂ are orthogonal and $\Sigma$ is diagonal as described above. The singular vectors of A are then U = U₁ U₂ and V = V₁ V₂.

The reduction to bidiagonal form is performed by the subroutine xGEBRD, or by xGBBRD for a band matrix.

The routine xGEBRD represents U₁ and V₁ in factored form as products of elementary reflectors, as described in section 5.4. If A is real, the matrices U₁ and V₁ may be computed explicitly using routine xORGBR, or multiplied by other matrices without forming U₁ and V₁ using routine xORMBR. If A is complex, one instead uses xUNGBR and xUNMBR, respectively.

If A is banded and xGBBRD is used to reduce it to bidiagonal form, U₁ and V₁ are determined as products of Givens rotations , rather than as products of elementary reflectors. If U₁ or V₁ is required, it must be formed explicitly by xGBBRD. xGBBRD uses a vectorizable algorithm, similar to that used by xSBTRD (see Kaufman [77]). xGBBRD may be much faster than xGEBRD when the bandwidth is narrow.

The SVD of the bidiagonal matrix is computed either by subroutine xBDSQR or by subroutine xBDSDC. xBDSQR uses QR iteration when singular vectors are desired [32,23], and otherwise uses the dqds algorithm [51]. xBDSQR is more accurate than its counterparts in LINPACK and EISPACK: barring underflow and overflow, it computes all the singular values of B to nearly full relative precision, independent of their magnitudes. It also computes the singular vectors much more accurately. See section 4.9 and [32,23,51] for details.

xBDSDC uses a variation of Cuppen's divide and conquer algorithm to find singular values and singular vectors [69,58]. It is much faster than xBDSQR if singular vectors of large matrices are desired. When singular values only are desired, it uses the same dqds algorithm as xBDSQR [51]. Divide-and-conquer is not guaranteed to compute singular values to nearly full relative precision, but in practice xBDSDC is often at least as accurate as xBDSQR. xBDSDC represents the singular vector matrices U₂ and V₂ in a compressed format requiring only $O(n \log n)$ space instead of n². U₂ and V₂ may subsequently be generated explicitly using routine xLASDQ, or multiplied by vectors without forming them explicitly using routine xLASD0.

If $m \gg n$, it may be more efficient to first perform a QR factorization of A, using the routine xGEQRF , and then to compute the SVD of the n-by-n matrix R, since if A = QR and $R = U \Sigma V^T$, then the SVD of A is given by $A = (QU) \Sigma V^T$. Similarly, if $m \ll n$, it may be more efficient to first perform an LQ factorization of A, using xGELQF. These preliminary QR and LQ factorizations are performed by the drivers xGESVD and xGESDD.

The SVD may be used to find a minimum norm solution to a (possibly) rank-deficient linear least squares problem (2.1). The effective rank, k, of A can be determined as the number of singular values which exceed a suitable threshold. Let $\hat{\Sigma}$ be the leading k-by-k submatrix of $\Sigma$, and $\hat{V}$ be the matrix consisting of the first k columns of V. Then the solution is given by:

$$x = \hat{V} \hat{\Sigma}^{-1} \hat{c}_1$$

where $\hat{c}_1$ consists of the first k elements of c = U^T b = U₂^T U₁^T b. U₁^T b can be computed using xORMBR, and xBDSQR has an option to multiply a vector by U₂^T.

Table 2.12: Computational routines for the singular value decomposition

Type of matrix	Operation	Single precision		Double precision
and storage scheme		real	complex	real	complex
general	bidiagonal reduction	SGEBRD	CGEBRD	DGEBRD	ZGEBRD
general band	bidiagonal reduction	SGBBRD	CGBBRD	DGBBRD	ZGBBRD
orthogonal/unitary	generate matrix after	SORGBR	CUNGBR	DORGBR	ZUNGBR
	bidiagonal reduction
	multiply matrix after	SORMBR	CUNMBR	DORMBR	ZUNMBR
	bidiagonal reduction
bidiagonal	SVD using	SBDSQR	CBDSQR	DBDSQR	ZBDSQR
	QR or dqds
	SVD using	SBDSDC		DBDSDC
	divide-and-conquer