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 7

$\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 7 - Practice final and one question

Please answer the following questions in complete sentences in a clearly prepared manuscript and submit the solution by midnight on Sunday, December 2nd, 2012. You should submit your PDF file to Blackboard http://mycourses.purdue.edu as well as your Matlab files.

Version 1: 2012-11-28: Added a much needed missing word to problem 5, fixed typo in 3 with the last coefficient.

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: Matrix operations

Suppose we have real $n\times n$ matrices $\mC, \mD, \mE, \mF$. Show how to compute real $n \times n$ matrices $\mA$ and $\mB$ so that

$$\mA + i \mB = (\mC + i \mD)(\mE + i \mF)$$

with just three real $n\times n$ matrix-matrix products.

Problem 2: Matrix operations part 2

Let $\mA$ be symmetric. Prove the polarization identity

$$\vp^T \mA \vq = \frac{1}{4} \left[ (\vp + \vq)^T \mA (\vp + \vq) - (\vp - \vq)^T \mA (\vp - \vq) \right].$$

Problem 3: LU

Let

$$\mT = \bmat{\alpha_1 & \beta_2 \\ \beta_2 & \alpha_2 & \ddots \\ & \ddots & \ddots & \beta_{n} \\ & & \beta_{n} & \alpha_n }.$$

Determine the LU factorization of $\mT$.

Problem 4: Property A

Consider the Jacobi iteration on a matrix with property A. Show that the eigenvalues of the iteration matrix $\mM^{-1} \mN$ come in positive, negative pairs. That is, if $\mu$ is an eigenvalue of $\mM^{-1} \mN$, then $-\mu$ is also an eigenvalue.

Problem 5: Eigenvectors

Consider the shifted power method for finding the largest eigenvalue of $\mA$, i.e.

$$\vx\itn{k+1} = (\mA - \mu \mI) \vx\itn{k} \text{ normalized after each iteration }.$$

Let $\mA$ have real-eigenvalues and suppose you know all the eigenvalues except for the largest. Find the value of a shift $\mu$ in $\mA - \mu \mI$ that will make the power method converge as fast as possible.

Problem 6: Iterative methods

Show that the Jacobi iteration converges for $2 \times 2$ symmetric positive definite systems.

Problem 7: Krylov-ish methods

(I probably wouldn’t have used this one on the exam, but I like it as a homework question that gets at material I might put on the exam.)

Consider the algorithm (10.2.1) in Golub and van Loan. Show that $\vr_i^T \vr_{i-1} = 0$.

Problem 8: Iterative methods & GMRES

This is not a final question. It’s based on an assignment from Michael Saunders’ iterative methods class.

Download the poisson2D.mat file from Tim Davis’s sparse matrix repository. If you type load poisson2D in Matlab, then you’ll get a new variable called Problem where the linear systems is: Problem.A and Problem.b. Use the following code to setup $\mA$ and $\vb$ for Matlab.

load poisson2D
A = Problem.A;
b = Problem.b;
A = (A+A')/2; % make it symmetric

The matrix A is symmetric-positive definite. But we are going to look at solving $(\mA - \sigma \mI) \vx = \vb$ where $\sigma = 1.7 \times 10^{-2}$. We’ll declare a method converged if the relative residual has $2$-norm less than $10^{-8}$.

1. Plot the relative residuals for 100 iterations of the MINRES method (minres in Matlab) and the GMRES method (gmres in Matlab) restarted every 20 iterations on a semilogy scale. Describe what you observe, in particular, do the methods converge?

2. If you use an incomplete cholesky ichol or incomplete LU ilu preconditioner with GMRES, it will converge. How many iterations does it take? Does MINRES benefit from an incomplete Cholesky preconditioner?

3. Show the eigenvalues of the matrix before and after preconditioning. Do you observe any clustering?