DimensionMismatch when taking hessian of vector function using DiffResults

[jacobmogic]

Hessian of vector function can be calculated by composing two jacobians:

f(x) = x.^2
g(x) = ForwardDiff.jacobian(f,x)
x = rand(3)
ForwardDiff.jacobian(g,x)

But if I want to retrieve the first jacobian using DiffResults,

result = DiffResults.JacobianResult(x)
result = ForwardDiff.jacobian!(result,g,x)

This would give

ERROR: DimensionMismatch("new dimensions (9, 3) must be consistent with array size 9")

Any workaround?

[lliill]

I have the same problem, waiting for an answer

[Jollywatt]

The reason is because DiffResults.JacobianResult(x) allocates a square result matrix, as the single-argument form assumes the Jacobian will be square. That is not the case for your function g, which has 3 inputs and 3×3 = 9 outputs, leading to a 9×3 Jacobian.

The 2-argument method DiffResults.JacobianResult(x, y) should be used when the Jacobian isn’t necessarily square. This works:

result = DiffResults.JacobianResult(x, g(x))
result = ForwardDiff.jacobian!(result, g, x)

Although it defeats the purpose of not evaluating g(x) first.

If you know the size that the Jacobian should be, you can preallocate the array with the correct size yourself:

result = Matrix{Float64}(undef, 3, 9)
result = ForwardDiff.jacobian!(result, g, x)

It’s annoying you can’t magically do figure out the resulting size automatically, though.