"Vectorized" Newton-Raphson vs SciPy.optimize Newton-Raphson

newton-raphson.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "import numpy as np\n",
    "from pandas import Series, DataFrame\n",
    "from scipy.optimize import newton\n",
    "from typing import Callable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "LRA = list()\n",
    "LRA.append(\"{levelname:s}\")\n",
    "LRA.append(\"{message:s}\")\n",
    "logging.basicConfig(level=logging.DEBUG, format=\" - \".join(LRA), style=\"{\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Function and Derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x: Series) -> Series:\n",
    "    return 6*x**5 - 5*x**4 - 4*x**3 + 3*x**2\n",
    "\n",
    "def df(x: Series) -> Series:\n",
    "    return 30*x**4 - 20*x**3 - 12*x**2 + 6*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vectorized Newton-Raphson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### while(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vectorized_newton(f: Callable, df: Callable, x0: Series, epsilon: float=1.48e-08, kappa: int=50):\n",
    "    logging.debug(\"expect convergence within ε={} after κ={} iterations\".format(epsilon, kappa))\n",
    "    i = 0\n",
    "    x = x0.copy()  # retain x0 values\n",
    "    delta = np.abs(0 - f(x)) # initial state\n",
    "    mask = delta > epsilon  # index mask (updated as converence reached)\n",
    "    try:\n",
    "        while any(mask):\n",
    "            state = DataFrame({\"Δ\": delta, \">ε\": mask})\n",
    "            logging.debug(\"iter {}:\\n{}\".format(i, state))\n",
    "            i += 1\n",
    "            assert i < kappa\n",
    "            x[mask] = x[mask] - f(x[mask])/df(x[mask])\n",
    "            delta[mask] = np.abs(0 - f(x))\n",
    "            mask = delta > epsilon\n",
    "    except AssertionError:\n",
    "        failed = DataFrame({\"x0\": x0[mask], \"x\": x[mask], \"Δ\": delta[mask]})\n",
    "        logging.error(\"convergence failed within {} iterations:\\n{}\".format(kappa, failed))\n",
    "        return x\n",
    "    else:\n",
    "        logging.info(\"all converged within {} iterations\".format(i+1))\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### for(range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vectorized_newton(f: Callable, df: Callable, x0: Series, epsilon: float=1.48e-08, kappa: int=50):\n",
    "    logging.debug(\"expect convergence within ε={} after κ={} iterations\".format(epsilon, kappa))\n",
    "    x = x0.copy()  # retain x0 values\n",
    "    delta = np.abs(0 - f(x))  # initial state\n",
    "    mask = delta > epsilon  # index mask (updated as converence reached)\n",
    "    for i in range(kappa):\n",
    "        state = DataFrame({\"Δ\": delta, \">ε\": mask})\n",
    "        logging.debug(\"iter {}:\\n{}\".format(i, state))\n",
    "        if not any(mask):\n",
    "            break\n",
    "        x[mask] = x[mask] - f(x[mask])/df(x[mask])\n",
    "        delta[mask] = np.abs(0 - f(x))\n",
    "        mask = delta > epsilon\n",
    "    if any(mask):\n",
    "        state = DataFrame({\"x0\": x0[mask], \"x\": x[mask], \"Δ\": delta[mask]})\n",
    "        logging.error(\"convergence failed within {} iterations:\\n{}\".format(kappa, state))\n",
    "    else:\n",
    "        logging.info(\"all converged within {} iterations\".format(i+1))\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SciPy Newton-Raphson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scipy_newton(x0: Series, **kwargs) -> Series:\n",
    "    x = x0.copy()\n",
    "    for i in x0.index:\n",
    "        x[i] = newton(f, x0[i], df ,**kwargs)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = Series([0.01, .5, .9], index=list(\"abc\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "DEBUG - expecting convergence within ε=1.48e-08 after κ=10 iterations\n",
      "DEBUG - iter 0:\n",
      "     >ε         Δ\n",
      "a  True  0.000296\n",
      "b  True  0.125000\n",
      "c  True  0.223560\n",
      "DEBUG - iter 1:\n",
      "     >ε         Δ\n",
      "a  True  0.000073\n",
      "b  True  0.094080\n",
      "c  True  1.725605\n",
      "DEBUG - iter 2:\n",
      "     >ε         Δ\n",
      "a  True  0.000018\n",
      "b  True  0.001213\n",
      "c  True  0.464877\n",
      "DEBUG - iter 3:\n",
      "     >ε         Δ\n",
      "a  True  0.000005\n",
      "b  True  0.000001\n",
      "c  True  0.093756\n",
      "DEBUG - iter 4:\n",
      "      >ε             Δ\n",
      "a   True  1.141123e-06\n",
      "b  False  1.804334e-12\n",
      "c   True  8.175036e-03\n",
      "DEBUG - iter 5:\n",
      "      >ε             Δ\n",
      "a   True  2.851629e-07\n",
      "b  False  1.804334e-12\n",
      "c   True  8.467543e-05\n",
      "DEBUG - iter 6:\n",
      "      >ε             Δ\n",
      "a   True  7.127603e-08\n",
      "b  False  1.804334e-12\n",
      "c  False  9.407033e-09\n",
      "DEBUG - iter 7:\n",
      "      >ε             Δ\n",
      "a   True  1.781718e-08\n",
      "b  False  1.804334e-12\n",
      "c  False  9.407033e-09\n",
      "DEBUG - iter 8:\n",
      "      >ε             Δ\n",
      "a  False  4.454065e-09\n",
      "b  False  1.804334e-12\n",
      "c  False  9.407033e-09\n",
      "INFO - all converged within 9 iterations\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "a    0.000039\n",
       "b    0.628667\n",
       "c    1.000000\n",
       "dtype: float64"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectorized_newton(f, df, x0, kappa=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "RuntimeError",
     "evalue": "Failed to converge after 10 iterations, value is 9.632789518030636e-06",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-ebe1fb7cc604>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mscipy_newton\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxiter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-10-f1072022db72>\u001b[0m in \u001b[0;36mscipy_newton\u001b[0;34m(x0, **kwargs)\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m         \u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnewton\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdf\u001b[0m \u001b[0;34m,\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Developer/smcveigh-phunware/utm/venv/lib/python3.6/site-packages/scipy/optimize/zeros.py\u001b[0m in \u001b[0;36mnewton\u001b[0;34m(func, x0, fprime, args, tol, maxiter, fprime2)\u001b[0m\n\u001b[1;32m    213\u001b[0m             \u001b[0mq1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    214\u001b[0m     \u001b[0mmsg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"Failed to converge after %d iterations, value is %s\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmaxiter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 215\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    216\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    217\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: Failed to converge after 10 iterations, value is 9.632789518030636e-06"
     ]
    }
   ],
   "source": [
    "scipy_newton(x0, maxiter=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}