1: Implement Classical Gram-Schmidt, Modified Gram-Schmidt, and Householder algorithms for the QR decomposition of a square matrix A. Your functions should return the Q and R matrices like this (in Matlab syntax)

 
[Q,R] = qr_cgs(A);
[Q,R] = qr_mgs(A);
[Q,R] = qr_house(A);

Implement a backsubstitution algorithm for solving R x = b where R is a square upper-triangular matrix, with function signature

 
x = backsub(R,b)

There is an indeterminacy in the QR decomposition that we'll need to resolve in order to make comparisons of different QR algorithms. Namely, you can change the sign of any column of Q and the sign of the same row of R without changing the product QR. Write a qrfix that makes appropriate sign changes so that all diagonal elements of R are nonnegative. The function signature should be like this

 
[Q,R] = qrfix(Q,R); 

(Note: production-quality code should check for things like non-square input matrices, encountering linear dependence among columns of A and hence inability to produce a new direction q by the core algorithm, etc. But feel free to gloss over such issues for this homework.)

2: Investigate the behavior of your three QR algorithms as follows:

(Note: Performing the following operations with a smaller (say 10 x 10) well-behaved (condition number of 100 or so) matrix might help you debug your QR algorithm codes)

Construct a random 100 x 100 matrix A with a known QR decomposition and a condition number of approximately 10^7

m = 100;                 % investigate QR on m x m matrices
[U,tmp] = qr(randn(m));  % get a random unitary U
[V,tmp] = qr(randn(m));  % get a random unitary V

% Construct a set of singular values in range 10^-3 to 10^4
s = sort(rand(1,m).* 10.^linspace(-3,4,m), 'descend');
S = diag(s);
Atmp = U*S*V';  

[tmp, R] = qr(Atmp);      % get R from the QR decomp of Atmp
[Q, tmp] = qr(randn(m));  % get a random unitary Q

[Q,R] = qrfix(Q,R);       % remove sign indeterminacy in Q,R

A = Q*R;                  % assign A = QR, so A has known QR dec

From this A, construct an Ax=b problem with a known solution x

x = randn(m,1);
b = A*x;

Compute the QR decomposition via Classical Gram-Schmidt and fix indeterminacy

[Qc,Rc] = qr_cgs(A); 
[Qc,Rc] = qrfix(Qc,Rc);

Now see how well the computed QR decomp does what it should by computing five error measures!

% How far off is Qc* Qc from identity?
unitary_error = norm(Qc'*Qc-eye(m))

% How far off is Rc from upper triangularity?
uptrian_error = norm(Rc - triu(Rc))

% How far off is Qc from the known Q?
q_error = norm(Q - Qc)

% How far off is Rc from the known R?
r_error = norm(R - Rc)

% How far off is Qc Rc from A?
qr_error = norm(A - Qc*Rc)

And look at some 2d plots of Qc* Qc to see exactly where it deviates from the identity

% View the matrix of inner products of Qc. 
% Should be an identity matrix (but won't be!)
imagesc(abs(Qc'*Qc))
colorbar
title('| Q* Q |')
xlabel('i')
ylabel('j')

% View how far off the matrix of inner products of Qc 
% is from the identity, on a log scale
imagesc(log10(abs(Qc'*Qc-eye(m))))
caxis([-16 0])
colorbar
title('log10 | Q* Q - I |')
xlabel('i')
ylabel('j')

Now, keeping A,x,b constant, run that same series of calculations using your qr_cgs, qr_mgs, and qr_house codes for the QR decomposition.

Use script files to automate the above calculations. Put the generation of A,x,b in one script file (or maybe function) and the others in another script so that you can rerun the whole series of calculations for different QR algorithms just by changing qr_cgs to qr_mgs or qr_house.

Run the script a few times for each algorithm. Note that the computed errors might vary over one or two orders of magnitude based on the particular random A,x,b. Again, running the script a few times for each algorithm, record the nearest usual power of ten for each error, and make a table of the 5 error measures for the 3 algorithms.

Based on your results, how well do each of the three QR algorithms work? Make some general observations based on the plots and the error table. Keep in mind that errors of order 10^-16 are superb, 10^-8 are ok, and 10^0 (one) are very bad, and that the order of magnitude (exponent of 10) is all that matters.

3: For the ambitious: To be more precise about the averaging over random A,x,b, run the tests repeatedly and take an average of the errors. Put the entire sequence of above commands (minus the plots) in a for loop over, say, 100 trials, and reduce the dimension of the matrices to say m=30 so they run faster. Then compute the geometric mean errors as follows

q_error = 0;

for n = 1:Ntrials
  ...
  q_error = q_error + log(norm(Q - Qc));
end

q_error = exp(q_error/Ntrials)

Put the various calculated geometric-mean errors in a table and base your discussion of the behavior of the three algorithms on the geometic means rather than the estimated order-of-magnitude errors of 3.

4: For those seeking world domination: calculate the geometric-mean errors for each algorithm as a function of condition number, and produce log-log plots for each error type with three lines, one each for CGS, MGS, and Householder, with condition number ranging from 1 to 10^16.

Turn in your codes, the error table, one inner-product matrix plot for each algorithm, and your discussion of the behavior of the three algorithms.