Cubic-Newton.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Super \"Newton\" method\n",
    "\n",
    "Newton's method corresponds to approximately solving $f(x)=0$ by doing a first-order Taylor approximation $f(x+\\delta)\\approx f(x)+\\delta f'(x)$ and solving for the root of this *line* to get our step $\\delta$.\n",
    "\n",
    "Suppose that instead we do a Taylor expansion to *second* order, finding the root of a *quadratic*:\n",
    "\n",
    "$$\n",
    "f(x+\\delta) \\approx f(x) + f'(x) \\delta + \\frac{f''(x)}{2} \\delta^2 = 0\n",
    "$$\n",
    "\n",
    "This is a quadratic equation, so we can use the quadratic formula to obtain:\n",
    "$$\n",
    "\\delta =  - \\frac{2f}{f' \\pm \\sqrt{(f')^2 - 2ff''}}\n",
    "$$\n",
    "where we've applied a couple of tricks to get the \"right\" solution:\n",
    "\n",
    "1. We want to make sure our equation is reasonable as $f'' \\to 0$, in which case we should get an ordinary Newton step.   This is a problem for the \"ordinary\" quadratic formula because the coefficient of $\\delta^2$ is $f''$, and this goes in the denominator.  We can get a \"better\" quadratic formula for this case by dividing by $1/\\delta^2$ and writing the quadratic equation for $1/\\delta$.\n",
    "\n",
    "2. The quadratic formula gives two roots.  We only want the *smaller* $\\delta$ root, because we want to take small steps, and that determines the sign.  This means that we should choose the $\\pm$ with the **same sign** as $f'$.\n",
    "\n",
    "Note also that the quadratic formula can go \"crazy\" (give complex steps for real functions) if $(f')^2 < 2ff''$.   However, this will never happen if we have a good enough initial guess, so that $f(x)$ is sufficiently small.   In general, higher-order Newton methods tend to be more sensitive to the initial guess!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "using ForwardDiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "newton (generic function with 1 method)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function newton(f, x, n)\n",
    "    for i = 1:n\n",
    "        fx = f(x)\n",
    "        fx′ = ForwardDiff.derivative(f, x)\n",
    "        x = x - fx / fx′\n",
    "    end\n",
    "    return x\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "supernewton (generic function with 1 method)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function supernewton(f, x, n)\n",
    "    for i = 1:n\n",
    "        fx = f(x)\n",
    "        fx′ = ForwardDiff.derivative(f, x)\n",
    "        fx″ = ForwardDiff.derivative(x -> ForwardDiff.derivative(f, x), x)\n",
    "        x = x - 2fx / (fx′ + copysign(sqrt(fx′^2 - 2fx*fx″), fx′))\n",
    "    end\n",
    "    return x\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2.6305807789924294\n",
      "-9.512351136373354\n",
      "-30.157759191560796\n",
      "-92.09398335713593\n",
      "-277.90265585386135\n",
      "-835.3286733440375\n",
      "-2507.6067258145663\n",
      "-7524.4408832261515\n",
      "-12041.199826559248\n",
      "-12041.199826559248\n"
     ]
    }
   ],
   "source": [
    "f(x) = cos(x) - x\n",
    "setprecision(BigFloat, 40000)\n",
    "for i = 1:10\n",
    "    println(Float64(log10(abs(f(supernewton(f, big(1), i))))))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We roughly **triple** the number of digits (**cube** the error) on each step.\n",
    "\n",
    "Basically, this is because expanding to second order neglects the $O(\\delta^3)$ term in the Taylor expansion, so it solves $f(x+\\delta)=0$ up to an $O(\\delta^3)$ error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.0",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}