I think the solver is supposed to be an improvement on a Newton iteration (Algebra solver details?). It’s computing the Jacobian internally using autodiff I think.
Something that I think I’d like would be able to be able to program custom Jacobian’s so that the ODE solver n’ such didn’t have to compute them itself. edit: the Hessian of the likelihood term is the Jacobian of the optimization problem here (just to be clear on the words).
There are plenty of downsides to this though in terms of development/testing. Even if this could be exposed through a C++ hack though…