On the variance reduction of Hamiltonian Monte Carlo via an approximation scheme

Abstract

Let $\pi$ be a probability distribution in $\mathbb{R}^d$ and f a test function, and consider the problem of variance reduction in estimating $\mathbb{E}_\pi(f)$. We first construct a sequence of estimators for $\mathbb{E}_\pi (f)$, say $({1}/{k})\sum_{i=0}^{k-1} g_n(X_i)$, where the $X_i$ are samples from $\pi$ generated by the Metropolized Hamiltonian Monte Carlo algorithm and $g_n$ is the approximate solution of the Poisson equation through the weak approximate scheme recently invented by Mijatović and Vogrinc (2018). Then we prove under some regularity assumptions that the estimation error variance $\sigma_\pi^2(g_n)$ can be as arbitrarily small as the approximation order parameter $n\rightarrow\infty$. To illustrate, we confirm that the assumptions are satisfied by two typical concrete models, a Bayesian linear inverse problem and a two-component mixture of Gaussian distributions.

1. Introduction and main results

Let $\pi$ be a target distribution in $\mathbb{R}^d$ and $f\colon \mathbb{R}^d\mapsto \mathbb{R}$ a test function. We are mainly concerned with the issue of variance reduction in estimating $\mathbb{E}_\pi(f)$ based on the samples generated by the Hamiltonian Monte Carlo algorithm. The variance reduction is a widely used method in statistical inferences. In fact, if $\pi$ is comparatively simple so that independent samples are easily drawn, say $(X_i, i\ge 1)$ , then we can use $({1}/{k})\sum_{i=1}^k f(X_i)$ as an estimator of $\mathbb{E}_\pi(f)$ by the well-known law of large numbers. The estimation error is asymptotically given by $({1}/{k})\mathbb{E}_\pi(f-\mathbb{E}_\pi(f))^2$ . However, if $\mathbb{E}_\pi(f-\mathbb{E}_\pi(f))^2$ is itself not negligible, then sufficiently many samples must be drawn in order to make the error as small as possible. Thus, in order to upgrade the efficiency of estimation we need to resort to the variance reduction method. That is, to find a function g with known $\mathbb{E}_\pi (g)$ , say $\mathbb{E}_\pi(g)=0$ and $\mathbb{E}_\pi(f+g-\mathbb{E}_\pi(f))^2 < \mathbb{E}_\pi(f-\mathbb{E}_\pi(f))^2$ , so that we might prefer to use $({1}/{k})\sum_{i=1}^k(f(X_i)+g(X_i))$ as an alternative estimator of $\mathbb{E}_\pi(f)$ by the law of large numbers again. A natural problem now arises: how do we construct such a g(x)? This is not an easy task in general.

On the other hand, we are also required to solve the issue of simulating the target distribution $\pi$ . In fact, this latter issue is more challenging and often appears in diverse research fields like high-dimensional data analysis, Bayesian statistical inference, high-performance numeric computation, and machine learning.

1.1. Metropolis–Hastings Monte Carlo

Markov chain theory turns out to offer a powerful tool for managing an issue like simulation, which resulted in the Markov chain Monte Carlo (MCMC) method. Assume that $(X_i,\,i\ge 0)$ is an ergodic Markov chain with $\pi$ as its stationary distribution. Then it follows by Birkhoff’s theorem that $\lim_{k\rightarrow\infty}\big[({1}/{k})\sum_{i=0}^{k-1}f(X_i)\big] = \mathbb{E}_\pi(f)$ $\pi$ -almost everywhere (a.e.). Therefore, the time mean $({1}/{k})\sum_{i=0}^{k-1}f(X_i)$ may be used as a good estimator for $\mathbb{E}_\pi(f)$ when the Markov chain is run for a sufficiently long time.

Constructing such a Markov chain that targets the desired distribution, however, is itself a nontrivial problem. Fortunately, various procedures have been outlined in the literature for automatically constructing appropriate transitions for any given target distribution, with the foremost among these the Metropolis–Hastings algorithm.

Let $\pi$ be a fixed probability measure, and assume it possesses a density function, denoted as $\pi(x)$ (with a minor abuse of notation). Start with a proposal transition kernel $Q(x,\cdot)$ with density function q(x, y) and run the Metropolis–Hastings algorithm to generate a chain $(X_{i},\,i\ge 0)$ as follows. Given the current state $X_i=x$ , sample a candidate $Y=y$ from the proposal kernel $Q(x,\cdot)$ , and then the calculate the acceptance rate:

(1.1) \begin{equation} \alpha(x, y)\;:\!=\;\min\bigg\{1,\dfrac{\pi(y)q(y,x)}{\pi(x)q(x, y)}\bigg\}.\end{equation}

Then independently draw a uniform random variable $Z\sim U[0,1]$ . If $Z\le \alpha(x,y)$ then set $X_{i+1}=y$ ; otherwise, set $X_{i+1}=x$ . The associated Markov kernel can be written as follows. For any set $A\in\mathcal{B}(\mathbb{R}^d)$ ,

\begin{equation*} P(x,A) = \int_{A}\alpha(x,y)q(x,y)\,\textrm{d}y + \delta_{x}(A)\int_{\mathbb{R}^d}(1-\alpha(x,y))q(x,y)\,\textrm{d}y.\end{equation*}

The Markov kernel $P(x,\textrm{d}y)$ is remarkably reversible with respect to $\pi$ , i.e. $\pi(x)P(x,\textrm{d}y) = \pi(y)P(y,\textrm{d}x)$ , and so the stationary distribution of the chain $(X_i,\,i\ge 0)$ is exactly the target distribution $\pi$ [5, Proposition 2.3.1]. It is worth noting that only the ratio ${\pi(y)}/{\pi(x)}$ is used in the acceptance/rejection probability (1.1), and thus it easily extends to a larger class of distributions, particularly like the Bayesian posterior distributions whose normalization constants are hardly computable.

As the reader may realize, the Metropolis–Hastings framework is quite flexible by instantiating it with different choices for the proposal transition kernel $Q(x, \textrm{d}y)$ , which in turn have a significant influence upon sampling. A simple and widely used choice is the so-called Metropolized random walk (MRW). This corresponds to simply taking a random walk with transition kernel $q(x, \cdot)\sim \mathcal{N}(x, hI_d)$ around the state space, where some steps are occasionally rejected. Note that the proposal is independent of the target $\pi$ , and so uses only a zeroth-order oracle for $\pi$ .

If $\pi$ has a smooth enough density function, particularly $\pi(x)\propto \textrm{e}^{-U(x)}$ , where $U(x)\in\mathcal{C}^2(\mathbb{R}^d)$ , we can exploit the information from the gradient. Consider the stochastic Langevin diffusion process

(1.2) \begin{equation} \textrm{d}X_t=-\nabla U(X_t)\,\textrm{d}t + \sqrt{2}\,\textrm{d}B_t,\end{equation}

where $B_t$ is the standard Brownian motion on $\mathbb{R}^d$ . Under certain smoothness conditions, the distribution of $X_t$ converges in probability to $\pi$ as $t\rightarrow\infty$ regardless of the initial value $X_0$ [22]. In practice, numerical machine computation requires finite precision and updates to X. We adapt the Euler–Maruyama discretization of (1.2) to offer a proper choice for transition,

(1.3) \begin{equation} X_{k+1}=X_k-h\nabla U(X_k) + \sqrt{2h}\xi_{k+1},\end{equation}

where $h>0$ is the step size and $(\xi_k,\, k \ge 1)$ is a sequence of independent and identically distributed (i.i.d.) $\mathcal{N}(0, 1)$ random variables. In other words, the proposal transition kernel $q(x, \cdot)$ is just $\mathcal{N}(x-h \nabla U(x) , 2hI_d)$ . The algorithm corresponding to (1.3) is often referred to as the unadjusted Langevin algorithm. Note that the invariant measure, say $\pi_h$ , induced by (1.3) is not identical to the target $\pi$ , so a bias must be corrected. The corresponding Metropolis–Hastings algorithm is called the Metropolized adjusted Langevin algorithm (MALA).

There have been plenty of works investigating the efficiency of MRW and MALA; see [10, 18–20, 22, 24] and references therein for more details. MRW is still popular in many applications because of its conceptual simplicity and the ease of implementation. Unfortunately, it is typically the slowest in terms of the total number of iterations, and performs poorly with increasing dimension and complexity of the target distribution.

To avoid the slow exploration of the state space that results from the diffusive behavior of simple random walk proposals, we resort to the Hamiltonian Monte Carlo (HMC), which automatically generate distant and coherent exploration for sufficiently well-behaved target distributions by carefully exploiting the differential structure of the target probability density.

1.2. Metropolized Hamiltonian Monte Carlo

Metropolized Hamiltonian Monte Carlo, abbreviated to MHMC, was introduced in [6] in computational physics, and came to the statistics community two decades later, quickly gaining popularity. The reader is referred to [2–4, 23] for nice introductory reviews.

