Hessian corrections to Hybrid Monte Carlo

TL;DR

This paper introduces a novel approach to Hybrid Monte Carlo that incorporates second-order derivatives of the log likelihood, improving efficiency by computing the Hessian only at trajectory endpoints.

Contribution

The method allows Hessian information to be used in HMC without the need for frequent Hessian evaluations during trajectories, enhancing computational efficiency.

Findings

Hessian corrections improve sampling accuracy in HMC.

The method reduces computational cost compared to traditional Hessian-based approaches.

Hessian evaluations are limited to start and end points of trajectories.

Abstract

A method for the introduction of second-order derivatives of the log likelihood into HMC algorithms is introduced, which does not require the Hessian to be evaluated at each leapfrog step but only at the start and end of trajectories.

Equations

$$\mathcal{H}=\frac{1}{2}{\mathbf{p}}^2-l\text{ ,}$$

$$\frac{\mathrm{d}\mathbf{\theta }}{\mathrm{d}t}=\frac{\partial \mathcal{H}}{\partial \mathbf{p}}=\mathbf{p}\text{ ,}\frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t}=-\frac{\partial \mathcal{H}}{\partial \mathbf{\theta }}=\partial l\text{ .}$$

$${\mathbf{\theta }}_0$$

$${\mathbf{p}}_{1/2}$$

$${\mathbf{\theta }}_i$$

$${\mathbf{p}}_{i+(1/2)}$$

$${\mathbf{\theta }}^{\ast }$$

$${\mathbf{p}}^{\ast }$$

$$l({\mathbf{\theta }}^n+\mathbf{x})\approx l({\mathbf{\theta }}^n)+{\mathbf{v}}^{\top }\mathbf{x}+\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}\text{ ,}$$

$$v_i:={\frac{\partial l}{\partial {\theta }^i}|}_{{\mathbf{\theta }}^n}\text{ ,}{H}_{ij}:={\frac{{\partial }^{2}l}{\partial {\theta }^i\partial {\theta }^j}|}_{{\mathbf{\theta }}^n}\text{ ,}\mathbf{v}:=(v_i)\text{ ,}\mathbf{H}:=(H_{ij})\text{ .}$$

$$\mathrm{d}\mathbf{x}=\left(\mathbf{Ax}+\mathbf{b}\right)\mathrm{d}t\text{ ,}$$

$$\mathbf{x}=\left(\begin{array}{c}\mathbf{\theta }\\ \mathbf{p}\end{array}\right)\text{ ,}\mathbf{A}=\left(\begin{array}{cc}0& \text{mathbbm}1\\ \mathbf{H}& 0\end{array}\right)\text{ ,}\mathbf{b}=\left(\begin{array}{c}0\\ \mathbf{v}\end{array}\right)\text{ .}$$

$$\frac{\mathrm{d}\mathbf{m}}{\mathrm{d}t}=\mathbf{Am}+\mathbf{b}\text{ ,}\frac{\mathrm{d}\mathbf{S}}{\mathrm{d}t}=\mathbf{A}\mathbf{S}+\mathbf{S}{\mathbf{A}}^{\top }\text{ ,}\mathbf{m}(0)=\left(\begin{array}{c}{\mathbf{\theta }}^n\\ \mathbf{q}\end{array}\right)\text{ ,}\mathbf{S}(0)=\left(\begin{array}{cc}0& 0\\ 0& \mathbf{Q}\end{array}\right)\text{ .}$$

$$\mathrm{Mean}({\mathbf{\theta }}^{\ast }(\delta ))$$

$$\mathrm{Cov}({\mathbf{\theta }}^{\ast }(\delta ))$$

$$\mathbf{q}({\mathbf{\theta }}^n)=\cosh \left({\mathbf{H}}^{1/2}\delta \right){\mathbf{H}}^{-1/2}{\left(\sinh \left({\mathbf{H}}^{1/2}\delta \right)\right)}^{-1}\partial l\text{ ;}\mathbf{Q}({\mathbf{\theta }}^n)=-{\left(\sinh \left({\mathbf{H}}^{1/2}\delta \right)\right)}^{-2}\text{ .}$$

$$\frac{\rho ({\mathbf{\theta }}^n,{\mathbf{p}}^n|{\mathbf{\theta }}^{\ast },{\mathbf{p}}^{\ast })}{\rho ({\mathbf{\theta }}^{\ast },{\mathbf{p}}^{\ast }|{\mathbf{\theta }}^n,{\mathbf{p}}^n)}=\frac{\mathcal{N}(-{\mathbf{p}}^{\ast }|\mathbf{q}({\mathbf{\theta }}^{\ast }),\mathbf{Q}({\mathbf{\theta }}^{\ast }))}{\mathcal{N}({\mathbf{p}}^n|\mathbf{q}({\mathbf{\theta }}^n),\mathbf{Q}({\mathbf{\theta }}^n))}$$

$$\pi (\mathbf{\theta })=\mathcal{N}(\mathbf{\theta }|0;\mathbf{V})\text{,}$$

Full text

1 Introduction

