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 3

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

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 27th, 2012. You should submit your PDF file to Blackboard http://mycourses.purdue.edu as well as your Matlab files. Given the shorter due-date, feel free to inquire about extensions, although, realize you have a midterm on October 4th. I’ve also made the homework shorter.

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 Cholesky Factorization

Through this problem, assume that the matrices are all symmetric. You may assume this in your code as well.

1. Implement a Matlab function for the Cholesky factorization.

2. Research the Cholesky factorization and answer the following question, either providing a proof or a reference. Is the Cholesky factorization unique?

3. Compare your Cholesky factorization code to Matlab’s chol function on the following matrix.

    n = 10;
    Q = qr(randn(n));
    D = abs(randn(n,1));
    A = Q*diag(D)*Q';

4. Your professor boldy asserted that using the Cholesky factorization was faster than using the LU factorization. Make sure he’s still right. Show the performance of Matlab’s LU and Cholesky factorization for matrices with size that varies from 10 to 500.

5. Your professor also boldy asserted that we preserve positive definiteness after one step of the Cholesky factorization. Recall that this was:

   $$\mA = \bmat{\alpha & \vb^T \\ \vb & \mC} = \mF \mF^T = \bmat{\gamma & 0 \\ \vf & \mF_1} \bmat{\gamma & \vf^T \\ 0 & \mF_1^T}.$$

   So we set $\gamma = \sqrt{\alpha}$ and $\vf = \vb/\gamma$. Then, we recursively compute the Cholesky factorization of

   $$\mC - \vb \vb^T/\alpha = \mF_1 \mF_1^T.$$

   This only works if $\mC - \vb \vb^T/\alpha$ is positive definite. Prove that this is true. (Hint: you can do it from the definition of positive definiteness: $\vx^T \mA \vx > 0$ for all $\vx$.)

6. Briefly explain why step 5 justifies that a Cholesky factorization exists for any positive definite matrix. (This is really like a one or two sentence answer.)

7. One of the most useful aspects of the Choleskfy factorization is as a way to check if a matrix is positive definite or negative definite. A matrix is negative definite if

   $$\vx^T \mA \vx < 0$$

   for all $\vx$. Change your Matlab implementation to report if a matrix is not positive definite. (Hint: this relates to why a Cholesky factorization always exists for a positive definite matrix in the previous problem.) Use this to test if the matrix for Poisson’s equation from Homework 2, Problem 7 is positive or negative definite.

Problem 2: Stability analysis

1. Suppose that $x_1, \ldots, x_n$ are a sequence of samples from a random variable. The sample variance of $x_1, \ldots, x_n$ is:

   $$V = \frac{1}{n-1} \sum_{i=1}^n (x_i - \frac{1}{n} \sum_{i=1}^n x_i )^2.$$

   Compute the condition number of the sample variance.

2. This is a real problem that I ran into this summer! I was using a program written by someone else to compute a particular statistical estimate. The details don’t really matter, but at one step, the code needed to compute the Hurwitz zeta function:

   $$H(s,q) = \sum_{n=0}^\infty \frac{1}{(q+n)^s}$$

   where $s$ ranged from 1 to 7 and where $q$ ranged from 1 to 500. They were using Matlab, which doesn’t provide a function for the Hurtwitz zeta function, but does provide a function for the Riemann zeta function:

   $$R(s) = \sum_{n=1}^\infty \frac{1}{n^s}.$$

   Note that the only difference between these two functions is that the first term in the Riemann zeta function is $1/(1^s)$ and the first term in the Hurwitz zeta function is $1/(q^s)$. To account for this difference, the code I was using did the following:

     function h=hzeta(s,q)
     z = zeta(s)
     h = z - sum((1:(q-1)).^(-s));

   Was I happy when I realized this? You should use what you’ve learned about finite precision computations to answer this question. You can evaluate the Hurwitz zeta function to arbitrary precision via Wolfram Alpha. There is an answer that’s more correct on this problem, but the most critical part of your answer is how you justify it.