Perturbation theory for Markov chains via Wasserstein distance

Introduction

Markov chain Monte Carlo (MCMC) algorithms are one of the key tools in computational statistics. They are used for the approximation of expectations with respect to probability measures given by unnormalized densities. For almost all classical MCMC methods it is essential to evaluate the target density. In many cases, this requirement is not an issue, but there are also important applications where it is a problem. This includes applications where the density is not available in closed form, see , or where an exact evaluation is computationally too demanding, see . Problems of this kind lead to the approximation of Markov chains and to the question of how small differences in the transitions of two Markov chains affect the differences between their distributions.

In Bayesian inference when big data sets are involved an exact evaluation of the target density is typically very expensive. For instance, in each step of a Metropolis-Hastings algorithm the likelihood of a proposed state must be computed. Every observation in the underlying data set contributes to the likelihood and must be taken into account in the calculation. This may result in evaluating several terabytes of data in each step of the algorithm. These are the reasons for the recent interest in numerically cheaper approximations of classical MCMC methods, see . A reduction of the computational costs can, e.g., be achieved by relying on a moderately sized random subsample of the data in each step of the algorithm. The function value of the target density is thus replaced by an approximation. Naturally, subsampling and alternative attempts at “cutting the Metropolis-Hastings budget” induce additional biases. These biases can lead to dramatic changes in the properties of the algorithms as discussed in .

We provide perturbation bounds based on Wasserstein distances, which lead to flexible quantitative estimates of the biases of approximate MCMC methods. Our first main result is the Wasserstein perturbation bound of Theorem 3.1. Under a Wasserstein ergodicity assumption, explained in Section 2, it provides an upper bound on the distance of the $n$th step distribution between an ideal and an approximating Markov chain in terms of the difference between their one-step transition probabilities. The result is well-suited for applications on a non-compact state space, since the difference of the one-step transition probabilities is measured by a weighted supremum with respect to a suitable Lyapunov function. For an autoregressive model, we show in Section 4.1 that the resulting perturbation bound cannot be improved in general. As a consequence of the Wasserstein approach we also obtain perturbation estimates for geometrically ergodic Markov chains. We first adapt our Wasserstein perturbation bound to this setting. Then, as a second main result, Theorem 3.2, we prove a refined estimate for geometrically ergodic chains where the perturbation is measured by a weighted total variation distance. Our perturbation bounds, and earlier ones in , establish a direct connection between an exponential convergence property for Markov chains and their robustness to perturbations. In particular, fast convergence to stationarity implies insensitivity to perturbations in the transition probabilities. Geometric ergodicity has been studied extensively in the MCMC literature. Thus, our estimates can be used in combination with many existing convergence results for MCMC algorithms. In Section 4, we illustrate the applicability of both theorems by generalizing recent findings on approximate Metropolis-Hastings algorithms from and on noisy Langevin algorithms for Gibbs random fields from .

We refer to for an overview of the classical literature on perturbation theory for Markov chains. However, as Stuart and Shardlow observed in , the classical assumptions on the perturbation might be too restrictive for many interesting applications. As a consequence, they develop a perturbation theory for geometrically ergodic Markov chains which requires to control perturbations of iterated transition kernels in a weaker sense. In our bounds for geometrically ergodic Markov chains, we have similar flexibility in the perturbation due to the Lyapunov-type stability condition, and require only a control on the errors of one-step transition kernels.

Mitrophanov, in , considers uniformly ergodic Markov chains and provides the best estimates in those settings. In the geometrically ergodic case, there are further related results, see and the references therein. Compared to , our focus is on non-asymptotic estimates with explicit constants, while their main focus is on qualitative results such as inheritance of geometric ergodicity by the perturbation. Earlier related results on perturbations induced by floating-point roundoff errors are shown in .

Finally, let us point out that our paper is complementary to the work of Pillai and Smith who also present Wasserstein perturbation bounds for Markov chains. When moving beyond the uniformly ergodic Markov chain case, an important challenge is to handle the issue that in many applications suprema of relevant quantities over the whole state space are infinite. The authors of guarantee finiteness of supremum norms by restricting attention to subsets of the state space. Their bounds thus involve exit probabilities from these subsets. Our approach circumvents these issues by relying on Lyapunov-type stability conditions for the approximate algorithm.

Wasserstein ergodicity

Let $G$ be a Polish space and $\mathcal{B}(G)$ be the corresponding Borel $\sigma$-algebra. Let $d$ be a metric, possibly different from the one which makes the space Polish, which is assumed to be lower semi-continuous with respect to the product topology of $G$. Let $\mathcal{P}$ be the set of all Borel probability measures on $(G,\mathcal{B}(G))$. Then, we define the Wasserstein distance of $\nu,\mu\in\mathcal{P}$ by

which leads to the well-known duality formula

For details we refer to [45, Chapter 1.2]. By $\delta_{x}$ we denote the probability measure concentrated at $x$. Hence $W(\delta_{x},\delta_{y})=d(x,y)$ is finite for $x,y\in G$.

Let $P$ be a transition kernel on $(G,\mathcal{B}(G))$ which defines a linear operator $P\colon\mathcal{P}\to\mathcal{P}$ given by

with $Pf(x)=\int_{G}f(y)P(x,{\rm d}y)$ whenever one of the integrals exist, see for example [40, Lemma 3.6]. Now, by

