Solved–Computational Linear Algebra Assignment 3 –Solution

Description

This assignment contains four questions for a total of 50 marks (equal to 6% of the nal mark for this unit).

Late submissions will be subject to late penalties (see the unit guide for full details).

Please submit the hardcopy in the assignment box of your tutorial class. The code for the programming exercises as well as relevant output are to be submitted as part of your hardcopy submission. In addition, collect all your Matlab m- les in one zip or rar- le and submit this le through Moodle. The Moodle submission is due at the same time as the hardcopy submission.

1. Gradient Computation. (10 marks)

1. 1. Let ~x; ~y 2 Rⁿ. Show that

r_~x(~y^T ~x) = ~y;

where r_~x is the gradient with respect to ~x.

With

f(~x) = ~y^T ~x;

show this by direct computation of the vector components of r_~xf(~x), given by

@f(~x)

@x_i

2. Using the result from part (a), and the result from Tutorial 5, question 3, that, for A 2 R^{n n}, r_~x(~x^T A~x) = A~x + A^T ~x, show that

	~	2	= 2A	T	A~x	2A	T~
	r~xkb A~xk2						b:
~	2		~			T ~
Hint: Write f(~x) = kb A~xk2		as f(~x) = (b A~x)				(b A~x) and expand. Note that

the result of part (b) was also proven in Theorem 3.13, by direct computation of the vector components of r_~xf(~x). Here, you are asked to prove this in an alternative way, using the result from part (a), and the result from Tutorial 5, question 3.

School of Mathematical Sciences Monash University

2. SOR Preconditioner. (10 marks)

~

The iteration formula for the Successive Over-Relaxation (SOR) method for A~x = b is given by:

xinew = (1 !) xiold + ! aii		bi=1	aij xjnew j=i+1	aij xjold!;
1		i 1	n
		Xj	X

with ! a xed weight 2 (0; 2).

Using the decomposition A = A_D A_L A_U , write the formula in the standard form of the update formula for a one-step stationary iterative method,

~x_k+1 = ~x_k + P~r_k;

to derive the preconditioning matrix P for SOR.

3. Implementation of the ALS Algorithm for Movie Rec-ommendation. (15 marks)

Read the pdf notes on movie recommendation (sections 1.3 and 3.5), and make sure you understand the solutions of Questions 3, 4 and 5 from Tutorial 6.

(a) Download the les

system_movie_j.m

driver_ALS_test.m

from Moodle.

Write a Matlab function with header

function [Ai bi]=system_user_i(Q,M,lambda,i)

~

that computes the matrix A_i and RHS vector b_i for updating user vector ~u_i in the second phase of an ALS iteration (with M xed) according to Eq. (3.17) in the pdf notes, based on the transpose of the ratings matrix, Q = R^T , current approximation for the user matrix M, and regularisation parameter .

Note: This is the dual to function system movie j.m, and is, in fact, very similar to it.

Write a Matlab function with header function [U M]=myALS(R,U,M,lambda)

which performs one ALS iteration, by rst calling system movie j for each column m~_j of M and solving the system, with U xed, and by then calling system user i for each column ~u_i of U and solving the system.

School of Mathematical Sciences Monash University

Write a Matlab function with header function g=gValue(R,U,M,lambda)

to compute the value of the function we optimise, g(U; M) (as in Eq. (3.11) in the notes).

Test these implementations by executing the script driver ALS test.m, which

carries out 500 ALS iterations for the small toy problem of Question 3 of Tutorial 6, and displays the error and value of g(U; M) as a function of iterations. Hand in printouts of the error and iteration plots, along with printouts of your Matlab les.

2. In this question, we will apply the ALS recommendation algorithm to some real-life movie rating data. We will use the data gathered by the non-commercial MovieLens project, see https://grouplens.org/datasets/movielens/.

Download the data set

ml-latest-small.zip

from https://grouplens.org/datasets/movielens/. This data set contains 100,004 ratings from 671 users who have rated 9,066 movies between 1995 and 2016 with rating values between 0.5 and 5.0 (in increments of 0.5).

Read the README le on the website, which contains information about the data les.

Have a look a the data les (you can use xcel or any editor to read the csv for-mat). The ratings are contained in ratings.csv. The movie titles and tags are in movies.csv.

Download the les

movielens_load_data.m

getTitle.m

from Moodle. This is the same code you were given for Question 5 of Tuto-rial 6. If you have not done this as part of Tutorial 6, please take a close look at the code in movielens load data.m and analyse what each part of the code does. Run the code and investigate the output and the gures it generates. Make sure you can relate the output to the code parts that generate it. Note that movielens load data.m writes several les to disk that contains parts of the Movielens data processed into the format of Matlab arrays.

Write a Matlab function driver_ALS.m

that runs the ALS algorithm on the Movielens dataset. Model driver ALS.m after driver ALS test.m, and after parts of movielens load data.m (see below). At the beginning of the le, you should load the Movielens data that you will need:

School of Mathematical Sciences Monash University

function [x res steps]=myGMRES(A,x0,b,tol,maxit)

1. • Performs a sequence of GMRES iterations on

1. • A x = b using the initial estimate x0, until the 2-norm

1. • of the residual is smaller than tol (relative

1. • to the initial residual), or until the number of iterations

1. • reaches maxit. `steps’ is the number of steps

1. • that were performed, and `res’ is a vector with the

1. • residual size after every interation (the 2-norm of

1. • the residual vector).

2. Implement the following in m- le myGMRES SSOR.m:

function [x res steps]=myGMRES_SSOR(A,x0,b,tol,maxit)

1. • Performs a sequence of right-preconditioned GMRES iterations

1. • on A x = b with the initial estimate x0, using the SSOR preconditioner,

1. • until the 2-norm of the residual is smaller than tol (relative

1. • to the initial residual), or until the number of iterations

1. • reaches maxit. `steps’ is the number of steps

1. • that were performed, and `res’ is a vector with the

1. • residual size after every interation (the 2-norm of

1. • the residual vector).

3. Implement the following in m- le myCG SSOR.m:

function [x res steps]=myCG_SSOR(A,x0,b,tol,maxit)

1. • Performs a sequence of preconditioned CG iterations

1. • on A x = b with the initial estimate x0, using the SSOR preconditioner,

1. • until the 2-norm of the residual is smaller than tol (relative

1. • to the initial residual), or until the number of iterations

1. • reaches maxit. `steps’ is the number of steps

1. • that were performed, and `res’ is a vector with the

1. • residual size after every interation (the 2-norm of

1. • the residual vector).

Implementation notes:

k ^~k

{ for GMRES, you can solve the small least-squares problem min M~y f

n^~

by using Matlab’s backslash operator as in ~y = M f (which solves over-determined systems in the least-squares sense using QR decomposition)

{ to build the preconditioner, you can use tril and triu to extract the strictly lower and upper triangular parts of A; you can use spdiags to extract the diagonal of A, and again to build a sparse diagonal matrix containing the diagonal of A (see Question 3 of Tutorial 7)

{ take ! = 1:9 in SSOR

School of Mathematical Sciences Monash University

{ for doing the two triangular solves when applying the SSOR preconditioner, you can again use backslash in Matlab (which recognizes triangular matrices, and solves the system uses forward or backward substitution, as appropriate) (see Question 3 of Tutorial 7)

Submit your m- les to the Moodle drop box, and provide a printout in your as-signment package to be submitted in the assignment boxes.

3. Download test iterative.m to test the correctness of (a{c). Report on the errors and number of steps. (The errors should be close to 1e-14.)

3. Make a driver program driverPCG PGMRES.m that applies myGMRES.m, myGMRES SSOR.m,

~

and myCG SSOR.m to A~x = b with A being the 2D model problem matrix gener-

~

ated by build laplace 2D kron.m, and b a vector with all ones. Use maxit=500 and tol=1e-8. Generate a plot in which you compare, for N = 30, the residual convergence for the three methods (starting from a zero initial guess), as a func-tion of the iteration. On this plot, for each of the methods, plot the 10-log of the residuals as a function of iteration number. Also compare with CG, using your le conjGrad.m from assignment 2 question 4 (or, a simpli cation of myCG SSOR.m, without preconditioner).

(e) Extending driverPCG PGMRES.m, for problem sizes N equal to [80 90 100 110 120 130 140], make a loglog plot of the number of iterations required for con-vergence as a function of total problem size n = N², for each of the 4 methods (CG, CG+SSOR, GMRES, GMRES+SSOR), all on the same gure. Use the polyfit command to t a polynomial of degree 1 to the measured (n, iterations) points of

problem size versus iterations, for each of the methods. The tted slope for CG p

and GMRES should be close to 1/2, because the number of iterations is O( n) for these methods. It can be shown that SSOR is an optimal preconditioner for CG for the 2D model problem in terms of asymptotic order, with number of iterations O(n¹⁼⁴).

Note: This means that the total work is O(n⁵⁼⁴) when using the SSOR precon-ditioner, and O(n³⁼²) without it (both for CG and GMRES). GMRES gives very similar results for the 2D model problem than CG, in terms of the number of re-quired iterations. However, GMRES is much more versatile since it can also be applied to non-SPD systems.

6