The basic idea behind the MHMC method is as follows. Let $H(x, v)=U(x) + {\|v\|^2}/{2}$ , $(x, v)\in \mathbb{R}^d\times \mathbb{R}^d$ . We augment the target distribution $\pi\propto \textrm{e}^{-U(x)}$ to add a momentum variable v. Specifically, define the Boltzman–Gibbs probability measure $\mu_{\textrm{BG}}(x, v)\propto \textrm{e}^{-H(x, v)}$ . Obviously, the first marginal of $\mu_{\textrm{BG}}$ is $\pi$ , so if we obtain a sample from $\mu_{\textrm{BG}}$ , then upon projecting to the first coordinate we obtain a sample from $\pi$ .

Simulating $\mu_{\textrm{BG}}$ is in turn done using the Hamiltonian dynamics formulation, which is a reformulation of classical dynamics. Denote by $(x_t, v_t)$ the state at time t of a physical system, where $x_t$ is the position vector and $v_t$ the momentum vector. The evolution of the system through time is then given by the Hamiltonian equation

\begin{equation*} \frac{\textrm{d} x}{\textrm{d} t}=\nabla_vH(x_t, v_t) = v_t, \qquad \frac{\textrm{d} v}{\textrm{d} t}=-\nabla_xH(x_t, v_t) = -\nabla U(x_t).\end{equation*}

Denote by $\Phi_t\colon(x, v)\mapsto (x_t, v_t)$ the Hamiltonian flow, where (x, v) and $(x_t,v_t)$ are starting states and the states after time t, respectively. The key advantage of the Hamiltonian dynamics over other formulations of classical dynamics is that the Hamiltonian flow possesses the following fundamental properties.

Conservation of energy: Along the Hamiltonian dynamics, $(x_t, v_t)_{t\ge 0}$ satisfies $H(x_t, v_t)=H\big(x_0, v_0\big)$ , $t\ge 0$ . In fact, it is easy to check that ${\partial H(x_t, v_t)}/{\partial t} \equiv 0$ .

Conservation of volume: $\Phi_t$ is a volume-preserving map since $\det(\nabla \Phi_t(x, v))=1$ for all $t\ge 0$ and $x, v\in \mathbb{R}^d$ .

Time reversibility: Assume $(x_t, v_t)_{0\le t\le T}$ solve Hamilton’s equations, then $(x_{T-t}, -v_{T-t})_{0\le t\le T} $ also solve Hamilton’s equations.

Preservation of Boltzman–Gibbs measure: $\Phi_t$ leaves the augmented target $\mu_{\textrm{BG}}$ invariant, $\mu_{\textrm{BG}}\circ \Phi_t ^{-1}=\mu_{\textrm{BG}}$ .

However, simply running Hamilton’s equations does not yield a convergent sampling algorithm. To get around this issue we need to refresh the momentum periodically. The ideal HMC algorithm is as follows:

• Pick an integration time $T>0$ and draw $(X_0, V_0)\sim \mu_0$ , where $\mu_0$ is a fixed but arbitrary probability distribution. For each iteration $k=0, 1, 2,\ldots$

• Step 1: Refresh the momentum by drawing $\xi_{kT}\sim N(0, I_d)$ .

• Step 2: Integrate Hamilton’s equations, and set $(X_{(k+1)T}, V_{(k+1)T})=\Phi_T(X_{kT}, \xi_{kT})$ .

Since both steps of each iteration preserve $\mu_{\textrm{BG}}$ , the entire algorithm preserves $\mu_{\textrm{BG}}$ . At this stage, though, this algorithm is still idealized because it assumes the ability to exactly integrate Hamilton’s equations. This is generally not possible outside a few special cases. We approximately implement Hamilton’s equations through the use of a numerical integrator. The simplest and most well-known such integrator is the so-called leapfrog integrator.

Pick a step size $h>0$ and a total number of iterations K, corresponding to the total integration time via $T=Kh$ . Let $(x_0, v_0)$ be the initial point. For $i=0, 1, 2, \ldots, K-1$ ,

