%% Approximating a matrix with a low-rank matrix.
% Given a matrix A, we can approximate it with a matrix of lower rank using
% the SVD.  Let's say our original A is full rank -- for instance, a
% photograph of McGraw Tower.
A = double(imread('tower.png'))/255.0;
figure(1); subplot(1,2,1);
imshow(A);
title('original image');
%% Computing the approximation
% First, we compute the SVD of A.  We'll use the "thin" version -- though
% it hardly matters with a matrix this nearly square.
[U,S,V] = svd(A,0);
s = diag(S);
n = length(s);
figure(2); semilogy(s, '-o');
%% Reconstrucing A exactly
% Now U(:,1:r) * S(1:r,1:r) * U(1:r,:) is a rank-r matrix that approximates
% A. Note that when r = n we get back the exact A (well, up to rounding
% errors).
figure(1); subplot(1,2,2);
imshow(U * S * V');
title(sprintf('rank = %d', n));
%% Lower-rank approximations
% The orthogonality of U and V (the full U and V, not the ones we just
% computed) guarantees that setting parts of S to zero will induce the same
% F-norm error in the approximation.  So the approximation error is exactly
% the norm of the unused part of diag(S).
figure(1); subplot(1,2,2);
M = norm(A, 'fro');
for k = [1 2 3 10 50 100 200]
    A_apx = U(:,1:k) * S(1:k,1:k) * V(:,1:k)';
    imshow(A_apx);
    title(sprintf('rank = %d: err = %.1f%% | %.1f%%', k, ...
          100*norm(A - A_apx, 'fro')/M, 100*norm(s(k+1:n))/M));
    pause
end
%% Disclaimer
% This is not really a useful application of low-rank approximation,
% because an image is not usually close to low rank, viewed in terms of its
% natural rows and colums. (Indeed, to get a usable approximation we
% needed to include more than half the data, whereas a real image
% compression method (JPEG) can achieve similar quality at 10:1 or 20:1
% compression.) It's just a nice visual demo of the error properties.