we define the generalized ergodicity coefficient of transition kernel $P$. This coefficient can be understood as a generalized Dobrushin ergodicity coefficient, see . Dobrushin himself called $\tau(P)$ the Kantorovich norm of $P$, see [10, formula (14.34)]. Finally, $\tau(P)$ also provides a lower bound of the coarse Ricci curvature of $P$ introduced in .

Two essential properties of the ergodicity coefficient are submultiplicativity and contractivity, see [10, Proposition 14.3 and Proposition 14.4].

For two transition kernels $P$ and $\widetilde{P}$ on $(G,\mathcal{B}(G))$ and $\mu,\nu\in\mathcal{P}$, we have

As an immediate consequence of this contractivity, we obtain the following corollary.

Let $P$ be a transition kernel with stationary distribution $\pi$, i.e. $\pi P=\pi$, and assume for some (and hence any) $x_{0}\in G$ it holds that $\int_{G}d(x_{0},x)\,{\rm d}\pi(x)<\infty$. Then

Because of the assumption $\int_{G}d(x_{0},x)\,{\rm d}\pi(x)<\infty$ we have that $W(\delta_{x},\pi)$ is finite for any $x\in G$. Thus, the assertion follows by Proposition 2.1 and stationarity of $\pi$. ∎

For some special cases one also has an estimate of the form (2.2) in the other direction. To this end, consider the trivial metric $d(x,y)=2\cdot\mathbf{1}_{x\not=y}$ with indicator function

be the total variation norm of a signed measure $q$ on $G$. In this setting $W(\mu,\nu)=\left\|\mu-\nu\right\|_{\text{tv}}$. For $x,y\in G$ with $x\not=y$ we have $\left\|\delta_{x}-\delta_{y}\right\|_{\text{tv}}=d(x,y)=2$ so that

For the moment, let us assume that $P$ is uniformly ergodic, that is, there exist numbers $\rho\in[0,1)$ and $C\in(0,\infty)$ such that

An immediate consequence of the uniform ergodicity is that $\tau_{1}(P^{n})\leq C\rho^{n}$.

For the transition kernel $P$ there exist numbers $\rho\in[0,1)$ and $C\in(0,\infty)$ such that

For any probability measure $p_{0}\in\mathcal{P}$, a transition kernel $P$ with stationary distribution $\pi$ and $p_{n}=p_{0}P^{n}$ we have under the Wasserstein ergodicity condition that

Perturbation bounds

Similar as in [33, Theorem 3.1], we show quantitative bounds on the difference of $p_{n}$ and $\widetilde{p}_{n}$, but use the Wasserstein distance instead of total variation. Besides Assumption 2.1, the bounds depend on the difference of the initial distributions and on a suitably weighted one-step difference between $P$ and $\widetilde{P}$.

Let Assumption 2.1 be satisfied with the numbers $C\in(0,\infty)$ and $\rho\in[0,1)$, i.e., $\tau(P^{n})\leq C\rho^{n}$. Assume that there are numbers $\delta\in(0,1)$ and $L\in(0,\infty)$ and a measurable Lyapunov function $\widetilde{V}:G\rightarrow[1,\infty)$ of $\widetilde{P}$ such that

with ${\widetilde{p}_{0}}(\widetilde{V})=\int_{G}\widetilde{V}(x)\,{\rm d}{\widetilde{p}_{0}}(x)$. Then

so that we obtain $W(\widetilde{p}_{i}P,\widetilde{p}_{i}\widetilde{P})\leq\gamma\kappa$. By this fact we have

Then, by (3.3), (3.4) and the triangle inequality of the Wasserstein distance we have

Finally, by (2.4) we obtain $\sum_{i=0}^{n-1}\tau(P^{i})\leq\frac{C(1-\rho^{n})}{1-\rho},$ which allows us to complete the proof. ∎

The parameter $\kappa$ is an upper bound on $\widetilde{p}_{i}(\widetilde{V})$. It can be interpreted as a measure for the stability of the perturbed Markov chain. The parameter $\gamma$ quantifies with a weighted supremum norm the one-step difference between $P$ and $\widetilde{P}$. The use of the Lyapunov function increases the flexibility of the resulting estimate, since larger values of $\widetilde{V}$ compensate larger values of the Wasserstein distance between the kernels. Notice that the existence of a Lyapunov function satisfying (3.1) is weaker than assuming $\widetilde{V}$-uniform ergodicity of $\widetilde{P}$ since it is not associated with a small set condition. In particular, the condition is satisfied for any $\widetilde{P}$ with the trivial choice $\widetilde{V}(x)=1$ for all $x\in G$, see Corollary 3.2. As we will see in Section 4, allowing for non-trivial choices of $\widetilde{V}$ considerably increases the applicability of our results.

If $\widetilde{P}$ has a stationary distribution, say $\widetilde{\pi}\in\mathcal{P}$, as a consequence of the previous theorem, we obtain bounds on the difference between $\pi$ and $\widetilde{\pi}$.

Let the assumptions of Theorem 3.2 be satisfied. Assume that $\widetilde{P}$ has a stationary distribution $\widetilde{\pi}\in\mathcal{P}$ and let $W(\pi,\widetilde{\pi})$ be finite. Then

