% Matlab code QRStability.m
% For "Applied Numerical Linear Algebra",  Question 3.2
% Written by James Demmel, Aug 2, 1996;
%                Modified, Jun 2, 1997
%
% Compare the numerical stability of 3 algorithms for computing
% the QR decomposition: 
%   QR using Householder transformations
%   CGS (Classical Gram-Schmidt)
%   MGS (Modified Gram-Schmidt)
%
% Inputs:
%
%   m, n = numbers of rows, columns in test matrices
%          m should be at least n
%   samples = number of test matrices to generate
%             samples should be at least 1
%   cnd = condition number of test matrices to generate
%         (ratio of largest to smallest singular value)
%         cnd should be at least 1
%
% Outputs:
%
%   For each algorithm, the maximum value of the residual
%   norm(A-Q*R)/norm(A), which should be around 
%   macheps = 1e-16 for a stable algorithm
%
%   For each algorithm, a plot (versus sample number)
%   of the orthogonality norm(Q'*Q-I), which should 
%   also be around macheps if the computed Q is orthogonal
%   and the algorithm is stable
%
residual = [];
orth = [];
% Generate samples random matrices
for exnum = 1:samples,
% Generate random matrix A, starting with the SVD of a random matrix
  A=randn(m,n);
  [u,s,v]=svd(A);
% Let singular values range from 1 to cnd, with
% uniformly distributed logarithms
  sd = [1, cnd, exp(rand(1,n-2)*log(cnd))];
  s = diag(sd);
  A=u(:,1:n)*s*v';
% Perform CGS (A = Qcgs*Rcgs) and MGS (A = Qmgs*Rmgs)
  Rcgs=[]; Rmgs=[];
  Qcgs=[]; Qmgs=[];
  for i=1:n,
    Qcgs(:,i)=A(:,i);
    Qmgs(:,i)=A(:,i);
    for j=1:i-1,
      Rcgs(j,i) = (Qcgs(:,j))'*A(:,i);
      Qcgs(:,i) = Qcgs(:,i) - Rcgs(j,i)*Qcgs(:,j);
      Rmgs(j,i) = (Qmgs(:,j))'*Qmgs(:,i);
      Qmgs(:,i) = Qmgs(:,i) - Rmgs(j,i)*Qmgs(:,j);
    end
    Rcgs(i,i) = norm(Qcgs(:,i));
    Rmgs(i,i) = norm(Qmgs(:,i));
    Qcgs(:,i) = Qcgs(:,i) / Rcgs(i,i);
    Qmgs(:,i) = Qmgs(:,i) / Rmgs(i,i);
  end
% Perform Householder QR
  [Qhouse,Rhouse]=qr(A,0);
% Compute residuals for each algorithm (cnd = norm(A))
  residual(exnum,1:3) = ...
     [ norm(A-Qhouse*Rhouse), norm(A-Qcgs*Rcgs), norm(A-Qmgs*Rmgs) ]/cnd;
% Compute orthogonality of Q for each algorithm 
  ee = eye(n);
  orth(exnum,1:3) = ...
     [ norm(Qhouse'*Qhouse-ee), norm(Qcgs'*Qcgs-ee), norm(Qmgs'*Qmgs-ee) ];
end
% Limit orth to 1 (any larger error is as bad) so it shows up on plot
orth = min(orth,ones(size(orth)));
% Print max residuals
disp(['max( norm(A-Q*R)/norm(A) ) for Householder QR = ', ...
       num2str(max(residual(:,1)))])
disp(['max( norm(A-Q*R)/norm(A) ) for CGS            = ', ...
       num2str(max(residual(:,2)))])
disp(['max( norm(A-Q*R)/norm(A) ) for MGS            = ', ...
       num2str(max(residual(:,3)))])
% Plot orthogonalities
clf
% make sure orthogonalities are at least eps to avoid log(0)
orth = max(orth,eps);
semilogy((1:samples),orth(:,1),'b*', ...
         (1:samples),orth(:,2),'r*', ...
         (1:samples),orth(:,3),'g*')
axis([1 samples 1e-16 1])
grid
ylabel('orthogonality')
xlabel(['test matrix number, condition num =',num2str(cnd)])
title('Orthogonality of Q for CGS (red), MGS (green), Householder (blue)')