Nth (univariate) derivative using AD

General Usage
,
1

I need an n th derivative of a univariate function for testing purposes. I came up with

import ForwardDiff

function ddn(f, x, ::Val{N}) where {N}
    fn = foldl((f, _) -> x -> ForwardDiff.derivative(f, x),
               ntuple(_ -> nothing, Val(N)); init = f)
    fn(x)
end

@code_warntype ddn(sin, 0.0, Val(3)) # type stable

but I am interested in other solutions.

2 Likes
2

I’ve seen people use GitHub - JuliaDiff/TaylorSeries.jl: Taylor polynomial expansions in one and several independent variables. for higher-order scalar differentiation. I can’t tell you if this is the recommended approach today though.

3 Likes
3

I thought Diffractor was meant to be especially suited to higher-order derivatives? I’m not really an AD user myself but based on some of the discussions I’ve seen on Slack it sounds like people have got it to work on quite a few test cases at this point, although afaik it’s still rather experimental.

1 Like
4

Thanks for all the replies. I went with the original code above, finding it the simplest — I only need it for unit testing.

Incidentally, I found that implementing a quick & dirty AD system in <50 LOC is the best solution I could come up with for calculating higher-order derivatives of Chebyshev polynomials and their Smolyak combinations. It only supports +, -, * for the polynomials, and / in addition for transformations.

After many years of Julia, I still find it unbelievable how powerful the language is.

6 Likes
5

Wow. That’s a very neat use of foldl. It’s even faster than using TaylorSeries like @baggepinnen mentioned:

using BenchmarkTools, ForwardDiff, TaylorSeries

function ddn(f, x, ::Val{N}) where {N}
  fn = foldl((f, _) -> x -> ForwardDiff.derivative(f, x),
             ntuple(_ -> nothing, Val(N)); init = f)
  fn(x)
end

# also type stable:
ddn_ts(f, x, N) = f(x+Taylor1(Float64, N))[N]*factorial(N)

# making const as well to be extra paranoid about type stuff.
const x = 1.1
const v = 5
const val_v = Val(v)

# simple f function, but trying to avoid things that are literally in the lookup
# tables for ForwardDiff:
f(x) = exp(cos(log(x+1.1)+1.1)+1.1)

@btime ddn(f, $x, $val_v); # 329 ns, 7  alloc: 336 bytes
@btime ddn_ts(f, $x, $v);  # 389 ns  13 alloc: 1.36 KiB

Also, wow again about ForwardDiff being so hard to beat.

2 Likes
6

OT: [not using AD]

If you want you can calculate the derivatives with Symbolics

using Symbolics

function der_symbolics(f, a::Int) # f is a callable julia function
    x = Symbolics.variable("x")
    der = (Differential(x)^a)(f(x)) |> expand_derivatives
    eval(build_function(der, x))   # this is a callable julia function
end

sin_third_deriv = der_symbolics(sin, 3) 
sin_third_deriv(0.0)

and sin_third_deriv is compatible with AD like ForwardDiff.
der_symbolics is not fast, however sin_third_deriv is fast.