Access Hessian from optimize method in CmdStanR

Is it possible to access the Hessian from the output when using cmdstan’s optimize method in R or to specify a draws argument so I can get estimates of standard errors on the constrained parameter space?

I was using the optimize function from rstan, but when I moved over to using my universities’ HPC I ran into issues with rstan. I’ve found a workaround with CmdstanR but I can’t seem to figure out where in the output object the Hessian might be stored, or how to return it, or how to get draws from the constrained parameter space.

Here’s how I would get what I want from rstan::optimizing():

    library(cmdstanr)
    library(rstan)
    # Use example model from cmdstan
    file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan")
    mod <- stan_model(file)
    stan.data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
    # Get optimization results
    stan_output <- optimizing(mod,
                              data = stan.data,
                              algorithm = "LBFGS",
                              hessian = TRUE,
                              draws = 50000)
    # Extract parameter names
    par_names <- names(stan_output$par)
    # Get parameter estimates
    par_values <- stan_output$par
    # By grabbing many draws of the parameters, we can take the square root
    # of the variances in the cov matrix for these draws to estimate SE
    par_se <- sqrt(diag(cov(stan_output$theta_tilde)))
    # Get our t-values
    t_val <- par_values / par_se
    output <- data.frame(parameters = par_names,
                         est = par_values,
                         se = par_se,
                         t_val = t_val,
                         row.names = NULL)

Tagging @mitzimorris, because I honestly don’t know…

tldr: not yet; we’re working on it.

this came up at a recent Gelman group Stan meeting - my recollection is that @bgoodri said that there were problems with what RStan is doing. we agreed that this should be done in core Stan. as soon as it happens, it can then be exposed via CmdStan and the CmdStanX wrappers.

5 Likes

Hello @mitzimorris, @bgoodri,

Is it possible to extract the Hessian from a cmdstanMLE object now?
Thanks in anticipation of the response.

thanks for checking, sorry to report, no progress.

2 Likes

Hi @mitzimorris ,
I was wondering whether there might have been any further process. Is it now possible to extract the Hessian from a cmdstanMLE object? From the manual, there is an option as follows so I had a hope that I might be able to access hessian when using L-BFGS. Is it correct?
history_size (positive integer) The size of the history used when approximating the Hessian. Only available for L-BFGS.
Thanks.
P.S. Hessian please (if not supported yet)? 😊

1 Like

+1 to that. I’ve been using RStan just because it gives me access to the Hessian and cmdstanr does not.

Just note that it’s not necessarily giving you the Hessian you would like to get, see e.g. Visual illustration of Jacobian of parameter transformation

3 Likes

This is very helpful. Thanks, Aki.

Any updates on this?

It’s now available in the github version of cmdstanr, with background on usage in the PR description: Add optional methods for log_prob, grad_log_prob, hessian, un/constrain pars by andrjohns · Pull Request #701 · stan-dev/cmdstanr · GitHub

5 Likes

Rad. Thank you!

@andrjohns, @mitzimorris Would the exposed functions in cmdstanR be straightforwardly supported in cmdstanPy? Sorry I don’t know CmdStanX wrappers as much as you do. But this feature is so cool!

The log_prob method which is already implemented directly in cmdstan could be exposed fairly easily, but for the other methods (e.g., hessian() or exposing user-defined functions), these would need to be handled using python’s c++ bindings which is a little trickier. I believe @WardBrian has been experimenting with it though, so he might have more info

1 Like

The log_prob method as implemented in CmdStan doesn’t return the Hessian, just the gradients.

If you want the Hessian, you should use BridgeStan.

2 Likes

What’s the current status? I’m using the latest cmdstanr and the latest production release of cmdstan.

CmdStanR has support for getting Hessians: Calculate the log-probability , the gradient w.r.t. each input, and the hessian for a given vector of unconstrained parameters — fit-method-hessian • cmdstanr

CmdStan also supports outputting the Hessian used as part of Laplace approximation if you enable the diagnostic file output, it will be a JSON file with an entry called Hessian

4 Likes