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
by computing **five error measures**!

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


And look at some 2d plots of ''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.** 

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  

<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  
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.
channelflow.org

User Tools

Site Tools

Page Tools
