Computing the gradient of another gradient

General
1

Dear Stan Community,

As part of a Markov chain Monte Carlo sampling algorithm, I need to compute the following expression
\nabla |\nabla \xi(x)|
for some function \xi that maps a high dimensional real space to a real number. Using the reverse mode ‘gradient’ function, I can compute
|\nabla \xi(x)|,
in Stan, but I’m unable to implement the additional gradient without using the `hessian’ function. This last function, however, seems currently not fully implemented in Stan.

Would someone have an idea on how to compute the expression above using only the ‘gradient’ functionality in Stan? The main technical problem is that I can’t seem to use a

Eigen::Matrix<var, Eigen::Dynamic, 1> grad_fx

to store the gradient when calling

gradient(f, x, fx, grad_fx).

Only the numerical type ‘double’ is allowed.

Any tips on how to implement this ?

Hannes

1 Like
2

I was going to forward you here: Using the Stan Math C++ Library

But it seems that uses the hessian function.

What error are you getting?

The usual way to compute hessians is to nest forward mode in reverse mode. So your types will look like fvar<var> now.

Here is the implementation of hessian if that helps: https://github.com/stan-dev/math/blob/develop/stan/math/mix/functor/hessian.hpp

@Bob_Carpenter has a draft autodiff manual here you might look at: https://github.com/bob-carpenter/ad-handbook/blob/master/ad-handbook-draft.pdf

1 Like
3

Thank you for your reply, the implementation of the ‘hessian’ function is really helpful.

Since the numerical type of the functor changes to fvar<var>, my functor looks like this

fvar<var> operator()(const Eigen::Matrix<fvar<var>, Eigen::Dynamic, 1>& x) const {
	fvar<var> dab = sqrt((x[0]-x[3])*(x[0]-x[3]) + (x[1]-x[4])*(x[1]-x[4]) + (x[2]-x[5])*(x[2]-x[5]));
	fvar<var> dbc = sqrt((x[3]-x[6])*(x[3]-x[6]) + (x[4]-x[7])*(x[4]-x[7]) + (x[5]-x[8])*(x[5]-x[8]));
	fvar<var> inner = ((x[0]-x[3])*(x[6]-x[3]) + (x[1]-x[4])*(x[7]-x[4]) + (x[2]-x[5])*(x[8]-x[5]))/(dab*dbc);
	fvar<var> theta = acos(inner);
	return theta;
}

This function basically computes the angle between a set of three atoms. This code, however, runs into a compilation error

/Users/hannesvdc/StanHeaders/inst/include/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:45:

38: error: no member named 'RowsAtCompileTime' in 'stan::math::fvar<stan::math::var>'typedef 

Eigen::Matrix<scalar_t, T::RowsAtCompileTime, T::ColsAtCompileTime>

Either I’m implementing the functor wrongly, or I don’t have the latest version of StanHeaders (my current version of from 2019-09-05).

Any ideas?

Hannes

4

It’s not obvious to me what is breaking down. Can you share the whole code and I can mess with it over here?

Can you take the eigen stuff out and get Hessians of just simple expressions?