%
% Solve the linear least squares problem by sampling
%
% Note that the code is not written to optimize the cpu time
%
%-------------------------------------------------------------------
%
% generate A and b 
%rand('state', 2378); 
%randn('state', 23984); 
m = 20000; 
n = 100; 
[U,S,V] = svd(rand(m,n),0); 
S = diag(linspace(1,10^6,n)); 
A = U*S*V';    % A has prescribed condition number 
b = rand(m,1);

% compute the LS solution directly 
x = A\b; 

% compute ``randomized'' LS solution 
sampled_rows = find( rand(m,1) < 10*n*log(n)/m ); 
A1 = A(sampled_rows,:); 
b1 = b(sampled_rows); 
x1 = A1\b1; 

% compare the residual 
xres = norm(A*x - b); 
x1res = norm(A*x1 - b); 
residualratio = x1res/xres     % should be about 1.
