Class mcmc_stepper_hmc (o2scl)

O2scl : Class List

template<class func_t, class data_t, class vec_t, class grad_t = std::function<int(size_t, const vec_t&, func_t&, vec_t&, data_t&)>, class vec_bool_t = std::vector<bool>>
class mcmc_stepper_hmc : public o2scl::mcmc_stepper_base<func_t, data_t, vec_t>

Hamiltonian Monte Carlo for MCMC.

The vectors step_fac, mom_step, and auto_grad are provided to give different values for each of the parameters. However, if these vectors have a smaller size, then the vector index is wrapped back around to the beginning (using the modulus operator).

The Hamiltonian is

H(q,p)=12pTM1plogL(q)
where M is the mass matrix and L is the likelihood. (The mcmc_para_base class presumes flat prior distributions.) This class assumes the mass matrix is the identity matrix.

The step is then accepted with probability

min{1,exp[Hnew]exp[Hold]}
Because the user-specified function, f, computes logL(q), the step should be accepted with probability
min{1,exp[iNpi,new2μi/2+fnew+iNpi,old2μi/2fold]}
where i is an index over N parameters and μi is the inverse mass for parameter i.

This class can compute the gradients automatically by finite-differencing or can use a gradient function specified by the user. The vector auto_grad controls this behavior. If the value of

auto_grad[ i % auto_grad.size() ]
is false, then it is presumed that the derivative with respect to the parameter with index i should be computed automatically by this stepper class. If it is true, then it is presumed that then the value will be obtained from the user-specified gradient function.

The gradient function, if specified, should be of the form

int grad(size_t nv, const vec_t &x, func_t &f,
vec_t &g, data_t &dat);

If the initial gradient calculation fails, then the HMC cannot proceed and the random walk (RW) algorithm from mcmc_stepper_rw is used as a fallback. If this fallback method is required frequently over the course of a simulation, then this may mean the combined HMC plus RW method may not be sampling the target distribution. This class doesn’t yet have a method for tracking this.

The algorithm is taken from [Neal11].

See the class documentation for mcmc_stepper_base for more information.

Public Functions

inline virtual const char *step_type()

Stepper type, “HMC”.

inline mcmc_stepper_hmc()
inline virtual ~mcmc_stepper_hmc()
inline void allocate(size_t n_params, size_t n_threads)

Allocate cache.

inline virtual void write_params(o2scl_hdf::hdf_file &hf)

Write stepper parameters to the HDF5 file.

inline void set_gradients(std::vector<grad_t> &vg)

Set the vector of user-specified gradients.

Note that this function stores a pointer to the vector of gradient objects and thus the user must ensure that this object is in scope when the MCMC is performed.

inline int grad_pot(size_t n_params, const vec_t &x, func_t &f, vec_t &g, data_t &dat)

Automatically compute the gradient using finite-differencing.

Note

The potential energy is the negative of the log-likelihood. This function returns the gradient of the log-likelihood, which is the negative of the gradient of the potential energy.

inline virtual void step(size_t i_thread, size_t n_params, func_t &f, const vec_t &current, vec_t &next, double w_current, double &w_next, const vec_t &low, const vec_t &high, int &func_ret, bool &accept, data_t &dat, rng<> &r, int verbose)

Construct a step.

This function constructs next and w_next, the next point and log weight in parameter space. The objective function f is then evaluated at the new point, the return value is placed in func_ret, and the step acceptance or rejection is stored in accept.

The first half step is:

p1/2=p0(ϵ)/2Uq(q0)
Then for i[1,N],
qi=qi1+ϵpi1/2if (i<N) pi+1/2=pi1/2ϵUq(qi)
The last half step is:
pN=pN1/2(ϵ)/2Uq(qN)

Public Members

vec_t step_fac

The factor controlling the step size for the fallback random walk (default is a 1-element vector containing 10.0)

size_t traj_length

Trajectory length (default 20)

prob_dens_gaussian pdg

Standard Gaussian for kinetic energy.

vec_t mom_step

Initial distribution in momentum space (default is a one-element vector containing 1.0)

double epsilon

Scale for gradient (default 0.2)

vec_bool_t auto_grad

Indicate which elements of the gradient need to be computed automatically (default is a one-element vector containing true).

For parameters in which this vector has an entry of true, it is assumed that the user-specified gradient object (if present) cannot compute the gradient and thus a simple numerical finite-differencing gradient is required.

double epsrel

The relative stepsize for finite-differencing (default 106 )

double epsmin

The minimum stepsize (default 1015)

ubmatrix param_cache

Cache for parameters.

ubmatrix grad_cache

Cache for gradients.

Public Static Attributes

static const size_t grad_failed = 30

Error if gradient failed.

Protected Attributes

std::vector<grad_t> *grad_ptr

Pointer to user-specified gradients.