(1.4) \begin{equation}\left\{ \begin{aligned} v_{(i+{1}/{2})h} &= v_{ih} - \dfrac h2 \nabla U(x_{ih}), \\[5pt] x_{(i+1)h} &= x_{ih} + h v_{(i+{1}/{2})h}, \\[5pt] v_{(i+1)h} &= v_{(i+{1}/{2})h} - \dfrac h2 \nabla U(x_{(i+1)h}). \end{aligned}\right.\end{equation}

Let $\Phi_{h, T}(x, v)$ be the output of the leapfrog integrator with K steps started at (x, v). Once we apply the leapfrog integrator to HMC, we obtain a discrete-time sampling algorithm that is once again biased. We correct the bias through the use of the Metropolis–Hastings acceptance/rejection dynamic. The MHMC algorithm is now summarized as follows:

• Initialize at $X_0\sim \pi_0$ , where $\pi_0$ is a fixed but arbitrary probability distribution. For iterations $i=0, 1, 2, \ldots$

• Step 1: Refresh the momentum: draw $V_i\sim \mathcal{N}(0, I_d)$ .

• Step 2: Propose a trajectory: let $(X'_i, V'_i) = \Phi_{h, T}(X_i, V_i)$ .

• Step 3: Accept the trajectory with probability

  \begin{equation*} \alpha((X_i, V_i),(X_i', V_i')) = \min\bigg\{1,\frac{\exp(\!-\!H(X_i',V_i'))}{\exp(\!-\!H(X_i,V_i))}\bigg\}. \end{equation*}

If the trajectory is accepted, set $X_{i+1}=X_i'$ ; otherwise, we set $X_{i+1}=X_i$ . Iteratively, we can obtain a Markov chain $(X_i,\,i\ge0)$ with kernel $P_{h,T}$ as follows:

(1.5) \begin{equation} P_{h,T}(x, A) = \int_{\mathbb{R}^d}\boldsymbol{1}_{A\times\mathbb{R}^d}(\Phi_{h,T}(x,v))\alpha(x,v)\eta(v)\,\textrm{d}v + \delta_{x}(A)\bigg(1-\int_{\mathbb{R}^d}\alpha(x,v)\eta(v)\,\textrm{d}v\bigg),\end{equation}

with $\alpha(x,v) \;:\!=\;\alpha((x,v),\Phi_{h,T}(x,v)) = \min\{1,\exp(\!-\!H(\Phi_{h,T}(x,v)) + H(x,v))\}$ and $\eta(v) = (2\pi)^{-{d}/{2}}\exp\big({-}\frac12\|v\|^2\big)$ . It is easy to see that $\pi$ is an invariant probability measure with respect to kernel $P_{h,T}$ [3]. More properties of the HMC algorithm are detailed in Section 2.

The main purpose of this article is to construct an efficient estimator whose error variances are asymptotically negligible based on the samples generated by the MHMC algorithm.

1.3. Main results

Besides the invariance of $\pi$ , [9] showed that $P_{h,T}$ is ergodic under regular conditions, and hence, by Birkhoff’s ergodic theorem, for $\pi$ -almost every $x\in\mathbb{R}^d$ , the following limit holds for $\pi$ -integrable f: $\lim_{k\rightarrow\infty}({1}/{k})\sum_{i=0}^{k-1}f(X_i) = \mathbb{E}_\pi(f)$ , $P_x$ -a.e., where $P_x$ stands for the law with the initial state $X_0=x$ . Therefore, it is once again natural to take the mean value $({1}/{k})\sum_{i=0}^{k-1}f(X_i)$ as an approximation for $\mathbb{E}_\pi(f)$ .

In order to establish the central limit theorem (CLT) and variance reduction, we need to make some additional assumptions on the target density $\pi(x)$ , equivalently U(x), and the test function f. Motivated by [9], we introduce the following assumptions.

Assumption 1.1. For some fixed $l\in(1,2]$ , there exist $L_1>0$ and $A\in\mathbb{R}$ such that, for every $x\in\mathbb{R}^d$ , $\langle\nabla U(x), x\rangle \ge L_1\| x \|^l-A$ .

Assumption 1.2. For some fixed number $l \in (1,2]$ ,

1. (i) $U \in C^3(\mathbb{R}^d)$ and there exists a constant $L_2>0$ such that $\|\nabla^kU(x)\| \le L_2(1 + \| x \|^{l-k})$ for all $x\in\mathbb{R}^d$ and $k=2,3$ .

2. (ii) There exist $L_3>0$ and $R_0\ge 0$ such that $\langle\nabla^2 U(x)\nabla U(x),\nabla U(x)\rangle \ge L_3\|x\|^{3l-4}$ for all $x\in\mathbb{R}^d$ , $\|x\| \ge R_0$ .

Assumption 1.3. U is a perturbation of a quadratic function. In particular, $U(x) = \frac{1}{2}x^\top\Sigma x + \zeta(x)$ , where $\Sigma$ is a positive definite matrix, $\zeta\colon\mathbb{R}^d\rightarrow\mathbb{R}$ is a differentiable function, and there exist $L_4\ge 0$ and $\gamma\in[1,2)$ such that, for all $x,y\in\mathbb{R}^d$ ,

(1.6a) \begin{align} |\zeta(x)| & \le L_4(1+\|x\|)^\gamma, \\[1pt] \nonumber \end{align}

(1.6b) \begin{align} \|\nabla\zeta(x)\| & \le L_4(1+\|x\|)^{\gamma-1}, \\[1pt] \nonumber \end{align}

(1.6c) \begin{align} \|\nabla\zeta(x)-\nabla\zeta(y)\| & \le L_4\|x-y\|. \\[12pt] \nonumber \end{align}

Assumption 1.4. There exists an $r>0$ such that $|f(x)|\le V_r(x)$ , where $V_r(x)=\textrm{e}^{r\|x\|}$ , $x\in\mathbb{R}^d$ .

Remark 1.1. If Assumptions 1.1 and 1.2 hold for some $l \in (1,2]$ , then for all $T\ge 1$ and $h>0$ (we require further that $h < {M}/{T^{{3}/{2}}}$ for a certain constant $M>0$ when $l=2$ ), $P_{h,T}$ satisfies the drift condition $(V_r, \lambda, b, C)$ , i.e. there exists a set $C>0$ , $\lambda\in[0,1)$ , and $b\in(0,+\infty)$ such that

(1.7) \begin{equation} P_{h,T}V_r(x) \le \lambda V_r(x) + b {\boldsymbol{1}}_{x\in C} \quad\text{for all }x\in\mathbb{R}^d, \end{equation}

where $C=\{x\colon V_r(x)\le L_5\}$ for a constant $L_5$ . Consequently, $P_{h,T}$ is $V_r$ -uniformly geometrically ergodic.

Assumption 1.3 is stronger than Assumption 1.2, as mentioned in [9]. In fact, U in Assumption 1.3 has a special form, namely a quadratic function with a perturbation. If Assumption 1.3 holds, then (1.7) holds for all $T\ge 1$ and $0<h< {M}/{T}$ , where M is a constant. It is worth noting that the constant M can be taken small enough that the bound (2.1) holds as well.

We are now ready to state our main results.

Theorem 1.1. (Central limit theorem.)

1. (i) Fix $l \in (1,2)$ . Under Assumptions 1.1, 1.2 and 1.4, f is $\pi$ -integrable and, for each $T\ge 1$ , with $h>0$ ,

   (1.8) \begin{equation} \frac{1}{\sqrt{k}\,}\sum_{i=0}^{k-1}(f(X_i)-\mathbb{E}_\pi(f)) \stackrel{\textrm{d}}{\rightarrow} \sigma_{\pi}(f)\mathcal{N}(0,1), \quad k\rightarrow\infty, \end{equation}
   where $\sigma^2_{\pi}(f) = \textrm{Var}_\pi(f) + 2\sum_{k=1}^\infty\mathbb{E}_\pi\big[(f-\mathbb{E}_\pi(f))P_{h,T}^k(f-\mathbb{E}_\pi(f))\big] < \infty$ .

2. (ii) Fix $l=2$ . Under Assumptions 1.1, 1.2, and 1.4, f is $\pi$ -integrable and (1.8) holds for each $T\ge 1$ , $0<h<{M}/{T^{{3}/{2}}}$ with a constant $M>0$ .

3. (iii) Under Assumptions 1.1, 1.3, and 1.4, f is $\pi$ -integrable and (1.8) holds for each $T\ge 1$ , $0<h<{M}/{T}$ with a constant $M>0$ .

As a direct consequence, we can estimate the error variance $\mathbb{E}_\pi\big[({1}/{{k}})\sum_{i=0}^{k-1}(f(X_i)-\mathbb{E}_\pi(f))\big]^2$ by $k^{-1}\sigma_\pi^2(f)$ provided k is large enough, which characterizes the convergence rate of $({1}/{{k}})\sum_{i=0}^{k-1}f(X_i)$ . If $\sigma_{\pi}^2(f)$ is comparatively large, however, more iterations of MHMC are required to decrease the error of the estimation, which may reduce the efficiency of sampling. To solve such an issue, we introduce another function g(x) such that $\mathbb{E}_\pi(g)$ is known and $\mathbb{E}_\pi\big[({1}/{k})\sum_{i=0}^{k-1}((f(X_i)+g(X_i)-\mathbb{E}_\pi(g))-\mathbb{E}_\pi(f))\big]^2$ can be substantially smaller. For example, let $\hat{f}$ be the solution to the Poisson equation

(1.9) \begin{equation} \hat{f}-P_{h,T}\hat{f} = f-\mathbb{E}_\pi(f);\end{equation}

then it is easy to see that $({1}/{k})\sum_{i=0}^{k-1}(f+P_{h,T}\hat{f}-\hat{f})(X_i)$ is an unbiased estimator for $\mathbb{E}_\pi(f)$ with variance zero, which induces high accuracy. In practice, it is almost impossible to obtain a precise solution to (1.9), however. We need a proper approximation $\tilde{f}$ of $\hat{f}$ and to calculate $\mathbb\pi\big[({1}/{k})\sum_{i=0}^{k-1}(f+P_{h,T}\tilde{f}-\tilde{f})(X_i)-\mathbb{E}_\pi(f)\big]^2$ in order to determine the asymptotic variance.

Recently, [17] presented an approximation scheme for a solution of the Poisson equation of a geometrically ergodic Metropolis–Hastings chain, and further proved that the sequence of the asymptotic variance in the CLT for the control-variate estimators converged to zero. The major contribution of this article is to adapt that approximation scheme for a solution of the Poisson equation of MHMC chains to construct a sequence of control-variate estimators and to further prove the asymptotic variances converge to zero, and so realize the variance reduction. Specifically, split the whole space $\mathbb{R}^d$ into a number of subdomains, say $G_0^n, G_1^n, \ldots, G_{m_n}^n$ , and choose a suitable point $a_i^n$ inside each subdomain $G_i^n$ . Then, based on these points, construct a Markov chain with a finite number of states as a good approximation of the continuous-state Markov chain. In turn, there must exist a solution denoted by $\hat{f}_n$ to the Poisson equation induced by such a finite Markov chain. Finally, define $\tilde{f}_n(x)= \hat{f}_n(a_i^n)$ for $x\in G_i^n$ as in (3.4) and set $g_n=f+P_{h, T}\tilde{f}_n-\tilde{f}_n$ , which is a perturbation of the original test function f. Note that $E_\pi g_n=E_\pi f$ , and we have the following variance reduction theorem.

Theorem 1.2 (Variance reduction.) In the above setting and with the assumptions in Theorem 1.1, the central limit theorem holds for every $n\ge 1$ :

\begin{equation*} \frac{1}{\sqrt{k}\,}\sum_{i=0}^{k-1}(g_n(X_i)-\mathbb{E}_\pi(f)) \stackrel{\textrm{d}}{\longrightarrow} \sigma_{\pi}(g_n)\mathcal{N}(0,1), \quad k\rightarrow\infty, \end{equation*}

where $\sigma^2_{\pi}(g_n) = \textrm{Var}_\pi(g_n) + 2\sum_{k=1}^\infty\mathbb{E}_\pi\big[(g_n-\mathbb{E}_\pi(f))P_{h,T}^k(g_n-\mathbb{E}_\pi(f))\big]<\infty$ . Moreover, $\lim_{n\rightarrow\infty} \sigma^2_{\pi}(g_n)=0$ .

The rest of the paper is organized as follows. In Section 2, we first review some basic concepts and notions about Markov chains with general state spaces, and then give some more technical results about the MHMC Markov chains. In Section 3 we formulate the approximation scheme based on the idea of weak approximation following [17], and construct in a specific way $\tilde{f}_n$ , which is used in the variance reduction. Section 4 is devoted to the detailed proofs of the main results. The proof of Theorem 1.1 is actually standard through the use of Lyapunov’s drift conditions. However, Theorem 1.2 cannot be obtained by a simple combination of [9, 17]. In fact, there is an additional term $|\textrm{det}\,\textbf J_{\Phi_x}(\bar{x})|$ in the HMC kernel, which is so complex that it does not satisfy the assumption in [17]. Our proof makes delicate use of the specific approximate construction and moment controls. Section 5 is devoted to two concrete examples, a Bayesian linear inverse problem and a two-component mixture of Gaussians, which satisfy the conditions for variance reduction. In Section 6, we provide some numerical experiments to support our theoretical guarantees. The paper concludes with a discussion in Section 7.

1.4. Notation

Denote by $\mathbb{R}^{+}$ the set of positive real numbers. Denote by $\mathbb{N}^+$ the set of positive integers. Denote by $\|\cdot\|$ the Euclidean norm on $\mathbb{R}^{d}$ . Denote by $\mathcal{B}(\mathbb{R}^{d})$ the Borel $\sigma$ -field of $\mathbb{R}^{d}$ and $\mathsf{F}(\mathbb{R}^{d})$ the set of all Borel-measurable functions on $\mathbb{R}^{d}$ . Denote by $\mu(\cdot)$ the Lebesgue measure on $\mathbb{R}^{d}$ . For a Markov kernel P, denote by $(X_i)_{i\in\mathbb{N}}$ the corresponding canonical Markov chain, and denote by $(\Omega,\mathcal{F})$ the canonical space. For any probability measure $\nu$ on $\mathbb{R}^d$ and a Markov kernel P, denote by $\mathbb{P}_{\nu}$ the unique probability measure on the canonical space $(\Omega,\mathcal{F})$ such that the canonical process $X_i$ is a Markov chain with kernel P and initial distribution $\nu$ . For $f \in \mathsf{F}(\mathbb{R}^{d})$ , set $\|f\|_{\infty} =\sup _{x \in \mathbb{R}^{d}}|f(x)|$ . For $f \in \mathsf{F}(\mathbb{R}^{d})$ and $V\colon\mathbb{R}^d\rightarrow[1,\infty)$ , the V-norm of f is given by $\|f\|_{V}=\|f / V\|_{\infty}$ . Suppose M is a $p \times q$ matrix; denote by $M^\top$ and $\det (M)$ the transpose and the determinant of M, respectively. Denote by $C^{k}(\mathbb{R}^d)$ the set of all k-times continuously differentiable functions. Let $f \in C^{k}(\mathbb{R}^p, \mathbb{R}^{q})$ and denote by $\textbf{J}_{f}$ the Jacobian matrix of f. Let $g \in C^{k}(\mathbb{R}^p, \mathbb{R})$ and denote by $\nabla g$ , $\nabla^{2} g$ , and $\nabla^{3} g$ the first-, second-, and third-order derivatives of g, respectively.

2. Properties of HMC

In this section we introduce important properties of HMC that will be used later. The reader is also referred to [2, 3, 23] for more details. For the reader’s convenience, we briefly review some basic concepts and notions about Markov chains on $\mathbb{R}^d$ in Appendix A.

Suppose that the target density is $\pi(x)\propto \textrm{e}^{-U(x)}$ , with U satisfying Assumptions 1.1–1.3. Thanks to the leapfrog algorithm and Metropolis–Hastings step in the MHMC, the resulting Markov chain $(X_i,\,i\geq 0)$ with kernel $P_{h,T}$ is reversible with respect to $\pi(x)$ , and thus $\pi$ is an invariant probability measure [14].

As for the irreducibility and ergodicity of MHMC Markov chains, additional conditions are required. Recently, [9] provided some sufficient conditions for the irreducibility of MHMC Markov chains. Indeed, suppose that Assumption 1.2 holds and that both the step size h and step number T satisfy the inequality

(2.1) \begin{equation} \big[1+\vartheta\big(hL_2^{{1}/{2}}\big)\big]^{T}<2,\end{equation}

where $\vartheta(x) = x\big(1+\frac{1}{2}x+\frac{1}{4}x^2\big)$ ; then, for all $x\in\mathbb{R}^d$ there exists a $C^1(\mathbb{R}^d,\mathbb{R}^d)$ diffeomorphism $\Psi_x\colon\bar{x}\mapsto \Psi_x(\bar{x})$ such that $x_T=\textrm{proj}_1\circ\Phi_{h,T}(x,\Psi_x(x_T))$ , where $\textrm{proj}_1\colon(x,y)\mapsto x$ denotes coordinate projection and $\Phi_{h,T}$ is defined by (1.4).

Let $\textbf{J}_{\Psi_{x} }(\bar{x})$ be the Jacobian matrix of $\Psi_x$ at the point $\bar{x}$ . The following lemma is important when proving the main theorems.

Lemma 2.1. If Assumption 1.2 and condition (2.1) hold, then there exists a constant $\kappa(h,T)$ such that, for any $x,\bar{x}\in\mathbb{R}^d$ , $D(x,\bar{x}) \;:\!=\; |\det\textbf{J}_{\Psi_{x}}(\bar{x})| \le \kappa(h,T)$ .

Proof. By the structure of the leapfrog integrator in (1.4), for all $(x_0,v_0)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ and $t\in\{1,2,\ldots,T\}$ , the tth iteration $(x_t,v_t) = \Phi_{h,t}(x_0,v_0)$ takes the form

(2.2) \begin{align} x_t & = x_0 + thv_0 - \frac{1}{2}th^2\nabla U(x_0) - h^2\sum_{k=1}^{t-1}(t-k)\nabla U(x_k), \\[5pt] v_t & = v_0 - \frac{1}{2}h(\nabla U(x_{0}) + \nabla U(x_t)) - h\sum_{k=1}^{t-1}\nabla U(x_k). \notag \end{align}

Let $\Gamma_{h,t}(x,v) =\sum_{k = 1}^{t-1}(t-k)\nabla U(x_k)$ . Then, for any $t\ge 1$ and $h>0$ , Assumption 1.2 implies

\begin{equation*} \sup_{x,v,w\in\mathbb{R}^{d}}\frac{\|\Gamma_{h,t}(x,v)-\Gamma_{h,t}(x,w)\|}{\|v-w\|} \le \frac{t}{h}(\kappa_{h,t}-1), \end{equation*}

where $\kappa_{h,t} = \big(1+hL_2^{{1}/{2}}\vartheta\big(hL_2^{{1}/{2}}\big)\big)^t <2$ [9]. Therefore, by (2.1), we have

(2.3) \begin{equation} \sup_{x,v\in\mathbb{R}^{d}}\frac{h}{t}\big\|\textbf{J}_{\Gamma_{h,t}}(x,v)\big\| < 1, \end{equation}

where $\textbf{J}_{\Gamma_{h,t}}(x,v)$ is the Jacobian matrix of the function $v\mapsto \Gamma_{h,t}(x, v)$ .

As a direct consequence, the function $v_0 \mapsto x_t$ defined in (2.2) is a diffeomorphism and has a continuously differentiable inverse $\bar{x}\mapsto \Psi_x(\bar{x})$ for each $x\in\mathbb{R}^d$ . In addition, by (2.2) and (2.3), there exists a constant $\kappa_1\in (0,1)$ such that, for any $x,v,w\in\mathbb{R}^d$ ,

\begin{multline*} \|\textrm{proj}_1 \circ \Phi_{h,T}(x,v) - \textrm{proj}_1 \circ \Phi_{h,T}(x,w)\|\\= hT\bigg\|\bigg(v-\frac{h}{T}\Gamma_{h,T}(x,v)\bigg) - \bigg(w-\frac{h}{T}\Gamma_{h,T}(x,w)\bigg)\bigg\| \ge hT\kappa_1\|v-w\|, \end{multline*}

which in turn implies that there exists a constant $\kappa_2>0$ such that, for all $x,v,w\in\mathbb{R}^d$ ,

\begin{equation*} \| \Psi_{x}(v) - \Psi_x(w)\| \le \kappa_2 \| v -w \|. \end{equation*}

Therefore, for all $x,\bar{x}\in\mathbb{R}^d$ , it follows from Hadamard’s inequality that

\begin{equation*} |\det\textbf{J}_{\Psi_{x}}(\bar{x})| \le \|\textbf{J}_{\Psi_{x}}(\bar{x})\|^d \le \sup_{x,\bar{x}\in\mathbb{R}^{d}}\|\textbf{J}_{\Psi_{x}}(\bar{x})\|^d \le \kappa_2^d \;=\!:\; \kappa(h,T) > 0, \end{equation*}

where the last inequality is due to [7, Exercise 2.24].

Lemma 2.2. If Assumption 1.2 and the condition (2.1) hold, then for any $T\ge 1$ , $h\ge 0$ , and $(x,v)\in\mathbb{R}^d\times\mathbb{R}^d$ , $(x,v)\mapsto(\Phi_{h,T}(x,v),\Psi_x(v),D_{\Psi_x(\cdot)}(v))$ is continuous on $\mathbb{R}^d\times\mathbb{R}^d$ .

Proof. Continuity for the mapping $(x,v)\mapsto(\Phi_{h,T}(x,v),\Psi_x(v),D_{\Psi_x(\cdot)}(v))$ was presented in the proof of [9, Theorem 1].

Remark 2.1. The assumptions of Lemma 2.2 seem slightly different from those in Theorem 1.2. However, by adjusting the value of M to make h small enough, condition (2.1) can be satisfied and thus Lemma 2.2 remains valid.

It now follows from (1.5) that

(2.4) \begin{align} P_{h,T}(x,\textrm{d}\bar{x}) & = \alpha_{\textrm{H}}(x,\bar{x})\frac{\exp\{-\|\Psi_x(\bar{x})\|^2/2\}}{(2\pi)^{{d}/{2}}} D(x,\bar{x})\,\textrm{d}\bar{x} \nonumber \\[5pt] & \quad + \delta_{x}(\textrm{d}\bar{x})\int_{\mathbb{R}^{d}}(1-\alpha(x,v)) \frac{\exp\{-\|v\|^{2}/2\}}{(2\pi)^{\frac{d}{2}}}\,\textrm{d}v, \end{align}

where $\alpha_{\textrm{H}}(x,\bar{x})=\alpha(x,\Psi_x(\bar{x}))$ , $D(x,\bar{x})$ was defined in Lemma 2.1.

Furthermore, the Markov kernel $P_{h,T}$ is $\mu$ -irreducible, aperiodic, and Harris recurrent, and each compact set is 1-small. As a consequence, for all $x\in\mathbb{R}^d$ , $\lim_{k\rightarrow \infty}\big\|P^k_{h,T}(x,\cdot)-$ $\pi(\cdot)\big\|_{\textrm{TV}}=0$ . Note also that if a Markov chain has a unique invariant probability measure, then ergodicity is valid and the ergodic theorem can be obtained by [5, Theorem 5.2.6]. Therefore, by the irreducibility and [5, Theorem 9.2.15], $\pi$ is the unique invariant probability measure with respect to $P_{h,T}$ and hence, for all $f\in L^1(\pi)$ and $\pi$ -almost every $x\in\mathbb{R}^d$ , we have $\lim_{k\rightarrow\infty}({1}/{k})\sum_{i=0}^{k-1}f(X_i) = \mathbb{E}_\pi(f)$ , ${P}_x\text{-a.e.}$ This forms the starting point of our further study in this work.

3. Approximation scheme for variance reduction

Definition 3.1.

1. (i) Assume there exists a partition $\mathbb{G}$ of $(\mathbb{R}^d,\mu)$ into measurable subsets $G_0,G_1,\ldots,G_m$ such that $\mu(G_i)>0$ holds for $0\le i\le m$ and $\bigcup_{i=1}^m G_i$ is bounded. Choose $a_i\in G_i$ for all $0\le i \le m$ and let $Y =\{a_0,\ldots,a_m\}$ . We say that the pair $\mathbb{Y} =(\mathbb{G},Y)$ is an allotment, with m being the size of $\mathbb{Y}$ .

2. (ii) Let $W\colon\mathbb{R}^d\rightarrow[1,\infty)$ be a measurable function. The W-radius and W-mesh of an allotment $\mathbb{Y}$ are respectively defined by $\textrm{rad}(\mathbb{Y},W)=\inf_{y\in G_0} W(y)$ and

   \begin{equation*} \delta(\mathbb{Y},W) = \max\bigg\{\max_{1\le i\le m}\sup_{y\in G_i}|y-a_i|,\ \max_{0\le i\le m}\sup_{y\in G_i}\bigg(\frac{W(a_i)}{W(y)}-1\bigg)\bigg\}. \end{equation*}

3. (iii) If there exists a sequence of allotments $(\mathbb{Y}_n,n\ge 0)$ satisfying $\lim_{n\rightarrow\infty}\textrm{rad}(\mathbb{Y}_n,W) = \infty$ and $\lim_{n\rightarrow\infty}\delta(\mathbb{Y}_n,W)=0$ , then we say that $\mathbb{Y}_n$ is exhaustive with respect to W.

Consider $V_r(x)=\textrm{e}^{r\|x\|}$ , $x\in \mathbb{R}^d$ . Obviously, $V_r\colon\mathbb{R}^d\rightarrow[1,\infty)$ is a continuous function with bounded sublevel sets, i.e. for every $c\in \mathbb{R}$ , the pre-image $V_r^{-1}((\!-\!\infty, c])$ is a bounded ball. Proposition A.1 in [17] establishes the existence of an exhaustive sequence with respect to $V_r$ . For the sake of completeness, we briefly review the construction of such an allotment here.

Let $r_1>1$ , and let $(r_n, \,n\ge 1)$ be an increasing unbounded sequence of positive numbers. For each $n\ge 1$ , define sets $L_n=V_r^{-1}((\!-\!\infty, r_n])$ , $\tilde{L}_n=\{x\in\mathbb{R}^d\colon\text{there exists}\ y\in L_n\mbox{ such that }\|x-y\|<\sqrt d\}$ . The set $\tilde{L}_n$ is trivially a bounded and non-empty closed set. $V_r$ is uniformly continuous on $\tilde{L}_n$ . So there exists a decreasing sequence $\varepsilon_n<1$ , vanishing as n tends to infinity, such that $\|x-y\|<\varepsilon_n \sqrt{d}$ implies $|V_r(x)-V_r(y)|< 1/n$ for each n and all $x,y \in \tilde{L}_n$ .

Fix $n\ge 1$ . For $x=(x_1,x_2,\ldots,x_d)\in\mathbb{R}^d$ write $K_x^n \;:\!=\; \prod_{i=1}^d[x_i, x_i+\varepsilon_n)$ . Pick points $x^{(1)},x^{(2)},\ldots,x^{(m_n)}$ such that the sets $G_j^n\;:\!=\;K^n_{x^{(j)}}$ are disjoint and cover $L_n$ . Let $G_0^n$ be the closure of $\mathbb{R}^d\setminus \bigcup_{j=1}^{m_n}K^n_{x^{(j)}}$ . Finally, choose $a^n_0\in G_0^n$ such that $V_r(a^n_0)=\inf_{x\in G_0^n}V_r(x)$ , and pick $a^n_j\in G^n_j$ arbitrarily. It is now easy to verify that the allotment defined above is exhaustive with respect to $V_r$ .

Let $\mathbb{Y}_n=(\mathbb{G}^n,Y^n)$ , where $\mathbb{G}^n=(G_i^n,\,0\le i\le m_n)$ and $Y^n=(a_i^n,\, 0\le i\le m_n)$ . We will construct a sequence of approximating solutions $\tilde{f}_n$ based on the exhaustive allotments $\mathbb{Y}_n$ . Given the kernel $P_{h,T}$ with $h,T>0$ , define $Q^{(n)}_{ij} = P_{h,T}\big(a^n_i,G^n_j\big)$ and $Q^{(n)}= \big(Q^{(n)}_{ij}\big)_{(m_n+1)\times (m_n+1)}$ . Since each part $G^n_i$ has a positive measure the transition matrix $Q^{(n)}$ is well defined, and for every $i,j\in\{0,1,\ldots,m_n\}$ , $Q^{(n)}_{ij}>0$ by the $\mu$ -irreducibility of $P_{h,T}$ . Furthermore, $Q^{(n)}$ is irreducible, aperiodic, and recurrent, so there exists a unique invariant probability measure by Markov chain theory on a finite state space. Let $a_n(x)=\sum_{i=0}^{m_n}a_i^n \boldsymbol{1}_{x\in G_i^n}$ , $x\in\mathbb{R}^d$ . Then, by the definitions of $V_r$ -mesh and exhaustive allotment, for every $x\in\mathbb{R}^d$ we have $\lim_{n\rightarrow\infty}a_n(x)=x$ . Observe that $V_r(a_n(x))$ can be controlled by $V_r(x)$ . Indeed, for all $n\ge 1$ and $x\in\mathbb{R}^d$ , we have the useful inequality

(3.1) \begin{align} V_r(a_n(x)) & = V_r(x)\bigg(1+\frac{V_r(a_n(x))-V_r(x)}{V_r(x)}\bigg) \nonumber \\[5pt] & \le V_r(x)\bigg(1+\max_{0\le i\le m}\sup_{y\in G_i}\bigg(\frac{V_r(a_n(x))}{V_r(y)}-1\bigg)\bigg) \le V_r(x)(1+\delta_n), \end{align}

where $\delta_n=\delta(\mathbb{Y}_n, V_r)$ denotes the $V_r$ -mesh of the allotment $\mathbb{Y}_n$ , the first inequality following by the exhaustivity property.

Lemma 3.1

1. (i) Drift condition. There exist constants $\lambda_v\in(0,1)$ , $b_v>0$ satisfying

   (3.2) \begin{equation} Q^{(n)} V_r(a_i^n) \le \lambda_v V_r(a_i^n) + b_v\boldsymbol{1}_{a_i^n\in C} \end{equation}
   for all $n\ge 1$ and $a_i^n\in Y_n$ .

2. (ii) Minorization condition. There exist a compact set $C\subset\mathbb{R}^d$ , a constant $\varepsilon\in\mathbb{R}^+$ , and a probability measure $\nu_n$ on $(Y_n, \mathcal{P}\big(Y_n\big))$ such that

   (3.3) \begin{equation} Q^{(n)}_{ij} \ge \varepsilon\nu_n(\{a_j^n\}) \end{equation}
   for all $n\ge 1$ and $i, j \in \{0,1,\ldots,m_n\}$ satisfying $a_i^n\in Y_n\cap C$ .

3. (iii) Strong aperiodicity condition. There exists a constant $\bar{\varepsilon}\in (0,\infty)$ such that $\varepsilon \nu_n(C\cap Y_n)\ge \bar{\varepsilon}$ .

The proof of Lemma 3.1 can be found in Appendix B.

Let $\hat{\pi}_n$ be an invariant probability measure of $Q^{(n)}$ .

Proposition 3.1. The approximate Markov kernel $Q^{(n)}$ is $V_r$ -uniformly geometrically ergodic, i.e. there exist constants $M>0$ and $\rho\in(0,1)$ such that, for $k, n\ge 1$ ,

\begin{equation*} \big\|\big(Q^{(n)}\big)^k(y,\cdot)-\hat{\pi}_n(\cdot)\big\|_{V_r} \le M V_r(y)\rho^k \quad \text{for all }y \in Y_n. \end{equation*}

Proof. By Lemma 3.1, the existence of M and $\rho$ is a direct consequence of [1]. It is also implied that M and $\rho$ only depend on the parameters $\varepsilon$ , $\lambda_v$ , $b_v$ , and $\bar{\varepsilon}$ in Lemma 3.1.

We are now in a position to construct the approximation solution $\tilde{f}_n$ . To this end, we denote by $\hat{f}_n$ a solution to the Poisson equation $\hat{f}_n(y) - Q^{(n)}\hat{f}_n(y)=f(y) - \mathbb{E}_{\hat{\pi}_n}(f)$ , $y\in Y_n$ . Having $\hat{f}_n$ , the approximating solution $\tilde{f}_n$ is defined by

(3.4) \begin{equation} \tilde{f}_n(x)= \sum_{i=0}^{m} \hat{f}_n\big(a^n_i\big)\boldsymbol{1}_{x\in G^n_i},\quad x\in\mathbb{R}^d.\end{equation}

Obviously, the more precisely $Q^{(n)}$ approximates $P_{h,T}$ , the more closely $\tilde{f}_n$ approaches the exact solution of the Poisson equation.

4. Proofs of main results

Throughout this section, the parameter $r>0$ is fixed. We start with the proof of Theorem 1.1.

Proof of Theorem 1.1. We only prove (i) with $l\in(1, 2)$ , since the other cases are similar. Recall that $P_{h,T}$ satisfies the drift condition (1.7) with $V_r(x) = \textrm{e}^{r\|x\|}$ . By [5, Theorem 15.2.4], $P_{h,T}$ is both $V_r$ and $V_{2r}$ -uniformly geometrically ergodic, and so it follows by [21, Fact 10] that $V_r$ and $V_{2r}$ are $\pi$ -integrable.

In addition, we obviously have the following equivalence relations:

\begin{equation*} \begin{array}{lr@{\ }c@{\ }l} & P_{h,T}V_r & \le & \lambda V_r + b \boldsymbol{1}_C \\[5pt] \Leftrightarrow & P_{h,T}V_r + (1- \lambda)V_r & \le & V_r + b \boldsymbol{1}_C \\[5pt] \Leftrightarrow & P_{h,T}\dfrac{V_r}{1-\lambda} + V_r & \le & \dfrac{V_r}{1-\lambda} + \dfrac{b}{1-\lambda} \boldsymbol{1}_C, \end{array} \end{equation*}

where $C=\{x\colon V_r(x)\le L_2\}$ . Therefore, for every f such that $\|f\|_{V_r}<\infty$ , f is $\pi$ -integrable, and further the CLT (1.8) holds according to [5, Theorem 21.2.11].

Next, we turn to the proof of Theorem 1.2. We first give an upper bound for error variance in terms of the spectral radius of $P_{h,T}$ on $L^2_0(\pi)$ .

Proposition 4.1. Denote by $\rho$ the spectral radius of $P_{h,T}|_{L_0^2(\pi)}$ . For any $g\in L_0^2(\pi)$ ,

$$ \sigma_\pi^2(g)\le \frac{1+\rho}{1-\rho}\mathbb{E}_\pi(g^2). $$

The proof of Proposition 4.1 can be found in Appendix D.

For notational simplicity, set $g_n =f +P_{h,T}\tilde{f}_n-\tilde{f}_n$ , where $\tilde{f}_n$ is given by (3.4). Trivially, $\mathbb{E}_\pi(g_n)=\mathbb{E}_\pi(f)$ , and by Proposition 4.1,

$$\sigma_\pi^2(g_n)\le \frac{1+\rho}{1-\rho}\mathbb{E}_\pi(g_n-\mathbb{E}_\pi(f))^2.$$

Therefore, it suffices to prove that $\mathbb{E}_\pi(g_n-\mathbb{E}_\pi(f))^2\rightarrow 0$ as $n\rightarrow\infty$ . In turn, it is sufficient to verify the following statements due to the dominated convergence theorem:

• • For $\pi$ -almost every $x\in\mathbb{R}^d$ , $\lim_{n\rightarrow\infty} g_n(x) = \mathbb{E}_\pi(f)$ .

• • The $V_r$ -norm is uniformly bounded: $\sup_{n\ge 1}\|g_n-\mathbb{E}_\pi(f)\|_{V_r} < \infty$ .

For clarity, we restate these statements as Propositions 4.2 and 4.3.

Proposition 4.2 For $\pi$ -almost every $x\in\mathbb{R}^d$ , $\lim_{n\rightarrow\infty} g_n(x)=\mathbb{E}_\pi(f)$ .

Observe that for any $x\in\mathbb{R}^d$ , $\tilde{f}_n(a_n(x))-P_{h,T}\tilde{f}_n(a_n(x)) = f(a_n(x))-\mathbb{E}_{\hat{\pi}_n}(f)$ , and so

\begin{align*} |g_n(x)-\mathbb{E}_\pi f| & \le |f(x)-f(a_n(x))| + |\mathbb{E}_{\hat{\pi}_n}(f)-\mathbb{E}_\pi(f)| \\[5pt] & \quad + |(P_{h,T}\tilde{f}_n-\tilde{f}_n)(x) - (P_{h,T}\tilde{f}_n-\tilde{f}_n)(a_n(x))| \\[5pt] & \;=\!:\; M^{(1)}_n(x)+M_n^{(2)} +M_n^{(3)}(x).\end{align*}

Then it reduces to verifying that $M_n^{(1)}$ , $M_n^{(2)}$ , and $M_n^{(3)}$ converge $\pi$ -a.e. to zero as $n\rightarrow\infty$ , hence the following three lemmas.

Lemma 4.1. $\lim_{n\rightarrow\infty}M_n^{(1)}(x)=0$ for $\pi$ -almost every x.

Proof. Note that f is $\pi$ -a.e. continuous and $a_n(x)\rightarrow x$ as $n\rightarrow\infty$ . This concludes the proof.

Lemma 4.2. $\lim_{n\rightarrow\infty}M_n^{(2)}=0$ .

The proof of Lemma 4.2 can be found in Appendix C. The following lemma offers a uniform bound for the $V_r$ -norm of $\tilde{f}_n$ .

Lemma 4.3. There exist a constant $\beta>0$ and a sequence of real numbers $\{b_n\}_{n\in\mathbb{N}}$ such that $\sup_{n\ge 1}\|\tilde{f}_n+b_n\|_{V_r}\le \beta$ .

Proof. This is a direct consequence of Proposition 3.1 combined with [17, Proposition 3.5].

Lemma 4.4. $\lim_{n\rightarrow\infty}M_n^{(3)}(x)=0$ for $\pi$ -almost every x.

Proof. The form (2.4) and the structure of $P_{h,T}$ together yield

\begin{equation*} \begin{aligned} M_n^{(3)} & = |P_{h,T}\tilde{f}_n(x)-P_{h,T}\tilde{f}_n(a_n(x))| \\[5pt] & = \bigg|\int_{\mathbb{R}^{d}}(\tilde{f}_n(\bar{x})-\tilde{f}_n(x)) \big(\alpha_{\textrm{H}}(x,\bar{x})D(x,\bar{x})\textrm{e}^{-\|\Psi_x(\bar{x})\|^2/2} \\[5pt] & \qquad\quad - \alpha_{\textrm{H}}(a_n(x),\bar{x})D(a_n(x),\bar{x}) \textrm{e}^{-\|\Psi_{a_n(x)}(\bar{x})\|^2/2}\big)\,\textrm{d}\bar{x}\bigg|, \end{aligned} \end{equation*}

where in the first equation we used the fact that $\tilde{f}_n(a_n(x))=\tilde{f}_n(x)$ .

Since $(x,\bar{x})\mapsto(D(x,\bar{x}), \alpha_{\textrm{H}}(x,\bar{x}))$ is continuous by Lemma 2.2 and, for each $\bar{x}\in\mathbb{R}^d$ ,

(4.1) \begin{equation} \big|\tilde{f}_n(\bar{x}) - \tilde{f}_n(x)\big| \le \big|\tilde{f}_n(\bar{x}) + b_n - (\tilde{f}_n(x) + b_n)\big| \le \beta(V_r(\bar{x}) + V_r(x)), \end{equation}

where $b_n$ and $\beta$ were given in Lemma 4.3, the integrand above converges to zero.

Analogous to the proof of Lemma 4.2, we obtain, for each $\bar{x}\in\mathbb{R}^d$ ,

(4.2) \begin{multline} \big|\alpha_{\textrm{H}}(x,\bar{x})D(x,\bar{x})\textrm{e}^{-\|\Psi_x(\bar{x})\|^2/2} - \alpha_{\textrm{H}}(a_n(x),\bar{x})D(a_n(x),\bar{x})\textrm{e}^{-\|\Psi_{a_n(x)}(\bar{x})\|^2/2}\big| \\[5pt] \le \delta_x\pi(\bar{x}) + \alpha_{\textrm{H}}(x,\bar{x})k(x,\bar{x}). \end{multline}

Obviously, the product of the right-hand sides in (4.1) and (4.2) is integrable with respect to $\mu$ over $\mathbb{R}^d$ . Therefore, $M_n^{(3)}$ tends to zero by the dominated convergence theorem.

Proposition 4.3. $\sup_{n\ge 1}\|g_n-\mathbb{E}_\pi(f)\|_{V_r} < \infty$ .

Proof. By the drift condition of $P_{h,T}$ , the solution $\hat{f}$ to the Poisson equation (1.9) belongs to $L_{V_r}^\infty=\{f\colon\|f\|_{V_r}<\infty\}$ , and it is unique up to a constant [5]. Therefore, by Lemma 4.3, there exist a constant $\beta'$ and a sequence of real numbers $\{b_n\}_{n\ge 1}$ such that, for all $n\geq 1$ and $x\in\mathbb{R}^d$ , $\big|\tilde{f}_n(x)+b_n-\hat{f}(x)\big|\le\beta'V_r(x)$ . Observe that

$$ g_n - \mathbb{E}_\pi(f) = P_{h,T}\tilde{f}_n-\tilde{f}_n+f-\pi(f) = P_{h,T}\big(\tilde{f}_{n}+b_{n}-\hat{f}\big)-\big(\tilde{f}_{n}+b_{n}-\hat{f}\big). $$

By the definition of $P_{h,T}$ , it easily follows that, for all $n\ge 1$ and $x\in\mathbb{R}^d$ ,

\begin{equation*} \begin{aligned} |g_n - \mathbb{E}_\pi(f)| & = \bigg|\int_{\mathbb{R}^d}\big(\tilde{f}_n+b_n-\hat{f}\big)(\Phi_{h,T}(x,v))\alpha(x,v)\eta(v)\,\textrm{d}v \\[5pt] & \qquad\quad + \big(\tilde{f}_n+b_n-\hat{f}\big)(x)\cdot \bigg[1-\int_{\mathbb{R}^d}\alpha(x,v)\,\eta(\textrm{d}v)\bigg]-\big(\tilde{f}_{n}+b_{n}-\hat{f}\big)(x)\bigg| \\[5pt] & \le \int_{\mathbb{R}^d}\big[\big|\big(\tilde{f}_n+b_n-\hat{f}\big)(\Phi(x,v))\big| + \big|\big(\tilde{f}_n+b_n-\hat{f}\big)(x)\big|\big]\alpha(x,v)\eta(v)\,\textrm{d}v \\[5pt] & \le \int_{\mathbb{R}^d}\beta'V_r(\Phi_{h,T}(x,v))\alpha(x,v)\,\eta(\textrm{d}v) + \beta'V_r\int_{\mathbb{R}^d}\alpha(x,v)\eta(v)\,\textrm{d}v \\[5pt] & \le \beta'(P_{h,T}V_r(x)+V_r(x)) \leq (1+\lambda+b)\beta'V_r(x), \end{aligned} \end{equation*}

which together with the definition of the $V_r$ -norm implies the desired result.

Proof of Theorem 1.2. Combining Propositions 4.1, 4.2, and 4.3 concludes the proof.

5. Applications

In this section we aim to show how to verify the assumptions on the target distributions given in Section 1.3 through two concrete models.

5.1. Bayesian linear inverse problem

Let $\beta\in\big(\frac12,1\big)$ , $\lambda_1, \lambda_2, \delta>0$ , and $\textbf{A}$ a $p\times d$ matrix. Consider the linear model $b=\textbf{A}X+\varepsilon$ , where X has a prior distribution $\pi_X(x)$ given by

\[ \pi_X(x) \propto \exp\bigg({-}\lambda_1\big(x^\top x+\delta\big)^\beta - \frac{\lambda_2}{2}x^\top x\bigg), \quad x\in\mathbb{R}^d,\]

and $\varepsilon\sim N(0, I_d)$ .

The so-called Bayesian linear problem is to infer X from the observations b. In fact, the density of the posterior distribution of interest us given by $\pi(x)\propto\exp(\!-\!U(x))$ , where

$$U(x) = -\dfrac{1}{2}x^{\top}\big(\textbf{A}^{\top }\textbf{A}+\lambda_{2}I_{d}\big)x -\lambda_{1}\big(x^{\top}x+\delta\big)^{\beta} + \langle b, \textbf{A}x \rangle.$$

The exponential integrator version of MALA (EI-MALA) was devised to generate a Markov chain with invariant distribution $\pi$ . The Wasserstein convergence rate of a class of EI-MALA algorithms was studied at length in [8], and the CLT for the corresponding Markov chains established in [13].

We will confirm that the potential U satisfies Assumptions 1.1 and 1.2, and so both the CLT and the variance reduction theorem are applicable for the MHMC Markov chains as well. It is apparent that U(x) can be decomposed as

(5.1) \begin{equation} U(x) = U_1(x) + U_2(x) = \bigg[\frac{1}{2}x^{\top}\big(\textbf{A}^{\top}\textbf{A}+\lambda_{2}I_{d}\big)x\bigg] + \bigg[\lambda_{1}\big(x^{\top}x+\delta\big)^{\beta}-\langle b,\textbf{A}x\rangle\bigg]. \end{equation}

[9] offers a sufficient condition for Assumptions 1.1 and 1.2. The following is a slightly stronger version.

Lemma 5.1. Fix $l\in(1,2]$ . Suppose that the potential U(x) can be decomposed as $U(x) = U_1(x)+U_2(x)$ , where $U_1(x)$ and $U_2(x)$ satisfy the following conditions:

1. (i) $U_1$ and $U_2$ belong to $C^3(\mathbb{R}^d)$ .

2. (ii) For all $k\ge1$ and $x\in\mathbb{R}^d$ , $U_1(kx) = k^lU_1(x)$ and $\{y\in\mathbb{R}^d\colon U_1(y)\le U_1(x)\}$ is a convex set.

3. (iii) $\lim_{\|x\|\rightarrow\infty}U_1(x)=\infty$ .

4. (iv) For $k=2,3$ ,

   (5.2) \begin{equation} \lim\limits_{\|x\|\rightarrow\infty}\dfrac{\|\nabla^k U_2(x)\|}{\|x\|^{l-k}} = 0. \end{equation}

Then the potential U satisfies Assumptions 1.1 and 1.2.

Proposition 5.1. The potential U given by (5.1) satisfies Assumptions 1.1 and 1.2.

The proof of Proposition 5.1 can be found in Appendix E.

5.2. Two-component mixture of Gaussians

Consider the sampling from a mixture of Gaussians with two components

(5.3) \begin{equation} \pi \sim p\mathcal{N}(\mu,I_d) + (1-p)\mathcal{N}(\nu,I_d),\end{equation}

where $\mathcal{N}(\cdot, I_d)$ is a d-dimensional normal distribution with identity variance.

Without loss of generality, assume $p = {\textrm{e}^c}/({\textrm{e}^c+\textrm{e}^{-c}})$ for some $c\in\mathbb{R}$ , and let $U_1(x)=\frac{1}{2}\|x \|^2$ , $U_2(x) =- \log\big(\textrm{e}^{c-\|\mu\|^2/2}\textrm{e}^{-\langle x,\mu\rangle} +\textrm{e}^{-c-\|\nu\|^2/2}\textrm{e}^{-\langle x,\nu\rangle}\big)$ . Then

\begin{align*} \pi(x) & = \frac{1}{(2\pi)^{{d}/{2}}}\bigg[p\exp\bigg\{{-}\frac{1}{2}\|x - \mu\|^2\bigg\} + (1-p)\exp\bigg\{{-}\frac{1}{2}\|x-\nu\|^2\bigg\}\bigg] \\[5pt] & = \frac{1}{(2\pi)^{{d}/{2}}(\textrm{e}^c+\textrm{e}^{-c})}\textrm{e}^{-\|x\|^2/2} \big[\textrm{e}^{c-\|\mu\|^2/2}\textrm{e}^{-\langle x,\mu\rangle} + \textrm{e}^{-c-\|\nu\|^2/2}\textrm{e}^{-\langle x,\nu\rangle}\big] \\[5pt] & \propto \textrm{e}^{-U(x)} = \textrm{e}^{-(U_1(x) + U_2(x))}.\end{align*}

While there is a huge literature about this model, we mention only the recent work [15] related to ours. Note that the potential function U(x) is Lipschitz smooth but not strongly convex in $\mathbb{R}^d$ . In fact, it is double-well and strongly convex only when $|x|\ge r$ for a sufficiently large r. MALA and the classical expectation-maximization (EM) optimization algorithm were compared in [15], which claimed that MALA sampling can be faster than the EM algorithm. We will examine Assumption 1.3 and verify the conditions (1.6a)–(1.6c), so that we can apply MHMC sampling to the double-well potential case.

Proposition 5.2. Consider a two-component Gaussian mixture model as in (5.3). Then Assumptions 1.3 hold.

The proof of Proposition 5.2 can be found in Appendix F.

Remark 5.1. In fact, this result can be extended to the case where $\pi \sim p\mathcal{N}(\mu,\Sigma) + (1-p)\mathcal{N}(\nu,\Sigma)$ for some symmetric positive definite matrix $\Sigma$ . We can find that the linear terms $\langle \mu,x\rangle$ , $\langle \nu,x\rangle$ in $U_2$ are replaced by $\big\langle \Sigma^{-1}\mu,x\big\rangle$ , $\big\langle \Sigma^{-1}\nu,x\big\rangle$ , respectively, and the bounds of $\|\nabla U_2\|$ and $\|\nabla^2 U_2\|$ can be obtained similarly.

Remark 5.2. By the approximation scheme, the allotment only needs to be exhaustive with respect to the function $V_r(x) = \exp(r\Vert x\Vert)$ for some r such that $f(x)\leq V_r(x)$ for all $x\in\mathbb R^d$ . Consequently, the specific allotment construction depends solely on the test function f. Suppose a constant r is such that $f(x)\leq V_r(x)$ for all x; then we could construct an allotment according to Section 3. In particular, let $r_n = \textrm{e}^n$ and $\varepsilon_n = (n^2\sqrt{d} + nrd)^{-1}$ . This gives $L_n = \{x\in\mathbb R^d\colon \Vert x \Vert \leq n/r\}$ and $\tilde{L}_n = \{x\in\mathbb R^d\colon \Vert x \Vert \leq \sqrt{d} + n/r\}$ . Following the construction of $G_j^n$ in Section 3, we can generate an allotment such that $\Vert x - y\Vert<\varepsilon_n\sqrt{d}$ implies $|V_r(x) - V_r(y)|<n^{-1}$ for each n and all $x,y\in \tilde{L}_n$ , by Lagrange’s theorem. By [17, Proposition A.1], such an allotment forms an exhaustive sequence with respect to $V_r$ .

6. Simulations

In this section we explain how to implement the MHMC algorithm for variance reduction through two simple examples. In fact, we have to tackle the following two issues in practice:

• • The stochastic matrix $Q=\big(Q_{ij}\big)_{(m+1)\times (m+1)}$ cannot be computed analytically.

• • Once the approximate solution $\tilde{f}$ has been computed, the function $P_{h, T} \tilde{f}$ , and thus the control variate $P_{h,T}\tilde{f}-\tilde{f}$ , are not accessible in closed form.

Construct a partition $\{G_0,G_1,\ldots,G_m\}$ of $\mathbb{R}^d$ in such a way that the probability of $\pi(G_0)$ is small. We assume that it is easy to sample d-dimensional normal random points. Let $a_i\in G_i$ for $i>0$ be arbitrary and choose $a_0$ on the boundary of $G_0$ .

Let $Y=\{a_0,a_1,\ldots,a_m\}$ , $\mathbb{G}=\{G_0,G_1,\ldots,G_m\}$ . This induces an allotment $(Y, \mathbb{G})$ in $\mathbb{R}^d$ . Recall that $Q_{ij}=P_{h, T}(a_i, G_j)$ , where the transition kernel $P_{h, T}$ was given by (1.5). As the precise computation of each entry is not feasible, we construct an estimate of $Q_{ij}$ via another i.i.d. Monte Carlo. With this in mind, let $Z_1, Z_2, \ldots, Z_N$ be i.i.d. samples from $\mathcal{N}(0, I_d)$ , where N is sufficiently large. Define

(6.1) \begin{equation} \hat{Q}(a_i, a_j) \;:\!=\; \left\{ \begin{array}{l@{\quad}l} \dfrac{1}{N}\displaystyle\sum_{k=1}^{N}\boldsymbol{1}_{Z_k\in G_j}\alpha(a_i,Z_{k}) & \text{if } j\ne i, \\[5pt] 1-\displaystyle\sum_{k\in\{0,1,\dots,m\}\backslash i}\hat{Q}(a_i, a_k) & \text{if } j=i. \end{array}\right.\end{equation}

The approximating transition matrix $Q=\big(Q_{ij}\big)$ is now well estimated by $\hat{Q}=\big(\hat{Q}_{ij}\big) $ by the law of large numbers, so we use $\hat{Q}$ in the MHMC algorithm instead of Q.

Turning to the second issue, given a $\pi$ -integrable function f, we aim to estimate $\mathbb{E}_\pi(f)$ by $(1/k)\sum_{i=0}^{k-1}(f+P_{h,T}\tilde{f}-\tilde{f})(X_i)$ , where $(X_i,\, i\ge 0)$ is a Markov chain generated by the MHMC algorithm. However, $P_{h, T}\tilde{f}$ is not accessible in a closed form once again. Similarly to (6.1), define, for any $x \in \mathbb{R}^d$ ,

(6.2) \begin{equation} \hat{Q}(x, a_j) \;:\!=\; \left\{ \begin{array}{l@{\quad}l} \dfrac{1}{N}\displaystyle\sum_{k=1}^{N}\boldsymbol{1}_{Z_k\in G_j}\alpha(a_i,Z_{k}) & \text{if } j\ne i(x), \\[5pt] 1-\displaystyle\sum_{k\in\{0,1,\dots,m\}\backslash i}\hat{Q}(x,a_k) & \text{if } j=i(x), \end{array}\right.\end{equation}

where i(x) is the unique index $i\in \{0,1,\ldots, m\}$ such that $x\in G_{i(x)}$ ; then define, for any $x \in \mathbb{R}^d$ , $\hat{Q}\tilde{f}(x)=\sum_{j=0}^m \hat{f}(a_j)\hat{Q}(x, a_j)$ . Finally, we use $\hat{Q}\tilde{f}$ in place of $P_{h,T}\tilde{f}$ . We remark that the term $(1/k)\sum_{i=0}^{k-1}(f+\hat{Q}\tilde{f}-\tilde{f})(X_i)$ is unbiased to some extent; the interested reader is referred to [17, Remark 5.1] for details.

In the next subsection we provide a simple example to illustrate the efficiency of variance reduction in MHMC through the weak approximation scheme.

It is not obvious [how] to do so in a way which substantially decreases the variance of the simulation without substantially increasing the complexity of the simulations, and for which optimal values of the parameters can be approximated with a reasonable numerical cost. [12, p. 60].

6.1. Gaussian mixture distribution

Consider a one-dimensional Gaussian mixture distribution as mentioned in Section 5.2. Let $\mu=3$ , $\nu=-2$ , and $p=0.3$ . Choose $f(x) =x^2$ as the force function. For any $m\ge 1$ , $l>0$ , let $G_0^{m,l} =\mathbb{R}\backslash [\!-\!l,l]$ and $G_i^{m,l}$ , $i = 1,2,\dots,m$ , be intervals of length ${2l}/{m}$ partitioning $[{-}l,l]$ . The step size and length are chosen to be $h=0.1$ and $T=20$ . The construction of $\hat{Q}$ is derived from (6.1) using $N = 10^7$ , and the estimator $(1/k)\sum_{i=0}^{k-1}(f+\hat{Q}\tilde{f}-\tilde{f})(X_i)$ is then derived from (6.2) using $N=10^5$ and the approximating scheme. We study three allotments, where the parameters take values as shown in Table 1. We evaluate the associated quantities after 6000 iterations. The results are shown in Figure 1 for different iterations and different estimations from 50 replications. The red curve corresponds to conventional MHMC estimates, and the blue curve corresponds to the modified estimators $(1/k)\sum_{i=0}^{k-1}(f+\hat{Q}\tilde{f}-\tilde{f})(X_i)$ .

Table 1. The parameter setup for different allotments.

We observe that the mean square errors (MSE) of the modified estimator are consistently lower than those of the original MHMC estimator. Moreover, as the allotments become denser, the MSE decreases further. This result suggests that the variance reduction technique is effective and demonstrates that achieving improved convergence in sampling does not require excessively fine allotments. These findings are encouraging and point to the potential for applying this approach to a wider range of models, particularly those in high-dimensional settings.

Figure 1. The estimates for $\mathbb{E}_{\pi}(f)$ using different methods with different allotments.

7. Discussion

There are two interesting recent works in the literature: [17], in which the authors develop an approximation scheme for a solution of the Poisson equations of a geometrically ergodic Metropolis–Hastings chain, and construct a sequence of control-variate estimators to decrease the error variance in the mean estimate; and [9], in which the authors discussed the irreducibility and geometric ergodicity of the popular Hamiltonian Monte Carlo algorithm in a rigorous mathematical sense. Motivated by these two works, we have established the CLT and variance reduction theorem under mild regular conditions for the MHMC algorithm, which implies that the MHMC algorithm can be applied in a wide range of fields.

In addition, we offer several concrete examples satisfying the regular conditions, which suggests that our results can be applied to a broad class of potentials. The simulation results imply that the method works well in practice.

Some questions arise from our work. The simulations suggest that the approximation scheme takes effect even if m is not so large, though it is desirable to obtain non-asymptotic results for the convergence of variance reduction, like in [17, Section 4]. It seems difficult to analyze the convergence rate when using the diffeomorphism mentioned in (2.4), and alternative ways are needed. On the other hand, when implementing the approximation scheme, we use standard Monte Carlo (6.1) to approximate the kernel Q, which may affect the performance of the scheme. Lastly, an interesting question is how to execute an efficient implementation for the approximation scheme. We leave these for future work.