# Hessian of a vector-valued in-place function with ForwardDiff.jl

**URL:** <https://discourse.julialang.org/t/hessian-of-a-vector-valued-in-place-function-with-forwarddiff-jl/134904>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [January 6, 2026, 7:32am UTC](https://discourse.julialang.org/t/hessian-of-a-vector-valued-in-place-function-with-forwarddiff-jl/134904 "2026-01-06T07:32:07Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Yuricst](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuricst/32/207613_2.png) [@Yuricst](https://discourse.julialang.org/u/Yuricst)\
**Post date:** [January 6, 2026, 7:32am UTC](https://discourse.julialang.org/t/hessian-of-a-vector-valued-in-place-function-with-forwarddiff-jl/134904/1 "2026-01-06T07:32:07Z")

</div>

I am looking to compute the Hessian of a vector-valued in-place function that looks like

```julia
function eom!(dx, x, p, t)
    dx[1:3] = x[4:6]
    dx[4:6] = -p[1] / norm(x[1:3])^3 * x[1:3]
end

```

using `ForwardDiff.jl` but I am struggling to find the correct way to implement a `vector_hessian()` function, as suggested by the [ForwardDiff.jl tutorial page](https://juliadiff.org/ForwardDiff.jl/stable/user/advanced/#Hessian-of-a-vector-valued-function):

> you could write a version of `vector_hessian` which supports functions of the form `f!(y, x)` , or perhaps an in-place Jacobian with `ForwardDiff.jacobian!` .

So far, to compute the Jacobian & Hessian, I have something like

```julia
using LinearAlgebra
using ForwardDiff

p = [1.0]
function eom!(y, x, p, t)
    y[1:3] = x[4:6]
    y[4:6] = -p[1] / norm(x[1:3])^3 * x[1:3]
    return
end

t0 = 0.0
x0 = [1.0, 0.1, 0.2, 0.3, 1.0, -0.1]
nx = length(x0)

# Jacobian (works)
A1 = ForwardDiff.jacobian((y,x) -> eom!(y,x,p,t0), zeros(nx), x0)

# Hessian (doesn't work)
A2 = ForwardDiff.jacobian(x -> ForwardDiff.jacobian((y,x) -> eom!(y,x,p,t0), zeros(nx), x), x0)
A2 = reshape(A2, nx, nx, nx)

```

where I get the error with the second-to-last line

```julia-auto
ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#2124#2125", Float64}, Float64, 6})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

```

As I understand, the issue is with how I am passing the dummy output storage `y`, where specifying a `Float64` vector (as I am doing via `zeros(nx)`) works fine when calling `ForwardDiff.jacobian` once (as in when calculate `A1`), but it is not when I have it nested (as in when I calculate `A2`).

I’d appreciate any help/suggestion, including any tips to implement both the Jacobian & Hessian calculation in this context in the most “appropriate” way.

Thank you!

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [January 6, 2026, 8:04am UTC](https://discourse.julialang.org/t/hessian-of-a-vector-valued-in-place-function-with-forwarddiff-jl/134904/2 "2026-01-06T08:04:50Z")

</div>

To better understand your needs:

- Is performance important (is this routine a bottleneck in your code)?
- How big are the input and output vectors?

---

<div class="post-metadata">

**Author:** ![Yuricst](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuricst/32/207613_2.png) [@Yuricst](https://discourse.julialang.org/u/Yuricst)\
**Post date:** [January 6, 2026, 8:13am UTC](https://discourse.julialang.org/t/hessian-of-a-vector-valued-in-place-function-with-forwarddiff-jl/134904/3 "2026-01-06T08:13:07Z")

</div>

Hello @gdalle,

- Performance is actually important since I would make many calls to evaluate the Jacobian/Hessian, but I am not yet at the point where this is a critical bottleneck - I would be interested to know if there is both a “naive” fix and a more elaborate, better-performing approach, if possible (like leveraging the Hessian’s symmetry?).
- Input `x` and output `y` at the moment are length-6 vectors, so I am looking at a `6*6*6` Hessian. I would be potentially looking at increasing up to e.g. length-10 ~ length-100 vectors (wondering if on the higher end here, what I should use might change?)

Thanks!