function HW01p02_1

%%% 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.
acc = .05;

%%% Used normally distributed numbers as a fun test.
TestMatrix = test_matrix_fastdecay(m,n);
%TestMatrix = test_matrix_slowdecay(m,n);

%Compute partial SVD to given accuracy
[U,D,V] = SVD_given_tolerance(TestMatrix,acc);
V_star = transpose(V);

%Compute CPQR to given accuracy
[Q1,R1,ind1] = CPQR_given_tolerance(TestMatrix,acc);

%Check to ensure that remainder terms for each process behave.
ErrorCPQR = norm(TestMatrix(:,ind1) - (Q1*R1),'fro');
ErrorSVD = norm(TestMatrix - (U*D*V_star),'fro');
fprintf(1,'Error in CPQR: %8.6f\n',ErrorCPQR);
fprintf(1,'Error in SVD:  %8.6f\n',ErrorSVD);

return

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


function [U,D,V] = SVD_given_tolerance(A,acc)

%Note: Error term associated with the CPQR process will not change after
%computing the partial SVD of A.
    [Q,R,ind] = CPQR_given_tolerance(A,acc);
    
    %Construct the inverse permutation vector, P*.
    k = length(ind);
    P_star(ind) = 1:k;

    %Compute the SVD of RP*
    [U_hat,D,V] = svd(R(:,P_star));
    
    %Right multiply Q by U_hat to find U in our partial SVD of A
    U = Q*U_hat;
    
return

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

function [Q,R,ind] = CPQR_given_tolerance(A,acc)

%Modification of CPQR_given_rank to accept the user's desired accuracy,
%acc.

[m,n] = size(A);
p = min(m,n);
R   = zeros(p,n);
Q   = A;
ind = 1:n;

for j = 1:p
    [~, 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;
    %Compute residual in Frobenius norm.
    residual = norm(Q(:,(j+1):n),'fro');
    %Exit loop if accuracy is achieved.
    if residual <= acc
        k = j;
        break
    end
end

Q = Q(:,1:k);
R = R(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
