SVD: Singular Value Decomposition

The Singular Value Decomposition (SVD) solver object can be used for computing a partial SVD of a rectangular matrix, and other related problems. It provides uniform and efficient access to several specific SVD solvers included in SLEPc, and also gives the possibility to compute the decomposition via the eigensolvers provided in the EPS package.

In many aspects, the user interface of SVD resembles that of EPS. For this reason, this chapter and chapter EPS: Eigenvalue Problem Solver have a very similar structure.

Mathematical Background

This section provides some basic concepts about the singular value decomposition and other related problems. The objective is to set up the notation and also to justify some of the solution approaches, particularly those based on the EPS object. As in the case of eigensolvers, some of the implemented methods are described in detail in the SLEPc technical reports.

For background material about the SVD, see for instance [Bai et al., 2000, ch. 6]. Many other books such as [Björck, 1996] or [Hansen, 1998] present the SVD from the perspective of its application to the solution of least squares problems and other related linear algebra problems.

The (Standard) Singular Value Decomposition (SVD)

The singular value decomposition (SVD) of an \(m\times n\) matrix \(A\) can be written as

(1)\[A=U\Sigma V^*,\]

where \(U=[u_1,\ldots,u_m]\) is an \(m\times m\) unitary matrix (\(U^*U=I\)), \(V=[v_1,\ldots,v_n]\) is an \(n\times n\) unitary matrix (\(V^*V=I\)), and \(\Sigma\) is an \(m\times n\) diagonal matrix with real diagonal entries \(\Sigma_{ii}=\sigma_i\) for \(i=1,\ldots,\min\{m,n\}\). If \(A\) is real, \(U\) and \(V\) are real and orthogonal. The vectors \(u_i\) are called the left singular vectors, the \(v_i\) are the right singular vectors, and the \(\sigma_i\) are the singular values.

Scheme of the thin SVD of a rectangular matrix \(A\).

In the following, we will assume that \(m\geq n\). If \(m<n\) then \(A\) should be replaced by \(A^*\) (note that in SLEPc this is done transparently as described later in this chapter and the user need not worry about this). In the case that \(m\geq n\), the top \(n\) rows of \(\Sigma\) contain \(\mathrm{diag}(\sigma_1,\ldots,\sigma_n)\) and its bottom \(m-n\) rows are zero. The relation equation (1) may also be written as \(AV=U\Sigma\), or

(2)\[Av_i=u_i\sigma_i\;,\quad i=1,\ldots,n,\]

and also as \(A^*U=V\Sigma^*\), or

(3)\[\begin{split}\begin{aligned} A^*u_i&=v_i\sigma_i\;,\quad i=1,\ldots,n,\\ \end{aligned}\end{split}\]

(4)\[\begin{aligned} A^*u_i&=0\;,\quad i=n+1,\ldots,m. \end{aligned}\]

The last left singular vectors corresponding to equation (4) are often not computed, especially if \(m\gg n\). In that case, the resulting factorization is sometimes called the thin SVD, \(A=U_n\Sigma_n V_n^*\), and is depicted in figure Scheme of the thin SVD of a rectangular matrix A.. This factorization can also be written as

(5)\[A=\sum_{i=1}^{n}\sigma_iu_iv_i^*.\]

Each \((\sigma_i,u_i,v_i)\) is called a singular triplet.

The singular values are real and nonnegative, \(\sigma_1\geq\sigma_2\geq\ldots\geq\sigma_r>\sigma_{r+1}=\ldots=\sigma_n=0\), where \(r=\mathrm{rank}(A)\). It can be shown that \(\{u_1,\ldots,u_r\}\) span the range of \(A\), \(\mathcal{R}(A)\), whereas \(\{v_{r+1},\ldots,v_n\}\) span the null space of \(A\), \(\mathcal{N}(A)\).

If the zero singular values are dropped from the sum in equation (5), the resulting factorization, \(A=\sum_{i=1}^{r}\sigma_iu_iv_i^*\), is called the compact SVD, \(A=U_r\Sigma_r V_r^*\).

In the case of a very large and sparse \(A\), it is usual to compute only a subset of \(k\leq r\) singular triplets. We will refer to this decomposition as the truncated SVD of \(A\). It can be shown that the matrix \(A_k=U_k\Sigma_k V_k^*\) is the best rank-\(k\) approximation to matrix \(A\), in the least squares sense.

In general, one can take an arbitrary subset of the summands in equation (5), and the resulting factorization is called the partial SVD of \(A\). As described later in this chapter, SLEPc allows the computation of a partial SVD corresponding to either the \(k\) largest or smallest singular triplets.

Equivalent Eigenvalue Problems

It is possible to formulate the problem of computing the singular triplets of a matrix \(A\) as an eigenvalue problem involving a Hermitian matrix related to \(A\). There are two possible ways of achieving this:

1. With the cross product matrix, either \(A^*A\) or \(AA^*\).

2. With the cyclic matrix, \(H(A)=\left[\begin{smallmatrix}0&A\\A^*&0\end{smallmatrix}\right]\).

In SLEPc, the computation of the SVD is usually based on one of these two alternatives, either by passing one of these matrices to an EPS object or by performing the computation implicitly.

