Precomputing strategy for Hamiltonian Monte Carlo method based on regularity in parameter space

Abstract

Markov Chain Monte Carlo (MCMC) algorithms play an important role in statistical inference problems dealing with intractable probability distributions. Recently, many MCMC algorithms such as Hamiltonian Monte Carlo (HMC) and Riemannian Manifold HMC have been proposed to provide distant proposals with high acceptance rate. These algorithms, however, tend to be computationally intensive which could limit their usefulness, especially for big data problems due to repetitive evaluations of functions and statistical quantities that depend on the data. This issue occurs in many statistic computing problems. In this paper, we propose a novel strategy that exploits smoothness (regularity) in parameter space to improve computational efficiency of MCMC algorithms. When evaluation of functions or statistical quantities are needed at a point in parameter space, interpolation from precomputed values or previous computed values is used. More specifically, we focus on HMC algorithms that use geometric information for faster exploration of probability distributions. Our proposed method is based on precomputing the required geometric information on a set of grids before running sampling algorithm and approximating the geometric information for the current location of the sampler using the precomputed information at nearby grids at each iteration of HMC. Sparse grid interpolation method is used for high dimensional problems. Tests on computational examples are shown to illustrate the advantages of our method.

Appendices

Appendix 1: Convergence to the correct distribution

In order to prove that the equilibrium distribution remains the same, it suffices to show that the detailed balance condition still holds. Note that the alternative Hamiltonian $\tilde{H}(q,p)$ defines a surrogate-induced Hamiltonian flow, parameterized by the trajectory length t, which is a map $\tilde{\phi }_t:\;(q,p)\rightarrow (q^*,p^{*})$. Here, $(q^*,p^*)$ is the end-point of the trajectory governed by the following equations

$$\begin{aligned} \frac{dq}{dt} = \frac{\partial \tilde{H}}{\partial p} = M^{-1}p,\quad \frac{dp}{dt} = -\frac{\partial \tilde{H}}{\partial q} = -\frac{\partial \tilde{U}}{\partial q} = -\tilde{F} \end{aligned}$$

Denote $\theta =(q,p),\;\theta '=(q^*,p^*)=\tilde{\phi }_t(\theta )$. In the Metropolis-Hasting step, we use the original Hamiltonian to compute the acceptance probability

$$\begin{aligned} \alpha (\theta ,\theta ') = \min (1,\exp [-H(\theta ')+H(\theta )]) \end{aligned}$$

therefore,

$$\begin{aligned}&\alpha (\theta ,\theta '){\mathbb {P}}(d\theta ) = \alpha (\theta ,\theta ')\exp [-H(\theta )]d\theta \\&\quad \mathop {=}\limits ^{\theta =\tilde{\phi }_t^{-1}(\theta ')} \min (\exp [-H(\theta )],\exp [-H(\theta ')])\left| \frac{d\theta }{d\theta '}\right| d\theta '\\&\quad = \alpha (\theta ',\theta )\exp [-H(\theta ')]d\theta '\\&\quad = \alpha (\theta ',\theta ) {\mathbb {P}}(d\theta ') \end{aligned}$$

since $\left| \frac{d\theta }{d\theta '}\right| = 1$ due to the volume conservation property of the surrogate induced Hamiltonian flow $\tilde{\phi }_t$. Now that we showed the detailed balance condition is satisfied, along with the reversibility of the surrogate induced Hamiltonian flow, the modified Markov chain will converge to the correct target distribution.

Appendix 2: More on sparse grid

1.1 Construction of sparse grid and basis functions

The Clenshaw–Curtis type sparse grid $H^{CC}$ introduced in Sect. 3.4.2 is constructed from the following formulas. Here, the $x_j^i$ are defined as

$$\begin{aligned} x_j^i =\left\{ \begin{array}{ll}(j-1)/(m_i-1),\; &{}j = 1,\ldots ,m_i, m_i > 1\\ 0.5, &{} j=1, m_i =1 \end{array}\right. \end{aligned}$$

In order to obtain nested sets of points, the number of nodes is given by

$$\begin{aligned} m_1 = 1\quad \text { and }\quad m_i = 2^{i-1} + 1 \text { for } i > 1 \end{aligned}$$

Piecewise linear basis functions a can be used for the univariate interpolation formulas $U^i(f)$.

$$\begin{aligned} a_1^1(x) = 1,\quad a_j^i(x) = \left\{ \begin{array}{ll} 1-(m_i-1)\cdot |x-x_j^i|,\;&{} |x-x_j^i| < 1/(m_i-1),\\ 0, &{} \text {otherwise} \end{array}\right. \end{aligned}$$

for $i>1$ and $j =1,\ldots ,m_i$.

Fig. 15

Piecewise linear Nodal basis (a) and hierarchical functions (b) with support nodes $x_j^i\in X_{\Delta }^i,i=1,2,3$ for the Clenshaw–Curtis grid

Fig. 16

Piecewise linear interpolation: nodal versus Hierarchical. a Nodal basis. b Hierarchical basis

1.2 Derivation of the hierarchical form

With the selection of nested sets of points, we can easily transform the univariate nodal basis into the hierarchical one. By definition, we have

$$\begin{aligned} \Delta ^i(f)&= U^i(f) - U^{i-1}(f) \\&= \sum _{j=1}^{m_i}f(x_j^i)\cdot a_j^i - \sum _{j=1}^{m_i} U^{i-1}(f)(x_j^i)\cdot a_j^i\\&= \sum _{j=1}^{m_i}\big (f(x_j^i) - U^{i-1}(f)(x_j^i)\big ) \cdot a_j^i \end{aligned}$$

since $f(x_j^i) - U^{i-1}(f)(x_j^i) = 0,\; \forall \; x_j^i \in X^{i-1}$,

$$\begin{aligned} \Delta ^i(f) = \sum _{x_j^i\in X^i_{\Delta }}\big (f(x_j^i)-U^{i-1}(f)(x_j^i)\big )\cdot a_j^i \end{aligned}$$

(6.1)

From (6.1) we note that for all $\Delta ^i(f)$, only the basis functions belonging to the grid points that have not yet occurred in a previous set $X^{i-k},\;1\le k \le i-1$ are involved. Figure 15 gives a comparison of the nodal and the hierarchical basis functions and Fig. 16 shows the construction of the interpolation formula using nodal basis functions and function values versus using hierarchical basis functions and hierarchical surpluses for a univariate function f. Both figures are reproductions based on Klimke and Wohlmuth (2005).

Applying the tensor product formula (3.7) with $\Delta ^i$ given in (6.1), the hierarchical update in the Smolyak algorithm (3.8) now can be rewritten as

$$\begin{aligned} \Delta A_{q,d}(f) =&\sum _{|{\varvec{i}}|=q}(\Delta ^{i_1}\otimes \cdots \otimes \Delta ^{i_d})(f)\\ =&\sum _{|{\varvec{i}}|=q}\sum _{x_{j_1}^{i_1}\in X_{\Delta }^{i_1}}\cdots \sum _{x_{j_d}^{i_d}\in X_{\Delta }^{i_d}}\big (f(x_{j_1}^{i_1},\ldots ,x_{j_d}^{i_d})-A_{q-1,d}(f)(x_{j_1}^{i_1},\ldots ,x_{j_d})\big )\nonumber \\&\cdot (a_{j_1}^{i_1}\otimes \cdots \otimes a_{j_d}^{i_d}). \end{aligned}$$