svds
Subset of singular values and vectors
Syntax
Description
s = svds(A,k,sigma,Name,Value)svds(A,k,sigma,'Tolerance',1e-3) adjusts the
convergence tolerance for the algorithm.
Examples
Largest Singular Values
The matrix A = delsq(numgrid('C',15)) is a symmetric positive definite matrix with singular values reasonably well-distributed in the interval (0 8). Compute the six largest singular values.
A = delsq(numgrid('C',15));
s = svds(A)s = 6×1
7.8666
7.7324
7.6531
7.5213
7.4480
7.3517
Specify a second input to compute a specific number of the largest singular values.
s = svds(A,3)
s = 3×1
7.8666
7.7324
7.6531
Smallest Singular Values
The matrix A = delsq(numgrid('C',15)) is a symmetric positive definite matrix with singular values reasonably well-distributed in the interval (0 8). Compute the five smallest singular values.
A = delsq(numgrid('C',15)); s = svds(A,5,'smallest')
s = 5×1
0.5520
0.4787
0.3469
0.2676
0.1334
Smallest Nonzero Singular Values
Create a sparse 100-by-100 Neumann matrix.
C = gallery('neumann',100);Compute the ten smallest singular values.
ss = svds(C,10,'smallest')ss = 10×1
0.9828
0.9049
0.5625
0.5625
0.4541
0.4506
0.2256
0.1139
0.1139
0
Compute the 10 smallest nonzero singular values. Since the matrix has a singular value that is equal to zero, the 'smallestnz' option omits it.
snz = svds(C,10,'smallestnz')snz = 10×1
0.9828
0.9828
0.9049
0.5625
0.5625
0.4541
0.4506
0.2256
0.1139
0.1139
Largest Singular Values Using Function Handle
Create two matrices representing the upper-right and lower-left nonzero blocks in a sparse matrix.
n = 500; B = rand(500); C = rand(500);
Save Afun in your current directory so that it is available for use with svds.
function y = Afun(x,tflag,B,C,n) if strcmp(tflag,'notransp') y = [B*x(n+1:end); C*x(1:n)]; else y = [C'*x(n+1:end); B'*x(1:n)]; end
The function Afun uses B and C to compute either A*x or A'*x (depending on the specified flag) without actually forming the entire sparse matrix A = [zeros(n) B; C zeros(n)]. This exploits the sparsity pattern of the matrix to save memory in the computation of A*x and A'*x.
Use Afun to calculate the 10 largest singular values of A. Pass B, C, and n as additional inputs to Afun.
s = svds(@(x,tflag) Afun(x,tflag,B,C,n),[1000 1000],10)
s = 250.3248 249.9914 12.7627 12.7232 12.6988 12.6608 12.6166 12.5643 12.5419 12.4512
Directly compute the 10 largest singular values of A to compare the results.
A = [zeros(n) B; C zeros(n)]; s = svds(A,10)
s = 250.3248 249.9914 12.7627 12.7232 12.6988 12.6608 12.6166 12.5643 12.5419 12.4512
Singular Value Decomposition of Sparse Matrix
west0479 is a real-valued 479-by-479 sparse matrix. The matrix has a few large singular values, and many small singular values.
Load west0479 and store it as A.
load west0479
A = west0479;Compute the singular value decomposition of A, returning the six largest singular values and the corresponding singular vectors. Specify a fourth output argument to check convergence of the singular values.
[U,S,V,cflag] = svds(A); cflag
cflag = 0
cflag indicates that all of the singular values converged. The singular values are on the diagonal of the output matrix S.
s = diag(S)
s = 6×1
105 ×
3.1895
3.1725
3.1695
3.1685
3.1669
0.3038
Check the results by computing the full singular value decomposition of A. Convert A to a full matrix and use svd.
[U1,S1,V1] = svd(full(A));
Plot the six largest singular values of A computed by svd and svds using a logarithmic scale.
s2 = diag(S1); semilogy(s2(1:6),'r.') hold on semilogy(s,'ro','MarkerSize',10) title('Singular Values of west0479') legend('svd','svds')
Fix Convergence Problem
Create a sparse diagonal matrix and calculate the six largest singular values.
A = diag(sparse([1e4*ones(1, 8) 1e4:-1:1])); s = svds(A)
Warning: Only 2 of the 6 requested singular values converged. Singular values that did not converge are NaN.
s = 6×1
104 ×
1.0000
0.9999
NaN
NaN
NaN
NaN
The svds algorithm produces a warning since the maximum number of iterations were performed but the tolerance could not be met.
The most effective way to address convergence problems is to increase the maximum size of the Krylov subspace used in the calculation by using a larger value for 'SubspaceDimension'. Do this by passing in the name-value pair 'SubspaceDimension' with a value of 60.
s = svds(A,6,'largest','SubspaceDimension',60)
s = 6×1
104 ×
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
Use QR Decomposition to Compute SVD of Nearly Singular Matrix
Compute the 10 smallest singular values of a nearly singular matrix.
rng default format shortg B = spdiags([repelem([1; 1e-7], [198, 2]) ones(200, 1)], [0 1], 200, 200); s1 = svds(B,10,'smallest')
Warning: Large residual norm detected. This is likely due to bad condition of the input matrix (condition number 1.0008e+16).
s1 = 10×1
7.0945
7.0945
7.0945
7.0945
7.0945
7.0945
7.0945
7.0945
0.25927
7.0888e-16
The warning indicates that svds fails to calculate the proper singular values. The failure with svds is because of the gap between the smallest and second smallest singular values. svds(...,'smallest') needs to invert B, which leads to large numerical error.
For comparison, compute the exact singular values using svd.
s = svd(full(B)); s = s(end-9:end)
s = 10×1
0.14196
0.12621
0.11045
0.094686
0.078914
0.063137
0.047356
0.031572
0.015787
7.0888e-16
In order to reproduce this calculation with svds, do a QR decomposition of B. The singular values of the triangular matrix R are the same as for B.
[Q,R,p] = qr(B,0);
Plot the norm of each row of R.
rownormR = sqrt(diag(R*R')); semilogy(rownormR) hold on; semilogy(size(R, 1), rownormR(end), 'ro')
The last entry in R is nearly zero, which causes instability in the solution.
Prevent this entry from corrupting the good parts of the solution by setting the last row of R to be exactly zero.
R(end,:) = 0;
Use svds to find the 10 smallest singular values of R. The results are comparable to those obtained by svd.
sr = svds(R,10,'smallest')sr = 10×1
0.14196
0.12621
0.11045
0.094686
0.078914
0.063137
0.047356
0.031572
0.015787
0
To compute the singular vectors of B using this method, transform the left and right singular vectors using Q and the permutation vector p.
[U,S,V] = svds(R,20,'s');
U = Q*U;
V(p,:) = V;Input Arguments
Output Arguments
Tips
svdsketchis useful when you do not know ahead of time what rank to specify withsvds, but you know what tolerance the approximation of the SVD should satisfy.svdsgenerates the default starting vectors using a private random number stream to ensure reproducibility across runs. Setting the random number generator state usingrngbefore callingsvdsdoes not affect the output.Using
svdsis not the most efficient way to find a few singular values of small, dense matrices. For such problems, usingsvd(full(A))might be quicker. For example, finding three singular values in a 500-by-500 matrix is a relatively small problem thatsvdcan handle easily.If
svdsfails to converge for a given matrix, increase the size of the Krylov subspace by increasing the value of'SubspaceDimension'. As secondary options, adjusting the maximum number of iterations ('MaxIterations') and the convergence tolerance ('Tolerance') also can help with convergence behavior.Increasing
kcan sometimes improve performance, especially when the matrix has repeated singular values.
Compatibility Considerations
References
[1] Baglama, J. and L. Reichel, “Augmented Implicitly Restarted Lanczos Bidiagonalization Methods.” SIAM Journal on Scientific Computing. Vol. 27, 2005, pp. 19–42.
[2] Larsen, R. M. “Lanczos Bidiagonalization with partial reorthogonalization.” Dept. of Computer Science, Aarhus University. DAIMI PB-357, 1998.