By pre-multiplying equation (2) by \(A^*\) and then using equation (3), the following relation results

(6)\[A^*Av_i=\sigma_i^2v_i,\]

that is, the \(v_i\) are the eigenvectors of matrix \(A^*A\) with corresponding eigenvalues equal to \(\sigma_i^2\). Note that after computing \(v_i\) the corresponding left singular vector, \(u_i\), is readily available through equation (2) with just a matrix-vector product, \(u_i=\frac{1}{\sigma_i}Av_i\).

Alternatively, one could first compute the left vectors and then the right ones. For this, pre-multiply equation (3) by \(A\) and then use equation (2) to get

(7)\[AA^*u_i=\sigma_i^2u_i.\]

In this case, the right singular vectors are obtained as \(v_i=\frac{1}{\sigma_i}A^*u_i\).

The two approaches represented in equations (6) and (7) are very similar. Note however that \(A^*A\) is a square matrix of order \(n\) whereas \(AA^*\) is of order \(m\). In cases where \(m\gg n\), the computational effort will favor the \(A^*A\) approach. On the other hand, the eigenproblem equation (6) has \(n-r\) zero eigenvalues and the eigenproblem equation (7) has \(m-r\) zero eigenvalues. Therefore, continuing with the assumption that \(m\geq n\), even in the full rank case the \(AA^*\) approach may have a large null space resulting in difficulties if the smallest singular values are sought. In SLEPc, this will be referred to as the cross product approach and will use whichever matrix is smaller, either \(A^*A\) or \(AA^*\).

Computing the SVD via the cross product approach may be adequate for determining the largest singular triplets of \(A\), but the loss of accuracy can be severe for the smallest singular triplets. The cyclic matrix approach is an alternative that avoids this problem, but at the expense of significantly increasing the cost of the computation. Consider the eigendecomposition of

(8)\[\begin{split}H(A)=\begin{bmatrix}0&A\\A^*&0\end{bmatrix},\end{split}\]

which is a Hermitian matrix of order \((m+n)\). It can be shown that \(\pm\sigma_i\) is a pair of eigenvalues of \(H(A)\) for \(i=1,\ldots,r\) and the other \(m+n-2r\) eigenvalues are zero. The unit eigenvectors associated with \(\pm\sigma_i\) are \(\frac{1}{\sqrt{2}}\left[\begin{smallmatrix}\pm u_i\\v_i\end{smallmatrix}\right]\). Thus it is possible to extract the singular values and the left and right singular vectors of \(A\) directly from the eigenvalues and eigenvectors of \(H(A)\). Note that in this case the singular values are not squared, and therefore the computed values will be more accurate (especially the small ones). The drawback in this case is that small eigenvalues are located in the interior of the spectrum.

The Generalized Singular Value Decomposition (GSVD)

An extension of the SVD to the case of two matrices is the generalized singular value decomposition (GSVD), which can be applied in constrained least squares problems, among others. An overview of the problem can be found in [Golub and van Loan, 1996, 8.7.3].

Consider two matrices, \(A\in\mathbb{C}^{m\times n}\) with \(m\geq n\) and \(B\in\mathbb{C}^{p\times n}\). Note that both matrices must have the same column dimension. Then there exist two unitary matrices \(U\in\mathbb{C}^{m\times m}\) and \(V\in\mathbb{C}^{p\times p}\) and an invertible matrix \(X\in\mathbb{C}^{n\times n}\) such that

(9)\[U^*AX=C,\qquad V^*BX=S,\]

where \(C=\mathrm{diag}(c_1,\dots,c_n)\) and \(S=\mathrm{diag}(s_{n-q+1},\dots,s_n)\) with \(q=\min(p,n)\). The values \(c_i\) and \(s_i\) are real and nonnegative, and the ratios define the generalized singular values,

(10)\[\sigma(A,B)\equiv\{c_1/s_1,\dots,c_q/s_q\},\]

and if \(p<n\) we can consider that the first \(n-p\) generalized singular values are infinite, as if \(s_1=\dots=s_{n-p}=0\). Note that if \(B\) is the identity matrix, \(X\) can be taken to be unitary and then we recover the standard SVD, \(\sigma(A,I)=\sigma(A)\), that is why equation (9) is considered a generalization of the SVD.

The diagonal matrices \(C\) and \(S\) satisfy \(C^*C+S^*S=I\), and are related to the CS decomposition [Golub and van Loan, 1996, 2.6.4] associated with the orthogonal factor of the QR factorization of matrices \(A\) and \(B\) stacked, that is, if

(11)\[\begin{split}Z:=\begin{bmatrix}A\\B\end{bmatrix}=\begin{bmatrix}Q_1\\Q_2\end{bmatrix}R,\end{split}\]

then \(C\) and \(S\) can be obtained from the singular values of \(Q_1\) and \(Q_2\), respectively. The matrix \(Z\) is relevant for algorithms and is often built explicitly.

Scheme of the thin GSVD of two matrices \(A\) and \(B\), for the case \(p < n\).

The above description assumes that matrix \(Z\) has full column rank, and the rank of \(B\) is also \(q\). In the general case where these assumptions do not hold, the structure of matrices \(C\) and \(S\) is a bit more complicated. This includes also the case where \(m<n\). A detailed description of those cases can be found in [Anderson et al., 1999, 2.3.5.3].

