SVD_LAPACK
SVD_LAPACK computes the singular value decomposition of a matrix by calling LAPACK subroutines: ZGESVD or ZGESDD.
Contents
Syntax
[U, S, V] = svd_lapack (A)
[U, S, V] = svd_lapack (A, ECON)
[U, S, V] = svd_lapack (A, ECON, OPT)
S = svd_lapack (A, ...)
[U, S, V, err] = svd_lapack (A, ...)
Description
Compute the singular value decomposition of a matrix A = U*S*V' where A is complex or real (not necessarily square), U and V are unitary or orthogonal, and S is always non-negative real diagonal.
The singular values reveal the rank of the matrix. Moreover, the singular value decomposition may be used to compute a pseudo inverse.
ECON argument makes S square and unitary matrices compact, where the unnecessary rows or columns of S, U and V are eliminated.
OPT argument selects the lapack routine ZGESVD or ZGESDD by assigning 'gesvd' or 'gesdd', where ZGESVD is based on QR iteration, and ZGESDD is a divide-and-conquer approach.
When called with one return value, it returns the sinular values S as a column vector.
When called with 4 return values, it returns the residual error: err = norm(A-U*S*V')/norm(A).
The lapack interface routines come from the work of Tim Toolan. See http://www.mathworks.co.kr/matlabcentral/fileexchange/16777-lapack
Examples
Generate a matrix with sigular values of [1 3 5 7]:
A = randn(7,4)+i*randn(7,4); [U,S,V] = svd(A,0); A = U*diag([1 3 5 7])*V';
Change format:
format short g
Compute singular values by selecting 'gesvd' implicitly:
svd_lapack(A)
ans =
7
5
3
1
Compute the economical form of SVD:
[U,S,V] = svd_lapack(A,0)
U =
Columns 1 through 2
0.019284 - 0.0031143i -0.25328 - 0.47749i
-0.036434 + 0.30821i -0.36886 + 0.43067i
-0.14784 - 0.38791i -0.2006 - 0.023319i
0.13317 + 0.30514i -0.071609 - 0.017322i
-0.12755 - 0.13695i 0.27878 + 0.29131i
-0.27704 + 0.59255i -0.096337 + 0.037905i
0.15081 - 0.36672i -0.32304 + 0.24989i
Columns 3 through 4
0.18937 - 0.068451i -0.11218 + 0.27533i
0.48378 - 0.14251i -0.48427 + 0.10378i
0.095815 + 0.13215i 0.23688 + 0.30733i
0.072506 - 0.22784i 0.54001 - 0.35323i
0.35702 + 0.46046i 0.18032 - 0.017114i
-0.46968 + 0.23118i -0.022319 + 0.11216i
-0.034712 - 0.080997i 0.10152 - 0.2078i
S =
7 0 0 0
0 5 0 0
0 0 3 0
0 0 0 1
V =
Columns 1 through 2
0.71198 0.3682
0.089184 - 0.33512i 0.56755 + 0.56103i
0.24561 - 0.50293i -0.33691 - 0.0021006i
-0.23786 + 0.054688i 0.23647 + 0.24109i
Columns 3 through 4
0.53771 0.26149
-0.46569 - 0.0071743i -0.084365 + 0.13722i
0.11736 + 0.58223i -0.43568 + 0.17504i
0.35993 + 0.1078i -0.42549 - 0.71005i
Select the lapack subroutine ZGESVD:
[U,S,V,err] = svd_lapack(A,0,'gesvd')
U =
Columns 1 through 2
0.019284 - 0.0031143i -0.25328 - 0.47749i
-0.036434 + 0.30821i -0.36886 + 0.43067i
-0.14784 - 0.38791i -0.2006 - 0.023319i
0.13317 + 0.30514i -0.071609 - 0.017322i
-0.12755 - 0.13695i 0.27878 + 0.29131i
-0.27704 + 0.59255i -0.096337 + 0.037905i
0.15081 - 0.36672i -0.32304 + 0.24989i
Columns 3 through 4
0.18937 - 0.068451i -0.11218 + 0.27533i
0.48378 - 0.14251i -0.48427 + 0.10378i
0.095815 + 0.13215i 0.23688 + 0.30733i
0.072506 - 0.22784i 0.54001 - 0.35323i
0.35702 + 0.46046i 0.18032 - 0.017114i
-0.46968 + 0.23118i -0.022319 + 0.11216i
-0.034712 - 0.080997i 0.10152 - 0.2078i
S =
7 0 0 0
0 5 0 0
0 0 3 0
0 0 0 1
V =
Columns 1 through 2
0.71198 0.3682
0.089184 - 0.33512i 0.56755 + 0.56103i
0.24561 - 0.50293i -0.33691 - 0.0021006i
-0.23786 + 0.054688i 0.23647 + 0.24109i
Columns 3 through 4
0.53771 0.26149
-0.46569 - 0.0071743i -0.084365 + 0.13722i
0.11736 + 0.58223i -0.43568 + 0.17504i
0.35993 + 0.1078i -0.42549 - 0.71005i
err =
8.5392e-16
Select the lapack subroutine ZGESDD:
[U,S,V,err] = svd_lapack(A,0,'gesdd')
U =
Columns 1 through 2
0.019284 - 0.0031143i -0.25328 - 0.47749i
-0.036434 + 0.30821i -0.36886 + 0.43067i
-0.14784 - 0.38791i -0.2006 - 0.023319i
0.13317 + 0.30514i -0.071609 - 0.017322i
-0.12755 - 0.13695i 0.27878 + 0.29131i
-0.27704 + 0.59255i -0.096337 + 0.037905i
0.15081 - 0.36672i -0.32304 + 0.24989i
Columns 3 through 4
0.18937 - 0.068451i -0.11218 + 0.27533i
0.48378 - 0.14251i -0.48427 + 0.10378i
0.095815 + 0.13215i 0.23688 + 0.30733i
0.072506 - 0.22784i 0.54001 - 0.35323i
0.35702 + 0.46046i 0.18032 - 0.017114i
-0.46968 + 0.23118i -0.022319 + 0.11216i
-0.034712 - 0.080997i 0.10152 - 0.2078i
S =
7 0 0 0
0 5 0 0
0 0 3 0
0 0 0 1
V =
Columns 1 through 2
0.71198 0.3682
0.089184 - 0.33512i 0.56755 + 0.56103i
0.24561 - 0.50293i -0.33691 - 0.0021006i
-0.23786 + 0.054688i 0.23647 + 0.24109i
Columns 3 through 4
0.53771 0.26149
-0.46569 - 0.0071743i -0.084365 + 0.13722i
0.11736 + 0.58223i -0.43568 + 0.17504i
0.35993 + 0.1078i -0.42549 - 0.71005i
err =
1.0157e-15
Diagnostics
Lapack, a fortran computational library, has two different subroutines for the Singular Value Decompostion (SVD): xGESVD and xGESDD. xGESVD is based on an implicit QR iteration and xGESDD uses a divide-and-conquer approach. See http://www.netlib.org/lapack/lug/node32.html and http://www.netlib.org/lapack/lug/node53.html for Lapack subroutines.
The source of ZGESVD can be found at http://www.netlib.org/lapack/explore-html/d6/d42/zgesvd_8f.html. ZGESVD calls other lapack subroutines:
- ZBDSQR computes the singular values and, optionally, the right and/or left singular vectors from the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B using the implicit zero-shift QR algorithm. See http://www.netlib.org/lapack/lawnspdf/lawn03.pdf.
- DLARTG generates a plane rotation with real cosine and real sine.
- ZDROT applies a plane rotation, where the cos and sin (c and s) are real and the vectors cx and cy are complex.
- ZLASR applies a sequence of plane rotations to a general rectangular matrix. simplified version needed to compare with.
The source of ZGESDD can be found at http://www.netlib.org/lapack/explore-html/d8/d54/zgesdd_8f_source.html. ZGESDD calls a subroutine:
- DBDSDC computes the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B = U * S * VT, using a divide and conquer method, where S is a diagonal matrix with non-negative diagonal elements (the singular values of B), and U and VT are orthogonal matrices of left and right singular vectors, respectively.
Matlab's built-in function svd seems to use the lapack subroutine xGESVD. Meanwhile, Octave's built-in function svd support both algorithms by using svd_driver().
If people want to use xGESDD routines in Matlab, it is required to write a MEX-file to call the routine. See http://www.mathworks.co.kr/matlabcentral/newsreader/view_thread/250196 But, it is not easy to call the lapack properly, because it needs a good knowledge of fortran work space.
SVD_LAPACK function provides both routines and the lapack interfaces. It could be an example for calling SVD rountine from Lapack interface. Although it supports double precision routines, they can be replaced by single precision routines without any problem.
xGESVD vs. xGESDD
Speaking of performance, xGESDD is faster than xGESVD. If only the singular values are needed to calculate, both show similar speed. See the following examples:
Generate a big matrix:
A=randn(1000,1000)+i*randn(1000,1000);
Measure times when only sinular values are calculated:
tic, S=svd(A,0); toc tic, S=svd_lapack(A,0,'gesvd'); toc tic, S=svd_lapack(A,0,'gesdd'); toc
Elapsed time is 4.230231 seconds. Elapsed time is 4.416975 seconds. Elapsed time is 4.300102 seconds.
Measure times when the complete SVD is calculated:
tic, [U,S,V]=svd(A,0); toc tic, [U,S,V]=svd_lapack(A,0,'gesvd'); toc tic, [U,S,V]=svd_lapack(A,0,'gesdd'); toc
Elapsed time is 7.816645 seconds. Elapsed time is 14.802859 seconds. Elapsed time is 7.770088 seconds.
In Octave, we can use svd_driver() to select subroutines:
svd_driver('gesvd'), tic, [U,S,V]=svd(A,0); toc svd_driver('gesdd'), tic, [U,S,V]=svd(A,0); toc
Convergence test
SVD has implemented in many numerical libraries such as LINPACK, EISPACK, and Lapack. See http://www.alglib.net/matrixops/general/svd.php. Among them, Lapack's SVD focuses the smallest singular value. The following example is to see the small singular values as in http://www.netlib.org/lapack/lawnspdf/lawn03.pdf:
Generate a matrix A:
eta = eps(1)/2; x = eta; % the smallest sigma ~ eta^3 %x = 0; % the smallest sigma ~ eta^2/sqrt(2) % good test for deflation A = [eta^2 1 0 0; 0 1 x 0; 0 0 1 1; 0 0 0 eta^2];
Change format:
format long g
Calculate SVD using ZGESVD:
[U,S,V,err] = svd_lapack(A,[],'gesvd')
U =
Columns 1 through 3
0.707106781186547 -2.77555756156289e-17 0.707106781186547
0.707106781186547 2.77555756156289e-17 -0.707106781186547
0 1 3.92523114670944e-17
0 6.16297582203915e-33 1.57009245868378e-16
Column 4
-1.11022302462516e-16
1.11022302462516e-16
-1.23259516440783e-32
1
S =
Columns 1 through 3
1.4142135623731 0 0
0 1.41421356237309 0
0 0 5.55111512312578e-17
0 0 0
Column 4
0
0
0
1.3684555315672e-48
V =
Columns 1 through 3
1.23259516440783e-32 0 1.57009245868377e-16
1 0 -1.93528837224682e-48
0 0.707106781186547 -0.707106781186547
0 0.707106781186547 0.707106781186547
Column 4
-1
1.23259516440783e-32
-1.11022302462516e-16
1.11022302462516e-16
err =
2.25593997819691e-16
Calculate SVD using ZGESDD:
[U,S,V,err] = svd_lapack(A,[],'gesdd')
U =
Columns 1 through 3
0.707106781186547 -2.77555756156289e-17 0.707106781186547
0.707106781186547 2.77555756156289e-17 -0.707106781186547
0 1 3.92523114670944e-17
0 6.16297582203915e-33 1.57009245868378e-16
Column 4
-1.11022302462516e-16
1.11022302462516e-16
-1.23259516440783e-32
1
S =
Columns 1 through 3
1.4142135623731 0 0
0 1.41421356237309 0
0 0 5.55111512312578e-17
0 0 0
Column 4
0
0
0
1.3684555315672e-48
V =
Columns 1 through 3
1.23259516440783e-32 0 1.57009245868377e-16
1 0 -1.93528837224682e-48
0 0.707106781186547 -0.707106781186547
0 0.707106781186547 0.707106781186547
Column 4
-1
1.23259516440783e-32
-1.11022302462516e-16
1.11022302462516e-16
err =
2.25593997819691e-16
Pseudo inverse
SVD can be used to find a Pseudo inverse. For example,
Change format:
format short g
Generate A with sigular values of [1 3 5 7]:
A = randn(7,4)+i*randn(7,4); [U,S,V] = svd(A,0); A = U*diag([1 3 5 7])*V'
A =
Columns 1 through 2
-0.22876 - 0.19423i -0.011353 + 2.4587i
0.072388 - 2.3155i -1.1755 - 1.088i
-0.24808 + 0.25466i 0.35265 + 1.161i
1.1915 - 0.99036i 2.1921 + 0.51522i
-0.27298 - 2.8106i 0.2178 - 2.8519i
0.13228 + 1.2026i -0.22827 + 1.1462i
0.71953 + 0.62923i -0.057724 + 0.58474i
Columns 3 through 4
-0.87186 - 1.3685i 0.64031 - 0.61814i
-1.0071 - 0.46738i 0.60154 - 1.6568i
-0.63273 - 1.6053i 0.8606 + 0.17557i
-0.80766 + 0.21813i -1.2019 - 2.3228i
-0.072746 - 0.3348i -0.42885 + 2.9651i
-0.56978 + 1.206i -0.26363 - 1.3408i
0.34944 - 1.8386i 1.7471 + 1.5072i
Calculate SVD:
[U,S,V,err] = svd_lapack(A,0,'gesdd')
U =
Columns 1 through 2
0.0037729 + 0.33929i 0.10335 - 0.070657i
-0.051937 - 0.0094506i -0.12665 - 0.61861i
-0.063863 + 0.19991i 0.2212 - 0.049661i
0.51099 + 0.055765i 0.0084541 - 0.24813i
-0.20791 - 0.62146i 0.17541 - 0.31523i
0.15372 + 0.23565i -0.24697 + 0.12253i
-0.227 + 0.117i 0.50911 + 0.085028i
Columns 3 through 4
-0.39517 - 0.47651i -0.18225 + 0.21106i
0.36011 - 0.32907i -0.01546 - 0.15897i
-0.41357 - 0.16381i 0.16629 - 0.30977i
0.015274 - 0.19551i 0.382 + 0.17691i
-0.27042 - 0.06192i -0.37182 + 0.34768i
0.12709 + 0.16134i -0.55981 + 0.063587i
-0.1523 + 0.011652i -0.077497 - 0.11546i
S =
7 0 0 0
0 5 0 0
0 0 3 0
0 0 0 1
V =
Columns 1 through 2
0.36988 0.59485
0.61537 - 0.12674i 0.34248 - 0.020485i
-0.13808 - 0.20164i 0.14098 + 0.52547i
-0.47381 + 0.42922i 0.34299 - 0.33884i
Columns 3 through 4
0.53337 -0.4742
-0.44883 + 0.28757i 0.40477 + 0.1989i
0.44845 - 0.21283i 0.57356 + 0.26249i
0.096152 + 0.41919i 0.16883 + 0.38125i
err =
7.4966e-16
Find the inverse of singular values:
S_diag = diag(S); idx_non_zero_S_diag = (S_diag ~= 0); S_diag_inv = S_diag; S_diag_inv(idx_non_zero_S_diag) = S_diag_inv(idx_non_zero_S_diag).^(-1); S_inv = [diag(S_diag_inv); zeros(size(S,2)-size(S,1),size(S,1))]
S_inv =
0.14286 0 0 0
0 0.2 0 0
0 0 0.33333 0
0 0 0 1
Find the Pseudo inverse:
A_inv = V*S_inv*U'
A_inv =
Columns 1 through 2
0.028657 + 0.17528i 0.053544 + 0.057216i
-0.016788 - 0.25633i -0.13383 + 0.091224i
-0.088755 - 0.05019i -0.040707 + 0.11621i
0.0028751 - 0.12402i -0.061487 + 0.129i
Columns 3 through 4
-0.12944 - 0.12242i -0.15042 + 0.14523i
0.057987 + 0.080387i 0.21429 - 0.020607i
-0.039617 + 0.30558i 0.24417 + 0.02121i
-0.091049 + 0.061188i 0.091339 + 0.17573i
Columns 5 through 6
0.13812 + 0.24622i 0.2668 - 0.025561i
-0.040543 - 0.1706i -0.22567 - 0.13165i
-0.16421 - 0.24752i -0.30075 - 0.24574i
0.061822 - 0.28133i -0.064852 - 0.17787i
Column 7
0.058245 - 0.073122i
-0.017982 + 0.0043803i
-0.073953 + 0.1149i
-0.0086562 - 0.07804i
Measure inversion error:
err_inv = norm(A_inv * A - eye(size(A,2)))
err_inv = 1.7661e-15
