derivatives.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Derivatives\n",
    "\n",
    "This notebook was automatically generated from the Algorithms for Optimization source code. Each cell generates a figure from the original text. While this code is not optimized for use in lectures, we provide it here to be adapted for such projects. We hope you find it useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "include(\"support_code.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\n",
    "f = x->sin(x)\n",
    "\n",
    "x₀ = 0.5\n",
    "ymin = -1.1\n",
    "lo = -2\n",
    "hi = 3\n",
    "\n",
    "function get_axis(h)\n",
    "a, b = x₀, x₀ + h\n",
    "fa, fb = f(a), f(b)\n",
    "m = (fb - fa) / (b-a)\n",
    "\n",
    "p = Plots.Plot[]\n",
    "push!(p, Plots.Linear(f, (lo,hi), style=\"thick, black\"))\n",
    "push!(p, Plots.Linear([a,a],[fa, ymin], style=\"solid, black, mark=none\"))\n",
    "push!(p, Plots.Linear([b,b],[fb, ymin], style=\"solid, black, mark=none\"))\n",
    "push!(p, Plots.Linear([lo,hi],[ymin, ymin], style=\"solid, black, mark=none\"))\n",
    "push!(p, Plots.Linear([lo, hi],[fa + (lo-a)*m, fa + (hi-a)*m], style=\"ultra thick, solid, pastelBlue, mark=none\"))\n",
    "push!(p, Plots.Scatter([a, b],[fa, fb], style=\"solid, mark=*, mark options={draw=black, fill=black}\"))\n",
    "\n",
    "Axis(p, style=\"axis lines=left, axis line style = {transparent}, xtick={$(x₀ + h/2)}, xticklabels={\\$h\\$}, ytick=\\\\empty, xtick style={draw=none}, ytick style={draw=none}\") #$\n",
    "end\n",
    "\n",
    "g = GroupPlot(3, 1, groupStyle=\"horizontal sep=1cm\", style=\"ymin=-1.1, ymax=1.1, width=6cm\")\n",
    "push!(g, get_axis(2.0))\n",
    "push!(g, get_axis(1.0))\n",
    "push!(g, get_axis(0.5))\n",
    "g\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "% \tplot(Axis(\n",
    "% \t\t[Plots.Contour(x->exp(0-x[1]^2-x[2]^2)*(x[1]-2*x[2]), (-2,2), (-2,2), levels=[-0.1,-0.3, -0.5]),\n",
    "% \t\t Plots.Scatter([-0.75], [1.0], style=\"black, mark=*, mark options={draw=black}\"),\n",
    "% \t\t Plots.Command(\"\\\\draw [->, ultra thick] (axis cs:-0.75,1) -- (axis cs:-0.25,0.5);\"),\n",
    "% \t\t Plots.Command(\"\\\\node at (axis cs:-0.75,1) [anchor=north] {\\$\\\\vect x\\$};\"),\n",
    "% \t\t Plots.Command(\"\\\\node at (axis cs:-0.5,0.75) [anchor=south west] {\\$\\\\vect s\\$};\"),\n",
    "% \t\t],\n",
    "% \t          style=\"width=9cm, axis lines=left, xtick=\\\\empty, ytick=\\\\empty, xmin=-2, xmax=1.75, ymin=-0.5, ymax=2, view={0}{90}, contour/labels=false, contour/draw color={black}\"))\n",
    "% \\end{jlcode}\n",
    "% \\begin{marginfigure}[-9cm]\n",
    "% \t\\centering\n",
    "% \t\\plot{fig_derivatives_2d_derivative}\n",
    "% \t\\caption{The two-dimensional directional derivative}\n",
    "% \\end{marginfigure}\n",
    "\n",
    "\\begin{example}\n",
    "\tWe wish to compute the directional derivative of $f(\\vect x) = x_1 x_2$ at $\\vect x = [1,0]$ in the direction $\\vect s = [-1, -1]$:\n",
    "\t\\begin{align*}\n",
    "\t\t\\nabla f(\\vect x) & = \\begin{bmatrix}\\frac{\\partial f}{\\partial x_1}, & \\frac{\\partial f}{\\partial x_2}\\end{bmatrix} = \\brock*{x_2, x_1} \\\\\n",
    "\t\t\\nabla_{\\vect s}f(\\vect x) & = \\nabla f(\\vect x)^\\transpose \\vect s = \\begin{bmatrix}0 & 1\\end{bmatrix} \\begin{bmatrix}-1 \\\\ -1\\end{bmatrix} = -1\n",
    "\t\\end{align*}\n",
    "\tWe can also compute the directional derivative as follows:\n",
    "\t\\begin{align*}\n",
    "\t\tg(\\alpha) &= f(\\vect x + \\alpha \\vect s) = (1 - \\alpha)(-\\alpha) = \\alpha^2 - \\alpha \\\\\n",
    "\t\tg'(\\alpha) &= 2\\alpha - 1 \\\\\n",
    "\t\tg'(0) &= -1\n",
    "\t\\end{align*}\n",
    "\t\\caption{Computing a directional derivative.\\label{ex:directional-deriv}}\n",
    "\\end{example}\n",
    "\\begin{marginfigure}[-8cm]\n",
    "\t\\centering\n",
    "\t\\begin{tikzpicture}\n",
    "\t\\begin{axis}[\n",
    "\twidth=\\marginparwidth,\n",
    "\txlabel={$\\alpha$},\n",
    "\tylabel={$f(\\alpha)$},\n",
    "\taxis lines=middle, xtick={1}, ytick=\\empty,\n",
    "\t]\n",
    "\t\\addplot[mark=none, samples=100, domain=-0.25:1.25]{(1-x)*x};\n",
    "\t\\addplot[mark=none, samples=2, domain=-0.25:0.5]{x};\n",
    "\t\\end{axis}\n",
    "\t\\end{tikzpicture}\n",
    "\\end{marginfigure}\n",
    "\n",
    "\n",
    "\n",
    "\\begin{ignore}\n",
    "\tFor multi-input, multi-output functions\\sidenote{$\\vect f: \\mathbb{R}^n \\rightarrow \\mathbb{R}^m$}, we can compute the \\vocab{Jacobian}:\n",
    "\n",
    "\t\\begin{equation}\n",
    "\t\\nabla \\vect f(\\vect x) = \\begin{bmatrix}\n",
    "\t\\frac{\\partial f_1(\\vect x)}{\\partial x_1} & \\frac{\\partial f_1(\\vect x)}{\\partial x_2} & \\cdots & \\frac{\\partial f_1(\\vect x)}{\\partial x_n} \\\\\n",
    "\t& \\vdots & & \\\\\n",
    "\t\\frac{\\partial f_m(\\vect x)}{\\partial x_1} & \\frac{\\partial f_m(\\vect x)}{\\partial x_2} & \\cdots & \\frac{\\partial f_m(\\vect x)}{\\partial x_n} \\\\\n",
    "\t\\end{bmatrix}\n",
    "\t\\end{equation}\n",
    "\n",
    "\tThe Jacobian forms the best linear approximation for $\\vect f$ at $\\vect x$, because it contains all of the gradient components.\n",
    "\t\\tim{I don't think we use the Jacobian in the text. Should we leave it out?} \\mykel{Let's comment it out. Let's try to make the book as simple as possible (but no simpler). If a topic is not connected with anything, we should strongly consider dropping it. Also, the Jacobian gets weird when we have the convention where the gradient is a column vector.}\n",
    "\\end{ignore}\n",
    "\n",
    "\\section{Numerical Differentiation}\n",
    "\n",
    "The process of estimating derivatives numerically is referred to as \\vocab{numerical differentiation}.\n",
    "Estimates can be derived in different ways from function evaluations. This section discusses finite difference methods and the complex step method.\\sidenote{For a more comprehensive treatment of the topics discussed in the remainder of this chapter, see \\fullcite{Griewank2008}.}\n",
    "\n",
    "\n",
    "% The gradient $\\nabla f$ or the derivative along a search direction ${\\frac{\\partial}{\\partial \\alpha} f(\\vect x + \\alpha \\vect d)}$ are often not readily avaiable for direct use, but they are critical to many optimization methods.\n",
    "% Given a method for derivative estimation one also automatically has a method for gradient estimation - one need merely estimate the partial derivative along each component direction.\n",
    "\n",
    "\\subsection{Finite Difference Methods}\n",
    "\n",
    "As the name implies, \\vocab{finite difference methods} compute the difference between two values that differ by a finite step size.\n",
    "They approximate the derivative definitions in \\cref{eq:derivative_defintions} using small differences:\n",
    "\\begin{fullwidth}\n",
    "\\begin{equation}\n",
    "\tf'(x) \\approx \\underbrace{\\frac{f(x + h) - f(x)}{h}}_\\text{forward difference} \\approx \\underbrace{\\frac{f(x + h/2) - f(x - h/2)}{h}}_\\text{central difference} \\approx \\underbrace{\\frac{f(x) - f(x - h)}{h}}_\\text{backward difference}\n",
    "\\end{equation}\n",
    "\\end{fullwidth}\n",
    "\n",
    "Mathematically, the smaller the step size $h$, the better the derivative estimate.\n",
    "Practically, values for $h$ that are too small can result in numerical cancellation errors.\n",
    "This effect is shown later in \\cref{fig:derivatives_derivative_error}. \\Cref{alg:finite-diff} provides implementations for these methods.\n",
    "\n",
    "\\begin{algorithm}\n",
    "\\begin{juliaverbatim}\n",
    "diff_forward(f, x; h=sqrt(eps(Float64))) = (f(x+h) - f(x))/h\n",
    "diff_central(f, x; h=cbrt(eps(Float64))) = (f(x+h/2) - f(x-h/2))/h\n",
    "diff_backward(f, x; h=sqrt(eps(Float64))) = (f(x) - f(x-h))/h\n",
    "\\end{juliaverbatim}\n",
    "\t\\caption{\n",
    "\t\t\\label{alg:finite-diff}\n",
    "\t\tFinite difference methods for estimating the derivative of a function \\jlv{f} at \\jlv{x} with finite difference \\jlv{h}.\n",
    "\t\tThe default step sizes are the square root or cube root of the machine precision for floating point values.\n",
    "\t\tThese step sizes balance machine round-off error with step size error.\n",
    "\n",
    "\t\tThe \\jlv{eps} function provides the step size between \\jlv{1.0} and the next larger representable floating-point value.\n",
    "\t}\n",
    "\\end{algorithm}\n",
    "\n",
    "\n",
    "The finite difference methods can be derived using the Taylor expansion.\n",
    "We will derive the forward difference derivative estimate, beginning with the Taylor expansion of $f$ about $x$:\n",
    "\\begin{equation}\n",
    "\tf(x + h) = f(x) + \\frac{f'(x)}{1!} h + \\frac{f''(x)}{2!} h^2 + \\frac{f'''(x)}{3!} h^3 + \\cdots\n",
    "\\end{equation}\n",
    "\n",
    "We can rearrange and solve for the first derivative:\n",
    "\\begin{align}\n",
    "\t\tf'(x) h & = f(x+h) - f(x) - \\frac{f''(x)}{2!} h^2 - \\frac{f'''(x)}{3!} h^3 - \\cdots \\\\\n",
    "\t\tf'(x) & = \\frac{f(x+h) - f(x)}{h} - \\frac{f''(x)}{2!} h - \\frac{f'''(x)}{3!} h^2 - \\cdots \\\\\n",
    "\t\tf'(x) & \\approx \\frac{f(x+h) - f(x)}{h}\n",
    "\\end{align}\n",
    "\n",
    "The forward difference approximates the true derivative for small $h$ with error dependent on $\\frac{f''(x)}{2!} h + \\frac{f'''(x)}{3!} h^2 + \\cdots$.\n",
    "This error term is $O(h)$, meaning the forward difference has linear error as $h$ approaches zero.\\sidenote[][]{\n",
    "\tAsymptotic notation is covered in \\chref{the appendix}{ch:math_concepts}.\n",
    "}\n",
    "\n",
    "The central difference method has an error term of $O(h^2)$.\\cite{Mathews2004}\n",
    "We can derive this error term using the Taylor expansion.\n",
    "The Taylor expansions about $x$ for $f(x+h/2)$ and $f(x-h/2)$ are:\n",
    "\\begin{align}\n",
    "\t\tf(x+h/2) & = f(x) + f'(x) \\frac{h}{2} + \\frac{f''(x)}{2!} \\paren*{\\frac{h}{2}}^2 + \\frac{f'''(x)}{3!}\\paren*{\\frac{h}{2}}^3 + \\cdots \\\\\n",
    "\t\tf(x-h/2) & = f(x) - f'(x) \\frac{h}{2} + \\frac{f''(x)}{2!} \\paren*{\\frac{h}{2}}^2 - \\frac{f'''(x)}{3!}\\paren*{\\frac{h}{2}}^3 + \\cdots\n",
    "\\end{align}\n",
    "\n",
    "Subtracting these expansions produces:\n",
    "\\begin{equation}\n",
    "\tf(x+h/2) - f(x-h/2) \\approx 2f'(x)\\frac{h}{2} + \\frac{2}{3!} f'''(x) \\paren*{\\frac{h}{2}}^3\n",
    "\\end{equation}\n",
    "\n",
    "We rearrange to obtain:\n",
    "\\begin{equation}\n",
    "\tf'(x) \\approx \\frac{f(x+h/2) - f(x-h/2)}{h} - \\frac{f'''(x)h^2}{24}\n",
    "\\end{equation}\n",
    "which shows that the approximation has quadratic error.\n",
    "\n",
    "\n",
    "\\subsection{Complex Step Method}\n",
    "\n",
    "We often run into the problem of needing to choose a step size $h$ small enough to provide a good approximation but not too small so as to lead to numerical subtractive cancellation issues.\n",
    "The \\vocab{complex step method} bypasses the effect of subtractive cancellation by using a single function evaluation. We evaluate the function once after taking a step in the imaginary direction.\\sidenote{\\fullcite{Martins2003}. Special care must be taken to ensure that the implementation of $f$ properly supports complex numbers as input.}\n",
    "\n",
    "The Taylor expansion for an imaginary step is:\n",
    "\\begin{equation}\n",
    "\tf(x + i h) = f(x) + i h f'(x) - h^2 \\frac{f''(x)}{2!} - i h^3 \\frac{f'''(x)}{3!} + \\cdots\n",
    "\\end{equation}\n",
    "Taking only the imaginary component of each side produces a derivative approximation:\n",
    "\\begin{align}\n",
    "\t\t\\Imag(f(x + i h)) & = h f'(x)  - h^3 \\frac{f'''(x)}{3!} + \\cdots \\\\\n",
    "\t\t\\Rightarrow f'(x) & = \\frac{\\Imag (f(x + i h))}{h} + h^2 \\frac{f'''(x)}{3!} - \\cdots \\\\\n",
    "\t\t\t\t\t\t  & = \\frac{\\Imag (f(x + i h))}{h} + O(h^2) \\text{ as } h \\rightarrow 0\n",
    "\\end{align}\n",
    "An implementation is provided by \\cref{alg:complex-step}. The real part approximates $f(x)$ to within $O(h^2)$ as $h \\rightarrow 0$:\n",
    "\\begin{align}\n",
    "\t\t\\Real(f(x + i h)) = f(x) - h^2 \\frac{f''(x)}{2!} + \\hdots \\\\\n",
    "\t\t\\Rightarrow f(x) = \\Real(f(x + i h)) + h^2 \\frac{f''(x)}{2!} - \\cdots\n",
    "\\end{align}\n",
    "\n",
    "Thus, we can evaluate both $f(x)$ and $f'(x)$ using a single evaluation of $f$ with complex arguments. \\Cref{ex:complex-step} shows the calculations involved for estimating the derivative of a function at a particular point.\n",
    "\\Cref{alg:complex-step} implements the complex step method.\n",
    "\\Cref{fig:derivatives_derivative_error} compares the numerical error of the complex step method to the forward and central difference methods as the step size is varied.\n",
    "\n",
    "\n",
    "\\begin{algorithm}\n",
    "\\begin{juliaverbatim}\n",
    "diff_complex(f, x; h=1e-20) = imag(f(x + h*im)) / h\n",
    "\\end{juliaverbatim}\n",
    "\\caption{\n",
    "\tThe complex step method for estimating the derivative of a function \\jlv{f} at \\jlv{x} with finite difference \\jlv{h}.\n",
    "\t\\label{alg:complex-step}\n",
    "}\n",
    "\\end{algorithm}\n",
    "\n",
    "%Use of the complex step method requires that the function be able to accept complex values as input.\n",
    "%Existing packages sometimes require modifying the code that implements $f$, which can limit the method's direct application.\n",
    "\n",
    "\\begin{example}\n",
    "Consider $f(x) = \\sin(x^2)$.\n",
    "The function value at $x = \\pi/2$ is approximately \\num{0.624266} and the derivative is $\\pi \\cos(\\pi^2/4) \\approx -2.45425$.\n",
    "We can arrive at this using the complex step method:\n",
    "\\begin{juliaconsole}\n",
    "f = x -> sin(x^2);\n",
    "v = f(π/2 + 0.001im);\n",
    "real(v) # f(x)\n",
    "imag(v)/0.001 # f'(x)\n",
    "\\end{juliaconsole}\n",
    "\\caption{The complex step method for estimating derivatives.\\label{ex:complex-step}}\n",
    "\\end{example}\n",
    "\n",
    "\\begin{jlcode}\n",
    "\tdiff_forward(f, x; h=sqrt(eps(Float64))) = (f(x+h) - f(x))/h\n",
    "\tdiff_central(f, x; h=cbrt(eps(Float64))) = (f(x+h/2) - f(x-h/2))/h\n",
    "\tdiff_backward(f, x; h=sqrt(eps(Float64))) = (f(x) - f(x-h))/h\n",
    "\tdiff_complex(f, x; h=1e-20) = imag(f(x + h*im)) / h\n",
    "\n",
    "\tp = let\n",
    "\t\tx = 0.5\n",
    "\t\td_true = 1.0\n",
    "\n",
    "\t\tarr_h = collect(10 .^ range(-17, stop=1, length=101))\n",
    "\t\tarr_forward = [abs(d_true - diff_forward(sin, x, h=h)/cos(x))\n",
    "\t\t\t\t\t\tfor h in arr_h]\n",
    "\t\tarr_central = [abs(d_true - diff_central(sin, x, h=h)/cos(x))\n",
    "\t\t\t\t\t\tfor h in arr_h]\n",
    "\t\tarr_complex = [abs(d_true - diff_complex(sin, x, h=h)/cos(x))\n",
    "\t\t\t\t\t\tfor h in arr_h]\n",
    "\n",
    "\t\tAxis([\n",
    "\t\t\tPlots.Linear(arr_h, arr_complex, style=\"thick, solid,\n",
    "\t\t\t\tpastelGreen, mark=none\", legendentry=\"complex\"),\n",
    "\t\t\tPlots.Linear(arr_h, arr_forward, style=\"thick, solid,\n",
    "\t\t\t\tpastelBlue, mark=none\", legendentry=\"forward\"),\n",
    "\t\t\tPlots.Linear(arr_h, arr_central, style=\"thick, solid,\n",
    "\t\t\t\tpastelRed, mark=none\", legendentry=\"central\")],\n",
    "\t\t\txmode=\"log\", ymode=\"log\",\n",
    "\t\t\txlabel=\"step size \\$h\\$\",\n",
    "\t\t\tylabel=\"absolute relative error\",\n",
    "\t\t\twidth=\"10cm\", height=\"8cm\",\n",
    "\t\t\tstyle=\"legend pos=outer north east\"\n",
    "\t\t)\n",
    "\tend\n",
    "\tplot(p)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}