% November 16, 2020, P.G. Martinsson
%
% This code generates the plots showing the accuracy of randomized low rank
% approximation.
%
% DRIVER_single: Show the accuracy for a single instantiation.
%
% DRIVER_many: Show 100 error curves in one graph
%
% DRIVER_histogram: Generate histograms over the errors.
%
% DRIVER_histogram_oversamp: Use oversampling to improve the accuracy.
%
% In each file, the variable "flag_matrix" chooses a matrix.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function main_RSVD_accuracy

DRIVER_single
%DRIVER_many
%DRIVER_histogram
%DRIVER_histogram_oversamp

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_single

rng(7)

%%% Set a test matrix.
flag_matrix = 3;
if (flag_matrix == 1)
  m = 200;
  n = 160;
  A = LOCAL_fastdecay(m,n);
elseif (flag_matrix == 2)
  m = 200;
  n = 160;
  A = LOCAL_slowdecay(m,n);
elseif (flag_matrix == 3)
  m = 200;
  n = 160;
  A = LOCAL_S_curve(m,n,round(0.5*n),1e-4);
elseif (flag_matrix == 4)
  m = 200;
  n = 160;
  A = LOCAL_smooth(m,n);
elseif (flag_matrix == 5)
  m    = 300;
  n    = 300;
  diam = 40;
  A    = LOCAL_helmholtz(m,n,diam);
elseif (flag_matrix == 6)
  m = 200;
  n = 200;
  A = generate_BIE(n);
end
p  = min(m,n);
ss = svd(A);

kmax = 100;

