linear_solve_stability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Backward stability of solving $Ax=b$ by Householder QR factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we demonstrate the backward stability of solving $Ax=b$ by the following algorithm:\n",
    "\n",
    "1. Factor $A=QR$ using Householder triangularization.\n",
    "2. Compute $y=Q^* b$.\n",
    "3. Solve $Ry = x$ by backward substitution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's implement Householder (*note: we don't need to actually compute $Q$ in the algorithm above, so we could just apply $Q^*$ directly using algorithm 10.2 from Trefethen & Bau*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def householder(A):\n",
    "    \"\"\"QR factorization via Householder triangularization.\"\"\"\n",
    "    m, n = A.shape\n",
    "    V = np.zeros(A.shape)\n",
    "    R = A.copy()\n",
    "    for k in range(n-1):\n",
    "        x = R[k:,k].copy()\n",
    "        x[0] = x[0] + np.sign(x[0])*np.linalg.norm(x,2)\n",
    "        x = x/np.linalg.norm(x,2)\n",
    "        V[k:,k] = x.copy()\n",
    "        for j in range(k,n):\n",
    "            R[k:,j] = R[k:,j] - 2*V[k:,k]*np.dot(V[k:,k].T,R[k:,j])\n",
    "    return V,R[:n,:]\n",
    "\n",
    "def apply_Q(V,x):\n",
    "    \"\"\"Algorithm 10.3 of Trefethen & Bau.\"\"\"\n",
    "    m, n = V.shape\n",
    "    for k in range(n-1,-1,-1):\n",
    "        x[k:] = x[k:] - 2*np.dot(V[k:,k],x[k:])*V[k:,k]\n",
    "    return x\n",
    "\n",
    "def compute_Q(V):\n",
    "    \"\"\"Find Q given the Householder reflector vectors.\"\"\"\n",
    "    m, n = V.shape\n",
    "    Q = np.zeros((m,n))\n",
    "    for k in range(n):\n",
    "        x = np.zeros(m)\n",
    "        x[k] = 1.\n",
    "        Q[:,k] = apply_Q(V,x)\n",
    "    return Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to show stability, let's pick a really nasty matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Hilbert(n):\n",
    "    \"\"\"Return the n x n Hilbert matrix.\"\"\"\n",
    "    A = np.zeros([n,n])\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            A[i,j] = 1./(i+j+1)\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=10\n",
    "x = np.ones(N)\n",
    "A = Hilbert(N)\n",
    "b = np.dot(A,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kappa = np.linalg.cond(A)\n",
    "np.log10(kappa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The condition number is about $10^{13}$, which means we could lose 13 digits of accuracy just due to rounding of the inputs.  Hence a backward stable algorithm can be expected to give us $16-13 = 3$ accurate digits in double precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V, R = householder(A)\n",
    "Q = compute_Q(V)\n",
    "y = np.dot(Q.T,b)\n",
    "x_tilde = np.linalg.solve(R,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.linalg.norm(x_tilde-x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, we are left with about 3 significant accurate digits; in other words, the forward error is about $10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(x_tilde)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How large is the backward error?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm = np.linalg.norm\n",
    "r = np.dot(A,x_tilde)-b\n",
    "print(norm(r)/norm(A)*norm(x_tilde))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Still on the order of unit roundoff, as we expect."
   ]
  }
 ],
 "metadata": {
  "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.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}