By Theorem 3.2 we obtain with $p_{0}=\pi$, $\widetilde{p}_{0}=\widetilde{\pi}$, the stationarity of the distributions $\pi$, $\widetilde{\pi}$ and by letting $n\to\infty$ that

By the Lyapunov condition and [16, Proposition 4.24], it holds that

which leads to $\kappa\leq L/(1-\delta)$ and finishes the proof. ∎

It may seem artificial to assume $W(\pi,\widetilde{\pi})<\infty$ but this is needed for the limit argument in the proof. This condition is often satisfied a priori. For example, it holds if the metric is bounded, i.e., $\sup_{x,y\in G}d(x,y)$ is finite, or, more generally, if the distributions $\pi$ and $\widetilde{\pi}$ possess a first moment in the sense that there exist $x_{0},\widetilde{x}_{0}\in G$ such that

As pointed out in Remark 3.1, we do not need to impose condition (3.1) to obtain a non-trivial perturbation bound:

Assume that Assumption 2.1 holds with the numbers $C\in(0,\infty)$ and $\rho\in[0,1)$, i.e., $\tau(P^{n})\leq C\rho^{n}$, and let

The statement follows by Theorem 3.1 with $\widetilde{V}(x)=1$ and $L=1-\delta$. ∎

For the trivial metric $d(x,y)=2\cdot\mathbf{1}_{x\not=y}$ the last corollary states essentially the result of [33, Theorem 3.1], where instead of the general Wasserstein distance the total variation distance is used. There, the bound’s dependence on $C$ and $\rho$ can be further improved by using the a priori bound $\tau_{1}(P^{n})\leq 1$ in addition to uniform ergodicity. For another metric $d$ such an a priori bound is in general not available.

Table 1 provides a detailed comparison between our Theorem 3.1 and the related Wasserstein perturbation result of Pillai and Smith, [35, Lemma 3.3]. An important ingredient in their estimate is a set $\widehat{G}\subseteq G$ which can be interpreted as the part of $G$ where both Markov chains remain with high probability. When a good uniform upper bound on $W(\delta_{x}P,\delta_{x}\widetilde{P})$ for all $x\in G$ is available, we can choose $\widehat{G}=G$ in [35, Lemma 3.3] and $\widetilde{V}(x)=1$ in Theorem 3.1. In that case, both results essentially simplify to Corollary 3.2. The results become entirely different when such a bound is not available or too rough. In our estimate, one then needs a non-trivial Lyapunov function for $\widetilde{P}$ and a uniform upper bound on $W(\delta_{x}P,\delta_{x}\widetilde{P})/\widetilde{V}(x)$. To apply their estimate, one needs a uniform bound on $W(\delta_{x}P,\delta_{x}\widetilde{P})$ for all $x\in\widehat{G}$. In addition, a bound on $\pi(G\setminus\widehat{G})$, Lyapunov functions and estimates of the exit probabilities from $\widehat{G}$ of both Markov chains need to be available. Finally, while [35, Lemma 3.3] requires slightly more regularity on the Lyapunov function, contractivity of the unperturbed transition kernel $P$ (with $C=1$) is not needed on the whole state space but only on $\widehat{G}$.

2 Perturbation bounds for geometrically ergodic Markov chains

In this section, we derive general perturbation bounds for geometrically ergodic Markov chains. First, we recall some results from , and which are helpful to apply our Wasserstein perturbation bounds in the geometrically ergodic case. Then we present the new estimates:

Corollary 3.3 is an application of Theorem 3.1 with Wasserstein distances replaced by $V$-norms of differences between measures.

In Corollary 3.4, we show that having a Lyapunov function $V$ for $P$ is sufficient for our bounds if the transition kernels $P$ and $\widetilde{P}$ are sufficiently close (in a suitable sense).

In Theorem 3.2, we provide a quantitative perturbation bound which still applies if we can only control the total variation distance between $P(x,\cdot)$ and $\widetilde{P}(x,\cdot)$. To measure the perturbation in such a weak sense is new for geometrically ergodic Markov chains.

A transition kernel $P$ with stationary distribution $\pi$ is called geometrically ergodic if there is a constant $\rho\in[0,1)$ and a measurable function $C\colon G\to(0,\infty)$ such that for $\pi$-a.e. $x\in G$ we have

For $\phi$-irreducible and aperiodic Markov chains, it is well known that geometric ergodicity is equivalent to $V$-uniform ergodicity, see [36, Proposition 2.1]. Namely, if $P$ is geometrically ergodic, then there exists a $\pi$-a.e. finite measurable function $V\colon G\to[1,\infty]$ with finite moments with respect to $\pi$ and there are constants $\rho\in[0,1)$ and $C\in(0,\infty)$ such that

The following result establishes the connection between $V$-norms and certain Wasserstein distances. It is basically due to Hairer and Mattingly , see also .

Assume that $V$ is lower semi-continuous on $G$. For $x,y\in G$, let us define the metric

Then, for any $\mu,\nu\in\mathcal{P}$ we have

where $W_{d_{V}}$ denotes the Wasserstein distance based on the metric $d_{V}$.

