to compute the second derivative but that sounds inefficient, thinking about the algebraic structure of the computation. I’d like to get the first derivative as well using the DiffResults API, and this looks like a mess to do in this format with nested calls.

Alternatively, I can wrap everything in a vector and use

ForwardDiff.hessian(x -> f(x[1]), [x])

but I am not sure about the impact of extra vector allocations.

Is there a better builtin method already implemented in the package? The source of DiffResults suggests that it supports higher derivatives natively, but I don’t see anything corresponding in ForwardDiff, so maybe the API is there but not the implementation.

If the scalar/vector distinction is causing issues, you can always use SVector(x) from the StaticArrays package to create a 1-element vector (or small-ish vectors, in general) with virtually no overhead. These play well with these autodiff packages.

julia> using DiffResults, ForwardDiff, StaticArrays
julia> x = SVector(1.0)
1-element SVector{1, Float64} with indices SOneTo(1):
1.0
julia> res = DiffResults.HessianResult(x)
ImmutableDiffResult(1.0, ([1.0], [0.0;;]))
julia> ForwardDiff.hessian!(res,z->sin(only(z)),x)
ImmutableDiffResult(0.8414709848078965, ([0.5403023058681398], [-0.8414709848078965;;]))

Note that I used sin(only(z)) instead of sin(z) to extract the one element from the vector input.

Many thanks; this works and it seems like a good solution!

In the meantime I have also found an open issue on Github from 2017 for an API to compute higher scalar derivatives, so I assume there is still nothing implemented natively in ForwardDiff.