Algorithms: SVD (singular value decomposition)

The SVD of a matrix provides with a powerful tool to reduce the dimensionality of a dataspace. It is used in the DMRG (density matrix renormalization group) to perform the decimation of the wavefunction basis to keep it from growing as nodes/spins are added to the system. It is used in image compression, and in PCA (principle component analysis) to reduce the dataspace to the most relevant dimensions.

SVD (singular value decomposition) is a decomposition of a $n\times p$ matrix into $U$ , an $n\times n$ , $S$ a “diagonal” $n\times p$ , and $V$ , a $p\times p$ matrix $M=USV^T$ . Note that $MM^T$ is a square $n\times n$ . The columns of $U$ are the eigenvectors of $MM^T$ . $M^TM$ is a square $p\times p$ and the columns of $V$ are it’s eigenvectors. Both $MM^T$ and $M^TM$ are positive semi-definite, their eigenvalues are non-negative. A matrix $A$ is positive semi-definite if for any real vector $v$ : $v^TAv\ge 0$

$v^T(MM^T)v=|M^T v|^2\ge 0$

If $M$ is a square matrix then $M=U\Lambda U^T$ is its eigenvalue decomposition: $U$ is an orthogonal matrix made of the eigenvectors such that $MU=U\Lambda$ ie $Mu_j=\lambda_j u_j$ , where $u_j$ is a column vector (and $u_j^T$ is a row vector)

$M[u_1,u_2,\cdots, u_n]=[\lambda_1u_1,\lambda_2 u_2,\cdots, \lambda_n u_n]=[u_1,u_2,\cdots, u_n]\Lambda$

The SVD generalizes this no non-square arrays. Another useful decomposition:

$M=\sum_js_j u_jv_j^T, \qquad 1=\sum_j v_jv_j^T=\sum_j P_j$

(the eigendecomposition and decomposition of unity into projection operators)

$u_ju_j^T\cdot u_k=u_j\delta_{jk}=u_k$

$M^TM=\Big\sum_js_j v_ju_j^T(\Big)\Big(\sum_k s_k u_k v_k^T\Big)=\sum_{jk}s_j s_k \delta_{jk}v_j v_k^T=\sum_j s_j^2 v_j v^T_j$

Note that $(M^T M)^T=M^T M$ so $M^T M$ is symmetric and real hence hermitian, so $\{v_k\}$ form a basis of orthogonal vectors, a property which we have used.
The Lanczos algorithm works very well to compute the largest singular value $S[1,1]$ ,
consider that
an arbitrary vector $w$ has an expansion

$w=\sum_j b_j v_j$

$(M^T M)^n w=\sum_j b_j v_j=\sum_j b_j s_j^{2n} v_j\rightarrow b_{1}s_j^{2n}n v_1+\cdots$

where $1$ denotes the index of the maximal eigenvalue/vector of $M^TM$

${(M^T M)^n w\over |(M^T M)^n w|}=v_1+\sum_{j\ne 1}({s_j\over s_1})^{2n} {b_j\over b_1} v_j\rightarrow v_1$

as $n\rightarrow \infty$ . Then since $Mv_1=s_1 u_1$

$M{(M^T M)^n w\over |(M^T M)^n w|}\rightarrow s_1 u_1$

$\bullet$ Extract the most relevant (largest) singular value/generalized eigenvalue $s_1$ by the Lanczos algorithm , $\bullet$ then subtract its contribution from $M$ to get

$M'=M-s_1 u_1 v_1^T, \qquad M'^T M'=M^TM-s_1^2 v_1 v_1^T$

(since $s_1 v_1 u_1^T M=s_1^2 v_1 v_1^T$ and $M^T s_1 u_1 v_1^T=s_1^2 v_1 u_1^T u_1 v_1^T=s_1^2 v_1 v_1^T$ ). Next apply the Lanczos algorithm to $M'1$ to get $s_2$ , the second most relevant singular value.
$\bullet$ repeat $n$ times in total to get the $n$ most relevant singular values.

Let’s have it in julia

""" returns U,S,V  for maximal singular value"""
maxsingularval=function(M,prec)
    n,p=size(M)
    v=ones(p)
    v=v ./ sqrt(sum(v .* v))
    w=M'*M*v
    val=sqrt(sum(M*w .* M*w)/sum(w .* w))
    iters=0
    while true
        w=M'*M*w
        w=w ./ sqrt(sum(w .* w))
        newval=sqrt(sum(M*w .* M*w)/sum(w .* w))
        if(abs(newval-val) < prec)
            break
        end
        val=newval
    end
return val, M*w ./ sqrt(sum(M*w .* M*w)), w
end

""" svd decomposition, NN<= min(ncol(M),nrow(M))"""
mysvd=function(M,NN,prec)
   n,p=size(M)
   U=zeros(n,n)
   S=zeros(n,p)
   V=zeros(p,p)
S[1,1], U[:,1],V[:,1]=maxsingularval(M,prec)
Mp=M
for j in 2:NN
Mp=Mp - S[j-1,j-1]*U[:,j-1]*V[:,j-1]'
S[j,j], U[:,j],V[:,j]=maxsingularval(Mp,prec)
end
return S,U,V'
end

Let’s test it with a $7\times 9$ random array. We can use our decomposition to create a pseudoinverse of a nonsquare array and check its integrity.

psinv=function(M,prec)
   n,p=size(M)
   k=minimum([n,p])
   ms,mu,mv=mysvd(M,k,prec)
for j in 1:k
if(ms[j,j]!=0.0)
ms[j,j]=1.0/ms[j,j]
end
end
mv'*ms'*mu'
end

m=random(7,9)
mS,mU,mV=mysvd(m,7,0.00001)
mU*mS*mV # should equal m, and it is VERY close
sum(abs.(m*psinv(m,0.0000001))) # should be a 7x9 with ones down the "diagonal" so sum
# of all entries should be 7, it returns 7.003860424

Next: applications.

Home 2.0

Algorithms: SVD (singular value decomposition)

Related Posts

Algorithms: SVD (singular value decomposition)

Published by Jeff on December 20, 2021December 20, 2021

Related Posts

Protected: AI/ML-3. Flux loss functions

ML/algorithms: PCA (principal component analysis)

Protected: AI/ML-6. ConvNets