05 Linear Equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\"> Computation for Physicists </h1>\n",
    "<h2 align=\"center\"> <em> Solving Linear Equations</em> </h2>\n",
    "<h2 align=\"center\" > <a href=\"mailto:duan@unm.edu\">Dr. Duan</a> (UNM) </h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Error Analysis\n",
    "- Total error = <span style='color:blue'> computational error </span> + <span style='color:red'> propagated data error </span>\n",
    "$$\\hat{f}(\\hat{x}) - f(x) = \\color{blue}{[\\hat{f}(\\hat{x}) - f(\\hat{x})]} + \\color{red}{[f(\\hat{x}) - f(x)]},$$\n",
    "where $f$ is the exact algorithm approximated by $\\hat{f}$.\n",
    "- Computational error = truncation error + rounding error.\n",
    "    - Rounding error: due to the limited precision.\n",
    "    - Truncation error: due to the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def mycos(x, tol=1e-6):\n",
    "    # x = x % (2.0*pi) # move x to [0, 2*pi)\n",
    "    res = 0. # result\n",
    "    term = 1. # 1st term in the Taylor series\n",
    "    n = 0 # power of x for the term\n",
    "    while abs(term) >= tol: \n",
    "        res += term # add the term to the result\n",
    "        term *= -x*x / ((n+1)*(n+2)) # compute the new term\n",
    "        n += 2 # skip the odd power of x    \n",
    "    # res += term # add the last term\n",
    "    print(f\"(-1)^{n//2} x^{n} / {n}! = {term}\")\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Is `tol` the tolerance of the absolute error or the relative error?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import math\n",
    "math.cos(3.) - mycos(3.) # dominated by truncation error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "math.cos(100.) - mycos(100.) # dominated by rounding error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- The _condition number_ \n",
    "$$\\mathrm{cond} =  \\left|\\frac{[f(\\hat{x}) - f(x)]/f(x)}{(\\hat{x}-x)/x}\\right| = \\left|\\frac{x f'(x)}{f(x)}\\right|$$\n",
    "measures the amplification of the data error in the original problem $f$.\n",
    "- A problem is _ill-conditioned_ if $\\mathrm{cond}\\gg1$, which always amplifies the error no matter which algorithm you choose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "math.tan(1.57078), math.tan(1.57079) # infinite cond at pi/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- A reliable numerical solution is obtained for a <span style='color:red'>well-conditioned problem</span> (with a small propagated data error) using a <span style='color:blue'>stable algorithm</span> (with a small computational error)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Homework 4\n",
    "- The standard deviation of a dataset $\\{x_1, x_2, \\cdots, x_n\\}$ is defined as\n",
    "$$\\sigma = \\sqrt{\\frac{1}{n}\\sum_{i=1}^n (x_i - \\bar{x})^2},$$\n",
    "where $\\bar{x}$ is the average value of $x_i$. It is straightforward to show that an equivalent definition of the standard deviation is\n",
    "$$\\sigma = \\sqrt{\\left(\\frac{1}{n}\\sum_{i=1}^n x_i^2\\right)- \\bar{x}^2}.$$\n",
    "Which of the two algorithms is potentially faster to compute, and which is more accurate?\n",
    "\n",
    "- Verify your conjecture with some example(s). You may find the NumPy functions `numpy.random.random()`, `numpy.sum()`, `numpy.mean()`, and `numpy.std()` useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import random\n",
    "\n",
    "n = 1000 # number of values\n",
    "data = np.empty((n,), dtype=np.float32)\n",
    "data[:] = 10+random(n) # dataset\n",
    "sigma = np.std(data)\n",
    "mean = np.mean(data)\n",
    "print(f\"mean = {mean}, sigma = {sigma}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "std1=np.sqrt(np.sum((data-mean)**2)/n)\n",
    "std2=np.sqrt(np.sum(data**2)/n - mean**2)\n",
    "print(f\"std1 has relative error {(std1-sigma)/sigma}\")\n",
    "print(f\"std2 has relative error {(std2-sigma)/sigma}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%timeit np.sqrt(np.sum((data-mean)**2)/n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%timeit np.sqrt(np.sum(data**2)/n - mean**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%timeit np.std(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import homework.hw4 as hw4\n",
    "\n",
    "help(hw4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "(hw4.std1(data)-sigma)/sigma, (hw4.std2(data)-sigma)/sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%timeit hw4.std1(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%timeit hw4.std2(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Systems of Linear Equations\n",
    "- Most of the numerical problems are reduced to solving a system of linear equations,\n",
    "$$\\mathbf{A} \\boldsymbol{x} \n",
    "=\\begin{bmatrix}\n",
    "a_{00} & a_{01} & \\cdots & a_{0,n_1} \\\\\n",
    "a_{10} & a_{11} & \\cdots & a_{1,n_1} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "a_{m-1,0} & a_{m-1,1} & \\cdots & a_{m-1,n-1}\n",
    "\\end{bmatrix}\\begin{bmatrix}\n",
    "x_0 \\\\ x_1 \\\\ \\vdots \\\\ x_{n-1}\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "b_0 \\\\ b_1 \\\\ \\vdots \\\\ b_{m-1}\n",
    "\\end{bmatrix}\n",
    "= \\boldsymbol{b}.$$\n",
    "- For $\\mathbf{A}\\in\\mathbb{R}^{m\\times n}$, the solution exists if $\\boldsymbol{b}\\in\\mathrm{span}(\\mathbf{A})=\\{\\mathbf{A}\\boldsymbol{x}: \\boldsymbol{x}\\in\\mathbb{R}^n\\}$.\n",
    "- If $m=n$, the solution exists and is unique if $\\mathbf{A}$ is _nonsingular_, which implies that $\\mathbf{A}^{-1}$ exists, $\\det(\\mathbf{A})\\neq0$, $\\mathbf{A}$ has rank $n$, and $\\mathbf{A}\\boldsymbol{z}\\neq\\mathbf{0}$ for any nonzero vector $\\boldsymbol{z}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- The [`linalg`](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html) module in SciPy (and NumPy) provides the linear algebra tools for solving linear equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.linalg as la\n",
    "n = 100 # problem size\n",
    "A = random(n*n).reshape((n, n)) # random matrix\n",
    "x = random(n) # expected solution\n",
    "b = A @ x # same as np.dot(A, x)\n",
    "\n",
    "x1 = la.inv(A) @ b # x = A^-1 b\n",
    "x2 = la.solve(A, b) # solve x from linear eqs\n",
    "print('||x1-x|| = ', la.norm(x1-x)) # absolute error\n",
    "print('||x2-x|| = ', la.norm(x2-x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Solving linear equations directly is faster than inverting the matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%timeit la.inv(A) @ b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit la.solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Vector Norms\n",
    "- The $p$-norm of a vector $\\boldsymbol{x}$ is defined as \n",
    "$$\\Vert\\boldsymbol{x}\\Vert_p = \\left(\\sum_{i=0}^{n-1} |x_i|^p \\right)^{1/p}.$$\n",
    "For example, $\\Vert\\boldsymbol{x}\\Vert_1 = \\sum_i |x_i|$, $\\Vert\\boldsymbol{x}\\Vert_2 = \\sqrt{\\sum_i |x_i|^2}$, and $\\Vert\\boldsymbol{x}\\Vert_\\infty = \\max |x_i|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Unit Circles](fig/vecnorm.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([1, -2, 3])\n",
    "print(f'Vector {x} has 2-norm {la.norm(x)}.')\n",
    "print(f'It also has 1-norm {la.norm(x, 1)} and ∞-norm {la.norm(x, np.inf)}.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- All norms are equivalent up to a constant factor.\n",
    "- For any vector norm:\n",
    "    - $\\Vert\\boldsymbol{x}\\Vert >0 $ if $\\boldsymbol{x}\\neq\\mathbf{0}$.\n",
    "    - $\\Vert \\gamma \\boldsymbol{x}\\Vert = |\\gamma|\\cdot \\Vert\\boldsymbol{x}\\Vert $ for any scalar $\\gamma$.\n",
    "    - $\\Vert\\boldsymbol{x} + \\boldsymbol{y}\\Vert \\leq \\Vert\\boldsymbol{x}\\Vert + \\Vert\\boldsymbol{y}\\Vert $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Matrix Norms\n",
    "- A vector-induced matrix norm is defined as the largest stretching it can do to a nonzero vector:\n",
    "$$\\Vert\\mathbf{A}\\Vert = \\max_{\\boldsymbol{x}\\neq\\mathbf{0}} \n",
    "\\frac{\\Vert \\mathbf{A}\\boldsymbol{x}\\Vert}{\\Vert\\boldsymbol{x}\\Vert}.$$\n",
    "- It can be shown that\n",
    "$$\\Vert\\mathbf{A}\\Vert_1 = \\max_j \\sum_i |a_{ij}|\n",
    "\\qquad\\mathrm{and}\\qquad\n",
    "\\Vert\\mathbf{A}\\Vert_\\infty = \\max_i \\sum_j |a_{i j}|.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[1, 2], [3, 4]])\n",
    "print(f\"{A} has a Frobenius norm {la.norm(A)}.\") # default is not 2-norm\n",
    "print(f\"It also has 1-norm {la.norm(A, 1)} and ∞-norm {la.norm(A, np.inf)}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- A vector-induced matrix norm has the following properties:\n",
    "    - $\\Vert\\mathbf{A}\\Vert > 0$ if $\\mathbf{A}\\neq\\mathbf{0}$.\n",
    "    - $\\Vert\\gamma\\mathbf{A}\\Vert = |\\gamma|\\cdot\\Vert\\mathbf{A}\\Vert$ for any scalar $\\gamma$.\n",
    "    - $\\Vert\\mathbf{A} + \\mathbf{B}\\Vert\\leq \\Vert\\mathbf{A}\\Vert + \\Vert\\mathbf{B}\\Vert$.\n",
    "    - $\\Vert\\mathbf{A}\\mathbf{B}\\Vert \\leq \\Vert\\mathbf{A}\\Vert \\cdot \\Vert\\mathbf{B}\\Vert$.\n",
    "    - $\\Vert\\mathbf{A} \\boldsymbol{x}\\Vert \\leq \\Vert\\mathbf{A}\\Vert \\cdot \\Vert\\boldsymbol{x}\\Vert$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Matrix Condition Number\n",
    "- The condition number of a square nonsingular matrix $\\mathbf{A}$ is defined by\n",
    "$$\\mathrm{cond}(\\mathbf{A}) = \\Vert\\mathbf{A}\\Vert \\cdot \\Vert\\mathbf{A}^{-1}\\Vert\n",
    "= \\left(\\max_{\\boldsymbol{x}\\neq\\mathbf{0}} \\frac{\\Vert\\mathbf{A}\\boldsymbol{x}\\Vert}{\\Vert\\boldsymbol{x}\\Vert}\\right)\n",
    "\\left(\\min_{\\boldsymbol{x}\\neq\\mathbf{0}} \\frac{\\Vert\\mathbf{A}\\boldsymbol{x}\\Vert}{\\Vert\\boldsymbol{x}\\Vert}\\right)^{-1}.$$\n",
    "It is the ratio of the maximum stretching to the maximum shrinking a matrix does to any nonzero vector.\n",
    "- $\\mathrm{cond}(\\mathbf{A})=\\mathrm{cond}(\\mathbf{A}^{-1})\\geq 1$.\n",
    "- $\\mathrm{cond}(\\mathbf{A})=\\infty$ if $\\mathbf{A}$ is singular.\n",
    "- $\\mathrm{cond}(\\gamma \\mathbf{A}) = \\mathrm{cond}(\\mathbf{A})$.\n",
    "- $\\mathrm{cond}(\\mathbf{D})=(\\max |d_{ii}|)/(\\min |d_{ii}|)$ if $\\mathbf{D}$ is diagonal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Suppose $\\mathbf{A}\\boldsymbol{x}=\\boldsymbol{b}$ and $\\mathbf{A}(\\boldsymbol{x}+\\Delta\\boldsymbol{x})\n",
    "=\\boldsymbol{b} + \\Delta\\boldsymbol{b}$, then\n",
    "$$\\frac{\\Vert\\Delta\\boldsymbol{x}\\Vert}{\\Vert\\boldsymbol{x}\\Vert}\n",
    "\\leq \\mathrm{cond}(\\mathbf{A})\n",
    "\\frac{\\Vert\\Delta\\boldsymbol{b}\\Vert}{\\Vert\\boldsymbol{b}\\Vert}.$$\n",
    "- An _ill-conditioned_ problem (with a very large condition number) can produce a poor solution because of the error amplification."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Ill-conditioning can be because of the poor scaling.\n",
    "- This can be because of the poor choice of units. Choose the right units so that the magnitudes of the matrix elements are comparable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from numpy.linalg import cond\n",
    "\n",
    "A = np.array([[1, 1e-14], [2, -1e-14]])\n",
    "x = np.array([1, 2e14])\n",
    "b = A @ x\n",
    "\n",
    "x1 = la.solve(A, b) # solver equilibrates A\n",
    "print(f\"cond(A) = {cond(A)}, x1 = {x1}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Ill-conditioning may also come from a near singularity condition, e.g. a column or row vector is approximately a linear combination of other vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "A = np.array([[1-1e-12, 1], [1, 1]])\n",
    "x = np.array([1, 2])\n",
    "b = A @ x\n",
    "\n",
    "x1 = la.solve(A, b) \n",
    "cond(A), la.norm(x1-x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Gassian Elimination\n",
    "- The basic technique of solving $\\mathbf{A}\\boldsymbol{x}=\\boldsymbol{b}$ is called _Gaussian elimination_. In the first step, one uses a series of matrices $\\mathbf{M}^{(k)}$ to make $\\mathbf{A}$ upper triangular:\n",
    "$$\\mathbf{M}^{(n-2)}\\cdots\\mathbf{M}^{(0)} \\mathbf{A}\n",
    "=\\mathbf{U}\n",
    "=\\begin{bmatrix} u_{00} & u_{01} & \\cdots & u_{0,n-1} \\\\\n",
    "0 & u_{11} & \\cdots & u_{1,n-1} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "0 & 0 & \\cdots & u_{n-1,n-1}\\end{bmatrix}.$$\n",
    "\n",
    "- In the second step, one solves $\\mathbf{U}\\boldsymbol{x}=\\boldsymbol{y} =\\mathbf{M}^{(n-2)}\\cdots\\mathbf{M}^{(0)}\\boldsymbol{b}$ with _backward substitution_:\n",
    "$$x_i = \\begin{cases}\n",
    "y_{n-1}/u_{n-1,n-1} & \\quad\\mathrm{if}\\, i=n-1, \\\\\n",
    "\\left(y_i - \\sum_{j=i+1}^{n-1}u_{ij}x_j\\right)u_{ii}^{-1} & \\quad\\mathrm{otherwise}.\n",
    "\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Let $\\mathbf{A}^{(k)}=\\mathbf{M}^{(k-1)}\\cdots\\mathbf{M}^{(0)} \\mathbf{A}$. Then\n",
    "$\\mathbf{M}^{(k)} = \\mathbf{1}-\\boldsymbol{m}^{(k)}\\boldsymbol{e}_k^T$ annhilates the elements of the $k$th column vector of $\\mathbf{A}^{(k)}$ below $a^{(k)}_{kk}$, where $a^{(k)}_{kk}$ is called the _pivot_, and\n",
    "$$\\boldsymbol{m}^{(k)} \n",
    "=\\begin{bmatrix}\n",
    "m^{(k)}_0 \\\\ \\vdots \\\\ m^{(k)}_k \\\\ m^{(k)}_{k+1} \\\\ \\vdots \\\\ m^{(k)}_{n-1}\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "0 \\\\ \\vdots \\\\ 0  \\\\ a^{(k)}_{k+1,k}/a^{(k)}_{k,k} \\\\ \\vdots \\\\ a^{(k)}_{n-1,k}/a^{(k)}_{k,k}\n",
    "\\end{bmatrix}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- It can be shown that\n",
    "\\begin{align*}\n",
    "\\mathbf{L} &= (\\mathbf{M}^{(k-1)}\\cdots\\mathbf{M}^{(0)})^{-1}\n",
    "= \\mathbf{1} + \\sum_{k=0}^{n-2}\\boldsymbol{m}^{(k)}\\boldsymbol{e}_k^T\\\\\n",
    "& = \\begin{bmatrix}\n",
    "1 & 0 & \\cdots & 0 \\\\\n",
    "l_{10} & 1 & \\cdots &0 \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "l_{n-1,0} & l_{n-1,1} & \\cdots & 1\n",
    "\\end{bmatrix}\n",
    "\\end{align*}\n",
    "is lower triangular. Therefore, \n",
    "$$\\mathbf{A}=\\mathbf{L}\\mathbf{U}.$$ \n",
    "This is known as the _LU factorization_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Explicit LU factorization is useful when many linear systems with the same $\\mathbf{A}$ but different $\\boldsymbol{b}$'s are to be solved.\n",
    "- Given the LU factorization of $\\mathbf{A}$, the first step of Gassian elimination can be achieved by solving $\\mathbf{L}\\boldsymbol{y}=\\boldsymbol{b}$ with _forward substitution_:\n",
    "$$y_i = \\begin{cases}\n",
    "b_0 & \\quad\\mathrm{if}\\, i=0, \\\\\n",
    "b_i - \\sum_{j=0}^{i-1}l_{ij}y_j & \\quad\\mathrm{otherwise}.\n",
    "\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Gassian elimination breaks down if the pivot $a^{(k)}_{kk}$ is 0. This can be fixed by permuting the rows of $\\mathbf{A}^{(k)}$, which is known as (partial) pivoting.\n",
    "\n",
    "- Gaussian elimination becomes more stable when the largest pivot is chosen at each step.\n",
    "\n",
    "- $\\mathbf{A}$ is singular if there is no nonzero pivot, but the LU factorization can still be completed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])\n",
    "lu, piv = la.lu_factor(A) # LU factorization with pivoting\n",
    "b1 = np.array([1, 1, 1, 1])\n",
    "b2 = np.array([1, 2, 3, 4])\n",
    "x1 = la.lu_solve((lu, piv), b1) # solve with LU and pivoting\n",
    "x2 = la.lu_solve((lu, piv), b2)\n",
    "la.norm(A @ x1 - b1), la.norm(A @ x2 - b2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Homework Problem\n",
    "- Think how the LU factorization can be used to compute the determinent and inverse of a nonsingular matrix.\n",
    "- Implement a simple function to compute the inverse of an arbitrary square matrix using  `lu_factor()` and `lu_solve()` in `scipy.linalg`. Document and validate the function.\n",
    "- Hint: For $\\mathbf{A}\\mathbf{X}=\\mathbf{1}$, you can treat $\\mathbf{X}$ and $\\mathbf{1}$ as collections of column vectors. `numpy.eye(n)` returns a $n\\times n$ identity matrix."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}