CSC338. Homework 5
Due Date: Wednesday Feburary 12, 9pm
Please see the guidelines at https://www.cs.toronto.edu/~lczhang/338/homework.html
What to Hand In
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 [ ]:
import math
import numpy as np
Question 1
Part (a) — 2 pt
Are permutation matrices positive definite? Include your explanation in your pdf writeup.
Part (b) — 2 pt
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.
Part (c) — 8 pt
Complete the function cholesky_factorize that returns the Cholesky factorization of a matrix, according to its docstring.
In [ ]:
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.]])
“””
Question 2
Part (a) — 4 pts
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 [ ]:
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.])
“””
Part (b) — 2 pt
Explain why using solve_rank_one_update does not give us accurate results in the below example:
In [ ]:
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)
Part (c) — 4 pt
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 [ ]:
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.
“””
Question 3
Part (a) — 4 pt
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.
Part (b) — 4 pt
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!