notes/Is-Gaussian-Elimination-Unstable.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b4738d86",
   "metadata": {},
   "source": [
    "Gaussian elimination for $A = LU$ may not be backward stable when $||U|| >> ||A||$. Here is a special matrix for which this occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6103f719",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "93919aaa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8×8 Matrix{Float64}:\n",
       "  1.0   0.0   0.0   0.0   0.0   0.0   0.0  1.0\n",
       " -1.0   1.0   0.0   0.0   0.0   0.0   0.0  1.0\n",
       " -1.0  -1.0   1.0   0.0   0.0   0.0   0.0  1.0\n",
       " -1.0  -1.0  -1.0   1.0   0.0   0.0   0.0  1.0\n",
       " -1.0  -1.0  -1.0  -1.0   1.0   0.0   0.0  1.0\n",
       " -1.0  -1.0  -1.0  -1.0  -1.0   1.0   0.0  1.0\n",
       " -1.0  -1.0  -1.0  -1.0  -1.0  -1.0   1.0  1.0\n",
       " -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  1.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m=8\n",
    "A = tril(2*I - ones(m,m))\n",
    "A[:,end] = ones(m)\n",
    "display(A)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aabda139",
   "metadata": {},
   "source": [
    "Let's compute $A = LU$ and compare $||U||$ and $||A||$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c055f548",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8×8 Matrix{Float64}:\n",
       "  1.0   0.0   0.0   0.0   0.0   0.0   0.0  0.0\n",
       " -1.0   1.0   0.0   0.0   0.0   0.0   0.0  0.0\n",
       " -1.0  -1.0   1.0   0.0   0.0   0.0   0.0  0.0\n",
       " -1.0  -1.0  -1.0   1.0   0.0   0.0   0.0  0.0\n",
       " -1.0  -1.0  -1.0  -1.0   1.0   0.0   0.0  0.0\n",
       " -1.0  -1.0  -1.0  -1.0  -1.0   1.0   0.0  0.0\n",
       " -1.0  -1.0  -1.0  -1.0  -1.0  -1.0   1.0  0.0\n",
       " -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  1.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "8×8 Matrix{Float64}:\n",
       " 1.0  0.0  0.0  0.0  0.0  0.0  0.0    1.0\n",
       " 0.0  1.0  0.0  0.0  0.0  0.0  0.0    2.0\n",
       " 0.0  0.0  1.0  0.0  0.0  0.0  0.0    4.0\n",
       " 0.0  0.0  0.0  1.0  0.0  0.0  0.0    8.0\n",
       " 0.0  0.0  0.0  0.0  1.0  0.0  0.0   16.0\n",
       " 0.0  0.0  0.0  0.0  0.0  1.0  0.0   32.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0  1.0   64.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0  0.0  128.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "F = lu(A)\n",
    "display(F.L)\n",
    "display(F.U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "642a2b3b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.557438524302"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "147.82421993705904"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(norm(A))\n",
    "display(norm(F.U))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "914617dc",
   "metadata": {},
   "source": [
    "Because $||U|| / ||A||$ grows exponentially with $m$, the backward error $||A - \\tilde L \\tilde U||$ may be exponentially large! Solving $Ax = b$ with Gaussian elimination sounds like a bad idea.\n",
    "\n",
    "But what actually happens on a computer?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c8473ff4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(norm(A - F.L*F.U))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d55ebb2b",
   "metadata": {},
   "source": [
    "It looks like the backward error for the LU factorization is actually perfect. Why?\n",
    "\n",
    "Can you modify the problem to bring the expected exponentially bad errors to life? (See Trefethen, problem 22.2.)\n",
    "\n",
    "Note that we DO lose accuracy rapidly (as dimension $m$ increases) when using the LU factors to solve Ax = b. Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d0fae337",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.1181705310112556e-15"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = rand(m)\n",
    "y = F.L \\ b\n",
    "x = F.U \\ y\n",
    "norm(A*x - b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7d74597",
   "metadata": {},
   "source": [
    "In practice, it is extremely rare to encounter matrices for which $||U|| >> ||A||$. We can get a sense of this by examining the ratio $\\rho = ||U||/||A||$ for various random matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "95e9b192",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.23189077910945"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_iter = 10000\n",
    "ρ = zeros(n_iter)\n",
    "for i = 1:n_iter\n",
    "    A = randn(m,m)\n",
    "    F = lu(A)\n",
    "    ρ[i] = norm(F.U) / norm(A)\n",
    "end\n",
    "maximum(ρ)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.3",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}