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 2

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

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

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.

Section 1: The SVD

Problem 1: Compute it!

Produce the analytic SVDs of the following matrices. (That is, no numerical approximations, but feel free to let Matlab give you a good guess!). It’s important to think about these questions because I may give similar questions on a midterm or final and you’ll be expected to remember these. It’s also handy to think about how to construct these, even though we haven’t seen any algorithms yet. You should be able to work them all out directly from the definition.

1. $\bmat{1 & 0 \\ 0 & 0}$

2. $\bmat{5 & 0 \\ 2 & 0}$

3. $\bmat{5 & -5 \\ 2 & -2 \\ 0 & 0}$

Problem 2: Fun with the SVD.

The following quesitons are some warm up questions about the theory of the SVD. Each of these is designed to teach you something about the SVD. So there is usually a nice fact that I want you to know.

1. What is $||\mQ||$ where $\mQ$ is an orthogonal matrix?

2. Find a matrix with two different SVD factorizations. (Remember, someone asked about this in class!) What is the same in both cases?

3. Show that $||\mA^{-1}|| = 1/\sigma_{\min}$.

Problem 3: Extracting rank from the SVD

In this problem, we’ll work through the details of how to show that the SVD encodes the rank of a matrix.

1. Show that the rank of a diagonal matrix is the number of non-zeros.

2. Show that the rank of $\mA \mB = \min(\rank \mA, \rank \mB)$ when $\mA$ and $\mB$ are full-rank.

3. Briefly explain that $\rank \mQ$ for a square orthogonal $\mQ$ is full rank.

4. Explain why we’ve determined that the rank of a matrix is given by the number of non-zero singular values.

Problem 4: The celebrated Eckart-Young theorem

One of the early communities to use the SVD for data analyais were the psychometricians. They were trying to understand how to analyze matrices of test scores!