%%%
G      = randn(n,kmax);
[Q0,~] = qr(A*G,  0);
[U0,~] = qr(A'*Q0,0);
[Q1,~] = qr(A*U0, 0);
[U1,~] = qr(A'*Q1,0);
[Q2,~] = qr(A*U1, 0);

err0 = zeros(1,kmax);
err1 = zeros(1,kmax);
err2 = zeros(1,kmax);
for k = 1:kmax
  err0(k) = norm(A - Q0(:,1:k)*(Q0(:,1:k)'*A));
  err1(k) = norm(A - Q1(:,1:k)*(Q1(:,1:k)'*A));
  err2(k) = norm(A - Q2(:,1:k)*(Q2(:,1:k)'*A));
end

%figure(1)
%LW = 'LineWidth';
%semilogy(1:p,ss,'k',...
%         1:(p-1),R11_err0,'b',...
%         2:p,R22_err0,'b--',LW,1.2), hold on
%semilogy(1:(p-1),R11_err1,'r',...
%				 2:p,R22_err1,'r--',LW,1.2)
%semilogy(1:(p-1),R11_err2,'g',...
%				 2:p,R22_err2,'g--',LW,1.2)
%semilogy(1:(p-1),R11_errs,'m',...
%				 2:p,R22_errs,'m--',LW,1.2)
%hold off
%legend({'singular values',...
	%'$\sigma_k(R_{11})$, CPQR',...
	%'$\sigma_1(R_{22})$, CPQR',...
	%'$\sigma_k(R_{11})$, PowerURV $q = 0$',...
	%'$\sigma_1(R_{22})$, PowerURV $q = 0$',...
	%'$\sigma_k(R_{11})$, PowerURV $q = 1$',...
	%'$\sigma_1(R_{22})$, PowerURV $q = 1$',...
	%'$\sigma_k(R_{11})$, PowerURV $q = 2$',...
	%'$\sigma_1(R_{22})$, PowerURV $q = 2$',...
	%},'Interpreter','LaTeX')

figure(1)
hold off
han = semilogy(1:kmax,err0,'r',...
               1:kmax,err1,'g',...
               1:kmax,err2,'b',...
               0:kmax,ss(1:(kmax+1)),'k',...
               'LineWidth',2);
legend('Exact SVD (optimal)',...
	   'RSVD q = 0',...
	   'RSVD q = 1',...
	   'RSVD q = 2')
axis tight
xlabel('$k$','Interpreter','LaTeX')
%ylabel('$\|A - Q_k\, Q_k^{*}\,A\|$','Interpreter','LaTeX')
ylabel('Spectral norm error')
title('Accuracy in rank k approximation','FontWeight','normal','FontAngle','italic')
legend([han(4),han(1),han(2),han(3)],...
       {'Exact SVD (optimal)',...
        'Basic RSVD (q=0)',...
        'RSVD with one step of power iteration (q=1)',...
        'RSVD with two steps of power iteration (q=2)'})
keyboard
%export_fig -trans -pdf fig_RSVD_accuracy_helm.pdf
%export_fig -trans -pdf fig_RSVD_accuracy_slow.pdf

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_many

rng(0)
nexp = 20;

%%% Set a test matrix.
flag_matrix = 5;
if (flag_matrix == 1)
  m = 200;
  n = 160;
  A = LOCAL_fastdecay(m,n);
elseif (flag_matrix == 2)
  m = 200;
  n = 160;
  A = LOCAL_slowdecay(m,n);
elseif (flag_matrix == 3)
  m = 200;
  n = 160;
  A = LOCAL_S_curve(m,n,round(0.5*n),1e-4);
elseif (flag_matrix == 4)
  m = 200;
  n = 160;
  A = LOCAL_smooth(m,n);
elseif (flag_matrix == 5)
  m    = 300;
  n    = 300;
  diam = 40;
  A    = LOCAL_helmholtz(m,n,diam);
elseif (flag_matrix == 6)
  m = 200;
  n = 200;
  A = generate_BIE(n);
end
ss = svd(A);

%kk   = 10:10:100;
kk   = 1:100;
ERR0 = zeros(length(kk),nexp);
ERR1 = zeros(length(kk),nexp);
ERR2 = zeros(length(kk),nexp);

for j = 1:nexp

  fprintf(1,'Experiment %3d of %3d\n',j,nexp)
  G      = randn(n,max(kk));
  [Q0,~] = qr(A*G,  0);
  [U0,~] = qr(A'*Q0,0);
  [Q1,~] = qr(A*U0, 0);
  [U1,~] = qr(A'*Q1,0);
  [Q2,~] = qr(A*U1, 0);

  for i = 1:length(kk)
    k         = kk(i);
    ERR0(i,j) = norm(A - Q0(:,1:k)*(Q0(:,1:k)'*A));
    ERR1(i,j) = norm(A - Q1(:,1:k)*(Q1(:,1:k)'*A));
    ERR2(i,j) = norm(A - Q2(:,1:k)*(Q2(:,1:k)'*A));
  end
end

mean_err0 = (1/nexp)*sum(ERR0,2);
mean_err1 = (1/nexp)*sum(ERR1,2);
mean_err2 = (1/nexp)*sum(ERR2,2);

figure(1)
hold off
han = semilogy(kk',ERR0,'r',...
               kk',ERR1,'g',...
               kk',ERR2,'b',...
               kk',mean_err0,'r',...
               kk',mean_err1,'g',...
               kk',mean_err2,'b',...
               0:max(kk),ss(1:(max(kk)+1)),'k',...
               'LineWidth',1);
set(han(0*nexp + (1:nexp)),'Color',[1,0.7,0.7])
set(han(1*nexp + (1:nexp)),'Color',[0.7,1,0.7])
set(han(2*nexp + (1:nexp)),'Color',[0.7,0.7,1])
axis tight
xlabel('$k$','Interpreter','LaTeX')
%ylabel('$\|A - Q_k\, Q_k^{*}\,A\|$','Interpreter','LaTeX')
ylabel('Spectral norm error')
%title(sprintf('Accuracy in rank k approximation - showing %d instantiations',nexp),...
%      'FontWeight','normal','FontAngle','italic')
legend([han(3*nexp+4),han(3*nexp+1),han(3*nexp+2),han(3*nexp+3)],...
       {'Exact SVD (optimal)',...
        'Basic RSVD (average of all runs)',...
        'RSVD with one step of power iteration (average of all runs)',...
        'RSVD with two steps of power iteration (average of all runs)'})
keyboard
%export_fig -trans -pdf fig_RSVD_accuracy_slow_many.pdf

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_histogram

rng(0)
nexp = 5000;

%%% Set a test matrix.
flag_matrix = 5;
if (flag_matrix == 1)
  m = 200;
  n = 160;
  A = LOCAL_fastdecay(m,n);
elseif (flag_matrix == 2)
  m = 200;
  n = 160;
  A = LOCAL_slowdecay(m,n);
elseif (flag_matrix == 3)
  m = 200;
  n = 160;
  A = LOCAL_S_curve(m,n,round(0.5*n),1e-4);
elseif (flag_matrix == 4)
  m = 200;
  n = 160;
  A = LOCAL_smooth(m,n);
elseif (flag_matrix == 5)
  m    = 300;
  n    = 300;
  diam = 40;
  A    = LOCAL_helmholtz(m,n,diam);
elseif (flag_matrix == 6)
  m = 200;
  n = 200;
  A = generate_BIE(n);
end
ss = svd(A);

k    = 50;
err0 = zeros(1,nexp);
err1 = zeros(1,nexp);
err2 = zeros(1,nexp);

for j = 1:nexp

  fprintf(1,'Experiment %3d of %3d\n',j,nexp)
  G      = randn(n,k);
  [Q0,~] = qr(A*G,  0);
  [U0,~] = qr(A'*Q0,0);
  [Q1,~] = qr(A*U0, 0);
  [U1,~] = qr(A'*Q1,0);
  [Q2,~] = qr(A*U1, 0);

  err0(j) = norm(A - Q0(:,1:k)*(Q0(:,1:k)'*A));
  err1(j) = norm(A - Q1(:,1:k)*(Q1(:,1:k)'*A));
  err2(j) = norm(A - Q2(:,1:k)*(Q2(:,1:k)'*A));
end

figure(1)
subplot(1,3,1)
hold off
plot(ss(k+1),0,'k.','MarkerSize',25)
hold on
histogram(err0)
subplot(1,3,2)
hold off
histogram(err1)
hold on
plot(ss(k+1),0,'k.','MarkerSize',25)
subplot(1,3,3)
hold off
histogram(err2)
hold on
plot(ss(k+1),0,'k.','MarkerSize',25)

figure(2)
subplot(1,1,1)
hold off
plot(ss(k+1),0,'k.','MarkerSize',25)
hold on
han0 = histogram(err0,10);
han1 = histogram(err1,10);
han2 = histogram(err2,10);
set(han0,'FaceColor','r')
set(han1,'FaceColor','g')
set(han2,'FaceColor','b')
export_fig -png -trans fig_RSVD_accuracy_histogram_fast.png

figure(3)
errmax = max([err0,err1,err2]);
xbins  = linspace(ss(k+1),errmax,50);
subplot(1,3,1)
hold off
plot(ss(k+1),0,'k.','MarkerSize',25)
hold on
histogram(err0,xbins)
subplot(1,3,2)
hold off
histogram(err1,xbins)
hold on
plot(ss(k+1),0,'k.','MarkerSize',25)
subplot(1,3,3)
hold off
histogram(err2,xbins)
hold on
plot(ss(k+1),0,'k.','MarkerSize',25)


keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_histogram_oversamp

rng(0)
nexp = 500;

%%% Set a test matrix.
flag_matrix = 2;
if (flag_matrix == 1)
  m = 200;
  n = 160;
  A = LOCAL_fastdecay(m,n);
elseif (flag_matrix == 2)
  m = 200;
  n = 160;
  A = LOCAL_slowdecay(m,n);
elseif (flag_matrix == 3)
  m = 200;
  n = 160;
  A = LOCAL_S_curve(m,n,round(0.5*n),1e-4);
elseif (flag_matrix == 4)
  m = 200;
  n = 160;
  A = LOCAL_smooth(m,n);
elseif (flag_matrix == 5)
  m    = 300;
  n    = 300;
  diam = 40;
  A    = LOCAL_helmholtz(m,n,diam);
elseif (flag_matrix == 6)
  m = 200;
  n = 200;
  A = generate_BIE(n);
end
ss = svd(A);

k    = 50;
p    = 10;
err0 = zeros(1,nexp);
err1 = zeros(1,nexp);
err2 = zeros(1,nexp);

for j = 1:nexp

  fprintf(1,'Experiment %3d of %3d\n',j,nexp)
  G      = randn(n,k+p);
  [Q0,~] = qr(A*G,  0);
  [U0,~] = qr(A'*Q0,0);
  [Q1,~] = qr(A*U0, 0);
  [U1,~] = qr(A'*Q1,0);
  [Q2,~] = qr(A*U1, 0);

  [Uhat,D,V] = svd(Q0'*A,'econ');
  WW0 = Q0*Uhat;
  [Uhat,D,V] = svd(Q1'*A,'econ');
  WW1 = Q1*Uhat;
  [Uhat,D,V] = svd(Q2'*A,'econ');
  WW2 = Q2*Uhat;
  
  err0(j) = norm(A - WW0(:,1:k)*(WW0(:,1:k)'*A));
  err1(j) = norm(A - WW1(:,1:k)*(WW1(:,1:k)'*A));
  err2(j) = norm(A - WW2(:,1:k)*(WW2(:,1:k)'*A));
end

figure(1)
subplot(1,3,1)
hold off
plot(ss(k+1),0,'k.','MarkerSize',25)
hold on
histogram(err0)
subplot(1,3,2)
hold off
histogram(err1)
hold on
plot(ss(k+1),0,'k.','MarkerSize',25)
subplot(1,3,3)
hold off
histogram(err2)
hold on
plot(ss(k+1),0,'k.','MarkerSize',25)

figure(2)
subplot(1,1,1)
hold off
plot(ss(k+1),0,'k.','MarkerSize',25)
hold on
han0 = histogram(err0,10);
han1 = histogram(err1,10);
han2 = histogram(err2,10);
set(han0,'FaceColor','r')
set(han1,'FaceColor','g')
set(han2,'FaceColor','b')
export_fig -png -trans fig_RSVD_accuracy_histogram_fast.png

figure(3)
errmax = max([err0,err1,err2]);
xbins  = linspace(ss(k+1),errmax,50);
subplot(1,3,1)
hold off
plot(ss(k+1),0,'k.','MarkerSize',25)
hold on
histogram(err0,xbins)
subplot(1,3,2)
hold off
histogram(err1,xbins)
hold on
plot(ss(k+1),0,'k.','MarkerSize',25)
subplot(1,3,3)
hold off
histogram(err2,xbins)
hold on
plot(ss(k+1),0,'k.','MarkerSize',25)


keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_fastdecay(m,n)

acc   = 1e-10;  
p     = min([m,n]);
[U,~] = qr(randn(m,p),0);
[V,~] = qr(randn(n,p),0);
ss    = acc.^linspace(0,1,p);
A     = U*diag(ss)*V';

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_slowdecay(m,n)

p     = min([m,n]);
[U,~] = qr(randn(m,p),0);
[V,~] = qr(randn(n,p),0);
ss    = 1./(1:p);
A     = U*diag(ss)*V';

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_S_curve(m,n,ktarg,acc)

% Determine the "inner dimension"
%r = min([m,n,2000]);
r = min([m,n]);

tt    = linspace(-1,1,ktarg);
uu    = [0.5*(1+tanh(5*tt)),ones(1,r-ktarg)];
ss    = 10.^(log10(acc)*uu);
[U,~] = qr(randn(m,r),0);
[V,~] = qr(randn(n,r),0);
A     = (U.*(ones(m,1)*ss))*V';

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A  = LOCAL_smooth(m,n)

[JJ,II]  = meshgrid(1:n,1:m);
A        = 1./(II+JJ);
[~,indi] = sort(rand(1,m));
[~,indj] = sort(rand(1,n));
A        = A(indi,indj);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_helmholtz(m,n,diam)

xxs = rand(2,n);
xxt = [2+rand(1,m);...
         rand(1,m)];

DD1 = xxt(1,:)'*ones(1,n) - ones(m,1)*xxs(1,:);
DD2 = xxt(2,:)'*ones(1,n) - ones(m,1)*xxs(2,:);

A   = besselh(0,2*pi*diam*sqrt(DD1.*DD1 + DD2.*DD2));

return