Module2/AD_with_dual_numbers_Julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An intro to Julia with Autodiff in 15 minutes\n",
    "\n",
    "Largely inspired from [Automatic Differentiation in 10 minutes with Julia](https://youtu.be/vAp6nUMrKYg)\n",
    "\n",
    "## 1- Dual numbers\n",
    "\n",
    "[source](https://en.wikipedia.org/wiki/Automatic_differentiation#Automatic_differentiation_using_dual_numbers)\n",
    "\n",
    "In linear algebra, the dual numbers extend the real numbers by adjoining one new element $\\varepsilon$ (epsilon) with the property $\\varepsilon^2 = 0$. Thus the multiplication of dual numbers is given by:\n",
    "$$(a+b\\varepsilon )(c+d\\varepsilon )=ac+(ad+bc)\\varepsilon,$$\n",
    "and addition is done componentwise:\n",
    "$$(a+b\\varepsilon )+(c+d\\varepsilon )=(a+c)+(b+d)\\varepsilon.$$\n",
    "\n",
    "Using [Julia](https://docs.julialang.org/en/v1/), we can create a data type for dual numbers as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct Dual <: Number\n",
    "    value::Float64\n",
    "    epsilon::Float64\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define addition, substraction and multiplication for dual numbers. To do this, we define a new [method](https://docs.julialang.org/en/v1/manual/methods/) to the [Julia Base](https://docs.julialang.org/en/v1/base/base/) operators `+,-,*`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Base: +, -, *\n",
    "+(z::Dual, w::Dual) = Dual(z.value+w.value, z.epsilon+w.epsilon)\n",
    "-(z::Dual, w::Dual) = Dual(z.value-w.value, z.epsilon-w.epsilon)\n",
    "*(z::Dual, w::Dual) = Dual(z.value*w.value, z.value*w.epsilon+z.epsilon*w.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience, we also define a new display function for dual numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Base.show(io::IO,x::Dual) = print(io,x.value,\" + \",x.epsilon,\" ε\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=Dual(2,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to start playing with dual numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a*a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a^7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This last two commands are rather surprising as we never define the power of a dual numbe! But it turns out that in Julia power is defined for any `x` supporting `*` as can be seen [here](https://github.com/JuliaLang/julia/blob/44fa15b1502a45eac76c9017af94332d4557b251/base/intfuncs.jl#L188) or by following the link given by the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@which a^7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2*a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To correct this problem, we need to:\n",
    "- convert the real `2` into the `Dual(2,0)`\n",
    "- tell Julia to make this conversion each time there is an expression implying dual numbers and reals\n",
    "\n",
    "This is what is called [conversion and promotion](https://docs.julialang.org/en/v1/manual/conversion-and-promotion/#conversion-and-promotion) and here it can be done as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Base: convert, promote_rule\n",
    "convert(::Type{Dual}, x::Real) = Dual(x,zero(x))\n",
    "promote_rule(::Type{Dual}, ::Type{<:Number}) = Dual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2*a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b=Dual(1,1)\n",
    "3*(a+b)^2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2- Automatic differentiation for polynomials\n",
    "\n",
    "We can now get derivatives for polynomials. Consider $P(x) = p_0+p_1x+p_2x^2+\\dots +p_n x^n$, note that since $\\varepsilon^2 =0$ ($\\varepsilon$ is nilpotent), we have $(x+\\varepsilon y)^k = x^k + kx^{k-1}\\varepsilon y$. Hence, we have\n",
    "$$\n",
    "P(x+\\varepsilon y) = P(x) +  \\left(p_1 + 2p_2 x+ \\dots np_n x^{n-1}\\right)y\\varepsilon = P(x) + P'(x)y\\varepsilon \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "3*(1+a)^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "value(z::Dual) = z.value\n",
    "epsilon(z::Dual) = z.epsilon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epsilon(3*(1+a)^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define a simple [function](https://docs.julialang.org/en/v1/manual/functions/) to compute the derivative of a polynomials at a given point as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function derivative(f,x::Real)\n",
    "    epsilon(f(Dual(x,1)))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "derivative(x->1+x+3x^2,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3- Going further and dealing with $\\sqrt{}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a^(1/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To deal with this problem, we follow the Babylonian as described in this nice [tutorial](https://github.com/JuliaAcademy/JuliaTutorials/blob/master/introductory-tutorials/intro-to-julia/AutoDiff.ipynb)\n",
    "> Repeat $ t \\leftarrow  \\frac{1}{2}\\left(t+\\frac{x}{t}\\right)$ until $t$ converges to $\\sqrt{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function Babylonian(x; N = 10) \n",
    "    t = (1+x)/2\n",
    "    for i = 2:N; t=(t+x/t)/2  end    \n",
    "    t\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Babylonian(2), √2 # Type \\sqrt+<tab> to get the symbol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Babylonian algorithm uses only addition and division, hence we need to define a division for dual numbers:\n",
    "$$\\frac{a+b\\varepsilon}{c+d\\varepsilon}= \\frac{a}{c}\\left(1+\\frac{b}{a}\\varepsilon\\right)\\left(1-\\frac{d}{c}\\varepsilon\\right)= \\frac{a}{c} +\\frac{bc-ad}{c^2}\\varepsilon.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Base:/\n",
    "/(x::Dual, y::Dual) = Dual(x.value/y.value, (y.value*x.epsilon - x.value*y.epsilon)/y.value^2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1/(1+a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Babylonian(a), √2, 0.5/√2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "derivative(x->Babylonian(x),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function dBabylonian(x; N = 10) \n",
    "    t = (1+x)/2\n",
    "    dt = 1/2\n",
    "    for i = 1:N;  \n",
    "        t = (t+x/t)/2; \n",
    "        dt = (dt+(t-x*dt)/t^2)/2; \n",
    "    end    \n",
    "    dt\n",
    "end  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 2; dBabylonian(x), .5/√x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "derivative(x->Babylonian(x)^4,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "derivative(x->(1+3*Babylonian(x))^3/Babylonian(x),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using Pkg; Pkg.add(\"ForwardDiff\")\n",
    "using ForwardDiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ForwardDiff.derivative(sqrt, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ForwardDiff.derivative(Babylonian, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ForwardDiff.derivative(x->(1+3*sqrt(x))^3/sqrt(x),2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more on dual numbers, have a look at the (not active anymore) Julia package [DualNumbers.jl](https://github.com/JuliaDiff/DualNumbers.jl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.2",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}