Matrix computations &
Numerical linear algebra

David Gleich

Purdue University

Fall 2012

Course number CS-51500

Tuesday and Thursday, 10:30-11:45am

Location Lawson 1106

Homework 5

$\newcommand{\mat}[1]{\boldsymbol{#1}} \renewcommand{\vec}[1]{\boldsymbol{\mathrm{#1}}} \newcommand{\vecalt}[1]{\boldsymbol{#1}} \newcommand{\conj}[1]{\overline{#1}} \newcommand{\normof}[1]{\|#1\|} \newcommand{\onormof}[2]{\|#1\|_{#2}} \newcommand{\itr}[2]{#1^{(#2)}} \newcommand{\itn}[1]{^{(#1)}} \newcommand{\eps}{\varepsilon} \newcommand{\kron}{\otimes} \DeclareMathOperator{\diag}{diag} \DeclareMathOperator{\trace}{trace} \DeclareMathOperator{\tvec}{vec} \newcommand{\prob}{\mathbb{P}} \newcommand{\probof}[1]{\prob\left\{ #1 \right\}} \newcommand{\pmat}[1]{\begin{pmatrix} #1 \end{pmatrix}} \newcommand{\bmat}[1]{\begin{bmatrix} #1 \end{bmatrix}} \newcommand{\spmat}[1]{\left(\begin{smallmatrix} #1 \end{smallmatrix}\right)} \newcommand{\sbmat}[1]{\left[\begin{smallmatrix} #1 \end{smallmatrix}\right]} \newcommand{\RR}{\mathbb{R}} \newcommand{\CC}{\mathbb{C}} \newcommand{\eye}{\mat{I}} \newcommand{\mA}{\mat{A}} \newcommand{\mB}{\mat{B}} \newcommand{\mC}{\mat{C}} \newcommand{\mD}{\mat{D}} \newcommand{\mE}{\mat{E}} \newcommand{\mF}{\mat{F}} \newcommand{\mG}{\mat{G}} \newcommand{\mH}{\mat{H}} \newcommand{\mI}{\mat{I}} \newcommand{\mJ}{\mat{J}} \newcommand{\mK}{\mat{K}} \newcommand{\mL}{\mat{L}} \newcommand{\mM}{\mat{M}} \newcommand{\mN}{\mat{N}} \newcommand{\mO}{\mat{O}} \newcommand{\mP}{\mat{P}} \newcommand{\mQ}{\mat{Q}} \newcommand{\mR}{\mat{R}} \newcommand{\mS}{\mat{S}} \newcommand{\mT}{\mat{T}} \newcommand{\mU}{\mat{U}} \newcommand{\mV}{\mat{V}} \newcommand{\mW}{\mat{W}} \newcommand{\mX}{\mat{X}} \newcommand{\mY}{\mat{Y}} \newcommand{\mZ}{\mat{Z}} \newcommand{\mLambda}{\mat{\Lambda}} \newcommand{\mPbar}{\bar{\mP}} \newcommand{\ones}{\vec{e}} \newcommand{\va}{\vec{a}} \newcommand{\vb}{\vec{b}} \newcommand{\vc}{\vec{c}} \newcommand{\vd}{\vec{d}} \newcommand{\ve}{\vec{e}} \newcommand{\vf}{\vec{f}} \newcommand{\vg}{\vec{g}} \newcommand{\vh}{\vec{h}} \newcommand{\vi}{\vec{i}} \newcommand{\vj}{\vec{j}} \newcommand{\vk}{\vec{k}} \newcommand{\vl}{\vec{l}} \newcommand{\vm}{\vec{l}} \newcommand{\vn}{\vec{n}} \newcommand{\vo}{\vec{o}} \newcommand{\vp}{\vec{p}} \newcommand{\vq}{\vec{q}} \newcommand{\vr}{\vec{r}} \newcommand{\vs}{\vec{s}} \newcommand{\vt}{\vec{t}} \newcommand{\vu}{\vec{u}} \newcommand{\vv}{\vec{v}} \newcommand{\vw}{\vec{w}} \newcommand{\vx}{\vec{x}} \newcommand{\vy}{\vec{y}} \newcommand{\vz}{\vec{z}} \newcommand{\vpi}{\vecalt{\pi}}$

Homework 5

Please answer the following questions in complete sentences in a clearly prepared manuscript and submit the solution by 5pm on Friday, November 9th, 2012. You should submit your PDF file to Blackboard http://mycourses.purdue.edu as well as your Matlab files.

Version 2: Posted 11-5 at 6:26pm, fixed typo in linsys_sor.m and removed unneeded line from problem 1-4.

Version 1: Posted 10-27 at 11:15am, fixed indices in problem 1-4; fixed “write implement” typo.

Problem 0: List your collaborators.

Please identify anyone, whether or not they are in the class, with whom you discussed your homework. This problem is worth 1 point, but on a multiplicative scale.

Problem 1: The power method, and beyond!

We’ll show that the humble power-method is really rather more interested than I alluded to in class! It’s a flexible starting point that serves as the basis for almost all of the eigenvalue algorithms.

1. (Warm up) Consider the power method as a Matlab code:

   maxit = 1000;
   for i=1:maxit
    x = A*x;
    x = x/norm(x);
    lam = x'*A*x;
   end

   Let $\vx\itn{i}$ be the value of x at the start of the $i$th instance of this loop. And let $\lambda\itn{i}$ be the value of lam at the start of the $i$th iteration. Suppose that $\vx$ and $\lambda$ are the true eigenvector, eigenvalue pair that we are converging too. If $\normof{\vx - \vx\itn{i}} = \eps$, show that: $|\lambda - \lambda^{i}|$ is $O(\eps^2)$ (where the constant may depend on the matrix $A$).

2. (Warm up) Show that if $\vx$ is an eigenvector of $\mA$ with eigenvalue $\lambda$, then $\vx$ is an eigenvector of $\mA^{-1}$. What is the eigenvalue?

3. Let $\mA$ have a complete set of eigenvalues $\lambda_i \not= \lambda_j$, for all $i \not= j$ and associated eigenvectors $\vx_i$. For the sake of simplicity, let’s suppose that all the eigenvalues are real. Suppose we run the following Matlab code:

   maxit = 1000;
   [L,U,P] = lu(A);
   for i=1:maxit
    x = U\(L\(P*x));
    x = x/norm(x);
    lam = x'*A*x;
   end

   Does this converge? If so, what does it converge to, and how fast does it converge?

4. Suppose you were in a parallel linear algebra class, taught by another professor, and you were answering the following question:

   Make the power method converge faster to any eigenvalue, eigenvector pair!
   Dense linear algebra is cheap. Work in pairs.

   Your friend says that you should look at this algorithm:

   maxit = 1000;
   for i=1:maxit
    x = (A - lam*eye(size(A,1)))\x);
    x = x/norm(x);
    lam = x'*A*x;
   end

   This algorithm is called “Rayleigh Quotient Iteration” and the standard result is that $x\itn{k}$ converges to an eigenvector cubically. That is, of $\normof{\vx\itn{k} - \vx} = \eps$, then $\normof{\vx\itn{k+1} - \vx} = O(\eps^3)$. Look up this result in a textbook, or on the web, and explain the proof in your own words.

5. Implement this method. Do you observe cubic convergence if you try to find an eigenvalue, eigenvector pair for a symmetric matrix?

Problem 2: Gershgorin disks

Theorem 7.2.1 in Golub and van Loan states the Gershgorin circle theorem. This theorem provides an easy check to find regions containing the eigenvalues of a matrix.

1. Read this theorem. Explain it in your own words.

2. Write a function to plot the Gershgorin circles in Matlab.

3. Plot the Gershgorin circles for the matrix:

   n = 10; 
   on = ones(n,1); 
   A = spdiags([-2*on 4*on -2*on],-1:1,n,n);

Problem 3: PageRank and the power method

In class, I derived PageRank as the solution of the linear system:

$$(\eye - \alpha \mP) \vx = (1-\alpha) \vv$$

where $\mP$ is a column-stochastic matrix, $0 \le \alpha < 1$ and $\vv$ is a non-negative vector whose elements sum to 1.

1. (Warm up) Use Gershgorin circle’s to show that $(\eye - \alpha \mP)$ is a non-singular matrix.

2. The standard definition of PageRank is as the largest eigenvector of the matrix:

   $$\mM = \alpha \mP + (1-\alpha) \vv \ve^T$$

   normalized to have sum one (instead of the standard 2-norm scaling). Consider the power method, without normalization, applied to $\mM$ starting from the vector $\vv$:

   $$\vx\itn{k+1} = \mM \vx\itn{k} = \mM^{k+1} \vv.$$

   Show that the iterates from this method are always non-negative and sum to 1.

3. Use this fact to simplify the iteration and show that the power method will converge to the solution of the linear system $(\eye - \alpha \mP) \vx = (1-\alpha) \vv$. You should get the same iteration that I discussed in class.

4. Download the data purdue_web.mat from http://www.cs.purdue.edu/homes/dgleich/cs515/homeworks/purdue_web.mat Load the data. It contains the matrix $\mP$ for Purdue’s website from 2001. The url variable contains the URL for each row/column of the matrix. Use url(1) to show the first url.

5. Use the power method to compute the PageRank vector $\vx$ for $\mP$ using $\alpha = 0.85$, $v_i = 1/n$ where $n$ is the size of the matrix. In Matlab

   n = size(P,1);
   v = ones(n,1)./n;
   alpha = 0.85;

   What is x(1)? Use the following code to sort the vector and output the top 27 entries:

   [ignore p] = sort(x,'descend');
   urls(p(1:27))

   Report the output from this code.

Problem 4: Changing SOR’s omega

I botched the in-class demo of SOR because I was in a hurry and used the difference between iterations to check for convergence instead of the relative residual. My apologies for this problem!
But we’ll fix it in this homework problem.

Modify the SOR method that your professor wrote for class to check the residual for convergence instead of the difference in iterates.
We’ll look at the same Poisson problem we saw in class.

A = delsq(numgrid('S', 22));

Except this time, we’ll do it correctly – without that issue that arose during the class-demo. Prepare a plot on how many iterations it takes the SOR method to compute a solution with residual $10^{-6}$ for various values of $\omega$. Find the fastest $\omega$ you can!

Problem 5: Compare CG as described in class to the standard CG iteration

One of the most remarkable results in iterative methods is how straightforward many algorithms can be derived, and how these straightforward derivations give rise to beautifully simple implementations. In this problem you’ll explore this property with the conjugate gradient method.

1. Implement a CG method based on the Lanczos method. That is, implement the Lanczos method to compute the matrix $\mT_k$ and $\mV_k$. Then solve $\mT_k \vy\itn{k} = \ve_1$ and compute $\vx = \normof{\vb} \mV_k \vy_k$. Show the code for your implementation. At each step, you should compute the normalized residuals.

2. Compare the first 25 residuals from your Lanczos-based CG code with the a standard implementation of CG from: http://www.cs.purdue.edu/homes/dgleich/cs515/matlab/cg.m for the linear system

   n = 100; 
   on = ones(n,1); 
   A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
   b = ones(n,1);

3. Using the cg.m function, look at how many iterations are required for CG to converge to a tolerance of $10^{-8}$ for the matrix in the last part. Determine how this scales with $n$.