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 4

$\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 4

Please answer the following questions in complete sentences in a clearly prepared manuscript and submit the solution by the end of class on Thursday, October 25th, 2012. You should submit your PDF file to Blackboard http://mycourses.purdue.edu as well as your Matlab files.

Version 3: Fixed the $A_1, L_1$ in problem 1, fixed the negative signs in Problem 4, clarified Problem 5.

Version 2: Fixed the problem numbers (5, 6), fixed the $\vx_2$ symbol in problem 3.

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: Tridiagonal systems

At this point in the class, you should be able to solve the following problems. These will return to haunt us when we study iterative methods. You’ll see!

Let $\mA$ be a symmetric, tridiagonal, positive definite matrix:

$$\mA = \bmat{\alpha_1 & \beta_1 \\ \beta_1 & \alpha_2 & \ddots \\ & \ddots & \ddots & \ddots \\ & & \beta_{n-2} & \alpha_{n-1} & \beta_{n-1} \\ & & & \beta_{n-1} & \alpha_n }$$

1. Let $\mA_{n-1}$ is the $n-1 \times n-1$ leading prinicipal minor of $\mA$. Suppose you are given the Choleksy factorization of $\mA_{n-1} = \mL_{n-1} \mL_{n-1}^T$. Determine the structure of $\mL_{n-1}$ (hint, $\mL_{n-1} \mL_{n-1}^T$ must be tridiagonal!).

2. Now, show how to compute the Cholesky factorization of $\mA = \mL \mL^T$ from $\mL_{n-1}$ in a constant number of operations.

3. Use the result of the first two parts to write down an algorithm that computes the Cholesky factorization of $\mA$ in a linear amount of work (starting from scratch!)

Problem 2: Tridiagonal least-squares problems

Let $\mA$ be a tridiagonal, $n+1 \times n$, “almost-symmetric” matrix:

$$\mA = \bmat{\alpha_1 & \beta_1 \\ \beta_1 & \alpha_2 & \ddots \\ & \ddots & \ddots & \ddots \\ & & \beta_{n-2} & \alpha_{n-1} & \beta_{n-1} \\ & & & \beta_{n-1} & \alpha_n \\ & & & & \beta_n}.$$

(For those who don’t appreciate mathematically imprecise terms, “almost-symmetric” means that it’s symmetric, except for the fact that it’s a rectangular matrix!)

1. Implement a Matlab algorithm to solve a least-squares problem with such a matrix using Given’s rotations.

2. Compare the performance of Matlab’s backslash operator to your Givens function on matrices of size $10,20,50,100,200,$ and $500$. Repeat each timing measurement at least 5 times.

Problem 3: The Schur Complement

Suppose we wish to solve

$$\mM \vx = \vb$$

and further suppose that you KNOW some of the values of $\vx$.

Let us permute and partition $\mM$ to be a block system:

$$\mM \vx = \bmat{ \mA & \mB \\ \mC & \mD } \bmat{\vx_1 \\ \vx_2} = \bmat{ \vb_1 \\ \vb_2 }$$

where $\vx_1$ is what you know.

1. Show how to solve for $\vx_2$ given $\vx_1$. Under what conditions is this possible?

2. Now, suppose that you have a very kind, but very confused dog that happened to eat your flash memory stick holding the piece of $\vx_1$ that you knew. However, you had saved your computed $\vx_2$ on your Purdue account, and so you have a backup. (This means you can assume that computing $vx_2$ from $\vx_1$ is possible for this problem if you determined it wasn’t always possible above.) Can you get $\vx_1$ back?

3. Combine these two parts to derive a single linear system to compute $\vx_1$ without computing $\vx_2$. The system you’ll derive is called the Schur complement.

Problem 4: Sparse matrices in Matlab

If you wish to use another language, please see me

In this problem, we’ll return to our favorite linear system, Poisson’s equation. (I told you that we’d see a lot of this system.)

We’ll assume that the boundary conditions are zero everywhere. Consequently, we’ll use the form of the system we derived on the exam (without all the extra rows!). Some of you did this already, so you may have an easier time on this question. That is, for the $N=3$ case, you’d have the $4 \times 4$ linear system from the exam. For the $N=4$ case, we have:

$$\bmat{-4 & 1 & & 1 & & & & & \\ 1 & -4 & 1 & & 1 & & & & \\ & 1 & -4 & & & 1 & & & \\ 1 & & & -4 & 1 & & 1 & & \\ & 1 & & 1 & -4 & 1 & & 1 & \\ & & 1 & & 1 & -4 & & & 1 \\ & & & 1 & & & -4 & 1 & \\ & & & & 1 & & 1 & -4 & 1 \\ & & & & & 1 & & 1 & -4 \\ } \bmat{ U(1,1) \\ U(1,2) \\ U(1,3) \\ U(2,1) \\ U(2,2) \\ U(2,3) \\ U(3,1) \\ U(3,2) \\ U(3,3) } = \bmat{ F(1,1) \\ F(1,2) \\ F(1,3) \\ F(2,1) \\ F(2,2) \\ F(2,3) \\ F(3,1) \\ F(3,2) \\ F(3,3) \\ }$$

where $U(i,j) = u(x_i,y_j)$ and $F(i,j) = f(x_i,y_j)$.

1. You can change an entry of a sparse matrix in Matlab with $A(i,j) = 1$. Write a Matlab routine to create the matrix above for a general $N$ by changing each element one at a time.

   A = sparse((N-1)*(N-1),N-1)*(N-1)); % this creates an all-zero sparse matrix
   <fill in the rest>

2. As shown in class, another way to create a sparse matrix is to use the coordinate form. Create the arrays i,j,andv for the matrix above. Hint there are $(N-1)^2 + 4(N-1)(N-2)$ non-zeros. Use this information to initialize the arrays i,j, and v. Write a routine to create the matrix above for a general $N$ using the coordinate form

   nz = (N-1)^2 + 4*(N-1)*(N-2);
   i = zeros(nz,1);
   j = zeros(nz,1);
   v = zeros(nz,1);
   <fill in the rest>
   A = sparse(i,j,v,(N-1)*(N-1),(N-1)*(N-1));

3. Which way is faster? Try creating matrices with $N = 10, 50, 100, 500$.

4. Is solving the linear system using backslash faster or slower with a sparse matrix? Does it depend on the value of $N$? Investigate and report.

5. Implement the Jacobi method for this problem. How many iterations does it take to converge to a relative residual of $10^{-4}$?

Problem 5: Prove or disprove

People seem to like these questions!

1. The eigenvalues of an $n \times n$ real-valued matrix are always real.

2. The solution to $(\mA^T \mA + \gamma I) \vx = \vb$ is unique for any $\gamma > 0$.

3. Every $n \times n$ matrix can be written as the sum of two non-singular matrices.

4. An $n \times n$ real-valued matrix has an orthogonal set of eigenvectors.

Problem 6: An experiment

As I’ve tried to emphasize in this class, matrix methods have broad applicability. In this problem, I’d like you to use what you’ve learned in this class to describe a problem in your own research in matrix terms.

If you don’t have one, spend 20 minutes searching on Google for others who have worked in your area with matrix methods and write about what they do.

If you cannot find anything, please document what you’ve done in an attempt to find something. But better yet, drop by office hours and ask for help finding something. Please try not to email us asking for help on this question as we’d like to discuss your area with you to help guide your search.

In any of these three cases, I’d expect somewhere between 2-3 paragraphs with equations describing what you’ve found (or failed to find … ) Some great things to mention are terms that you haven’t seen in class yet, or ideas that you’d like to learn more about, or anything (matrix related) you find interesting!