Add rev functor to compute hessian-times-vector by finite differences

[WardBrian]

Summary

This adds a function analogous to finite_diff_hessian_auto for computing Hessian-vector products of functions which only support rev mode types. This can be done with only two calls to gradient, compared to the 2N calls needed for the naive approach of using finite_diff_hessian_auto and then multiplying.

The approach is described by @justindomke here: https://justindomke.wordpress.com/2009/01/17/hessian-vector-products/

In short,

If $g(x) = \frac{\partial f}{\partial x} f(x)$, then $g(x+\Delta x) \approx g(x) + H(x)\Delta x$
So
$H(x)v \approx \frac{g(x + rv) - g(x - rv)}{2r}$

See also: discussion in roualdes/bridgestan#126

Tests

I have duplicated the existing test used for autodiff Hessian-vector products. There is only one, so more might be nice.

Side Effects

None

Release notes

New function in stan::math::internal, finite_diff_hessian_vector_product_auto

Checklist

• Math issue: closes Hessian-vector product with central finite differences #2913

• Copyright holder: Simons Foundation

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

[aseyboldt]

I don't think it is valid to choose a different epsilon for each dimension, wouldn't that just change the vector we multiply the hessian with?

[WardBrian]

I don't think it is valid to choose a different epsilon for each dimension, wouldn't that just change the vector we multiply the hessian with?

I think you're right - this just so happened to have worked because the finite_diff_stepsize function uses
stepsize(u) = cbrt(epsilon) * max(1, abs(u)), and because v is normalized max(1,u) is one always, so the cube root of epsilon is used.

Do you think there is a better way of automatically deciding an epsilon, or should I just skip the call to finite_diff_stepsize and use cbrt(epsilon)?

[aseyboldt]

Do you think there is a better way of automatically deciding an epsilon, or should I just skip the call to finite_diff_stepsize and use cbrt(epsilon)?

I don't know of any better way, but that doesn't mean there isn't one... :-)

[andrjohns]

LGTM!

[WardBrian]

@andrjohns mind re-approving? The auto formatting commit dismissed your last review