% 
% The prototype randomized SVD algorithm 
% with acceleration (the power method) 
%
% 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
figure(1) 
subplot(2,2,1) 
colormap(gray(64));
image(X)
axis image
axis off
title('Magic Sqaure in Melencolia, Durer 1514')
A = X;    % 359-by-371

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

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

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

figure(1) 
subplot(2,2,2) 
colormap(gray(64)); 
image(Akopt) 
axis image
axis off
title('Optimal rank-k')
figure(2) 
semilogy(diag(S(1:k+5,1:k+5)),'r+'); 
hold on

% --- the prototype randomized algorithm 

randn('state',31456);
Omega = randn(m,s); 

C = A*Omega;    % sketching of A 
Q = orth(C); 
[Ua,Sa,Va] = svd(Q'*A); 
Ak = (Q*Ua(:,1:k))*Sa(1:k,1:k)*Va(:,1:k)';  

figure(1) 
subplot(2,2,3) 
colormap(gray(64)); 
image(Ak) 
axis image
axis off
title('Randomized rank-k')
figure(2) 
semilogy(diag(Sa(1:k,1:k)),'bo'); 

% --- the accelerated prototype randomized algorithm 

C = A*Omega; 
Q = orth(C); 
for i = 1:q 
   C = A'*Q; 
   Q = orth(C); 
   C = A*Q; 
   Q = orth(C); 
end 
[Ua2,Sa2,Va2] = svd(Q'*A); 
Ak2 = (Q*Ua2(:,1:k))*Sa2(1:k,1:k)*Va2(:,1:k)';  

figure(1) 
subplot(2,2,4) 
colormap(gray(64)); 
image(Ak) 
axis image
axis off
title('Accelerated randomized rank-k')
figure(2) 
semilogy(diag(Sa2(1:k,1:k)),'gx'); 
hold off

diff1 = norm(A - Akopt)/norm(A);
diff2 = norm(A - Ak)/norm(A);
diff3 = norm(A - Ak2)/norm(A);
diff_with_org = [diff1 diff2 diff3 ]

diff4 = norm(Akopt - Ak)/norm(Akopt);
diff5 = norm(Akopt - Ak2)/norm(Akopt);
diff_with_opt = [diff4 diff5]
