SCHOOL OF MATHEMATICAL SCIENCES
MTH3320
Computational Linear Algebra
Assignment 3
Due date: Thursday 21 May, 2020, 9pm (submit the
electronic copy of your assignment and the matlab
code via Moodle)
This assignment contains four questions for a total of 50 marks (equal to 7.5% of the
final mark for this unit).
Late submissions will be subject to late penalties (see the unit guide for full
details).
The relevant output of the programming exercises are to be submitted as part of your
electronic submission. You can either include the outputs as a part of your pdf file (if
you can edit the pdf document), or include them in a separate word document (you need
to label the figures and outputs consistently with the question number). In addition,
collect all your Matlab m-files in one zip file and submit this file through Moodle.
1. Unitary Similarity Transformation to a Lower Trian-
gular Matrix. (10 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 ∈ 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.
2. Convergence of the Inverse Iteration on Symmetric
Matrices. (15 marks)
Consider the inverse iteration (Algorithm 6.32) applied to a real symmetric matrix A ∈
Rn×n. Suppose λK is the closest eigenvalue to µ and λL is the second closest, that is,
|λK − µ| < |λL − µ| ≤ |λi − µ|, for each i 6= 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 6= 0. Prove that, at iteration k, the vector
~b(k) in the inverse iteration converges as
‖~b(k) − ~qK‖ = O
(∣∣∣∣λK − µλL − µ
∣∣∣∣k
)
.
You may consider taking the following steps:
(i) Suppose all the eigenvectors ~q1, ~q2, . . . , ~qn are orthonormal, write down the eigen-
decomposition of (A− µI)−1.
(ii) Use induction to show that the inverse iteration method produces a sequence of
vectors ~b(k), k = 1, 2, . . . in the form of
~b(k) =
(A− µI)−k~b(0)
‖(A− µI)−k~b(0)‖ ,
(iii) Use the results in (i) and (ii) to complete the rest of the proof.
3. 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 or build laplace 2D kron.m to build the 2D Laplacian matrix).
(a) Implement the following in m-file 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
% 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).
(b) Implement the following in m-file 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
2020 2
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).
Submit your m-files steepest.m and conjGrad.m, and provide a printout in your
main assignment package.
(c) 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.)
(d) Make a driver program CG driver.m that applies steepest.m and conjGrad.m to
A~x = ~b withA 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?
(e) 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, briefly
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. Implementation of the Inverse Iteration and the Rayleigh
Quotient Iteration. (10 marks)
(a) 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.
2020 3
School of Mathematical Sciences Monash University
% 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
(b) 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.
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
(c) Check the correctness of your implementation of RayleighIteration.m as follows.
Download the Matlab files flip vec.m, test Rayleigh.m and matrix data Q4.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
2020 4
School of Mathematical Sciences Monash University
Step 1 in test Rayleigh.m). Comment on the result plot if your implementation
produces a desirable answer.
(d) Run the same test for your implementation of InverseIteration.m using the
initial shift value given in my test Rayleigh.m and matrix data Q4.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 main assignment
package. Comment on the performance of your implementation of the Rayleigh
quotient iteration, compared with the inverse iteration.
(e) You should submit the code InverseIteration.m, RayleighIteration.m and
my test Rayleigh.m to the Moodle. Also, include a printout of your code in the
main assignment package.
2020 5