$29
Computational Linear Algebra Assignment 2 Solution
This assignment contains four questions for a total of 50 marks (equal to 6% of the nal mark for this unit).
Late submissions will be subject to late penalties (see the unit guide for full details).
Please submit the hardcopy in the assignment box of your tutorial class. The code for the programming exercises as well as relevant output are to be submitted as part of your hardcopy submission. In addition, collect all your Matlab m- les in one zip or rar- le and submit this le through Moodle. The Moodle submission is due at the same time as the hardcopy submission.
1. Householder re ection in R2. (10 marks)
Let
A = ~a1
~a2
=
6
20
:
8
10
Construct the rst Householder re ection matrix, Q1, which re ects the rst col-umn of A, ~a1, to a vector
Q1~a1 =
k~a1k
;
0
i.e., choose the sign according to the rule used to ensure numerical stability, de-termine vector ~v1 and its normalised version ~u1, and then the matrix Q1.
Verify that Q1 is an orthogonal matrix.
Compute
2
Q1 6 ;
the image of the vector (2; 6)T under the re ection.
School of Mathematical Sciences
Monash University
(d) Compute
4
Q1
3
;
the image of the vector ( 3; 4)T
under the re ection.
(e)
Make a sketch (by hand) in the
R
2 plane indicating the vectors ~x = (x; y)T
2
R
2
that arise in parts (a)-(d):
the original vector ~a1, its image Q1~a1 under the
re ection, the vectors ~v1 and ~u1, and the vectors (2; 6)T and ( 3; 4)T and their images under the re ection. Also indicate the line about which the vectors are re ected.
(f) Compute Q1A and write down the QR decomposition of A.
Note: Do all computations by hand. (You can verify the results by Matlab if you want.)
Eigenvalues of a Householder re ector matrix. (10 marks)
Determine the m eigenvalues of an m m Householder re ector matrix
Q = I 2u~T ;
where ~u 2 Rm with k~uk2 = 1, and nd m corresponding eigenvectors.
(Hint: Some of the m eigenvalues may occur multiple times. Rather than trying to compute the eigenvalues by the determinant formula, it is more straightforward to think geometrically. Consider that a Householder re ection is a re ection of vectors in Rm about a hyperplane that is orthogonal to ~u. What does the re ection do to a vector in the hyperplane? (I.e., what does the re ection do to any vector orthogonal to ~u?) What is the dimension of the hyperplane? What does the re ection do to ~u? Sketching a picture of the situation in R3 may help.)
Implementation of QR decomposition by Householder re ection. (15 marks)
Implement the QR algorithm using Householder re ections for A 2 Rm n in Matlab according to the pseudocode discussed in class. Your implementation myHouseholder.m should return the full versions of the factors, i.e. Q 2 Rm m and R 2 Rm n. You should use the function header given below and submit the code to the Moodle drop box on the unit webpage. Also, include a printout of your code in the assignment package submitted to the assignment boxes.
function [Q,R] = myHouseholder(A)
Usage: [Q,R] = myHouseholder(A)
Compute the QR factorization of matrix A using
Householder reflections and return the resulting
factors Q and R.
2
School of Mathematical Sciences Monash University
Check the correctness of your implementation myHouseholder.m as follows. Down-load from the solutions of Tutorial 4 the Matlab les test grams 2.m, myGramSchmidt.m and myGramSchmidtMod.m (read Question 6b and c from Tutorial 4 and make sure you understand the results of running test grams 2.m).
Make a new le test grams 3.m that extends test grams 2.m by adding a test
that computes the QR factorisation of the ill-conditioned matrix using myHouseholder.m. Compute and report the norm of QT Q I and A QR. You will see that QT Q I
is almost as accurate as for Matlab's built-in QR decomposition, and much more accurate than for any of the versions of Gram-Schmidt orthogonalisation. Submit your le test grams 3.m to Moodle.
myHouseholder.m builds QT as R is updated, according to
QT = Qn(: : : (Q2(Q1I)));
and then takes the transpose at the end to obtain Q. However, as discussed in class and in the notes, this does unnecessary work. A more e cient implementation to compute Q is to use the expression
Q = Q1(: : : (Qn 1(QnI))):
This computation is carried out in a separate loop after the loop that computes R has nished, and it requires that the Householder vectors ~uk are stored as R is being computed, for use in the second loop that computes Q. You are asked to implement this second, more e cient mechanism to compute Q, in a new le myHouseholder mod.m. (Note: don't build the actual Qk matrices, but apply I 2~uk~uTk to the relevant part of the intermediate result matrix.)
Make a new le test grams 4.m that modi es test grams 3.m by adding a test that computes the Q factors of the ill-conditioned matrix using the two di erent methods: myHouseholder.m and myHouseholder mod.m. Compute and report the 2-norm of the di erence of the two Q factors.
Implementation of the steepest descent and conjugate gradient methods. (15 marks)
You are asked to implement the steepest descent and conjugate gradient methods for
~
A ~x = b in Matlab, and apply the methods to the 2D model problem (download build laplace 2D.m from Question 9 of Tutorial 2 or build laplace 2D kron.m from Assignment 1 to build the 2D Laplacian matrix).
Implement the following in m- le steepest.m:
function [x res steps]=steepest(A,x0,b,tol,maxit)
Performs a sequence of steepest descent iterations on
A x = b using the initial estimate x0, until the 2-norm
of the residual is smaller than tol (relative
3
School of Mathematical Sciences Monash University
to the initial residual), or until the number of iterations
reaches maxit. `steps' is the number of steps
that were performed, and `res' is a vector with the
residual size after every interation (the 2-norm of
the residual vector).
Implement the following in m- le conjGrad.m:
function [x res steps]=conjGrad(A,x0,b,tol,maxit)
Performs a sequence of conjugate gradient iterations
on A x = b using the initial estimate u0, until the 2-norm
of the residual is smaller than tol (relative
to the initial residual), or until the number of iterations
reaches maxit. `steps' is the number of steps
that were performed, and `res' is a vector with the
residual size after every interation (the 2-norm of
the residual vector).
Submit your m- les steepest.m and conjGrad.m to the Moodle drop box, and provide a printout in your assignment package to be submitted in the assignment boxes.
Download test steepest cg.m to test the correctness of (a) and (b). Report on the errors and number of steps. (The errors should be smaller than 10 10.)
Make a driver program CG driver.m that applies steepest.m and conjGrad.m to
~
A~x = b with A being the 2D model problem matrix generated by build laplace 2D.m
~
or build laplace 2D kron.m, and b a vector with all twos. Use maxit=400 and tol=1e-6. Generate a plot in which you compare, for N = 30, the residual con-vergence for the steepest descent and conjugate gradient methods (starting from a zero initial guess), as a function of the iteration. For each of the methods, plot the 10-log of the residuals as a function of iteration number. Which method requires the least iterations to obtain an accurate solution, and is this as you expected?
Provide a table in which you list, for steepest descent and conjugate gradient, how
many iterations you need to reduce the residual by 1e 4, for N = 16; 32; 64: (You can use maxit = 500.) What are the total number of unknowns and the condition number for each problem size? (You can use cond(full(A)); no need to do this for N = 64 though, because it may take too long.) For each method, do you see the expected behaviour in the number of required iterations as a function of the total number of unknowns? (Explain.) Using these numerical results, brie y discuss the computational cost/complexity of each of the methods as a function of the total number of unknowns (discuss the cost per iteration, the number of iterations required, and the total cost, as a function of total problem size). Which method is best for large problem size?
4