As in the case of the SVD, one can consider thin, compact, truncated, and partial versions of the GSVD. A pictorial representation of the thin GSVD is shown in figure Scheme of the thin GSVD of two matrices A and B, for the case p < n..

The columns of \(X\), \(x_i\), are called the (right) generalized singular vectors. The left vectors in this case would correspond to the columns of \(U\) and \(V\). In SLEPc, the user interface will provide these vectors stacked on top of each other, as a single \((m+p)\)-vector \(\begin{bmatrix}u_i\\v_i\end{bmatrix}\).

Equivalent Eigenvalue Problems

In the GSVD it is also possible to formulate the problem as an eigenvalue problem, which opens the door to approach its solution via EPS. The columns of \(X\) satisfy

(12)\[s_i^2A^*Ax_i=c_i^2B^*Bx_i,\]

and so if \(s_i\neq 0\) then \(A^*Ax_i=\sigma_i^2B^*Bx_i\), a generalized eigenvalue problem for the matrix pair \((A^*A,B^*B)\). This is the analog of the cross product matrix eigenproblem of equation (6).

The formulation that is analog to the eigenproblem associated with the cyclic matrix equation (8) is to solve the generalized eigenvalue problem defined by any of the matrix pairs

(13)\[\begin{split}\left( \begin{bmatrix}0&A\\A^*&0\end{bmatrix}, \begin{bmatrix}I&0\\0&B^*B\end{bmatrix} \right), \qquad\text{or}\qquad \left( \begin{bmatrix}0&B\\B^*&0\end{bmatrix}, \begin{bmatrix}I&0\\0&A^*A\end{bmatrix} \right).\end{split}\]

The Hyperbolic Singular Value Decomposition (HSVD)

The hyperbolic singular value decomposition (HSVD) was introduced in [Onn et al., 1991], motivated by some signal processing applications such as the so-called covariance differencing problem. The formulation of the HSVD is similar to that of the SVD, except that \(U\) is orthogonal with respect to a signature matrix,

(14)\[A=U\Sigma V^*,\qquad U^*\Omega U=\tilde\Omega,\]

where \(\Omega=\mathrm{diag}(\pm 1)\) is an \(m\times m\) signature matrix provided by the user, while \(\tilde\Omega\) is another signature matrix obtained as part of the solution. Sometimes \(U\) is said to be a hyperexchange matrix, or also an \((\Omega,\tilde\Omega)\)-orthogonal matrix. Note that in the problem definition normally found in the literature it is \(V\) that is \((\Omega,\tilde\Omega)\)-orthogonal and not \(U\). We choose this definition for consistency with respect to the generalized HSVD of two matrices. If the user wants to compute the HSVD according to the alternative definition, then it suffices to (conjugate) transpose the input matrix \(A\), using for instance MatHermitianTranspose.

As in the case of the SVD, the solution of the problem consists in singular triplets \((\sigma_i,u_i,v_i)\), with \(\sigma_i\) real and nonnegative and sorted in nonincreasing order. Note that these quantities are different from those of section The (Standard) Singular Value Decomposition (SVD), even though we use the same notation here. With each singular triplet, there is an associated sign \(\tilde\omega_i\) (either 1 or \(-1\)), the corresponding diagonal element of \(\tilde\Omega\). In SLEPc, this value is not returned by the user interface, but if required it can be easily computed as \(u_i^*\Omega u_i\).

The relations between left and right singular vectors are slightly different from those of the standard SVD. We have \(AV=U\Sigma\) and \(A^*\Omega U=V\Sigma^*\tilde\Omega\), so for \(m\geq n\) the following relations hold, together with equation (2):

(15)\[\begin{split}\begin{aligned} A^*\Omega u_i&=v_i\sigma_i\tilde\omega_i\;,\quad i=1,\ldots,n,\\ \end{aligned}\end{split}\]

(16)\[\begin{aligned} A^*\Omega u_i&=0\;,\quad i=n+1,\ldots,m. \end{aligned}\]

In SLEPc we will compute a partial HSVD consisting of either the largest or smallest hyperbolic singular triplets. Note that the sign \(\tilde\omega_i\) is not used when sorting for largest or smallest \(\sigma_i\).

Equivalent Eigenvalue Problems

Once again, we can derive cross and cyclic schemes to compute the decomposition by solving an eigenvalue problem. The cross product matrix approach has two forms, to be selected depending on whether \(m\geq n\) or not, as discussed in section The (Standard) Singular Value Decomposition (SVD). The first form,

(17)\[A^*\Omega Av_i=\sigma_i^2\tilde\omega_iv_i,\]

is derived by pre-multiplying equation (2) by \(A^*\Omega\) and then using equation (15). This eigenproblem can be solved as a HEP (cf. section Defining the Problem) and may have both positive and negative eigenvalues, corresponding to \(\tilde\omega_i=1\) and \(\tilde\omega_i=-1\), respectively. Once the right vector \(v_i\) has been computed, the corresponding left vector can be obtained using equation (2) with just a matrix-vector product, \(u_i=\sigma_i^{-1}Av_i\).