Lower semi-continuity of $V$ implies lower semi-continuity of $d_{V}$, which leads to the duality formula (2.1) by [45, Theorem 1.14]. We thus impose the standing assumption of lower semi-continuity of $V$ whenever we speak of $V$-uniform ergodicity in the following. In principle, this requirement can be removed and (3.8) remains true, but we do not go into further detail in that direction. In applications, this is typically not restrictive since $V$ is continuous anyway.

By similar arguments as in the proof of [26, Theorem 1.1] we observe that (3.7) implies a suitable upper bound on

If (3.7) is satisfied for the transition kernel $P$, then $\tau_{V}(P^{n})\leq C\rho^{n}.$

For any positive real numbers $a_{1},a_{2},b_{1},b_{2}$ we have the following elementary inequality

Now, by using (3.7) we obtain the assertion. ∎The lemmas above and Theorem 3.1 lead to the following new perturbation bound for geometrically ergodic Markov chains.

Let $P$ be $V$-uniformly ergodic, i.e., there are constants $\rho\in[0,1)$ and $C\in(0,\infty)$ such that

We also assume that there are numbers $\delta\in(0,1)$ and $L\in(0,\infty)$ and a measurable Lyapunov function $\widetilde{V}:G\rightarrow[1,\infty)$ of $\widetilde{P}$ such that

with $\widetilde{p_{0}}(\widetilde{V})=\int_{G}\widetilde{V}(x)\,{\rm d}{\widetilde{p}_{0}}(x)$. Then

In [41, Theorem 3.1], a related perturbation bound is proven. The convergence property of the unperturbed transition kernel is slightly weaker than our $V$-uniform ergodicity, but also based on a kind of Lyapunov function. More restrictively, there it is assumed that the difference of $P^{n}$ and $\widetilde{P}^{n}$ for all $n>0$ can be controlled. In addition, the perturbation error is measured with a weight given by the same Lyapunov function as in the convergence property of $P$, but by taking a supremum over a subset of test functions. With our approach we can take the supremum over all test functions and obtain similar estimates by setting $p_{0}=\pi$.

The next corollary demonstrates how the Lyapunov function of $\widetilde{P}$ can be replaced by a Lyapunov function of $P$, provided that the distance between the transition kernels is sufficiently small. Notice that assuming the existence of a Lyapunov function of $P$ in addition to the $V$-uniform ergodicity is a definition of constants rather than an additional requirement, see, e.g., .

Let $P$ be $V$-uniformly ergodic, i.e., there are constants $\rho\in[0,1)$ and $C\in(0,\infty)$ such that

Moreover, $V\colon G\to[1,\infty)$ is a measurable Lyapunov function of $P$, such that

with constants $\delta\in(0,1)$ and $L\in(0,\infty)$. Let

with $\widetilde{p_{0}}(V)=\int_{G}V(x)\,{\rm d}{\widetilde{p}_{0}}(x)$. If $\gamma+\delta<1$, then

which implies (3.14). The assertion follows by the assumption that $\delta+\gamma<1$ and an application of Corollary 3.3. ∎

For discrete state spaces and under the requirement $p_{0}=\widetilde{p}_{0}$, a result similar to the previous corollary is obtained in [21, Theorem 3, Corollary 3]. The authors of replace our constant $\kappa$ by $\max_{0\leq i\leq n}\widetilde{p}_{i}(V)$. This we could do as well, see the proof of Theorem 3.1.

It is easily seen that $(B_{V},\left|\cdot\right|_{V})$ is a normed linear space. In the setting of Corollary 3.3, we have

In Corollary 3.4, the more restrictive case $V=\widetilde{V}$ is considered. The corresponding operator norm $\VERT P-\widetilde{P}\VERT_{B_{V}\to B_{V}}$ appears in classical perturbation theory for Markov chains, see . But as discussed in [41, p. 1126] and it might be too restrictive to measure the perturbation with this operator norm for $V=\widetilde{V}$.

By relying, e.g., on [28, Proposition 2] we have some flexibility in the choice of $V$. There it is shown that, for $r\in(0,1)$, $V$-uniform ergodicity implies $V^{r}$-uniform ergodicity. This leads to less favorable constants in the $V^{r}$-uniform ergodicity of $P$, but can relax the requirements on the similarity of $P$ and $\widetilde{P}$. Namely, with a Lyapunov function $\widetilde{V}$ of $\widetilde{P}$ we can apply Corollary 3.3 with a $V^{r}$-uniformly ergodic $P$ and $\gamma=\VERT P-\widetilde{P}\VERT_{B_{V^{r}}\to B_{\widetilde{V}}}$.

Unfortunately, this approach breaks down for $r=0$. To see this, notice that $V^{r}$-uniform ergodicity with $r=0$ is just uniform ergodicity which is not implied by geometric ergodicity. The next theorem overcomes this limitation by separating the two roles of the function $V$ in the previous perturbation bounds. Roughly, we set $V=1$ in the sense that we measure the distances between $P$ and $\widetilde{P}$ as well as between $p_{n}$ and $\widetilde{p}_{n}$ in the total variation distance. At the same time, we set $V=\widetilde{V}$ in the sense that we assume $P$ is $\widetilde{V}$-uniformly ergodic with Lyapunov function $\widetilde{V}$.

Let $P$ be $\widetilde{V}$-uniformly ergodic, i.e., there are constants $\rho\in[0,1)$ and $C\in(0,\infty)$ such that

Moreover, $\widetilde{V}\colon G\to[1,\infty)$ is a measurable Lyapunov function of $\widetilde{P}$ and $P$, such that

with constants $\delta\in(0,1)$ and $L\in(0,\infty)$. Let

with $\widetilde{p_{0}}(\widetilde{V})=\int_{G}\widetilde{V}(x)\,{\rm d}{\widetilde{p}_{0}}(x)$. Then, for $\gamma\in(0,\exp(-1))$ we have

From the proof of Theorem 3.2 we know that

Fix a real number $r\in(0,1)$ and let $s=1-r$. By considering (2.3) one can see that $\tau_{1}(P)\leq 1$. This leads to

For $\gamma\in(0,\exp(-1))$, we can choose the numbers $r=1+\log(\gamma)^{-1}$ and $s=\log(\gamma^{-1})^{-1}$. This yields $\gamma^{r}=\exp(1)\gamma$ and the proof is complete. ∎

Let $\widetilde{\pi}\in\mathcal{P}$ be a stationary distribution of $\widetilde{P}$. Notice that by the assumption that $\widetilde{V}$ is Lyapunov function of $\widetilde{P}$ and [16, Proposition 4.24] it follows that $\widetilde{\pi}(\widetilde{V})\leq L/(1-\delta)$. Further, by the $\widetilde{V}$-uniform ergodicity of $P$ we also know that $\pi(\widetilde{V})$ is finite. Thus,

Now, by Theorem 3.2 we can bound $\left\|\pi-\widetilde{\pi}\right\|_{\text{tv}}$ with $p_{0}=\pi$, $\widetilde{p}_{0}=\widetilde{\pi}$ and by letting $n\to\infty$. We obtain

In the setting of Theorem 3.2, we can also interpret $\gamma$ as an operator norm. Namely,

Applications

We illustrate our perturbation bounds in three different settings. We begin with studying an autoregressive process also considered in . After this, we show quantitative perturbation bounds for approximate versions of two prominent MCMC algorithms, namely the Metropolis-Hastings and stochastic Langevin algorithms.

and it is well known that there exists a stationary distribution, say $\pi_{\alpha}$, of $P_{\alpha}$.

leads to $\tau(P_{\alpha}^{n})\leq\left|\alpha\right|^{n}$. Similarly, one obtains

and $p_{\alpha,n}=p_{0}P^{n}_{\alpha}$, $\widetilde{p}_{\widetilde{\alpha},n}=\widetilde{p}_{0}P^{n}_{\widetilde{\alpha}}$. Then, inequality (3.2) of Theorem 3.1 gives

and for $p_{0}=\widetilde{p}_{0}$ we have

From the previous two inequalities one can see that if $\widetilde{\alpha}$ is sufficiently close to $\alpha$, then the distance of the distribution $p_{\alpha,n}$ and $\widetilde{p}_{\widetilde{\alpha},n}$ is small. Let us emphasize here that we provide an explicit estimate rather than an asymptotic statement.

The dependence on $\left|\alpha-\widetilde{\alpha}\right|$ in the previous inequality cannot be improved in general. To see this, let us assume that $X_{0,\alpha}$ and $X_{0,\widetilde{\alpha}}$ are real-valued random variables with distribution $\pi_{\alpha}$ and $\pi_{\widetilde{\alpha}}$, respectively. Then, because of the stationarity we have that $X_{1,\alpha}=\alpha X_{0,\alpha}+Z_{1}$ and $X_{1,\widetilde{\alpha}}=\widetilde{\alpha}X_{0,\widetilde{\alpha}}+Z_{1}$ are also distributed according to $\pi_{\alpha}$ and $\pi_{\widetilde{\alpha}}$, respectively. Thus

Let us now discuss the application of Corollary 3.4 and Theorem 3.2. Under the additional assumption that $\mu$, the distribution of $Z_{1}$, has a Lebesgue density $h$, it is shown in [15, Section 4] that the autoregressive model (4.1) is also $\widetilde{V}$-uniformly ergodic. Precisely, there is a constant $C\geq 1$ such that

Moreover, from [13, Example 1] we know that

does not go to $0$ when $\widetilde{\alpha}\downarrow\alpha$. Hence, Corollary 3.4 cannot quantify for small $|\widetilde{\alpha}-\alpha|$ whether the $n$th step distributions are close to each other. However, also in [13, Example 1] it is proven that

The first summand on the right hand side we can bound by

and similarly for the second summand. Using that $F(a)=F(-a)$, we obtain $F(a)\leq 2|a|\,h_{\max}$. Finally, by substitution we can write

For simplicity set $p_{0}=\widetilde{p}_{0}$ and assume that $h_{\max}\leq 1$ as well as $\left|\alpha-\widetilde{\alpha}\right|\in(0,\exp(-1)/2)$. Then, Theorem 3.2 implies

2 Approximate Metropolis-Hastings algorithms

We apply our perturbation results to the approximate (or noisy) Metropolis-Hastings algorithms analyzed in . We assume either that the unperturbed transition kernel of the Metropolis-Hastings algorithm satisfies the Wasserstein ergodicity condition stated in Assumption 2.1 or is geometrically ergodic. In particular, we do not assume that the transition kernel is uniformly ergodic. Let $\pi$ be a probability distribution on $(G,\mathcal{B}(G))$ and assume that we are interested in sampling realizations from this distribution. Let $Q$ be a transition kernel which serves as the proposal for the Metropolis-Hastings algorithm. From [44, Proposition 1] we know that there exists a set $S\subset G\times G$ such that we can define the “acceptance ratio” for $(x,y)\in G\times G$ as

Then, let the acceptance probability be $\alpha(x,y)=\min\{1,r(x,y)\}$. With this notation the Metropolis-Hastings algorithm defines a transition kernel

A single transition from $X_{n}$ to $X_{n+1}$ of the Metropolis-Hastings algorithm works as follows:

Draw a sample $Y\sim Q(X_{n},\cdot)$ and $U\sim\mbox{Unif}$ independently, call the result $y$ and $u$;

Set $r:=r(X_{n},y)$, with the ratio $r(\cdot,\cdot)$ defined in (4.5);

If $u<r$, then accept the proposal, and set $X_{n+1}:=y$, else reject the proposal and set $X_{n+1}:=X_{n}$.

A single transition from $\widetilde{X}_{n}$ to $\widetilde{X}_{n+1}$ works as follows:

Draw a sample $Y\sim Q(\widetilde{X}_{n},\cdot)$ and $U\sim\mbox{Unif}$ independently, call the result $y$ and $u$;

Draw a sample $R\sim\mu_{\widetilde{X}_{n},y,u}$, call the result $\widetilde{r}$;

If $u<\widetilde{r}$, then accept the proposal, and set $\widetilde{X}_{n+1}:=y$, else reject the proposal and set $\widetilde{X}_{n+1}:=\widetilde{X}_{n}$.

and the transition kernel of such a Markov chain is still of the form (4.6) with $\alpha(x,y)$ substituted by $\widetilde{\alpha}(x,y)$, i.e., it is given by $P_{\widetilde{\alpha}}$. The following results hold in the slightly more general case where $\widetilde{\alpha}(x,y)$ is any approximation of the acceptance probability $\alpha(x,y)$.

The next lemma provides an estimate for the Wasserstein distance between transition kernels of the form (4.6) in terms of the acceptance probabilities.

Let $Q$ be a transition kernel on $(G,\mathcal{B}(G))$ and let $\alpha\colon G\times G\rightarrow$ and $\widetilde{\alpha}\colon G\times G\rightarrow$ be measurable functions. By $P_{\alpha}$ and $P_{\widetilde{\alpha}}$ we denote the transition kernels of the form (4.6) with acceptance probabilities $\alpha$ and $\widetilde{\alpha}$. Then, for all $x\in G$, we have

with $\mathcal{E}(x,y)=|\alpha(x,y)-\widetilde{\alpha}(x,y)|$.

By the use of the dual representation of the Wasserstein distance it follows that

∎By the previous lemma and Theorem 3.1, we obtain the following Wasserstein perturbation bound for the approximate Metropolis-Hastings algorithm.

Let $Q$ be a transition kernel on $(G,\mathcal{B}(G))$ and let $\alpha\colon G\times G\rightarrow$ and $\widetilde{\alpha}\colon G\times G\rightarrow$ be measurable functions. By $P_{\alpha}$ and $P_{\widetilde{\alpha}}$ we denote the transition kernels of the form (4.6) with acceptance probabilities $\alpha$ and $\widetilde{\alpha}$. Let the following conditions be satisfied:

Assumption 2.1 holds for the transition kernel $P_{\alpha}$, i.e., $\tau(P_{\alpha}^{n})\leq C\rho^{n}$ for $\rho\in[0,1)$ and $C\in(0,\infty)$.

There are numbers $\delta\in(0,1)$, $L\in(0,\infty)$ and a measurable Lyapunov function $\widetilde{V}:G\rightarrow[1,\infty)$ of $P_{\widetilde{\alpha}}$, i.e.,

Let $\mathcal{E}(x,y)=|\alpha(x,y)-\widetilde{\alpha}(x,y)|$ and assume that

Then, for any $p_{0}\in\mathcal{P}$ and finite $p_{0}(\widetilde{V})=\int_{G}\widetilde{V}(x){\rm d}p_{0}(x)$ we have

where $\kappa=\max\left\{p_{0}(\widetilde{V}),\frac{L}{1-\delta}\right\}$.

Let us point out several aspects of condition (4.7). Recall that (4.7) is always satisfied with $\widetilde{V}(x)=1$ for all $x\in G$. However, in this case it seems more difficult to control $\gamma$. If some additional knowledge in form of a Lyapunov function $V\colon G\to[1,\infty)$ of $P_{\alpha}$, i.e., $P_{\alpha}V(x)\leq\delta V(x)+L$ for some $\delta\in(0,1)$ and $L\in(0,\infty)$, is available, then a non-trivial candidate for $\widetilde{V}$ is $V$. For sufficiently small

Then, $P_{\widetilde{\alpha}}V(x)\leq(\delta+\delta_{V})V(x)+L$ and whenever $\delta+\delta_{V}<1$ it is clear that condition (4.7) is verified.

To highlight the usefulness of a non-trivial Lyapunov function, we consider the following scenario which is related to a local perturbation of an independent Metropolis-Hastings algorithm.

Let us assume that for $P_{\alpha}$ Assumption 2.1, as formulated in Corollary 4.1, is satisfied. For some probability measure $\mu$ on $(G,\mathcal{B}(G))$ define $Q(x,\cdot)=\mu$ and $p_{0}=\widetilde{p}_{0}=\mu$. For $\widetilde{G}\subseteq G$ let

Hence, for $x\in\widetilde{G}$ the transition kernel $P_{\widetilde{\alpha}}(x,\cdot)$ accepts any proposed state and for $x\not\in\widetilde{G}$ we have $P_{\widetilde{\alpha}}(x,\cdot)=P_{\alpha}(x,\cdot)$. It is easily seen that $\mathcal{E}(x,y)\leq\mathbf{1}_{\widetilde{G}}(x)$. For arbitrary $R>0$ and $r\in(0,1)$ set $\widetilde{V}(x)=1+R\mathbf{1}_{\widetilde{G}}(x)$ and note that

The last inequality of the previous formula follows by distinguishing the cases $x\in\widetilde{G}$ and $x\not\in\widetilde{G}$. Define $D(\widetilde{G})=\sup_{x\in\widetilde{G}}\int_{G}d(x,y)\mu({\rm d}y)$ and observe

for arbitrary $R\in(0,\infty)$ and $r\in(0,1)$. Under the assumption that $D(\widetilde{G})$ is finite and letting $R\to\infty$ as well as $r\downarrow 0$ we obtain

which tells us that basically $\mu(\widetilde{G})$ measures the difference of the distributions. A small perturbation set $\widetilde{G}$ with respect to $\mu$, thus implies a small bias. In contrast, with the trivial Lyapunov function $\widetilde{V}=1$, and if there is $(x,y)\in\widetilde{G}\times G$ such that $\alpha(x,y)=0$, we only obtain

The resulting upper bound on $W(p_{0}P_{\alpha}^{n},p_{0}P_{\widetilde{\alpha}}^{n})$ will typically be bounded away from zero regardless of the set $\widetilde{G}$.

The constant $\gamma$ essentially depends on the distance $d(x,y)$ and the difference of the acceptance probabilities $\mathcal{E}(x,y)$. By applying the Cauchy-Schwarz inequality to the numerator of $\gamma$, we can separate the two parts, i.e.,

If both integrals remain finite we see that an appropriate control of $\mathcal{E}(x,y)$ suffices for making the constant $\gamma$ small.

By using a Hoeffding-type bound, in Bardenet et al. [3, Lemma 3.1.] it is shown that for their version of the approximate Metropolis-Hastings algorithm with adaptive subsampling the approximation error $\mathcal{E}(x,y)$ is bounded uniformly in $x$ and $y$ by a constant $s>0$. Moreover, $s$ can be chosen arbitrarily small for the implementation of the algorithm.

Now we consider the case where the unperturbed transition kernel $P_{\alpha}$ is geometrically ergodic. Motivated by Remark 4.2, we also assume that $\mathcal{E}(x,y)\leq s$ for a sufficiently small number $s>0$. The following corollary generalizes a main result of Bardenet et al. [3, Proposition 3.2] to the geometrically ergodic case.

Let $Q$ be a transition kernel on $(G,\mathcal{B}(G))$ and let $\alpha\colon G\times G\rightarrow$ and $\widetilde{\alpha}\colon G\times G\rightarrow$ be measurable functions. By $P_{\alpha}$ and $P_{\widetilde{\alpha}}$ we denote the transition kernels of the form (4.6) with acceptance probabilities $\alpha$ and $\widetilde{\alpha}$. Let the following conditions be satisfied:

The unperturbed transition kernel $P_{\alpha}$ is $V$-uniformly ergodic, that is,

for numbers $\rho\in[0,1)$, $C\in(0,\infty)$ and a measurable function $V\colon G\rightarrow[1,\infty)$. Moreover, $V$ is a Lyapunov function of $P_{\alpha}$, i.e.,

for numbers $\delta\in(0,1)$ and $L\in(0,\infty)$.

A uniform bound $s>0$ on the difference of the acceptance probabilities is given, that is, for all $x,y\in G$, we have

If $s<(1-\delta)/\lambda$, then, for any $p_{0}\in\mathcal{P}$ with finite $\kappa=\max\left\{p_{0}(V),\frac{L}{1-\delta-\lambda s}\right\}$ we have

We consider the metric $d_{V}$, defined in Lemma 3.1, set $V=\widetilde{V}$ and use $\mathcal{E}(x,y)\leq s$ so that it is easily seen that the constant $\gamma$ from Corollary 4.1 satisfies $\gamma\leq s\lambda$. From the proof of Corollary 3.4, we know that $V$ is a Lyapunov function of $P_{\widetilde{\alpha}}$ provided that $\gamma+\delta<1$. Thus, we have

Now if $s<(1-\delta)/\lambda$, then $\delta+\lambda s<1$ and the assertion follows from Corollary 4.1 by writing the Wasserstein distances in terms of $V$-norms as in Section 3.2. ∎

Without $V(x)$ in the denominator, i.e., if we had relied on Corollary 3.2 instead of Theorem 3.1, the constant $\lambda$ would often be infinite. Consider the following toy example: Let $\pi$ be the exponential distribution with density $\exp(-x)$ on $G=[0,\infty)$ and assume that $Q(x,{\rm d}y)$ is a uniform proposal with support $[x-1,x+1]$. With $V(x)=\exp(x)$ it is well known that the Metropolis-Hastings algorithm is $V$-uniformly ergodic, see or [37, Example 4]. In this example

whereas $\int_{x-1}^{x+1}\exp(y){\rm d}y$ is unbounded in $x$. Notice that $\lambda$ only depends on the unperturbed Markov chain so that a bound on $\lambda$ can be combined with any approximation.

Let $P_{\widetilde{\alpha}}$ and $P_{\alpha}$ be $\phi$-irreducible and aperiodic. Then, one can prove under the assumptions of Corollary 4.2 that $P_{\widetilde{\alpha}}$ is $V$-uniformly ergodic if $s$ is sufficiently small. To see this, note that by [31, Theorem 16.0.1] the $V$-uniform ergodicity of $P_{\alpha}$ implies that $P_{\alpha}$ satisfies their drift condition (V4). By the arguments stated in the proof of Corollary 3.4, one obtains that $P_{\widetilde{\alpha}}$ also satisfies (V4) for sufficiently small $s$ and this implies $V$-uniform ergodicity. In this case, clearly $P_{\widetilde{\alpha}}$ possesses a stationary distribution, say $\widetilde{\pi}$, and

The previous inequality follows by (3.5) and the fact that

Here the finiteness of $\pi(V)$ follows by the $V$-uniform ergodicity of $P$ and $\widetilde{\pi}(V)\leq L/(1-\delta-\lambda s)$ follows by (4.10) and [16, Proposition 4.24].

3 Noisy Langevin algorithm for Gibbs random fields

An alternative to the Metropolis-Hastings algorithm is the Langevin algorithm, see . Unfortunately, in its implementation one needs the gradient of the density of the target distribution. To overcome this problem, different approximate Langevin algorithms have been proposed and studied, see .

where the prior density $p(\theta)$ is the Lebesgue density of the normal distribution $\mathcal{N}(0,\sigma_{p}^{2})$ with $\sigma_{p}>0$.

In general $\pi_{y}$ is not a stationary distribution of $P_{\sigma}$, but there exists a stationary distribution (see Proposition 4.1 below), say $\pi_{\sigma}$, which is close to $\pi_{y}$ depending on $\sigma$. Let $z(\theta)=\sum_{y\in\mathcal{Y}^{M}}\exp(\theta\,s(y))$ then, by the definition of $\pi_{y}$ we have

A single transition from $\widetilde{X}_{n}$ to $\widetilde{X}_{n+1}$ works as follows:

Draw $Z_{n}\sim\mathcal{N}(0,\sigma^{2})$, independent from step 1., call the result $z_{n}$. Set

From [2, Lemma 3] and by applying arguments of , we obtain the following facts about the noisy Langevin algorithm.

the function $V$ is a Lyapunov function for $P_{\sigma}$ and $P_{\sigma,N}$. We have

with $\delta=1-\frac{\sigma^{2}}{4\sigma_{p}^{2}}$, $L=\sigma+\sigma^{2}\left\|s\right\|_{\infty}+\frac{\sigma^{2}}{2\sigma_{p}^{2}}$ and the interval

the transition kernels $P_{\sigma}$ and $P_{\sigma,N}$ are $V$-uniformly ergodic.

for $N>4\max\left\{\|s\|_{\infty}^{2}\sigma^{4},\|s\|_{\infty}^{-3}\sigma^{-6}\right\}$ we have

Thus, the assertions from 1. to 3. are proven. The statement of 4. is a consequence of [2, Lemma 3]. There it is shown that for $N>4\|s\|_{\infty}^{2}\sigma^{4}$ it holds that

By using $\exp(\theta)-1\leq\theta\exp(\theta)$ and $N>4$ we further estimate the right-hand side by

Since $\log(N)\cdot N^{-1/3}<2$, we have the bound $K_{N,s,\sigma}\leq\exp(1)$ provided that $4N^{2/3}\|s\|_{\infty}^{2}\sigma^{4}\geq 2$ which follows from $N\geq\|s\|_{\infty}^{-3}\sigma^{-6}$. The assertion of (4.13) follows now by a simple calculation. ∎

By using the facts collected in the previous proposition, we can apply the perturbation bound of Theorem 3.2 and obtain a quantitative perturbation bound for the noisy Langevin algorithm.

We have by Proposition 4.1 that $P_{\sigma}$ is $V$-uniformly ergodic with $V(\theta)=1+\left|\theta\right|$, i.e., there are numbers $\rho\in[0,1)$ and $C\in(0,\infty)$ such that

Now, by combining Theorem 3.2 and Remark 3.8 with the results from Proposition 4.1 we obtain the result. ∎

We want to point out that the assumptions imposed are the same as in [2, Theorem 3.2], but instead of the asymptotic result we provide an explicit estimate. The numbers $\rho\in[0,1)$ and $C\in(0,\infty)$ are not stated in terms of the model parameters. In principle, these values can be derived from the drift condition (4.12) through [5, Theorem 1.1].

Acknowledgements

We thank Alexander Mitrophanov and the referees for their valuable comments which helped to improve the paper. D.R. was supported by the DFG Research Training Group 2088.

References