Adding logdensity_and_hessian_vector_product to the API

[sethaxen]

Sometimes one only needs the action of the Hessian on a relatively small number of vectors (e.g. this paper, which uses Lanczos algorithm to estimate the largest eigenvalue of the Hessian using just Hessian-vector-products). In such cases, for high-dimensional models, computing the full Hessian can be very wasteful.

For these cases, it would be useful to add logdensity_and_hessian_vector_product to the API. Ideally one would like to be able to use capabilities to distinguish between a model that just implements logdensity_gradient_and_hessian and logdensity_and_hessian_vector_product.

Notes

• Bridgestan exposes a similar function: https://roualdes.github.io/bridgestan/latest/languages/julia.html#BridgeStan.log_density_hessian_vector_product
• DifferentiationInterface exposes hvp and gradient_and_hvp: https://juliadiff.org/DifferentiationInterface.jl/DifferentiationInterface/stable/api/#DifferentiationInterface.hvp
• Some optimizers supported by Optimization.jl can use HVPs: https://docs.sciml.ai/Optimization/stable/optimization_packages/optim/#Hessian-Free-Second-Order

[gdalle]

While we're at it, sparse Hessians could also be helpful in probabilistic programming I suppose?

[sethaxen]

I don't have experience with them, but I imagine they could be useful for memory-efficient Laplace approximations for high-dimensional models. I haven't seen this done yet though, so I don't know how useful it would be to add it to the API here.

FWIW, searching GitHub, I couldn't find any public code where someone is using logdensity_gradient_and_hessian, so any use of Hessian in the API is probably niche functionality.

[gdalle]

I agree that it's probably niche right now, but one of the reasons might be the lack of support for sparse Hessians in the ecosystem, until now?

[tpapp]

I am happy to add API for anything if someone actually needs it. If there is a project, however experimental, that needs this, please comment here or make a PR.

I am skeptical about the value of expanding the API just for the sake of completeness.

[sethaxen]

@tpapp great, thanks! I'm going to run some experiments with the above metric adaptation approach. If it performs well, I'll be wanting to add it to AdvancedHMC (open issue here: TuringLang/AdvancedHMC.jl#277) and will need HVPs in the API here, so I'll open a PR here to propose API changes.