%% 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')

% 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)

%
% 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')

%
% 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;

% Copyright 2012 The MathWorks, Inc.