{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CSC338. Homework 5\n",
    "\n",
    "Due Date: Wednesday Feburary 12, 9pm\n",
    "\n",
    "Please see the guidelines at https://www.cs.toronto.edu/~lczhang/338/homework.html\n",
    "\n",
    "### What to Hand In\n",
    "\n",
    "Please hand in 2 files:\n",
    "\n",
    "- Python File containing all your code, named `hw5.py`.\n",
    "- 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.\n",
    "\n",
    "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.\n",
    "\n",
    "**Make sure to remove or comment out all matplotlib or other expensive code before submitting your homework!**\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1\n",
    "\n",
    "### Part (a) -- 2 pt\n",
    "\n",
    "Are permutation matrices positive definite?\n",
    "Include your explanation in your pdf writeup.\n",
    "\n",
    "### Part (b) -- 2 pt\n",
    "\n",
    "Suppose that $A$ is a symmetric positive definite matrix. Show that the function\n",
    "\n",
    "$$ ||x||_A = (x^T A x)^{\\frac{1}{2}} $$\n",
    "\n",
    "on a vector $x$ satisfies the three properties of a vector norm.\n",
    "\n",
    "Include your solution in your pdf writeup.\n",
    "\n",
    "### Part (c) -- 8 pt\n",
    "\n",
    "Complete the function `cholesky_factorize` that returns the Cholesky\n",
    "factorization of a matrix, according to its docstring."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cholesky_factorize(A):\n",
    "    \"\"\"Return the Cholesky Factorization L of A, where\n",
    "        * A is an nxn symmetric, positive definite matrix\n",
    "        * L is lower triangular, with positive diagonal entries\n",
    "        * $A = LL^T$\n",
    "        \n",
    "    >>> M = np.array([[8., 3., 2.],\n",
    "                      [3., 5., 1.],\n",
    "                      [2., 1., 3.]])\n",
    "    >>> L = cholesky_factorize(M)\n",
    "    >>> np.matmul(L, L.T)\n",
    "    array([[8., 3., 2.],\n",
    "           [3., 5., 1.],\n",
    "           [2., 1., 3.]])\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2\n",
    "\n",
    "### Part (a) -- 4 pts\n",
    "\n",
    "Complete the function `solve_rank_one_update` that solves the system\n",
    "$(A - {\\bf u}{\\bf v}^T){\\bf x} = {\\bf b}$, assuming that the factorization\n",
    "$A = LU$ has already been done for you. You are welcome to add any helper\n",
    "functions that you wish, including functions that you wrote in\n",
    "homeworks 3 and 4. Just make sure that you include the helper functions\n",
    "in your python script submission."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_rank_one_update(L, U, b, u, v):\n",
    "    \"\"\"Return the solution x to the system (A - u v^T)x = b, where\n",
    "    A = LU, using the approach we derived in class using\n",
    "    the Sherman Morrison formula. You may assume that\n",
    "    the LU factorization of A has already been computed for you, and\n",
    "    that the parameters of the function have:\n",
    "        * L is an invertible nxn lower triangular matrix\n",
    "        * U is an invertible nxn upper triangular matrix\n",
    "        * b is a vector of size n\n",
    "        * u and b are also vectors of size n\n",
    "\n",
    "    >>> A = np.array([[2., 0., 1.],\n",
    "                      [1., 1., 0.],\n",
    "                      [2., 1., 2.]])\n",
    "    >>> L, U = lu_factorize(A) # from homework 3\n",
    "    >>> L\n",
    "    array([[1. , 0. , 0. ],\n",
    "           [0.5, 1. , 0. ],\n",
    "           [1. , 1. , 1. ]])\n",
    "    >>> U\n",
    "    array([[ 2. ,  0. ,  1. ],\n",
    "           [ 0. ,  1. , -0.5],\n",
    "           [ 0. ,  0. ,  1.5]])\n",
    "    >>> b = np.array([1., 1., 0.])\n",
    "    >>> u = np.array([1., 0., 0.])\n",
    "    >>> v = np.array([0., 2., 0.])\n",
    "    >>> x = solve_rank_one_update(L, U, b, u, v)\n",
    "    >>> x\n",
    "    array([1. , 0. , -1.])\n",
    "    >>> np.matmul((A - np.outer(u, v)), x)\n",
    "    array([1. , 1. , 0.])\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (b) -- 2 pt\n",
    "\n",
    "Explain why using `solve_rank_one_update` does not give us accurate results\n",
    "in the below example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_example():\n",
    "    A = np.array([[2., 0., 1.],\n",
    "                  [1., 1., 0.],\n",
    "                  [1., 1., 1.]])\n",
    "    L = np.array([[1., 0., 0.],\n",
    "                  [0.5, 1., 0.],\n",
    "                  [0.5, 1., 1.]])\n",
    "    U = np.array([[2., 0., 1.],\n",
    "                  [0., 1., -0.5],\n",
    "                  [0., 0., 1.]])\n",
    "    b = np.array([1, 1, -1])\n",
    "    u = np.array([0, 0, 0.9999999999999999])\n",
    "    v = np.array([0, 0, 0.9999999999999999])\n",
    "    x = solve_rank_one_update(L, U, b, u, v)\n",
    "    print(np.matmul((A - np.outer(u, v)), x) - b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (c) -- 4 pt\n",
    "\n",
    "Write function `solve_rank_one_update_iterative` that, given an\n",
    "approximation $x$ to the solution $(A - {\\bf u}{\\bf v}^T){\\bf x} = {\\bf b}$,\n",
    "performs one iteration of iterative refinement to obtain a better\n",
    "estimate $x^\\star$. You may use `solve_rank_one_update` as a helper\n",
    "function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_rank_one_update_iterative(L, U, b, u, v, x):\n",
    "    \"\"\"Return a better solution x* to the system (A - u v^T)x = b,\n",
    "    where A = LU. The first 5 parameters are the same as those of the \n",
    "    function `solve_rank_one_update`. The last parameter is an \n",
    "    estimate `x` of the solution.\n",
    "\n",
    "    This function should perform exactly *one* iterative refinement\n",
    "    iteration.\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3\n",
    "\n",
    "### Part (a) -- 4 pt\n",
    "\n",
    "We proved the Sherman-Morrison formula in Tutorial 5:\n",
    "    \n",
    "$$(A - uv^T)^{-1} = A^{-1} + A^{-1}u(1-v^TA^{-1}u)^{-1}v^T A^{-1}$$\n",
    "\n",
    "A related formula is the Woodbury formula for a rank-k update of the matrix $A$.\n",
    "Specifically, if $U$ and $V$ are $n \\times k$ matrices, then\n",
    "\n",
    "$$(A - UV^T)^{-1} = A^{-1} + A^{-1}U(I-V^TA^{-1}U)^{-1}V^T A^{-1}$$\n",
    "\n",
    "Prove that the Woodbury formula is correct.\n",
    "Include your solution in your pdf writeup.\n",
    "\n",
    "### Part (b) -- 4 pt\n",
    "\n",
    "A rank 1 matrix has the form $xy^T$ where $x$ and $y$ are column vectors.\n",
    "Suppose $A$ and $B$ are non-sigular matrices. Show that $A - B$ is rank 1\n",
    "if and only if $A^{-1} - B^{-1}$ is also rank 1. \n",
    "\n",
    "Include your solution in your pdf writeup.\n",
    "\n",
    "Hint: Use the Sherman-Morrison formula. This question is not supposed to be\n",
    "easy, so leave aside time to think!"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}