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