Expose model hessian in Rstan

Hi Stan-dev,

I’d like to expose model Hessians to R through Rstan with something like grad_log_prob. I have been told by Mike B that the Stan team doesn’t want to do this in a full release because the autodiff Hessians aren’t ready for all models, so I’m happy to hack it up myself in a custom branch. But I think I need some additional guidance.

Looking at the rstan source, it appears that stan_fit.hpp is the first step in exposing C++ functions to R. So I would like to write my own version of

SEXP grad_log_prob(SEXP upar, SEXP jacobian_adjust_transform)

However, in the Stan model C++ code, stan::model::log_prob_grad appears to have a different structure than stan::model::hessian

In particular, I’m not quite sure how to pass the arguments (par_r and par_i in the log_prob_grad function) into the Hessian. The unit test appears to use a stan::io::dump object:

Upon seeing this, I thought maybe I had better just ask for help. :P Can you someone give me a hand figuring this out?

I don’t know the details of the R interface, but from the Stan side, what you need to do is get the functor corrsponding to the model log density then pass that into the Hessian functional. That’s what the algorithm code does with the model to compute gradients.

Ok, as of now I’m actually stuck on the C++ part, so maybe you could help me see how to do that. Specifically, what should I put in the “hessian here” comment of the below function? (I mostly copied it from the current Rstan grad_log_prob) upar is a vector of the real-valued model parameters (the things we’re taking the Hessian with respect to).

    SEXP hessian_log_prob(SEXP upar, SEXP jacobian_adjust_transform) {
      std::vector<double> par_r = Rcpp::as<std::vector<double> >(upar);
      if (par_r.size() != model_.num_params_r()) {
        std::stringstream msg;
        msg << "Number of unconstrained parameters does not match "
               "that of the model ("
            << par_r.size() << " vs "
            << model_.num_params_r()
            << ").";
        throw std::domain_error(msg.str());
      std::vector<int> par_i(model_.num_params_i(), 0);
      std::vector<double> gradient;
      SOMETHING hessian;
      double lp;
      if (Rcpp::as<bool>(jacobian_adjust_transform))
        // Previously:
        // lp = stan::model::log_prob_grad<true,true>(model_, par_r, par_i, gradient, &rstan::io::rcout);
        // Hessian here!
        // Previously:
        // lp = stan::model::log_prob_grad<true,false>(model_, par_r, par_i, gradient, &rstan::io::rcout);
        // Hessian here!

      // ... then Rcpp wrapping stuff...

Look at stan/math/mix/mat/functor/hessian.hpp and you’ll see that it’s an Eigen matrix that gets computed. The function’s well documented.

Cool, and what goes in here?

It’s C++, you just have to look at the signature for the function you are calling. If you ask about producing specific pieces for the call you might get more helpful answers.

Let me try to be more specific.

I can see that the hessian function takes as arguments a model and memory that has been pre-allocated to store the value of the log prob, the gradient, and the Hessian.

hessian(const M& model,
                 const Eigen::Matrix<double, Eigen::Dynamic, 1>& x,
                 double& f,
                 Eigen::Matrix<double, Eigen::Dynamic, 1>& grad_f,
                 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& hess_f,
                 std::ostream* msgs = 0)

I assume that the model object that is passed in here has had its parameter values set at some point, since of course the Hessian is evaluated at some particular set of parameter values and they do not appear here. My question is how to set them. That is, suppose I give you a model object, a params_r vector, and a params_i vector. Howe to I get a version of the model that I can pass into the hessian function so I get the hessian at those parameter values?

Note that, unlike the hessian function, log_prob_grad takes par_r and par_i as arguments. I have looked at how they are used inside log_prob_grad. They are passed as arguments to the construction of a var type which is then differentiated.

var adLogProb
          = model.template log_prob<propto, jacobian_adjust_transform>
          (ad_params_r, params_i, msgs);

But I don’t see how to use this to produce something I could pass to the hessian function.

Does this make sense?

I don’t recall syntax off hand but if you compile a really simple model with CmdStan and look at the .hpp file you’ll see the methods that are available for setting parameters.

@rgiordan If you look at the code in stan/model/ you’ll see that the hessian functional defined there does not take in a model object directly but rather something called a model_functional (which is actually a functor – I think that one was my fault from a long time ago). Anyways, in stan/model/model_functional.hpp you’ll see that model_functional just wraps a model instance and hides the somewhat ungainly log_prob function in an operator() so that we can call all of the nice higher-order functionals. Following the files we then see that the input vector x is the vector of unconstrained parameter values at which you want to evaluate the Hessian.

Perfect, I understand now. It’s simpler than I thought it was! I guess params_i is not actually used for anything. And it looks like if I want to make jacobian_adjust_transform an option I could just make the functor myself and call stan::math::hessian on it rather than go through the stan::model::hessain function.

Thanks! I’ll follow up if I have any more problems.

Correct – params_i is not used.