The second form of cross computes the left vectors first, by pre-multiplying equation (15) by \(A\) and then using equation (2),

(18)\[AA^*\Omega u_i=\sigma_i^2\tilde\omega_iu_i.\]

In this case, the right singular vectors are obtained as \(v_i=(\sigma_i\tilde\omega_i)^{-1}A^*\Omega u_i\). The coefficient matrix of (18) is non-Hermitian, so the eigenproblem has to be solved as non-Hermitian, or alternatively it can be formulated as a generalized eigenvalue problem of GHIEP type (cf. section Defining the Problem) for the indefinite pencil \((AA^*,\Omega)\),

(19)\[AA^*\hat{u}_i=\sigma_i^2\tilde\omega_i\Omega \hat{u}_i,\]

with \(\hat{u}_i=\Omega u_i\). The eigenvectors obtained from equation (18) or equation (19) must be normalized so that \(U^*\Omega U=\tilde\Omega\) holds.

Finally, in the cyclic matrix approach for the HSVD we must solve a generalized eigenvalue problem defined by the matrices of order \((m+n)\)

(20)\[\begin{split}H(A)=\begin{bmatrix}0&A\\A^*&0\end{bmatrix},\qquad \hat\Omega=\begin{bmatrix}\Omega&0\\0&I\end{bmatrix}.\end{split}\]

As in the case of equation (19), this pencil is Hermitian-indefinite and hence it may have complex eigenvalues. However, it can be shown that nonzero eigenvalues of the pencil \((H(A),\hat\Omega)\) are either real (equal to \(\pm\sigma_i\) for a certain hyperbolic singular value \(\sigma_i\)) or purely imaginary (equal to \(\pm\sigma_ij\) for a certain \(\sigma_i\) with \(j=\sqrt{-1}\)). The associated eigenvectors are \(\left[\begin{smallmatrix}\varsigma_i u_i\\v_i\end{smallmatrix}\right]\), where \(\varsigma_i\) is either \(\pm 1\) or \(\pm j\).

Basic Usage

From the perspective of the user interface, the SVD package is very similar to EPS, with some differences that will be highlighted shortly.

Example code for basic solution with SVD.

SVD       svd;       /*  SVD solver context  */
Mat       A;         /*  problem matrix      */
Vec       u, v;      /*  singular vectors    */
PetscReal sigma;     /*  singular value      */
PetscInt  j, nconv;
PetscReal error;

SVDCreate( PETSC_COMM_WORLD, &svd );
SVDSetOperators( svd, A, NULL );
SVDSetProblemType( svd, SVD_STANDARD );
SVDSetFromOptions( svd );
SVDSolve( svd );
SVDGetConverged( svd, &nconv );
for (j=0; j<nconv; j++) {
  SVDGetSingularTriplet( svd, j, &sigma, u, v );
  SVDComputeError( svd, j, SVD_ERROR_RELATIVE, &error );
}
SVDDestroy( &svd );

The basic steps for computing a partial SVD with SLEPc are illustrated in figure Example code for basic solution with SVD.. The steps are more or less the same as those described in chapter EPS: Eigenvalue Problem Solver for the eigenvalue problem. First, the solver context is created with SVDCreate. Then the problem matrices have to be specified with SVDSetOperators and the type of problem can be selected via SVDSetProblemType. Then, a call to SVDSolve invokes the actual solver. After that, SVDGetConverged is used to determine how many solutions have been computed, which are retrieved with SVDGetSingularTriplet. Finally, SVDDestroy gets rid of the object.

If one compares this example code with the EPS example in figure Example code for basic solution with EPS, the most outstanding differences are the following:

• The singular value is a PetscReal, not a PetscScalar.

• Each singular vector is defined with a single Vec object, not two as was the case for eigenvectors.

Defining the Problem

Defining the problem consists in specifying the problem matrix \(A\) (and also a second matrix \(B\) in case of the GSVD), the problem type, and the portion of the spectrum to be computed.

The problem matrices are provided with the following function SVDSetOperators

SVDSetOperators(SVD svd,Mat A,Mat B);

where A can be any matrix, not necessarily square, stored in any allowed PETSc format including the matrix-free mechanism (see section Supported Matrix Types for a detailed discussion), and B is generally set to NULL except in the case of computing the GSVD. If the problem is hyperbolic (section The Hyperbolic Singular Value Decomposition (HSVD)) then in addition a signature matrix \(\Omega\) must be provided with SVDSetSignature

SVDSetSignature(SVD svd,Vec omega);

Note that \(\Omega\) is represented as a vector containing values equal to 1 or \(-1\).

The problem type can be specified with the function SVDSetProblemType

SVDSetProblemType(SVD svd,SVDProblemType type);

Note that in SVD calling it is currently not strictly required, since the problem type can be deduced from the information passed by the user: if only one matrix is passed to SVDSetOperators the problem type will default to a standard SVD (or hyperbolic SVD if a signature has been provided), while if two matrices are passed then it defaults to GSVD. Still, it is recommended to always call SVDSetProblemType. Table Problem types considered in SVD. lists the supported problem types.

Problem types considered in SVD.

