2. Wasserstein Convergence

We start our discussion with the Wasserstein convergence. Following [Reference Rudolf and Schweizer38], we first introduce the metric we use and the notion of ergodicity. For a lower semi-continuous function $V\;:\;\mathbb{R}^d \rightarrow [1,\infty]$ , we define

(6) \begin{equation}d_V(x,y)=(V(x)+V(y))1_{x\neq y}.\end{equation}

For two probability measures $\mu$ and $\nu$ on $\mathbb{R}^d$ , we define

\begin{align*}\|\mu-\nu\|_V=\sup_{|\,f|\leq V}\left |\int f(x)(\mu(dx)-\nu(dx))\right |.\end{align*}

It can be shown that $\|\mu-\nu\|_V=W_{d_V}(\mu,\nu)$ where $W$ denotes the Wasserstein distance (Lemma 3.1 in [Reference Rudolf and Schweizer38]). If we use the constant function $V(x)=1$ , this gives the well-known total variation distance, i.e.

\begin{align*}\|\mu-\nu\|_{TV} = \sup_{|\,f|\leq 1}\left |\int f(x)(\mu(dx)-\nu(dx))\right |.\end{align*}

Note that using $V(x)=1$ neglects the location information of $x$ . This location information can be crucial for problems with unbounded domains.

In general, for problems with unbounded domains, one often chooses $V$ to be a Lyapunov function. Given a Markov chain $X_n$ with transition kernel $P$ , which we will write together as $(X_n, P)$ for short, we say $V\;:\;\mathbb{R}^d\to [1,\infty)$ is a Lyapunov function if there exist $\lambda\in(0,1)$ and $L>0$ such that

(7) \begin{equation}PV (x)=\int P(x,dy) V(y) \leq \lambda V(x) + L,\end{equation}

and the sublevel sets of $V$ are compact. The choice of $V$ depends on the Markov chain and the target density. It can often be chosen as $\|x\|^2$ or $-\log \pi(x)$ or as functions of certain moments; see, for example, [Reference Hairer and Mattingly20, Reference Hairer, Stuart and Vollmer21, Reference Rudolf and Schweizer38].

If $\pi$ is the invariant measure of $X_n$ , we say $X_n$ is geometrically ergodic under $\|\,\cdot\,\|_V$ (see Theorem 16.1 in [Reference Meyn and Tweedie31]) if there are constants $\rho\in(0,1)$ and $C_0\in(0,\infty)$ , such that for any $n\in \mathbb{Z}^+$ ,

(8) \begin{equation} \|\delta_x P^n - \pi\|_V \leq C_0 \rho^n V(x).\end{equation}

We refer to $\rho$ as the ergodicity coefficient. Note that the smaller the value of $\rho$ , the faster the convergence to stationarity. From (8), using the triangle inequality, we obtain an equivalent definition of geometric ergodicity, which requires that for any $x$ and $y$ ,

(9) \begin{equation} \|\delta_x P^n - \delta_yP^n\|_V \leq C_0^{\prime}\rho^n d_V(x,y).\end{equation}

The equivalence can be seen from

(10) \begin{equation}\|\delta_x P^n - \pi\|_V\leq \int \pi(dy)\|\delta_x P^n - \delta_yP^n\|_V \leq C_0^{\prime}\rho^n(V(x)+ \pi V)\leq C_0\rho^nV(x).\end{equation}

where $C_0=C_0^{\prime}(1+ \pi V )$ .

The approximation problem under Wasserstein distance has been studied in [Reference Rudolf and Schweizer38]. We present one of their main results here, which is related to our subsequent development. Interested readers can find more general discussion in the original work.

Theorem 1. (Corollary 3.3 in [Reference Rudolf and Schweizer38]) Suppose $(X_n, P)$ is geometrically ergodic, i.e., as in (9). Suppose $\widehat{V}$ is a Lyapunov function for $(\widehat{X}_n,\widehat{P})$ in the sense of (7) and that

(11) \begin{equation}\|\delta_x P - \delta_x\widehat{P} \|_V \leq \epsilon \widehat{V}(x).\end{equation}

Then, for some constant $C$ , we have

(12) \begin{equation}\|\delta_x P^n - \delta_x\widehat{P}^n\|_V \leq C\epsilon \frac{1-\rho^n}{1-\rho}\left(\widehat{V}(x) + \frac{L}{1-\lambda}\right).\end{equation}

The bound in (12) and the triangular inequality give us an approximation error bound

\begin{align*}\|\delta_x \widehat{P}^n - \pi\|_V\leq \|\delta_x P^n - \delta_x\widehat{P}^n\|_V+\|\delta_x P^n-\pi\|_V\leq C'(\epsilon+\rho^n) (\widehat{V}(x)+V(x))\end{align*}

for some constant $C'$ .

Note that the right-hand side of (12) is not converging to zero as $n\to 0$ . Thus, it cannot help us learn the ergodicity of $\widehat{X}_n$ or whether $\widehat{X}_n$ has a unique invariant measure. The next result shows that ergodicity can be obtained with essentially the same conditions as Theorem 1 (note that condition (8) leads to (13) through (10)).

Theorem 2. Suppose $V$ is a Lyapunov function for $P$ in the sense of (7). In addition, assume there exist $N\in \mathbb{Z}^+$ and $\rho\in(0,1)$ , such that for any $n\geq N$ ,

(13) \begin{equation}\|\delta_xP^n-\delta_y P^n\|_V\leq \rho^nd_V(x,y).\end{equation}

Lastly, there is an $\epsilon_0>0$ , and the following holds for some $\epsilon\in(0,\epsilon_0]$ ,

(14) \begin{equation}\|\delta_x P - \delta_x \widehat{P}\|_V\leq \epsilon V (x).\end{equation}

Then, $V$ is a Lyapunov function for $\widehat{P}$ as well with

\begin{align*}\widehat{P} V(x)\leq (\lambda+\epsilon) V(x) + L.\end{align*}

Moreover, $\widehat{X}_n$ has a unique invariant measure $\widehat{\pi}$ and there exist $C_1,D_1\in(0,\infty)$ independent of $\epsilon$ , such that

\begin{align*}\|\delta_x\widehat{P}^n-\delta_y \widehat{P}^n\|_V\leq C_1(\rho+D_1\epsilon)^{n}d_V(x,y).\end{align*}

In Theorem 2, $\epsilon_0$ is chosen such that $\lambda+\epsilon_0<1$ and $\rho+D_1\epsilon_0<1$ . Theorem 2 indicates that if $P$ is geometrically ergodic with ergodicity coefficient $\rho$ and if $\widehat{P}$ is $\epsilon$ -close to $P$ as characterized by (14), $\widehat{P}$ is also geometrically ergodic. Moreover, the ergodicity coefficient of $\widehat{P}$ is bounded above by $\rho+D_1\epsilon$ .

In statistical applications, we are more interested in turning convergence results into error bounds for the Monte Carlo estimators. The central limit theorem of ergodic Markov chains was studied in [Reference Jones26, Reference Tierney44], which provide asymptotic error quantifications. In practice, non-asymptotic bounds for finite values of $n$ may be more desirable. The following proposition appeared in [Reference Joulin and Ollivier27, Reference Rudolf37]. We provide an explicit statement here to show the variance bound along with a simple proof for self-completeness. For simplicity, we assume the Markov chain is initialized with the invariant measure, i.e. $\widehat{X}_0\sim\widehat{\pi}$ , so a burn-in period is not necessary.

Proposition 1. Suppose $\|\delta_x\widehat{P}^n-\delta_y \widehat{P}^n\|_V\leq \widehat{\rho}^n d_V(x,y)$ for some $\widehat{\rho}\in (0,1)$ . Then, for any $f$ that is 1-Lipschitz under $\|\,\cdot\,\|_V$ ,

\begin{align*}\left|\delta_x \widehat{P}^n f -\widehat{\pi} f\right|\leq \widehat{\rho}^n(V(x)+\widehat{\pi} V).\end{align*}

In addition, if we use $\hat{f}_M=\frac{1}{M}\sum_{k=1}^M f(\widehat{X}_k)$ as an estimator of $\widehat{\pi} f$ starting from $\widehat{X}_0\sim \widehat{\pi}$ ,

\begin{align*}\mathbb{E}_{\widehat{\pi}}\left[(\hat{f}_M-\widehat{\pi} f)^2\right]\leq \frac{2}{ (1-\widehat{\rho})M}\mathbb{E}_{\widehat{\pi}}\Big[|\,f(\widehat{X}_0)|(V(\widehat{X}_0)+\widehat{\pi} V)\Big].\end{align*}

In many applications, we are interested in the properties of $\widehat{\pi}$ on compact regions. In these scenarios, the associated test functions will be bounded, and $1$ -Lipschitz will be under $\|\,\cdot\,\|_V$ . To learn tail properties of $\widehat{\pi}$ such as intermittency, the test functions often need to grow with $\|x\|$ . In order to apply Proposition 2, we would need the Lyapunov function $V$ to grow at a similar scale.

3. $\chi^2$ convergence

In this section, we discuss convergence under the $\chi^2$ divergence. We first introduce some notations. For a transition kernel $P$ and density $\mu$ , we define

\begin{align*}\|P\|_{\mu} = \max_{f:0<\|\,f\|_{\mu}<\infty} \frac{\|Pf\|_{\mu}}{\|\,f\|_{\mu}},\end{align*}

where $\|\,f\|_{\mu}^2=\langle f, f\rangle_{\mu}$ . For two probability measures $\mu$ and $\nu$ on $\mathbb{R}^d$ , where $\nu$ is absolutely continuous with respect to $\mu$ , we define the $\chi^2$ divergence of $\nu$ from $\mu$ as:

\begin{align*}D_{\chi^2}(\nu\|\mu)=\int \left(\frac{\nu(x)}{\mu(x)}-1\right)^2\mu(x)dx=\int\frac{\nu(x)^2}{\mu(x)}dx-1.\end{align*}

For a transition kernel $P$ that is reversible with respect to $\pi$ , the spectral gap of $P$ is defined as [Reference Hairer, Stuart and Vollmer21]

(15) \begin{equation} \kappa(P)=1-\sup\left\{\frac{\|Pf-\pi f\|_{\pi}^2}{\|\,f-\pi f\|_{\pi}^2}:\ \; f\in L^2(\pi), \textrm{var}_\pi f \neq 0\right\}.\end{equation}

Note that by repeatedly applying (15), we have for any $f\in L^2(\pi)$

\begin{align*}\|P^n f -\pi f\|_{\pi}^2 \leq (1-\kappa(P))^n \|\,f-\pi f\|_{\pi}^2.\end{align*}

Thus, the larger the spectral gap, the more quickly $X_n$ converges to its invariant measure.

Remark 2. An alternative definition of the spectral gap takes the form [Reference Andrieu and Vihola2, Reference Bakry, Gentil and Ledoux3]

\begin{align*}\kappa_a(P)=\inf\left\{\frac{\langle f, (I-P)f\rangle_{\pi}}{\textrm{var}_{\pi} f}\;:\; f\in L^2(\pi), \textrm{var}_\pi f \neq 0\right\}.\end{align*}

Note that these two spectral gaps are related through $1-\kappa(P)=(1-\kappa_a(P))^2$ .

3.1. General $\chi^2$ approximation and convergence

Theorem 3. Our first result assumes that $(X_n,P)$ has a spectral gap and that $\widehat{P}$ is a close approximation of $P$ : Suppose $P$ is a reversible transition kernel with invariant measure $\pi$ and a spectral gap $\kappa(P)>0$ in the sense of (15). There is an $\epsilon_0>0$ such that for any transition kernel $\widehat{P}$ satisfying $\|P-\widehat{P}\|_\pi\leq \epsilon<\epsilon_0$ and any $a\in(0,1)$ , there exists a constant $C$ that may depend on $\kappa(P)$ such that the following holds with $\widehat{\kappa}=(1-a)\kappa(P)-\frac{C\epsilon^2}{a}\in(0,1)$ :

1. (i) For any $f\in L^2(\pi)$ ,

   \begin{align*}\|\widehat{P}^nf-\pi\widehat{P}^nf\|_{\pi}^2 \leq (1-\widehat{\kappa})^{n}\textrm{var}_\pi f \end{align*}

2. (ii) $\widehat{P}$ has an invariant measure $\widehat{\pi}$ , which satisfies

   \begin{align*}|(\widehat{\pi}-\pi \widehat{P}^n)f|^2 \leq C\epsilon^2(1-\widehat{\kappa})^{n}\textrm{var}_\pi f.\end{align*}

   Moreover, $D_{\chi^2}(\widehat{\pi}\|\pi)\leq C\epsilon^2$ .

In Theorem 3, $\epsilon_0$ is chosen such that $\hat\kappa\in(0,1)$ . Theorem 3 indicates that if $\widehat{P}$ and $P$ are $\epsilon$ -close to each other as quantified by $\|\widehat{P}-P\|_{\pi}\leq \epsilon$ , then $\widehat{P}$ has a stationary distribution $\widehat{\pi}$ . Moreover, $\widehat{\pi}$ and $\pi$ are $\epsilon$ -close to each other as quantified by $D_{\chi^2}(\widehat{\pi}\|\pi)\leq C\epsilon^2$ . We also note that showing that $\|\widehat{P}^nf-\pi\widehat{P}^nf\|_{\pi}^2 \leq (1-\widehat{\kappa})^{n} \textrm{var}_\pi f$ is different than finding the spectral gap of $\widehat{P}$ , since the latter would need a similar inequality but with $\pi$ replaced by $\widehat{\pi}$ . In other words, Theorem 3 does not provide a spectral gap for $\widehat{P}$ . On the other hand, we can obtain error bounds for Monte Carlo estimators using the bounds established in Theorem 3:

Proposition 2. Under the same conditions as those in Theorem 3, for any $f\in L^2(\pi)$ and any initial distribution $\widehat{X}_0\sim \nu\ll\pi$ , there exists a constant $C$ such that

\begin{align*}\left|\nu \widehat{P}^n f -\widehat{\pi} f\right|^2\leq (1-\widehat{\kappa})^{n}\textrm{var}_\pi(f)\left(\sqrt{D_{\chi^2}(\nu\|\pi)+1}+C\epsilon\right)^2.\end{align*}

In addition, if $f$ is bounded, there exists a constant $C$ such that

\begin{align*}\mathbb{E}_{\widehat{\pi}}[(\hat{f}_M-\widehat{\pi} f)^2]\leq \frac{C}{M(1-(1-\widehat{\kappa})^{1/4})}\sqrt{\textrm{var}_{\widehat{\pi}}(f)\textrm{var}_\pi (f)},\end{align*}

where $\hat{f}_M=\frac{1}{M}\sum_{k=1}^M f(\widehat{X}_k)$ .

Remark 3. [Reference Negrea and Rosenthal32] provides a result similar to the first claim in Theorem 3 (see Lemma A.6 in [Reference Negrea and Rosenthal32]). But it requires stronger assumptions on $\widehat{P}$ , namely, it requires that $\widehat{P}$ is ergodic and aperiodic and that it is a mapping from $L_2(\pi)$ to $L_2(\pi)$ . Our result does not require these assumptions.

3.2. Spectral gap with density ratio bounds

In this section, we show that stronger results can be established if we can bound the ratio between the invariant densities $\pi$ and $\widehat{\pi}$ . Such a bound is assessable if we have an explicit characterization of $\widehat{\pi}$ . For example, in Bayesian inverse problems, $\pi(x)\propto p_0(x)\exp\!({-}\frac12\|G(x)-y\|^2)$ , while $\widehat{\pi}(x)\propto p_0(x)\exp\!({-}\frac12\|\widehat{G}(x)-y\|^2)$ . In this case, a density ratio bound can be obtained if $\|G(x)-\widehat{G}(x)\|$ is bounded. This is practically feasible by using an accurate numerical approximation of $G$ ; and the approximation error can be estimated by the grid size or Galerkin truncation used in the numerical scheme (see, for example, [Reference Hoang, Schwab and Stuart23]). Similar assumptions have also been imposed in existing Bayesian computation literature; see, for example, [Reference Beskos, Jasra, Law, Marzouk and Zhou5].

Theorem 4. Suppose $P$ and $\widehat{P}$ are two reversible transition kernels with invariant densities $\pi$ and $\widehat{\pi}$ respectively. We further assume that $\pi(x)/\widehat{\pi}(x)\in [(1+\epsilon)^{-1}, 1+\epsilon]$ , and that $\|P-\widehat{P}\|_{\pi}\leq\epsilon$ . Then there exists a universal constant $C$ such that

\begin{align*}\kappa(\widehat{P})\geq \kappa(P)-C\epsilon.\end{align*}

Note that for $\epsilon$ small enough, $\kappa(P)-C\epsilon>0$ . Based on the spectral gap, we have the following non-asymptotic Monte Carlo error bound:

Proposition 3. Suppose $(\widehat{X}_n,\widehat{P})$ has a spectral gap $\widehat{\kappa}$ . Suppose the initial distribution is $\nu$ , i.e. $\widehat{X}_0\sim \nu$ . Then

\begin{align*}\left|\mathbb{E} f(\widehat{X}_n) -\widehat{\pi} f\right|^2\leq (1-\widehat{\kappa})^{n}\textrm{var}_{\widehat{\pi}}f(D_{\chi^2}(\nu\|\widehat{\pi})+1).\end{align*}

In addition, if $\nu=\widehat{\pi}$ , then for $\hat{f}_M=\frac{1}{M}\sum_{k=1}^M$ , we have

\begin{align*}\mathbb{E}_{\widehat{\pi}}[(\hat{f}_M-\widehat{\pi} f)^2]\leq { \frac{2}{M(1-(1-\widehat{\kappa})^{1/2})}}\textrm{var}_{\widehat{\pi}}[f].\end{align*}

Note that Proposition 3 is not a new result. A more delicate central limit theorem version of it can be found as Theorem 4.4 of [Reference Hairer, Stuart and Vollmer21]. We provide a short proof of the proposition in the Appendix B for self-completeness.

Before we conclude our discussion of the $\chi^2$ convergence, we remark that even though the condition $\|P-\widehat{P}\|_\pi\leq \epsilon$ is reasonable for the spectral gap analysis, it can be hard to verify directly in some applications. To remedy this issue, the next proposition shows that we can bound $\|P-\widehat{P}\|_\pi$ through a bound for $\|\delta_x P - \delta_x\widehat{P}\|_{TV}$ , which can be easier to obtain using coupling techniques.

Proposition 4. Suppose there exists a $\pi$ -measurable function $V\;:\;\mathbb{R}^d\rightarrow[1,\infty)$ such that $\|\delta_x P - \delta_x\widehat{P}\|_{TV}\leq \epsilon V(x)$ . In addition, suppose $\frac1a\leq \pi(x)/\hat\pi(x)\leq a$ for some constant $a>0$ . Then

\begin{align*}\|P-\widehat{P}\|_{\pi}\leq \sqrt{2(1+a^2)}\sqrt{\epsilon}\|V\|_{\pi}^{1/2}.\end{align*}

4. Application: Metropolis–Hastings MCMC on perturbed densities

Random walk Metropolis (RWM) and Metropolis adjusted Langevin algorithm (MALA) are two popular MCMC samplers when it comes to sampling a generic density $\pi$ . Many existing works have already studied their spectral gap under suitable conditions on $\pi$ [Reference Dwivedi, Chen, Wainwright and Yu15, Reference Hairer, Stuart and Vollmer21, Reference Roberts and Tweedie36]. When implementing these samplers, it is often the case that we have access only to an approximation of $\pi$ , which we denote as $\widehat{\pi}$ . In this section, we will demonstrate how to apply our analysis framework to establish proper bounds for the spectral gap of the “numerical” RWM and MALA.

In fact, we can develop some general results for the Metropolis–Hastings (MH) type of Monte Carlo algorithm. Assume the proposals are given by some smooth transition density $R(x,x^{\prime})$ . Due to the possibility of rejection, MH Monte Carlo transition densities w.r.t. Lebesgue measure can be written as $P(x,x^{\prime})=\gamma(x)\delta_x(x^{\prime})+\beta(x,x^{\prime})$ with

(16) \begin{equation} \beta(x,x^{\prime})=\min\left\{\frac{\pi(x^{\prime})R(x^{\prime},x)}{\pi(x)}, R(x,x^{\prime})\right\},\quad \gamma(x)=1-\int \beta(x,x^{\prime})dx^{\prime}.\end{equation}

The perturbed transition density can be written as $\widehat{P}(x,x^{\prime})=\hat{\gamma}(x)\delta_x(x^{\prime})+\widehat{\beta}(x,x^{\prime})$ . We provide some sufficient conditions under which the difference between $P$ and $\widehat{P}$ is of order  $\epsilon$ .

Lemma 1. If the transition density is of the form $P(x,x^{\prime})=\gamma(x)\delta_x(x^{\prime})+\beta(x,x^{\prime}) $ with $\nu(x) P(x,x^{\prime})=\nu(x^{\prime}) P(x^{\prime},x)$ , suppose that $\widehat{P}(x,x^{\prime})=\hat{\gamma}(x)\delta_x(x^{\prime})+\hat{\beta}(x,x^{\prime})$ with

\begin{align*}|\hat{\gamma}(x)-\gamma(x)|\leq C\epsilon\ \ and \ \ (1-C\epsilon )\beta(x,x^{\prime})\leq \hat{\beta}(x,x^{\prime})\leq (1+C\epsilon) \beta(x,x^{\prime}).\end{align*}

for some constant $C\in(0,\infty)$ . Then there exists a constant $C_1\in(0,\infty)$ such that $\|P-\widehat{P}\|_\nu\leq C_1\epsilon.$

4.1. Random walk Metropolis

RWM considers implementing the MH procedure on random walk proposals. That is, we use

\begin{align*}R(x,x^{\prime})=\frac{1}{(2\pi h )^{d/2}}\exp\left({-}\frac1{4h}\|x^{\prime}-x\|^2\right)\end{align*}

in (16). It is worth noting that using a perturbed density $\widehat{\pi}$ does not affect this proposal.

Proposition 5. For RWM, there is an $\epsilon_0 > 0$ so that for any $\epsilon<\epsilon_0$ and $\sup_x|\log\pi(x)-\log\hat\pi (x)|\leq C\epsilon$ , there is a constant $C_1$ so that

\begin{align*}\|P_{RWM} - \widehat{P}_{RWM}\|_{\pi}\leq C_1\epsilon.\end{align*}

If the original RWM has a spectral gap and if $\sup_x|\log\pi(x)-\log\hat\pi (x)|\leq C\epsilon$ , then Proposition 5, together with Theorem 4, implies that the perturbed RWM has a proper spectral gap as well. In practice, an estimate of $\epsilon$ can be obtained by analyzing the numerical scheme used; see Section 3.2 for more details.

4.2. Metropolis adjusted Langevin algorithm

MALA considers implementing the MH procedure on proposals following the Langevin diffusion. That is, we use

\begin{align*}R(x,x^{\prime})=\frac{1}{(4\pi h )^{d/2}}\exp\left({-}\frac1{4h}\|x^{\prime}-x-h\nabla \log \pi(x)\|^2\right)\end{align*}

in (16). Using a perturbed density $\widehat{\pi}$ does change this proposal. We discuss the perturbation in two separate cases. In particular, we shall verify that the condition $\|P_{MALA}-\widehat{P}_{MALA}\|_{\pi}\leq \epsilon$ holds under appropriate assumptions on $\widehat{\pi}$ in the two cases. Then, if $P_{MALA}$ has a spectral gap, the numerical sampler $\widehat{P}_{MALA}$ has a proper spectral gap as well.

4.2.1. Bounded domain

When the support of $\pi$ and $\widehat{\pi}$ is bounded, the analysis is quite straightforward with Lemma 1.

Proposition 6. For MALA, there is an $\epsilon_0>0$ so that for any $\epsilon<\epsilon_0$ , if $\sup_x|\log\pi(x)-\log\hat\pi(x)|\leq C\epsilon$ , $\sup_x\|\nabla\log\pi(x)-\nabla\log\hat\pi(x)\|\leq C\epsilon$ and if the support of $\pi$ and $\widehat{\pi}$ is bounded, then

\begin{align*}\|P_{MALA}-\widehat{P}_{MALA}\|_{\pi}=O(\epsilon).\end{align*}

4.2.2. Unbounded support

When the support of the density is unbounded, directly bounding $\|P_{MALA}-\widehat{P}_{MALA}\|_{\pi}$ becomes difficult. Instead, we consider establishing $\| \delta_xP- \delta_x\widehat{P}\|_{TV}=O(\epsilon)$ .

Proposition 7. For MALA, if $\log \pi$ is Lipschitz, $\sup_x|\log\pi(x)-\log\hat\pi(x)|\leq L_\pi\epsilon$ and, moreover, $\sup_x\|\nabla\log\pi(x)-\nabla\log\widehat{\pi}(x)\|\leq L_\pi\epsilon$ , for any $\delta>0$ , there exists $C_{\delta}\in(0,\infty)$ such that for $h<(\frac{5L_{\pi}}{\delta} + 20L_{\pi})^{-1}$ ,

\begin{align*}\| \delta_xP- \delta_x\widehat{P}\|_{TV}\leq C_{\delta}\epsilon\exp\!(\delta \|x\|^2).\end{align*}

When $\pi(x)$ is sub-Gaussian, we can find a $\delta>0$ such that $V(x)=\exp\!(\delta \|x\|^2)$ is $L_2$ -integrable under $\pi$ . Then Proposition 4 indicates that $\|P- \widehat{P}\|_{\pi}=O(\sqrt{\epsilon})$ .

5. Application: Parallel Tempering with Perturbed Densities

In this section, we demonstrate how to apply our framework to parallel tempering (PT) algorithms [Reference Earl and Deem17, Reference Tawn and Roberts42, Reference Tawn, Roberts and Rosenthal43]. These algorithms are also referred to as the replica exchange methods [Reference Dong and Tong11, Reference Dupuis, Liu, Plattner and Doll13, Reference Sugita and Okamoto41]. Compared with regular MCMC samplers like RWM and MALA, PT tries to sample a multiple-tempered version of the target density. Such a design can improve the convergence rate on densities with multiple isolated modes.

To implement PT, a sequence of distributions $\pi_0,\ldots, \pi_K$ is considered, where the last one is the target density $\pi_K =\pi$ . The first density, $\pi_0$ , is usually a distribution that is easy to draw samples from. The intermediate distributions, $\pi_k$ ’s $1\leq k\leq K-1$ , are set up so that the two neighboring densities are similar to each other. A common choice for the intermediate distributions is to consider interpolations between $\pi_K$ and $\pi_0$ :

\begin{align*}\pi_k(x)\propto\pi^{\beta_k}(x)\pi_0^{1-\beta_k}(x),\end{align*}

where $0=\beta_0<\beta_1<\ldots <\beta_K=1$ is a sequence of parameters. PT intends to generate samples from the product density

\begin{align*}\Pi=\pi_0\times \pi_1\times\cdots\times \pi_K \mbox{ on $\mathbb{R}^{d(K+1)}$.}\end{align*}

To do so, its iterations consist of $K+1$ parts, i.e. $X_n=(X^0_{n},\ldots,X^K_{n})$ , and the updating rule is given by the following two steps:

1. (i) Update each $X^k_{n}$ to $X^{k}_{n+1}$ according to a transition kernel $M_k$ , whose stationary distribution is $\pi_k$ . In practice, $M_k$ is often taken as the transition kernel obtained by repeating the RWM or MALA update for $t_k$ steps. That is, $M_k=P_{RWM}^{t_k}$ or $M_k=P_{MALA}^{t_k}$ .

2. (ii) Pick an index $k\in \{0,\ldots,K-1\}$ uniformly at random, and swap the values of $X^k_{n+1}$ and $X^{k+1}_{n+1}$ with probability $\alpha_k(X_{n+1}^k,X^{k+1}_{n+1})$ , where

\begin{align*}\alpha_k(x,x^{\prime})=\min\left\{1, \frac{\pi_k(x^{\prime})\pi_{k+1}(x)}{\pi_{k}(x)\pi_{k+1}(x^{\prime})}\right\}.\end{align*}

The pseudo-code of PT is given in Algorithm 1.

The exchange procedure can be described by the transition probability on $R^{(K+1)d}\times R^{(K+1)d}$ :

\begin{align*}Q_k(x,x)=1-\alpha_k(x^k, x^{k+1}),\quad Q_k(x,S_k(x))=\alpha_k(x^k, x^{k+1}),\end{align*}

where $S_k(x)=(x^0,\ldots, x^{k-1},x^{k+1},x^k, x^{k+2},\ldots, x^K)$ . With a little change of notation, we write the transition kernel as $Q_k$ as well; i.e. $Q_kf(x)=Q_k(x,x)f(x)+Q_k(x,S_k(x))f(S_k(x))$ . The transition kernel of PT can then be written as

(17) \begin{equation}P=\left(M_{0}\otimes \cdots \otimes M_K\right) \left(\frac{1}{K}\sum_{\mathbf{k}\in \{0,\ldots,K-1\}} Q_{k}\right),\end{equation}

where the direct product of two transition kernels is given by

\begin{align*}M_0\otimes M_1f(x^0,x^1)=\int\int M_0(x^0,y^0)M_1(x^1,y^1) f(y^0,y^1) dy^0dy^1.\end{align*}

Algorithm 1

Parallel Tempering

The spectral gap of $P$ in (17) has been studied in [Reference Woodard, Schmidler and Huber46]. Assuming the state space can be partitioned into $\mathbb{R}^d=\cup_{i=1}^J A_j$ , it is shown that $\kappa(P)$ can be seen as the product of three elements: (1) the maximal spectral gap when sampling $\pi_k$ , $k\geq 1$ , constrained on one piece, $A_j$ ; (2) the spectral gap when sampling $\pi_0$ using $M_0$ ; and (3) the density ratio $\pi_k(A_j)/\pi_{k+1}(A_j)$ . In particular, if $\pi_0$ is easy to sample, $\pi_k$ is not very different from $\pi_{k+1}$ , and the sampling of $\pi_k$ constrained on $A_j$ is efficient, then PT can be highly efficient. When implementing PT numerically, we may not have access to the exact values of $\pi_k$ but only an $\epsilon$ -approximation, which we denote as $\widehat{\pi}_k$ . Then the corresponding PT uses a sampler $\widehat{M}_k$ , with invariant measure $\widehat{\pi}_k$ at each replica, while the exchange probability is given by

\begin{align*}\widehat{\alpha}_k(x,x^{\prime})=\min\left\{1, \frac{\widehat{\pi}_k(x^{\prime})\widehat{\pi}_{k+1}(x)}{\widehat{\pi}_{k}(x)\widehat{\pi}_{k+1}(x^{\prime})}\right\}.\end{align*}

The corresponding transition kernel can be written as

\begin{align*}\widehat{P}=\left(\widehat{M}_{0}\otimes \cdots \otimes \widehat{M}_K\right) \left(\frac{1}{K}\sum_{\mathbf{k}\in \{0,\ldots,K-1\}} \widehat{Q}_{k}\right).\end{align*}

It is natural to ask whether this numerical PT will inherit the spectral gap of $P$ . The next result, together with Theorem 4, indicate that under appropriate regularity conditions on $\widehat{\pi}_k$ ’s, the numerical PT also has a proper spectral gap.

Proposition 8. Suppose that for each replica, if the target distribution satisfies $\sup_x|\log\widehat{\pi}_k(x)-\log\pi_k(x)|\leq \epsilon$ and the transition kernel satisfies $\|P_k-\widehat{P}_k\|_{\pi_k}\leq \epsilon$ , then the transition kernel of PT satisfies the following for some constant $C$ :

\begin{align*}\|P-\widehat{P}\|_{\Pi}\leq C\epsilon.\end{align*}

Before we prove Proposition 8, we first prove two auxiliary lemmas. The first lemma shows that different compositions of approximated transition kernels yield approximation kernels of similar accuracy. In particular, it helps us establish the condition $\|P_k-\widehat{P}_k\|_{\pi_k}\leq \epsilon$ in Proposition 8 if we use $M_k=P_{RWM}^{t_k}$ or $M_k=P_{MALA}^{t_k}$ .

The second lemma establishes proper bounds for the swapping transition.

Lemma 3. Let $Q$ be a transition probability of form

\begin{align*}Q(x,S(x))=a(x,S(x)),\quad Q(x,x)=1-a(x,S(x)),\end{align*}

where $S(x)$ is some given map. Suppose $Q$ is reversible with a density $\nu$ , i.e.

\begin{align*}\nu(x) Q(x,S(x))=\nu(S(x)) Q(S(x),x).\end{align*}

Similarly, let $\widehat{Q}$ denote the transition probability of the form

\begin{align*}\widehat{Q}(x,S(x))=\widehat{a}(x,S(x)),\quad \widehat{Q}(x,x)=1-\widehat{a}(x,S(x)),\end{align*}

reversible with $\widehat{\nu}$ . If for some constant $C$ , $(1-C\epsilon)a(x,S(x))\leq \widehat{a}(x,S(x))\leq (1+C\epsilon)a(x,S(x))$ , then

\begin{align*}\|(Q-\widehat{Q}) f\|_{\nu}\leq 2C \epsilon \|\,f\|_\nu.\end{align*}

6. Numerical examples

In this section, we present some numerical examples based on the predator–prey system to illustrate the theoretical results developed in the preceding sections.

6.1. Predator–prey system

We consider inferring the parameters of a system of ordinary differential equations (ODEs) that models the predator–prey system [Reference Lotka29]. Denoting the populations of prey and predator by $( \gamma_p, \gamma_q )$ , the populations change over time according to the pair of coupled ODEs is

(18) \begin{align}\frac{d \gamma_p}{dt} & = r \gamma_p \left(1 - \frac{\gamma_p}{K}\right) - s \left(\frac{\gamma_p \,\gamma_q }{w + \gamma_p}\right), \nonumber \\\frac{d \gamma_q}{dt} & = u \left(\frac{\gamma_p \,\gamma_q }{w + \gamma_p}\right) - v \gamma_q, \end{align}

with initial conditions $\gamma_p(0)$ and $\gamma_q(0)$ . $r$ , $K$ , $a$ , $s$ , $u$ , and $v$ are model parameters that control the dynamics of the populations of prey and predator. In the absence of the predator, the population of prey evolves according to the logistic equation, which is characterized by $r$ and $K$ . In the absence of the prey, the population of the predator has an exponential decay rate $v$ . The additional parameters $s,w$ , and $u$ characterize the interaction between the predator population and the prey population.

In the inference problem, we want to estimate both the model parameters and the initial conditions. In this case, we have $d = 8$ and denote

\begin{align*}\theta = ( \gamma_p(0), \gamma_q(0), r, K, a, s, u, v ).\end{align*}

A commonly used prior for this problem is a uniform distribution over a hypercube $(a_1,b_1) \times \cdots \times (a_d,b_d) $ (see e.g. [Reference Parno and Marzouk33]). Here we set $a_i = 10^{-3}$ and $b_i = 2\times 10^2$ for all $i$ . Noisy observations of both $\gamma_p(t;\; \theta)$ and $\gamma_q(t;\; \theta)$ at times regularly spaced at $m = 20$ time points in $t \in [2, 40]$ are used to infer $\theta$ . This defines a so-called forward model,

\begin{align*}F(\theta) = [\gamma_p(t_1;\; \theta), \gamma_q(t_1;\; \theta), \ldots, \gamma_p(t_m;\; \theta), \gamma_q(t_m;\; \theta)],\end{align*}

that maps a given parameter $\theta$ to the observables. The observables are perturbed with independent Gaussian observational errors with mean zero and variance $4$ . A “true” parameter

\begin{align*}\theta_{\rm true} = [50, 5, 0.6, 100, 1.2, 25, 0.5, 0.3]^\top\end{align*}

is used to generate the synthetic observed data set, which is denoted by $y$ . The trajectories of $\gamma_p(t;\; \theta_{\rm true})$ and $\gamma_q(t;\; \theta_{\rm true})$ together with the synthetic data set are shown in Figure 1.

Figure 1.

