Implementation of an implace Hessian of a vector-valued function?

HI, i’m looking to implement an implace hessian of a vector-valued function, for example:

f(x) = [x[1],x[2],sin(x[1]),sum(x)]

The ForwardDiff Documentation gives an out of place example (slightly wrong, as size(f(x)) !=size(x) in the general case:

julia> function vector_hessian(f, x)
       n = length(x)
       out = ForwardDiff.jacobian(x -> ForwardDiff.jacobian(f, x), x)
       return reshape(out, n, n, n)

the problem arrives when you use a inplace version of the same function:

function f!(res,x) 
res .= [x[1],x[2],sin(x[1]),sum(x)]

the problem with nested inplace jacobians is that the result of the intermediate jacobian is a dual, and its very difficult to create a cache storage.
Anybody has encountered something like that? how can i create an efficient cache for this type of function?

See Nested ForwardDiff.jacobian calls with inplace function for one solution.


Although I didn’t use the same method, the idea is here, to pass a implace Jacobian to another implace Jacobian, the cache needs to be build dinamically, and the inner jacobian should be passed as an out of place cached function.

And for reference, FiniteDiff.jl has a fully cached version of Hessian computations through finite differencing:


Nice! I also need to implement this on FiniteDiff (as an user selected option)

But the Hessian computations in FiniteDiff.jl do no support vector-valued functions, right?