24 Computing derivatives in Julia

This section uses these add-on packages:

using CalculusWithJulia
using Plots
plotly()
using ForwardDiff
using SymPy
using Roots

SymPy returns symbolic derivatives. Up to choices of simplification, these answers match those that would be derived by hand. This is useful when comparing with known answers and for seeing the structure of the answer. However, there are times we just want to work with the answer numerically. For that we have other options within Julia. We discuss approximate derivatives and automatic derivatives in this section. The latter will find wide usage in these notes.

24.1 Approximate derivatives

By approximating the limit of the secant line with a value for a small, but positive, \(h\), we get an approximation to the derivative. That is

\[ f'(x) \approx \frac{f(x+h) - f(x)}{h}. \]

This is the forward-difference approximation. The central difference approximation looks to both sides of \(x\):

\[ f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}. \]

Though in general they are different, both are easy and fast to compute, useful approximations to the derivative. The central difference is usually more accurate for the same size \(h\). However, both are susceptible to round-off errors. The numerator is a subtraction of like-size numbers—a perfect opportunity to lose precision.

Due to numeric issues there is a balancing act:

if \(h\) is too big the approximation to the limit is not good.
if \(h\) is too small the round-off errors are problematic,

For the forward difference \(h\) values around \(10^{-8}\) are typically good, for the central difference, values around \(10^{-6}\) are typically good, but these ranges aren’t always the case.

Example

Let’s verify that the forward difference isn’t too far off.

f(x) = exp(-x^2/2)
c = 1
h = 1e-8
fapprox = (f(c+h) - f(c)) / h

-0.6065306479285937

We can compare to the actual with:

@syms x
df = diff(f(x), x)
factual = convert(Float64, df(c))
abs(factual - fapprox)

1.1784039744000552e-8

The error is about \(1\) part in \(100\) million.

The central difference is better here, even with a bigger \(h\):¹

h = 1e-6
cdapprox = (f(c+h) - f(c-h)) / (2h)
abs(factual - cdapprox)

1.5675682973892435e-11

24.1.1 Automatic derivatives

Roughly speaking, symbolic derivatives are exact, but can be slow to compute, whereas approximate derivatives are not exact, but fast to compute. However, the widely used automatic derivatives are both exact and fast to compute.

Automatic differentiation (AD) is the general name for a few different approaches. We utilize forward mode automatic differentiation. The ForwardDiff package provides one of several ways for Julia to compute automatic derivatives. ForwardDiff is well suited for functions encountered in these notes, which depend on at most a few variables and output no more than a few values at once.

The ForwardDiff package was loaded in this section; in general its features are available when the CalculusWithJulia package is loaded, as that package provides a more convenient interface. The derivative function is not exported by ForwardDiff, so its usage requires qualification. To illustrate, to find the derivative of \(f(x)\) at a point we have this syntax for ForwardDiff:

ForwardDiff.derivative(f, c)   # derivative is qualified by a module name

-0.6065306597126334

The CalculusWithJulia package defines a method for ' (a postfix operator in Julia) so that when applied to a Function object a function for computing values of the derivative, as above, is returned.²

To be clear, this usage returns a function that computes the derivative of f:

f'

#D##0 (generic function with 1 method)

And this usage calls the derivative of f at the value stored by c:

f'(c)

-0.6065306597126334

This is the same mathematical notation that is commonly used.

