Singular Value Decomposition - Algorithms and Error Analysis

We study only algorithms for real matrices, which are most commonly used in the applications described in this course.

For more details, see A. Kaylor Cline and I. Dhillon, Computation of the Singular Value Decomposition and the references therein.

Prerequisites

The reader should be familiar with facts about the singular value decomposition and perturbation theory and algorithms for the symmetric eigenvalue decomposition.

Competences

The reader should be able to apply an adequate algorithm to a given problem, and to assess the accuracy of the solution.

Basics

Definitions

The singular value decomposition (SVD) of $A\in\mathbb{R}^{m\times n}$ is $A=U\Sigma V^T$, where $U\in\mathbb{R}^{m\times m}$ is orthogonal, $U^TU=UU^T=I_m$, $V\in\mathbb{R}^{n\times n}$ is orthogonal, $V^TV=VV^T=I_n$, and $\Sigma \in\mathbb{R}^{m\times n}$ is diagonal with singular values $\sigma_1,\ldots,\sigma_{\min\{m,n\}}$ on the diagonal.

If $m>n$, the thin SVD of $A$ is $A=U_{1:m,1:n} \Sigma_{1:n,1:n} V^T$.

Facts

  1. Algorithms for computing SVD of $A$ are modifications of algorithms for the symmetric eigenvalue decomposition of the matrices $AA^T$, $A^TA$ and $\begin{bmatrix} 0 & A\\ A^T & 0 \end{bmatrix}$.

  2. Most commonly used approach is the three-step algorithm:

    1. Reduce $A$ to bidiagonal matrix $B$ by orthogonal transformations, $X^TAY=B$.
    2. Compute the SVD of $B$ with QR iterations, $B=W\Sigma Z^T$.
    3. Multiply $U=XW$ and $V=YZ$.
  3. If $m\geq n$, the overall operation count for this algorithm is $O(mn^2)$ operations.

  4. Error bounds. Let $U\Sigma V^T$ and $\tilde U \tilde \Sigma \tilde V^T$ be the exact and the computed SVDs of $A$, respectively. The algorithms generally compute the SVD with errors bounded by $$ |\sigma_i-\tilde \sigma_i|\leq \phi \epsilon\|A\|_2, \qquad \|u_i-\tilde u_i\|_2, \| v_i-\tilde v_i\|_2 \leq \psi\epsilon \frac{\|A\|_2} {\min_{j\neq i} |\sigma_i-\tilde \sigma_j|}, $$ where $\epsilon$ is machine precision and $\phi$ and $\psi$ are slowly growing polynomial functions of $n$ which depend upon the algorithm used (typically $O(n)$ or $O(n^2)$). These bounds are obtained by combining perturbation bounds with the floating-point error analysis of the algorithms.

Bidiagonalization

Facts

  1. The reduction of $A$ to bidiagonal matrix can be performed by applying $\min\{m-1,n\}$ Householder reflections $H_L$ from the left and $n-2$ Householder reflections $H_R$ from the right. In the first step, $H_L$ is chosen to annihilate all elements of the first column below the diagonal, and $H_R$ is chosen to annihilate all elements of the first row right of the first super-diagonal. Applying this procedure recursively yields the bidiagonal matrix.

  2. $H_L$ and $H_R$ do not depend on the normalization of the respective Householder vectors $v_L$ and $v_R$. With the normalization $[v_L]_1=[V_R]_1=1$, the vectors $v_L$ are stored in the lower-triangular part of $A$, and the vectors $v_R$ are stored in the upper-triangular part of $A$ above the super-diagonal.

  3. The matrices $H_L$ and $H_R$ are not formed explicitly - given $v_L$ and $v_R$, $A$ is overwritten with $H_L A H_R$ in $O(mn)$ operations by using matrix-vector multiplication and rank-one updates.

  4. Instead of performing rank-one updates, $p$ transformations can be accumulated, and then applied. This block algorithm is rich in matrix-matrix multiplications (roughly one half of the operations is performed using BLAS 3 routines), but it requires extra workspace.

  5. If the matrices $X$ and $Y$ are needed explicitly, they can be computed from the stored Householder vectors. In order to minimize the operation count, the computation starts from the smallest matrix and the size is gradually increased.

  6. The backward error bounds for the bidiagonalization are as follows: The computed matrix $\tilde B$ is equal to the matrix which would be obtained by exact bidiagonalization of some perturbed matrix $A+\Delta A$, where $\|\Delta A\|_2 \leq \psi \varepsilon \|A\|_2$ and $\psi$ is a slowly increasing function of $n$. The computed matrices $\tilde X$ and $\tilde Y$ satisfy $\tilde X=X+\Delta X$ and $\tilde Y=Y+\Delta Y$, where $\|\Delta X \|_2,\|\Delta Y\|_2\leq \phi \varepsilon$ and $\phi$ is a slowly increasing function of $n$.

  7. The bidiagonal reduction is implemented in the LAPACK subroutine DGEBRD. The computation of $X$ and $Y$ is implemented in the subroutine DORGBR, which is not yet wrapped in Julia.

  8. Bidiagonalization can also be performed using Givens rotations. Givens rotations act more selectively than Householder reflectors, and are useful if $A$ has some special structure, for example, if $A$ is a banded matrix. Error bounds for function myBidiagG() are the same as above, but with slightly different functions $\psi$ and $\phi$.


