Computational Linear Algebra Assignment 3 Solution

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.







 
Gradient Computation. (10 marks)




 
Let ~x; ~y 2 Rn. Show that

r~x(~yT ~x) = ~y;

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

With

f(~x) = ~yT ~x;




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




@f(~x)

@xi







 
Using the result from part (a), and the result from Tutorial 5, question 3, that, for A 2 Rn n, r~x(~xT A~x) = A~x + AT ~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 = AD AL AU , write the formula in the standard form of the update formula for a one-step stationary iterative method,




~xk+1 = ~xk + P~rk;




to derive the preconditioning matrix P for SOR.







 
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 Ai and RHS vector bi for updating user vector ~ui 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 = RT , 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 ~ui of U and solving the system.
















2
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.




 
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:
















3
School of Mathematical Sciences
Monash University




load Rsp


load
titles


load
mov_index_from_column_to_global





Use the parameters




lambda=0.01;




nits=200;




f=20;




and the initialisation




U = ones(f, nUsers);




Then driver ALS.m should run 200 iterations of ALS on the Movielens data. Hand in a plot of the value of g(U; M) as a function of the number of iterations, along with printouts of all the Matlab code you write.







Complete driver ALS.m with the following additional tasks:




{ For user 123, list the title, categories and rating of the 10 movies that the user ranked highest in the original ratings matrix. (Use getTitle.m as in movielens load data.m.)




{ For user 123, list the title and categories of the 10 movies that are most highly recommended for the user.










Hand in printouts of these lists. (Note that the highest predicted ratings are a bit higher than 5, which is not surprising because the model does not explicitly limit the predicted ratings between 0.5 and 5.)







 
Implementation of the preconditioned CG and GM-RES methods. (15 marks)




~

You are asked to implement the preconditioned CG and GMRES methods for A ~x = b in Matlab, and apply the methods to the 2D model problem (download build laplace 2D kron.m from Moodle).







Please carefully look over the answer of Question 3 of Tutorial 7 on preconditioning, and try out and experiment with the solution code, before starting this assignment question. Make sure you understand well what's going on in Question 3 of Tutorial 7.




(a) Implement the following in m- le myGMRES.m:
























4
School of Mathematical Sciences Monash University










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




 
Performs a sequence of GMRES iterations on




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




 
of the residual is smaller than tol (relative




 
to the initial residual), or until the number of iterations




 
reaches maxit. `steps' is the number of steps




 
that were performed, and `res' is a vector with the




 
residual size after every interation (the 2-norm of




 
the residual vector).




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




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




 
Performs a sequence of right-preconditioned GMRES iterations




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




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




 
to the initial residual), or until the number of iterations




 
reaches maxit. `steps' is the number of steps




 
that were performed, and `res' is a vector with the




 
residual size after every interation (the 2-norm of




 
the residual vector).




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




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




 
Performs a sequence of preconditioned CG iterations




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




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




 
to the initial residual), or until the number of iterations




 
reaches maxit. `steps' is the number of steps




 
that were performed, and `res' is a vector with the




 
residual size after every interation (the 2-norm of




 
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






















5
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.




 
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.)




 
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 = N2, 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(n1=4).




Note: This means that the total work is O(n5=4) when using the SSOR precon-ditioner, and O(n3=2) 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