If $N$ individuals are each subjected to $n$ tests, it is a fundamental postulate of factor theory that the resulting $n \times N$ score matrix a can be adequately approximated by another matrix whose rank $r$ is less than the smaller of $n$ or $N$.
Closely associated to this postulate is the purely mathematical problem of finding that matrix of rank $r$ which most closely approximates a given matrix of higher rank $R$. It will be shown that if the least-squares criterion of approximation be adopted, this problem has a general solution which is relatively simple in a theoretical sense, though the amount of numerical work involved in applications may be prohibitive. Certain conclusions can be drawn from the theoretical solution which may be of importance in practical work. – [Eckart and Young](http://dx.doi.org/10.1007/BF02288367)

In this question, we will work through a proof of the Eckart-Young theorem. You are welcome to go look at the paper for the original proof. However, I think the more modern steps below will be easier to follow. Note that you can also find this same proof in many references. Again, you are free to look at these and even copy them. However, I will ask that you explain everything! If you find yourself reading these and thinking, “hmm… why is that true…?” Those are exactly the steps we are expecting you to explain.

The Eckart-Young theorem states:

Let $\mA \in \RR^{m \times n}$, without loss of generality, assume $m \ge n$. Let $\mA = \mU \mat{\Sigma} \mV^T$ be the SVD. Let $\mU = \bmat{\vu_1, \ldots, \vu_m}$, $\mat{\Sigma} = \diag(\sigma_1,\ldots,\sigma_n)$, and $\mV = \bmat{\vv_1, \ldots, \vv_n}$. Define

$$\mA_k = \sum_{i=1}^k \sigma_i \vu_i \vv_i^T.$$

Then:

$$\normof{\mA - \mA_k} \le \normof{\mA - \mB}$$

for any matrix $\mB$ of rank $k$.

It’s helpful to translate this into real language. Suppose I give you a matrix $\mA$, and ask for the closest matrix to $\mA$ but with rank $k$ instead of rank $\mA$ instead. Eckart-Young tells us that we can construct such as matrix from the SVD.

In this problem, we’ll use the matrix 2-norm.

1. Explain why we don’t lose any generality with the assumption $m \ge n$.

2. Warmup Suppose that $\mA$ is a diagonal matrix with non-negative entries. What is the best diagonal rank $k$ approximation to $\mA$? (Yes, you have to prove it. No, you cannot assume the Eckart-Young theorem is true!)

3. Warmup Show that $\normof{\mA - \mA_n} = 0$.

4. Warmup Show that $\normof{\mA - \mA_R} = 0$ where $R$ is the rank of $\mA$. Why does this mean we can assume that $k < R$ without losing any further generality?

5. Show that $\normof{\mA - \mA_k} = \sigma_{k+1}$.

6. We’ll proceed by contradiction now. Suppose that a rank $k$ matrix $\mB$ exists with $\normof{\mA - \mB} < \sigma_{k+1}$. Briefly justify why:

   $$\normof{(\mA - \mB) \vx} < \sigma_{k+1} \normof{\vx}.$$

7. Using the last result, we can choose $\vx$ to be anything we want. Note that $\mB$ is rank $k$. This means it must have a null-space of dimension $n - k$. Show that, for any vector $\vw$ in the null-space of $\mB$ that:

   $$\normof{\mA \vw} < \sigma_{k+1} \normof{\vw}.$$

8. Now we have a subspace of $\RR^{n}$ such that $\normof{\mA \vw}$ is bounded above. We’ll now show that there’s another subspace where it’s bounded below! Let $\vz \in \text{span}(\vv_1, \ldots, \vv_{k+1})$. Show that $\normof{\mA \vz} \ge \sigma_{k+1} \normof{\vz}$. (Hint, note that $\mV_{k+1} = \bmat{\vv_1, \ldots, \vv_{k+1}}$ is an orthogonal basis for that span.)

9. What’s the dimension of the space where $\mA \vw$ is bounded above? What’s the dimension of the space where $\mA \vz$ is bounded below? Why does this lead to a contradiction?

10. Give yourself a pat on the back!

Problem 5: Eigenfaces

Download the matrix ‘Yale_32x32.mat’ from http://www.cad.zju.edu.cn/home/dengcai/Data/Yale/Yale_32x32.mat Recall, you can use

urlwrite('http://www.cad.zju.edu.cn/home/dengcai/Data/Yale/Yale_32x32.mat','Yale_32x32.mat');

as we did in the matlab/cosine_similarity.m file

1. Compute the singular value decomposition of the matrix of images as a “pixels $\times$ faces” matrix. What is the largest singular value?

2. Plot the leading singular vector $\vu_1$ as a $32\times 32$ pixel image. What do you see?

3. If you’re curious about this, this website has more info http://en.wikipedia.org/wiki/Eigenface. Note that it describes something slightly different from what we are doing in the first two parts. The key difference is that our matrix isn’t centered with respect to the images.

Section 2: Linear systems and QR

Problem 6: Backsolve

1. Implement a backsolve as a function in Matlab. Show your code.

2. Construct a random upper-triangular linear system via:

   A = triu(randn(n)); b = randn(n,1);

Compare the performance of your backsolve to Matlab’s \ method to solve a linear system.

Problem 7: Poisson’s equation

In this problem, we’ll meet one of the most common matrices studied in numerical linear algebra: the $2d$-Laplacian. We arrive at this matrix by discretizing a partial differential equation. Poisson’s equation is:

$$\Delta u = f$$

where $u(x,y)$ is a continuous function defined over the unit-plane (i.e. $0 \le x \le 1, 0 \le y \le 1$), $f(x,y)$ is a continuous function defined over the same region, and $\Delta$ is the Laplacian operator:

$$\Delta u = \partial^2 u/\partial x^2 + \partial^2 u/\partial y^2.$$

Given a function $f$, we want to find a function $u$ that satifies this equation. There are many approaches to solve this problem theoeretically and numerically. We’ll take a numerical approach here.

Suppose we discretize the function $u$ at regular points $x_0, \ldots, x_n$, and $y_1, \ldots, y_n$ where $x_i = y_i = i/n$ so that we have:

$$u(x,y) \approx \begin{matrix} u(x_0, y_0) & \cdots & u(x_n, y_0) \\ \vdots & \ddots & \vdots \\ u(x_0, y_n) & \cdots & u(x_n, y_n). \end{matrix}$$

For this discretization, note that

$$\begin{aligned} \Delta u(x_i, y_j) & \approx {n^2} \left( u(x_{i-1},y_{j}) - 2 u(x_i,y_j) + u(x_{i+1},y_j) \right) \\ & \qquad+ {n^2} \left( u(x_{i},y_{j-1}) - 2 u(x_i,y_j) + u(x_{i},y_{j+1}) \right) \\ & = {n^2} \left( u(x_{i-1},y_{j})+ u(x_{i},y_{j-1}) - 4 u(x_i,y_j) + u(x_{i+1},y_j) + u(x_{i},y_{j+1}) \right) \\ & = f(x_i, x_j). \end{aligned}$$

What we’ve done here is use the approximation:

$$\partial^2 u / \partial x^2 \approx \frac{1}{h^2} \left( u(x-h) - 2 u(x) + u(x+h) \right)$$

for both partial terms.
We need this equation to hold at each point $x_i, y_j$. But note that there are some issues with this equation at the boundary values (where x=0 or 1, or where y=0 or 1).
For this problem, we’ll make it very simple and set:

$$u(0,y_j) = u(1,y_j) = u(x_i,0) = u(x_i,1) = 0.$$

Now, we’ll do what we always do! Turn this into some type of matrix equation!

Let $\mU$ be an $n+1 \times n+1$ matrix that we’ll index from zero instead of one:

$$\mU = \bmat{U_{0,0} & \cdots & U_{0,n} \\ \vdots & \ddots & \vdots \\ U_{n,0} & \cdots & U_{n,n} }.$$

where $U_{i,j} = u(x_i, y_j)$. At this point, we are nearly done. What we are going to do is turn Poisson’s equation into a linear system. This will be somewhat like how we turned image resampling into a matrix vector equation in the last homework.

In order to write $\mU$ as a vector, we’ll keep the convention from last time:

$$\mU = \bmat{ u_1 & \cdots & u_{n+1} \\ u_{n+2} & \cdots & u_{2(n+1)} \\ \vdots & \ddots & \vdots \\ u_{n(n+1) + 1} & \cdots & u_{(n+1)(n+1)} }.$$

Let $\vu$ be the vector of elements here. Note that our approximation to $\Delta u$, just involved a linear combination of the elements of $\vu$. This means we have a linear system:

$$\mA \vu = \vf$$

where $\mA$ and $\vf$ correspond to all the rows in:

$$\frac{1}{h^2} \left( u(x_{i-1},y_{j})+ u(x_{i},y_{j-1}) - 4 u(x_i,y_j) + u(x_{i+1},y_j) + u(x_{i},y_{j+1} \right) = f(x_i, y_j).$$

1. Let $n=3$. Write down the $16 \times 16$ linear equation for $\vu$ including all the boundary conditions. Note that you can encode the boundary conditions by adding a row of $\mA$ where: $u_i = 0$.

2. Write a Matlab code to construct the matrix $\mA$ and vector $\vf$ when $n=10$ and $f(x,y) = 1$. Here’s some pseudo-code to help out:

   n = 10;
   A = zeros((n+1)^2,(n+1)^2);
   f = zeros((n+1)^2,1);
   G = reshape(1:((n+1)^2), n+1, n+1)'; 
   for i=0:n
    for j=0:n
        row = G(i+1,j+1);
        if i==0 || j == 0 || i==n || j==n
            % we are on a boundary
            f(row) = 0;
            <fill in statement to set A(row,:)>
        else
            % we are NOT on a boundary
            f(row) = n^2;
            <fill in statement to set A(row,:)>
        end
    end
   end

3. Solve for $\vu$ using Matlab’s backslash solver show the result using the mesh function.

(REMOVED) Problem 8

I’ve removed problem 8

Problem 9: Bonus! (5 points)

1. There is no collaboration allowed on bonus problems. Please answer the following question:

I collaborated with ________ while working on this problem.

1. In class, we showed that the first column of a QR factorization of a column vector $\va$ is:

   $$\va = \underbrace{(\va/\normof{\va})}_{\mQ} \underbrace{\normof{\vq}}_{\mR}.$$

   Next, we showed that the full QR factorization is given by:

   $$\va = \mH \normof{\va} \ve_1,$$

   where $\mH$ is a Householder reflector. This requires that $\mH \ve_1 = \va/\normof{\va}$.

I want you to show that $\mH \ve_1 = \va/\normof{\va}$ directly from the definition of $\mH$, that is:

$$\mH = \eye - 2 \frac{\vv \vv^T}{\vv^T \vv}$$

where $\vv = \normof{\va}\ve_1 - \va$. I gave this one a few stabs myself and couldn’t come up with anything that wasn’t a tedious derivation. But I didn’t think about it that long. I’d love to see if there’s a slick proof of this result, again, directly from the defintion – and hence, this bonus question. To get full points, you can just slog it out based on the definition.