MTH3320 Computational Linear Algebra
Computational Linear Algebra
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
MTH3320
Computational Linear Algebra
Assignment
This assignment contains four questions for a total of 100 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 Shifts. (20 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. (20 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 (5 marks each):
School of Mathematics 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. (10 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 AlgorithmWithout Shifts.
(30 marks)
You are asked to implement the QR algorithm without shifts 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
2024 2
School of Mathematics 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 shifts in Matlab according to the pseudocode
discussed in class. Your implementation QRWithoutShifts.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] = QRWithoutShifts(H, k)
% Usage: [Q,T] = QRWithoutShifts(H, k) produces an orthogonal matrix
% Q and an upper triangular matrix T so that H = Q*T*Q’ and
% Q’*Q = EYE(SIZE(Q)). The QR algorithm without shifts is used.
% The input H 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(V)).
% Columns of V are eigenvectors and diagonal entries of D are
% eigenvalues.
% The input A is a square matrix
% 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 QRWithoutShifts.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.
2024 3
School of Mathematics Monash University
Now the next task is to use the code you implemented to find the eigenvalues of the
symmetric matrix A provided in test QREig.m.
(d) Download flip vec.m and test QREig.m to test the correctness of (a), (b), and
(c). For the test of (c), report on the convergence versus the number of steps. In-
clude 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.
5. X-Ray Imaging. (20 marks)
Implement the truncated SVD method for reconstructing X-ray images. Download the
data file xray data.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, and the scalar m defines the image size (m×m).
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 = 500, 510, 520, . . . , 1850 to compute the goodness of fit (of the reconstructed
image) to the data, ∥Fx⃗+k − d⃗n∥, 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
measure ∥x⃗+k ∥). Plot the corner you find on the same plot.
(iv) Save your reconstructed image to a vector xr. 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.
Submit your figures in your assignment package and your code in a MATLAB script file
xray driver.m.