%
% Use randomized preconditioning to solve the normal equation
% of the linear least squares problem 
%
% Reference: 
% Coakley, Rokhlin, Tygert, SIAM J. Sci. Comput. 33:849-868, 2011 
%
% Note: the code is not written to optimize the cpu time 
%
% ----------------------------------------------------------------
%
% generate matrix m by n A (short, fat) with prescribed singular 
% value distribution (=condition number)
%
m = 20;
n = 100; 
U = orth(randn(m)); 
V = orth(randn(n)); 
S = zeros(m,n); 
smax = 10^10; 
sigma = linspace(1,smax,m); 
for i = 1:m,
   S(i,i) = sigma(i); 
end 
A = U*S*V';
Acondnum  = cond(A) 

%
% Construct randomized preconditioner 
%
ell = m+4; 
G = randn(n,ell); 
S = A*G;     % sketching of A 
[Q,R,E] = qr(S');  % QR decomposition with column pivoting
P = E*R(1:m,1:m)'; % preconditioner P 
condP = cond(P) 

% preconditioned matrix 
B = P\A; 
PAcondnum  = cond(B)

% solver the normal equation associated with min||A^T*x - b||
b = rand(n,1); 

xls = A'\b;   % solve LS directly 

xrand = P'\((B*B')\(B*b));  % solve LS by solving preconditioned normal eq.

%  check 
exact_res = norm(A'*xls - b)
rand_res = norm(A'*xrand - b)
residualratio = exact_res/rand_res     % should be around 1