Problem Type	SVDProblemType	Command line key
Standard SVD	SVD_STANDARD	-svd_standard
Generalized SVD (GSVD)	SVD_GENERALIZED	-svd_generalized
Hyperbolic SVD (HSVD)	SVD_HYPERBOLIC	-svd_hyperbolic

It is important to note that all SVD solvers in SLEPc make use of both \(A\) and \(A^*\), as suggested by the description in section The (Standard) Singular Value Decomposition (SVD). \(A^*\) is not explicitly passed as an argument to SVDSetOperators, therefore it will have to stem from \(A\). There are two possibilities for this: either \(A\) is transposed explicitly and \(A^*\) is created as a distinct matrix, or \(A^*\) is handled implicitly via MatMultTranspose (or MatMultHermitianTranspose in the complex case) operations whenever a matrix-vector product is required in the algorithm. The default is to build \(A^*\) explicitly, but this behavior can be changed with SVDSetImplicitTranspose

SVDSetImplicitTranspose(SVD svd,PetscBool impl);

In section The (Standard) Singular Value Decomposition (SVD), it was mentioned that in SLEPc the cross product approach chooses the smallest of the two possible cases \(A^*A\) or \(AA^*\), that is, \(A^*A\) is used if \(A\) is a tall, thin matrix (\(m\geq n\)), and \(AA^*\) is used if \(A\) is a fat, short matrix (\(m<n\)). In fact, what SLEPc does internally is that if \(m<n\) the roles of \(A\) and \(A^*\) are reversed. This is equivalent to transposing all the SVD factorization, so left singular vectors become right singular vectors and vice versa. This is actually done in all singular value solvers, not only the cross product approach. The objective is to simplify the number of cases to be treated internally by SLEPc, as well as to reduce the computational cost in some situations. Note that this is done transparently and the user need not worry about transposing the matrix, only to indicate how the transpose has to be handled, as explained above.

The user can specify how many singular values and vectors to compute. The default is to compute only one singular triplet. The function SVDSetDimensions

SVDSetDimensions(SVD svd,PetscInt nsv,PetscInt ncv,PetscInt mpd);

allows the specification of the number of singular values to compute, nsv. The next argument can be set to prescribe the number of column vectors to be used by the solution algorithm, ncv, that is, the largest dimension of the working subspace. These two parameters can also be set at run time with the options -svd_nsv and -svd_ncv. For example, the command line

$ ./program -svd_nsv 10 -svd_ncv 24

requests 10 singular values and instructs to use 24 column vectors. Note that ncv must be at least equal to nsv, although in general it is recommended (depending on the method) to work with a larger subspace, for instance \(\mathtt{ncv}\geq2\cdot\mathtt{nsv}\) or even more. As in the case of the EPS object, the last argument, mpd, can be used to limit the maximum dimension of the projected problem, as discussed in section Computing a Large Portion of the Spectrum. Using this parameter is especially important in the case that a large number of singular values are requested.

Available possibilities for selection of the singular values of interest.

SVDWhich	Command line key	Sorting criterion
SVD_LARGEST	-svd_largest	Largest \(\sigma\)
SVD_SMALLEST	-svd_smallest	Smallest \(\sigma\)

For the selection of the portion of the spectrum of interest, there are only two possibilities in the case of SVD: largest and smallest singular values, see table Available possibilities for selection of the singular values of interest.. The default is to compute the largest ones, but this can be changed with SVDSetWhichSingularTriplets

SVDSetWhichSingularTriplets(SVD svd,SVDWhich which);

which can also be specified at the command line. This criterion is used both for configuring how the algorithm seeks singular values and also for sorting the computed values. In contrast to the case of EPS, computing singular values located in the interior part of the spectrum is difficult, the only possibility is to use an EPS object combined with a spectral transformation (this possibility is explained in detail in the next section). Note that in this case, the value of which applies to the transformed spectrum.

Using a threshold to specify wanted singular values

In some applications, the number of wanted singular values is not known a priori. For instance, suppose that matrix \(A\) is known to have low rank, and we want to approximate it by a truncated SVD containing the leading singular values. The goal is to have a good approximation, that is, discard only small singular values. The numerical rank \(k\) might be difficult to estimate, due to discretization error or other reasons. But since we assume that singular values decay abruptly around some unknown value \(k\), we can configure SLEPc to detect this decay.

The threshold \(\delta\) can be set with SVDSetThreshold

SVDSetThreshold(SVD svd,PetscReal thres,PetscBool rel);

When rel is PETSC_TRUE, the solver will compute \(k\) singular values, with \(\sigma_j\geq\delta\cdot\sigma_1\), for \(j=1,\dots,k-1\), where \(k\) is computed internally. In this case, the threshold is relative to the largest singular value \(\sigma_1\), i.e., the matrix norm. It can be interpreted as a percentage, for instance \(\delta=0.8\) means compute singular values that are at least 80% of the matrix norm. Alternatively, an absolute threshold can be used by setting rel to PETSC_FALSE

Illustration of threshold usage with the singular values of the rdb200 matrix. The red line represents a threshold of \(\delta=0.8\).

