{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import sys, os\n",
    "import numpy as np\n",
    "import scipy\n",
    "import skimage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# MAX_ITERS = 10\n",
    "TOL = 1.0e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## LU Factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Gaussian Elimination"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This section is based on lectures 20-22 in Trefethen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "If you are unfamiliar with Gaussian elimination or need a refresher, watch [this Khan Academy video](https://www.khanacademy.org/math/precalculus/precalc-matrices/row-echelon-and-gaussian-elimination/v/matrices-reduced-row-echelon-form-2)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's use Gaussian Elimination by hand to review:\n",
    "\n",
    "A =\n",
    " \\begin{pmatrix}\n",
    "  1 & -2 & -2 & -3 \\\\\n",
    "  3 & -9 & 0 & -9 \\\\\n",
    "  -1 & 2 & 4 & 7  \\\\\n",
    "  -3 & -6 & 26 & 2\n",
    " \\end{pmatrix}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Answer:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$$ LU = \\begin{bmatrix} 1 & 0 & 0 & 0\\\\ 3 & 1 & 0 & 0 \\\\ -1 & 0 & 1 & 0 \\\\ -3 & 4 & -2 & 1\\end{bmatrix} \\cdot \\begin{bmatrix} 1 & -2 & -2 & -3 \\\\ 0 & -3 & 6 & 0 \\\\ 0 & 0 & 2 & 4 \\\\ 0 & 0 & 0 & 1 \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Above example is from Lectures 20, 21 of Trefethen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**Gaussian Elimination** transform a linear system into an upper triangular one by applying linear transformations on the left.  It is *triangular triangularization*.\n",
    "\n",
    "$ L_{m-1} \\dots L_2 L_1 A = U $\n",
    "\n",
    "L is *unit lower-triangular*: all diagonal entries are 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def LU(A):\n",
    "    U = np.copy(A)\n",
    "    m, n = A.shape\n",
    "    L = np.eye(n)\n",
    "    for k in range(n-1):\n",
    "        for j in range(k+1,n):\n",
    "            L[j,k] = U[j,k]/U[k,k]\n",
    "            U[j,k:n] -= L[j,k] * U[k,k:n]\n",
    "    return L, U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "A = np.array([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]).astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "L, U = LU(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.allclose(A, L @ U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The LU factorization is useful!\n",
    "\n",
    "Solving Ax = b becomes LUx = b:\n",
    "1. find A = LU\n",
    "2. solve Ly = b\n",
    "3. solve Ux = y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**Work**\n",
    "\n",
    "Work for Gaussian Elimination: $2\\cdot\\frac{1}{3} n^3$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**Memory**\n",
    "\n",
    "Above, we created two new matrices, $L$ and $U$.  However, we can store the values of $L$ and $U$ in our matrix A (overwriting the original matrix).  Since the diagonal of $L$ is all $1$s, it doesn't need to be stored.  Doing factorizations or computations **in-place** is a common technique in numerical linear algebra to save memory.  Note: you wouldn't want to do this if you needed to use your original matrix $A$ again in the future.  One of the homework questions is to rewrite the LU method to operate in place."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Consider the matrix $$ A = \\begin{bmatrix} 10^{-20} & 1 \\\\ 1 & 1 \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "A = np.array([[1e-20, 1], [1,1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "By hand, use Gaussian Elimination to calculate what L and U are:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "#Exercise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "#Exercise:\n",
    "L1, U1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "L2, U2 = LU(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "L2, U2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.allclose(L1, L2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.allclose(U1, U2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "np.allclose(A, L2 @ U2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This is the motivation for $LU$ factorization **with pivoting**.\n",
    "\n",
    "This also illustrates that LU factorization is *stable*, but not *backward stable*. (spoiler alert: even with partial pivoting, LU is \"explosively unstable\" for certain matrices, yet stable in practice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Stability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "An algorithm $\\hat{f}$ for a problem $f$ is **stable** if for each $x$,\n",
    "$$ \\frac{\\lVert \\hat{f}(x) - f(y) \\rVert}{ \\lVert f(y) \\rVert } = \\mathcal{O}(\\varepsilon_{machine}) $$\n",
    "\n",
    "for some $y$ with\n",
    "$$ \\frac{\\lVert y - x \\rVert }{\\lVert x \\rVert} = \\mathcal{O}(\\varepsilon_{machine}) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**A stable algorithm gives nearly the right answer to nearly the right question** (Trefethen, pg 104)\n",
    "\n",
    "To translate that:\n",
    "- right question: $x$\n",
    "- nearly the right question: $y$\n",
    "- right answer: $f$\n",
    "- right answer to nearly the right question: $f(y)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Backwards Stability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Backwards stability is both **stronger** and **simpler** than stability. \n",
    "\n",
    "An algorithm $\\hat{f}$ for a problem $f$ is **backwards stable** if for each $x$,\n",
    "$$ \\hat{f}(x) = f(y) $$\n",
    "\n",
    "for some $y$ with\n",
    "$$ \\frac{\\lVert y - x \\rVert }{\\lVert x \\rVert} = \\mathcal{O}(\\varepsilon_{machine}) $$\n",
    "\n",
    "**A backwards stable algorithm gives exactly the right answer to nearly the right question** (Trefethen, pg 104)\n",
    "\n",
    "Translation:\n",
    "\n",
    "- right question: $x$\n",
    "- nearly the right question: $y$\n",
    "- right answer: $f$\n",
    "- right answer to nearly the right question: $f(y)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### LU factorization with Partial Pivoting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's now look at the matrix $$ \\hat{A} = \\begin{bmatrix} 1 & 1 \\\\ 10^{-20} & 1 \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "A = np.array([[1,1], [1e-20, 1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "By hand, use Gaussian Elimination to calculate what L and U are:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "#Exercise:\n",
    "L1, U1 = np.array([[1, 0],[10**-20, 1]]), np.array([[1, 1],[0, 1 - 10**-20]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "L2, U2 = LU(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.allclose(A, L2 @ U2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Idea: We can switch the order of the rows around to get more stable answers! This is equivalent to multiplying by a permutation matrix $P$.  For instance,\n",
    "\n",
    "$$\\begin{bmatrix} 0 & 1 \\\\ 1 & 0 \\end{bmatrix} \\cdot \\begin{bmatrix} 10^{-20} & 1 \\\\ 1 & 1 \\end{bmatrix} =  \\begin{bmatrix} 1 & 1 \\\\ 10^{-20} & 1 \\end{bmatrix} $$\n",
    "\n",
    "$$ PA = \\hat{A} $$\n",
    "\n",
    "Apply Gaussian elimination to $PA$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "At each step, choose the largest value in column k, and move that row to be row k."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Homework"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def swap(a,b):\n",
    "    temp = np.copy(a)\n",
    "    a[:] = b\n",
    "    b[:] = temp\n",
    "\n",
    "a=np.array([1,2,3])\n",
    "b=np.array([3,2,1])\n",
    "swap(a,b)\n",
    "a,b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "#Exercise: re-write the LU factorization above to use pivoting\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "A = np.array([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]).astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "L, U, P = LU_pivot(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Can compare below to answers in Trefethen, page 159:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**Partial pivoting** permutes the rows. It is such a universal practice, that this is usually what is meant by *LU factorization*.\n",
    "\n",
    "**Complete pivoting** permutes the rows and columns.  Complete pivoting is significantly time-consuming and rarely used in practice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Consider the system of equations:\n",
    "\n",
    "$$ \\begin{bmatrix} 1 & 0 & 0  & 0 & 0 & 1 \\\\ -1 & 1 & 0  & 0 & 0 & 1 \\\\ -1 & -1 & 1  & 0 & 0 & 1 \\\\ -1 & -1 & -1  & 1 & 0 & 1  \\\\  -1 & -1 & -1  & -1 & 1 & 1 \\\\ -1 & -1 & -1  & -1 & -1 & 1 \\end{bmatrix} \\mathbf{x} = \\begin{bmatrix} 1 \\\\ 1 \\\\ 1  \\\\ 1 \\\\ 2 \\\\ 1 \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def make_matrix(n):\n",
    "    A = np.eye(n)\n",
    "    for i in range(n):\n",
    "        A[i,-1] = 1\n",
    "        for j in range(i):\n",
    "            A[i,j] = -1\n",
    "    return A "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def make_vector(n):\n",
    "    b = np.ones(n)\n",
    "    b[-2] = 2\n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "make_vector(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Exercise: Let's use Gaussian Elimination on the $5 \\times 5$ system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Scipy has this funtionality as well.  Let's look at the solution for the last 5 equations with $n=10,\\,20,\\,30,\\,40,\\,50,\\,60$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "?scipy.linalg.solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "for n, ls in zip(range(10, 70, 10), ['--', ':', '-', '-.', '--', ':']):\n",
    "    soln = scipy.linalg.lu_solve(scipy.linalg.lu_factor(make_matrix(n)), make_vector(n))\n",
    "    plt.plot(soln[-5:], ls)\n",
    "    print(soln)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "What is going on when $n=60$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**Theorem**: Let the factorization $PA = LU$ of a matrix A be computed by Gaussian Elimination with partial pivoting.  The *computed* (by a computer with Floating Point Arithmetic) matrices $\\hat{P}$, $\\hat{L}$, and $\\hat{U}$ satisfy\n",
    "\n",
    "$$\\hat{L}\\hat{U} = \\hat{P} A + \\delta A, \\quad \\frac{\\delta A}{A} = \\mathcal{O}(\\rho \\varepsilon_{machine}) $$\n",
    "\n",
    "where $\\rho$ is the *growth factor*, \n",
    "\n",
    "$$\\rho = \\frac{max_{i,j} \\lvert u_{ij} \\rvert }{max_{i,j} \\lvert a_{ij} \\rvert } $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "For our matrix above, $\\rho = 2^{m-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Unstable in theory, stable in practice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Stability of most algorithms (such as QR) is straightforward.  Not the case for Gaussian Elimination with partial pivoting.  Instability in Gaussian elimination (with or without pivoting) arises only if L and/or U is large relative to the size of A.\n",
    "\n",
    "Trefethen: \"Despite examples like (22.4), Gaussian elimination with partial pivoting is utterly stable in practice... In fifty years of computing, no matrix problems that excite an explosive instability are known to have arisen under natural circumstances.\" [although can easily be constructed as contrived examples]\n",
    "\n",
    "Although some matrices cause instability, but extraordinarily small proportion of all matrices so \"never\" arise in practice for statistical reasons.  \"If you pick a billion matrices at random, you will almost certainly not find one for which Gaussian elimination is unstable.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Further Reading\n",
    "- Gaussian Elimination/LU factorization-- Trefethn Lecture 20\n",
    "- Pivoting -- Trefethn Lecture 21\n",
    "- Stability of Gaussian Elimination -- Trefethn Lecture 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}