14_LA_iterative.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": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Iterative Methods\n",
    "\n",
    "In this lecture we will consider a number of classical and more modern methods for solving sparse linear systems like those we found from our consideration of boundary value problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Motivation:  Direct Methods\n",
    "\n",
    "Gaussian elimination or other direct methods for solving linear systems of equations but are there other, better methods for solving these systems when there is structure that can be exploited?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this lecture we will be considering an important type of linear system of equations that arise from discretizing the Poisson problem in one-dimension\n",
    "$$\n",
    "    u''(x) = f(x) \\quad x \\in [a,b] \\quad u(a) = \\alpha \\quad u(b) = \\beta.\n",
    "$$\n",
    "How might you discretize the above ODE?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Following our finite difference discussion we can setup the following set of algebraic equations:\n",
    "$$\n",
    "    \\frac{U_{i-1} - 2 U_i + U_{i+1}}{\\Delta x^2} = f(x_i).\n",
    "$$\n",
    "\n",
    "This forms the $i$th row of the matrix.  See if you can write down this matrix in the form of $Ax = b$ ignoring the boundary conditions for now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    \\frac{1}{\\Delta x^2} \\begin{bmatrix}\n",
    "    -2 &  1 &    &    &    \\\\\n",
    "     1 & -2 &  1 &    &    \\\\\n",
    "       &  1 & -2 &  1 &    \\\\\n",
    "       &    &  1 & -2 &  1 \\\\\n",
    "       &    &    &  1 & -2 \\\\\n",
    "    \\end{bmatrix} \\begin{bmatrix}\n",
    "        U_1 \\\\ U_2 \\\\ U_3 \\\\ U_4 \\\\ U_5\n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix}\n",
    "        f(x_1) - \\frac{\\alpha}{\\Delta x^2} \\\\ f(x_2) \\\\ f(x_3) \\\\ f(x_4) \\\\ f(x_5) - \\frac{\\beta}{\\Delta x^2} \\\\\n",
    "    \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Computational Load\n",
    "\n",
    "Now consider using Gaussian elimination on the above matrix.  For good measure let us consider a 3D problem and discretize each dimension with a $m = 100$ leading to $100 \\times 100 \\times 100 = 10^6$ unknowns.\n",
    "\n",
    "Gaussian Elimination - $\\mathcal{O}(N^3)$ operations to solve, $(10^6)^3 = 10^{18}$ operations.\n",
    "\n",
    "Suppose you have a machine that can perform 100 gigaflops (floating point operations per second):\n",
    "$$\n",
    "    \\frac{10^{18}~ [\\text{flop}]}{10^{11}~ [\\text{flop / s}]} = 10^7~\\text{s} \\approx 115~\\text{days}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Memory Load\n",
    "\n",
    "What about memory?\n",
    "\n",
    "We require $N^2$ to store entire array.  In double precision floating point we would require 8-bytes per entry leading to\n",
    "$$\n",
    "    (10^6)^2 ~[\\text{entries}] \\times 8 ~[\\text{bytes / entry}] = 8 \\times 10^{12} ~[\\text{bytes}] = 8 ~[\\text{terabytes}].\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The situation really is not as bad as we are making it out to be as long as we take advantage of the sparse nature of the matrices.  In fact for 1 dimensional problems direct methods can be reduced to $\\mathcal{O}(N)$ in the case for a tridiagonal system.  The situation is not so great for higher-dimensional problems however unless more structure can be leveraged.  Examples of these types of solvers include fast Fourier methods such as fast Poisson solvers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi and Gauss-Seidel\n",
    "\n",
    "The Jacobi and Gauss-Seidel methods are simple approaches to introducing an iterative means for solving the problem $Ax = b$ when $A$ is sparse.  Consider the general equation derived from the Poisson problem\n",
    "$$\n",
    "    \\frac{U_{i-1} - 2 U_i + U_{i+1}}{\\Delta x^2} = f(x_i).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we rearrange this expression to solve for $U_i$ we have\n",
    "$$\n",
    "    U_i = \\frac{1}{2} (U_{i+1} + U_{i-1}) - f(x_i) \\frac{\\Delta x^2}{2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For a direct method we would simultaneously find the values of $U_i$, $U_{i+1}$ and $U_{i-1}$ but instead consider the iterative scheme computes an update to the equation above by using the past iterate (values we already know)\n",
    "$$\n",
    "    U_i^{(k+1)} = \\frac{1}{2} (U_{i+1}^{(k)} + U_{i-1}^{(k)}) - f(x_i) \\frac{\\Delta x^2}{2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since this allows us to evaluate $U_i^{(k + 1)}$ without knowing the values of $U_{i+1}^{(k)} + U_{i-1}^{(k)}$ we directly evaluate this expression!  This process is called **Jacobi iteration**.  It can be shown that for this particular problem Jacobi iteration will converge from any initial guess $U^{(0)}$ although slowly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Advantages\n",
    " - Matrix $A$ is never stored or created\n",
    " - Storage is optimal\n",
    " - $\\mathcal{O}(m)$ are required per iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "Let's try to solve the problem before in the BVP section but use Jacobi iterations to replace the direct solve\n",
    "$$\n",
    "    u_{xx} = e^x, \\quad x \\in [0, 1] \\quad \\text{with} \\quad u(0) = 0, \\text{ and } u(1) = 3\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Problem setup\n",
    "a = 0.0\n",
    "b = 1.0\n",
    "u_a = 0.0\n",
    "u_b = 3.0\n",
    "f = lambda x: numpy.exp(x)\n",
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)\n",
    "\n",
    "# Descretization\n",
    "N = 100\n",
    "x_bc = numpy.linspace(a, b, N + 2)\n",
    "x = x_bc[1:-1]\n",
    "delta_x = (b - a) / (N + 1)\n",
    "\n",
    "# Expected iterations needed\n",
    "iterations_J = int(2.0 * numpy.log(delta_x) / numpy.log(1.0 - 0.5 * numpy.pi**2 * delta_x**2))\n",
    "\n",
    "# Solve system\n",
    "# Initial guess for iterations\n",
    "U_new = numpy.zeros(N + 2)\n",
    "U_new[0] = u_a\n",
    "U_new[-1] = u_b\n",
    "convergence_J = numpy.zeros(iterations_J)\n",
    "for k in range(iterations_J):\n",
    "    U = U_new.copy()\n",
    "    for i in range(1, N + 1):\n",
    "        U_new[i] = 0.5 * (U[i+1] + U[i-1]) - f(x_bc[i]) * delta_x**2 / 2.0\n",
    "\n",
    "    convergence_J[k] = numpy.linalg.norm(u_true(x_bc) - U_new, ord=2)\n",
    "        \n",
    "# Plot result\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x_bc, U, 'o', label=\"Computed\")\n",
    "axes.plot(x_bc, u_true(x_bc), 'k', label=\"True\")\n",
    "axes.set_title(\"Solution to $u_{xx} = e^x$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"u(x)\")\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(range(iterations_J), convergence_J, 'o')\n",
    "axes.set_title(\"Convergence of Jacobi Iterations for Poisson Problem\")\n",
    "axes.set_xlabel(\"Iteration\")\n",
    "axes.set_ylabel(\"$||U^{(k)} - U^{(k-1)}||_2$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A slight modification to the above leads also to the Gauss-Seidel method.  Programmtically it is easy to see the modification but in the iteration above we now will have\n",
    "$$\n",
    "    U_i^{(k+1)} = \\frac{1}{2} (U_{i+1}^{(k)} + U_{i-1}^{(k+1)}) - f(x_i) \\frac{\\Delta x^2}{2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Problem setup\n",
    "a = 0.0\n",
    "b = 1.0\n",
    "u_a = 0.0\n",
    "u_b = 3.0\n",
    "f = lambda x: numpy.exp(x)\n",
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)\n",
    "\n",
    "# Descretization\n",
    "N = 100\n",
    "x_bc = numpy.linspace(a, b, N + 2)\n",
    "x = x_bc[1:-1]\n",
    "delta_x = (b - a) / (N + 1)\n",
    "\n",
    "# Expected iterations needed\n",
    "iterations_GS = int(2.0 * numpy.log(delta_x) / numpy.log(1.0 - numpy.pi**2 * delta_x**2))\n",
    "\n",
    "# Solve system\n",
    "# Initial guess for iterations\n",
    "U = numpy.zeros(N + 2)\n",
    "U[0] = u_a\n",
    "U[-1] = u_b\n",
    "convergence_GS = numpy.zeros(iterations_GS)\n",
    "success = False\n",
    "for k in range(iterations_GS):\n",
    "    for i in range(1, N + 1):\n",
    "        U[i] = 0.5 * (U[i+1] + U[i-1]) - f(x_bc[i]) * delta_x**2 / 2.0\n",
    "\n",
    "    convergence_GS[k] = numpy.linalg.norm(u_true(x_bc) - U, ord=2)\n",
    "\n",
    "# Plot result\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x_bc, U, 'o', label=\"Computed\")\n",
    "axes.plot(x_bc, u_true(x_bc), 'k', label=\"True\")\n",
    "axes.set_title(\"Solution to $u_{xx} = e^x$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"u(x)\")\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(range(iterations_GS), convergence_GS, 'o')\n",
    "axes.set_title(\"Convergence of Gauss-Seidel Iterations for Poisson Problem\")\n",
    "axes.set_xlabel(\"Iteration\")\n",
    "axes.set_ylabel(\"$||U^{(k)} - U^{(k-1)}||_2$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Matrix Splitting Methods\n",
    "\n",
    "One way to view Jacobi and Gauss-Seidel is as a splitting of the matrix $A$ so that\n",
    "$$\n",
    "    A = M - N.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Then the system $A U = b$ can be viewed as\n",
    "$$\n",
    "    M U - N U = b \\Rightarrow MU = NU + b.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Viewing this instead as an iteration we have then\n",
    "$$\n",
    "    M U^{(k+1)} = N U^{(k)} + b.\n",
    "$$\n",
    "The goal then would be to pick $M$ and $N$ such that $M$ contains as much of $A$ as possible while remaining easier to solve than the original system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The resulting update for each of these then becomes\n",
    "$$\\begin{aligned}\n",
    "    U^{(k+1)} &= M^{-1} N U^{(k)} &+& M^{-1} b \\\\\n",
    "    &= G U^{(k)} &+&  c\n",
    "\\end{aligned}$$\n",
    "where $G$ is called the **iteration matrix** and $c = M^{-1} b$.  We also want\n",
    "$$\n",
    "    u = G u + c\n",
    "$$\n",
    "where $u$ is the true solution of the original $A u = b$, in other words $u$ is the fixed point of the iteration.  Is this fixed point stable though?  If the spectral radius $\\rho(G) < 1$ we can show that in fact the iteration is stable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For Jacobi the splitting is \n",
    "$$\n",
    "    M = -\\frac{2}{\\Delta x^2} I, \\quad \\text{and} \\quad N = \\frac{1}{\\Delta x^2} \\begin{bmatrix}\n",
    "        0 & 1 & \\\\\n",
    "        1 & 0 & 1 \\\\\n",
    "          & \\ddots & \\ddots & \\ddots \\\\\n",
    "          & & 1 & 0 & 1 \\\\\n",
    "          & &   & 1 & 0\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "(sticking to the Poisson problem).  $M$ is now a diagonal matrix and easy to solve.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For Gauss-Seidel we have\n",
    "$$\n",
    "    M = \\frac{1}{\\Delta x^2} \\begin{bmatrix}\n",
    "        -2 &  & \\\\\n",
    "         1 & -2 &  \\\\\n",
    "           & \\ddots & \\ddots \\\\\n",
    "           & & 1 & -2 & \\\\\n",
    "           & &   & 1 & -2\n",
    "    \\end{bmatrix} \\quad \\text{and} \\quad \n",
    "    N = \\frac{1}{\\Delta x^2} \\begin{bmatrix}\n",
    "         0 & 1 & \\\\\n",
    "          & 0 & 1 \\\\\n",
    "          & & \\ddots & \\ddots \\\\\n",
    "           & &  & 0 & 1\\\\\n",
    "           & &   &  & 0\n",
    "    \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Stopping Criteria\n",
    "\n",
    "How many iterations should we perform?  Let $E^{(k)}$ represent the error present at step $k$.  If we want to reduce the error at the first step $E^{(0)}$ to order $\\epsilon$ then we have\n",
    "$$\n",
    "    ||E^{(k)}|| \\approx \\epsilon ||E^{(0)}||.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Under suitable assumption we can bound the error in the 2-norm as\n",
    "$$\n",
    "    ||E^{(k)}|| \\leq \\rho(G)^k ||E^{(0)}||.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Moving back to our estimate of the number of iterations we can combine our two expressions involving the error $E$ by taking $\\Delta x \\rightarrow 0$ which allows us to write\n",
    "$$\n",
    "    k \\approx \\frac{\\log \\epsilon}{\\log \\rho(G)}\n",
    "$$\n",
    "taking into account error convergence.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Picking $\\epsilon$ is a bit tricky but one natural criteria to use would be $\\epsilon = \\mathcal{O}(\\Delta x^2)$ since our original discretization was 2nd-order accurate.  This leads to\n",
    "$$\n",
    "    k = \\frac{2 \\cdot \\log \\Delta x}{\\log \\rho}.\n",
    "$$\n",
    "This also allows us to estimate the total number of operations that need to be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For Jacobi we have the spectral radius of $G$ as\n",
    "$$\n",
    "    \\rho_J \\approx 1 - \\frac{1}{2} \\pi^2 \\Delta x^2.\n",
    "$$\n",
    "so that\n",
    "$$\n",
    "    k = \\mathcal{O}(m^2 \\log m) \\quad \\text{as} \\quad m \\rightarrow \\infty\n",
    "$$\n",
    "where $m$ here is now the number of points used.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Combining this with the previous operation count per iteration we find that Jacobi would lead to $\\mathcal{O}(m^3 \\log m)$ work which is not very promising.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For two dimensions we have $\\mathcal{O}(m^4 \\log m)$ so even compared to Gaussian elimination this approach is not ideal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What about Gauss-Seidel?  Here the spectral radius is approximately\n",
    "$$\n",
    "    \\rho_{GS} \\approx 1 - \\pi^2 \\Delta x^2\n",
    "$$\n",
    "so that\n",
    "$$\n",
    "    k = \\frac{2 \\cdot \\log \\Delta x}{\\log (1 - \\pi^2 \\Delta x^2)}\n",
    "$$\n",
    "which still does not lead to any advantage over direct solvers.  It does show that Gauss-Seidel does actually converge faster due to the factor of 2 difference between $\\rho_J$ and $\\rho_{GS}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Successive Overrelaxation (SOR)\n",
    "\n",
    "Well that's a bit dissapointing isn't it?  These iterative schemes do not seem to be worth much but it turns out we can do better with a slight modification to Gauss-Seidel.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If you look at Gauss-Seidel iteration it turns out it moves $U$ in the correct direction to $u$ but is very conservative in the amount.  If instead we do\n",
    "$$\\begin{aligned}\n",
    "    U^{GS}_i &= \\frac{1}{2} \\left(U^{(k+1)}_{i-1} + U^{(k)}_{i+1} - \\Delta x^2 f_i\\right) \\\\\n",
    "    U^{(k+1)}_i &= U_i^{(k)} + \\omega \\left( U_i^{GS} - U_i^{(k)}\\right )\n",
    "\\end{aligned}$$\n",
    "where we get to pick $\\omega$ we can do much better.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    U^{GS}_i &= \\frac{1}{2} \\left(U^{(k+1)}_{i-1} + U^{(k)}_{i+1} - \\Delta x^2 f_i\\right) \\\\\n",
    "    U^{(k+1)}_i &= U_i^{(k)} + \\omega \\left( U_i^{GS} - U_i^{(k)}\\right )\n",
    "\\end{aligned}$$\n",
    "\n",
    "If $\\omega = 1$ then we get back Gauss-Seidel.  \n",
    "\n",
    "If $\\omega < 1$ we move even less and converges even more slowly (although is sometimes used for multigrid under the name underrelaxation).  \n",
    "\n",
    "If $\\omega > 1$ then we move further than Gauss-Seidel suggests and any method where $\\omega > 1$ is known as **successive overrelaxation** (SOR)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can write this as a matrix splitting method as well.  We can combine the two-step formula above to find\n",
    "$$\n",
    "    U^{(k+1)}_i = \\frac{\\omega}{2} \\left( U^{(k+1)}_{i-1} + U^{(k)}_{i+1} - \\Delta x^2 f_i \\right ) + (1 - \\omega) U_i^{(k)}\n",
    "$$\n",
    "corresponding to a matrix splitting of\n",
    "$$\n",
    "    M = \\frac{1}{\\omega} (D - \\omega L) \\quad \\text{and} \\quad  N = \\frac{1}{\\omega} ((1-\\omega) D + \\omega U)\n",
    "$$\n",
    "where $D$ is the diagonal of the matrix $A$, and $L$ and $U$ are the lower and upper triangular parts without the diagonal of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It can be shown that picking an $\\omega$ such that $0 < \\omega < 2$ the SOR method converges.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It turns out we can also find an optimal $\\omega$ for a wide class of problems.  For Poisson problems in any number of space dimensions for instance it can be shown that SOR method converges optimaly if\n",
    "$$\n",
    "    \\omega_{opt} = \\frac{2}{1 + \\sin(\\pi \\Delta x)} \\approx 2 - 2 \\pi \\Delta x.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What about the number of iterations?  We can follow the same tactic as before with the spectral radius of $G_{SOR}$ now\n",
    "$$\n",
    "    \\rho = \\omega_{opt} - 1 \\approx 1 - 2 \\pi \\Delta x.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This leads to an iteration count of\n",
    "$$\n",
    "    k = \\mathcal{O}(m \\log m)\n",
    "$$\n",
    "an order of magnitude better than Gauss-Seidel alone!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Problem setup\n",
    "a = 0.0\n",
    "b = 1.0\n",
    "u_a = 0.0\n",
    "u_b = 3.0\n",
    "f = lambda x: numpy.exp(x)\n",
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)\n",
    "\n",
    "# Descretization\n",
    "N = 100\n",
    "x_bc = numpy.linspace(a, b, N + 2)\n",
    "x = x_bc[1:-1]\n",
    "delta_x = (b - a) / (N + 1)\n",
    "\n",
    "# SOR parameter\n",
    "omega = 2.0 / (1.0 + numpy.sin(numpy.pi * delta_x))\n",
    "\n",
    "# Expected iterations needed\n",
    "iterations_SOR = int(2.0 * numpy.log(delta_x) / numpy.log(omega - 1.0))\n",
    "\n",
    "# Solve system\n",
    "# Initial guess for iterations\n",
    "U = numpy.zeros(N + 2)\n",
    "U[0] = u_a\n",
    "U[-1] = u_b\n",
    "convergence_SOR = numpy.zeros(iterations_SOR)\n",
    "for k in range(iterations_SOR):\n",
    "    for i in range(1, N + 1):\n",
    "        U_gs = 0.5 * (U[i-1] + U[i+1] - delta_x**2 * f(x_bc[i]))\n",
    "        U[i] += omega * (U_gs - U[i])\n",
    "\n",
    "    convergence_SOR[k] = numpy.linalg.norm(u_true(x_bc) - U, ord=2)\n",
    "        \n",
    "# Plot result\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x_bc, U, 'o', label=\"Computed\")\n",
    "axes.plot(x_bc, u_true(x_bc), 'k', label=\"True\")\n",
    "axes.set_title(\"Solution to $u_{xx} = e^x$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"u(x)\")\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(range(iterations_SOR), convergence_SOR, 'o')\n",
    "axes.set_title(\"Convergence of SOR Iterations for Poisson Problem\")\n",
    "axes.set_xlabel(\"Iteration\")\n",
    "axes.set_ylabel(\"$||U^{(k)} - U^{(k-1)}||_2$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Plotting all the convergence rates\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(range(iterations_SOR), convergence_J[:iterations_SOR], 'r', label=\"Jacobi\")\n",
    "axes.semilogy(range(iterations_SOR), convergence_GS[:iterations_SOR], 'b', label=\"Gauss-Seidel\")\n",
    "axes.semilogy(range(iterations_SOR), convergence_SOR, 'k', label=\"SOR\")\n",
    "axes.legend(loc=3)\n",
    "axes.set_title(\"Comparison of Convergence Rates\")\n",
    "axes.set_xlabel(\"Iteration\")\n",
    "axes.set_ylabel(\"$||U^{(k)} - U^{(k-1)}||_2$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Descent Methods\n",
    "\n",
    "One special case of matrices are amenable to another powerful way to iterate to the solution.  A matrix is said to be **symmetric positive definite** (SPD) if \n",
    "$$\n",
    "    x^T A x > 0 \\quad \\forall \\quad x \\neq 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As an example check that the matrix\n",
    "$$\n",
    "    A = \\begin{bmatrix}\n",
    "        2 &-1 &0 &0 \\\\\n",
    "        -1 & 2 & -1 & 0 \\\\\n",
    "        0 & -1 & 2 & -1 \\\\\n",
    "        0 & 0 & -1 & 2\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "is symmetric positive definite."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now define a function $\\phi: \\mathbb R^m \\rightarrow \\mathbb R$ such that\n",
    "$$\n",
    "    \\phi(u) = \\frac{1}{2} u^T A u - u^T f.\n",
    "$$\n",
    "This is a quadratic function in the variables $u_i$ and in the case where $m = 2$ forms a parabolic bowl.  Since this is a quadratic function there is a unique minimum, $u^\\ast$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def phi(X, Y, A, f):\n",
    "    return 0.5 * (A[0, 0] * X**2 + A[0, 1] * X * Y + A[1, 0] * X * Y + A[1, 1] * Y**2) - X * f[0] - Y * f[1]\n",
    "\n",
    "x = numpy.linspace(-15, 15, 100)\n",
    "y = numpy.linspace(-15, 15, 100)\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "f = numpy.array([1.0, 2.0])\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 3)\n",
    "axes = fig.add_subplot(1, 3, 1, aspect='equal')\n",
    "A = numpy.identity(2)\n",
    "axes.contour(X, Y, phi(X, Y, A, f), 25)\n",
    "\n",
    "axes = fig.add_subplot(1, 3, 2, aspect='equal')\n",
    "A = numpy.array([[2, -1], [-1, 2]])\n",
    "axes.contour(X, Y, phi(X, Y, A, f), 25)\n",
    "\n",
    "axes = fig.add_subplot(1, 3, 3, aspect='equal')\n",
    "A = numpy.array([[2, -1.8], [-1.7, 2]])\n",
    "axes.contour(X, Y, phi(X, Y, A, f), 25)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Lets see how approaching the problem like this helps us:\n",
    "\n",
    "$$\n",
    "    \\phi(u) = \\frac{1}{2} u^T A u - u^T f.\n",
    "$$\n",
    "For the $m = 2$ case write the function $\\phi(u)$ out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    \\phi(u) = \\frac{1}{2} (A_{11} u_1^2 + A_{12} u_1 u_2 + A_{21} u_1 u_2 + A_{22} u^2_2) - u_1 f_1 - u_2 f_2\n",
    "$$\n",
    "\n",
    "What property of the matrix $A$ simplifies the expression above?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Symmetry!  This implies that $A_{21} = A_{12}$ and the expression above simplifies to \n",
    "$$\n",
    "   \\phi(u) = \\frac{1}{2} (A_{11} u_1^2 + 2 A_{12} u_1 u_2 + A_{22} u^2_2) - u_1 f_1 - u_2 f_2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now write two expressions that when evaluated at $u^\\ast$ are identically 0 that express that $u^\\ast$ minimizes $\\phi(u)$.\n",
    "\n",
    "$$\n",
    "   \\phi(u) = \\frac{1}{2} (A_{11} u_1^2 + 2 A_{12} u_1 u_2 + A_{22} u^2_2) - u_1 f_1 - u_2 f_2\n",
    "$$\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since $u^\\ast$ minimizes $\\phi(u)$ we know that the first derivatives should be zero at the minimum:\n",
    "   $$\\begin{aligned}\n",
    "       \\frac{\\partial \\phi}{\\partial u_1} &= A_{11} u_1 + A_{12} u_2 - f_1 = 0 \\\\\n",
    "       \\frac{\\partial \\phi}{\\partial u_1} &= A_{21} u_1 + A_{22} u_2 - f_2 = 0\n",
    "   \\end{aligned}$$\n",
    "   Note that these equations can be rewritten as\n",
    "   $$\n",
    "       A u = f.\n",
    "   $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The conclusion here then is that finding the minimum of $\\phi$ is equivalent to solving the system $A u = f$!  This is a common type of reformulation for many problems where it may be easier to treat a given equation as a minimization problem rather than directly solve it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that this is not quite the matrix that we have been using for our Poisson problem so far which is actually symmetric negative definite although these same methods work as well.  In this case we actually want to find the maximum of $\\phi$ instead, other than that everything is the same."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Also note that if $A$ is indefinite then the eigenvalues of $A$ will change sign and instead of a stable minimum or maximum we have a saddle point which are much more difficult to handle (GMRES can for instance)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Method of Steepest Descent\n",
    "\n",
    "So now we turn to finding the $u^\\ast$ that minimizes the function $\\phi(u)$.  The simplest approach to this is called the **method of steepest descent** which finds the direction of the largest gradient of $\\phi(u)$ and goes in that direction.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Mathematically we then have\n",
    "$$\n",
    "    u^{(k+1)} = u^{(k)} - \\alpha^{(k)} \\nabla \\phi(u^{(k)})\n",
    "$$\n",
    "where $\\alpha^{(k)}$ will be the step size chosen in the direction we want to go.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can find $\\alpha$ by\n",
    "$$\n",
    "    \\alpha^{(k)} = \\min_{\\alpha \\in \\mathbb R} \\phi\\left(u^{(k)} - \\alpha \\nabla \\phi(u^{(k)}\\right),\n",
    "$$\n",
    "i.e. the $\\alpha$ that takes us just far enough so that if we went any further $\\phi$ would increase.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This implies that $\\alpha^{(k)} \\geq 0$ and $\\alpha^{(k)} = 0$ only if we are at the minimum of $\\phi$.  We can compute the gradient of $\\phi$ as\n",
    "$$\n",
    "    \\nabla \\phi(u^{(k)}) = A u^{(k)} - f \\equiv -r^{(k)}\n",
    "$$\n",
    "where $r^{(k)}$ is the residual defined as\n",
    "$$\n",
    "    r^{(k)} = f - A u^{(k)}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Looking back at the definition of $\\alpha^{(k)}$ then leads to the conclusion that the $\\alpha$ that would minimize the expression would be the one that satisfies\n",
    "$$\n",
    "    \\frac{\\text{d} \\phi(\\alpha)}{\\text{d} \\alpha} = 0. \n",
    "$$\n",
    "\n",
    "To find this note that\n",
    "$$\n",
    "    \\phi(u + \\alpha r) = \\left(\\frac{1}{2} u^T A u - u^T f \\right) + \\alpha(r^T A u - r^T f) + \\frac{1}{2} \\alpha^2 r^T A r\n",
    "$$\n",
    "so that the derivative becomes\n",
    "$$\n",
    "    \\frac{\\text{d} \\phi(\\alpha)}{\\text{d} \\alpha} = r^T A u - r^T f + \\alpha r^T A r\n",
    "$$\n",
    "\n",
    "Setting this to zero than leads to\n",
    "$$\n",
    "    \\alpha = \\frac{r^T r}{r^T A r}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Problem setup\n",
    "a = 0.0\n",
    "b = 1.0\n",
    "u_a = 0.0\n",
    "u_b = 3.0\n",
    "f = lambda x: numpy.exp(x)\n",
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)\n",
    "\n",
    "# Descretization\n",
    "N = 50\n",
    "x_bc = numpy.linspace(a, b, N + 2)\n",
    "x = x_bc[1:-1]\n",
    "delta_x = (b - a) / (N + 1)\n",
    "\n",
    "# Construct matrix A\n",
    "A = numpy.zeros((N, N))\n",
    "diagonal = numpy.ones(N) / delta_x**2\n",
    "A += numpy.diag(diagonal * 2.0, 0)\n",
    "A += numpy.diag(-diagonal[:-1], 1)\n",
    "A += numpy.diag(-diagonal[:-1], -1)\n",
    "\n",
    "# Construct right hand side\n",
    "b = -f(x)\n",
    "b[0] += u_a / delta_x**2\n",
    "b[-1] += u_b / delta_x**2\n",
    "\n",
    "# Algorithm parameters\n",
    "MAX_ITERATIONS = 10000\n",
    "tolerance = 1e-3\n",
    "\n",
    "# Solve system\n",
    "U = numpy.empty(N)\n",
    "convergence_SD = numpy.zeros(MAX_ITERATIONS)\n",
    "success = False\n",
    "for k in range(MAX_ITERATIONS):\n",
    "    r = b - numpy.dot(A, U)\n",
    "    if numpy.linalg.norm(r, ord=2) < tolerance:\n",
    "        success = True\n",
    "        break\n",
    "        \n",
    "    alpha = numpy.dot(r, r) / numpy.dot(r, numpy.dot(A, r))\n",
    "    U = U + alpha * r\n",
    "\n",
    "    convergence_SD[k] = numpy.linalg.norm(u_true(x) - U, ord=2)\n",
    "        \n",
    "if not success:\n",
    "    print(\"Iteration failed to converge!\")\n",
    "    print(convergence_SD[-1])\n",
    "else:\n",
    "    # Plot result\n",
    "    fig = plt.figure()\n",
    "    axes = fig.add_subplot(1, 1, 1)\n",
    "    axes.plot(x, U, 'o', label=\"Computed\")\n",
    "    axes.plot(x_bc, u_true(x_bc), 'k', label=\"True\")\n",
    "    axes.set_title(\"Solution to $u_{xx} = e^x$\")\n",
    "    axes.set_xlabel(\"x\")\n",
    "    axes.set_ylabel(\"u(x)\")\n",
    "\n",
    "    fig = plt.figure()\n",
    "    axes = fig.add_subplot(1, 1, 1)\n",
    "    axes.semilogy(range(k), convergence_SD[:k], 'o')\n",
    "    axes.set_title(\"Convergence of Steepest Descent for Poisson Problem\")\n",
    "    axes.set_xlabel(\"Iteration\")\n",
    "    axes.set_ylabel(\"$||U^{(k)} - U^{(k-1)}||_2$\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Convergence of Steepest Descent\n",
    "\n",
    "What controls the convergence of steepest descent?  It turns out that the shape of the parabolic bowl formed by $\\phi$ is the major factor determining the convergence of steepest descent.  For example, if $A$ is a scalar multiple of the identity than these ellipses are actually circles and steepest descent converges in $m$ steps.  If $A$ does not lead to circles, the convergence is based on the ratio between the semi-major and semi-minor axis of the resulting ellipses $m$ dimensional ellipses.  This is controlled by the smallest and largest eigenvalues of the matrix $A$ hence why steepest descent grows increasingly difficult as $m$ increases for the Poisson problem.  Note that this also relates to the condition number of the matrix in the $\\ell_2$ norm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Ellipses](./images/ellipses.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def steepest_descent(A, U, b, axes):\n",
    "    MAX_ITERATIONS = 10000\n",
    "    tolerance = 1e-3\n",
    "    success = False\n",
    "    for k in range(MAX_ITERATIONS):\n",
    "        axes.text(U[0] + 0.1, U[1] + 0.1, str(k), fontsize=12)\n",
    "        axes.plot(U[0], U[1], 'ro')\n",
    "        r = b - numpy.dot(A, U)\n",
    "        if numpy.linalg.norm(r, ord=2) < tolerance:\n",
    "            success = True\n",
    "            break\n",
    "\n",
    "        alpha = numpy.dot(r, r) / numpy.dot(r, numpy.dot(A, r))\n",
    "        U = U + alpha * r\n",
    "        \n",
    "    if success:\n",
    "        return k\n",
    "    else:\n",
    "        raise Exception(\"Iteration did not converge.\")\n",
    "\n",
    "phi = lambda X, Y, A: 0.5 * (A[0, 0] * X**2 + A[0, 1] * X * Y + A[1, 0] * X * Y + A[1, 1] * Y**2) - X * f[0] - Y * f[1]\n",
    "\n",
    "x = numpy.linspace(-15, 15, 100)\n",
    "y = numpy.linspace(-15, 15, 100)\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "# fig.set_figwidth(fig.get_figwidth() * 3)\n",
    "# axes = fig.add_subplot(1, 3, 1, aspect='equal')\n",
    "axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "A = numpy.identity(2)\n",
    "f = numpy.array([0.0, 0.0])\n",
    "axes.contour(X, Y, phi(X, Y, A), 25)\n",
    "print(\"Iteration count: %s\" % steepest_descent(A, numpy.array([-10.0, -12.0]), f, axes))\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "# axes = fig.add_subplot(1, 3, 2, aspect='equal')\n",
    "axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "A = numpy.array([[2, -1], [-1, 2]])\n",
    "f = numpy.array([1.0, 2.0])\n",
    "axes.contour(X, Y, phi(X, Y, A), 25)\n",
    "print(\"Iteration count: %s\" % steepest_descent(A, numpy.array([-10.0, -12.0]), f, axes))\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "# axes = fig.add_subplot(1, 3, 3, aspect='equal')\n",
    "A = numpy.array([[2, -1.8], [-1.7, 2]])\n",
    "f = numpy.array([1.0, 2.0])\n",
    "axes.contour(X, Y, phi(X, Y, A), 25)\n",
    "print(\"Iteration count: %s\" % steepest_descent(A, numpy.array([-10.0, -12.0]), f, axes))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### $A$-Conjugate Search Directions and Conjugate Gradient\n",
    "\n",
    "An alternative to steepest descent is to choose a slightly different direction to descend down.  Generalizing our step from above let the iterative scheme be  \n",
    "$$\n",
    "    u^{(k+1)} = u^{(k)} - \\alpha^{(k)} p^{(k)}\n",
    "$$\n",
    "and as before we want to pick an $\\alpha$ such that\n",
    "$$\n",
    "    \\min_{\\alpha} \\phi(u^{(k)} - \\alpha p^{(k)})\n",
    "$$\n",
    "leading again to the choice of $\\alpha$ of\n",
    "$$\n",
    "    \\alpha^{(k)} = \\frac{(p^{(k)})^T p^{(k)}}{(p^{(k)})^T A p^{(k)}}\n",
    "$$\n",
    "accept now we are also allowed to pick the search direction $p^{(k)}$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Ways to choose $p^{(k)}$:\n",
    " - One bad choice for $p$ would be orthogonal to $r$ since this would then be tangent to the level set (ellipse) of $\\phi(u^{(k)})$ and would cause it to only increase so we want to make sure that $p^T r \\neq 0$ (the inner-product).\n",
    " - We also want to still move downwards so require that $\\phi(u^{(k+1)}) < \\phi(u^{(k)})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We know that $r^{(k)}$ is not always the best direction to go in but what might be better?  We could head directly for the minimum but how do we do that without first knowing $u^\\ast$?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It turns out when $m=2$ we can do this from any initial guess $u^{(0)}$ and initial direction $p^{(k)}$ we will arrive at the minimum in 2 steps if we choose the second search direction dependent on\n",
    "$$\n",
    "    (p^{(1)})^T A p^{(0)} = 0.\n",
    "$$\n",
    "In general if two vectors satisfy this property they are said to be $A$-conjugate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that if $A = I$ then these two vectors would be orthogonal to each other and in this sense $A$-conjugacy is a natural extension from orthogonality and the simple case from before where the ellipses are circles to the case where we can have very distorted ellipses.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In fact the vector $p^{(0)}$ is tangent to the level set that $u^{(1)}$ lies on and therefore choosing $p^{(1)}$ so that it is $A$-conjugate to $p^{(0)}$ then always heads to the center of the ellipse.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In other words, once we know a tangent to one of the ellipses we can always choose a direction that minimizes in one of the dimensions of the search space.  Choosing the $p^{(k)}$ iteratively this way forms the basis of the **conjugate gradient** method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Ellipses CG](./images/ellipses_CG.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to generalize beyond $m = 2$ consider the $m=3$ case.  As stated before we are now in a three-dimensional space where the level-sets are concentric ellipsoids.  Taking a slice through this space will lead to an ellipse on the slice.\n",
    "\n",
    "1. Start with an initial guess $u^{(0)}$ and choose a search direction $p^{(0)}$.\n",
    "1. Minimize $\\phi(u)$ in the direction $u^{(0)} + \\alpha p^{(0)}$ resulting in the choice\n",
    "$$\n",
    "    \\alpha^{(0)} = \\frac{(p^{(0)})^T p^{(0)}}{(p^{(0)})^T A p^{(0)}}\n",
    "$$\n",
    "as we saw before.  Now set $u^{(1)} = u^{(0)} + \\alpha^{(0)} p^{(0)}$.\n",
    "1. Choose $p^{(1)}$ to be $A$-conjugate to $p^{(0)}$.  In this case there are an infinite set of vectors that are possible that satisfy $(p^{(1)})^T A p^{(0)} = 0$.  Beyond requiring $p^{(1)}$ to be $A$-conjugate we also want it to be linearly-independent to $p^{(0)}$.\n",
    "1. Again choose an $\\alpha^{(1)}$ that minimizes the residual (again tangent to the level sets of $\\phi$) in the direction $p^{(1)}$ and repeat the process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Problem setup\n",
    "a = 0.0\n",
    "b = 1.0\n",
    "u_a = 0.0\n",
    "u_b = 3.0\n",
    "f = lambda x: numpy.exp(x)\n",
    "u_true = lambda x: (4.0 - numpy.exp(1.0)) * x - 1.0 + numpy.exp(x)\n",
    "\n",
    "# Descretization\n",
    "N = 50\n",
    "x_bc = numpy.linspace(a, b, N + 2)\n",
    "x = x_bc[1:-1]\n",
    "delta_x = (b - a) / (N + 1)\n",
    "\n",
    "# Construct matrix A\n",
    "A = numpy.zeros((N, N))\n",
    "diagonal = numpy.ones(N) / delta_x**2\n",
    "A += numpy.diag(diagonal * 2.0, 0)\n",
    "A += numpy.diag(-diagonal[:-1], 1)\n",
    "A += numpy.diag(-diagonal[:-1], -1)\n",
    "\n",
    "# Construct right hand side\n",
    "b = -f(x)\n",
    "b[0] += u_a / delta_x**2\n",
    "b[-1] += u_b / delta_x**2\n",
    "\n",
    "# Algorithm parameters\n",
    "MAX_ITERATIONS = 2 * N\n",
    "tolerance = delta_x**2\n",
    "\n",
    "# Solve system\n",
    "U = numpy.zeros(N)\n",
    "convergence_CG = numpy.zeros(MAX_ITERATIONS)\n",
    "success = False\n",
    "r = numpy.dot(A, U) - b\n",
    "p = -r\n",
    "r_prod_prev = numpy.dot(r, r)\n",
    "for k in range(MAX_ITERATIONS):\n",
    "    w = numpy.dot(A, p)\n",
    "    alpha = r_prod_prev / numpy.dot(p, w)\n",
    "    U += alpha * p\n",
    "    r += alpha * w\n",
    "    beta = numpy.dot(r, r) / r_prod_prev\n",
    "    r_prod_prev = numpy.dot(r, r)\n",
    "    if numpy.linalg.norm(r, ord=2) < tolerance:\n",
    "        success = True\n",
    "        break\n",
    "    p = beta * p - r\n",
    "    convergence_CG[k] = numpy.linalg.norm(u_true(x) - U, ord=2)\n",
    "\n",
    "if not success:\n",
    "    print(\"Iteration failed to converge!\")\n",
    "else:\n",
    "    # Plot result\n",
    "    fig = plt.figure()\n",
    "    axes = fig.add_subplot(1, 1, 1)\n",
    "    axes.plot(x, U, 'o', label=\"Computed\")\n",
    "    axes.plot(x_bc, u_true(x_bc), 'k', label=\"True\")\n",
    "    axes.set_title(\"Solution to $u_{xx} = e^x$\")\n",
    "    axes.set_xlabel(\"x\")\n",
    "    axes.set_ylabel(\"u(x)\")\n",
    "\n",
    "    fig = plt.figure()\n",
    "    axes = fig.add_subplot(1, 1, 1)\n",
    "    axes.semilogy(convergence_CG[:k+2], 'o')\n",
    "    axes.set_title(\"Convergence of Conjugate Gradient for Poisson Problem\")\n",
    "    axes.set_xlabel(\"Iteration\")\n",
    "    axes.set_ylabel(\"$||U^{(k)} - U^{(k-1)}||_2$\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def CG(A, U, b, axes):\n",
    "    MAX_ITERATIONS = 10000\n",
    "    tolerance = 1e-3\n",
    "    success = False        \n",
    "    r = numpy.dot(A, U) - b\n",
    "    p = -r\n",
    "    r_prod_prev = numpy.dot(r, r)\n",
    "    axes.text(U[0] + 0.1, U[1] + 0.1, str(0), fontsize=12)\n",
    "    axes.plot(U[0], U[1], 'ro')\n",
    "    for k in range(MAX_ITERATIONS):\n",
    "        w = numpy.dot(A, p)\n",
    "        alpha = r_prod_prev / numpy.dot(p, w)\n",
    "        U += alpha * p\n",
    "        axes.text(U[0] + 0.1, U[1] + 0.1, str(k + 1), fontsize=12)\n",
    "        axes.plot(U[0], U[1], 'ro')\n",
    "        r += alpha * w\n",
    "        beta = numpy.dot(r, r) / r_prod_prev\n",
    "        r_prod_prev = numpy.dot(r, r)\n",
    "        if numpy.linalg.norm(r, ord=2) < tolerance:\n",
    "            success = True\n",
    "            break\n",
    "        p = beta * p - r\n",
    "        \n",
    "    if success:\n",
    "        return k\n",
    "    else:\n",
    "        raise Exception(\"Iteration did not converge.\")\n",
    "\n",
    "phi = lambda X, Y, A: 0.5 * (A[0, 0] * X**2 + A[0, 1] * X * Y + A[1, 0] * X * Y + A[1, 1] * Y**2) - X * f[0] - Y * f[1]\n",
    "\n",
    "x = numpy.linspace(-15, 15, 100)\n",
    "y = numpy.linspace(-15, 15, 100)\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "# fig.set_figwidth(fig.get_figwidth() * 3)\n",
    "# axes = fig.add_subplot(1, 3, 1, aspect='equal')\n",
    "axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "A = numpy.identity(2)\n",
    "f = numpy.array([0.0, 0.0])\n",
    "axes.contour(X, Y, phi(X, Y, A), 25)\n",
    "print(\"Iteration count: %s\" % CG(A, numpy.array([-10.0, -12.0]), f, axes))\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "# axes = fig.add_subplot(1, 3, 2, aspect='equal')\n",
    "axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "A = numpy.array([[2, -1], [-1, 2]])\n",
    "f = numpy.array([1.0, 2.0])\n",
    "axes.contour(X, Y, phi(X, Y, A), 25)\n",
    "print(\"Iteration count: %s\" % CG(A, numpy.array([-10.0, -12.0]), f, axes))\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "# axes = fig.add_subplot(1, 3, 3, aspect='equal')\n",
    "A = numpy.array([[2, -1.8], [-1.7, 2]])\n",
    "f = numpy.array([1.0, 2.0])\n",
    "axes.contour(X, Y, phi(X, Y, A), 25)\n",
    "print(\"Iteration count: %s\" % CG(A, numpy.array([-10.0, -12.0]), f, axes))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Preconditioning\n",
    "\n",
    "One way to get around the difficulties with these types of methods due to the distortion of the ellipses (and consequently the conditioning of the matrix) is to precondition the matrix.  The basic idea is that we take our original problem $A u = f$ and instead solve\n",
    "$$\n",
    "    M^{-1} A u = M^{-1} f.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that since we need to find the inverse of $M$, this matrix should be nice.  A couple of illustrative examples may help to illustrate why this might be a good idea:\n",
    "\n",
    " - If $M = A$ then we essentially have solved our problem already although that does not help us much\n",
    " - If $M = \\text{diag}(A)$, then $M^{-1}$ is easily computed and it turns out for some problems this can decrease the condition number of $M^{-1} A$ significantly.  Note though that this is not actually helpful in the case of the Poisson problem.\n",
    " - If $M$ is based on another iterative method used on $A$, for instance Gauss-Seidel, these can be effective general use preconditioners for many problems.\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So the next question then becomes how to choose a preconditioner.  This is usually very problem specific and a number of papers suggest strategies for particular problems."
   ]
  }
 ],
 "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.8.8"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}