In the low-rank example suggested above, the singular values will decay fast for a relatively small \(k\), so a relative threshold \(\delta=0.5\) or less could do the job. But care must be taken in matrices whose singular values decay progressively, as in the example of figure Illustration of threshold usage with the singular values of the rdb200 matrix, since a small \(\delta\) would imply computing many singular triplets and hence a very high cost, both computationally and in memory. Since the number of computed singular values is not known a priori, the solver will need to reallocate the basis of vectors internally, to have enough room to accommodate all the singular vectors, so this option must be used with caution to avoid out-of-memory problems. The recommendation is to set the value of ncv to be larger than the estimated number of singular values, to minimize the number of reallocations. You can also use the nsv parameter in combination with the threshold to stop in case the number of computed singular triplets exceeds that value.

An absolute threshold is also available when computing smallest singular values, in which case the solver will compute the \(k\) smallest singular values, where \(\sigma_j<\delta\), \(j=n-k+1,\dots,n\).

Selecting the SVD Solver

The available methods for computing the partial SVD are shown in table List of solvers available in the SVD module. In the column of supported problems, ‘G’ means GSVD and ‘H’ means HSVD (the standard SVD is supported by all solvers).. These methods can be classified in the following categories:

• Solvers based on EPS. These solvers set up an EPS object internally, thus using the available eigensolvers for solving the SVD problem. The two possible approaches in this case are the cross product matrix and the cyclic matrix, as described in section Mathematical Background.

• Specific SVD solvers. These are typically eigensolvers that have been adapted algorithmically to exploit the structure of the SVD or GSVD problems. There are two Lanczos-type solvers in this category: Lanczos and thick-restart Lanczos, see Hernandez et al. [2007] for a detailed description of these methods. In this category, we could also add the randomized SVD (RSVD), a special solver that does not compute individual singular vectors accurately, but rather a low-rank approximation of \(A\) by means of randomized techniques.

• The LAPACK solver. This is an interface to some LAPACK routines, analog of those in the case of eigenproblems. These routines operate in dense mode with only one processor and therefore are suitable only for moderate size problems. This solver should be used only for debugging purposes.

• External packages: ScaLAPACK and ELEMENTAL are dense packages and compute the complete SVD, while PRIMME offers Davidson-type methods to compute only a few singular triplets.

The default solver is the one that uses the cross product matrix (cross), usually the fastest and most memory-efficient approach, at least for the standard SVD. See a more detailed explanation below.

List of solvers available in the SVD module. In the column of supported problems, ‘G’ means GSVD and ‘H’ means HSVD (the standard SVD is supported by all solvers).

Method	SVDType	Options Database Name	Supported Problems	Default
Cross Product	SVDCROSS	cross	G,H	\(\star\)
Cyclic Matrix	SVDCYCLIC	cyclic	G,H
Lanczos	SVDLANCZOS	lanczos
Thick-restart Lanczos	SVDTRLANCZOS	trlanczos	G,H
Randomized SVD (RSVD)	SVDRANDOMIZED	randomized
LAPACK solver	SVDLAPACK	lapack	G
Wrapper to ScaLAPACK	SVDSCALAPACK	scalapack
Wrapper to KSVD	SVDKSVD	ksvd
Wrapper to ELEMENTAL	SVDELEMENTAL	elemental
Wrapper to PRIMME	SVDPRIMME	primme

The solution method can be specified procedurally or via the command line. The application programmer can set it by means of the command SVDSetType

SVDSetType(SVD svd,SVDType method);

while the user writes the options database command -svd_type followed by the name of the method (see table List of solvers available in the SVD module. In the column of supported problems, ‘G’ means GSVD and ‘H’ means HSVD (the standard SVD is supported by all solvers).).

Note that not all solvers in the SVD module support non-standard problems (sections The Generalized Singular Value Decomposition (GSVD)-The Hyperbolic Singular Value Decomposition (HSVD)). Table List of solvers available in the SVD module. In the column of supported problems, ‘G’ means GSVD and ‘H’ means HSVD (the standard SVD is supported by all solvers). indicates which solvers can be used for these.

The EPS-based solvers deserve some additional comments. These SVD solvers work by creating an EPS object internally and setting up an eigenproblem of type EPS_HEP (or EPS_GHEP in the case of the GSVD, or EPS_GHIEP in the case of the HSVD). These solvers implement the cross product matrix and the cyclic matrix approaches as described in section Mathematical Background. Therefore, the operator matrix associated with the EPS object will be \(A^*A\) in the case of the cross solver and \(H(A)\) in the case of the cyclic solver (with variations in the case of non-standard problems).

In the case of the cross solver, the matrix \(A^*A\) is not built explicitly by default, since the resulting matrix may be much less sparse than the original matrix. By default, a shell matrix is created internally in the SVD object and passed to the EPS object. Still, the user may choose to force the computation of \(A^*A\) explicitly, by means of PETSc’s sparse matrix-matrix product subroutines. This is set with SVDCrossSetExplicitMatrix

SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicit);

In the cyclic solver the user can also choose between handling \(H(A)\) implicitly as a shell matrix (the default), or forming it explicitly, that is, storing its elements in a distinct matrix. The function for setting this option is SVDCyclicSetExplicitMatrix

SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicit);

