{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CSC338. Homework 3\n",
    "\n",
    "Due Date: Wednesday January 29, 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 `hw3.py`.\n",
    "- PDF file named `hw3_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) -- 4 pt\n",
    "\n",
    "Solve the following system $A{\\bf x} = {\\bf b}$ using Gauss Elimination **by hand**. Show all your steps.\n",
    "Show all your steps in your pdf writeup.\n",
    "\n",
    "\\begin{align*}\n",
    "A = \\begin{bmatrix}\n",
    "1 & 0 & -2 \\\\\n",
    "-1 & 1 & 1 \\\\\n",
    "0 & 2 & -1\n",
    "\\end{bmatrix}, \n",
    "b = \\begin{bmatrix}\n",
    "-1 \\\\\n",
    "2 \\\\\n",
    "3 \\\\\n",
    "\\end{bmatrix}, \n",
    "\\end{align*}\n",
    "\n",
    "### Part (b) -- 2 pt\n",
    "\n",
    "Using your work from part (a), find the LU factorization of $A$.\n",
    "Show all your steps in your pdf writeup.\n",
    "\n",
    "## Question 2\n",
    "\n",
    "For this question, you can refer to and use any of the code that we wrote together in tutorial 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code from tutorial 3\n",
    "\n",
    "def backward_substitution(A, b):\n",
    "    \"\"\"Return a vector x with np.matmul(A, x) == b, where \n",
    "        * A is an nxn numpy matrix that is upper-triangular and non-singular\n",
    "        * b is an nx1 numpy vector\n",
    "    \"\"\"\n",
    "    n = A.shape[0]\n",
    "    x = np.zeros_like(b, dtype=np.float)\n",
    "    for i in range(n-1, -1, -1):\n",
    "        s = 0\n",
    "        for j in range(n-1, i, -1):\n",
    "            s += A[i,j] * x[j]\n",
    "        x[i] = (b[i] - s) / A[i,i]\n",
    "    return x\n",
    "\n",
    "def eliminate(A, b, k):\n",
    "    \"\"\"Eliminate the k-th row of A, in the system np.matmul(A, x) == b,\n",
    "    so that A[i, k] = 0 for i < k. The elimination is done in place.\"\"\"\n",
    "    n = A.shape[0]\n",
    "    for i in range(k + 1, n):\n",
    "        m = A[i, k] / A[k, k]\n",
    "        for j in range(k, n):\n",
    "            A[i, j] = A[i, j] - m * A[k, j]\n",
    "        b[i] = b[i] - m * b[k]\n",
    "\n",
    "def gauss_elimination(A, b):\n",
    "    \"\"\"Return a vector x with np.matmul(A, x) == b using\n",
    "    the Gauss Elimination algorithm, without partial pivoting.\"\"\"\n",
    "    for k in range(A.shape[0] - 1):\n",
    "        eliminate(A, b, k)\n",
    "    x = backward_substitution(A, b)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (a)  -- 3 pt\n",
    "\n",
    "Write a function `forward_substitution` that solves the lower-triangular\n",
    "system $Ax = b$.\n",
    "\n",
    "Hint: This function should be very similar to the `backward_substitution` function\n",
    "from the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_substitution(A, b):\n",
    "    \"\"\"Return a vector x with np.matmul(A, x) == b, where \n",
    "        * A is an nxn numpy matrix that is lower-triangular and non-singular\n",
    "        * b is a nx1 numpy vector\n",
    "        \n",
    "    >>> A = np.array([[2., 0.],\n",
    "                      [1., -2.]])\n",
    "    >>> b = np.array([1., 2.])\n",
    "    >>> forward_substitution(A, b)\n",
    "    array([ 0.5 , -0.75])\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (b) -- 4 pt \n",
    "\n",
    "Write a function `elementary_elimination_matrix` that returns the \n",
    "$k$-th elementary elimination matrix ($M_k$ in your notes).\n",
    "\n",
    "You may assume that `A[i,j] = 0 for i > j, j < k - 1`. (The subtraction\n",
    "in $k-1$ is because in Python, indices begin at 0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def elementary_elimination_matrix(A, k):\n",
    "    \"\"\"Return the elements below the k-th diagonal of the\n",
    "    k-th elementary elimination matrix.    \n",
    "\n",
    "    (Do not use partial pivoting, since we haven't\n",
    "    introduced the idea yet.)\n",
    "    \n",
    "    Precondition: A is an nxn numpy matrix, non-singular\n",
    "                  A[i,j] = 0 for i > j, j < k - 1\n",
    "\n",
    "    As always, these examples are for your understanding only.\n",
    "    The actual Python output might differ slightly.\n",
    "    \n",
    "    >>> A = np.array([[2., 0., 1.],\n",
    "                      [1., 1., 0.],\n",
    "                      [2., 1., 2.]])\n",
    "    >>> elementary_elimination_matrix(A, 1)\n",
    "    np.array([[1., 0., 0.],\n",
    "              [-0.5, 1., 0.],\n",
    "              [-1., 0., 1.]])\n",
    "    >>> A = np.array([[2., 0., 1.],\n",
    "                      [0., 1., 0.],\n",
    "                      [0., 1., 2.]])\n",
    "    >>> elementary_elimination_matrix(A, 2)\n",
    "    np.array([[1., 0., 0.],\n",
    "              [0., 1., 0.],\n",
    "              [0., -1., 1.]])\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (c) -- 4pt\n",
    "\n",
    "Write a function `lu_factorize` that factors a matrix A into its upper and lower\n",
    "triangular components. Use the function `elementary_elimination_matrix` as a helper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lu_factorize(A):\n",
    "    \"\"\"Return two matrices L and U, where\n",
    "            * L is lower triangular\n",
    "            * U is upper triangular\n",
    "            * and np.matmul(L, U) == A\n",
    "    >>> A = np.array([[2., 0., 1.],\n",
    "                      [1., 1., 0.],\n",
    "                      [2., 1., 2.]])\n",
    "    >>> L, U = lu_factorize(A)\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",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (d) -- 2pt\n",
    "\n",
    "Write a function `solve_lu` that solves a linear system Ax=b by\n",
    "    * factoring $A = LU$ (using the `lu_factorize` function)\n",
    "    * solving $Ly = b$ using forward substitution (using the `forward_substitution` function)\n",
    "    * solving $Ux = y$ using backward substitution (using the `backward_substitution` function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_lu(A, b):\n",
    "    \"\"\"Return a vector x with np.matmul(A, x) == b using\n",
    "    LU factorization. (Do not use partial pivoting, since we haven't\n",
    "    introduced the idea yet.)\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part (e) -- 5 pt\n",
    "\n",
    "Write a function `invert_matrix` that takes an $n \\times n$ matrix `A`, and computes its inverse by solving $n$\n",
    "systems of linear equations of the form $Ax = b$.\n",
    "\n",
    "You can use the functions that you wrote earlier, but note that you **will** be graded on efficiency.\n",
    "In particular, avoid writing code that repeats a computation unnecessarily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def invert_matrix(A):\n",
    "    \"\"\"Return the inverse of the nxn matrix A by \n",
    "    solving n systems of linear equations of the form \n",
    "    Ax = b.\n",
    "    \n",
    "    >>> A = np.array([[0.5, 0.],\n",
    "                      [-1., 2.]])\n",
    "    >>> invert_matrix(A)\n",
    "    array([[ 2. , -0. ],\n",
    "           [ 1. ,  0.5]])\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3\n",
    "\n",
    "### Part (a) -- 4 pt\n",
    "\n",
    "Prove that the matrix\n",
    "\n",
    "\\begin{align*}\n",
    "A = \\begin{bmatrix}\n",
    "0 & 1 \\\\\n",
    "1 & 0 \\\\\n",
    "\\end{bmatrix}\n",
    "\\end{align*}\n",
    "     \n",
    "has no $LU$ facotirzation.\n",
    "That is, there are no lower triangular matrix $L$ and upper\n",
    "triangular matrix $U$ such that $A = LU$.\n",
    "\n",
    "Include your proof in your pdf writeup.\n",
    "\n",
    "### Part (b) -- 2 pt\n",
    "\n",
    "Show that for an elementary matrix $M_k = I - {\\bf m}{\\bf e_k^T}$, we have $M_k^{-1} = I + {\\bf m}{\\bf e_k^T}$.\n",
    "(Property 4 of the elementary elimination matrix from your lecture notes)\n",
    "\n",
    "Include your solution in your pdf writeup."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}