function hw01p2

rng(0)

%%% Set basic problem parameters.
m    = 120;   % The matrix A to be compressed is m x n.
n    = 100;   % The matrix A to be compressed is m x n.
p    = min(m,n);

%%% Create a test matrix
A         = test_matrix_fastdecay(m,n);
%A         = test_matrix_slowdecay(m,n);

%%% Compute the column pivoted QR factorization.
[Q,R,ind] = CPQR_given_rank(A,p);

%%% Next we check the quality of the computed factorization
ss    = svd(A);     % Compute singular values of A.
qrerr = zeros(1,p-1); % Vector holding truncation errors in QR factorization.
for j = 1:(p-1)
  qrerr(j) = norm(A(:,ind) - Q(:,1:j)*R(1:j,:));
end

figure(1)
semilogy(0:(p-1),ss,'rx-',...
         1:(p-1),qrerr,'bx-')
xlabel('Rank k')
ylabel('Error ||A - A_k||')
legend('Error in SVD','Error in QR')

return

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

function [Q,R,ind] = CPQR_given_rank(A,k)

n   = size(A,2);
R   = zeros(k,n);
Q   = A;
ind = 1:n;

for j = 1:k
    [~, j_max]     = max(sum(Q(:,j:n).*Q(:,j:n),1));
    j_max          = j_max + j - 1;
    Q(:,[j,j_max]) = Q(:,[j_max,j]);
    R(:,[j,j_max]) = R(:,[j_max,j]);
    ind([j,j_max]) = ind([j_max,j]);
    r_jj   = norm(Q(:,j));
    Q(:,j) = Q(:,j)/r_jj;
    Q(:,j) = Q(:,j) - Q(:,1:(j-1))*(Q(:,1:(j-1))'*Q(:,j));
    Q(:,j) = Q(:,j)/norm(Q(:,j));
    R(j,j) = r_jj;
    rr     = Q(:,j)'*Q(:,(j+1):n);
    R(j,(j+1):n) = rr;
    Q(:,(j+1):n) = Q(:,(j+1):n) - Q(:,j)*rr;
end

Q = Q(:,1:k);

return


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

function A = test_matrix_fastdecay(m,n)

acc   = 1e-20;  
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 = test_matrix_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