In [36]:
m=8
n=5
import Random
Random.seed!(421)
A=rand(-9.0:9,m,n)


Out[36]:
8×5 Array{Float64,2}:
  9.0  -8.0  -6.0  -8.0   6.0
 -7.0  -3.0  -3.0  -5.0   1.0
 -6.0   5.0   4.0  -4.0   9.0
 -6.0  -9.0   9.0   4.0   2.0
  6.0  -3.0   8.0   2.0  -9.0
 -5.0   3.0   4.0  -9.0   9.0
  5.0   4.0   6.0  -6.0   4.0
 -3.0  -1.0   5.0   4.0   3.0

In [39]:
function H(x)
    v=copy(x)
    v[1]+=sign(x[1])*norm(x)
    display(v/v[1])
    I-(2/(v⋅v))*v*v'
end


Out[39]:
H (generic function with 1 method)

In [40]:
H1=H(A[:,1])


8-element Array{Float64,1}:
  1.0
 -0.26683247952453054
 -0.22871355387816905
 -0.22871355387816905
  0.22871355387816905
 -0.19059462823180753
  0.19059462823180753
 -0.11435677693908453
Out[40]:
8×8 Array{Float64,2}:
 -0.522233   0.406181    0.348155   …   0.290129   -0.290129    0.174078
  0.406181   0.891618   -0.0928991     -0.077416    0.077416   -0.0464496
  0.348155  -0.0928991   0.920372      -0.0663565   0.0663565  -0.0398139
  0.348155  -0.0928991  -0.0796278     -0.0663565   0.0663565  -0.0398139
 -0.348155   0.0928991   0.0796278      0.0663565  -0.0663565   0.0398139
  0.290129  -0.077416   -0.0663565  …   0.944703    0.0552971  -0.0331783
 -0.290129   0.077416    0.0663565      0.0552971   0.944703    0.0331783
  0.174078  -0.0464496  -0.0398139     -0.0331783   0.0331783   0.980093

In [16]:
H1'*H1


Out[16]:
8×8 Array{Float64,2}:
  1.0          -5.92989e-17  -4.22513e-17  …  -1.41681e-17  -4.90355e-17
 -5.92989e-17   1.0           8.95151e-18     -2.85312e-17   1.40079e-17
 -4.22513e-17   8.95151e-18   1.0             -1.56306e-17  -2.82284e-18
 -4.22513e-17   1.58904e-17  -1.12064e-17     -1.56306e-17   4.11605e-18
  4.91902e-17  -2.4564e-17    1.64106e-17      1.56306e-17  -4.11605e-18
 -3.17913e-18   2.7447e-17    1.32454e-17  …  -3.04723e-17   7.96389e-18
 -1.41681e-17  -2.85312e-17  -1.56306e-17      1.0          -7.96389e-18
 -4.90355e-17   1.40079e-17  -2.82284e-18     -7.96389e-18   1.0

In [19]:
B=H1*A


Out[19]:
8×5 Array{Float64,2}:
 -17.2337         2.14696    3.94576    1.27657   6.20877
   1.04083e-15   -5.70754   -5.65385   -7.47529   0.944293
   5.55112e-16    2.67925    1.72527   -6.12168   8.95225
   1.11022e-16  -11.3207     6.72527    1.87832   1.95225
   3.33067e-16   -0.679253  10.2747     4.12168  -8.95225
   4.09395e-16    1.06604    2.10439  -10.7681    8.96021
   6.93889e-18    5.93396    7.89561   -4.23194   4.03979
  -3.33067e-16   -2.16037    3.86263    2.93916   2.97613

In [35]:
H2=H(B[1,2:5])


Out[35]:
4×4 SparseArrays.SparseMatrixCSC{Float64,Int64} with 4 stored entries:
  [1, 1]  =  -1.0
  [2, 2]  =  1.0
  [3, 3]  =  1.0
  [4, 4]  =  1.0

In [23]:
H3=[1 zeros(1,4);zeros(4) H2]


Out[23]:
5×5 Array{Float64,2}:
 1.0   0.0        0.0        0.0        0.0
 0.0  -0.27635   -0.507887  -0.164316  -0.799175
 0.0  -0.507887   0.797901  -0.065385  -0.318009
 0.0  -0.164316  -0.065385   0.978846  -0.102885
 0.0  -0.799175  -0.318009  -0.102885   0.499604

