05_root_finding_optimization.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>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Root Finding and Optimization\n",
    "\n",
    "**GOAL:** Find where $f(x) = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  Future Time Annuity\n",
    "\n",
    "When can I retire?\n",
    "\n",
    "$$ A = \\frac{P}{(r / m)} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ] $$\n",
    "\n",
    "$A$ total value after $n$ years\n",
    "\n",
    "$P$ is payment amount per compounding period\n",
    "\n",
    "$m$ number of compounding periods per year\n",
    "\n",
    "$r$ annual interest rate\n",
    "\n",
    "$n$ number of years to retirement\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If I want to retire in 20 years what does the annual interest rate $r$ need to be?\n",
    "\n",
    "Set $P = \\frac{\\$18,000}{12} = \\$1500, \\quad m=12, \\quad n=20$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def total_value(P, m, r, n):\n",
    "    \"\"\"Total value of portfolio given parameters\n",
    "    \n",
    "    Based on following formula:\n",
    "    \n",
    "    A = \\frac{P}{(r / m)} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n}\n",
    "                - 1 \\right ] \n",
    "    \n",
    "    :Input:\n",
    "     - *P* (float) - Payment amount per compounding period\n",
    "     - *m* (int) - number of compounding periods per year\n",
    "     - *r* (float) - annual interest rate\n",
    "     - *n* (float) - number of years to retirement\n",
    "     \n",
    "     :Returns:\n",
    "     (float) - total value of portfolio\n",
    "     \n",
    "    \"\"\"\n",
    "    \n",
    "    return P / (r / float(m)) * ( (1.0 + r / float(m))**(float(m) * n)\n",
    "                                 - 1.0)\n",
    "\n",
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "    \n",
    "r = numpy.linspace(0.05, 0.1, 100)\n",
    "goal = 1e6\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, total_value(P, m, r, n))\n",
    "axes.plot(r, numpy.ones(r.shape) * goal, 'r--')\n",
    "axes.set_xlabel(\"r (interest rate)\")\n",
    "axes.set_ylabel(\"A (total value)\")\n",
    "axes.set_title(\"When can I retire?\")\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fixed Point Iteration\n",
    "\n",
    "How do we go about solving this?\n",
    "\n",
    "Could try to solve at least partially for $r$:\n",
    "\n",
    "$$ A = \\frac{P}{(r / m)} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ] ~~~~ \\Rightarrow ~~~~~$$\n",
    "\n",
    "$$ r = \\frac{P \\cdot m}{A} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ] ~~~~ \\Rightarrow ~~~~~$$\n",
    "\n",
    "$$ r = g(r)$$\n",
    "or \n",
    "$$ g(r) - r = 0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def g(P, m, r, n, A):\n",
    "    \"\"\"Reformulated minimization problem\n",
    "    \n",
    "    Based on following formula:\n",
    "    \n",
    "    g(r) = \\frac{P \\cdot m}{A} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ]\n",
    "    \n",
    "    :Input:\n",
    "     - *P* (float) - Payment amount per compounding period\n",
    "     - *m* (int) - number of compounding periods per year\n",
    "     - *r* (float) - annual interest rate\n",
    "     - *n* (float) - number of years to retirement\n",
    "     - *A* (float) - total value after $n$ years\n",
    "     \n",
    "     :Returns:\n",
    "     (float) - value of g(r)\n",
    "     \n",
    "    \"\"\"\n",
    "    \n",
    "    return P * m / A * ( (1.0 + r / float(m))**(float(m) * n)\n",
    "                                 - 1.0)\n",
    "\n",
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "    \n",
    "r = numpy.linspace(0.00, 0.1, 100)\n",
    "goal = 1e6\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, g(P, m, r, n, goal))\n",
    "axes.plot(r, r, 'r--')\n",
    "axes.set_xlabel(\"r (interest rate)\")\n",
    "axes.set_ylabel(\"$g(r)$\")\n",
    "axes.set_title(\"When can I retire?\")\n",
    "axes.set_ylim([0, 0.12])\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Guess at $r_0$ and check to see what direction we need to go...\n",
    "\n",
    "1. $r_0 = 0.0800, \\quad g(r_0) - r_0 = -0.009317550125425428$\n",
    "1. $r_1 = 0.0850, \\quad g(r_1) - r_1 = -0.00505763375972$\n",
    "1. $r_2 = 0.0875, \\quad g(r_2) - r_2 = -0.00257275331014$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A bit tedious, we can also make this algorithmic:\n",
    "```python\n",
    "r_values = numpy.linspace(0.08, 0.09, 10)\n",
    "for r in r_values:\n",
    "    print(\"r = \", r, \"g(r) =\", g(P, m, r, n, goal))\n",
    "    print(\"Difference = \", numpy.abs(g(P, m, r, n, goal) - r))\n",
    "    r = g(P, m, r, n, goal)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "r_values = numpy.linspace(0.08, 0.09, 11)\n",
    "for r in r_values:\n",
    "    print(\"r = \", r, \"g(r) =\", g(P, m, r, n, goal))\n",
    "    print(\"Difference = \", numpy.abs(g(P, m, r, n, goal) - r))\n",
    "    r = g(P, m, r, n, goal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 2:\n",
    "\n",
    "Let $f(x) = x - e^{-x}$, solve $f(x) = 0$\n",
    "\n",
    "Equivalent to $x = e^{-x}$ or $x = g(x)$ where $g(x) = e^{-x}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.2, 1.0, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, numpy.exp(-x), 'r')\n",
    "axes.plot(x, x, 'b')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"f(x)\")\n",
    "\n",
    "x = 0.4\n",
    "for steps in range(3):\n",
    "    print(\"x = \", x, \"Residual = \", numpy.abs(numpy.exp(-x) - x))\n",
    "    x = numpy.exp(-x)\n",
    "    axes.plot(x, numpy.exp(-x),'kx')\n",
    "    axes.text(x, numpy.exp(-x), steps+1, fontsize=\"15\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 3:\n",
    "\n",
    "Let $f(x) = \\ln x + x$ and solve $f(x) = 0$ or $x = -\\ln x$.\n",
    "\n",
    "Note that this problem is equivalent to $x = e^{-x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.1, 1.0, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, -numpy.log(x), 'r')\n",
    "axes.plot(x, x, 'b')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"f(x)\")\n",
    "axes.set_ylim([0, 1.5])\n",
    "\n",
    "x = 0.5\n",
    "for steps in range(3):\n",
    "    print(\"x = \", x, \"Residual = \", numpy.abs(numpy.log(x) + x))\n",
    "    x = -numpy.log(x)\n",
    "    axes.plot(x, -numpy.log(x),'kx')\n",
    "    axes.text(x, -numpy.log(x), steps+1, fontsize=\"15\")\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "These are equivalent problems!  Something is awry..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Analysis of Fixed Point Iteration\n",
    "\n",
    "Existence and uniqueness of fixed point problems\n",
    "\n",
    "*Existence:*\n",
    "\n",
    "Assume $g \\in C[a, b]$, if the range of the mapping $y = g(x)$ satisfies $y \\in [a, b] \\quad \\forall \\quad x \\in [a, b]$ then $g$ has a fixed point in $[a, b]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.0, 1.0, 100)\n",
    "\n",
    "# Plot function and intercept\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, numpy.exp(-x), 'r')\n",
    "axes.plot(x, x, 'b')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"f(x)\")\n",
    "\n",
    "# Plot domain and range\n",
    "axes.plot(numpy.ones(x.shape) * 0.4, x, '--k')\n",
    "axes.plot(numpy.ones(x.shape) * 0.8, x, '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * numpy.exp(-0.4), '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * numpy.exp(-0.8), '--k')\n",
    "\n",
    "axes.set_xlim((0.0, 1.0))\n",
    "axes.set_ylim((0.0, 1.0))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.1, 1.0, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, -numpy.log(x), 'r')\n",
    "axes.plot(x, x, 'b')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"f(x)\")\n",
    "axes.set_xlim([0.1, 1.0])\n",
    "axes.set_ylim([0.1, 1.0])\n",
    "\n",
    "# Plot domain and range\n",
    "axes.plot(numpy.ones(x.shape) * 0.4, x, '--k')\n",
    "axes.plot(numpy.ones(x.shape) * 0.8, x, '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * -numpy.log(0.4), '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * -numpy.log(0.8), '--k')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "r = numpy.linspace(0.06, 0.1, 100)\n",
    "goal = 1e6\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, g(P, m, r, n, goal))\n",
    "axes.plot(r, r, 'r--')\n",
    "axes.set_xlabel(\"r\")\n",
    "axes.set_ylabel(\"$g(r)$\")\n",
    "axes.set_xlim([0.06, 0.1])\n",
    "axes.set_ylim([g(P, m, 0.06, n, goal), g(P, m, 0.1, n, goal)])\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "\n",
    "axes.plot([0.08, 0.08], [g(P, m, 0.06, n, goal), g(P, m, 0.1, n, goal)], '--k')\n",
    "axes.plot([0.09, 0.09], [g(P, m, 0.06, n, goal), g(P, m, 0.1, n, goal)], '--k')\n",
    "axes.plot(r, numpy.ones(r.shape) * g(P, m, 0.08, n, goal), '--k')\n",
    "axes.plot(r, numpy.ones(r.shape) * g(P, m, 0.09, n, goal), '--k')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "*Uniqueness:*\n",
    "\n",
    "Additionally, suppose $g'(x)$ is defined on $x \\in [a, b]$ and $\\exists K < 1$ such that\n",
    "\n",
    "$$\n",
    "    |g'(x)| \\leq K < 1 \\quad \\forall \\quad x \\in (a,b)\n",
    "$$\n",
    "\n",
    "then $g$ has a unique fixed point $P \\in [a,b]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.4, 0.8, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, numpy.abs(-numpy.exp(-x)), 'r')\n",
    "axes.plot(x, numpy.ones(x.shape), 'k--')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"f(x)\")\n",
    "axes.set_ylim((0.0, 1.1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "*Asymptotic convergence*: Behavior of fixed point iterations\n",
    "\n",
    "$$x_{k+1} = g(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Assume that $\\exists ~ x^\\ast$ s.t. $x^\\ast = g(x^\\ast)$ (i.e. $x^\\ast$ is the fixed point), then define\n",
    "\n",
    "$$\n",
    "    x_k = x^\\ast + e_k \\quad \\quad x_{k+1} = x^\\ast + e_{k+1}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "    x^\\ast + e_{k+1} = g(x^\\ast + e_k)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Taylor expand the function $g$ about $x^\\ast$:\n",
    "\n",
    "$$\n",
    "    g(x) = g(x^\\ast) + g'(x^\\ast) (x - x^\\ast) + \\frac{g''(x^\\ast)}{2!} (x - x^\\ast)^2 + \\mathcal{O}((x - x^\\ast)^3)\n",
    "$$\n",
    "\n",
    "Evaluate this series at $x_k = x^\\ast + e_k$ to find\n",
    "\n",
    "$$\n",
    "    g(x^\\ast + e_k) = g(x^\\ast) + g'(x^\\ast) e_k + \\frac{g''(x^\\ast) e_k^2}{2} + \\mathcal{O}(e_k^3)\n",
    "$$\n",
    "\n",
    "therefore from our definition from before that $x^\\ast + e_{k+1} = g(x^\\ast + e_k)$ we have\n",
    "\n",
    "$$\n",
    "    x^\\ast + e_{k+1} = g(x^\\ast) + g'(x^\\ast) e_k + \\frac{g''(x^\\ast) e_k^2}{2} + \\mathcal{O}(e_k^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that because $x^* = g(x^*)$ these terms cancel leaving\n",
    "\n",
    "$$e_{k+1} = g'(x^*) e_k + \\frac{g''(x^*) e_k^2}{2}$$\n",
    "\n",
    "So if $|g'(x^*)| \\leq K < 1$ we can conclude that\n",
    "\n",
    "$$|e_{k+1}| = K |e_k|$$\n",
    "\n",
    "which shows convergence.  Also note that $K$ is related to $|g'(x^*)|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convergence of iterative schemes\n",
    "\n",
    "Given any iterative scheme where\n",
    "\n",
    "$$|e_{k+1}| = C |e_k|^n$$\n",
    "\n",
    "If $C < 1$ and:\n",
    " - $n=1$ then the scheme is **linearly convergent**\n",
    " - $n=2$ then the scheme is **quadratically convergent**\n",
    " - $n > 1$ the scheme can also be called **superlinearly convergent**\n",
    "\n",
    "If $C > 1$ then the scheme is **divergent**\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Examples Revisited\n",
    "$g(x) = e^{-x}$ with $x^* \\approx 0.56$\n",
    " \n",
    "   $$|g'(x^*)| = |-e^{-x^*}| \\approx 0.56$$\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$g(x) = - \\ln x \\quad \\text{with} \\quad x^* \\approx 0.56$\n",
    "\n",
    "   $$|g'(x^*)| = \\frac{1}{|x^*|} \\approx 1.79$$\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    r = g(r) = \\frac{P \\cdot m}{A} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import sympy\n",
    "r, P, m, A, n = sympy.symbols('r P m A n')\n",
    "g = P * m / A * ((1 + r /m)**(m * n) - 1)\n",
    "g_prime = g.diff(r)\n",
    "r_star = 0.08985602484084668\n",
    "print(\"g'(r) = \", g_prime)\n",
    "print(\"g'(r*) = \", g_prime.subs({P: 1500.0, m: 12, n:20, A: 1e6, r: r_star}))\n",
    "\n",
    "f = sympy.lambdify(r, g_prime.subs({P: 1500.0, m: 12, n:20, A: 1e6}))\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "r = numpy.linspace(-0.01, 0.1, 100)\n",
    "axes.plot(r, f(r))\n",
    "axes.plot(r, numpy.ones(r.shape), 'k--')\n",
    "axes.plot(r_star, f(r_star), 'ro')\n",
    "axes.plot(0.0, f(0.0), 'ro')\n",
    "axes.set_xlim((-0.01, 0.1))\n",
    "axes.set_xlabel(\"$r$\")\n",
    "axes.set_ylabel(\"$g'(r)$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Better ways for root-finding/optimization\n",
    "\n",
    "If $x^*$ is a fixed point of $g(x)$ then $x^*$ is also a *root* of $f(x^*) = g(x^*) - x^*$ s.t. $f(x^*) = 0$.\n",
    "\n",
    "For instance:\n",
    "\n",
    "$$f(r) = r - \\frac{m P}{A} \\left [ \\left (1 + \\frac{r}{m} \\right)^{m n} - 1 \\right ] =0 $$\n",
    "\n",
    "or\n",
    "\n",
    "$$f(r) = A - \\frac{m P}{r} \\left [ \\left (1 + \\frac{r}{m} \\right)^{m n} - 1 \\right ] =0 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Classical Methods\n",
    " - Bisection (linear convergence)\n",
    " - Newton's Method (quadratic convergence)\n",
    " - Secant Method (super-linear)\n",
    " \n",
    "## Combined Methods\n",
    " - RootSafe (Newton + Bisection)\n",
    " - Brent's Method (Secant + Bisection)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bracketing and Bisection\n",
    "\n",
    "A **bracket** is an interval $[a,b]$ that contains exactly one zero or minima/maxima of interest.  \n",
    "\n",
    "In the case of a zero the bracket should satisfy \n",
    "$$\n",
    "    \\text{sign}(f(a)) \\neq \\text{sign}(f(b)).\n",
    "$$\n",
    "\n",
    "In the case of minima or maxima we need \n",
    "$$\n",
    "    \\text{sign}(f'(a)) \\neq \\text{sign}(f'(b))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Theorem**:  \n",
    "\n",
    "Let\n",
    "$$\n",
    "    f(x) \\in C[a,b] \\quad \\text{and} \\quad \\text{sign}(f(a)) \\neq \\text{sign}(f(b))\n",
    "$$\n",
    "\n",
    "then there exists a number \n",
    "$$\n",
    "    c \\in (a,b) \\quad \\text{s.t.} \\quad f(c) = 0.\n",
    "$$\n",
    "(proof uses intermediate value theorem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.1, 100)\n",
    "f = lambda r, A, m, P, n: A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r, A, m, P, n), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "axes.set_xlabel(\"r (%)\")\n",
    "axes.set_ylabel(\"f(r)\")\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "\n",
    "a = 0.075\n",
    "b = 0.095\n",
    "axes.plot(a, f(a, A, m, P, n), 'ko')\n",
    "axes.plot([a, a], [0.0, f(a, A, m, P, n)], 'k--')\n",
    "axes.plot(b, f(b, A, m, P, n), 'ko')\n",
    "axes.plot([b, b], [f(b, A, m, P, n), 0.0], 'k--')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Basic bracketing algorithms shrink the bracket while ensuring that the root/extrema remains within the bracket.\n",
    "\n",
    "What ways could we \"shrink\" the bracket so that the end points converge to the root/extrema?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Bisection Algorithm\n",
    "\n",
    "Given a bracket $[a,b]$ and a function $f(x)$ - \n",
    "1. Initialize with bracket\n",
    "2. Iterate\n",
    "   1. Cut bracket in half and check to see where the zero is\n",
    "   2. Set bracket to new bracket based on what direction we went"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "f = lambda r, A=A, m=m, P=P, n=n: A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "\n",
    "# Initialize bracket\n",
    "a = 0.07\n",
    "b = 0.10\n",
    "\n",
    "# Setup figure to plot convergence\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r, A, m, P, n), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "axes.set_xlabel(\"r (%)\")\n",
    "axes.set_ylabel(\"f(r)\")\n",
    "# axes.set_xlim([0.085, 0.091])\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.plot(a, f(a, A, m, P, n), 'ko')\n",
    "axes.plot([a, a], [0.0, f(a, A, m, P, n)], 'k--')\n",
    "axes.plot(b, f(b, A, m, P, n), 'ko')\n",
    "axes.plot([b, b], [f(b, A, m, P, n), 0.0], 'k--')\n",
    "\n",
    "# Algorithm parameters\n",
    "TOLERANCE = 1e-4\n",
    "MAX_STEPS = 5\n",
    "\n",
    "# Initialize loop\n",
    "delta_x = b - a\n",
    "c = a + delta_x / 2.0\n",
    "f_a = f(a)\n",
    "f_b = f(b)\n",
    "f_c = f(c)\n",
    "\n",
    "# Loop until we reach the TOLERANCE or we take MAX_STEPS\n",
    "for step in range(1, MAX_STEPS + 1):\n",
    "    \n",
    "    # Plot iteration\n",
    "    axes.plot(c, f_c,'kx')\n",
    "    axes.text(c, f_c, str(step + 1), fontsize=\"15\")\n",
    "\n",
    "    # Check tolerance - Could also check the size of delta_x\n",
    "    # We check this first as we have already initialized the values\n",
    "    # in c and f_c\n",
    "    if numpy.abs(f_c) < TOLERANCE:\n",
    "        break\n",
    "\n",
    "    if numpy.sign(f_a) != numpy.sign(f_c):\n",
    "        b = c\n",
    "        f_b = f_c\n",
    "    else:\n",
    "        a = c\n",
    "        f_a = f_c\n",
    "    delta_x = b - a\n",
    "    c = a + delta_x / 2.0\n",
    "    f_c = f(c)\n",
    "        \n",
    "if step == MAX_STEPS:\n",
    "    print(\"Reached maximum number of steps!\")\n",
    "else:\n",
    "    print(\"Success!\")\n",
    "    print(\"  x* = %s\" % c)\n",
    "    print(\"  f(x*) = %s\" % f(c))\n",
    "    print(\"  number of steps = %s\" % step)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Convergence of Bisection\n",
    "\n",
    "Generally have\n",
    "$$\n",
    "    |e_{k+1}| = C |e_k|^n\n",
    "$$\n",
    "where we need $C < 1$ and $n > 0$.\n",
    "\n",
    "Letting $\\Delta x_k$ be the width of the $k$th bracket we can then estimate the error with\n",
    "$$\n",
    "    e_k \\approx \\Delta x_k\n",
    "$$\n",
    "and therefore\n",
    "$$\n",
    "    e_{k+1} \\approx \\frac{1}{2} \\Delta x_k.\n",
    "$$\n",
    "Due to the relationship then between $x_k$ and $e_k$ we then know\n",
    "$$\n",
    "    |e_{k+1}| = \\frac{1}{2} |e_k|\n",
    "$$\n",
    "so therefore the method is linearly convergent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Newton's Method (Newton-Raphson)\n",
    " - Given a bracket, bisection is guaranteed to converge linearly to a root\n",
    " - However bisection uses almost no information about $f(x)$ beyond its sign at a point\n",
    " \n",
    "**Basic Idea**: Given $f(x)$ and $f'(x)$ use a linear approximation to $f(x)$ \"locally\" and use the x-intercept of the resulting line to predict where $x^*$ might be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given current location $x_k$, we have $f(x_k)$ and $f'(x_k)$ and form a line through the point $(x_k, f(x_k))$:\n",
    "\n",
    "Form equation for the line:\n",
    "\n",
    "$$y = f'(x_k) x + b$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Solve for the y-intercept value $b$\n",
    "\n",
    "$$f(x_k) = f'(x_k) x_k + b$$\n",
    "\n",
    "$$b = f(x_k) - f'(x_k) x_k$$\n",
    "\n",
    "and simplify.\n",
    "\n",
    "$$y = f'(x_k) x + f(x_k) - f'(x_k) x_k$$\n",
    "\n",
    "$$y = f'(x_k) (x - x_k) + f(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now find the intersection of our line and the x-axis (i.e. when $y = 0$) and use the resulting value of $x$ to set $x_{k+1}$ \n",
    "\n",
    "$$\n",
    "    0 = f'(x_k) (x_{k+1}-x_k) + f(x_k)\n",
    "$$\n",
    "\n",
    "$$\n",
    "    x_{k+1} = x_k-\\frac{f(x_k)}{f'(x_k)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "An alternative method of derivation for Newton-Raphson (and more in line with our methods) uses Taylor series.  Expand the function $f(x)$ in a Taylor series about the current Newton-Raphson iteration $x_k$:\n",
    "\n",
    "$$\n",
    "    f(x) = f(x_k) + f'(x_k) (x - x_k) + \\frac{f''(x_k)}{2!} (x - x_k)^2 + \\mathcal{O}((x-x_k)^3)\n",
    "$$\n",
    "\n",
    "Let $\\delta_k$ be the update to the $x_{k+1}$ iteration such that\n",
    "$$\n",
    "    x_{k+1} = x_k + \\Delta x_k\n",
    "$$\n",
    "and evaluate our expression for $f(x)$ at $x_{k+1}$:\n",
    "$$\n",
    "    f(x_{k+1}) = f(x_k) + f'(x_k) \\Delta x_k + \\frac{f''(x_k)}{2!} \\Delta x_k^2 + \\mathcal{O}(\\Delta x_k^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now assume that $x_{k+1} = x^\\ast$, if this is the case the above simplifies to\n",
    "$$\n",
    "    0 = f(x_k) + f'(x_k) \\Delta x_k + \\frac{f''(x_k)}{2!} \\Delta x_k^2 + \\mathcal{O}(\\Delta x_k^3)\n",
    "$$\n",
    "and dropping the higher order terms leads to\n",
    "$$\n",
    "    \\Delta x_k = - \\frac{f(x_k)}{f'(x_k)}\n",
    "$$\n",
    "assuming that $f \\in \\mathbb R$ leading to the update\n",
    "$$\n",
    "    x_{k+1} = x_k - \\frac{f(x_k)}{f'(x_k)}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "f = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "        A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "f_prime = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "            -P*m*n*(1.0 + r/m)**(m*n)/(r*(1.0 + r/m)) \\\n",
    "                + P*m*((1.0 + r/m)**(m*n) - 1.0)/r**2\n",
    "\n",
    "# Initial guess\n",
    "x_k = 0.06\n",
    "\n",
    "# Setup figure to plot convergence\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "\n",
    "# Plot x_k point\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "axes.plot(x_k, f(x_k), 'ko')\n",
    "axes.text(x_k, -5e4, \"$x_k$\", fontsize=16)\n",
    "axes.plot(x_k, 0.0, 'xk')\n",
    "axes.text(x_k, f(x_k) + 2e4, \"$f(x_k)$\", fontsize=16)\n",
    "axes.plot(r, f_prime(x_k) * (r - x_k) + f(x_k), 'k')\n",
    "\n",
    "# Plot x_{k+1} point\n",
    "x_k = x_k - f(x_k) / f_prime(x_k)\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "axes.plot(x_k, f(x_k), 'ko')\n",
    "axes.text(x_k, 1e4, \"$x_{k+1}$\", fontsize=16)\n",
    "axes.plot(x_k, 0.0, 'xk')\n",
    "axes.text(0.0873, f(x_k) - 2e4, \"$f(x_{k+1})$\", fontsize=16)\n",
    "\n",
    "axes.set_xlabel(\"r (%)\")\n",
    "axes.set_ylabel(\"f(r)\")\n",
    "axes.set_title(\"Newton-Raphson Steps\")\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What does the alogrithm look like for Newton-Raphson?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Algorithm\n",
    "\n",
    "1. Initialize $x_k$\n",
    "1. Begin loop\n",
    "  1. Compute $f(x_k)$ and $f'(x_k)$\n",
    "  1. Use these to compute new $x_{k+1}$\n",
    "  1. Check stopping criteria"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "f = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "        A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "f_prime = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "            -P*m*n*(1.0 + r/m)**(m*n)/(r*(1.0 + r/m)) \\\n",
    "                + P*m*((1.0 + r/m)**(m*n) - 1.0)/r**2\n",
    "\n",
    "# Algorithm parameters\n",
    "MAX_STEPS = 200\n",
    "TOLERANCE = 1e-4\n",
    "        \n",
    "# Initial guess\n",
    "x_k = 0.06\n",
    "\n",
    "# Setup figure to plot convergence\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "\n",
    "for n in range(1, MAX_STEPS + 1):\n",
    "    axes.plot(x_k, f(x_k),'kx')\n",
    "    axes.text(x_k, f(x_k), str(n), fontsize=\"15\")\n",
    "    x_k = x_k - f(x_k) / f_prime(x_k)\n",
    "    if numpy.abs(f(x_k)) < TOLERANCE:\n",
    "        break\n",
    "        \n",
    "if n == MAX_STEPS:\n",
    "    print(\"Reached maximum number of steps!\")\n",
    "else:\n",
    "    print(\"Success!\")\n",
    "    print(\"  x* = %s\" % x_k)\n",
    "    print(\"  f(x*) = %s\" % f(x_k))\n",
    "    print(\"  number of steps = %s\" % n)\n",
    "\n",
    "axes.set_xlabel(\"r (%)\")\n",
    "axes.set_ylabel(\"f(r)\")\n",
    "axes.set_title(\"Newton-Raphson Steps\")\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:\n",
    "\n",
    "$$f(x) = x - e^{-x}$$\n",
    "\n",
    "$$f'(x) = 1 + e^{-x}$$\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{f(x_k)}{f'(x_k)} = x_k - \\frac{x_k - e^{-x_k}}{1 + e^{-x_k}}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Asymptotic Convergence of Newton's Method\n",
    "\n",
    "For a simple root (non-multiplicative) - Let $g(x) = x - \\frac{f(x)}{f'(x)}$, then\n",
    "\n",
    "$$x_{k+1} = g(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Definitions of errors and iteration:\n",
    "\n",
    "$$x_{k+1} = x^* + e_{k+1} \\quad \\quad x_k = x^* + e_k$$\n",
    "\n",
    "General Taylor expansion:\n",
    "\n",
    "$$\n",
    "    x^* + e_{k+1} = g(x^* + e_k) = g(x^*) + g'(x^*) e_k + \\frac{g''(x^*) e_k^2}{2!} + \\mathcal{O}(e_k^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that as before $x^*$ and $g(x^*)$ cancel:\n",
    "\n",
    "$$e_{k+1} = g'(x^*) e_k + \\frac{g''(x^*) e_k^2}{2!} + \\ldots$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What about $g'(x^*)$ though?  \n",
    "\n",
    "$$\\begin{aligned}\n",
    "    g(x) &= x - \\frac{f(x)}{f'(x)} \\\\\n",
    "    g'(x) & = 1 - \\frac{f'(x)}{f'(x)} + \\frac{f(x) f''(x)}{(f'(x))^2} = \\frac{f(x) f''(x)}{(f'(x))^2}\n",
    "\\end{aligned}$$\n",
    "\n",
    "which evaluated at $x = x^*$ becomes\n",
    "\n",
    "$$\n",
    "    g'(x^*) = \\frac{f(x^*)f''(x^*)}{f'(x^*)^2} = 0\n",
    "$$\n",
    "\n",
    "since $f(x^\\ast) = 0$ by definition (assuming $f''(x^\\ast)$ and $f'(x^\\ast)$ are appropriately behaved)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Back to our expansion we have again\n",
    "\n",
    "$$\n",
    "    e_{k+1} = g'(x^*) e_k + \\frac{g''(x^*) e_k^2}{2!} + \\ldots\n",
    "$$\n",
    "\n",
    "which simplifies to \n",
    "\n",
    "$$\n",
    "    e_{k+1} = \\frac{g''(x^*) e_k^2}{2!} + \\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    e_{k+1} = \\frac{g''(x^*) e_k^2}{2!} + \\ldots\n",
    "$$\n",
    "leads to \n",
    "$$\n",
    "    |e_{k+1}| = \\left | \\frac{g''(x^*)}{2!} \\right | |e_k|^2\n",
    "$$\n",
    "\n",
    "Newton's method is therefore quadratically convergent where the the constant is controlled by the second derivative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For a multiple root (e.g. $f(x) = (x-1)^2$) the case is not particularly rosey unfortunately.  Why might this be?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:\n",
    "$f(x) = \\sin (2 \\pi x)$\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{\\sin (2 \\pi x)}{2 \\pi \\cos (2 \\pi x)}= x_k - \\frac{1}{2 \\pi} \\tan (2 \\pi x)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0, 2, 1000)\n",
    "f = lambda x: numpy.sin(2.0 * numpy.pi * x)\n",
    "f_prime = lambda x: 2.0 * numpy.pi * numpy.cos(2.0 * numpy.pi * x)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, f(x),'b')\n",
    "axes.plot(x, f_prime(x), 'r')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.set_title(\"Comparison of $f(x)$ and $f'(x)$\")\n",
    "axes.set_ylim((-2,2))\n",
    "axes.set_xlim((0,2))\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "\n",
    "x_k = 0.3\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'ko')\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "axes.plot(x, f_prime(x_k) * (x - x_k) + f(x_k), 'k')\n",
    "\n",
    "\n",
    "x_k = x_k - f(x_k) / f_prime(x_k)\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'ko')\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0, 2, 1000)\n",
    "f = lambda x: numpy.sin(2.0 * numpy.pi * x)\n",
    "x_kp = lambda x: x - 1.0 / (2.0 * numpy.pi) * numpy.tan(2.0 * numpy.pi * x)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, f(x),'b')\n",
    "axes.plot(x, x_kp(x), 'r')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.set_title(\"Comparison of $f(x)$ and $f'(x)$\")\n",
    "axes.set_ylim((-2,2))\n",
    "axes.set_xlim((0,2))\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Basins of Attraction\n",
    "\n",
    "Given a point $x_0$ can we determine if Newton-Raphson converges?\n",
    "\n",
    "A *basin of attraction* $X$ for Newton's methods is defined as the set such that $\\forall x \\in X$ Newton iterations converges.  Unfortunately this is far from a trivial thing to determine and even for simple functions can lead to regions that are fractal.  \n",
    "\n",
    "Plotted below are two fairly simple equations which demonstrate the problem:\n",
    "1. $f(x) = x^3 - 1$\n",
    "2. Kepler's equation $\\theta - e \\sin \\theta = M$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: x**3 - 1\n",
    "f_prime = lambda x: 3 * x**2\n",
    "\n",
    "N = 1001\n",
    "x = numpy.linspace(-2, 2, N)\n",
    "X, Y = numpy.meshgrid(x, x)\n",
    "R = X + 1j * Y\n",
    "\n",
    "for i in range(30):\n",
    "    R = R - f(R) / f_prime(R)\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.contour(X, Y, R)\n",
    "axes.set_xlabel(\"Real\")\n",
    "axes.set_ylabel(\"Imaginary\")\n",
    "axes.set_title(\"Basin of Attraction for $f(x) = x^3 - 1$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def f(theta, e=0.083, M=1):\n",
    "    return theta - e * numpy.sin(theta) - M\n",
    "def f_prime(theta, e=0.083):\n",
    "    return 1 - e * numpy.cos(theta)\n",
    "\n",
    "N = 1001\n",
    "x = numpy.linspace(-30.5, -29.5, N)\n",
    "y = numpy.linspace(-17.5, -16.5, N)\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "R = X + 1j * Y\n",
    "\n",
    "for i in range(30):\n",
    "    R = R - f(R) / f_prime(R)\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.contour(X, Y, R)\n",
    "axes.set_xlabel(\"Real\")\n",
    "axes.set_ylabel(\"Imaginary\")\n",
    "axes.set_title(\"Basin of Attraction for $f(x) = x - e \\sin x - M$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Other Issues\n",
    "\n",
    "Need to supply both $f(x)$ and $f'(x)$, could be expensive\n",
    " \n",
    "Example:  FTV equation $f(r) = A - \\frac{m P}{r} \\left[ \\left(1 + \\frac{r}{m} \\right )^{m n} - 1\\right]$\n",
    "\n",
    "Can use symbolic differentiation (`sympy`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Secant Methods\n",
    "\n",
    "Is there a method with the convergence of Newton's method but without the extra derivatives?  What way would you modify Newton's method so that you would not need $f'(x)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given $x_k$ and $x_{k-1}$ represent the derivative as the approximation\n",
    "\n",
    "$$f'(x) \\approx \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}}$$\n",
    "\n",
    "Combining this with the Newton approach leads to\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{f(x_k) (x_k - x_{k-1}) }{f(x_k) - f(x_{k-1})}$$\n",
    "\n",
    "This leads to superlinear convergence and not quite quadratic as the exponent on the convergence is $\\approx 1.7$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Alternative interpretation, fit a line through two points and see where they intersect the x-axis.\n",
    "\n",
    "$$(x_k, f(x_k)) ~~~~~ (x_{k-1}, f(x_{k-1})$$\n",
    "\n",
    "$$y = \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} (x - x_k) + b$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$b = f(x_{k-1}) - \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} (x_{k-1} - x_k)$$\n",
    "\n",
    "$$ y = \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} (x - x_k) + f(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now solve for $x_{k+1}$ which is where the line intersects the x-axies ($y=0$)\n",
    "\n",
    "$$0 = \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} (x_{k+1} - x_k) + f(x_k)$$\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{f(x_k)  (x_k - x_{k-1})}{f(x_k) - f(x_{k-1})}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "f = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "        A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "\n",
    "# Initial guess\n",
    "x_k = 0.07\n",
    "x_km = 0.06\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "\n",
    "axes.plot(x_k, 0.0, 'ko')\n",
    "axes.plot(x_k, f(x_k), 'ko')\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "axes.plot(x_km, 0.0, 'ko')\n",
    "axes.plot(x_km, f(x_km), 'ko')\n",
    "axes.plot([x_km, x_km], [0.0, f(x_km)], 'k--')\n",
    "\n",
    "axes.plot(r, (f(x_k) - f(x_km)) / (x_k - x_km) * (r - x_k) + f(x_k), 'k')\n",
    "x_kp = x_k - (f(x_k) * (x_k - x_km) / (f(x_k) - f(x_km)))\n",
    "axes.plot(x_kp, 0.0, 'ro')\n",
    "axes.plot([x_kp, x_kp], [0.0, f(x_kp)], 'r--')\n",
    "axes.plot(x_kp, f(x_kp), 'ro')\n",
    "\n",
    "axes.set_xlabel(\"r (%)\")\n",
    "axes.set_ylabel(\"f(r)\")\n",
    "axes.set_title(\"Secant Method\")\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What would the algorithm look like for such a method?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Algorithm\n",
    "\n",
    "Given $f(x)$, given bracket $[a,b]$, a `TOLERANCE`, and a `MAX_STEPS` (note we need two points to start).\n",
    "\n",
    "1. Initialize $x_1 = a$, $x_2 = b$, $f_1 = f(x_1)$, and $f_2 = f(x_2)$\n",
    "2. Loop until either `MAX_STEPS` is reached or `TOLERANCE` is achieved\n",
    "   1. Calculate new update $x_{k+1}$ by update formula\n",
    "   2. Check for convergence and break if reached\n",
    "   3. Update parameters $x_1$, $x_2$, $f_1 = f(x_1)$ and $f_2(x_2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "f = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "        A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "f_prime = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "            -P*m*n*(1.0 + r/m)**(m*n)/(r*(1.0 + r/m)) \\\n",
    "                + P*m*((1.0 + r/m)**(m*n) - 1.0)/r**2\n",
    "\n",
    "# Algorithm parameters\n",
    "MAX_STEPS = 50\n",
    "TOLERANCE = 1e-4\n",
    "        \n",
    "# Initial bracket\n",
    "x_k = 0.07\n",
    "x_km = 0.06\n",
    "\n",
    "# Setup figure to plot convergence\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "\n",
    "for n in range(1, MAX_STEPS + 1):\n",
    "    axes.plot(x_k, f(x_k), 'o')\n",
    "    axes.text(x_k + 0.0025, f(x_k), n, fontsize=\"15\")\n",
    "    x_kp = x_k - f(x_k) * (x_k - x_km) / (f(x_k) - f(x_km))\n",
    "    x_km = x_k\n",
    "    x_k = x_kp\n",
    "    print(\"Residual = \", numpy.abs(f(x_k)))\n",
    "    if numpy.abs(f(x_k)) < TOLERANCE:\n",
    "        break\n",
    "        \n",
    "if n == MAX_STEPS:\n",
    "    print(\"Reached maximum number of steps!\")\n",
    "else:\n",
    "    print(\"Success!\")\n",
    "    print(\"  x* = %s\" % x_k)\n",
    "    print(\"  f(x*) = %s\" % f(x_k))\n",
    "    print(\"  number of steps = %s\" % n)\n",
    "\n",
    "axes.set_xlabel(\"r (%)\")\n",
    "axes.set_ylabel(\"f(r)\")\n",
    "axes.set_title(\"Secant Method\")\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Comments\n",
    "\n",
    " - Secant method as shown is equivalent to linear interpolation\n",
    " - Can use higher order interpolation for higher order secant methods\n",
    " - Convergence is not quite quadratic\n",
    " - Not guaranteed to converge\n",
    " - Do not preserve brackets\n",
    " - Almost as good as Newton's method if your initial guess is good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Hybrid Methods\n",
    "\n",
    "Combine attributes of methods with others to make one great algorithm to rule them all (not really)\n",
    "\n",
    "#### Goals\n",
    "1. Robustness:  Given a bracket $[a,b]$, maintain bracket\n",
    "1. Efficiency:  Use superlinear convergent methods when possible\n",
    "\n",
    "#### Options\n",
    " - Methods requiring $f'(x)$\n",
    "   - NewtSafe (RootSafe, Numerical Recipes)\n",
    "   - Newton's Method within a bracket, Bisection otherwise\n",
    " - Methods not requiring $f'(x)$\n",
    "   - Brent's Algorithm (zbrent, Numerical Recipes)\n",
    "     - Combination of bisection, secant and inverse quadratic interpolation\n",
    "   - `scipy.optimize` package"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimization (finding extrema)\n",
    "\n",
    "I want to find the extrema of a function $f(x)$ on a given interval $[a,b]$.\n",
    "\n",
    "A few approaches:\n",
    " - Bracketing Algorithms:  Golden-Section Search (linear)\n",
    " - Interpolation Algorithms: Repeated parabolic interpolation\n",
    " - Hybrid Algorithms "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Bracketing Algorithm (Golden Section Search)\n",
    "\n",
    "Given $f(x) \\in C[x_0,x_3]$ that is convex (concave) over an interval $x \\in [x_0,x_3]$ reduce the interval size until it brackets the minimum (maximum).\n",
    "\n",
    "Note that we no longer have the $x=0$ help we had before so bracketing and doing bisection is a bit more tricky in this case.  In particular choosing your initial bracket is important!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Bracket Picking\n",
    "\n",
    "Say we start with a bracket $[x_0, x_3]$ and pick to new points $x_1 < x_2 \\in [x_0, x_3]$.  We want to pick a new bracket that guarantees that the extrema exists in it.  We then can pick this new bracket with the following rules:\n",
    " - If $f(x_1) < f(x_2)$ then we know the minimum is between $x_0$ and $x_2$.\n",
    " - If $f(x_1) > f(x_2)$ then we know the minimum is between $x_1$ and $x_3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: x**2\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "\n",
    "search_points = [-1.0, -0.5, 0.75, 1.0]\n",
    "axes = fig.add_subplot(2, 2, 1)\n",
    "x = numpy.linspace(search_points[0] - 0.1, search_points[-1] + 0.1, 100)\n",
    "axes.plot(x, f(x), 'b')\n",
    "for (i, point) in enumerate(search_points):\n",
    "    axes.plot(point, f(point),'or')\n",
    "    axes.text(point + 0.05, f(point), str(i))\n",
    "axes.plot(0, 0, 'sk')\n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_title(\"$f(x_1) < f(x_2) \\Rightarrow [x_0, x_2]$\")\n",
    "\n",
    "search_points = [-1.0, -0.75, 0.5, 1.0]\n",
    "axes = fig.add_subplot(2, 2, 2)\n",
    "x = numpy.linspace(search_points[0] - 0.1, search_points[-1] + 0.1, 100)\n",
    "axes.plot(x, f(x), 'b')\n",
    "for (i, point) in enumerate(search_points):\n",
    "    axes.plot(point, f(point),'or')\n",
    "    axes.text(point + 0.05, f(point), str(i))\n",
    "axes.plot(0, 0, 'sk')\n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_title(\"$f(x_1) > f(x_2) \\Rightarrow [x_1, x_3]$\")\n",
    "\n",
    "search_points = [-1.0, 0.25, 0.75, 1.0]\n",
    "axes = fig.add_subplot(2, 2, 3)\n",
    "x = numpy.linspace(search_points[0] - 0.1, search_points[-1] + 0.1, 100)\n",
    "axes.plot(x, f(x), 'b')\n",
    "for (i, point) in enumerate(search_points):\n",
    "    axes.plot(point, f(point),'or')\n",
    "    axes.text(point + 0.05, f(point), str(i))\n",
    "axes.plot(0, 0, 'sk')\n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_title(\"$f(x_1) < f(x_2) \\Rightarrow [x_0, x_2]$\")\n",
    "\n",
    "search_points = [-1.0, -0.75, -0.25, 1.0]\n",
    "axes = fig.add_subplot(2, 2, 4)\n",
    "x = numpy.linspace(search_points[0] - 0.1, search_points[-1] + 0.1, 100)\n",
    "axes.plot(x, f(x), 'b')\n",
    "for (i, point) in enumerate(search_points):\n",
    "    axes.plot(point, f(point),'or')\n",
    "    axes.text(point + 0.05, f(point), str(i))\n",
    "axes.plot(0, 0, 'sk')\n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_title(\"$f(x_1) > f(x_2) \\Rightarrow [x_1, x_3]$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Picking Brackets and Points\n",
    "\n",
    "Again say we have a bracket $[x_0,x_3]$ and suppose we have two new search points $x_1$ and $x_2$ that separates $[x_0,x_3]$ into two new overlapping brackets.\n",
    "\n",
    "Define \n",
    "$$\\begin{aligned}\n",
    "    a &= x_1 - x_0, \\\\\n",
    "    b &= x_3 - x_1,\\\\\n",
    "    c &= x_2 - x_1 \\quad \\text{and} \\\\\n",
    "    d &= x_3 - x_2.\n",
    "\\end{aligned}$$\n",
    "\n",
    "For **Golden Section Search** we require two conditions:\n",
    " - The two new possible brackets are of equal length.  If we pick the left bracket $[x_0, x_2]$ then\n",
    " $$\n",
    "     a+c = b\n",
    " $$\n",
    " and the right bracket $[x_1, x_3]$\n",
    " $$\n",
    "     d + c = b.\n",
    " $$\n",
    " - The distances between subsequent triplets is proportional.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: (x - 0.25)**2 + 0.5\n",
    "phi = (1.0 + numpy.sqrt(5.0)) / 2.0\n",
    "\n",
    "x = [-1.0, None, None, 1.0]\n",
    "x[1] = x[3] - 1.0 / phi * (x[3] - x[0])\n",
    "x[2] = x[0] + 1.0 / phi * (x[3] - x[0])\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "axes = []\n",
    "axes.append(fig.add_subplot(1, 2, 1))\n",
    "axes.append(fig.add_subplot(1, 2, 2))\n",
    "t = numpy.linspace(-2.0, 2.0, 100)\n",
    "for i in range(2):\n",
    "    axes[i].plot(t, f(t), 'k')\n",
    "\n",
    "    # First set of intervals\n",
    "    axes[i].plot([x[0], x[2]], [0.0, 0.0], 'g')\n",
    "    axes[i].plot([x[1], x[3]], [-0.2, -0.2], 'r')\n",
    "    axes[i].plot([x[0], x[0]], [0.0, f(x[0])], 'g--')\n",
    "    axes[i].plot([x[2], x[2]], [0.0, f(x[2])], 'g--')\n",
    "    axes[i].plot([x[1], x[1]], [-0.2, f(x[2])], 'r--')\n",
    "    axes[i].plot([x[3], x[3]], [-0.2, f(x[3])], 'r--')\n",
    "    for (n, point) in enumerate(x):\n",
    "        axes[i].plot(point, f(point), 'ok')\n",
    "        axes[i].text(point, f(point)+0.1, n, fontsize='15')\n",
    "\n",
    "    axes[i].set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "    axes[i].set_ylim((-1.0, 3.0))\n",
    "\n",
    "# Left new interval\n",
    "x_new = [x[0], None, x[1], x[2]]\n",
    "x_new[1] = 1.0 / phi * (x[1] - x[0]) + x[0]\n",
    "\n",
    "axes[0].plot([x_new[0], x_new[2]], [1.5, 1.5], 'b')\n",
    "axes[0].plot([x_new[1], x_new[3]], [1.75, 1.75], 'c')\n",
    "axes[0].plot([x_new[0], x_new[0]], [1.5, f(x_new[0])], 'b--')\n",
    "axes[0].plot([x_new[2], x_new[2]], [1.5, f(x_new[2])], 'b--')\n",
    "axes[0].plot([x_new[1], x_new[1]], [1.75, f(x_new[1])], 'c--')\n",
    "axes[0].plot([x_new[3], x_new[3]], [1.75, f(x_new[3])], 'c--')\n",
    "axes[0].plot(x_new[1], f(x_new[1]), 'ko')\n",
    "axes[0].text(x_new[1] + 0.05, f(x_new[1]) + 0.1, \"*\", fontsize='15')\n",
    "    \n",
    "# Right new interval\n",
    "x_new = [x[1], x[2], None, x[3]]\n",
    "x_new[2] = (x[2] - x[1]) / phi + x[2]\n",
    "\n",
    "axes[1].plot([x_new[0], x_new[2]], [1.25, 1.25], 'b')\n",
    "axes[1].plot([x_new[1], x_new[3]], [1.5, 1.5], 'c')\n",
    "axes[1].plot([x_new[0], x_new[0]], [1.25, f(x_new[0])], 'b--')\n",
    "axes[1].plot([x_new[2], x_new[2]], [1.25, f(x_new[2])], 'b--')\n",
    "axes[1].plot([x_new[1], x_new[1]], [1.5, f(x_new[2])], 'c--')\n",
    "axes[1].plot([x_new[3], x_new[3]], [1.5, f(x_new[3])], 'c--')\n",
    "axes[1].plot(x_new[2], f(x_new[2]), 'ko')\n",
    "axes[1].text(x_new[2] + 0.05, f(x_new[2]) + 0.1, \"*\", fontsize='15')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The first rule implies:\n",
    "$$\\begin{aligned}\n",
    "    a + c &= b \\\\\n",
    "    x_1 - x_0 + x_2 - x_1 &= x_3 - x_1 \\\\\n",
    "    x_2 - x_0 &= x_3 - x_1.\n",
    "\\end{aligned}$$\n",
    "\n",
    "Assume that this allows us to pick $x_2$ (we need to figure out how to choose $x_1$).  We then know\n",
    "$$\n",
    "    x_2 = x_3 - x_1 + x_0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Subsequent proportionality implies that the distances between the 4 points at one iteration is proportional to the next.  Since we have two choices for our new interval we write down many proportionality constraints however let us focus on the two defined by the distances $a$, $b$, and $c$.\n",
    "\n",
    "If $f(x_1) < f(x_2)$ then we choose $(x_0, x_1, x_2)$ as our new triplet meaning\n",
    "$$\n",
    "    \\frac{a}{b} = \\frac{c}{a}\n",
    "$$\n",
    "\n",
    "If $f(x_1) > f(x_2)$ then we choose $(x_1, x_2, x_3)$ as our new triplet meaning\n",
    "$$\n",
    "    \\frac{a}{b} = \\frac{c}{b-c}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using these relations we can solve for the ratio $b / a$ via the following.  Take\n",
    "$$\n",
    "    \\frac{a}{b} = \\frac{c}{a} \\quad \\text{and} \\quad \\frac{a}{b} = \\frac{c}{b-c}\n",
    "$$\n",
    "and eliminate $c$ to find\n",
    "$$\\begin{aligned}\n",
    "    c &= \\frac{a^2}{b} \\Rightarrow  \\\\\n",
    "    \\frac{a}{b} &= \\frac{a^2}{b^2-a^2} \\\\\n",
    "    ab^2 - a^3 &= a^2 b \\\\\n",
    "    \\frac{b^2}{a^2} - \\frac{b}{a} - 1 &= 0\n",
    "\\end{aligned}$$\n",
    "whose solution is\n",
    "$$\n",
    "    \\frac{b}{a} = \\frac{1 \\pm \\sqrt{5}}{2} = \\varphi\n",
    "$$\n",
    "where $\\varphi$ is the well known \"golden ratio\" (note that there are two values here, the most common definition of $\\varphi$ uses the $+$ branch but in fact you can use either depending on the application)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Back to the problem at hand, we now need to pick our new set of points.  Note that we only need one new point as the the other three are left-overs from the previous iteration.  Let us concentrate on the case where the extrema is between $[x_0, x_2]$.  Denote the new bracket values with $\\hat{\\quad}$ and identify\n",
    "$$\n",
    "    \\hat{x_0} = x_0, \\quad \\hat{x_2} = x_1, \\quad \\text{and} \\quad \\hat{x_3} = x_2.\n",
    "$$\n",
    "In this case we need to find $\\hat{x_1}$, in that case use the subsequent intervals $a$ and $\\hat{a_~}$ and equate\n",
    "$$\n",
    "    \\varphi \\hat{a~} = a \\Rightarrow \\varphi (\\hat{x_1} - \\hat{x_0}) = x_1 - x_0\n",
    "$$\n",
    "which in terms of the previous values can be solved for $\\hat{x_1}$ to lead to\n",
    "$$\n",
    "    \\hat{x_1} = \\frac{x_1 - x_0}{\\varphi} + x_0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In the alternative case we have the bracket $[x_1, x_3]$ and\n",
    "$$\n",
    "    \\hat{x_0} = x_1, \\quad \\hat{x_1} = x_2, \\quad \\text{and} \\quad \\hat{x_3} = x_3\n",
    "$$\n",
    "where we now need to find $\\hat{x_2}$.  Instead of using $\\hat{a~}$ we can use $\\hat{b~}$ and the relationship\n",
    "$$\n",
    "    \\varphi \\hat{c~} = c \\Rightarrow \\varphi (\\hat{x_2} - \\hat{x_1}) = x_2 - x_1\n",
    "$$\n",
    "which again can be manipulated to lead to the value of $\\hat{x_2}$ as\n",
    "$$\n",
    "    \\hat{x_2} = \\frac{x_2 - x_1}{\\varphi} + x_0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Algorithm\n",
    "1. Initialize bracket $[x_0,x_3]$\n",
    "1. Initialize points $x_1 = x_3 - \\frac{1}{\\varphi} \\cdot (x_3 - x_0)$ and $x_2 = x_0 + \\frac{1}{\\varphi} \\cdot (x_3 - x_0)$\n",
    "1. Loop\n",
    "   1. Evaluate $f_1$ and $f_2$\n",
    "   1. If $f_1 < f_2$ then we pick the left interval for the next iteration\n",
    "   1. and otherwise pick the right interval\n",
    "   1. Check size of bracket for convergence $x_3 - x_0 <$ `TOLERANCE`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# New Test Function!\n",
    "def f(t):\n",
    "    \"\"\"Simple function for minimization demos\"\"\"\n",
    "    return -3.0 * numpy.exp(-(t - 0.3)**2 / (0.1)**2) \\\n",
    "           +      numpy.exp(-(t - 0.6)**2 / (0.2)**2) \\\n",
    "           +      numpy.exp(-(t - 1.0)**2 / (0.2)**2) \\\n",
    "           +      numpy.sin(t)                        \\\n",
    "           -      2.0\n",
    "\n",
    "t = numpy.linspace(0, 2, 200)\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, f(t))\n",
    "axes.set_xlabel(\"t (days)\")\n",
    "axes.set_ylabel(\"People (N)\")\n",
    "axes.set_title(\"Decrease in Population due to SPAM Poisoning\")\n",
    "axes.set_xlim((0.0, 2.0))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def f(t):\n",
    "    \"\"\"Simple function for minimization demos\"\"\"\n",
    "    return -3.0 * numpy.exp(-(t - 0.3)**2 / (0.1)**2) \\\n",
    "           +      numpy.exp(-(t - 0.6)**2 / (0.2)**2) \\\n",
    "           +      numpy.exp(-(t - 1.0)**2 / (0.2)**2) \\\n",
    "           +      numpy.sin(t)                        \\\n",
    "           -      2.0\n",
    "\n",
    "phi = (1.0 + numpy.sqrt(5.0)) / 2.0\n",
    "\n",
    "TOLERANCE = 1e-4\n",
    "MAX_STEPS = 100\n",
    "\n",
    "x = [0.2, None, None, 0.5]\n",
    "x[1] = x[3] - 1.0 / phi * (x[3] - x[0])\n",
    "x[2] = x[0] + 1.0 / phi * (x[3] - x[0])\n",
    "\n",
    "t = numpy.linspace(0, 2, 200)\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, f(t))\n",
    "axes.set_xlabel(\"t (days)\")\n",
    "axes.set_ylabel(\"People (N)\")\n",
    "axes.set_title(\"Decrease in Population due to SPAM Poisoning\")\n",
    "\n",
    "success = False\n",
    "for n in range(1, MAX_STEPS + 1):\n",
    "    axes.plot(x[0], f(x[0]),'ko')\n",
    "    axes.plot(x[3], f(x[3]),'ko')\n",
    "    f_1 = f(x[1])\n",
    "    f_2 = f(x[2])\n",
    "    \n",
    "    if f_1 < f_2:\n",
    "        # Pick the left bracket\n",
    "        x_new = [x[0], None, x[1], x[2]]\n",
    "        x_new[1] = 1.0 / phi * (x[1] - x[0]) + x[0]\n",
    "    else:\n",
    "        # Pick the right bracket\n",
    "        x_new = [x[1], x[2], None, x[3]]\n",
    "        x_new[2] = (x[2] - x[1]) / phi + x[2]\n",
    "    x = x_new\n",
    "    \n",
    "    if numpy.abs(x[3] - x[0]) < TOLERANCE:\n",
    "        success = True\n",
    "        break\n",
    "        \n",
    "if success:\n",
    "    print(\"Success!\")\n",
    "    print(\"  t* = %s\" % str((x[3] + x[0]) / 2.0))\n",
    "    print(\"  f(t*) = %s\" % f((x[3] + x[0]) / 2.0))\n",
    "    print(\"  number of steps = %s\" % n)\n",
    "else:\n",
    "    print(\"Reached maximum number of steps!\")\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Interpolation Approach\n",
    "\n",
    "Successive parabolic interpolation - similar to secant method\n",
    "\n",
    "Basic idea:  Fit polynomial to function using three points, find it's minima, and guess new points based on that minima"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "1. What do we need to fit a polynomial $p_n(x)$ of degree $n \\geq 2$?\n",
    "\n",
    "2. How do we construct the polynomial $p_2(x)$?\n",
    "\n",
    "3. Once we have constructed $p_2(x)$ how would we find the minimum?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Algorithm\n",
    "\n",
    "Given $f(x)$ and $[x_0,x_1]$ - Note that unlike a bracket these will be a sequence of better approximations to the minimum.\n",
    "1. Initialize $x = [x_0, x_1, (x_0+x_1)/2]$\n",
    "1. Loop\n",
    "   1. Evaluate function $f(x)$\n",
    "   1. Use a polynomial fit to the function: \n",
    "   \n",
    "      $$p(x) = p_0 x^2 + p_1 x + p_2$$\n",
    "\n",
    "   1. Calculate the minimum: \n",
    "      \n",
    "      $$p'(x) = 2 p_0 x + p_1 = 0 \\quad \\Rightarrow \\quad x^\\ast = -p_1 / (2 p_0)$$\n",
    "\n",
    "   1. New set of points $x = [x_1, (x_0+x_1)/2, x^\\ast]$\n",
    "   1. Check tolerance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def f(t):\n",
    "    \"\"\"Simple function for minimization demos\"\"\"\n",
    "    return -3.0 * numpy.exp(-(t - 0.3)**2 / (0.1)**2) \\\n",
    "           +      numpy.exp(-(t - 0.6)**2 / (0.2)**2) \\\n",
    "           +      numpy.exp(-(t - 1.0)**2 / (0.2)**2) \\\n",
    "           +      numpy.sin(t)                        \\\n",
    "           -      2.0\n",
    "\n",
    "MAX_STEPS = 100\n",
    "TOLERANCE = 1e-4\n",
    "\n",
    "x = numpy.array([0.5, 0.2, (0.7) / 2.0])\n",
    "\n",
    "t = numpy.linspace(0, 2, 200)\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, f(t))\n",
    "axes.set_xlabel(\"t (days)\")\n",
    "axes.set_ylabel(\"People (N)\")\n",
    "axes.set_title(\"Decrease in Population due to SPAM Poisoning\")\n",
    "axes.plot(x[0], f(x[0]), 'ko')\n",
    "axes.plot(x[1], f(x[1]), 'ko')\n",
    "\n",
    "success = False\n",
    "for n in range(1, MAX_STEPS + 1):\n",
    "    axes.plot(x[2], f(x[2]), 'ko')\n",
    "    poly = numpy.polyfit(x, f(x), 2)\n",
    "    axes.plot(t, poly[0] * t**2 + poly[1] * t + poly[2], 'r--')\n",
    "    x[0] = x[1]\n",
    "    x[1] = x[2]\n",
    "    x[2] = -poly[1] / (2.0 * poly[0])\n",
    "    if numpy.abs(x[2] - x[1]) / numpy.abs(x[2]) < TOLERANCE:\n",
    "        success = True\n",
    "        break\n",
    "    \n",
    "if success:\n",
    "    print(\"Success!\")\n",
    "    print(\"  t* = %s\" % x[2])\n",
    "    print(\"  f(t*) = %s\" % f(x[2]))\n",
    "    print(\"  number of steps = %s\" % n)\n",
    "else:\n",
    "    print(\"Reached maximum number of steps!\")\n",
    "    \n",
    "axes.set_ylim((-5, 0.0))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Scipy Optimization\n",
    "\n",
    "Scipy contains a lot of ways for optimization!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.optimize as optimize\n",
    "print(optimize.golden(f, brack=(0.2, 0.25, 0.5)))"
   ]
  }
 ],
 "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
}