Computational Linear Algebra Assignment 4 Solution

This assignment contains ve 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.







 
Unitary Similarity Transformation to a Lower Trian-gular Matrix. (5 marks)




Use induction to show that every square matrix can be transformed to a lower triangular matrix via a unitary similarity transformation. That is, for a matrix A 2 Rn n, we can always identify a unitary matrix Q and a lower triangular matrix L such that




A = QLQ :




Hint: use the same procedure in the proof of the existence of Schur factorisation, and start from the lower-rightmost entry.







 
Convergence of the Inverse Iteration on Symmetric Matrices. (10 marks)




Consider the inverse iteration (Algorithm 6.32) applied to a real symmetric matrix A 2 Rn n. Suppose K is the closest eigenvalue to and L is the second closest, that is,




j K j < j L j j i j; for each i ≠ K:

School of Mathematical Sciences Monash University










Let ⃗q1; : : : ; ⃗qn denote eigenvectors corresponding the eigenvalues of A. Suppose further we have an initial vector ⃗x such that ⃗x ⃗qK ≠ 0. Prove that, at iteration k, the vector

⃗(k)

b in the inverse iteration converges as










K




k














∥⃗b(k)
⃗qK ∥ = O
(
L




) :














You may consider taking the following steps:




 
Suppose all the eigenvectors ⃗q1; ⃗q2; : : : ; ⃗qn are orthonormal, write down the eigen-decomposition of (A I) 1.




 
Use induction to show that the inverse iteration method produces a sequence of

⃗(k)

vectors b , k = 1; 2; : : : in the form of

⃗(k)


(A
I) k⃗b(0)
b
=




;
∥(A






I) k⃗b(0)∥
 
Use the results in (i) and (ii) to complete the rest of the proof.







 
Properties of Eigenvalues of a Special Symmetric Ma-trix. (optional)




Consider a real symmetric tridiagonal matrix A 2 Rn n in the form of



2
a1
b1










3




b1
a2 b2
. . .










A =
6


b2 a3








7
:
6










7


6














7




6


. . . . . .








7




6












7




4








1
bn 1
5




6




an


7




6




bn




1
an
7


Suppose all the entries on the subdiagonal and superdiagonal are nonzero, that is, bi ≠ 0 for i = 1; : : : ; n 1. Prove that all eigenvalues of A are distinct. You may take the following steps.










1. First show that the geometric multiplicity of each eigenvalue of A is 1. Hint: show the rank of the matrix A I is at least n 1 for any 2 R.

 
Show that the symmetric A is not defective.




 
Use the relation between geometric multiplicity and algebraic multiplicity to com-plete the proof.






















2
School of Mathematical Sciences Monash University










4. QR Algorithm Without Shift. (10 marks)




Consider applying one step of the QR Algorithm without shifts (Algorithm 7.11) to a real symmetric tridiagonal matrix A 2 Rn n. Suppose that we only want to compute eigenvalues.







 
In the QR factorisation QR = A(k 1), explain the following: which entries of R are generally nonzero and which entries of Q are generally nonzero?




 
Show that the tridiagonal structure is preserved in forming the product A(k) = RQ.




3. Determine the order of total work (in ops) required to get from A(k 1) to A(k).







 
Implementation of the Inverse Iteration and the Rayleigh Quotient Iteration. (10 marks)




 
Implement the inverse iterations in Matlab according to the pseudocode discussed in class. Your implementation InverseIteration.m should store estimated eigen-




values and eigenvectors in each step, except the step 0. You should use the function header given below.




function [B, lamb] = InverseIteration(A, mu, x, k)




 
 

 
Usage: [B, lamb] = InverseIteration(A, mu, x, k) carries k steps




 
of the inverse iteration for estimating the eigenvalue of A,




 
with an initial shift mu, and an initial vector x.




 
It generates outputs B and lamb.




 
-- B is an (n x k) matrix that stores the estimated eigenvector at




 
i-th step as its i-th column, B(:,i)




 
-- lamb is a (1 x k) vector that stores the estimated eigenvalue at




 
i-th step as its i-th element, lamb(i)




 
Note that the initial vector and initial eigenvalue estimate are




 
not stored.




n = size(A, 1); % size of the matrix




B = zeros(n, k); % matrix B




lamb = zeros(1, k); % vector lamb




% your code




end




 
Implement the Rayleigh quotient iteration in Matlab according to the pseudocode discussed in class. Your implementation RayleighIteration.m should store esti-mated eigenvalues and eigenvectors in each step, except the step 0. You should use the function header given below.













3
School of Mathematical Sciences Monash University










function [B, lamb] = RayleighIteration(A, x, k)