The EPS object associated with the cross and cyclic SVD solvers is created with a set of reasonable default parameters. However, it may sometimes be necessary to change some of the EPS options such as the eigensolver. To allow application programmers to set any of the EPS options directly within the code, the following routines are provided to extract the EPS context from the SVD object, SVDCrossGetEPS SVDCyclicGetEPS

SVDCrossGetEPS(SVD svd,EPS *eps);
SVDCyclicGetEPS(SVD svd,EPS *eps);

A more convenient way of changing EPS options is through the command-line. This is achieved simply by prefixing the EPS options with -svd_cross_ or -svd_cyclic_ as in the following example:

$ ./program -svd_type cross -svd_cross_eps_type gd

At this point, one may consider changing also the options of the ST object associated with the EPS object in cross and cyclic SVD solvers, for example to compute singular values located at the interior of the spectrum via a shift-and-invert transformation. This is indeed possible, but some considerations have to be taken into account. When \(A^*A\) or \(H(A)\) are managed as shell matrices, then the potential of the spectral transformation is limited seriously, because some of the required operations will not be defined if working with implicit matrices (this is discussed briefly in sections Supported Matrix Types and Explicit Computation of Coefficient Matrix). The following example illustrates the computation of interior singular values with the cyclic solver with explicit \(H(A)\) matrix:

$ ./program -svd_type cyclic -svd_cyclic_explicitmatrix
            -svd_cyclic_st_type sinvert -svd_cyclic_eps_target 12.0
            -svd_cyclic_st_ksp_type preonly -svd_cyclic_st_pc_type lu

Similarly, in the case of GSVD the thick-restart Lanczos solver uses a KSP solver internally, that can be configured by accessing it with SVDTRLanczosGetKSP

SVDTRLanczosGetKSP(SVD svd,KSP *ksp);

or with the corresponding command-line options prefixed with -svd_trlanczos_. This KSP object is used to solve a linear least squares problem at each Lanczos step with coefficient matrix \(Z\) equation (11), which by default is a shell matrix but the user can choose to create it explicitly with the function SVDTRLanczosSetExplicitMatrix.

Retrieving the Solution

Once the call to SVDSolve is complete, all the data associated with the computed partial SVD is kept internally in the SVD object. This information can be obtained by the calling program by means of a set of functions described below.

As in the case of eigenproblems, the number of computed singular triplets depends on the convergence and, therefore, it may be different from the number of solutions requested by the user. So the first task is to find out how many solutions are available, with SVDGetConverged

SVDGetConverged(SVD svd,PetscInt *nconv);

Usually, the number of converged solutions, nconv, will be equal to nsv, but in general it can be a number ranging from 0 to ncv (here, nsv and ncv are the arguments of function SVDSetDimensions).

Normally, the user is interested in the singular values only, or the complete singular triplets. The function SVDGetSingularTriplet

SVDGetSingularTriplet(SVD svd,PetscInt j,PetscReal *sigma,Vec u,Vec v);

returns the \(j\)-th computed singular triplet, \((\sigma_j,u_j,v_j)\), where both \(u_j\) and \(v_j\) are normalized to have unit norm[1]. Typically, this function is called inside a loop for each value of j from 0 to nconv–1. Note that singular values are ordered according to the same criterion specified with function SVDSetWhichSingularTriplets for selecting the portion of the spectrum of interest.

In some applications, it may be enough to compute only the right singular vectors. This is especially important in cases in which memory requirements are critical (remember that both \(U_k\) and \(V_k\) are dense matrices, and \(U_k\) may require much more storage than \(V_k\), see figure Scheme of the thin SVD of a rectangular matrix A.). In SLEPc, there is no general option for specifying this, but the default behavior of some solvers is to compute only right vectors and allocate/compute left vectors only in the case that the user requests them. This is done in the cross solver and in some special variants of other solvers such as one-sided Lanczos (consult the Hernandez et al. [2007] technical report for specific solver options).

In the case of the GSVD, the sigma argument of SVDGetSingularTriplet contains \(\sigma_i=c_i/s_i\) and the second Vec argument (v) contains the right singular vectors (\(x_i\)), while the first Vec argument (u) contains the other vectors of the decomposition stacked on top of each other, as a single \((m+p)\)-vector: \(\left[\begin{smallmatrix}u_i\\v_i\end{smallmatrix}\right]\).

Reliability of the Computed Solution

In SVD computations, a-posteriori error bounds are much the same as in the case of Hermitian eigenproblems, due to the equivalence discussed in section The (Standard) Singular Value Decomposition (SVD). The residual vector is defined in terms of the cyclic matrix, \(H(A)\), so its norm is

(21)\[\|r_\mathrm{SVD}\|_2=\left(\|A\tilde{v}-\tilde{\sigma}\tilde{u}\|_2^2+\|A^*\tilde{u}-\tilde{\sigma}\tilde{v}\|_2^2\right)^{\frac{1}{2}},\]

where \(\tilde{\sigma}\), \(\tilde{u}\) and \(\tilde{v}\) represent any of the nconv computed singular triplets delivered by SVDGetSingularTriplet.

Given the above definition, the following relation holds

(22)\[|\sigma-\tilde{\sigma}|\leq \|r_\mathrm{SVD}\|_2,\]

