% 
% The prototype randomized SVD algorithm 
%
% Note that the code is not rewritten to optimize the cpu time
%
%-----------------------------------
%
% small matrix example 
%n = 10; 
%m = 8; 
%s = 5; 
%k = 2; 
%A = randn(m,n);  

% image matrix example 
load durer % Albrecht Durer (1471-1528)'s Melencolia
subplot(2,2,1) 
colormap(gray(128));
image(X)
axis image
axis off
title('Melencolia, Durer (1514)')

load detail % the detail of the magic square
subplot(2,2,2) 
colormap(gray(64));
image(X)
axis image
axis off
title('Magic Sqaure')
A = X;    % 359-by-371

%---------------------------------

[n,m]=size(A); 
s = 50;   % sketching size
k = 40;   % the rank of approximation 

% --- rank-k truncated SVD 
[U,S,V]=svd(A); 
Akopt = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';  

% --- randomized algorithm for rank-k approximation

Omega = randn(m,s); 
C = A*Omega;     % sampling 
Q = orth(C); 
[Ua,Sa,Va] = svd(Q'*A); 
Ak = (Q*Ua(:,1:k))*Sa(1:k,1:k)*Va(:,1:k)';  

%----------------------

diff_with_opt = norm(Akopt - Ak)/norm(Akopt)

subplot(2,2,3) 
colormap(gray(64)); 
image(Akopt) 
axis image
axis off
title('Optimal rank-k')

subplot(2,2,4) 
colormap(gray(64)); 
image(Ak) 
axis image
axis off
title('Randomized rank-k')
