function hw03p02

%DRIVER_symmetric
DRIVER_nonsymmetric

return

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

function DRIVER_symmetric

rng('default')
rng(0)

%%% Set problem parameters.
n = 100;
p = 0; % Over sampling parameter

%%% Construct a test matrix.
A = LOCAL_fast_decay_sym(n,200);
%A = LOCAL_slow_decay_sym(n);

%%% Run accuracy tests for different values of k.
%%% Observe that we use the SAME random matrix in every experiment.
kvec  = 5:5:80;
err1  = zeros(1,length(kvec));
err2  = zeros(1,length(kvec));
err3  = zeros(1,length(kvec));
err4  = zeros(1,length(kvec));
err1f = zeros(1,length(kvec));
err2f = zeros(1,length(kvec));
err3f = zeros(1,length(kvec));
err4f = zeros(1,length(kvec));
G     = randn(n,2*n);
for j = 1:length(kvec)
  k = kvec(j);
  [U,D,V]  = LOCAL_rsvd(A,k,p,G);
  err1(j)  = norm(A - U*D*V');
  err1f(j) = norm(A - U*D*V','fro');
  [U,D,V]  = LOCAL_rsvd_power(A,k,p,1,G);
  err2(j)  = norm(A - U*D*V');
  err2f(j) = norm(A - U*D*V','fro');
  [U,D,V]  = LOCAL_rsvd_power(A,k,p,2,G);
  err3(j)  = norm(A - U*D*V');
  err3f(j) = norm(A - U*D*V','fro');
  [U,D]    = LOCAL_lanczos(A,randn(n,1),k);
  err4(j)  = norm(A - U*D*U');
  err4f(j) = norm(A - U*D*U','fro');
end

%%% Compute the errors from truncating the SVD.
ss  = svd(A);
ssf = sqrt(triu(ones(length(ss)))*(ss.*ss));

%%% Create error plots.
figure(1)
subplot(1,2,1)
hold off
semilogy(0:(length(ss)-1),ss,'k-',...
         kvec,err1,'r.-',...
         kvec,err2,'b.-',...
         kvec,err3,'m.-',...
         kvec,err4,'g.-','LineWidth',2)
axis([0,kvec(end),ss(kvec(end)+1),ssf(1)])
xlabel('k')
ylabel('||A - A_k||')
legend('svd','rsvd','rsvdpower q=1','rsvdpower q=2','lanczos')
title('Spectral norm errors')
subplot(1,2,2)
hold off
semilogy(0:(length(ss)-1),ssf,'k-',...
         kvec,err1f,'r.-',...
         kvec,err2f,'b.-',...
         kvec,err3f,'m.-',...
         kvec,err4f,'g.-','LineWidth',2)
xlabel('k')
ylabel('||A - A_k||')
axis([0,kvec(end),ss(kvec(end)+1),ssf(1)])
legend('svd','rsvd','rsvdpower q=1','rsvdpower q=2','lanczos')
title('Frobenius norm errors')
keyboard
return

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

function DRIVER_nonsymmetric

rng('default')
rng(0)

%%% Set problem parameters.
n = 400;
p = 0; % Over sampling parameter

%%% Construct a test matrix.
%A = LOCAL_fast_decay(n,n,200);
A = LOCAL_fast_decay_normal(n,200);
%A = LOCAL_slow_decay_normal(n);

%%% Run accuracy tests for different values of k.
%%% Observe that we use the SAME random matrix in every experiment.
kvec  = 5:5:160;
err1  = zeros(1,length(kvec));
err2  = zeros(1,length(kvec));
err3  = zeros(1,length(kvec));
err4  = zeros(1,length(kvec));
err1f = zeros(1,length(kvec));
err2f = zeros(1,length(kvec));
err3f = zeros(1,length(kvec));
err4f = zeros(1,length(kvec));
G     = randn(n,2*n);
for j = 1:length(kvec)
  k = kvec(j);
  [U,D,V]  = LOCAL_rsvd(A,k,p,G);
  err1(j)  = norm(A - U*D*V');
  err1f(j) = norm(A - U*D*V','fro');
  [U,D,V]  = LOCAL_rsvd_power(A,k,p,1,G);
  err2(j)  = norm(A - U*D*V');
  err2f(j) = norm(A - U*D*V','fro');
  [U,D,V]  = LOCAL_rsvd_power(A,k,p,2,G);
  err3(j)  = norm(A - U*D*V');
  err3f(j) = norm(A - U*D*V','fro');
  [U,D,V]  = LOCAL_arnoldi(A,randn(n,1),k);
  err4(j)  = norm(A - U*D*V');
  err4f(j) = norm(A - U*D*V','fro');
end

%%% Compute the errors from truncating the SVD.
ss  = svd(A);
ssf = sqrt(triu(ones(length(ss)))*(ss.*ss));

%%% Create error plots.
figure(1)
subplot(1,2,1)
hold off
semilogy(0:(length(ss)-1),ss,'k-',...
         kvec,err1,'r.-',...
         kvec,err2,'b.-',...
         kvec,err3,'m.-',...
         kvec,err4,'g.-', 'LineWidth',2)
axis([0,kvec(end),ss(kvec(end)+1),ssf(1)])
xlabel('k')
ylabel('||A - A_k||')
legend('svd','rsvd','rsvdpower q=1','rsvdpower q=2','lanczos')
title('Spectral norm errors')
subplot(1,2,2)
hold off
semilogy(0:(length(ss)-1),ssf,'k-',...
         kvec,err1f,'r.-',...
         kvec,err2f,'b.-',...
         kvec,err3f,'m.-',...
         kvec,err4f,'g.-',...
         'LineWidth',2)
xlabel('k')
ylabel('||A - A_k||')
axis([0,kvec(end),ss(kvec(end)+1),ssf(1)])
legend('svd','rsvd','rsvdpower q=1','rsvdpower q=2','lanczos')
title('Frobenius norm errors')
keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function build a test matrix A whose singular values decay
% exponentially. To be precise, it builds a matrix A via
%   A = U * D * V'
% where U and V are randomly drawn ON matrices, and D is diagonal. The
% entries of D are taken to be D(i,i) = beta^(i-1) where beta is chosen so
% that D(k,k) = 1e-15, where k is an input parameter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_fast_decay(m,n,k)

p       = min(m,n);
[U,~,~] = qr(randn(m,p),0);
[V,~,~] = qr(randn(n,p),0);
beta    = (1e-15)^(1/(k-1));
ss      = beta.^(0:(p-1));
A       = U * diag(ss) * V';
    
return

function A = LOCAL_fast_decay_sym(n,k)

[U,~,~] = qr(randn(n),0);
beta    = (1e-15)^(1/(k-1));
ss      = beta.^(0:(n-1));
A       = U * diag(ss) * U';
A       = 0.5*(A + A');
    
return

function A = LOCAL_fast_decay_normal(n,k)

nh            = round(n/2);
beta          = (1e-15)^(1/k);
ss            = beta.^(2*(0:(nh-1)));
th            = 2*pi*rand(1,nh);
dd            = zeros(1,n);
dd(1:2:(n-1)) = ss.*cos(th);
dd(2:2:n)     = ss.*cos(th);
od            = zeros(1,n-1);
od(1:2:end)   = ss.*sin(th);
X             = diag(dd) - diag(od,1) + diag(od,-1);
[U,~,~]       = qr(randn(n),0);
A             = U*X*U';
    
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function build a test matrix A whose singular values decay
% slowly. To be precise, it builds a matrix A via
%    A = U * D * V'
% where U and V are randomly drawn ON matrices, and D is diagonal. The
% entries of D are taken to be D(i,i) = 1/i.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_slow_decay(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_slow_decay_sym(n)

[U,~,~] = qr(randn(n),0);
ss      = 1./(1:n);
A       = U * diag(ss) * U';
A       = 0.5*(A + A');
    
return

function A = LOCAL_slow_decay_normal(n)

nh            = round(n/2);
ss            = 1./(1:nh);
th            = 2*pi*rand(1,nh);
dd            = zeros(1,n);
dd(1:2:(n-1)) = ss.*cos(th);
dd(2:2:n)     = ss.*cos(th);
od            = zeros(1,n-1);
od(1:2:end)   = ss.*sin(th);
X             = diag(dd) - diag(od,1) + diag(od,-1);
[U,~,~]       = qr(randn(n),0);
A             = U*X*U';
    
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the single pass randomized algorithm to compute an
% approximate low rank SVD of a given matrix A. In other words, it computes
% an approximate factorization
%    A \approx U * D * V'
% where U and V and ON matrices with k columns, and D is diagonal.
% 
% INPUTS:  A       Given matrix to be factorized.
%          k       Target rank.
%          p       Over-sampling parameter.
%
% OUTPUTS: U,D,V   Factors in approximate SVD.
%
% NOTE: In this particular code, the random matrices Gc and Gr are taken as
% input parameters. This is artificial, and is not how the code would
% normally work - it is done simply to make testing across several k
% consistent.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [U,D,V] = LOCAL_single_pass(A,k,p,Gc,Gr)

%%% Execute Stage A to build the sampling and basis matrices.
Gc       = Gc(:,1:(k+p));
Gr       = Gr(:,1:(k+p));
Yc       = A *Gc;
Yr       = A'*Gr;
[Qc,~,~] = svd(Yc,'econ');
[Qr,~,~] = svd(Yr,'econ');
Qc       = Qc(:,1:k);
Qr       = Qr(:,1:k);

%%% Build an approximation to the projected matrix C = Qc' * A * Qr
C = (Gr'*Qc)\(Yr'*Qr);
%C = (Qc'*Yc)/(Qr'*Gc);

%%% Form the SVD of C, and then project back to the "big" space.
[Uhat,D,Vhat] = svd(C);
U             = Qc*Uhat;
V             = Qr*Vhat;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the basic randomized algorithm to compute an
% approximate low rank SVD of a given matrix A. In other words, it computes
% an approximate factorization
%    A \approx U * D * V'
% where U and V and ON matrices with k columns, and D is diagonal.
% 
% INPUTS:  A       Given matrix to be factorized.
%          k       Target rank.
%          p       Over-sampling parameter.
%
% OUTPUTS: U,D,V   Factors in approximate SVD.
%
% NOTE: In this particular code, the random matrix G is taken as an
% input parameter. This is artificial, and is not how the code would
% normally work - it is done simply to make testing across several k
% consistent.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [U,D,V] = LOCAL_rsvd(A,k,p,G)

Y          = A*G(:,1:(k+p));
[Q,~,~]    = qr(Y,0);
B          = Q'*A;
[Uhat,D,V] = svd(B,'econ');
U          = Q*Uhat(:,1:k);
D          = D(1:k,1:k);
V          = V(:,1:k);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the basic randomized algorithm to compute an
% approximate low rank SVD of a given matrix A. In other words, it computes
% an approximate factorization
%    A \approx U * D * V'
% where U and V and ON matrices with k columns, and D is diagonal.
% 
% INPUTS:  A       Given matrix to be factorized.
%          k       Target rank.
%          p       Over-sampling parameter.
%          q       Power scheme parameter.
%
% OUTPUTS: U,D,V   Factors in approximate SVD.
%
% NOTE: In this particular code, the random matrix G is taken as an
% input parameter. This is artificial, and is not how the code would
% normally work - it is done simply to make testing across several k
% consistent.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [U,D,V] = LOCAL_rsvd_power(A,k,p,q,G)

Y          = A*G(:,1:(k+p));
[Q,~,~]    = qr(Y,0);
for j = 1:q
  [Q,~,~] = qr(A'*Q,0);
  [Q,~,~] = qr(A*Q,0);
end
B          = Q'*A;
[Uhat,D,V] = svd(B,'econ');
U          = Q*Uhat(:,1:k);
D          = D(1:k,1:k);
V          = V(:,1:k);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the basic randomized algorithm to compute an
% approximate low rank SVD of a given matrix A. In other words, it computes
% an approximate factorization
%    A \approx U * D * V'
% where U and V and ON matrices with k columns, and D is diagonal.
% 
% INPUTS:  A       Given matrix to be factorized.
%          k       Target rank.
%          p       Over-sampling parameter.
%
% OUTPUTS: U,D,V   Factors in approximate SVD.
%
% NOTE: The output matrices U and V will be COMPLEX even if A is real.
%
% NOTE: In this code, we apply the FULL fft, as opposed to the subsampled
% FFT. This ruins the asymptotic flop count. In other words, this code is
% fine for testing ACCURACY, but should not under any circumstances be used
% to test the speed of the SRFT algorithm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [U,D,V] = LOCAL_rsft(A,k,p)

n          = size(A,2);
dd         = exp(1i*2*pi*rand(1,n));
Y          = A*diag(dd);
Y          = fft(Y,[],2);
[~,ind]    = sort(rand(1,n));
J          = ind(1:(k+p));
Y          = Y(:,J);
[Q,~,~]    = qr(Y,0);
B          = Q'*A;
[Uhat,D,V] = svd(B,'econ');
U          = Q*Uhat(:,1:k);
D          = D(1:k,1:k);
V          = V(:,1:k);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the Lanczos iteration to compute an approximate low 
% rank EVD of a given matrix A. In other words, it computes an approximate 
% factorization
%    A \approx U * D * U'
% where U is an ON matrices with k columns, and D is diagonal.
% 
% INPUTS:  A         Given matrix to be factorized.
%          vec_init  Starting vector.
%          n_steps   Number of steps to be taken.
%
% OUTPUTS: U,D   Factors in approximate EVD.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [U,D] = LOCAL_lanczos(A,vec_init,n_steps)

n_tot = size(A,1);

Q     = zeros(n_tot,n_steps);
a_vec = zeros(1,n_steps);
b_vec = zeros(1,n_steps-1);

Q(:,1)    = vec_init/norm(vec_init);
A_times_q = A*Q(:,1);
a_vec(1)  = Q(:,1)' * A_times_q;
vec_res   = A_times_q - a_vec(1) * Q(:,1);

for j = 2:n_steps
  b_vec(j-1) = norm(vec_res);
  Q(:,j)     = vec_res / b_vec(j-1);
  A_times_q  = A * Q(:,j);
  a_vec(j)   = Q(:,j)' * A_times_q;
  vec_res    = A_times_q - Q(:,j) * a_vec(j) - Q(:,j-1) * b_vec(j-1);
end

T = diag(a_vec) + diag(b_vec,1) + diag(b_vec,-1);
[Uhat,D] = eig(T);
U = Q*Uhat;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the Arnoldi iteration to compute an approximate low 
% rank SVD of a given matrix A. In other words, it computes an approximate 
% factorization
%    A \approx U * D * V'
% where U and V are an ON matrices with k columns, and D is diagonal.
% 
% INPUTS:  A         Given matrix to be factorized.
%          vec_init  Starting vector.
%          n_steps   Number of steps to be taken.
%
% OUTPUTS: U,D,V     Factors in approximate SVD.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [U,D,V] = LOCAL_arnoldi(A,vec_init,n_steps)

n_tot = size(A,1);

Q = zeros(n_tot,n_steps);
H = zeros(n_steps);

Q(:,1)    = vec_init/norm(vec_init);
A_times_q = A*Q(:,1);
H(1,1)    = Q(:,1)' * A_times_q;
vec_res   = A_times_q - H(1,1) * Q(:,1);

for j = 2:n_steps
  H(j,j-1)  = norm(vec_res);
  Q(:,j)    = vec_res / H(j,j-1);
  A_times_q = A * Q(:,j);
  H(1:j,j)  = Q(:,1:j)' * A_times_q;
  vec_res   = A_times_q - Q(:,1:j) * H(1:j,j);
end

[Uhat,D,Vhat] = svd(H);
U = Q*Uhat;
V = Q*Vhat;

return
