function hw02p6

rng('default')
rng(0)

%%% Set problem parameters.
m = 100;
n = 140;
p = 10;

%%% Construct a test matrix.
%A  = LOCAL_fast_decay(m,n,200);
%A  = LOCAL_slow_decay(m,n);
A  = LOCAL_helmholtz(m,n,10);

%%% 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));
err1f = 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');
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.-')
axis([0,kvec(end),ss(kvec(end)+1),ssf(1)])
legend('svd','rsvd')
title('Spectral norm errors')
subplot(1,2,2)
hold off
semilogy(0:(length(ss)-1),ssf,'k-',...
         kvec,err1f,'r.-')
axis([0,kvec(end),ss(kvec(end)+1),ssf(1)])
legend('svd','rsvd')
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate a test matrix with decaying spectrum. We emulate a problem
% that arises when solving the Helmholtz equation. The input parameter "kh"
% controls the rank of the matrix. The higher kh, the higher the numerical
% rank.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_helmholtz(m,n,kh)

xxt = [2+randn(1,m);randn(1,m)];
xxs = [  randn(1,n);randn(1,n)];
DD  = sqrt((xxt(1,:)'*ones(1,n) - ones(m,1)*xxs(1,:)).^2 + ...
           (xxt(2,:)'*ones(1,n) - ones(m,1)*xxs(2,:)).^2);
A   = besselj(0,kh*DD);

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);
[Uhat,D,V] = svd(Q'*A,'econ');
U          = Q*Uhat(:,1:k);
D          = D(1:k,1:k);
V          = V(:,1:k);

return

