Computing Hessian when the cost function consists of elementary functions

[Kumar Shaurya Shankar]

Hi all

I am trying to use the stan C++ autodiff library for my visual pose tracking algorithm. Essentially I have a (reasonably complex) non linear cost function that warps a reference image and computes a difference measure between the warped and reference image. The parameters to my warp are defined by a 6 element vector (3 for translation, and 3 for rotation)

I have been able to get the Jacobian of my cost function using stan, but I'm stuck at trying to obtain the hessian. The crux of the problem comes down to elementary functions like sin and cos trying to operate on the stan::math::fvar<stan::math::var> template parameters. If I modify my code to explicitly use the .val() method to obtain the stan::math::var element, although the code compiles, I think the gradients aren't propagated, and so I get incorrect (0) jacobians and hessians.

Could someone please point me how to go about solving this issue? I tried to look at the sample code in earlier posts and unit test code, but they all use Eigen matrix elements, and not the simpler methods such as sin, cos, log. Is it the case that the master branch hasn't implemented the nested templatized math functions?

Thanks a lot!

Regards

Shaurya

[Bob Carpenter]

See:

stan/math/mix/mat/functor/hessian.hpp

The individual fvar<T> functions have not been deeply tested
in most case, so anything involving fvar should be taken as
alpha or beta code. Having said that, you can use

stan/math/prim/mat/functor/finite_diff_hessian.hpp

to test.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Andrew Gelman]

Hi, I have no idea but I do seem to recall that 2nd and 3rd derivatives will be in Stan soon?

A

On May 6, 2016, at 6:44 PM, Kumar Shaurya Shankar <kumar....@gmail.com> wrote:

[Kumar Shaurya Shankar]

As a minimal example

  template<typename Derived> Derived operator()(

      const Eigen::Matrix<Derived, Eigen::Dynamic, 1>& theta) const {

    return sin(theta(0)*theta(1));

  }

This works with numeric differentiation
stan::math::finite_diff_hessian(f, theta, cost, G, H );

as you recommended, but doesn't compile with hessian

stan::math::hessian(f, theta, cost, G, H);

error: no matching function for call to ‘sin(stan::math::fvar<stan::math::var>)’

I just want to confirm that this is expected. 

Thanks!

[Daniel Lee]

Hi Kumar,

1) it helps to provide the code -- I was trying to piece together your example, but it isn't complete.

2) is the class f defined in the stan::math namespace? If not, you'll need to change it to:

template<typename Derived>
Derived operator()(const Eigen::Matrix<Derived, Eigen::Dynamic, 1>& theta) const {

  using stan::math::sin;

  return sin(theta(0)*theta(1));
}

Daniel

[Kumar Shaurya Shankar]

Apologies for not providing the full example code. Adding the sin namespace didn't help,

error: invalid initialization of reference of type ‘const stan::math::var&’ from expression of type ‘stan::math::fvar<stan::math::var>’

---

#include <cmath>

#include <stan/math.hpp>

#include <iostream>

#include <Eigen/Dense>

#include <stan/math/mix/mat/functor/hessian.hpp>

#include <stan/math/prim/mat/functor/finite_diff_hessian.hpp>

struct test_functor {

  template<typename Derived> Derived operator()(

      const Eigen::Matrix<Derived, Eigen::Dynamic, 1>& theta) const {

    using stan::math::sin;

    return sin(theta(0)*theta(1));

  }

};

int main() {

  Eigen::Matrix<double, Eigen::Dynamic, 1> theta = Eigen::Matrix<double,

      Eigen::Dynamic, 1>::Constant(6, 0), G = Eigen::Matrix<double,

      Eigen::Dynamic, 1>::Zero(6);

  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> H = Eigen::Matrix<

      double, Eigen::Dynamic, Eigen::Dynamic>::Zero(6, 6);

  test_functor f;

  double cost;

  stan::math::hessian(f, theta, cost, G, H); // <-- This fails

  std::cout << "Gradient computed is \n" << G << "\n and Hessian is \n" << H;

  return 0;

}

[Daniel Lee]

Instead of your current includes, use this:

#include <stan/math/mix/mat.hpp>
#include <cmath>
#include <iostream>

I was able to get it to run by doing that. In order to get your code to work with the finite_diff_hessian() functor, you'll also need to add "using std::sin;" to the operator() you've definited.

The stan/math.hpp include only includes reverse mode (that should be clear in the autodiff paper, but it's a long read and we don't describe how to use forward mode there).

We've now moved to having these includes:
- stan/math/prim/{scal, arr, mat}.hpp for primitives without using autodiff. scal is for scalar functions (no std::vector or Eigen), arr is for vector (std::vector, but no Eigen), and mat is for Eigen and std::vector.

- stan/math/rev/{scal, arr, mat}.hpp for reverse mode autodiff. stan/math.hpp includes stan/math/rev/mat.hpp
- stan/math/fwd/{scal, arr, mat}.hpp for forward mode autodiff.
- stan/math/mix/{scal, arr, mat}.hpp for mixed mode autodiff. This is what you need. Since you're using Eigen, you'll need to bring in the mat.hpp file.

We really need to add instructions on how to use forward mode and mixed forward / backward mode autodiff.

Daniel

[Kumar Shaurya Shankar]

a. Oh! So THAT's what all the subfolder hierarchy is for. It all makes sense now.

b. So that means I checkout the develop branch? stan/math/mix/mat.hpp is not in release (2.9.0)

[Daniel Lee]

Yes, you'll need the develop version for that to work. In v2.9.0, you'll need to include the headers in a particular order to get the template metaprograms instantiated properly. We included those headers for convenience. (We should have done it long before.)

[Kumar Shaurya Shankar]

OK, great, thanks a lot! I was able to compile with that. I'll try integrating it with the rest of my code now

I didn't have the cvodes library(Is it this - http://computation.llnl.gov/projects/sundials-suite-nonlinear-differential-algebraic-equation-solvers/cvodes?), so I just commented out the includes in math/prim/arr.hpp and math/rev/arr.hpp

[Kumar Shaurya Shankar]

So I finally got the autodiff hessian to work, and I expected some spike in memory usage, but it turns out that my code requires ~12GB of RAM to run, haha. Ouch.

[Daniel Lee]

Is that for one call to hessian? There is some management to do with regards to the autodiff stack. We haven't really pushed on the forward mode, so it's possible we haven't correctly reused the allocated memory or we weren't clear about what needs to happen between iterations to reset memory. 

If it's multiple calls, can you see if the memory footprint grows and doesn't shrink?

[Kumar Shaurya Shankar]

Yes, just the one call to the hessian. The memory growth ramps up the first time my cost function is called, and then remains more or less constant through the subsequent calls to the cost function.

Also, the hessian I obtain by autodiff vs the one obtained using finite differences is drastically different (i.e. the Hessian goes from negative definite (as expected) with the autodiff to a hessian with mixed eigen values). The gradients are pretty close, however.

[Bob Carpenter]

How many dimensions do you have? At the very least,
it's going to require N^2 * 20 bytes plus about 40
bytes per operator/function call in computing the function
and its gradient.

The reason to use autodiff is that finite differences
have low precision and you need to nest one inside the other
to calculate Hessians. We only use the finite diffs to
validate our higher-order autodiff code---we don't recommend
using them for real work.

The hessian() code does call recover_memory() after each row
is computed. We should probably modify it to define H(i,j) = H(j,i)
just to make sure everything's symmetric.

- Bob

[Kumar Shaurya Shankar]

Well, that makes sense. I'm operating on 640x480 images, so..yeah. Just the storage requirement for the image comes down to ~4GB based on what you mentioned.