% demo low-rank approximation of an image matrix. Barnett 11/11/16

clear;
A = double(imread('mystery.jpg'));   % if you want to read JPG directly
% A = load('image.mat');                 % one already in Matlab format
size(A)

[U S V] = svd(A);

sigma = diag(S);  % singular vals

figure; plot(sigma,'+-'); title('singular values {\sigma_j}');
xlabel('j'); set(gca','yscale','log')

vtot = norm(A(:))^2    % sum of squares of entries of A = sum of sigma_j^2
% aka "total variance"
vtot = sum(sigma.^2)   % another way to get 

klist = [1 2 3 4 5 7 10 20 50 100];  % ranks to explore

for k=klist
  % use just first k singular values & vectors to build best rank-k approx...
  Ak = U(:,1:k) * ( S(1:k,1:k) * V(:,1:k)' );
  vk = sum(sigma(1:k).^2);   % variance explained by just sigma_1,..,sigma_k.
  figure; imagesc(Ak); axis equal tight; colormap(gray(256));
  title(sprintf('rank k = %d (%.3g of variance)',k,vk/vtot)); drawnow
  pause
end

figure; imagesc(A); axis equal tight; colormap(gray(256));
title('exact'); drawnow

sum(sigma(5:end).^2) / vtot   % frac of variance explained by best rank-4 approx

figure; n=numel(klist);
for i=1:n, k=klist(i);
  % use just first k singular values & vectors to build best rank-k approx...
  Ak = U(:,1:k) * ( S(1:k,1:k) * V(:,1:k)' );
  vk = sum(sigma(1:k).^2);   % variance explained by just sigma_1,..,sigma_k.
  subplot(2,n,i); imagesc(Ak); axis equal tight off; colormap(gray(256));
  title(sprintf('rank k=%d (%.3g of var)',k,vk/vtot)); c = caxis;
  subplot(2,n,i+n); imagesc(Ak-A); caxis(c); % use same grayscale
  axis equal tight off; colormap(gray(256));
end

% Warning: this is not how image compression is done in real life!
% One reason is that it's highly rotation-dependent.
% There are better ways
% (eg, based on sparsification of Fourier series in blocks)