Markov chain Monte Carlo (MCMC) is a highly influential computationally intensive method for performing Bayesian inference, with a large variety of applications (Brooks et al., 2011). While earlier MCMC algorithms made use of random walks in parameter space (Gilks et al., 1995), as highlighted in a recent review by Green et al. (2015), the use of derivatives can lead to improved algorithms.

One derivative-based approach is HMC, standing for either Hybrid (Duane et al., 1987) or Hamiltonian (Neal, 2011) Monte Carlo, which requires first derivatives of the log likelihood to be available. More recently, second derivatives have been included via the use of geometric approaches (Girolami and Calderhead, 2011, Betancourt, 2013). This working paper introduces a different route to inclusion of second-order derivatives through truncated Taylor expansion of the log-likelihood, after which Hamilton’s equations can be solved exactly without further approximation. This algorithm is called HHMC (for Hessian-corrected HMC) and is able to sample accurately from distributions with different scales for each parameter, which is challenging for standard HMC.

2 A Hessian HMC algorithm

2.1 Local solution of Hamilton’s equations

The idea behind HMC is to propose a new set of parameters $\boldsymbol{\theta}^{*}$ starting from $\boldsymbol{\theta}^{n}$ by making use of additional ‘momentum’ variables $\mathbf{p}$. We will consider the general case where $\mathbf{p}^{n}\sim\mathcal{N}(\mathbf{q},\mathbf{Q})$ in general, although in standard HMC, $\mathbf{q}=\mathbf{0}$ and $\mathbf{Q}=\mathbbm{1}$.

Supposing our aim is to sample according to $\pi$, and we let $l=\mathrm{ln}(\pi)$, then the basis for HMC algorithm design is Hamilton’s equations. These make use of the Hamiltonian

[TABLE]

and take the form

[TABLE]

In general, these cannot be solved analytically and so proposals are made on the basis of a numerical approximation, for example the leapfrog method, in which $L$ steps of length $\varepsilon$ are

[TABLE]

This algorithm, together with the initial random choice of $\mathbf{p}^{n}$, defines a marginal proposal density $\rho(\boldsymbol{\theta}^{*}|\boldsymbol{\theta}^{n})$. This will be close to the solution of (2) over a time period $\delta:=\varepsilon L$ for large $L$ and small $\varepsilon$.

Now suppose that we approximate $l=\mathrm{ln}(\pi)$ in the neighbourhood of some value $\boldsymbol{\theta}^{n}$ through Taylor expansion

[TABLE]

where

[TABLE]

Then we can approximate Hamilton’s equations in the region of $\boldsymbol{\theta}^{n}$ through the linear SDE

[TABLE]

where

[TABLE]

Note that this SDE does not have a white-noise term, but is nevertheless a special case of the results of Archambeau et al. (2007), and therefore has Gaussian solution with mean $\mathbf{m}$ and covariance matrix $\mathbf{S}$ obeying

[TABLE]

Solving these ODEs over the interval $[0,\delta]$ gives, after some analytic work,

[TABLE]

If we then impose that this should correspond to a situaion where the marginal proposal $\rho(\boldsymbol{\theta}^{*}|\boldsymbol{\theta}^{n})=\pi(\boldsymbol{\theta}^{*})$ under the approximation (4), then we obtain a solution for the proposal distribution after additional analytic work:

[TABLE]

In terms of numerical computation of (10), despite the seeming complexity of the matrix functions, the general approach of Davies and Higham (2003) is applicable. In particular, despite the appearance of fractional powers in the compact expression (10), all of the functions involved have only integer matrix powers in their Taylor series and so is expressible as matrix polynomials.

Note that we can no longer write the accept-reject probabilities in this algorithm in terms of the Hamiltonian, however since the leapfrog integrator (3) is reversible it is still possible to calculate these probabilities. In particular, if we start with $\boldsymbol{\theta}^{n},\mathbf{p}^{n}$ and propose $\boldsymbol{\theta}^{*},\mathbf{p}^{*}$,

[TABLE]

The MCMC algorighm based on (10) and (11) is called HHMC for Hessian-corrected HMC. This algorithm does not require the Hessian $\mathbf{H}$ at each leapfrog step, just at the start and end of trajectories, and does not involve third-order derivatives $\partial_{i}\mathbf{H}$ making it much less computationally costly than Riemannian approaches (Girolami and Calderhead, 2011, Betancourt, 2013).

3 A target distribution with heterogeneous scales

Following Neal (2012), consider the following multivariate Gaussian target density:

[TABLE]

where $\mathbf{V}$ is a diagonal matrix with entries equal to the squares of: 110; 100; twenty-six evenly spaced standard deviations between 16 and 8; 1.1; and 1.0. Neal suggests this distribution as a diagnostic for HMC because the variable scales associated with each parameter create difficulties for the algorithm.

Running MCMC chains of length $10^{3}$ using $L=10$ leapfrog steps each with a step size of $\epsilon=0.2$ for HMC and HHMC gives the results shown in Figures 1, 2, 3 and 4. The results show that HMC is out-performed by HHMC.

Acknowledgements

Work supported by the UK Engineering and Physical Sciences Research Council.