Assignment 4 Scripts

P1

% P1
% Checks out QRTLS
clc
randn('seed',0)
m = 20; n = 10; A = randn(m,n); b = randn(m,1);

[U,S,V] = svd([A b]);
xTLS = -V(1:n,n+1)/V(n+1,n+1);

disp(sprintf('[m,n] = [%1d,%1d]',m,n))
disp(sprintf('sigma(n+1)/sigma(n+1) = %10.3e\n',S(n+1,n+1)/S(n,n)))
disp('    tol       #Iterations        Relative Error')
disp('-----------------------------------------------------')
for tol = [.01  .0001  .000001   .00000001]
   [x,j] = QRTLS(A,b,tol,100);
   disp(sprintf('%10.1e      %4d         %16.3e',tol,j,norm(x-xTLS)/norm(xTLS)))
end

disp(sprintf('\n\n\n'))
disp('Behavior of QRTLS when ratio = sigma(n+1)/sigma(n) approaches 1.')
tol = .000001;
disp(sprintf('tol = %10.3e',tol))
disp(' ')
disp('    ratio       #Iterations        Relative Error')
disp('-----------------------------------------------------')
for ratio = [.9  .99   .999   .9999]
    C = randn(m,n+1);
    [U,S,V] = svd(C);
    S(n+1,n+1) = ratio*S(n,n);
    C = U*S*V';
    A = C(:,1:n);
    b = C(:,n+1);
    xTLS = -V(1:n,n+1)/V(n+1,n+1);
    [x,j] = QRTLS(A,b,tol,100000);
    disp(sprintf('%10.6f      %6d         %16.3e',ratio,j,norm(x-xTLS)/norm(xTLS)))
end
    
  function v = MyCond(R)
% R is n-by-n upper triangular
% v is a unit 2-norm vector that approximates V(:,n) where U'*R*V = S
% is the SVD of R. 

[n,n] = size(R);
v = zeros(n,1);
% Solve Rv = d where d is vector of plus or minus 1's
v(n) = 1/R(n,n);
for k=n-1:-1:1
    s = R(k,k+1:n)*v(k+1:n);
    if s>0
        % Set d(k) = -1 to maximize growth
        v(k) = (-1 - s)/R(k,k);
    else
        % Set d(k) = +1 to maximize growth
        v(k) = ( 1 - s)/R(k,k);
    end
end
v = v/norm(v); 


   

P2

% P2 Tests LowRankLS
%
% Your implementation of LowrankLS should make efficient use of the QR
% factorization. Type "help qr" to find out just what it can do for you.
% 
randn('seed',0)
m = 1000; n = 500; r = 10;
F = randn(m,r); G = randn(n,r); b = randn(m,1);
x = LowrankLS(F,G,b);
clc
disp('First r components of the solution =')
format long
x1r = x(1:r)

P3

%P3 Checks CS

clc
randn('seed',0)

disp(' n1  n2   p   r     errU1      errU2      errV       errC       errS    diagC   diagS')
disp('--------------------------------------------------------------------------------------------')
n=100;p=10;
for n1=[10 80 95 99]
    n2 = n-n1;
    % generate soem random orthogonal matrices
    [U1,R] = qr(randn(n1,n1));
    [U2,R] = qr(randn(n2,n2));
    [V,R]  = qr(randn(p,p));
    C = zeros(n1,p);
    S = zeros(n2,p);
    for r=0:min(p,n2)
        % Generate a problem where Q1 has r small singular values
        c = [.6*rand(r,1);.8+.2*rand(p-r,1)];
        if n2<p
            c(n2+1:p) = ones(p-n2,1);
        end
        C(1:p,1:p) = diag(c);
        Q1 = U1*C*V';
        s = sqrt(1-c.^2);
        if n2<p
            S(1:n2,1:n2) = diag(s(1:n2));
        else
            S(1:p,1:p) = diag(s(1:p));
        end
        Q2 = U2*S*V';
        Q = [Q1;Q2];
        
        % Compute the CS decomposition
        [U1,U2,V,C,S] = CS(Q,n1);
        
        % Make sure U1, U2, V are orthogonal and that C and S are diagonal.
        % And make sure U1'*Q1*V = C and U2'*Q2*V = S
        errU1 = norm(U1'*U1 - eye(n1,n1));
        errU2 = norm(U2'*U2 - eye(n2,n2));
        errV  = norm(V'*V - eye(p,p));
        errC  = norm(U1'*Q1*V - C);
        errS  = norm(U2'*Q2*V - S);
        c = diag(C);
        if sum(~((c==sort(c)) & (norm(c)==norm(C,'fro')) & (c>=0)))==0
            C_Check = '    OK   ';
        else
            C_Check = '    Not OK   ';
        end
        s = diag(S);
        if sum(~((s==-sort(-s)) & (norm(s)==norm(S,'fro')) &(s>=0)))==0
            S_Check = '   OK     ';
        else
            S_Check = '   Not OK     ';
        end
        errors = [errU1 errU2 errV errC errS];
        first = sprintf('%3d %3d %3d %3d',n1,n2,p,r);
        second = sprintf('   %8.1e',errors);
        disp([first second  C_Check  S_Check]) 
    end
end