SCHOOL OF MATHEMATICAL SCIENCES
MTH3320 – 2020 S1
Computational Linear Algebra
Assignment 4
Due date: Friday 12 June, 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. 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 ∈ Rn×n. Suppose that we only want to compute
eigenvalues.
1. 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?
2. Show that the tridiagonal structure is preserved in forming the product A(k) = RQ.
3. Determine the order of total work (in flops) required to get from A(k−1) to A(k).
2. Bidiagonalisation. (10 marks)
Suppose we have a matrix A ∈ Rm×n. Recall the Golub-Kahan bidiagonalisation pro-
cedure and the Lawson-Hanson-Chan (LHC) bidiagonalisation procedure (Section 8.2).
Answer the following questions:
School of Mathematical Sciences Monash University
(i) Workout the operation counts required by the Golub-Kahan bidiagonalisation.
(ii) Workout the operation counts required by the LHC bidiagonalisation.
(iii) Using the ratio m
n
, derive and explain under what circumstances the LHC is com-
putationally more advantageous than the Golub-Kahan.
(iv) Suppose we have a bidiagonal matrix B ∈ Rn×n, show that both B>B and BB>
are tridiagonal matrices.
Hint: recall that the operation counts of the QR factorisation (using Householder reflec-
tion) is about 2mn2− 2
3
n3. You can relate those two bidiagonalisation procedures to the
QR factorisation to work out their operation counts.
3. Matrix reduction. (5 marks)
Consider a matrix A ∈ Rm×n has a full QR factorisation
A = QR, with R =
[
Rˆ
0
]
where Q is an orthogonal matrix and Rˆ is an upper-triangular square matrix. Consid-
ering that the matrix Rˆ has an SVD Rˆ = UΣV >, express the SVD of A in terms of Q,
U , Σ, and V .
4. Implementation of the QR Algorithm Without Shift.
(15 marks)
You are asked to implement the QR algorithm without shift for finding eigenvalues of a
matrix A ∈ Rn×n. This takes three steps:
(a) Implement a function myHess.m that transforms a square matrix to a Hessenberg
matrix, using the Householder reflection. 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
2020 2
School of Mathematical Sciences Monash University
Hint: you need to apply the Householder reflection 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.
(b) 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)
to produce the Hessenberg matrix that can be used in the step (b).
(c) 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.
2020 3
School of Mathematical Sciences Monash University
Now the next task is to use the code you implemented to find the eigenvalues of the
symmetric matrix A provided in matrix data A4.mat.
(d) 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.
(e) Submit your m-files 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.
5. X-Ray Imaging. (10 marks)
Implement the truncated SVD method for reconstructing X-ray images. Download the
data file xray data 64.mat from Moodle. An X-ray imaging problem with the resolution
64×64 is given in the data file. In the data file, the matrix F is the model, the vector d is
the noisy data, the scalar m defines the image size (m×m), and the vector xr represents
a random image used for demonstrating how to plot the reconstructed image.
Complete the following tasks:
(i) Compute the SVD of F , use the truncated SVD with rank-k to reconstruct the
image, i.e., ~x+k = VkΣ
−1
k U
>
k
~d.
(ii) Use k = 100, 101, . . . , 930 to compute the goodness of fit (of the reconstructed
image) to the data, ‖F~x+k − ~dn‖, and the 2-norm of the reconstructed image, ‖~x+k ‖.
(iii) Plot the L-curve and determine the index of the corner of the L-curve (which
indicates the balance point between the representation error and the robustness).
Plot the corner you find on the same plot.
(iv) Plot the reconstructed image using the truncation index found in step (iii). Hint:
the MATLAB command imagesc(reshape(xr, m, m), [-0.01, 0.01]); plots
the reconstructed image defined by the vector xr.
Submit your figures and your code in a MATLAB script file xray driver.m.
2020 4