Skip to content

Instantly share code, notes, and snippets.

@z-m-k

z-m-k/Julia's AutoDiff.ipynb

Created August 7, 2015 13:00
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save z-m-k/03e1f3f1452d3d9e9595 to your computer and use it in GitHub Desktop.
Save z-m-k/03e1f3f1452d3d9e9595 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
Julia's AutoDiff.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"using DataFrames"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"numerical_derivatives (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function numerical_derivatives(fun::Function, x::Array)\n",
" function derivative(f0::Array, f, x)\n",
" grad=fill(0.0, size(v)[1], length(x))\n",
" for i in 1:length(x)\n",
" h=sqrt(eps())*maximum((abs(x[i]),1))\n",
" xph = copy(x)\n",
" xph[i]+=h\n",
" dx = (xph[i] - x[i])#::Float64\n",
" f1 = fun(xph...)#::Float64\n",
" grad[:,i]=(f1 - f0) / dx\n",
" end\n",
" return [f0, grad...]\n",
" end\n",
" function derivative(f0::Number, f, x)\n",
" grad=fill(0.0, length(x))\n",
" for i in 1:length(x)\n",
" h=sqrt(eps())*maximum((abs(x[i]),1))\n",
" xph = copy(x)\n",
" xph[i]+=h\n",
" dx = (xph[i] - x[i])#::Float64\n",
" f1 = fun(xph...)#::Float64\n",
" grad[i]=(f1 - f0) / dx\n",
" end\n",
" return [f0, grad...]\n",
" end\n",
" f0=fun(x...)\n",
" return derivative(f0, fun, x)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"autodiff_derivatives (generic function with 1 method)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using DualNumbers\n",
"function autodiff_derivatives(f, x)\n",
" function derivative(v::Array, f, x)\n",
" grad=fill(0.0, size(v)[1], length(x))\n",
" for i=1:length(x)\n",
" xp=copy(xp0)\n",
" xp[i]=dual(x[i], 1)\n",
" grad[:,i]=epsilon(f(xp...))\n",
" end\n",
" return [v, grad...]\n",
" end\n",
" function derivative(v::Number, f, x)\n",
" grad=fill(0.0, length(x))\n",
" grad=fill(0.0, length(x))\n",
" for i=1:length(x)\n",
" xp=copy(xp0)\n",
" xp[i]=dual(x[i], 1)\n",
" grad[i]=epsilon(f(xp...))\n",
" end\n",
" return [v, grad...]\n",
" end\n",
" xp0=dual(x, 0)\n",
" v=f(x...)\n",
" return derivative(v, f, x)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"#Example\n",
"Consider $f(x,y)=4.25x^2+xy^x$\n",
"\n",
"Then:\n",
"\n",
"$\\frac{\\partial f}{\\partial x}=8.5x+y^x+xy^x\\log{y}$\n",
"\n",
"$\\frac{\\partial f}{\\partial y}=x^2y^{x-1}$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f(x,y) = 4.25*x^2 + x*y^x\n",
"∂f∂x(x,y) = 8.5x + y^x + x*y^x*log(y)\n",
"∂f∂y(x,y) = x^2*y^(x-1)\n",
"\n",
"x=2.5\n",
"y=3\n",
"hand_written=[f(x,y), ∂f∂x(x,y), ∂f∂y(x,y)]\n",
"auto_diff =autodiff_derivatives(f, [x, y])\n",
"nume_diff =numerical_derivatives(f, [x,y]);"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"data-frame\"><tr><th></th><th>Index</th><th>Hand</th><th>Auto</th><th>Numerical</th></tr><tr><th>1</th><td>$f()$</td><td>65.53364317029974</td><td>65.53364317029974</td><td>65.53364317029974</td></tr><tr><th>2</th><td>$\\frac{∂f}{∂x}$</td><td>79.65263405845546</td><td>79.65263405845546</td><td>79.65263595581055</td></tr><tr><th>3</th><td>$\\frac{∂f}{∂y}$</td><td>32.47595264191645</td><td>32.47595264191645</td><td>32.475952784220375</td></tr></table>"
],
"text/plain": [
"3x4 DataFrame\n",
"| Row | Index | Hand | Auto | Numerical |\n",
"|-----|-----------------|---------|---------|-----------|\n",
"| 1 | \"\\$f()\\$\" | 65.5336 | 65.5336 | 65.5336 |\n",
"| 2 | \"\\$\\\\frac{\\u2202f}{\\u2202x}\\$\" | 79.6526 | 79.6526 | 79.6526 |\n",
"| 3 | \"\\$\\\\frac{\\u2202f}{\\u2202y}\\$\" | 32.476 | 32.476 | 32.476 |"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data=DataFrame(Index=[\"\\$f()\\$\", \"\\$\\\\frac{∂f}{∂x}\\$\", \"\\$\\\\frac{∂f}{∂y}\\$\"], \n",
" Hand=hand_written, \n",
" Auto=auto_diff, \n",
" Numerical=nume_diff)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Differences between them:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"data-frame\"><tr><th></th><th>Index</th><th>Auto</th><th>Numerical</th></tr><tr><th>1</th><td>$f()$</td><td>0.0</td><td>0.0</td></tr><tr><th>2</th><td>$\\frac{∂f}{∂x}$</td><td>0.0</td><td>-1.8973550908185643e-6</td></tr><tr><th>3</th><td>$\\frac{∂f}{∂y}$</td><td>0.0</td><td>-1.42303925088072e-7</td></tr></table>"
],
"text/plain": [
"3x3 DataFrame\n",
"| Row | Index | Auto | Numerical |\n",
"|-----|-----------------|------|-------------|\n",
"| 1 | \"\\$f()\\$\" | 0.0 | 0.0 |\n",
"| 2 | \"\\$\\\\frac{\\u2202f}{\\u2202x}\\$\" | 0.0 | -1.89736e-6 |\n",
"| 3 | \"\\$\\\\frac{\\u2202f}{\\u2202y}\\$\" | 0.0 | -1.42304e-7 |"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"DataFrame(Index=[\"\\$f()\\$\", \"\\$\\\\frac{∂f}{∂x}\\$\", \"\\$\\\\frac{∂f}{∂y}\\$\"], \n",
" Auto=data[:Hand]-data[:Auto], \n",
" Numerical=data[:Hand]-data[:Numerical])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So AutoDiff is __exact__."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#Works with more complex functions \n",
"\n",
"Let $f(x_{1..N},y_{1..N})=4.5a\\sum_{i=1}^N{x^{ia}_i}+ab\\sum_{i=1}^N{i x_i y_i^{\\frac{a}{b}x_i}}$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(\n",
"1x3 Array{Float64,2}:\n",
" 0.355383 1.25013 0.376598,\n",
"\n",
"1x3 Array{Float64,2}:\n",
" 0.512889 0.467054 3.45291)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a=2.5; b=3;\n",
"X=exp(randn(3))\n",
"Y=exp(randn(3))\n",
"X', Y'"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"data-frame\"><tr><th></th><th>Index</th><th>Auto</th><th>Numerical</th></tr><tr><th>1</th><td>$f$</td><td>58.37780401230243</td><td>58.37780401230243</td></tr><tr><th>2</th><td>$\\frac{∂f}{∂a}$</td><td>36.87103006799046</td><td>36.87103042602539</td></tr><tr><th>3</th><td>$\\frac{∂f}{∂b}$</td><td>8.490680931708148</td><td>8.490681012471518</td></tr></table>"
],
"text/plain": [
"3x3 DataFrame\n",
"| Row | Index | Auto | Numerical |\n",
"|-----|-----------------|---------|-----------|\n",
"| 1 | \"\\$f\\$\" | 58.3778 | 58.3778 |\n",
"| 2 | \"\\$\\\\frac{\\u2202f}{\\u2202a}\\$\" | 36.871 | 36.871 |\n",
"| 3 | \"\\$\\\\frac{\\u2202f}{\\u2202b}\\$\" | 8.49068 | 8.49068 |"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#One line...\n",
"f(a,b)=4.5*a*sum([x^(a*i) for (i,x) in enumerate(X)])+sum([a*b*i*xy[1]*xy[2]^(a/b*xy[1]) for (i,xy) in enumerate(zip(X,Y))])\n",
"one=DataFrame(Index=[\"\\$f\\$\", \"\\$\\\\frac{∂f}{∂a}\\$\", \"\\$\\\\frac{∂f}{∂b}\\$\"], \n",
" Auto=autodiff_derivatives(f, [a, b]), \n",
" Numerical=numerical_derivatives(f, [a, b]))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"data-frame\"><tr><th></th><th>Index</th><th>Auto</th><th>Numerical</th></tr><tr><th>1</th><td>$f$</td><td>58.37780401230242</td><td>58.37780401230242</td></tr><tr><th>2</th><td>$\\frac{∂f}{∂a}$</td><td>36.87103006799046</td><td>36.87103061676025</td></tr><tr><th>3</th><td>$\\frac{∂f}{∂b}$</td><td>8.490680931708148</td><td>8.490681171417236</td></tr></table>"
],
"text/plain": [
"3x3 DataFrame\n",
"| Row | Index | Auto | Numerical |\n",
"|-----|-----------------|---------|-----------|\n",
"| 1 | \"\\$f\\$\" | 58.3778 | 58.3778 |\n",
"| 2 | \"\\$\\\\frac{\\u2202f}{\\u2202a}\\$\" | 36.871 | 36.871 |\n",
"| 3 | \"\\$\\\\frac{\\u2202f}{\\u2202b}\\$\" | 8.49068 | 8.49068 |"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#Of course this could be also written...\n",
"function f(a,b)\n",
" out=0.0\n",
" for i=1:endof(X)\n",
" @inbounds out+=4.5*a*X[i]^(a*i)+a*b*i*X[i]*Y[i]^(a/b*X[i])\n",
" end\n",
" out\n",
"end\n",
"long=DataFrame(Index=[\"\\$f\\$\", \"\\$\\\\frac{∂f}{∂a}\\$\", \"\\$\\\\frac{∂f}{∂b}\\$\"], \n",
" Auto=autodiff_derivatives(f, [a, b]), \n",
" Numerical=numerical_derivatives(f, [a, b]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Differences (due to order of operations on floats)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"data-frame\"><tr><th></th><th>Index</th><th>Auto</th><th>Numerical</th></tr><tr><th>1</th><td>$f$</td><td>7.105427357601002e-15</td><td>7.105427357601002e-15</td></tr><tr><th>2</th><td>$\\frac{∂f}{∂a}$</td><td>0.0</td><td>-1.9073485901799359e-7</td></tr><tr><th>3</th><td>$\\frac{∂f}{∂b}$</td><td>0.0</td><td>-1.5894571880892272e-7</td></tr></table>"
],
"text/plain": [
"3x3 DataFrame\n",
"| Row | Index | Auto | Numerical |\n",
"|-----|-----------------|-------------|-------------|\n",
"| 1 | \"\\$f\\$\" | 7.10543e-15 | 7.10543e-15 |\n",
"| 2 | \"\\$\\\\frac{\\u2202f}{\\u2202a}\\$\" | 0.0 | -1.90735e-7 |\n",
"| 3 | \"\\$\\\\frac{\\u2202f}{\\u2202b}\\$\" | 0.0 | -1.58946e-7 |"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"DataFrame(Index=[\"\\$f\\$\", \"\\$\\\\frac{∂f}{∂a}\\$\", \"\\$\\\\frac{∂f}{∂b}\\$\"], \n",
" Auto=one[:Auto]-long[:Auto], \n",
" Numerical=one[:Numerical]-long[:Numerical])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.3.9",
"language": "julia",
"name": "julia-0.3"
},
"language_info": {
"name": "julia",
"version": "0.3.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment