Hi,
I am interested in calculating the hessian of a function that requires calls to multiple other built-in and local functions. What is otherwise a run-able function quickly becomes incompatible with the Hessian type. Here is a simple example illustrating the issue:

function TestFun(x::Vector)
Sigma = [x[1]^2 x[1] * x[2] * x[3]; x[1] * x[2] * x[3] x[2]^2];
lambda = eigvals(Sigma);
ldV = sum(log(lambda));
return ldV
end

Generate an arbitrary input vector:

W = rand(3)

Compute the hessian:

hess = ForwardDiff.hessian(TestFun, W)

Which produces the error,

LoadError: MethodError: eigvals! has no method matching eigvals!(::Array{ForwardDiff.HessianNumber{3,Float64,Tuple{Float64,Float64,Float64}},2})
while loading In[331], in expression starting on line 1

in eigvals at linalg/eigen.jl:92
in TestFun at In[329]:3
in _calc_hessian at C:\Users....julia\v0.4\ForwardDiff\src\api\hessian.jl:98
in hessian at C:\Users....julia\v0.4\ForwardDiff\src\api\hessian.jl:27

In this case, it is easy enough to avoid the call to eigvals by hardcoding the equations for the eigenvalues since Sigma is just a 2x2 matrix, but this is far from ideal for larger matrices and doesn't address the underlying difficulty of passing along Hessian type variables.

I am curious if anyone has thoughts on solutions/best practices. Any help is greatly appreciated.
Thanks!