I have a function f(x,y) where x and y are both vectors, and I would like to calculate the matrix of mixed partials. The following works but is not type stable
using ForwardDiff
f(x,y) = sum(x)*sum(y)
x0 = [1.0, 2.0, 3.0]
y0 = [-1.0, -2.0, -3.0, -4.0]
fx = y -> ForwardDiff.gradient(x -> f(x, y), x0)
fxy = ForwardDiff.jacobian(fx, y0)
Normally, the way I make gradients and Jacobians type stable with ForwardDiff is to preallocate the output vector or matrix, but since the Jacobian computation needs to differentiate throught the gradient computation I cannot figure out the correct way to do it.
The function f above is just an example, in my actual code I need to get this to work for any function (within reason) of two vectors x and y. Regrettably, some instances will be too large for StaticArrays to be useful, which otherwise works fine for small systems.