Homework #2 of Advanced Programing Methodology at SNU

AMC_homework2.m

%STABILITY OF QR
N = 100;
C = [];
S0 = []; S1 = []; S2 = []; T0 = []; T1 = []; T2 = [];
for k = 1:N
    c = 10^(10*rand); C = [C c];
    % Create m by n (m >= n) matrix with condition number c
    m = 1 + round(100*rand);
    n = min(1+round(m*rand),m);
    A = givcond(m,n,c);
    [Q,R] = clgs(A); % Classical Gram-Schmidt
    S0 = [S0 norm(Q*R-A)];
    T0 = [T0 norm(Q'*Q-eye(n))];
    [Q,R] = mgs(A); % Modified Gram-Schmidt
    S1 = [S1 norm(Q*R-A)]
    T1 = [T1 norm(Q'*Q-eye(n))];
    disp(['k=' int2str(k)]);
end
figure(1);
loglog(C,S0,'*r*', 'MarkerSize',10);
hold on
loglog(C,S1,'*r*', 'MarkerSize',10);
xlabel('condition number', 'Fontweight', 'bold','Fontsize', 18);
ylabel('||QR-A||_2','Fontweight','bold','Fontsize', 18);
title('stability in QR decomposition','Fontweight','bold' ,'Fontsize', 18);
legend('Gram-Schmidt','Modified Gram-Schmidt');
hold off

figure(2);
loglog(C,T0,'*r','MarkerSize',10);
hold on
loglog(C,T1,'+y','MarkerSize',10);
xlabel('condition number', 'Fontweight', 'bold','Fontsize', 18);
ylabel('||Q^*Q-I||_2','Fontweight','bold','Fontsize', 18);
title('stability in QR decomposition','Fontweight','bold' ,'Fontsize', 18);
legend('Gram-Schmidt','Modified Gram-Schmidt');

hold off

function A = givcond(m,n,c)
c = max(abs(c),1);
p = min(m,n);
[U,X] = qr(randn(m,p),0);
[V,X] = qr(randn(n,p),0);
S = diag(logspace(0,-log10(c),p));
A = U*S*V';


% The Classical style of Gram Schmidt orthogonalization
function [Q,R] = clgs(A)
% Keep here the length of the matrix A, for iteration counters
len = size(A,2);
for j = 1:len
  if j>1
    % For elements after the first one, we should run these calculations
    R(1:j-1,j) = A(:,1:j-1)'*A(:,j);
    A(:,j) = A(:,j)-A(:,1:j-1)*R(1:j-1,j);
  end
  % We should keep normalized vectors in Q.
  R(j,j) = norm(A(:,j),2);
  A(:,j) = A(:,j)/R(j,j);
end
Q = A;

% The Modified style of Gram Schmidt orthogonalization
function [Q,R] = mgs(A)
% Keep here the length of the matrix A, for iteration counters
len = size(A,2);
R = zeros(len);
for j = 1:len
  R(j,j) = norm(A(:,j),2);
  A(:,j) = A(:,j)/R(j,j);
  if j==len
    break
  end
  % If we have not reached the final element, we should
  % run these calculations.
  R(j,j+1:len) = A(:,j)'*A(:,j+1:len);
  A(:,j+1:len) = A(:,j+1:len)-A(:,j)*R(j,j+1:len);
end
Q = A;