Higher dimensionnal Hessians

[Magalame]

Hi!

I was wondering if there was a practical way to produce an array (of size m ) of n*n Hessians for a function f: R^n -> R^m.

I checked out vector_hessian from the docs, but it seems to give a weird output. For example, I'd get

[:, :, 1] =
 10  10
  2  -8

[:, :, 2] =
  0  -8
 -8   0

[:, :, 3] =
 -8  0
  4  4

Instead of

[:, :, 1] =
 10  -8
  -8  4

[:, :, 2] =
 2   0
 0   0

[:, :, 3] =
 10  -8
  -8  4

I suspect I'm just poorly using the reshape function: I pass out to reshape(n,n,m)

Thanks!

[Magalame]

I think I found a way to get what I want:

function vector_hessian(f, x)
       n = length(x)
       out = ForwardDiff.jacobian(x -> ForwardDiff.jacobian(f, x)', x)
       m = fld(length(out),n^2)
       return reshape(out',n,n,m)
end

Although I'm not absolutely certain it's the most efficient way to do it.

[longemen3000]

I have an (somewhat) efficient implementation of a vector_hessian, used by autodiff twiceDifferentiableConstraints, but for in place functions, waiting in a pull request to NLSolversBase. Hopefully in can be ported to TwiceDifferentiable.