Due Date: Wednesday Feburary 12, 9pm

Please see the guidelines at https://www.cs.toronto.edu/~lczhang/338/homework.html

Please hand in 2 files:

- Python File containing all your code, named
`hw5.py`

. - PDF file named
`hw5_written.pdf`

containing your solutions to the written parts of the assignment. Your solution can be hand-written, but must be legible. Graders may deduct marks for illegible or poorly presented solutions.

If you are using Jupyter Notebook to complete the work, your notebook can be exported as a .py file (File -> Download As -> Python). Your code will be auto-graded using Python 3.6, so please make sure that your code runs. There will be a 20% penalty if you need a remark due to small issues that renders your code untestable.

**Make sure to remove or comment out all matplotlib or other expensive code before submitting your homework!**

Submit the assignment on **MarkUs** by 9pm on the due date. See the syllabus for the course policy regarding late assignments. All assignments must be done individually.

In [1]:

```
import math
import numpy as np
```

Are permutation matrices positive definite? Include your explanation in your pdf writeup.

Suppose that $A$ is a symmetric positive definite matrix. Show that the function

$$ ||x||_A = (x^T A x)^{\frac{1}{2}} $$

on a vector $x$ satisfies the three properties of a vector norm.

Include your solution in your pdf writeup.

Complete the function `cholesky_factorize`

that returns the Cholesky
factorization of a matrix, according to its docstring.

In [2]:

```
def cholesky_factorize(A):
"""Return the Cholesky Factorization L of A, where
* A is an nxn symmetric, positive definite matrix
* L is lower triangular, with positive diagonal entries
* $A = LL^T$
>>> M = np.array([[8., 3., 2.],
[3., 5., 1.],
[2., 1., 3.]])
>>> L = cholesky_factorize(M)
>>> np.matmul(L, L.T)
array([[8., 3., 2.],
[3., 5., 1.],
[2., 1., 3.]])
"""
```

Complete the function `solve_rank_one_update`

that solves the system
$(A - {\bf u}{\bf v}^T){\bf x} = {\bf b}$, assuming that the factorization
$A = LU$ has already been done for you. You are welcome to add any helper
functions that you wish, including functions that you wrote in
homeworks 3 and 4. Just make sure that you include the helper functions
in your python script submission.

In [3]:

```
def solve_rank_one_update(L, U, b, u, v):
"""Return the solution x to the system (A - u v^T)x = b, where
A = LU, using the approach we derived in class using
the Sherman Morrison formula. You may assume that
the LU factorization of A has already been computed for you, and
that the parameters of the function have:
* L is an invertible nxn lower triangular matrix
* U is an invertible nxn upper triangular matrix
* b is a vector of size n
* u and b are also vectors of size n
>>> A = np.array([[2., 0., 1.],
[1., 1., 0.],
[2., 1., 2.]])
>>> L, U = lu_factorize(A) # from homework 3
>>> L
array([[1. , 0. , 0. ],
[0.5, 1. , 0. ],
[1. , 1. , 1. ]])
>>> U
array([[ 2. , 0. , 1. ],
[ 0. , 1. , -0.5],
[ 0. , 0. , 1.5]])
>>> b = np.array([1., 1., 0.])
>>> u = np.array([1., 0., 0.])
>>> v = np.array([0., 2., 0.])
>>> x = solve_rank_one_update(L, U, b, u, v)
>>> x
array([1. , 0. , -1.])
>>> np.matmul((A - np.outer(u, v)), x)
array([1. , 1. , 0.])
"""
```

Explain why using `solve_rank_one_update`

does not give us accurate results
in the below example:

In [4]:

```
def run_example():
A = np.array([[2., 0., 1.],
[1., 1., 0.],
[1., 1., 1.]])
L = np.array([[1., 0., 0.],
[0.5, 1., 0.],
[0.5, 1., 1.]])
U = np.array([[2., 0., 1.],
[0., 1., -0.5],
[0., 0., 1.]])
b = np.array([1, 1, -1])
u = np.array([0, 0, 0.9999999999999999])
v = np.array([0, 0, 0.9999999999999999])
x = solve_rank_one_update(L, U, b, u, v)
print(np.matmul((A - np.outer(u, v)), x) - b)
```

Write function `solve_rank_one_update_iterative`

that, given an
approximation $x$ to the solution $(A - {\bf u}{\bf v}^T){\bf x} = {\bf b}$,
performs one iteration of iterative refinement to obtain a better
estimate $x^\star$. You may use `solve_rank_one_update`

as a helper
function.

In [5]:

```
def solve_rank_one_update_iterative(L, U, b, u, v, x):
"""Return a better solution x* to the system (A - u v^T)x = b,
where A = LU. The first 5 parameters are the same as those of the
function `solve_rank_one_update`. The last parameter is an
estimate `x` of the solution.
This function should perform exactly *one* iterative refinement
iteration.
"""
```

We proved the Sherman-Morrison formula in Tutorial 5:

$$(A - uv^T)^{-1} = A^{-1} + A^{-1}u(1-v^TA^{-1}u)^{-1}v^T A^{-1}$$

A related formula is the Woodbury formula for a rank-k update of the matrix $A$. Specifically, if $U$ and $V$ are $n \times k$ matrices, then

$$(A - UV^T)^{-1} = A^{-1} + A^{-1}U(I-V^TA^{-1}U)^{-1}V^T A^{-1}$$

Prove that the Woodbury formula is correct. Include your solution in your pdf writeup.

A rank 1 matrix has the form $xy^T$ where $x$ and $y$ are column vectors. Suppose $A$ and $B$ are non-sigular matrices. Show that $A - B$ is rank 1 if and only if $A^{-1} - B^{-1}$ is also rank 1.

Include your solution in your pdf writeup.

Hint: Use the Sherman-Morrison formula. This question is not supposed to be easy, so leave aside time to think!