where \(\sigma\) is an exact singular value. The associated error can be obtained in terms of \(\|r_\mathrm{SVD}\|_2\) with the following function: SVDComputeError

SVDComputeError(SVD svd,PetscInt j,SVDErrorType type,PetscReal *error);

In the case of the GSVD, the function SVDComputeError will compute a residual norm based on the two relations equation (9),

(23)\[\|r_\mathrm{GSVD}\|_2=\left(\|\tilde{s}^2A^*\tilde{u}-\tilde{c}B^*B\tilde{x}\|_2^2+\|\tilde{c}^2B^*\tilde{v}-\tilde{s}A^*A\tilde{x}\|_2^2\right)^{\frac{1}{2}},\]

where \(\tilde{x}\), \(\tilde{u}\), \(\tilde{v}\) are the computed singular vectors corresponding to \(\tilde{\sigma}\), and \(\tilde{c}\), \(\tilde{s}\) are obtained from \(\tilde{\sigma}\) as \(\tilde{s}=1/\sqrt{1+\tilde{\sigma}^2}\) and \(\tilde{c}=\tilde{\sigma}\tilde{s}\). See [Alvarruiz et al., 2024] for details.

Similarly, in the HSVD we employ a modified residual

(24)\[\|r_\mathrm{HSVD}\|_2=\left(\|A\tilde{v}-\tilde{\sigma}\tilde{u}\|_2^2+\|A^*\Omega\tilde{u}-\tilde{\sigma}\tilde{\omega}\tilde{v}\|_2^2\right)^{\frac{1}{2}},\]

where \(\tilde\omega\) is the corresponding element of the signature \(\tilde\Omega\) of the definition equation (14).

Controlling and Monitoring Convergence

Similarly to the case of eigensolvers, in SVD the number of iterations carried out by the solver can be determined with SVDGetIterationNumber, and the tolerance and maximum number of iterations can be set with SVDSetTolerances. Also, convergence can be monitored with command-line keys -svd_monitor, -svd_monitor_all, -svd_monitor_conv, or graphically with -svd_monitor draw::draw_lg, or alternatively with -svd_monitor_all draw::draw_lg. See section Controlling and Monitoring Convergence for additional details.

Available possibilities for the convergence criterion, with \(Z=A\) in the standard SVD or as defined in equation (11) for the GSVD.

Convergence criterion	SVDConv	Command line key	Error bound
Absolute	SVD_CONV_ABS	-svd_conv_abs	\(|r|\)
Relative to eigenvalue	SVD_CONV_REL	-svd_conv_rel	\(|r|/|\lambda|\)
Relative to matrix norms	SVD_CONV_NORM	-svd_conv_norm	\(|r|/|Z|\)
Fixed number of iterations	SVD_CONV_MAXIT	-svd_conv_maxit	\(\infty\)
User-defined	SVD_CONV_USER	-svd_conv_user	user function

As in the case of eigensolvers, the user can choose different convergence tests, based on an error bound obtained from the computed residual norm, \(\|r\|\). Table Available possibilities for the convergence criterion, with Z=A in the standard SVD or as defined in equation eq:qr for the GSVD. lists the available options. It is worth mentioning the SVD_CONV_MAXIT convergence criterion, which is a bit special. With this criterion, the solver will not compute any residual norms and will stop with a successful status when the maximum number of iterations is reached. This is intended for the SVDRANDOMIZED solver in cases when a low-rank approximation of a matrix needs to be computed instead of accurate singular vectors.

Viewing the Solution

There is support for different kinds of viewers for the solution, as in the case of eigensolvers. One can for instance use -svd_view_values, -svd_view_vectors, -svd_error_relative, or -svd_converged_reason. See the description in section Viewing the Solution.

Footnotes

[Alvarruiz:2024:TLB]

F. Alvarruiz, C. Campos, and J. E. Roman. Thick-restarted joint Lanczos bidiagonalization for the GSVD. Journal of Computational and Applied Mathematics, 440:115506, 2024. doi:10.1016/j.cam.2023.115506.

[Anderson:1999:LUG]

E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users' Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999.

[Bai:2000:TSA]

Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, editors. Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2000.

[Bjorck:1996:NML]

Å. Björck. Numerical Methods for Least Squares Problems. Society for Industrial and Applied Mathematics, Philadelphia, 1996.

[Golub:1996:MC] (1,2)

G. H. Golub and C. F. van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, MD, third edition, 1996.

[Hansen:1998:RDI]

Per Christian Hansen. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1998.

[str-8] (1,2)

V. Hernandez, J. E. Roman, and A. Tomas. Restarted Lanczos Bidiagonalization for the SVD in SLEPc. Technical Report STR-8, Universitat Politècnica de València, 2007. Available at https://slepc.upv.es.

[Onn:1991:HSV]

R. Onn, A. O. Steinhardt, and A. Bojanczyk. The hyperbolic singular value decomposition and applications. IEEE Transactions on Signal Processing, 39(7):1575–1588, 1991. doi:10.1109/78.134396.

---

[1]

The exception is in the HSVD, where \(u_j\) is normalized so that \(U^*\Omega U=\tilde\Omega\).