Backward error of the LU decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "## Backward error of the LU decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## An almost backward error result\n",
    "\n",
    "**Theorem**\n",
    "Let the factorization $A=LU$ of a nonsingular matrix $A\\in\\mathbb{R}^{m\\times m}$ be computed by Gaussian elimination without pivoting. If $A$ has an LU factorization, then for all sufficiently small $\\epsilon_{mach}$, the factorization completes successfully in floating point arithmetic, and the computed matrices $\\tilde{L}$ and $\\tilde{U}$ satisfy\n",
    "\n",
    "$$\\tilde{L}\\tilde{U} = A+\\delta A,\\quad \\frac{\\|\\delta A\\|}{\\|L\\|\\|U\\|} = \\mathcal{O}(\\epsilon_{mach})$$\n",
    "\n",
    "for some $\\delta A\\in\\mathbb{R}^{m\\times m}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Does the above theorem tell us that LU decomposition is backward stable? No, it doesn't!! The crucial fact is that the error is measured relative to $\\|L\\|\\|U\\|$. But the backward error is a measure relative to $\\|A\\|$. Hence, backward stability follows only if $\\|L\\|\\|U\\|\\approx \\|A\\|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Consider again the example with\n",
    "\n",
    "$$\n",
    "A_{\\epsilon} = \\begin{bmatrix}\\epsilon & 1\\\\ 1 & 1\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "We have\n",
    "$$\n",
    "A_{\\epsilon} = \\begin{bmatrix}1 & 0 \\\\ \\epsilon^{-1} & 1\\end{bmatrix}\\begin{bmatrix}\\epsilon & 1 \\\\ 0 & 1-\\epsilon^{-1}\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Hence, $\\|L\\|\\|U\\|\\approx \\epsilon^{-2}$ for very small $\\epsilon$. But $\\|A_\\epsilon\\|=\\mathcal{O}(1)$, and we cannot expect backward stability. Indeed, we showed that for small $\\epsilon$ the LU is not backward stable in this example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Growth Factors\n",
    "\n",
    "**Definition (Growth factor)**\n",
    "The growth factor $\\rho$ is defined as the ratio\n",
    "$$\n",
    "\\rho:=\\frac{\\max_{ij}|u_{ij}|}{\\max_{ij}|a_{ij}|}.\n",
    "$$\n",
    "\n",
    "In Gaussian elimination with partial pivoting we have $\\|L\\| = \\mathcal{O}(1)$ since $|\\ell_{ij}|\\leq 1$ for all multipliers $\\ell_{ij}$. Hence, the main contribution to the product $\\|L\\|\\|U\\|$ is in $\\|U\\|$, which is measured by the growth factor. Let us go back to the previous result. We have that $\\|L\\|\\approx 1$ and $\\|U\\|\\approx \\rho\\|A\\|$. Hence, the backward error result\n",
    "\n",
    "$$\n",
    "\\frac{\\|\\delta A\\|}{\\|L\\|\\|U\\|} = \\mathcal{O}(\\epsilon_{mach})\n",
    "$$\n",
    "becomes\n",
    "$$\n",
    "\\frac{\\|\\delta A\\|}{\\|A\\|} = \\mathcal{O}(\\rho\\epsilon_{mach}).\n",
    "$$\n",
    "\n",
    "\n",
    "The following theorem makes this precise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "**Theorem (Backward stability of LU factorization)**\n",
    "Let the factorization $PA=LU$ of a matrix $A\\in\\mathbb{R}^{m\\times m}$ be computed by Gaussian elimination with partial pivoting. Then the computed matrices $\\tilde{P}$, $\\tilde{L}$, $\\tilde{U}$ satisfy\n",
    "\n",
    "$$\n",
    "\\tilde{L}\\tilde{U}=\\tilde{P}A+\\delta A,\\quad \\frac{\\|\\delta A\\|}{\\|A\\|}=\\mathcal{O}(\\rho\\epsilon_{mach}).\n",
    "$$\n",
    "If $\\ell_{ij}<1$ for each $i>j$ (there are no ties in the selection of pivots in exact arithmetic), then $\\tilde{P}=P$ for all sufficiently small $\\epsilon_{mach}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Worst-Case Instability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  0.  0.  0.  0.  0.  0.  1.]\n",
      " [-1.  1.  0.  0.  0.  0.  0.  1.]\n",
      " [-1. -1.  1.  0.  0.  0.  0.  1.]\n",
      " [-1. -1. -1.  1.  0.  0.  0.  1.]\n",
      " [-1. -1. -1. -1.  1.  0.  0.  1.]\n",
      " [-1. -1. -1. -1. -1.  1.  0.  1.]\n",
      " [-1. -1. -1. -1. -1. -1.  1.  1.]\n",
      " [-1. -1. -1. -1. -1. -1. -1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def matrix_with_large_rho(n):\n",
    "    A = np.zeros((n,n),dtype=np.float64)\n",
    "    for j in range(n-1):\n",
    "        A += np.diag(-np.ones(n-1-j,dtype=np.float64),-j-1)\n",
    "    A +=np.diag(np.ones(n,dtype=np.float64))\n",
    "    A[:,-1]=np.ones(n)\n",
    "    return A\n",
    "print(matrix_with_large_rho(8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Its LU decomposition has the following shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.linalg import lu\n",
    "P,L,U=lu(matrix_with_large_rho(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.   0.   0.   0.   0.   0.   0.   0.   0.   1.]\n",
      " [  0.   1.   0.   0.   0.   0.   0.   0.   0.   2.]\n",
      " [  0.   0.   1.   0.   0.   0.   0.   0.   0.   4.]\n",
      " [  0.   0.   0.   1.   0.   0.   0.   0.   0.   8.]\n",
      " [  0.   0.   0.   0.   1.   0.   0.   0.   0.  16.]\n",
      " [  0.   0.   0.   0.   0.   1.   0.   0.   0.  32.]\n",
      " [  0.   0.   0.   0.   0.   0.   1.   0.   0.  64.]\n",
      " [  0.   0.   0.   0.   0.   0.   0.   1.   0. 128.]\n",
      " [  0.   0.   0.   0.   0.   0.   0.   0.   1. 256.]\n",
      " [  0.   0.   0.   0.   0.   0.   0.   0.   0. 512.]]\n"
     ]
    }
   ],
   "source": [
    "print(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let us plot the growth factor with respect to n."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'rho')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmW0lEQVR4nO3deXhV1dn+8e9DMFj9iRNqLdRixSIhaMAUwQGDimBtlCpSUKpACIOmAo44a6EgFBFFZAikDC/gBMiUQEogjEdMAkoYZBCh5gUFh4KiAknW74+EvJEmkEBO9tnJ/bmuXOTsc87O7ebIk7XW3s825xwiIiIANbwOICIioUNFQUREiqgoiIhIERUFEREpoqIgIiJFanod4FTUqVPH1a9f3+sYIiK+kpWV9ZVz7oKSnvN1Uahfvz6ZmZlexxAR8RUz21Xac5o+EhGRIioKIiJSxJdFwcxizWz8/v37vY4iIlKl+LIoOOfmOed6nn322V5HERGpUnxZFEREJDhUFEREpIgvi4LWFESkOgsEAgwZMoRAIFDh+/bldQrOuXnAvOjo6Hivs4iIVKa0tDTatWuHc47w8HDS0tJo2bJlhe3flyMFEZHqKDU1lQ4dOpCbm0teXh6HDx8mPT29Qn+GioKISIj75ptv6NatG23btuXss8+mVq1ahIWFER4eTkxMTIX+LF9OH4mIVBczZ87koYce4quvvuKZZ57h2WefZd26daSnpxMTE1OhU0egoiAiEpL27NlDQkICs2bNomnTpixcuJCoqCgAWrZsWeHF4ChNH4mIhBDnHJMmTSIiIoIFCxbw8ssv8+GHHxYVhGDz5UjBzGKB2AYNGngdRUSkwuzcuZOePXvyr3/9i+uvv54JEybQsGHDSs3gy5GC2lyISFWSn5/PqFGjiIyMJBAIMHr0aJYtW1bpBQF8OlIQEakqNm/eTI8ePVi9ejXt2rVj3LhxXHLJJZ7l8eVIQUTE744cOcLf//53oqKi+OSTT5gyZQrJycmeFgTQSEFEpFIFAgGmTZvGokWL2L59Ox07duT111/noosu8joaoKIgIlJp0tPTadOmDbm5uQAMGTKEAQMGeJzq5zR9JCJSCVasWME999xTVBDCwsJwznmc6r+FTFEwsxgzW2FmY80sxus8IiIV4cCBAzz00EO0atWK0047LagtKipCUIuCmSWZ2V4z23DM9nZmtsXMtpvZ0bGTA74HTgdygplLRKQypKSkEBkZyZgxY+jXrx/btm1j6dKlDBw4sMK7m1YUC+bwxcxaUfAP/RTnXGThtjBgK9CGgn/8M4DOwCfOuXwzuwgY4Zy770T7j46OdpmZmUHLLyJyMr7++mv69+/P1KlTiYiIYOLEibRo0cLrWEXMLMs5F13Sc0EdKTjnlgPfHLO5ObDdObfDOXcYeAu40zmXX/j8t0CtYOYSEQkG5xzvvPMOjRo1YsaMGTz33HOsXbs2pArCiXhx9lFd4PNij3OAa8zsLqAtcA7wRmlvNrOeQE/A8/N5RUSO2r17Nw8++CBz5szh6quvZvHixVx55ZVexyo3L4qClbDNOedmAbNO9Gbn3Hgz2wPEhoeHX13h6UREysE5R1JSEo8++iiHDh1i2LBh9O/fn5o1/XnGvxdnH+UAvy72uB6wuzw7UO8jEQkFO3bs4JZbbqFHjx5ERUWxfv16Hn/8cd8WBPCmKGQAl5vZpWYWDnQC5pZnB2YWa2bj9+/fH5SAIiLHk5eXx8iRI2nSpAkZGRmMHTuWJUuWcPnll3sd7ZQF+5TUGUAAaGhmOWYW55zLBRKARcBm4B3n3Mby7FcjBRHxQiAQoF+/flx55ZX079+f1q1bs2nTJnr16kWNGiFz2dcpCeoYxznXuZTtyUDyye5X91MQkcq2fPlybr755qIrkl988UWef/55zEpaJvUvX5Y2jRREpDJlZGTQsWPHn7WoCA8Pr3IFAXxaFEREKsMPP/zA448/TosWLcjLyyM8PDykW1RUBF8ukWv6SESCLT09nfj4eLZv307Pnj0ZNmwYmzZtIj09nZiYmJBsUVERgtrmItjU5kJEKtr+/ft58sknGTduHJdddhmJiYm0bt3a61gVyrM2F8GiU1JFJBjmz59P48aNSUxM5NFHH2X9+vVVriCciC+LghaaRaQi7du3j3vvvZfY2FjOPfdcAoEAw4cP54wzzvA6WqXzZVEQEakIzjlmzJhBREQE7733Hi+++CJZWVk0b97c62ie8eVCs4jIqcrJyaFPnz7Mnz+f5s2bM3HiRCIjI72O5TlfjhS0piAiJys/P5/x48fTuHFj0tLSGDFiBKtXr1ZBKOTLkYJzbh4wLzo6Ot7rLCLiD4FAgJkzZ7JkyRLWrVtH69atSUxM5LLLLvM6WkjxZVEQESmPlStXctNNN3HkyBEABgwYwODBg6vkFcmnypfTRyIiZZWdnU2nTp2KCkJYWBi1a9dWQSiFioKIVEmHDh3ihRdeoFmzZhw8eLBatKioCL6cPlKbCxE5njVr1hAXF8fGjRvp0qULr776Ktu2bavyLSoqgtpciEiVcfDgQZ577jlGjhxJ3bp1GTt2LLfffrvXsULO8dpc+HKkICJyrCVLlhAfH8+OHTvo06cPL7/8MrVr1/Y6lu9oTUFEfO0///kP8fHx3HzzzdSoUYP09HTefPNNFYSTpKIgIr41d+5cGjduTFJSEk888QTr16/nxhtv9DqWr6koiIjv7N27l06dOnHnnXdSp04d1qxZw9ChQ/nFL37hdTTf82VRUJsLkerJOce0adOIiIhg9uzZDBw4kMzMTKKjS1wzlZPgy4VmtbkQqV4CgQDvv/8+K1asIBAI0KJFCyZOnEhERITX0aocXxYFEak+Vq1aRevWrYuuSO7bty+vvPIKYWFhHiermnw5fSQi1cO2bdvo0qXLz1pUXHTRRSoIQaSiICIhJzc3l2HDhnHllVeyb98+TjvtNLWoqCSaPhKRkPLxxx8TFxdHVlYW7du3Z/To0ezatUstKiqJioKIhIRDhw4xaNAgXn75Zc477zzeffdd7r77bsyMX/3qVyoGlURFQUQ8FwgEiIuLY/Pmzdx///2MGDGC888/3+tY1VJIrSmY2ZlmlmVmf/Q6i4gE3/fff0+/fv247rrrOHjwICkpKUyePFkFwUNBLQpmlmRme81swzHb25nZFjPbbmYDij31JPBOMDOJSGj417/+RZMmTXjttdd46KGH2LBhA+3atfM6VrUX7JHCJOBnf8tmFgaMBm4DIoDOZhZhZrcAm4Avg5xJRDz07bffEhcXx6233kqtWrVYsWIFo0aN4qyzzvI6mhDkNQXn3HIzq3/M5ubAdufcDgAzewu4E/h/wJkUFIofzSzZOZd/7D7NrCfQE+CSSy4JYnoRqWizZ8/mwQcfZN++fTz11FM8//zznH766V7HkmK8WGiuC3xe7HEOcI1zLgHAzLoCX5VUEACcc+OB8VBwk53gRhWRUxUIBJg3bx4ffPABS5cuJSoqigULFtCsWTOvo0kJvCgKJd0tu+gfd+fcpBPuQLfjFPGF1atX07p1aw4fPgxA7969ef311znttNM8Tial8eLsoxzg18Ue1wN2l2cHzrl5zrmeZ599doUGE5GKs2vXLrp27VpUEMLCwrjkkktUEEKcF0UhA7jczC41s3CgEzC3PDtQ62yR0JWfn8/o0aOJjIzk888/V4sKnwn2KakzgADQ0MxyzCzOOZcLJACLgM3AO865jeXZr0YKIqFpy5YttGrVioSEBK677jo++eQTli1bxsCBA0lLS9NVyT5gzvlvrbbYmkL8tm3bvI4jUu0dOXKE4cOH89JLL3HGGWfw6quvcv/992NW0hKieM3MspxzJd6ZKKSuaC4rjRREQse6deu45pprePrpp4mNjWXTpk088MADKgg+5cuioDUFEe/99NNPPP300/z+979n9+7dzJw5k3fffZdf/vKXXkeTU+DLoqCRgoi3Vq1aRVRUFEOGDOH+++9n8+bN3HXXXV7Hkgrgy6IgIt747rvv+Otf/8oNN9zAoUOHSE1NJSkpiXPPPdfraFJBfFkUNH0kUvkWLVpEZGQko0eP5q9//SvZ2dm0adPG61hSwXx5PwXn3DxgXnR0dLzXWUSqskAgQHJyMllZWaSkpHDFFVewcuVKrr32Wq+jSZD4siiISPAFAgFiYmKKrkh+4IEHGDt2rBrYVXG+nD4SkeDas2cP8fHxP2tR0bBhQxWEasCXRUFrCiLB4Zzjn//8JxEREWzbtk0tKqohXxYFnZIqUvF27txJ27Zt6d69O02aNCE7O1stKqohrSmIVHN5eXmMHj2ap59+GjPjzTffpFevXtSoUfA7o4pB9aKiIFKNbd68mbi4OAKBALfddhtjx47VHQ2rOV9OH2lNQeTUHDlyhL///e9ERUWxZcsWpk6dyoIFC1QQxJ9FQWsKIicvKyuL6Ohonn32Wdq3b8/mzZvp0qWLGtgJ4NOiICLl9+OPPzJgwACuueYa9u3bx+zZs3n77be58MILvY4mIURrCiLVwPLly+nRowfbtm0jLi6O4cOHc84553gdS0KQioJIFRUIBFi4cCEbNmxg1qxZXHrppSxevJibb77Z62gSwlQURKqgY1tUdOzYkaSkJM4880yPk0mo8+Wags4+Eind119/zYMPPlhUEGrUqEFUVJQKgpSJL4uCzj4S+W/OOd555x0aNWpEdnY2NWvWJCwsjFq1aqlFhZSZpo9EqoDdu3fz4IMPMmfOHKKjo1m8eDEHDx4kPT2dmJgYXZUsZaaiIOJjzjmSkpJ49NFHOXToEP/4xz/o168fNWsW/K+tYiDlpaIg4lM7duwgPj6eJUuWcOONNzJhwgQaNGjgdSzxOV+uKYhUZ3l5ebz66qs0adKEjIwMxo0bx5IlS1QQpEJopCDiIxs3biQuLo41a9Zw++23M3bsWOrVq+d1LKlCNFIQ8YHDhw/zt7/9jaZNm/Lpp58yffp05s2bp4IgFS5kRgpm1gjoC9QB0pxzYzyOJBISMjIyiIuLIzs7m3vvvZeRI0dywQUXeB1LqqigFgUzSwL+COx1zkUW294OeA0IAyY45152zm0GeptZDSAxmLlEQl0gECA1NZWtW7fy1ltvcfHFFzN37lxiY2O9jiZVXLBHCpOAN4ApRzeYWRgwGmgD5AAZZjbXObfJzO4ABhS+R6RaCgQCtG7dmkOHDgFw5513MnnyZHSxplSGoK4pOOeWA98cs7k5sN05t8M5dxh4C7iz8PVznXPXAveVtk8z62lmmWaWuW/fvmBFF/HE/v376du3b1FBqFGjBtdcc40KglQaLxaa6wKfF3ucA9Q1sxgze93MxgHJpb3ZOTfeORftnIvWvKpUJfPnz6dx48ZkZmaqRYV4xouF5pJu7+Scc+lAepl2YBYLxOq8bKkK9u3bR9++fZkxYwaRkZHMmjWLvLw8tagQT3hRFHKAXxd7XA/YXZ4dOOfmAfOio6PjKzKYSGVyzjFjxgwefvhhDhw4wEsvvcSAAQMIDw8H1KJCvOHF9FEGcLmZXWpm4UAnYG55dqDW2eJ3OTk53HHHHdx3331cdtllrF27lueff76oIIh4JahFwcxmAAGgoZnlmFmccy4XSAAWAZuBd5xzG8uzX7XOFr/Kz89n3LhxREREkJaWxogRI1i9ejWRkZEnfrNIJSjT9JGZnQ28CNxQuGkZ8Dfn3HF/VXfOdS5lezLHWUwuQx6tKYjvbN++nfj4eNLT02ndujWJiYlcdtllXscS+ZmyjhSSgANAx8KvA8A/gxXqRDRSED/Jzc1l+PDhNGnShLVr15KYmEhaWpoKgoSksi40X+acu7vY45fM7KMg5BGpUtavX09cXByZmZnExsYyZswY6tat63UskVKVtSj8aGbXO+dWApjZdcCPwYt1fJo+klAWCARYvHgxn332GVOnTuXcc8/lrbfeomPHjpiVdEa2SOgoa1HoDUwpXFswCq5S7hqsUCeiU1IlVB3bouLWW29l2rRp1KlTx+NkImVTpqLgnPsYuMrMahc+PhDUVCegkYKEooMHD/LYY4/9rEVFTEyMCoL4SlnPPqoF3A3UB2oeHQI75/4WtGTHoZGChJq0tDTi4+P57LPPCAsLAyA8PFwtKsR3yjp9NAfYD2QBh4IXR8Rf/vOf//DYY48xceJEGjRoQHp6OuHh4WpRIb5V1qJQzznXLqhJRHxmzpw59OnThy+//JInnniCF198kV/84heAWlSIf5X1OoXVZtYkqEnKQW0uxEtffvklf/7zn2nfvj116tRhzZo1DB06tKggiPjZcYuCmWWb2QbgZmCtmW0xs/WF29dXTsT/povXxAvOOf7nf/6HiIgI3n//fQYOHEhmZibR0dFeRxOpMCeaPvojBYUjC9CpPlJt/fvf/6Z3796kpKTQokULJk6cSEREhNexRCrccYuCc24XFDW2u9A5l1EpqURCxNEGdk888QT5+fmMHDmShISEojOMRKqasi40twZ6mdku4CAFF7A559yVQUsm4rGtW7fSo0cPVqxYwS233ML48eO59NJLvY4lElRlLQq3BTVFOeniNQmmFStWMHjwYJYsWcIZZ5xBUlISXbt2VYsKqRbMOed1hpMWHR3tMjMzvY4hVciUKVPo2rUrzjnCwsKYPXs2sbGxXscSqVBmluWcK/EMCS/uvCYScg4dOsRzzz1Ht27dKP6L0oYNGzxMJVL5VBSk2gsEAkRFRTFo0CDatm3L6aefTlhYmNpUSLWkoiDV1vfff0/fvn257rrr+OGHH0hJSSE5OZklS5YwcOBA0tLSdGWyVDtlXWgWqVJSU1Pp2bMnu3btIiEhgcGDB3PWWWcBBS0qVAykuvLlSEFtLuRkffvtt3Tr1q1ommjFihWMGjWqqCCIVHe+LApqcyEnY9asWURERDB16lSeeuopPvroI66//nqvY4mEFE0fSZX3xRdfkJCQwMyZM4mKimLBggU0a9bM61giIcmXIwWRsnDOMWnSJCIiIpg/fz6DBw/mww8/VEEQOQ6NFKRK2rlzJ7169SI1NZXrrruOCRMmcMUVV3gdSyTkaaQgVcqqVauIjY2lUaNGrFq1ilGjRrF8+XIVBJEy0khBqowZM2bQpUsX8vPzqVGjBu+++y533XWX17FEfEUjBfG9I0eOMHjwYP7yl7+Qn58PgJmxZcsWj5OJ+E9IFQUza29miWY2x8xu9TqPhL5169bRvHlznnnmGVq1aqUWFSKnKOhFwcySzGxv4W09i29vV3h7z+1mNgDAOfe+cy4e6Ar8OdjZxL9++uknnnrqKX7/+9/zxRdfMHPmTJYsWaIWFSKnKOits82sFfA9MMU5F1m4LQzYCrQBcoAMoLNzblPh868A05xza4+3b7XOrp5WrlxJXFwcW7dupVu3brzyyiuce+65XscS8Q1PW2c755YD3xyzuTmw3Tm3wzl3GHgLuNMKDAVSSisIZtbTzDLNLHPfvn3BDS8h5bvvviMhIYEbbriBw4cPk5qaSlJSkgqCSAXyak2hLvB5scc5hdv+CtwCdDCz3iW90Tk33jkX7ZyLvuCCC4KfVELCwoULiYyM5M033+Thhx8mOzubNm3aeB1LpMrx6pTUku5r6JxzrwOvn/DNuh1ntfH111/zyCOPMGXKFK644gpWrlzJtdde63UskSrLq5FCDvDrYo/rAbs9yiIhyDnHe++9R0REBNOnT+eZZ55h3bp1KggiQeZVUcgALjezS80sHOgEzC3rm9UltWrbs2cPd999N/fccw/16tUjIyODQYMGcfrpp3sdTaTKq4xTUmcAAaChmeWYWZxzLhdIABYBm4F3nHMbg51FQtvq1avp0KEDv/vd70hJSWHo0KGsWbOGqKgor6OJVBtBX1NwznUuZXsykHwy+9SaQtUzc+ZMOnbsWNSiYvr06fz5z7pURaSyhdQVzWWl6aOqIy8vj9dee43OnTv/rEXFjh07PE4mUj35sijodpxVw6ZNm7jhhhvo168fV199tVpUiIQAXxYFjRT87ciRIwwaNIimTZuydetWpk6dyurVq9WiQiQEqHW2VKqsrCy6d+/O+vXr6dixI6NGjeLCCy8EoGXLlioGIh7z5UhB00f+8+OPP/Lkk0/SvHlz9u3bx+zZs3n77beLCoKIhAZfFgVNH/nLsmXLuOqqqxg2bBjdu3dn06ZNtG/f3utYIlICXxYF8YcDBw7Qp08fYmJiyM3NZfHixSQmJnLOOed4HU1ESuHLoqDpo9CXnJxM48aNGTduHP379yc7O5ubb77Z61gicgK+LAqaPgpdX331FV26dOH222+ndu3arF69mhEjRnDmmWd6HU1EysCXRUFCz+rVq+ncuTOXX345b7/9Ns8//zxr166lRYsWXkcTkXLQKalyyubOncuf/vQn8vPzMTMmT57MX/7yF69jichJ0EhBTppzjgkTJhT1LAKoUaMGOTk5HicTkZPly6KghWbvffrpp9xyyy3Ex8fTqFEjtagQqSJ8WRS00OydvLw8RowYQZMmTcjIyGDcuHFkZWWpRYVIFaE1BSmzDRs2EBcXx4cffsjtt9/O2LFjqVevHqAWFSJVhS9HClK5Dh8+zEsvvUSzZs3YsWMH06dPZ968eUUFQUSqDo0U5LgyMjLo3r07GzZs4N5772XkyJFccMEFXscSkSDRSEFK9MMPP/DYY4/RokULvv32W+bOncu0adNUEESqOF8WBZ19FFxLly7lyiuv5JVXXiE+Pp6NGzcSGxvrdSwRqQS+LAo6+yg49u/fT69evbjpppsAWLJkCWPHjkXHWaT60JqCEAgEGDNmDCkpKXzzzTc89thjvPTSS5xxxhleRxORSqaiUM0lJydzxx13kJeXh5mRmJhIXFyc17FExCO+nD6SU+ecY/r06dxzzz3k5eUBBS0q9u7d63EyEfGSikI1lJOTwx133MF9991H/fr11aJCRIpo+qgayc/PJzExkccff5zc3FxGjBjBww8/zIcffkh6ejoxMTG6KlmkmlNRqCa2b99OfHw86enp3HTTTSQmJvLb3/4WUIsKEfk/ITN9ZGa/NbOJZvae11mqktzcXIYPH06TJk1Yu3YtiYmJLF68uKggiIgUF9SiYGZJZrbXzDYcs72dmW0xs+1mNgDAObfDOafTXipQdnY21157LY8//ji33normzZtokePHpiZ19FEJEQFe6QwCWhXfIOZhQGjgduACKCzmUUEOUe1cujQIV544QWaNWvGzp07efvtt3n//fepW7eu19FEJMQFdU3BObfczOofs7k5sN05twPAzN4C7gQ2lWWfZtYT6AlwySWXVFzYKuKDDz4gLi6OTZs20aVLF0aOHMn555/vdSwR8Qkv1hTqAp8Xe5wD1DWz881sLNDUzJ4q7c3OufHOuWjnXLSas/2fgwcP8sgjj3Dttddy4MABFixYwNSpU1UQRKRcvDj7qKQJbeec+xroXaYdmMUCsQ0aNKjQYH4UCASYMGECKSkp7Nmzhz59+vDyyy9Tu3Ztr6OJiA95URRygF8Xe1wP2F2eHTjn5gHzoqOj4ysymN+kpqbyhz/8oahFxejRo3nwwQe9jiUiPubF9FEGcLmZXWpm4UAnYG55dqDW2TBnzhw6dOjwsxYV1fl4iEjFCPYpqTOAANDQzHLMLM45lwskAIuAzcA7zrmN5dlvdW6dvXfvXjp16kT79u256KKLqFWrllpUiEiFCfbZR51L2Z4MJJ/sfqvjmoJzjmnTptG3b1++//57Bg0axBNPPEFmZqZaVIhIhTHnnNcZTlp0dLTLzMz0OkbQ/fvf/6Z3796kpKTQsmVLJk6cSKNGjbyOJSI+ZWZZzrnokp4LmTYX8t/y8/MZM2YMjRs3ZtmyZbz22musWLFCBUFEgsaXDfGqw/TR1q1b6dGjBytWrKBNmzaMHz+e+vXrex1LRKo4X44UqvJCc25uLsOGDeOqq64iOzubSZMmsWjRIhUEEakUGimEkI8//pju3buzdu1a/vSnPzF69Gguvvhir2OJSDWikUII+Omnn3j22WeJjo7mf//3f3nvvfeYNWuWCoKIVDpfjhSqikAgwOTJk1m4cCG7du3igQceYMSIEZx33nleRxORakpFwSNpaWm0a9eO3NxczIwRI0bQv39/r2OJSDXny+kjv7e5SE1NpUOHDuTm5gIFLSp++uknj1OJiPi0KPh1TeHbb7+lW7dutG3bltq1a6tFhYiEHF8WBT+aNWsWERERTJ06laeeeootW7awdOlSBg4cSFpamlpUiEhI0JpCkH3xxRckJCQwc+ZMmjZtSnJyMk2bNgWgZcuWKgYiElI0UggS5xyTJ08mIiKC+fPnM2TIENasWVNUEEREQpEvi0KoLzTv3LmTdu3a0bVrVxo3bszHH3/MgAEDOO2007yOJiJyXL4sCqG60Jyfn8+oUaOIjIxk9erVvPHGGyxbtoyGDRt6HU1EpEy0plBBPvnkE3r06MGqVato27Yt48aN4ze/+Y3XsUREysWXI4VQcuTIEQYPHsxVV13Fpk2bmDx5MikpKSoIIuJLGimcpEAgwPTp01m0aBHbtm2jQ4cOvPHGG1x00UVeRxMROWkqCichPT2dNm3aFF2RPGTIEAYMGOBxKhGRU+fL6SMvzz5auXIl99xzT1FBCAsLw8+3NBURKc6XRcGLs4++++47EhISuOGGG6hZs6ZaVIhIlaTpozJYuHAhvXr14vPPP6dv374MGjSI7Oxs0tPTiYmJ0VXJIlJlqCgcx9dff80jjzzClClTaNSoEatWrSoqAGpRISJVkS+nj4LNOcd7771HREQE06dP59lnn2XdunUqAiJS5WmkcIw9e/bw0EMPMXv2bK6++mpSU1O56qqrvI4lIlIpNFIo5JwjKSmJRo0akZKSwtChQ/nggw9UEESkWtFIAfjss8/o2bMnixcvplWrViQmJvK73/3O61giIpUuZEYKZnammU02s0Qzu68yfmZeXh6vvfYakZGRrFmzhjFjxrB06VIVBBGptoJaFMwsycz2mtmGY7a3M7MtZrbdzI5eCnwX8J5zLh64I5i5AoEA/fv3Jyoqin79+nHjjTeyceNGevfuTY0aIVMnRUQqXbCnjyYBbwBTjm4wszBgNNAGyAEyzGwuUA/ILnxZXrACBQIBbrzxRo4cOQLACy+8wAsvvICZBetHioj4RlB/LXbOLQe+OWZzc2C7c26Hc+4w8BZwJwUFot6JcplZTzPLNLPMffv2lTtTeno6eXkFNScsLIxatWqpIIiIFPJirqQu8HmxxzmF22YBd5vZGGBeaW92zo13zkU756IvuOCCcv/wmJgYtagQESmFF2cflfRruXPOHQS6lWkHZrFAbIMGDcr9w1u2bElaWppaVIiIlMCLopAD/LrY43rA7soMoBYVIiIl82L6KAO43MwuNbNwoBMwtzw7CNV7NIuI+F2wT0mdAQSAhmaWY2ZxzrlcIAFYBGwG3nHObQxmDhERKZugTh855zqXsj0ZSD7Z/Z7KmoKIiJTOl1dqafpIRCQ4fFkUvLwdp4hIVebLoqCRgohIcJifbzpvZvuAXV7nOIE6wFdehygD5axYfskJ/smqnBXnN865Eq/+9XVR8AMzy3TORXud40SUs2L5JSf4J6tyVg5fTh+JiEhwqCiIiEgRFYXgG+91gDJSzorll5zgn6zKWQm0piAiIkU0UhARkSIqCiIiUkRF4RSZ2a/NbKmZbTazjWbWt4TXxJjZfjP7qPDreS+yFmbZaWbZhTkyS3jezOz1wvtnrzezZh5kbFjsWH1kZgfMrN8xr/HkmJZ033EzO8/M/mVm2wr/PLeU95Z0b/LKzvoPM/uk8O92tpmdU8p7j/s5qYScL5rZ/xb7+/1DKe+ttGNaSs63i2XcaWYflfLeSjuep8w5p69T+AIuBpoVfn8WsBWIOOY1McB8r7MWZtkJ1DnO838AUii4GVILYI3HecOALyi42MbzYwq0ApoBG4ptGwYMKPx+ADC0lP+OT4HfAuHAx8d+Tiop661AzcLvh5aUtSyfk0rI+SLwWBk+G5V2TEvKeczzrwDPe308T/VLI4VT5Jzb45xbW/j9dxS0A6/rbapTcicwxRX4ADjHzC72MM/NwKfOuZC4ct2VfN/xO4HJhd9PBtqX8NbS7k0eNCVldc6luoL29QAf8H/3RfdMKce0LCr1mB4vpxXc6L0jMCNYP7+yqChUIDOrDzQF1pTwdEsz+9jMUsysceUm+xkHpJpZlpn1LOH50u6h7ZVOlP4/Wqgc04ucc3ug4JcE4MISXhNqxxWgOwWjwpKc6HNSGRIKp7mSSpmSC6VjegPwpXNuWynPh8LxLBMVhQpiZv8PmAn0c84dOObptRRMf1wFjALer+R4xV3nnGsG3AY8ZGatjnm+xHtoBz/Wfyu8M98dwLslPB1Kx7QsQua4ApjZM0AuMK2Ul5zocxJsY4DLgChgDwVTM8cKpWPameOPErw+nmWmolABzOw0CgrCNOfcrGOfd84dcM59X/h9MnCamdWp5JhHs+wu/HMvMJuCIXhxnt9Du5jbgLXOuS+PfSKUjinw5dEptsI/95bwmpA5rmb2APBH4D5XOOF9rDJ8ToLKOfelcy7POZcPJJby80PimJpZTeAu4O3SXuP18SwPFYVTVDiXOBHY7JwbUcprfln4OsysOQXH/evKS1mU40wzO+vo9xQsOm445mVzgfsLz0JqAew/OjXigVJ/+wqVY1poLvBA4fcPAHNKeM0p35u8IphZO+BJ4A7n3A+lvKYsn5OgOmYd60+l/PyQOKbALcAnzrmckp4MheNZLl6vdPv9C7iegiHreuCjwq8/AL2B3oWvSQA2UnB2xAfAtR5l/W1hho8L8zxTuL14VgNGU3BWRzYQ7VHWMyj4R/7sYts8P6YUFKk9wBEKflONA84H0oBthX+eV/jaXwHJxd77BwrOTvv06LH3IOt2Cubhj35Wxx6btbTPSSXnnFr4+VtPwT/0F3t9TEvKWbh90tHPZbHXenY8T/VLbS5ERKSIpo9ERKSIioKIiBRRURARkSIqCiIiUkRFQUREiqgoiIhIERUFEREpoqIgUoHMrL4V3Fsj0Qrur5FqZr/wOpdIWakoiFS8y4HRzrnGwH+Au72NI1J2KgoiFe8z59xHhd9nAfW9iyJSPioKIhXvULHv84CaXgURKS8VBRERKaKiICIiRdQlVUREimikICIiRVQURESkiIqCiIgUUVEQEZEiKgoiIlJERUFERIqoKIiISJH/D4ujd9zobg5AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "nvals = 1+np.arange(19)\n",
    "rho_vals = np.zeros_like(nvals,dtype=np.float64)\n",
    "for i,n in enumerate(nvals):\n",
    "    P,L,U = lu(matrix_with_large_rho(n))\n",
    "    rho_vals[i] = np.max(np.abs(U))\n",
    "plt.semilogy(nvals,rho_vals,'k.-')\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('rho')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The growth factor $\\rho$ grows exponentially with $n$ for this matrix. What is the effect on the backward error for solving a linear system of equations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.05161206801240477\n"
     ]
    }
   ],
   "source": [
    "n = 100\n",
    "A = matrix_with_large_rho(n)\n",
    "b = np.random.rand(n)\n",
    "x = np.linalg.solve(A,b)\n",
    "print(np.linalg.norm(np.dot(A,x)-b)/(np.linalg.norm(A)*np.linalg.norm(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44.80225124630289"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Gaussian elimination with partial pivoting is not backward stable.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Growth factor of random matrices\n",
    "\n",
    "The following plot shows the probability density functions for the growth factor $\\rho$ for matrices of dimension $8$, $16$, and $32$ (see [Driscoll, Maki '07](http://www.math.udel.edu/~driscoll/papers/2007-DriscollMaki-673.pdf)). It shows that large growth factors are exponentially unlikely to occur randomly. Indeed, essentially no practical application is known, where growth factors cause difficulties. Hence, **in practice Gaussian elimination with partial pivoting is backward stable**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"./growth_factors.png\">"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python [conda env:dev]",
   "language": "python",
   "name": "conda-env-dev-py"
  },
  "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}