channelflow.org

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== IAM 961 HW 3 ======



**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)

<code> 
[Q,R] = qr_cgs(A);
[Q,R] = qr_mgs(A);
[Q,R] = qr_house(A);
</code>


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

<code> 
x = backsub(R,b)
</code>


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

<code> 
[Q,R] = qrfix(Q,R); 
</code>

(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

<code>
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
</code>

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

<code>
x = randn(m,1);
b = A*x;
</code>


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

<code>
[Qc,Rc] = qr_cgs(A); 
[Qc,Rc] = qrfix(Qc,Rc);
</code>


Now see how well the computed QR decomp does what it should!

<code>
% 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)
</code>


In particular, look at ''Qc* Qc'' to see exactly where it deviates from 
the identity

<code>
% 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')
</code>


**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.** 

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

Hint: 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''. Also, rerun each of the 


**3:** For the ambitious: Note that if you rerun the above commands repeatedly for 
different randomly constructed ''A,x,b'', you might get errors ranging over a couple 
different orders of magnitude. It's actually best to 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. Use a geometric mean rather than the arithmetic mean, i.e.  

<code>
q_error = 0;

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

q_error = exp(q_error/Ntrials)
</code>

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  
particular-case answers of **3**. 

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

User Tools

Site Tools

Page Tools
