12_LA_conditioning_stability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>\n",
    "\n",
    "Note:  This material largely follows the text \"Numerical Linear Algebra\" by Trefethen and Bau (SIAM, 1997) and is meant as a guide and supplement to the material presented there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "init_cell": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Conditioning and Stability\n",
    "\n",
    "Conditioning is the behavior of a perturbed problem mathematically (analytically).  Stability is concerned with how an algorithm behaves when perturbed (say with input)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Conditioning and Condition Numbers\n",
    "\n",
    "A **well-conditioned** problem is one where a small perturbation to the original problem leads to only small changes in the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Formally we can think of a function $f$ which maps $x$ to $y$\n",
    "\n",
    "$$\n",
    "    f(x) = y \\quad \\text{or} \\quad f: X \\rightarrow Y.\n",
    "$$\n",
    "\n",
    "Using $\\epsilon$ notation, let $x \\in X$ where we perturb $x$ with $\\epsilon$ and we ask how the result $y$ changes:\n",
    "\n",
    "$$\n",
    "    ||f(x) - f(x + \\epsilon)|| \\leq C ||x - (x+\\epsilon)||\n",
    "$$\n",
    "\n",
    "for some constant $C$ possible dependent on $\\epsilon$ depending on the type of conditioning we are considering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Absolute Condition Number\n",
    "\n",
    "If we let $\\delta x = x + \\epsilon$ be the small perturbation to the input and $\\delta f = f(x + \\delta x) - f(x)$ be the result the **absolute condition number** $\\hat{~\\kappa}$ can be defined as\n",
    "\n",
    "$$\n",
    "    \\hat{\\!\\kappa} = \\sup_{\\delta x} \\frac{||\\delta f||}{||\\delta x||}\n",
    "$$\n",
    "\n",
    "for most problems (assuming $\\delta f$ and $\\delta x$ are both infinitesimal).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "When $f$ is differentiable we can evaluate the condition number via the Jacobian.  Recall that the derivative of a multi-valued function can be termed in the form of a Jacobian $J(x)$ where\n",
    "$$\n",
    "    [J(x)]_{ij} = \\frac{\\partial f_i}{\\partial x_j}(x).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This allows us to write the infinitesimal $\\delta f$ as\n",
    "$$\n",
    "    \\delta f \\approx J(x) \\delta x\n",
    "$$\n",
    "with equality when $||\\delta x|| \\rightarrow 0$.  Then we can write the condition number as\n",
    "$$\n",
    "    \\hat{\\!\\kappa} = ||J(x)||\n",
    "$$\n",
    "where the norm is the one induced by the spaces $X$ and $Y$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Relative Condition Number\n",
    "\n",
    "The **relative condition number** is defined similarly and is related to the difference before between the absolute error and relative error as defined previously.  With the same caveats as before it can be defined as\n",
    "$$\n",
    "    \\kappa = \\sup_{\\delta x} \\left( \\frac{\\frac{||\\delta f||}{||f(x)||}}{\\frac{||\\delta x||}{||x||}} \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again if $f$ is differentiable we can use the Jacobian $J(x)$ to evaluate the relative condition number as\n",
    "$$\n",
    "    \\kappa = \\frac{||J(x)||}{||f(x)|| ~/ ~||x||}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Examples\n",
    "Calculate the following relative condition numbers of the following problems.\n",
    "\n",
    "$\\sqrt{x}$ for $x > 0$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    f(x) = \\sqrt{x}, \\quad J(x) = f'(x) = \\frac{1}{2 \\sqrt{x}} \\\\\n",
    "    \\kappa = \\frac{||J(x)||}{||f(x)|| / ||x||} = \\frac{1}{2 \\sqrt{x}} \\frac{x}{\\sqrt{x}} = \\frac{1}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Obtain the scalar $f(x) = x_1 - x_2$ from the vector $x = (x_1, x_2)^T \\in \\mathbb R^2$ using an $\\ell_\\infty$ norm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    f(x) = x_1 - x_2, \\quad J(x) = \\left [ \\frac{\\partial f}{\\partial x_1}, \\frac{\\partial f}{\\partial x_2}\\right ] = [1, -1] \\\\\n",
    "    \\kappa = \\frac{||J(x)||_\\infty}{||f(x)||_\\infty / ||x||_\\infty} = \\frac{2 \\max_{i=1,2} |x_i|}{|x_1 - x_2|}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "    ||J||_\\infty = 2.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Matrix-vector multiplication $Ax$ only perturbing $x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    \\kappa &= \\sup_{\\delta x} \\left ( \\frac{||A (x+\\delta x) - A x||}{||Ax||} \\frac{||x||}{||\\delta x||}\\right ) \\\\\n",
    "    &= \\sup_{\\delta x} \\frac{ ||A \\delta x||}{||\\delta x||} \\frac{||x||}{||Ax||} \\\\\n",
    "    &= ||A|| \\frac{||x||}{||A x||}\n",
    "\\end{aligned}$$\n",
    "\n",
    "If $A$ has an inverse, then we can use\n",
    "$$\n",
    "    \\frac{||x||}{||A x||} \\leq ||A^{-1}||\n",
    "$$\n",
    "and therefore\n",
    "$$\n",
    "    \\kappa \\leq ||A|| ||A^{-1}||.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Condition Number of a Matrix\n",
    "\n",
    "The condition number of a matrix is defined by the product\n",
    "$$\n",
    "    \\kappa(A) = ||A||~||A^{-1}||.\n",
    "$$\n",
    "where here we are thinking about the matrix rather than a problem.  If $\\kappa$ is small than $A$ is said to be **well-conditioned**.  If $A$ is singular we assign $\\kappa(A) = \\infty$ as the matrix's condition number."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "When we are considering the $\\ell_2$ norm then we can write the condition number as\n",
    "$$\n",
    "    \\kappa(A) = \\frac{\\sqrt{\\rho(A^\\ast A)}}{\\sqrt{\\rho((A^\\ast A)^{-1})}} = \\frac{\\sqrt{\\max |\\lambda|}}{\\sqrt{\\min |\\lambda|}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Condition Number of a System of Equations\n",
    "\n",
    "Another way to think about the conditioning of a problem we have looked at before is that the matrix $A$ itself is an input to the problem.  Consider than the system of equations $Ax = b$ where we will perturb both $A$ and $x$ resulting in\n",
    "$$\n",
    "    (A + \\delta A)(x + \\delta x) = b.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Assuming we solve the problem exactly we know that $Ax = b$ and that the infinitesimals multiplied $\\delta A \\delta x$ is small the above simplifies to\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    (A + \\delta A)(x + \\delta x) &= b \\\\\n",
    "    Ax + \\delta Ax + A \\delta x + \\delta A \\delta x &= b \\\\\n",
    "    \\delta Ax + A \\delta x & = 0\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Solving for $\\delta x$ leads to \n",
    "$$\n",
    "    \\delta x = -A^{-1} \\delta A x\n",
    "$$\n",
    "implying\n",
    "\n",
    "$$\n",
    "    ||\\delta x|| \\leq ||A^{-1}|| ~ ||\\delta A|| ~ ||x||\n",
    "$$\n",
    "\n",
    "and therefore\n",
    "$$\n",
    "    \\frac{\\frac{||\\delta x||}{||x||}}{\\frac{||\\delta A||}{||A||}} \\leq ||A^{-1}||~||A|| = \\kappa(A).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also say the following regarding the condition number of a system of equations then\n",
    "\n",
    "**Theorem:**  Let $b$ be fixed and consider the problem of computing $x$ in $Ax = b$ where $A$ is square and non-singular.  The condition number of this problem with respect to perturbations in $A$ is the condition number of the matrix $\\kappa(A)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Stability\n",
    "\n",
    "We now return to the consideration of the fact that we are interested not only in the well-conditioning of a mathematical problem but in how we might solve it on a finite precision machine.  In some sense conditioning describes how well we can solve a problem in exact arithmetic and stability how well we can solve the problem in finite arithmetic.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Accuracy and Stability\n",
    "\n",
    "As we have defined before we will consider **absolute error** as\n",
    "$$\n",
    "    ||\\hat{~f}(x) - f(x)||\n",
    "$$\n",
    "where $\\hat{~f(x)}$ is the approximation to the true solution $f(x)$.  Similarly we can define **relative error** as\n",
    "$$\n",
    "    \\frac{||\\hat{~f}(x) - f(x)||}{||f(x)||}.\n",
    "$$\n",
    "In the ideal case we would like the relative error to be $\\mathcal{O}(\\epsilon_{\\text{machine}})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Forwards Stability\n",
    "\n",
    "A **forward stable** algorithm for $x in \\in X$ has\n",
    "\n",
    "$$\n",
    "    \\frac{||\\hat{~f}(x) - f(x)||}{||f(x)||} = \\mathcal{O}(\\epsilon_{\\text{machine}})\n",
    "$$\n",
    "\n",
    "In other words\n",
    "> A forward stable algorithm gives almost the right answer to exactly the right question."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Backwards Stability\n",
    "\n",
    "A stronger notion of stability can also be defined which is satisfied by many approaches in numerical linear algebra.  We say that an algorithm $\\hat{~f}$ is **backward stable** if for $x \\in X$ we have\n",
    "\n",
    "$$\n",
    "    \\hat{~f}(x) = f(\\hat{\\!x})\n",
    "$$\n",
    "\n",
    "for some $\\hat{\\!x}$ with\n",
    "\n",
    "$$\n",
    "    \\frac{||\\hat{\\!x} - x||}{||x||} = \\mathcal{O}(\\epsilon_{\\text{machine}}).\n",
    "$$\n",
    "\n",
    "In other words\n",
    "> A backward stable algorithm gives exactly the right answer to nearly the right question."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Combining these ideas along with the idea that we should not expect to be able to accurately compute the solution to a poorly conditioned problem we can form the mixed forward-backward sense of stability as for $x \\in X$ if\n",
    "\n",
    "$$\n",
    "    \\frac{||\\hat{~f}(x) - f(\\hat{~x})||}{||f(\\hat{\\!x})||} = \\mathcal{O}(\\epsilon_{\\text{machine}})\n",
    "$$\n",
    "\n",
    "for some $\\hat{\\!x}$ with \n",
    "\n",
    "$$\n",
    "    \\frac{||\\hat{\\!x} - x||}{||x||} = \\mathcal{O}(\\epsilon_{\\text{machine}}).\n",
    "$$\n",
    "\n",
    "In other words\n",
    "> A stable algorithm gives nearly the right answer to nearly the right question."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "An important aspect of the above statement is that we can not necessarily guarantee an accurate result.  If the condition number $\\kappa(x)$ is small we would expect that a stable algorithm would give us an accurate result (by definition).  This is reflected in the following theorem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Theorem:**  Suppose a backward stable algorithm is applied to solve a problem $f: X \\rightarrow Y$ with condition number $\\kappa$ on a finite precision machine, then the relative errors satisfy\n",
    "$$\n",
    "    \\frac{||\\hat{~f}(x) - f(\\hat{\\!x})||}{||f(\\hat{\\!x})||} = \\mathcal{O}(\\kappa(x) ~ \\epsilon_{\\text{machine}}).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Proof:**  By the definition of the condition number of a problem we can write\n",
    "$$\n",
    "    \\frac{||\\hat{~f}(x) - f(\\hat{\\!x})||}{||f(\\hat{\\!x})||} \\leq (\\kappa(x) + \\mathcal{O}(\\epsilon_{\\text{machine}}))\\frac{||\\hat{\\!x} - x||}{||x||}.\n",
    "$$\n",
    "Combining this with the definition of backwards stability we can arrive at the statement of the theorem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The proof above demonstrates **backward error analysis**, in other words using the condition number of the problem and stability of the algorithm to determine the error.  A perhaps more obvious approach to determine eventual accuracy would be to consider the accrual of error at each step of an algorithm given slightly perturbed input.  This approach is known as **forward error analysis**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Stability of $Ax = b$ using Householder Triangularization\n",
    "\n",
    "As an example lets consider the conditioning and algorithm for solving $Ax = b$.  Here we will use a $QR$ factorization approach to solve $Ax = b$ given by a Householder triangularization.  First off lets discuss the $QR$ factorization itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Theorem:** Let the $QR$ factorization $A = QR$ of a matrix $A \\in \\mathbb C^{m \\times n}$ be computed using a Householder triangularization approach on a finite precision machine, then\n",
    "\n",
    "$$\n",
    "    \\hat{\\!Q}~\\hat{\\!R} = A + \\delta A \\quad \\frac{||\\delta A||}{||A||} = \\mathcal{O}(\\epsilon_{\\text{machine}})\n",
    "$$\n",
    "\n",
    "for some $\\delta A \\in \\mathbb C^{m \\times n}$ where $~\\hat{\\!Q}$ and $\\hat{~\\!R}$ are the finite arithmetic versions of $Q$ and $R$.  Householder triangularization is therefore backward stable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Solving $Ax = b$ with $QR$ Factorization\n",
    "\n",
    "So Householder triangularization is backwards stable but we also know that this does not guarantee accuracy if the problem itself is ill-conditioned.  Is backward stability enough to guarantee accurate results if we use it for $Ax = b$ for instance?  It turns out that the accuracy of the product of $QR$ is enough to guarantee accuracy of a larger algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Consider the steps to solving $A x = b$ using $QR$ factorization:\n",
    "1. Compute the $QR$ factorization of $A$\n",
    "1. Multiply the vector $b$ by $Q^\\ast$ so that $y = Q^\\ast b$.\n",
    "1. Solve using backward-substitution the triangular system $R x = y$ or $x = R^{-1} y$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We know that step (1) is backward stable, what about step (2), the matrix-vector multiplication?  We can write the estimate of the backwards stability as\n",
    "\n",
    "$$\n",
    "    (\\hat{\\!Q} + \\delta Q) ~\\hat{\\!y} = b \\quad \\text{with} \\quad ||\\delta Q|| = \\mathcal{O}(\\epsilon_{\\text{machine}})\n",
    "$$\n",
    "\n",
    "where we have inverted the matrix $~\\hat{\\!Q}$ since it is unitary.  Since this is exact we know also that the matrix-vector multiplication is also backwards stable since this is an equivalent statement to multiplying $b$ by a slightly perturbed matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Step (3) is backward substitution (or the computation of $R^{-1}$).  Writing the backwards stability estimate we have\n",
    "\n",
    "$$\n",
    "    (\\hat{\\!R} + \\delta R) ~\\hat{\\!x} = \\hat{\\!y} ~~~~\\text{with}~~~~ \\frac{||\\delta R||}{||\\hat{\\!R}||} = \\mathcal{O}(\\epsilon_{\\text{machine}})\n",
    "$$\n",
    "\n",
    "demonstrating that the results $\\hat{x}$ is the exact solution to a slight perturbation of the original problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "These results lead to the following two theorems:\n",
    "\n",
    "**Theorem:**  Using $QR$ factorization to solve $Ax=b$ as described above is backward stable, satisfying\n",
    "$$\n",
    "    (A + \\Delta A) \\hat{\\!x} = b, \\quad \\frac{||\\Delta A||}{||A||} = \\mathcal{O}(\\epsilon_{\\text{machine}})\n",
    "$$\n",
    "for some $\\Delta A \\in \\mathbb C^{m \\times n}$.\n",
    "\n",
    "**Theorem:** The solution $\\hat{x}$ computed by the above algorithm satisfies\n",
    "$$\n",
    "    \\frac{||\\hat{\\!x} - x||}{||x||} = \\mathcal{O}(\\kappa(x) ~ \\epsilon_{\\text{machine}}).\n",
    "$$"
   ]
  }
 ],
 "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.6.3"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}