Here we see the error in estimating \(f'(1)\):

abs(factual - f'(c))

0.0

In this case, the automatic derivative is exact.

Example

For \(f(x) = \sqrt{1 + \sin(\cos(x))}\) compare the difference between the forward derivative with \(h=1e-8\) and that computed by automatic differentiation at \(x=\pi/4\).

The forward derivative is found with:

f(x) = sqrt(1 + sin(cos(x)))
c, h = pi/4, 1e-8
fwd = (f(c+h) - f(c))/h

-0.20927346522370271

That given by automatic differentiation is:

ds_value = f'(c)
ds_value, fwd, ds_value - fwd

(-0.20927346371432803, -0.20927346522370271, 1.5093746807970376e-9)

Finally, SymPy gives an exact value we use to compare:

fp = diff(f(x), x)

\(- \frac{\sin{\left(x \right)} \cos{\left(\cos{\left(x \right)} \right)}}{2 \sqrt{\sin{\left(\cos{\left(x \right)} \right)} + 1}}\)

actual = float(fp(PI/4))
actual - ds_value, actual - fwd

(-5.551115123125783e-17, 1.5093746252858864e-9)

As expected, the automatic derivative is nearly exact and accurate up to possibly accumulated floating point differences; the forward difference is just pretty close.

Example

Suppose our task is to find a zero of the second derivative of \(k(x) = e^{-x^2/2}\) in \([0, 10]\), a known bracket. The second derivative is found by k'' below, with the derivative operation being applied twice behind the scenes. With this, we have:

k(x) = exp(-x^2/2)
find_zero(k'', (0, 10))

1.0

As with plotting and other uses of functions as arguments, we pass in the function object, k'', and not the function evaluated at a point.

24.2 Recap on derivatives in Julia

A quick summary of the \(3\) different ways for finding derivatives in Julia presented in these notes:

Symbolic derivatives are found using diff from SymPy
Automatic derivatives are found using the notation f' which utilizes ForwardDiff.derivative
approximate derivatives at a point, c, for a given h are found with (f(c+h)-f(c))/h.

For example, here all three are computed and compared:

f(x) = exp(-x) * sin(x)

c = pi
h = 1e-8

fp = diff(f(x),x)

Dict(:approximate=>(f(c+h) - f(c))/h, :automatic=>f'(c),
     :symbolic=>fp, :symbolic_evaluated => fp(x=>c))

Dict{Symbol, Number} with 4 entries:
  :approximate        => -0.0432139
  :automatic          => -0.0432139
  :symbolic           => -exp(-x)*sin(x) + exp(-x)*cos(x)
  :symbolic_evaluated => -exp(-pi)

24.3 Questions

Question

Find the derivative using a forward difference approximation of \(f(x) = x^x\) at the point \(x=2\) using h=0.1:

Use f' find the value using automatic differentiation

Question

Let \(f(x) = x^x\). Using automatic differentiation, find \(f'(3)\).

Question

Let \(f(x) = \lvert 1 - \sqrt{1 + x}\rvert\). Using automatic differentation, find \(f'(3)\).

Question

Let \(f(x) = e^{\sin(x)}\). Using automatic differentation, find \(f'(3)\).

Question

For Julia’s airyai function find a numeric derivative using the forward difference. For \(c=3\) and \(h=10^{-8}\) find the forward difference approximation to \(f'(3)\) for the airyai function.

Question

Find the rate of change with respect to time of the function \(f(t)= 64 - 16t^2\) at \(t=1\).

Question

Find the rate of change with respect to height, \(h\), of \(f(h) = 32h^3 - 62 h + 12\) at \(h=2\).

Question

Mathematically, as the value of h in the forward difference gets smaller the forward difference approximation gets better. On the computer, this is thwarted by floating point representation issues (in particular the error in subtracting two like-sized numbers in forming \(f(x+h)-f(x)\).)

For 1e-16 what is the error (in absolute value) in finding the forward difference approximation for the derivative of \(\sin(x)\) at \(x=0\)?

Repeat for \(x=\pi/4\):

Question

Let \(f(x) = e^{-x^2/2}\). At \(c=2\) we can compare the exact answer for the derivative to the forward difference approximation for different values of h. For example:

f(x) = exp(-x^2/2)
c = 2
@syms x
exact = float(diff(f(x), x)(x=>2))

-0.2706705664732254

Whereas,

hs = [1/10^i for i in 0:16]
fdiffs = [(f(c+h) - f(c))/h for h in hs]
error = fdiffs .- exact
[hs error]

17×2 Matrix{Float64}:
 1.0       0.146444
 0.1       0.019823
 0.01      0.00202549
 0.001     0.000202958
 0.0001    2.02998e-5
 1.0e-5    2.03002e-6
 1.0e-6    2.02962e-7
 1.0e-7    2.07743e-8
 1.0e-8    2.73314e-9
 1.0e-9   -3.05735e-8
 1.0e-10  -1.41596e-7
 1.0e-11  -1.80693e-6
 1.0e-12  -2.95625e-5
 1.0e-13   5.37042e-5
 1.0e-14  -0.00688519
 1.0e-15   0.0208704
 1.0e-16   0.270671

Which value of h provided the smallest error? Write the exponent as i in 1/10^i.

We repeat with the central difference approximation.

hs = [1/10^i for i in 0:16]
cdiffs = [(f(c+h) - f(c-h))/(2h) for h in hs]
error = cdiffs .- exact
[hs error]

17×2 Matrix{Float64}:
 1.0      -0.0270403
 0.1      -0.00044909
 0.01     -4.51097e-6
 0.001    -4.51117e-8
 0.0001   -4.51394e-10
 1.0e-5   -6.33293e-12
 1.0e-6   -7.81819e-13
 1.0e-7    2.35141e-10
 1.0e-8    1.34536e-9
 1.0e-9   -1.66958e-8
 1.0e-10  -2.81797e-9
 1.0e-11  -4.19152e-7
 1.0e-12  -1.56847e-5
 1.0e-13   0.000192482
 1.0e-14  -0.00272185
 1.0e-15   0.0069926
 1.0e-16   0.270671

Which value of h provided the smallest error? Write the exponent as i in 1/10^i.

The FiniteDifferences and FiniteDiff packages provide performant interfaces for differentiation based on finite differences.↩︎
The CalculusWithJulia package defines a method for Base.adjoint(f::Function). In Julia speak this is a form of type piracy, as the package modifies a method (adjoint) for a type (Function) neither of which belongs to the package. As well, the meaning given does not conform with the expected generic meaning of the operation elsewhere in the Julia ecosystem. Admittedly this is bad form, but a useful deviation from good practice for pedagogical reasons.↩︎