%




 
Usage: [B, lamb] = RayleighIteration(A, x, k) carries k steps




 
of the Rayleigh quotient iteration for estimating the eigenvalue




 
of A, with an initial vector x




 
It generates outputs B and lamb.




 
-- B is an (n x k) matrix that stores the estimated eigenvector at




 
i-th step as its i-th column, B(:,i)




 
-- lamb is a (1 x k) vector that stores the estimated eigenvalue at




 
i-th step as its i-th element, lamb(i)




 
Note that the initial vector and initial eigenvalue estimate are




 
not stored.




n = size(A, 1); % size of the matrix




B = zeros(n, k); % matrix B




lamb = zeros(1, k); % vector lamb




% your code end




 
Check the correctness of your implementation of RayleighIteration.m as follows. Download the Matlab les flip vec.m, test Rayleigh.m and matrix data A4.mat under Assignment 4 on Moodle. Save the test Rayleigh.m as my test Rayleigh.m. Then modify my test Rayleigh.m as the driver to call your RayleighIteration.m to solve an eigenvalue problem. The matrix and initial vector are provided (see Step 1 in test Rayleigh.m). Comment on the result plot if your implementation produces a desirable answer.




 
Run the same test for your implementation of InverseIteration.m using the initial shift value given in my test Rayleigh.m and matrix data A4.mat. Mod-ify the convergence plots in my test Rayleigh.m to also plot the convergence of eigenvalue and eigenvector estimates produced by the inverse iteration on the same graph showing the convergence of the Rayleigh quotient iteration. Include a printout of the plots produced by my test Rayleigh.m in the assignment package submitted to the assignment boxes. Comment on the performance of your imple-mentation of the Rayleigh quotient iteration, compared with the inverse iteration.




 
You should submit the code InverseIteration.m, RayleighIteration.m and my test Rayleigh.m 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.











































4
School of Mathematical Sciences Monash University










 
Implementation of the QR Algorithm Without Shift. (15 marks)










You are asked to implement the QR algorithm without shift for nding eigenvalues of a matrix A 2 Rn n. This takes three steps:




 
Implement a function myHess.m that transforms a square matrix to a Hessenberg matrix, using the Householder re ection. Your function should take the following form:




function [P,H] = myHess(A)




 
Usage: [P,H] = myHess(A) produces a unitary matrix P and a




 
Hessenberg matrix H so that A = P*H*P' and P'*P = EYE(SIZE(P)).




 
Householder reflections will be used.




 
your code




end










You can modify the Householder re ection code in Assignment 2 for this task. Hint: you need to apply the Householder re ection twice in each iteration for computing H. For computing the unitary matrix P , the exactly same technique used in Algorithm 3.10 can be applied.




 
Implement the QR algorithm without shift in Matlab according to the pseudocode discussed in class. Your implementation QRWithoutShift.m should produce an orthogonal matrix Q and a (nearly) upper triangular matrix T after k steps. Your function should take the following form:




function [Q,T] = QRWithoutShift(H, k)




 
Usage: [Q,T] = QRWithoutShift(H, k) produces an orthogonal matrix




 
Q and an upper triangular matrix T so that A = Q*T*Q' and




 
Q'*Q = EYE(SIZE(Q)). The QR algorithm without shift is used.




 
The input A is a Hessenberg matrix that is produced in Step 1.




 
The input k is the number of steps.




 
your code




end




If you do not manage to get the step (a) working. You can use the MATLAB function







[P,H] = hess(A)
















5
School of Mathematical Sciences Monash University













to produce the Hessenberg matrix that can be used in the step (b).




 
Implement the third function QREig.m that calls the above two functions to com-pute the eigenvectors and eigenvalues of a symmetric matrix. Your function should take the following form:




function [V,D] = QREig(A, k)




 
Usage: [V,D] = QREig(A, k) produces an orthogonal matrix V and




 
a diagonal matrix D so that A = V*D*V' and V'*V = EYE(SIZE(Q)).




 
Columns of V are eigenvectors and diagonal entries of D are




 
eigenvalues.




 
The input A is a Hessenberg matrix that is produced in Step 1.




 
The input k is the number of steps.




 
your code




end




Hint: If the input matrix is symmetric, the function myHess.m should produce a matrix H that is tridiagonal (some of the entries below the subdiagonal and above the superdiagonal can have small values close to zero). In the next step, the function QRWithoutShift.m should produce a matrix T that is diagonal (some of the off diagonal entries can have small values close to zero) if k is sufficiently large. Then eigenvalues and eigenvectors can be obtained.







Now the next task is to use the code you implemented to nd the eigenvalues of the symmetric matrix A provided in matrix data A4.mat.







 
Download flip vec.m, test QREig.m and matrix data A4.mat to test the cor-rectness of (a), (b), and (c). For the test of (c), report on the convergence versus the number of steps. Include a printout of the plots produced by test QREig.m in the assignment package submitted to the assignment boxes.




 
Submit your m- les QREig.m, QRWithoutShift.m, and myHess.m to the Moodle drop box, and provide a printout in your assignment package to be submitted in the assignment boxes.





























































6