In [24]:
B*H3


Out[24]:
8×5 Array{Float64,2}:
 -17.2337       -7.76897  -2.45528e-17    3.56275e-17  -4.02705e-16
   1.04083e-15   4.92246  -1.42395       -6.10679       7.60017
   5.55112e-16  -7.76518  -2.43079       -7.46629       2.41256
   1.11022e-16  -2.15602  10.3721         3.05818       7.69067
   3.33067e-16   1.44647  11.1206         4.39534      -7.62125
   4.09395e-16  -6.7548   -1.00769      -11.7749        4.06326
   6.93889e-18  -8.18305   2.27815       -6.04935      -4.79945
  -3.33067e-16  -4.22616   3.04061        2.67321       1.68265

In [25]:
H1*A*H3


Out[25]:
8×5 Array{Float64,2}:
 -17.2337       -7.76897  -2.45528e-17    3.56275e-17  -4.02705e-16
   1.04083e-15   4.92246  -1.42395       -6.10679       7.60017
   5.55112e-16  -7.76518  -2.43079       -7.46629       2.41256
   1.11022e-16  -2.15602  10.3721         3.05818       7.69067
   3.33067e-16   1.44647  11.1206         4.39534      -7.62125
   4.09395e-16  -6.7548   -1.00769      -11.7749        4.06326
   6.93889e-18  -8.18305   2.27815       -6.04935      -4.79945
  -3.33067e-16  -4.22616   3.04061        2.67321       1.68265

In [2]:
using LinearAlgebra

In [41]:
?LAPACK.gebrd!


Out[41]:
gebrd!(A) -> (A, d, e, tauq, taup)

Reduce A in-place to bidiagonal form A = QBP'. Returns A, containing the bidiagonal matrix B; d, containing the diagonal elements of B; e, containing the off-diagonal elements of B; tauq, containing the elementary reflectors representing Q; and taup, containing the elementary reflectors representing P.


In [42]:
# We need copy()
Outg=LAPACK.gebrd!(copy(A))


Out[42]:
([-17.233687939614086 -7.76897048243217 … 0.128739281508655 0.6261410509739129; -0.2668324795245306 -14.890247384185125 … -0.8769502559968347 0.025394350692863336; … ; 0.19059462823180756 -0.41302052973946735 … 0.36078939474593574 0.27938709255633576; -0.11435677693908453 -0.21330565047588582 … -0.14200886340449123 -0.27274236959087667], [-17.233687939614086, -14.890247384185125, 13.475190619553185, 13.673356118350501, 14.193966697750746], [-7.76897048243217, -9.85421673756534, -10.636159259139372, -1.2961631050111961, 1.58966657e-315], [1.5222329678670936, 1.3305825317403805, 1.4211703659081718, 1.1449365281707868, 1.7349788922897473], [1.2763503557700235, 1.130143592806911, 1.8042031786817188, 0.0, 0.0])

In [45]:
Outg[4]


Out[45]:
5-element Array{Float64,1}:
 1.5222329678670936
 1.3305825317403805
 1.4211703659081718
 1.1449365281707868
 1.7349788922897473

In [28]:
B=Bidiagonal(Outg[2],Outg[3][1:end-1],'U')


Out[28]:
5×5 Bidiagonal{Float64,Array{Float64,1}}:
 -17.2337   -7.76897    ⋅          ⋅        ⋅ 
    ⋅      -14.8902   -9.85422     ⋅        ⋅ 
    ⋅         ⋅       13.4752   -10.6362    ⋅ 
    ⋅         ⋅         ⋅        13.6734  -1.29616
    ⋅         ⋅         ⋅          ⋅      14.194

In [32]:
[svdvals(A) svdvals(B)]


Out[32]:
5×2 Array{Float64,2}:
 22.7144   22.7144
 19.4392   19.4392
 14.7417   14.7417
 13.9812   13.9812
  7.37431   7.37431

In [46]:
# Extract X
function myBidiagX(H::Matrix)
    m,n=size(H)
    T=typeof(H[1,1])
    X = Matrix{T}(I,m,n)
    v = Array{T}(undef,m)
    for j = n : -1 : 1
        v[j] = one(T)
        v[j+1:m] = H[j+1:m, j]
        γ = -2 / (v[j:m]v[j:m])
        w = γ * X[j:m, j:n]'*v[j:m]
        X[j:m, j:n] = X[j:m, j:n] + v[j:m]*w'
    end
    X
end

# Extract Y
function myBidiagY(H::AbstractMatrix)
    n,m=size(H)
    T=typeof(H[1,1])
    Y = Matrix{T}(I,n,n)
    v = Array{T}(undef,n)
    for j = n-2 : -1 : 1
        v[j+1] = one(T)
        v[j+2:n] = H[j+2:n, j]
        γ = -2 / (v[j+1:n]v[j+1:n])
        w = γ * Y[j+1:n, j+1:n]'*v[j+1:n]
        Y[j+1:n, j+1:n] = Y[j+1:n, j+1:n] + v[j+1:n]*w'
    end
    Y
