13_bisections_newton.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Foundations of Computational Economics #13\n",
    "\n",
    "by Fedor Iskhakov, ANU\n",
    "\n",
    "<img src=\"_static/img/dag3logo.png\" style=\"width:256px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Two very important algorithms for solving equations\n",
    "\n",
    "<img src=\"_static/img/lab.png\" style=\"width:64px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Classic algorithm for equation solving\n",
    "\n",
    "1. Bisection method  \n",
    "1. Newton-Raphson method  \n",
    "\n",
    "\n",
    "Solve equations of the form $ f(x) = 0 $\n",
    "\n",
    "Focus on the scalar case today."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Bisection method for solving equations\n",
    "\n",
    "Solve equation $ f(x)=0 $, conditional on $ x \\in [a,b] \\subset \\mathbb{R} $ such that $ f(a)f(b)<0 $\n",
    "\n",
    "Algorithm: similar to binary search, but in **continuous space**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def bisection(f,a=0,b=1,tol=1e-6,maxiter=100,callback=None):\n",
    "    '''Bisection method for solving equation f(x)=0\n",
    "    on the interval [a,b], with given tolerance and number of iterations.\n",
    "    Callback function is invoked at each iteration if given.\n",
    "    '''\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: -4*x**3+5*x+1\n",
    "a,b = -3,-.5  # upper and lower limits\n",
    "x = bisection(f,a,b)\n",
    "print('Solution is x=%1.3f, f(x)=%1.12f' % (x,f(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# make nice plot\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = [12, 8]\n",
    "xd = np.linspace(a,b,1000)  # x grid\n",
    "plt.plot(xd,f(xd),c='red')  # plot the function\n",
    "plt.plot([a,b],[0,0],c='black')  # plot zero line\n",
    "ylim=[f(a),min(f(b),0)]\n",
    "plt.plot([a,a],ylim,c='grey')  # plot lower bound\n",
    "plt.plot([b,b],ylim,c='grey')  # plot upper bound\n",
    "def plot_step(x,**kwargs):\n",
    "    plot_step.counter += 1\n",
    "    plt.plot([x,x],ylim,c='grey')\n",
    "plot_step.counter = 0  # new public attribute\n",
    "bisection(f,a,b,callback=plot_step)\n",
    "print('Converged in %d steps'%plot_step.counter)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def bisection(f,a=0,b=1,tol=1e-6,maxiter=100,callback=None):\n",
    "    '''Bisection method for solving equation f(x)=0\n",
    "    on the interval [a,b], with given tolerance and number of iterations.\n",
    "    Callback function is invoked at each iteration if given.\n",
    "    '''\n",
    "    if f(a)*f(b)>0:\n",
    "        raise ValueError('Function has the same sign at the bounds')\n",
    "    for i in range(maxiter):\n",
    "        err = abs(b-a)\n",
    "        if err<tol: break\n",
    "        x = (a+b)/2\n",
    "        a,b = (x,b) if f(a)*f(x)>0 else (a,x)\n",
    "        if callback != None: callback(err=err,x=x,iter=i)\n",
    "    else:\n",
    "        raise RuntimeError('Failed to converge in %d iterations'%maxiter)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Newton-Raphson (Newton) method\n",
    "\n",
    "General form $ f(x)=0 $\n",
    "\n",
    "- Equation solving  \n",
    "- Finding maximum/minimum based on FOC, then $ f(x)=Q'(x) $  \n",
    "\n",
    "\n",
    "Algorithm:\n",
    "1. Start with some good guess $ x_0 $ not too far from the solution\n",
    "2. Newton step: $ x_{i+1} = x_i - \\frac{f(x_i)}{f'(x_i)} $\n",
    "3. Iterate until convergence in some metric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Derivation for Newton method using Taylor series expansion\n",
    "\n",
    "$$\n",
    "f(x) = \\sum_{k=0}^{\\infty} \\frac{f^{(k)}(x_0)}{k!} (x-x_0)^k\n",
    "$$\n",
    "\n",
    "Take first two terms, assume $ f(x) $ is solution, and let\n",
    "$ x_0=x_i $ and $ x=x_{i+1} $\n",
    "\n",
    "$$\n",
    "0 = f(x) = f(x_i) + f'(x_i) (x_{i+1}-x_i) \\quad \\Rightarrow \\quad x_{i+1} = x_i - \\frac{f(x_i)}{f'(x_i)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def newton(fun,grad,x0,tol=1e-6,maxiter=100,callback=None):\n",
    "    '''Newton method for solving equation f(x)=0\n",
    "    with given tolerance and number of iterations.\n",
    "    Callback function is invoked at each iteration if given.\n",
    "    '''\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: -4*x**3+5*x+1\n",
    "g = lambda x: -12*x**2+5\n",
    "x = newton(f,g,x0=-2.5,maxiter=7)\n",
    "print('Solution is x=%1.3f, f(x)=%1.12f' % (x,f(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# make nice seriest of plots\n",
    "a,b = -3,-.5  # upper and lower limits\n",
    "xd = np.linspace(a,b,1000)  # x grid\n",
    "def plot_step(x0,x1,iter,**kwargs):\n",
    "    plot_step.counter += 1\n",
    "    if iter<10:\n",
    "        plt.plot(xd,f(xd),c='red')  # plot the function\n",
    "        plt.plot([a,b],[0,0],c='black')  # plot zero line\n",
    "        ylim = [min(f(b),0),f(a)]\n",
    "        plt.plot([x0,x0],ylim,c='grey') # plot x0\n",
    "        l = lambda z: g(x0)*(z - x1)\n",
    "        plt.plot(xd,l(xd),c='green')  # plot the function\n",
    "        plt.ylim(bottom=10*f(b))\n",
    "        plt.title('Iteration %d'%(iter+1))\n",
    "        plt.show()\n",
    "plot_step.counter = 0  # new public attribute\n",
    "newton(f,g,x0=-2.5,callback=plot_step)\n",
    "print('Converged in %d steps'%plot_step.counter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def newton(fun,grad,x0,tol=1e-6,maxiter=100,callback=None):\n",
    "    '''Newton method for solving equation f(x)=0\n",
    "    with given tolerance and number of iterations.\n",
    "    Callback function is invoked at each iteration if given.\n",
    "    '''\n",
    "    for i in range(maxiter):\n",
    "        x1 = x0 - fun(x0)/grad(x0)\n",
    "        err = abs(x1-x0)\n",
    "        if callback != None: callback(err=err,x0=x0,x1=x1,iter=i)\n",
    "        if err<tol: break\n",
    "        x0 = x1\n",
    "    else:\n",
    "        raise RuntimeError('Failed to converge in %d iterations'%maxiter)\n",
    "    return (x0+x1)/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Rate of convergence of the two methods\n",
    "\n",
    "- How fast does a solution method converge on the root of the equation?  \n",
    "- Rate of convergence = the rate of decrease of the bias (difference between current guess and the solution)  \n",
    "- Can be approximated by the rate of decrease of the error in the stopping criterion  \n",
    "\n",
    "\n",
    "Bisections: **linear convergence**\n",
    "\n",
    "Newton: **quadratic convergence**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def print_err(iter,err,**kwargs):\n",
    "    x = kwargs['x'] if 'x' in kwargs.keys() else kwargs['x0']\n",
    "    print('{:4d}:  x = {:17.14f}  err = {:8.6e}'.format(iter,x,err))\n",
    "\n",
    "print('Newton method')\n",
    "newton(f,g,x0=-2.5,callback=print_err,tol=1e-10)\n",
    "\n",
    "print('Bisection method')\n",
    "bisection(f,a=-3,b=-0.5,callback=print_err,tol=1e-10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Measuring complexity of Newton and bisection methods\n",
    "\n",
    "- What is the size of input $ n $?  \n",
    "- Desired precision of the solution!  \n",
    "- Thus, attention to the errors in the solution as algorithm proceeds  \n",
    "- Rate of convergence is part of the computational complexity of the algorithms  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Computational complexity of Newton method\n",
    "\n",
    "- Calculating a root of a function f(x) **with n-digit precision**  \n",
    "- Provided that a good initial approximation is known  \n",
    "- Is $ O((logn)F(n)) $, where $ F(n) $ is the cost of  \n",
    "- calculating $ f(x)/f'(x) $ with $ n $-digit precision  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Further learning resources\n",
    "\n",
    "- On computational complexity of Newton method\n",
    "  [https://m.tau.ac.il/~tsirel/dump/Static/knowino.org/wiki/Newton’s_method.html#Computational_complexity](https://m.tau.ac.il/~tsirel/dump/Static/knowino.org/wiki/Newton's_method.html#Computational_complexity)  \n",
    "- “Improved convergence and complexity analysis of Newton’s method for solving equations”\n",
    "  [https://www.tandfonline.com/doi/abs/10.1080/00207160601173431?journalCode=gcom20](https://www.tandfonline.com/doi/abs/10.1080/00207160601173431?journalCode=gcom20)  "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "date": 1597824775.202882,
  "filename": "13_bisections_newton.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Foundations of Computational Economics #13"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}