Left and middle: the trajectories of $\gamma_p(t;\; \theta_{\rm true})$ and $\gamma_q(t;\; \theta_{\rm true})$ computed using the second-order Runge–Kutta method, with different time step size $\eta$ . Right: the $L_2$ error of the model outputs with different time step sizes $\eta$ . Here $G(\theta_{\rm true})$ is computed using $\eta = \eta_0 \times 2^{-6}$ . The trajectories computed by the time step size $\eta = \eta_0 \times 2^{-6}$ are used to generate a synthetic data set. The observed data sets of the prey and the predator are shown as circles and squares, respectively.

To prevent rejections caused by proposal samples that fall outside of the hypercube, we further consider the prior distribution as the pushforward of the standard Gaussian measure, with the probability density function

\begin{align*}p_0(x) = (2\pi)^{-d/2} \exp\left( -\frac12 \|x\|^2 \right)\end{align*}

under a diffeomorphic transformation $T\; :\; \mathbb{R}^d \rightarrow \mathbb{R}^d$ that maps each coordinate

\begin{align*}\theta_i = T_i(x_i) = a_i + \frac{b_i - a_i}{\sqrt{2\pi} }\int_{-\infty}^{x_i} \exp\left( -\frac12 z^2 \right) d z.\end{align*}

In other words, $p_0(x)$ is the prior distribution for the transformed parameter $x = T^{-1}(\theta)$ . Writing $G(x) = F(T(x))$ , our goal is to characterize the posterior distribution

\begin{align*}\pi(x) \propto p_0(x) \exp\left( - \frac{1}{8} \|G(x) - y\|^2 \right).\end{align*}

The system of ODEs in (18), and hence the function $G(x)$ , has to be numerically solved by some ODE solvers. Here we use the second-order explicit Runge–Kutta method, with time step size $\eta$ , to solve (18). As shown in Figure 1, the trajectories of $\gamma_p(t;\; \theta_{\rm true})$ and $\gamma_q(t;\; \theta_{\rm true})$ converge as $\eta \rightarrow 0$ at a rate of $O(\eta^2)$ (see [Reference Butcher8] for a detailed analysis). The numerical solver, which is characterized by the step size $\eta$ , defines the approximate model $\widehat{G}(x)$ and the approximate posterior density $\widehat{\pi}(x)$ . Figure 2 shows the estimated marginal distributions (using Algorithm 1) of perturbed posteriors defined by various time step sizes. Here we observe that as $h$ decreases, the estimated marginal distributions almost overlap each other, which suggests that the perturbed distributions converge as the discretized model $\widehat{G}$ converges.

Figure 2.

Marginal distributions of perturbed posteriors defined by various time step sizes.

6.2. MCMC results

To validate the theoretical results of Metropolis–Hasting MCMC on perturbed densities in Section 4, we first simulate the RWM algorithm with invariant densities $\widehat{\pi}(x)$ defined by various time step sizes, as shown in Figure 1. All the Markov chains in this set of simulation experiments are generated using the same Gaussian random walk proposal distribution. The box plots of the integrated autocorrelation times (IACTs) of the resulting Markov chains are shown in Figure 3. Then we simulate MALA with invariant densities $\widehat{\pi}(x)$ defined by the same set of time step sizes. The box plots of the resulting IACTs are shown in Figure 4. Again, all the Markov chains are generated using the same proposal distribution. For both algorithms, we simulate each Markov chain for $10^{6}$ iterations after discarding burn-in samples and repeat the simulation for $20$ times with different initial states to produce the box plots. As established in our theoretical analysis, for both algorithms the resulting Markov chains targeting various approximate posterior densities produce similar IACTs.

Figure 3.

Box plots of integrated autocorrelation times (IACTs) of the first four parameter Markov chains simulated by the RWM algorithm. Here, $20$ realizations of Markov chains are used to estimate the IACT.

Figure 4.

Box plots of integrated autocorrelation times (IACTs) of the first four parameter Markov chains simulated by the MALA algorithm. Here, $20$ realizations of Markov chains are used to estimate the IACT.

To validate the theoretical results on the parallel tempering with perturbed densities in Section 5, we simulate Algorithm 1 with the same Gaussian random walk as in RWM. For each of the invariant densities $\widehat{\pi}(x)$ defined by various time step sizes, we set $K=4$ , and the intermediate distributions take the form

\begin{align*}\widehat{\pi}_k(x) \propto p_0(x) \exp\left( - \frac{\beta_k}{8} \|\widehat{G}(x) - y\|^2 \right),\end{align*}

where $\beta_k = 1 + \alpha^{-K} - \alpha^{-k}$ with $\alpha = 1.3$ and $k = \{0,1,2,3,4\}$ . Here $\beta_k$ is an increasing sequence such that $\beta_K = 1$ . The same Gaussian random walk is used across all replicas to simulate the Markov chain. The box plots of the IACTs of resulting Markov chains are shown in Figure 5. Similar to the previous experiments, the resulting Markov chains targeting various approximate posterior densities produce similar IACTs. Furthermore, we noticed that the IACTs of Algorithm 1 are smaller than those of the RMW and MALA algorithm, which suggests that Algorithm 1 has a better mixing rate.

Figure 5.

Box plots of integrated autocorrelation times (IACTs) of the first four parameters of Markov chains simulated by Algorithm 1. Here, $20$ realizations of Markov chains are used to estimate the IACT.