end


Out[46]:
myBidiagY (generic function with 1 method)

In [47]:
X=myBidiagX(Outg[1])


Out[47]:
8×5 Array{Float64,2}:
 -0.522233   0.153094   -0.389298    0.375871    0.602756
  0.406181  -0.371433   -0.589452   -0.19564     0.169813
  0.348155   0.48648    -0.0603962   0.0930414  -0.130866
  0.348155   0.109779    0.303295   -0.234362    0.726816
 -0.348155  -0.0621273   0.0720229  -0.799631    0.0912441
  0.290129   0.424461   -0.45901    -0.17818    -0.0117797
 -0.290129   0.578737   -0.129015   -0.268059   -0.143864
  0.174078   0.266313    0.415734    0.117688    0.182486

In [48]:
Y=myBidiagY(Outg[1]')


Out[48]:
5×5 Array{Float64,2}:
 1.0   0.0        0.0         0.0         0.0
 0.0  -0.27635   -0.0738167  -0.0301716  -0.957743
 0.0  -0.507887  -0.159517   -0.826093    0.184866
 0.0  -0.164316   0.981577   -0.0941737  -0.0252745
 0.0  -0.799175  -0.0749189   0.55479     0.218894

In [49]:
# Orthogonality
norm(X'*X-I), norm(Y'*Y-I)


Out[49]:
(1.1982936726328904e-15, 5.960014444977975e-16)

In [50]:
# Residual error
norm(A*Y-X*B)


Out[50]:
1.2134615767632303e-14

In [53]:
# Bidiagonalization using Givens rotations
function myBidiagG(A1::Matrix)
    A=deepcopy(A1)
    m,n=size(A)
    T=typeof(A[1,1])
    X=Matrix{T}(I,m,m)
    Y=Matrix{T}(I,n,n)
    for j = 1 : n        
        for i = j+1 : m
            G,r=givens(A,j,i,j)
            # Use the faster in-place variant
            # A=G*A
            lmul!(G,A)
            # X=G*X
            lmul!(G,X)
            # display(A)
        end
        for i=j+2:n
            G,r=givens(A',j+1,i,j)
            # A*=adjoint(G)
            rmul!(A,adjoint(G))
            # Y*=adjoint(G)
            rmul!(Y,adjoint(G))
            # display(A)
        end
    end
    X',Bidiagonal(diag(A),diag(A,1),'U'), Y
end


Out[53]:
myBidiagG (generic function with 1 method)

In [54]:
X₁, B₁, Y₁ = myBidiagG(A)


Out[54]:
([0.5222329678670936 -0.15309399920139824 … -0.15076651875883973 0.07009783931517549; -0.4061811972299618 0.37143298314761614 … 0.32864415280786885 0.32201295739000246; … ; 0.2901294265928297 -0.5787367959838046 … 0.6867321565310394 0.00948705221054524; -0.1740776559556978 -0.26631343612734604 … 0.0 0.823885833236553], 5×5 Bidiagonal{Float64,Array{Float64,1}}:
 diag: 17.233687939614086  14.890247384185123  …  -14.193966697750746
 super: 7.76897048243217  9.854216737565345  10.636159259139376  1.2961631050112015, [1.0 -0.0 … 0.0 -0.0; 0.0 -0.27635035577002326 … -0.030171588021228657 0.9577427912919487; … ; 0.0 -0.16431642775514896 … -0.09417365216726617 0.025274480905889372; 0.0 -0.7991753531727701 … 0.5547895992064668 -0.21889360819642828])

In [55]:
# Orthogonality
norm(X₁'*X₁-I), norm(Y₁'*Y₁-I)


Out[55]:
(1.3649803812779788e-15, 4.217501041969514e-16)

In [56]:
# Diagonalization
X₁'*A*Y₁


Out[56]:
8×5 Array{Float64,2}:
 17.2337        7.76897        8.9948e-16    1.19459e-15   -3.86745e-16
  1.11022e-15  14.8902         9.85422      -1.7583e-15     2.2406e-16
  1.11022e-15  -9.38165e-17  -13.4752       10.6362        -4.84943e-16
 -1.01308e-15  -2.70843e-15   -1.73276e-15  13.6734         1.29616
  5.55112e-17  -3.48249e-16    2.18673e-16  -1.91564e-15  -14.194
  7.77156e-16  -7.21419e-17    1.83272e-15   1.11989e-15   -1.58309e-15
  8.88178e-16  -1.05642e-15    1.69323e-15   4.62053e-16   -6.23447e-16
  1.44329e-15  -4.30364e-16   -5.11474e-16   1.86434e-16    1.90036e-16

In [57]:
# X, Y and B are not unique
B


Out[57]:
5×5 Bidiagonal{Float64,Array{Float64,1}}:
 -17.2337   -7.76897    ⋅          ⋅        ⋅ 
    ⋅      -14.8902   -9.85422     ⋅        ⋅ 
    ⋅         ⋅       13.4752   -10.6362    ⋅ 
    ⋅         ⋅         ⋅        13.6734  -1.29616
    ⋅         ⋅         ⋅          ⋅      14.194

In [58]:
B₁


Out[58]:
5×5 Bidiagonal{Float64,Array{Float64,1}}:
 17.2337   7.76897     ⋅         ⋅         ⋅ 
   ⋅      14.8902     9.85422    ⋅         ⋅ 
   ⋅        ⋅       -13.4752   10.6362     ⋅ 
   ⋅        ⋅          ⋅       13.6734    1.29616
   ⋅        ⋅          ⋅         ⋅      -14.194

In [61]:
Matrix(B)'*Matrix(B)


Out[61]:
5×5 Array{Float64,2}:
 297.0    133.888     0.0       0.0       0.0
 133.888  282.076   146.732     0.0       0.0
   0.0    146.732   278.686  -143.324     0.0
   0.0      0.0    -143.324   300.089   -17.7229
   0.0      0.0       0.0     -17.7229  203.149

Bidiagonal QR method

Let $B$ be a real upper-bidiagonal matrix of order $n$ and let $B=W\Sigma Z^T$ be its SVD.

All metods for computing the SVD of bidiagonal matrix are derived from the methods for computing the EVD of the tridiagonal matrix $T=B^T B$.

Facts

  1. The shift $\mu$ is the eigenvalue of the $2\times 2$ matrix $T_{n-1:n,n-1:n}$ which is closer to $T_{n,n}$. The first Givens rotation from the right is the one which annihilates the element $(1,2)$ of the shifted $2\times 2$ matrix $T_{1:2,1:2}-\mu I$. Applying this rotation to $B$ creates the bulge at the element $B_{2,1}$. This bulge is subsequently chased out by applying adequate Givens rotations alternating from the left and from the right. This is the Golub-Kahan algorithm.

  2. The computed SVD satisfes error bounds from the Fact 4 above.

  3. The special variant of zero-shift QR algorithm (the Demmel-Kahan algorithm) computes the singular values with high relative accuracy.

  4. The tridiagonal divide-and-conquer method, bisection and inverse iteration, and MRRR method can also be adapted for bidiagonal matrices.

  5. Zero shift QR algorithm for bidiagonal matrices is implemented in the LAPACK routine DBDSQR. It is also used in the function svdvals(). Divide-and-conquer algorithm for bidiagonal matrices is implemented in the LAPACK routine DBDSDC. However, this algorithm also calls zero-shift QR to compute singular values.

Examples


In [62]:
W,σ,Z=svd(B)


Out[62]:
SVD{Float64,Float64,Array{Float64,2}}
U factor:
5×5 Array{Float64,2}:
  0.471185   -0.699859    0.439019    0.274686   0.141392
  0.64595    -0.124171   -0.531454   -0.38439   -0.370313
 -0.548414   -0.554479   -0.0237941  -0.117773  -0.614295
  0.244482    0.430483    0.472585    0.270608  -0.677014
 -0.0143028  -0.0448941  -0.548558    0.830478  -0.0846805
singular values:
5-element Array{Float64,1}:
 22.71443872538939
 19.439175375504053
 14.741722158615167
 13.981197372822301
  7.374312302263113
Vt factor:
5×5 Array{Float64,2}:
 -0.357493  -0.584605  -0.605576   0.403968  -0.0228886
  0.620456   0.374817  -0.321418   0.606182  -0.0614842
 -0.513232   0.305443   0.333505   0.455503  -0.569728
 -0.338587   0.256748   0.157415   0.354245   0.818029
 -0.330432   0.598778  -0.627665  -0.369298  -0.0439948

In [65]:
@which svd(B)





In [66]:
σ₁=svdvals(B)


Out[66]:
5-element Array{Float64,1}:
 22.7144387253894
 19.43917537550406
 14.74172215861516
 13.9811973728223
  7.374312302263106

In [67]:
@which svdvals(B)


Out[67]:
svdvals(A::AbstractArray{#s664,2} where #s664<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64}) in LinearAlgebra at C:\Users\Ivan_Slapnicar\AppData\Local\Programs\Julia\Julia-1.4.1\share\julia\stdlib\v1.4\LinearAlgebra\src\svd.jl:194

In [68]:
σ-σ₁


Out[68]:
5-element Array{Float64,1}:
 -1.0658141036401503e-14
 -7.105427357601002e-15
  7.105427357601002e-15
  1.7763568394002505e-15
  7.105427357601002e-15

In [69]:
?LAPACK.bdsqr!


Out[69]:
bdsqr!(uplo, d, e_, Vt, U, C) -> (d, Vt, U, C)

Computes the singular value decomposition of a bidiagonal matrix with d on the diagonal and e_ on the off-diagonal. If uplo = U, e_ is the superdiagonal. If uplo = L, e_ is the subdiagonal. Can optionally also compute the product Q' * C.

Returns the singular values in d, and the matrix C overwritten with Q' * C.


In [70]:
B.dv


Out[70]:
5-element Array{Float64,1}:
 -17.233687939614086
 -14.890247384185125
  13.475190619553185
  13.673356118350501
  14.193966697750746

In [71]:
one(B)


Out[71]:
5×5 Bidiagonal{Float64,Array{Float64,1}}:
 1.0  0.0   ⋅    ⋅    ⋅ 
  ⋅   1.0  0.0   ⋅    ⋅ 
  ⋅    ⋅   1.0  0.0   ⋅ 
  ⋅    ⋅    ⋅   1.0  0.0
  ⋅    ⋅    ⋅    ⋅   1.0

In [77]:
BV=Matrix(one(B))
BU=Matrix(one(B))
BC=Matrix(one(B))
σ₂,Z₂,W₂,C = LAPACK.bdsqr!('U',copy(B.dv),copy(B.ev),BV,BU,BC)


Out[77]:
([22.71443872538939, 19.439175375504053, 14.741722158615167, 13.981197372822301, 7.374312302263113], [-0.3574931844947342 -0.5846050242892535 … 0.4039677234405411 -0.022888581149161945; 0.6204560208800932 0.3748165122085779 … 0.6061822062385038 -0.06148419671564419; … ; -0.33858691851450035 0.2567475263905849 … 0.3542448598305683 0.8180285726822231; -0.33043169301926517 0.5987778619961728 … -0.36929755632795147 -0.0439947614449794], [0.47118510340925945 -0.6998591041532833 … 0.2746858682943597 0.14139205185957202; 0.6459497679561967 -0.1241711534217206 … -0.38439013781932746 -0.3703125602362496; … ; 0.2444815803934529 0.4304832156808261 … 0.2706075571293588 -0.6770139926373585; -0.014302786105246401 -0.044894118385265115 … 0.8304775341367058 -0.08468046283229448], [0.47118510340925945 0.6459497679561967 … 0.2444815803934529 -0.014302786105246401; -0.6998591041532833 -0.1241711534217206 … 0.4304832156808261 -0.044894118385265115; … ; 0.2746858682943597 -0.38439013781932746 … 0.2706075571293588 0.8304775341367058; 0.14139205185957202 -0.3703125602362496 … -0.6770139926373585 -0.08468046283229448])

In [78]:
W₂'*B*Z₂'


Out[78]:
5×5 Array{Float64,2}:
 22.7144        2.60902e-15   1.44329e-15   9.4369e-16   -4.00374e-15
  2.95597e-15  19.4392       -4.55191e-15  -2.33147e-15   2.90046e-15
  4.96825e-15   6.32827e-15  14.7417       -2.66454e-15   5.16254e-15
  1.88738e-15   2.33147e-15   8.88178e-16  13.9812        6.66134e-16
  4.82253e-16   2.86229e-15   2.52576e-15   2.10942e-15   7.37431

In [79]:
?LAPACK.bdsdc!


Out[79]:
bdsdc!(uplo, compq, d, e_) -> (d, e, u, vt, q, iq)

Computes the singular value decomposition of a bidiagonal matrix with d on the diagonal and e_ on the off-diagonal using a divide and conqueq method. If uplo = U, e_ is the superdiagonal. If uplo = L, e_ is the subdiagonal. If compq = N, only the singular values are found. If compq = I, the singular values and vectors are found. If compq = P, the singular values and vectors are found in compact form. Only works for real types.

Returns the singular values in d, and if compq = P, the compact singular vectors in iq.


In [80]:
σ₃,ee,W₃,Z₃,rest=LAPACK.bdsdc!('U','I',copy(B.dv),copy(B.ev))


Out[80]:
([22.71443872538939, 19.439175375504053, 14.741722158615167, 13.981197372822301, 7.374312302263113], [0.0, 0.0, 0.0, 0.0], [0.47118510340925945 -0.6998591041532833 … 0.2746858682943597 0.14139205185957202; 0.6459497679561967 -0.1241711534217206 … -0.38439013781932746 -0.3703125602362496; … ; 0.2444815803934529 0.4304832156808261 … 0.2706075571293588 -0.6770139926373585; -0.014302786105246401 -0.044894118385265115 … 0.8304775341367058 -0.08468046283229448], [-0.3574931844947342 -0.5846050242892535 … 0.4039677234405411 -0.022888581149161945; 0.6204560208800932 0.3748165122085779 … 0.6061822062385038 -0.06148419671564419; … ; -0.33858691851450035 0.2567475263905849 … 0.3542448598305683 0.8180285726822231; -0.33043169301926517 0.5987778619961728 … -0.36929755632795147 -0.0439947614449794], [1.799709845e-315], [27174])

In [81]:
W₃'*B*Z₃'


Out[81]:
5×5 Array{Float64,2}:
 22.7144        2.60902e-15   1.44329e-15   9.4369e-16   -4.00374e-15
  2.95597e-15  19.4392       -4.55191e-15  -2.33147e-15   2.90046e-15
  4.96825e-15   6.32827e-15  14.7417       -2.66454e-15   5.16254e-15
  1.88738e-15   2.33147e-15   8.88178e-16  13.9812        6.66134e-16
  4.82253e-16   2.86229e-15   2.52576e-15   2.10942e-15   7.37431

Functions svd(), LAPACK.bdsqr!() and LAPACK.bdsdc!() use the same algorithm to compute singular values.


In [82]:
[σ₃-σ₂ σ₃-σ]


Out[82]:
5×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0

Let us compute some timings. We observe $O(n^2)$ operations.


In [84]:
n=1000
Abig=Bidiagonal(rand(n),rand(n-1),'U')
Bbig=Bidiagonal(rand(2*n),rand(2*n-1),'U')
@time svdvals(Abig);
@time svdvals(Bbig);
@time LAPACK.bdsdc!('U','N',copy(Abig.dv),copy(Abig.ev));
@time svd(Abig);
@time svd(Bbig);


  0.030342 seconds (11 allocations: 141.531 KiB)
  0.099849 seconds (11 allocations: 282.156 KiB)
  0.031539 seconds (11 allocations: 141.531 KiB)
  0.066507 seconds (16 allocations: 38.255 MiB)
  0.403456 seconds (16 allocations: 152.802 MiB, 34.88% gc time)

QR method

Final algorithm is obtained by combining bidiagonalization and bidiagonal SVD methods. Standard method is implemented in the LAPACK routine DGESVD. Divide-and-conquer method is implemented in the LAPACK routine DGESDD.

The functions svd(), svdvals(), and svdvecs() use DGESDD. Wrappers for DGESVD and DGESDD give more control about output of eigenvectors.


In [85]:
# The built-in algorithm
U,σA,V=svd(A)


Out[85]:
SVD{Float64,Float64,Array{Float64,2}}
U factor:
8×5 Array{Float64,2}:
  0.149592      0.697084    -0.454385    0.44584    -0.1969
  0.224464     -0.00315252   0.204138    0.411853    0.675146
  0.536028     -0.22465      0.0115001  -0.167755   -0.145731
  0.000933696  -0.55898     -0.422169    0.557899   -0.0806194
 -0.440475     -0.136885    -0.549489   -0.220845    0.463171
  0.619219     -0.0774185   -0.165031   -0.0874047   0.287435
  0.244406      0.0937872   -0.479637   -0.478975    0.0175803
  0.0522163    -0.342944    -0.119488    0.079884   -0.424519
singular values:
5-element Array{Float64,1}:
 22.71443872538939
 19.439175375504053
 14.741722158615167
 13.981197372822301
  7.374312302263113
Vt factor:
5×5 Array{Float64,2}:
 -0.357493   0.21599     0.0555667  -0.535823    0.731678
  0.620456  -0.0392581  -0.651222   -0.432618    0.0473813
 -0.513232   0.422882   -0.689941    0.248674   -0.14109
 -0.338587  -0.876721   -0.296922    0.0582911   0.158613
 -0.330432  -0.0658623   0.0929526  -0.678601   -0.646017

In [88]:
Sv=svd(A)


Out[88]:
SVD{Float64,Float64,Array{Float64,2}}
U factor:
8×5 Array{Float64,2}:
  0.149592      0.697084    -0.454385    0.44584    -0.1969
  0.224464     -0.00315252   0.204138    0.411853    0.675146
  0.536028     -0.22465      0.0115001  -0.167755   -0.145731
  0.000933696  -0.55898     -0.422169    0.557899   -0.0806194
 -0.440475     -0.136885    -0.549489   -0.220845    0.463171
  0.619219     -0.0774185   -0.165031   -0.0874047   0.287435
  0.244406      0.0937872   -0.479637   -0.478975    0.0175803
  0.0522163    -0.342944    -0.119488    0.079884   -0.424519
singular values:
5-element Array{Float64,1}:
 22.71443872538939
 19.439175375504053
 14.741722158615167
 13.981197372822301
  7.374312302263113
Vt factor:
5×5 Array{Float64,2}:
 -0.357493   0.21599     0.0555667  -0.535823    0.731678
  0.620456  -0.0392581  -0.651222   -0.432618    0.0473813
 -0.513232   0.422882   -0.689941    0.248674   -0.14109
 -0.338587  -0.876721   -0.296922    0.0582911   0.158613
 -0.330432  -0.0658623   0.0929526  -0.678601   -0.646017

In [92]:
Sv.Vt


Out[92]:
5×5 Array{Float64,2}:
 -0.357493   0.21599     0.0555667  -0.535823    0.731678
  0.620456  -0.0392581  -0.651222   -0.432618    0.0473813
 -0.513232   0.422882   -0.689941    0.248674   -0.14109
 -0.338587  -0.876721   -0.296922    0.0582911   0.158613
 -0.330432  -0.0658623   0.0929526  -0.678601   -0.646017

In [93]:
# With our building blocks
U₁=X*W
V₁=Y*Z
U₁'*A*V₁


Out[93]:
5×5 Array{Float64,2}:
 22.7144        5.25173e-15   2.61275e-16   2.45241e-15  -2.48326e-15
  9.82975e-16  19.4392       -2.94597e-15  -3.08156e-15   4.35454e-15
  2.75209e-15   7.27416e-15  14.7417        8.70922e-17   5.12142e-15
  2.117e-15    -6.18948e-16  -1.40526e-16  13.9812        2.83885e-16
  4.77623e-15   4.25303e-15   1.12624e-15   2.92931e-15   7.37431

In [94]:
?LAPACK.gesvd!


Out[94]:
gesvd!(jobu, jobvt, A) -> (U, S, VT)

Finds the singular value decomposition of A, A = U * S * V'. If jobu = A, all the columns of U are computed. If jobvt = A all the rows of V' are computed. If jobu = N, no columns of U are computed. If jobvt = N no rows of V' are computed. If jobu = O, A is overwritten with the columns of (thin) U. If jobvt = O, A is overwritten with the rows of (thin) V'. If jobu = S, the columns of (thin) U are computed and returned separately. If jobvt = S the rows of (thin) V' are computed and returned separately. jobu and jobvt can't both be O.

Returns U, S, and Vt, where S are the singular values of A.


In [95]:
# DGESVD
LAPACK.gesvd!('A','A',copy(A))


Out[95]:
([0.14959169523194776 0.697083564565451 … -0.18855944807540953 0.06394902354873996; 0.22446405928471988 -0.003152521629016692 … 0.2246752303874581 0.2771850637830364; … ; 0.24440607421233626 0.09378723173362023 … 0.6463125113361845 -0.026618264127683304; 0.052216273049930303 -0.34294390975336814 … 0.013928748176334986 0.8189122696488763], [22.71443872538939, 19.439175375504053, 14.741722158615156, 13.9811973728223, 7.374312302263106], [-0.35749318449473405 0.21599040938785102 … -0.5358232301323785 0.731677919060666; 0.6204560208800932 -0.03925807784782213 … -0.4326176228012191 0.04738127856925271; … ; -0.3385869185145026 -0.8767212111366486 … 0.05829113338948243 0.15861294372785595; 0.33043169301926517 0.06586234793778371 … 0.6786005823725091 0.646017120670051])

In [96]:
?LAPACK.gesdd!


Out[96]:
gesdd!(job, A) -> (U, S, VT)

Finds the singular value decomposition of A, A = U * S * V', using a divide and conquer approach. If job = A, all the columns of U and the rows of V' are computed. If job = N, no columns of U or rows of V' are computed. If job = O, A is overwritten with the columns of (thin) U and the rows of (thin) V'. If job = S, the columns of (thin) U and the rows of (thin) V' are computed and returned separately.


In [97]:
LAPACK.gesdd!('N',copy(A))


Out[97]:
(Array{Float64}(undef,8,0), [22.7144387253894, 19.43917537550406, 14.74172215861516, 13.9811973728223, 7.374312302263106], Array{Float64}(undef,5,0))

Let us perform some timings. We observe $O(n^3)$ operations.


In [98]:
n=1000
Abig=rand(n,n)
Bbig=rand(2*n,2*n)
@time Ubig,σbig,Vbig=svd(Abig);
@time svd(Bbig);
@time LAPACK.gesvd!('A','A',copy(Abig));
@time LAPACK.gesdd!('A',copy(Abig));
@time LAPACK.gesdd!('A',copy(Bbig));


  1.144909 seconds (17 allocations: 45.899 MiB, 1.67% gc time)
  7.258646 seconds (13 allocations: 183.350 MiB, 0.52% gc time)
  6.398566 seconds (10 allocations: 23.408 MiB, 2.25% gc time)
  1.123545 seconds (12 allocations: 45.899 MiB)
  7.074224 seconds (12 allocations: 183.350 MiB, 1.29% gc time)

In [99]:
# Residual
norm(Abig*Vbig-Ubig*Diagonal(σbig))


Out[99]:
9.385834663948128e-13

In [ ]: