Hessian matrix of log-likelhood function that uses structure-of-array var_value<VectorXd>

General
1

I have a log-likelihood function that uses not only stan::math::var’s, but also involves stan::math::var_value <<Eigen::VectorXd> > (and var_value < MatrixXd > as well) since it involves pixel operations with thousands of pixels. I would like to be able to calculate the Hessian of the log-likelihood (i.e. the Fisher matrix), but it seems this requires templating the log-likelihood for forward-mode autodiff as well as reverse-mode. Is there a forward-mode version of var_value<...> structures? Or if not, are there plans to implement such a thing? I don’t want to switch to using vectors of var/fvar because the autodiff is considerably slower using that approach. Any ideas would be appreciated here!

2

You are correct that currently you do need forward mode, and that this is not really set up to do struct of arrays work.

Do you need full hessians (which will be quite expensive no matter how you compute them) or can your use case get away with hessian-vector products? There are ways to approximate the latter using only a couple gradient evaluations

3

I need the full Hessian so that I can invert it to obtain the parameter covariance matrix. Is the best strategy to generalize my templating so that if it’s instantiated with fvar, it uses array of structs instead of struct of arrays? It’s a pain, but I think I could pull that off.

Is there a reason why there is no forward-mode struct of arrays, or is it just lack of demand? Just curious. Thanks for your input!

4

Are you writing your own analytic gradient? If so, it would be necessary to write a separate forward mode that also calculates the dual for forward mode

@stevebronder is the best person to ask about the history of fvar/SoA, but it probably is just a lack of demand, not a technical issue

You may also want to try math/stan/math/rev/functor/finite_diff_hessian_auto.hpp at develop · stan-dev/math · GitHub first, just to make sure you can achieve what you want when you do have the Hessian

5

At this time there is no fvar<var_value<Eigen::Matrix>> and fvar<Eigen::Matrix>. It is mostly that we did not want to rewrite the forward mode library for it. We would be open to a contributions for that work though!