This script demonstrates how BBSVDS can be used to compute singular triplets. To start it directly, use:
playshow bbdemo_svd2
First, we generate a small matrix with singular values that decrase as a hyperbole. Matrices of this type occur in data-modeling problems such as latent semantic indexing.
To make it a bit challenging, we permute the matrix and apply a sparse random Householder reflector on each side.
m=1600; n=500; sig0=.1./(.1+10*(0:499)'); A=spdiags(sig0,0,m,n); A=A(randperm(m),randperm(n)); qu=2*rand(m,1)-1; qu=qu.*(rand(m,1)>.95); qv=2*rand(n,1)-1; qv=qv.*(rand(n,1)>.95); qu=sparse(qu/norm(qu)); A=A-2*qu*(qu'*A); qv=sparse(qv/norm(qv)); A=A-2*(A*qv)*qv'; semilogy(sig0,'.b'); hold on; grid on; title('Typical data-modeling problem'); xlabel('index'); ylabel('Singular value');
Now compute the leading singular triplets of A with a call to BBSVDS. For later reference we wrap A in a counter and time the operation.
NOTE: triplets are not guaranteed to be "the" leading triplets. A few triplets may be paired with wrong exact triplets, so a few values may may show an unexpectedly large "absolute error".
Ac=bbcount(A); opt1=bbsvdopt(Ac,400,true); t=cputime; [U1,S1,V1]=bbsvds(Ac,400,opt1); T1=cputime-t; cstat1=bbcountstat(Ac,'bbsvds_count'); bbcountreset(Ac); sig1=diag(S1); l1=length(sig1); semilogy(abs(sig0(1:l1)-sig1),'.g'); title('Accuracy of singular values'); ylabel('Singular value and absolute error');
400 triplets converged.
We may limit the amount of memory by reducing the amount of deflation vectors (see BBDEMO_SVD1 for details).
opt2=opt1; opt2.dlead=50; opt2.dtail=50; t=cputime; [U2,S2,V2]=bbsvds(Ac,400,opt2); T2=cputime-t; cstat2=bbcountstat(Ac,'bbsvds_count'); bbcountreset(Ac); sig2=diag(S2); l2=length(sig2); semilogy(abs(sig0(1:l2)-sig2), '.r');
Total number of converged triplets: 242 Total number of converged triplets: 313 Halting: triplets stopped converging. Increase dlead or rmaxi if you need more triplets. 337 triplets converged.
As an alternative, we may resort to one-sided deflation. As we will show later this affects the accuracy and orthogonality of the singular triplets.
opt3=opt1; opt3.dsides=1; t=cputime; [U3,S3,V3]=bbsvds(Ac,400,opt3); T3=cputime-t; cstat3=bbcountstat(Ac,'bbsvds_count'); bbcountreset(Ac); sig3=diag(S3); l3=length(sig3); semilogy(abs(sig0(1:l3)-sig3),'.m');
Warning: One-sided orthogonalization forced. The singular triplets may not be accurate. 407 triplets converged.
Here is a comparison of the different methods. The estimated memory consumption does not include the returned matrices U and V.
mem1=bbsvdmem(Ac,opt1); mem2=bbsvdmem(Ac,opt2); mem3=bbsvdmem(Ac,opt3); bbclear Ac; fprintf('Converged: %5i, %5i, %5i\n', l1, l2, l3); fprintf('Computation times: %5.2f, %5.2f, %5.2f (s)\n', T1, T2, T3); fprintf('Memory usage: %5.2f, %5.2f, %5.2f (Mb)\n', mem1, mem2, mem3);
Converged: 400, 337, 400 Computation times: 6.95, 25.44, 5.67 (s) Memory usage: 16.66, 7.05, 7.39 (Mb)
Here we illustrate how much work each procedure used in terms of matrix-vector products: Green: Simple call to BBSVDS Red: Deflation limited to 100 triplets Black: One-sided deflation
hold off; plot([cstat1.mm],'.g'); hold on; grid on plot([cstat2.mm],'.r'); plot([cstat3.mm],'.k'); hold off; title('Effort of the SVD procedures'); xlabel('Converged triplets'); ylabel('Matrix-vector products');
We saw that both one- and two-sided deflation yielded accurate singular values. Another major usage for one-sided deflation is in low-rank approximations. According to the Schmidt-Mirsky theorem, the matrix B=U(:,1:N)*S(1:N,1:N)*V(:,1:N)' minimizes the norm of A-B among any matrix of rank-N.
Let us compare the difference between the one- and two-sided deflation for this purpose. We disregard the last few triplets, since there may be a few missed triplets in the end. That does not tell anything about the accuracy of the returned triplets.
N=390; B1=U1(:,1:N)*S1(1:N,1:N)*V1(:,1:N)'; B3=U3(:,1:N)*S3(1:N,1:N)*V3(:,1:N)'; fprintf('Frobenius norm of the residual:\n'); fprintf(' Theoretical minimum: %.8e\n', norm(sig0(N+1:end))); fprintf(' Two-sided deflation: %.8e\n', norm(A-B1,'fro')); fprintf(' One-sided deflation: %.8e\n', norm(A-B3,'fro'));
Frobenius norm of the residual: Theoretical minimum: 2.37774161e-04 Two-sided deflation: 2.37774161e-04 One-sided deflation: 2.37774161e-04
The downside of one-sided orthogonalization is that the individual triplets are less accurate. First let us compare the orthogonality levels.
In general, the orthogonality of the long singular space (in this case U) will only be eps*(sig/NORM(A)).
Top: two-sided deflation, Bottom: One-sided deflation. Left: Left singular vectors, Right: Right singular vectors.
MU1=abs(U1'*U1-eye(l1)); MV1=abs(V1'*V1-eye(l1)); MU3=abs(U3'*U3-eye(l3)); MV3=abs(V3'*V3-eye(l3)); o=ones(400,10); O=ones(10,810); G=[[MU1, o, MV1] ; O ; [MU3, o, MV3]]; imagesc(log10(realmin+G),[-16,0]); title('Cosine of angle between triplets'); xlabel('index 1-400'); ylabel('index 1-400'); set(gca,'XTick',[],'YTick',[]); hdl=get(colorbar,'ylabel'); set(hdl,'string','Log-magnitude (base 10)');
The norm of the residuals also increase for triplets in the long-space as the singular value decrease: Blue: residual for the left singular triplets. Green: residual for the right singular triplets.
EU3=A*V3-U3*S3; eu3=sqrt(sum(EU3.*EU3))'; EV3=A'*U3-V3*S3; ev3=sqrt(sum(EV3.*EV3))'; semilogy(eu3,'.b'); hold on; grid on; semilogy(ev3,'.g'); title('Accuracy of singular triplets'); xlabel('index'); ylabel('Norm of residual');
