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

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.

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

My current solution is as follows:

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

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

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

Thanks. The following code:

``````using ForwardDiff, DiffResults, StaticArrays, BenchmarkTools

function computehessian(f, x::AbstractArray)
result = DiffResults.HessianResult(x)
result = ForwardDiff.hessian!(result, f, x)
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. 2 Likes

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(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)
``````

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.