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
