Computing value and first & second derivatives of a scalar function with scalar input, using ForwardDiff

Specific Domains Numerics
,
1

I cannot find, either in the documentation of ForwardDiff, or on the web more widely, any information on how to compute the second derivative of a scalar function with a scalar input. How can I do this? (Given a function f that takes x<:Real as input).

Is there a way to compute the function output and first and second derivatives simultaneously, like calling hessian with a DiffResult? Note that this latter approach only works with array-like inputs.

2

I think just differentiating twice would work. Although you probably wanna take a look at GitHub - JuliaDiff/TaylorDiff.jl: Taylor-mode automatic differentiation for higher-order derivatives too

1 Like
3

My current solution is as follows:

function computehessian(f, x::AbstractArray)
    result = DiffResults.HessianResult(x)
    result = ForwardDiff.hessian!(result, f, x)
    return DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
end

computehessian(f, x::T) where T <: Number = map(first, computehessian(f ∘ first, StaticArrays.SVector{1, T}(x)))
4

You can differentiate twice:

df(x)  = ForwardDiff.derivative(f,x)
d2f(x) = ForwardDiff.derivative(df,x)

sought_result = (f(x), df(x), d2f(x))
1 Like
5

Thanks. The following code:

using ForwardDiff, DiffResults, StaticArrays, BenchmarkTools

function computehessian(f, x::AbstractArray)
    result = DiffResults.HessianResult(x)
    result = ForwardDiff.hessian!(result, f, x)
    return DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
end
computehessian(f, x::T) where T <: Number = map(first, computehessian(f ∘ first, SVector{1, T}(x)))

function computehessian2(f, x)
    df(x)  = ForwardDiff.derivative(f, x)
    d2f(x) = ForwardDiff.derivative(df, x)
    return (f(x), df(x), d2f(x))
end

x = rand()
@btime computehessian(exp, $x)
@btime computehessian2(exp, $x)

produces

  50.279 ns (1 allocation: 64 bytes)
  27.513 ns (0 allocations: 0 bytes)

so it seems your approach is better. :smiley:

2 Likes
6

However, with a more complicated f function it isn’t so clear cut:

x = rand()
w = rand()
s1 = rand()
s2 = rand()
f(w, s1, s2, x) = log(w * s1 * exp(-0.5 * x * s1 ^ 2) + (1-w) * s2 * exp(-0.5 * x * s2 ^ 2))
expandfunc(args, v) = args[1](args[2:end]..., v)
fixallbutlast(func, args...) = Base.Fix1(expandfunc, (func, args...))
g = fixallbutlast(f, w, s1, s2)
@btime computehessian($g, $x)
@btime computehessian2($g, $x)

outputs:

  69.701 ns (1 allocation: 64 bytes)
  83.289 ns (0 allocations: 0 bytes)
7

Given this slow-down (caused by calling f 4 times instead of once), I think it strange that the hessian! functionality doesn’t exist for scalar x, so I opened an issue about it.