{'fraction_non_alphanumeric': 0.05602552686257049, 'fraction_numerical': 0.021853665776194717, 'mean_word_length': 3.892448720275912, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving the Poisson equation using coupled Markov chains

This article draws connections between unbiased estimators constructed from coupled Markov chains that meet exactly after a random number of iterations, and solutions of the Poisson equation. We first show how such pairs of chains can be employed to obtain unbiased estimators of pointwise evaluations of solutions of the Poisson equation. We then propose new estimators of the asymptotic variance of Markov chain ergodic averages. We formally study the proposed estimators under realistic assumptions on the meeting times of the coupled chains and on the existence of moments of test functions under the target distribution. We illustrate their behaviour in toy examples and in a more challenging setting of high-dimensional Bayesian regression.


Central Limit Theorem and the Poisson equation
Markov chain Monte Carlo (MCMC) methods form a convenient family of simulation techniques with many applications in statistics. With (X, X ) a measurable space and π a probability measure of statistical interest, MCMC involves simulation of a discrete-time, time-homogeneous Markov chain X = (X t ) t≥0 , with a π-invariant Markov transition kernel P and some initial distribution π 0 . Letting L p (π) = {f : π(|f | p ) < ∞}, where π(f ) = f (x)π(dx), the interest is to approximate π(h) for some function h ∈ L 1 (π).
In particular, after simulating the chain until time t, one may approximate an integral of interest π(h) via the average t −1 t−1 s=0 h(X s ). Under weak assumptions, such averages converge almost surely to π(h) as t → ∞, and under stronger but still realistic assumptions on π 0 , P and h that hold in many applications, they satisfy central limit theorems (CLTs), 1) where v(P, h) is the asymptotic variance associated with the Markov kernel P and the function h (see, e.g., Meyn & Tweedie 2009, Douc et al. 2018, Jones 2004. One standard route to proving a CLT is via a solution of the Poisson equation for h associated with P , i.e. any function g such that 2) where P is viewed as a Markov operator and I is the identity. In particular, if h, g ∈ L 2 (π) then (1. 1) holds with initial distribution π by Douc et al. (2018, Theorem 21.2.5), and v(P, h) = E π {g(X 1 ) − P g(X 0 )} 2 = 2π(h 0 · g) − π(h 2 0 ), (1. 3) where E π indicates that X 0 ∼ π. In the sequel, we focus on solutions g = g +c, where g = ∞ t=0 P t h 0 , and this solution in (1.3) leads to the more familiar expression v(P, h) = var π (h(X 0 )) + 2 ∞ t=1 cov π (h(X 0 ), h(X t )), (1.4) where the subscript π indicates that X 0 ∼ π.

Monte Carlo methods using the Poisson equation
Despite its convenience for theoretical analysis, the Poisson equation is typically not analytically solvable for most Markov chains and functions of interest, and consistent approximations have been lacking. As a result, the asymptotic variance v(P, h) is commonly estimated using batch-means or spectral variance estimators (see, e.g., Flegal & Jones 2010, Vats et al. 2018, Chakraborty et al. 2019. On the other hand, inconsistent approximations of the solution of the Poisson equation have been proposed for the purpose of variance reduction via control variates. This involves replacing h by h − φ in an ergodic average, where φ ∈ L 1 0 (π), so that π(h − φ) = π(h) and the limit of the MCMC estimator is unchanged, but the variance may be smaller with a judicious choice of φ. A convenient family of φ is {(I − P )f : f ∈ L 1 (π)} since φ ∈ L 1 0 (π) by construction and indeed an optimal choice of f is g. Approximations of g have been considered for this purpose in, e.g., Andradóttir et al. (1993), Henderson (1997), Dellaportas & Kontoyiannis (2012), Mijatović & Vogrinc (2018), Alexopoulos et al. (2022).
We aim to contribute to methodological aspects of the Poisson equation, with a focus on asymptotic variance estimation via (1.3). We develop novel approximations of the solution g, which are unbiased under fairly weak conditions, using coupled Markov chains. Couplings have often been used as formal techniques to analyze the marginal convergence of Markov chains (e.g. Jerrum 1998), but many couplings are also implementable, with applications to exact sampling (Propp & Wilson 1996), unbiased estimation (Glynn & Rhee 2014), convergence diagnostics (Johnson 1996, Biswas et al. 2019, or variance reduction (Neal & Pinto 2001). Here we find new uses for couplings of MCMC algorithms: first we propose unbiased estimators of evaluations of the solution g, and then we use them to approximate (1.3) in combination with unbiased MCMC , leading to estimators of v(P, h), which are unbiased and with finite variance under fairly mild conditions.

Coupled chains and fishy functions 2.1 Coupled Markov chains
For a given π-invariant Markov kernel P we consider the time-homogeneous, discrete-time Markov chain (X, Y ) = (X t , Y t ) t≥0 with Markov kernelP , which is a "faithful" coupling of P with itself (Rosenthal 1997), i.e. it satisfies P (x, y; A × X) = P (x, A),P (x, y; X × A) = P (y, A), A ∈ X , (2. 1) andP (x, x; {(x , y ) : x = y }) = 1 for all x ∈ X. We observe that X = (X t ) t≥0 and Y = (Y t ) t≥0 are both time-homogeneous, discrete-time Markov chains with kernel P . We use subscripts to denote the distribution of (X 0 , Y 0 ). For example P x,y is the law of (X, Y ) when (X 0 , Y 0 ) = (x, y), and P ν indicates When only one chain is referenced, e.g. X, we may similarly write P x and P µ to indicate X 0 = x and X 0 ∼ µ, respectively. An important random variable when using coupled Markov chains is the meeting time, which is the first time at which X and Y take the same value. SinceP is faithful, X t = Y t with probability 1 for all times t ≥ τ . Johnson (1998 and others have shown how to construct the coupled kernelP for various MCMC algorithms, such that τ has finite expectation; see Appendix B for further information. The main assumption throughout this paper is that τ has a finite κth moment, κ > 1, when (X 0 , Y 0 ) ∼ π ⊗ π, the product measure of π and itself. Assumption 1. The Markov kernel P is π-irreducible and for some κ > 1, E π⊗π [τ κ ] < ∞.
The assumption implies a polynomially decaying survival function P x,y (τ > t) of order κ and, conversely, one may verify the assumption by showing that P x,y (τ > t) decays polynomially with order s > κ with a dependence on (x, y) that is not too strong. The proof is provided in Section 4.1.

Unbiased approximation of fishy functions
In the Markov chain setting, solutions of the Poisson equation for h, or fishy functions for brevity, have been studied extensively (see, e.g., Duflo 1970, Glynn & Meyn 1996, Glynn & Infanger 2022. It is known, for example, that fishy functions exist under fairly weak conditions. In this paper we focus on a particular, and more restrictive, core expression for the fishy function for h.

Definition 3.
Let h ∈ L 1 (π), define h 0 = h − π(h) and the function Although g is not necessarily well-defined π-almost everywhere in general, Assumption 1 implies that for p ≥ 1, g ∈ L p (π) for h ∈ L m (π) when κ and m are sufficiently large in relation to p; see Theorem 13. If g is well-defined, then it is straightforward to deduce that (I − P )g = h 0 , so that g is indeed fishy. In addition, g + c is fishy for any constant c ∈ R. Moreover, if g ∈ L 1 (π) then π(g ) = 0, and by Douc et al. (Lemma 21.2.2, 2018) any other fishy function g is equal to g up to an additive constant. We can now define a family of fishy functions that will play a central role in what follows.

Definition 4.
For y ∈ X, the function g y : x → g (x) − g (y), (2.4) is fishy if g is well-defined, since g (y) is a constant. When y is fixed, in the sequel we may write g instead of g y .
A priori, links between couplings of Markov chains and fishy functions may not appear as obvious.
However, the definition of g y lends itself naturally to unbiased approximation via coupled Markov chains.
Definition 5. For x ∈ X, the proposed estimator of g y (x) is, with (X 0 , Y 0 ) = (x, y), 5) where (X, Y ) and τ are defined in Section 2.1. In particular, the dependence of G y (x) on x and y is that the law of (X, Y ) is P x,y in (2. 5), and so we may omit the subscript (x, y) when there is no ambiguity. We will also denote G y by G below, in places where the explicit mention of y is unnecessary.
The simple intuition behind Definition 5 is that we can equivalently write and upon justifying the interchange of expectation and infinite sum as we do in Section 4.2, then The random variable G y (x) may be simulated using Algorithm Output: coupled chains and the meeting time defined as the smallest t at which X t = Y t−L .
Remark 6. The above reasoning demonstrates how one can approximate the difference between two fishy function evaluations. In contrast, to unbiasedly approximate g (x) we would need to estimate g ν (x) and π(g ν ) in an unbiased manner, for an arbitrary distribution ν, and taking their difference.
An alternative estimator to Definition 5 could employ a random "truncation" variable as in Glynn & Rhee (2014) or Agapiou et al. (2018). This would allow the use of coupled chains that do not exactly meet, as long as the distance between them goes to zero fast enough. We do not pursue this generalization here.

Recovering Glynn-Rhee estimators
Through the Poisson equation we can also recover the unbiased estimation techniques of Glynn & Rhee (2014. We observe that we may rearrange (1.2) as where the left-hand side does not depend on x and g is any fishy function for h. It then seems natural to estimate π(h) by estimating the terms on the right-hand side. With g = g we may write where the right-hand side is a familiar quantity in the light of Glynn & Rhee (2014). In particular, one may define a lag 1 process (X t , Y t ) t≥0 , starting from X 0 = Y 0 = x, then sampling X 1 ∼ P (X 0 , ·), and iteratively sampling (X t+1 , Y t ) ∼P ((X t , Y t−1 ), ·) for t ≥ 1, whereP is as defined in Section 2.1.
Such (X, Y ) process is characterized by a lag L = 1 such that (X t+L , Y t ) t≥0 is Markov and X t+L = Y t almost surely for all t large enough.
} may be unbiased for π(h) under suitable assumptions. This same perspective on (2.7) suggests that for any x ∈ X we may define the equivalent and The initialization of the chains can be generalized from a point mass at x ∈ X to any probability distribution π 0 . Indeed, a re-arranged and integrated Poisson equation is π(h) = π 0 (h) + π 0 P (g) − π 0 (g), (2.9) and this suggests the following estimator of π(h), with essentially the same justification as above.
Definition 7. For h ∈ L 1 (π), denote the approximation of π(h) by where marginally X 0 ∼ µ, Y 0 ∼ µ and X 1 ∼ µP for some probability measure µ. We denote by γ the joint probability measure for (X 1 , Y 0 ), since this features in our analysis, noting that this is a coupling of µP and µ.
The estimator in Definition 7 is identical to the estimator denoted by H 0 (X, Y ) in Jacob et al.
By changing the initial distributions, Definition 7 also admits the estimators in Jacob et al. (2020) denoted by H k (X, Y ) for some k ∈ N, and the estimators H k:m (X, Y ) are obtained as averages of in Section 3. We also show that finiteness of moments of |H| can be transferred also to k-offset and L-lagged estimators with L ≥ 1. Appendix A provides background on Glynn-Rhee estimators and their variants.
Remark 8. There were already known links between bias of MCMC estimators and fishy functions. Kontoyiannis & Dellaportas (2009, Section 4) observe that the average t −1 t−1 s=0 h(X s ) has expectation t −1 t−1 s=0 P s h(x 0 ) given X 0 = x, thus g (x) represents the leading term in the asymptotic bias: Hence, if g ∈ L 1 0 (π) and an approximation of g = g + c is available, we may use the identity g = g − π(g) to approximate g . On the other hand, even without accounting for c, the behaviour of g may be informative about the asymptotic bias, e.g. how g(x) changes as |x| increases, as illustrated in Section 5.2.

Ergodic Poisson asymptotic variance estimator
We consider the task of estimating the asymptotic variance in the CLT, v(P, h) Vats et al. 2018, Chakraborty et al. 2019, and existing estimators are consistent when the length of each chain goes to infinity. Here we employ coupled Markov chains, as generated by Algorithm 2.1,to define new estimators of v(P,h): we start with a consistent but biased estimator before introducing unbiased estimators in the next section. We first re-express (1. 3) as 1) where v(π, h) is the variance of h(X) under X ∼ π. Expectations with respect to π can be consistently estimated from long MCMC runs. Using t steps after burn-in we define the empirical measure π MC = t −1 t−1 s=0 δ Xs . We can thus approximate π(h) and v(π, h) using the empirical mean and variance, denoted by π MC (h) and v MC (h), and these are typically consistent as t → ∞.
The difficulty is in the term π((h − π(h)) · g) in (3.1), since g cannot be evaluated exactly. We employ unbiased estimators G(x) in Definition 5 in place of evaluations g(x). This leads to the following estimator of v(P, h), which we call the ergodic Poisson asymptotic variance estimator (EPAVE), where each G(X s ) is conditionally independent of all others given (X 0 , . . . , X t−1 ). In Section 4.6, we show that this estimator is strongly consistent and satisfies a √ t-CLT under Assumption 1 and moment assumptions on h. The proposed EPAVE might be practically relevant, but in this paper we view it as an intermediate step toward the estimators of Section 3.2, which have the advantage of being unbiased.
Remark 9. We can compute (3.2) online by keeping track of the sums We can also modulate the relative cost of estimating the fishy function evaluations in (3.2) by generating G(X s ) only for some of the times, e.g. when s mod D = 0 for some D ∈ N, which amounts to a thinning strategy and can be beneficial (Owen 2017).

Subsampled unbiased asymptotic variance estimator
Starting again from (3.1), we propose an unbiased estimator of v(P, h) by combining unbiased estimators G(x) of g(x) with unbiased Glynn-Rhee approximations of π (Glynn & Rhee 2014, reviewed in Appendix A. We thus assume that we can generate random signed measureŝ for a class of test functions h, where the expectation is with respect to coupled lagged Markov chains as generated by Algorithm 2.1, with X 0 , Y 0 initialized from an arbitrary coupling of π 0 with itself. Combining these measuresπ with G(x) in (2. 5), each term in (3.1) can be estimated without bias, as we next describe. First, unbiased estimators of v(π, h) = V π [h(X)] can be obtained using two independent unbiased measuresπ (1) Unbiased estimation of the term π((h − π(h)) · g) in (3.1) is more involved. We first provide an informal reasoning that motivates the proposed estimator given below in (3.5) and described in pseudo-code in Algorithm 3.1. A formal justification, including bounds on the moments of the proposed estimator, is postponed to Section 4.5. Consider the problem of estimating π(h · g) without bias, and assume that we can generate unbiased measuresπ = N n=1 ω n δ Zn of π and estimators G(x) with expectation equal to g(x) for all x. Then we can generate are conditionally independent of one another given (Z n ). Conditioning onπ, we have and then taking the expectation with respect toπ yields π(h · g), under adequate assumptions on h · g.
However, the variable N n=1 ω n h(Z n )·G(Z n ) requires estimators of the fishy function for all N locations, and N could be large. Alternatively, after generatingπ we can sample an index ∈ {1, . . . , N } according to a Categorical distribution with strictly positive probabilities ξ = (ξ 1 , . . . , ξ N ), and given Z we can generate G(Z ). Then we observe that, conditioning onπ, integrating out the randomness in G(Z ) given Z , and then the randomness in , and therefore ω h(Z )·G(Z )/ξ is an unbiased estimator of π(h·g) that requires only one estimation of g at Z . The estimator proposed below employs R ≥ 1 estimators of the fishy function for each signed measureπ, where R is a tuning parameter, and the choice of selection probabilities ξ is discussed in Section 3.3.
We can finally introduce the proposed estimator, which empirical performance is illustrated in Section 5. We writeπ (j) as 5) and its generation is described in Algorithm 3.1. The total cost ofv(P, h) will typically be dominated by the cost of obtaining two signed measuresπ (j) for j ∈ {1, 2} and 2R estimators of evaluations of g.
We callv(P, h) the subsampled unbiased asymptotic variance estimator (SUAVE). In Section 4.5, we show that under Assumption 1,v(P,h) is unbiased and has p finite moments whenever h has sufficiently many moments. For simplicity the proof assumes that ξ (j) k = 1/N (j) for all k ∈ {1, . . . , N (j) }.

Algorithm 3.1 Subsampled unbiased asymptotic variance estimator (SUAVE).
Input: generator of unbiased signed measuresπ, generator of unbiased fishy function estimators G, method to compute selection probabilities ξ, integer R.

Implementation and improvements
The proposed estimator relies on unbiased signed measuresπ of π. In the notation of Appendix A, for all experiments below we generate lagged chains, record the meeting times, and choose L and k as large quantiles of the meeting times. Then we choose m as a multiple of k such as 5k, following the suggestions in Jacob et al. (2020).
To implement SUAVE we need to choose selection probabilities ξ = (ξ 1 , . . . , ξ N (j) ) given N (j) , the number of atoms in the signed measureπ (j) . We set these probabilities to 1/N (j) as a default choice.
We can also try to minimize the variance of the resulting estimators with respect to ξ. This requires information on the variance of G. Indeed, if we condition on the realizations ofπ (j) , j ∈ {1, 2}, then the variance of the term is minimized over ξ as follows. Since its expectation is independent of ξ, we can equivalently minimize its second moment, thus we define √ α n , n = 1, . . . , N (j) .
The use of ξ leads to a second moment equal to ( √ α n ) 2 , but for any ξ such that N (j) n=1 ξ n = 1, the Cauchy-Schwarz inequality implies Therefore ξ in (3.7) results in the smallest variance of the term in (3.6). It might sometimes be possible to approximate E[G(z) 2 ] as a function of z using pilot runs. We investigate this with numerical experiments in Section 5.2.
We also need to choose R, where 2R is the number of states in X at which the fishy function g is estimated. Throughout Section 5 we observe gains in efficiency when using R greater than 1. We can guide the choice of R numerically by minimizing inefficiency defined as the product of expected cost and variance, which can be approximated using empirical averages of M independent copies (v (m) (P, h)) M m=1 . We note that when we run SUAVE with a given choice of R, we can easily output estimators corresponding to any smaller value of R at no extra cost, which is helpful for the assessment of the impact of R.
A naive implementation of SUAVE with Algorithm 3.1 could incur a large memory cost when each state in X is large, as in high-dimensional regression (Biswas et al. 2022) or phylogenetic inference (Kelly et al. 2021). Indeed storing all the atoms of the generated signed measures might be cumbersome.
However, for SUAVE we only need to select within each measure R atoms at which to evaluate h and to estimate g; see Line 3c in Algorithm 3.1. We can address the memory issue by setting ξ for all k and by using reservoir sampling (Vitter 1985). This technique allows to sample (r,j) uniformly in {1, . . . , N (j) }, R times independently, without knowing N (j) in advance and keeping only R objects in memory.
We mention methodological variations that we do not investigate further in this manuscript.
1. Instead of sampling R atoms from each signed measure with replacement, we can sample without replacement.
2. The 2R estimators G employed in (3.5) could be generated jointly instead of independently.
In particular, we can couple 2R chains starting from (Z (r,j) ) R r=1 , j ∈ {1, 2}, and a single chain starting from y: in other words, we can simulate a coupling of 2R + 1 chains instead of simulating a coupling of two chains 2R times.
3. As discussed around (2.6) we can either choose y ∈ X in the definition of g and G, and keep y constant throughout the simulations as we do in Section 5, or we can draw y from some distribution ν before each fishy function estimation.

Assumption on the meeting time
We first prove Proposition 2, which can be used to verify that Assumption 1 holds.
Proof of Proposition 2. The first part follows by Markov's inequality: For the second part, we have P Taylor's theorem gives and so The following sufficient condition for g ∈ L p (π) is used to ensure that g is well-defined under Assumption 1 when h has sufficiently many moments.
Proposition 11 (Douc et al. 2018, Prop. 21.2.3). Let P be a Markov kernel with unique invariant We consider here what Assumption 1 implies about the corresponding P and its fishy functions.
First, we observe that this assumption implies that P is aperiodic, and also ergodic of degree 2 (as defined e.g. in Nummelin 1984, Section 6.4), which implies e.g. that a CLT holds for ergodic averages of bounded functions.

Proposition 12. If Assumption 1 holds then
and P is aperiodic and ergodic of degree 2.
We now turn to the implication of Assumption 1 on properties of g . In particular, we see that κ > p with sufficiently many moments of h implies that g ∈ L p (π).
As is well known, it is not necessary that g ∈ L 2 0 (π) for a CLT to hold. Theorem 14 states that the CLT holds under Assumption 1 and that the asymptotic variance v(P, h) is given by the expression that we approximate in Section 3. The proof in Appendix D relies heavily on the strategy of Douc et al.
Remark 15. The results above rely only on the existence of meeting times with polynomial survival functions, so one can deduce that they hold for any Markov kernel P such that P t (x, ·) − π TV ≤ M (x)(t + 1) −s with s > κ and π(M ) < ∞, since then Conversely, for such a Markov kernel there exists a possibly non-Markovian coupling that would satisfy Assumption 1 (see, e.g., Griffeath 1975), but we do not pursue this here.

Unbiased approximation of fishy functions
The following technical lemma will be useful to obtain bounds on the moments of the fishy function estimator in Definition 5.
The following lemma employs dominated convergence to justify the interchange of expectation of infinite sum and thereby ensure that G y (x) is indeed an unbiased estimator of g y (x).

Lemma 17. Under Assumption 1, let
Since m > κ/(κ − 1), Theorem 13 provides that g ∈ L 1 0 (π), and hence for π-almost all x and y we have and so the assumed integrability of the right-hand side implies, by dominated convergence, that The following lemma is used several times to ensure that expectations of functions of X t are uniformly bounded in t under reasonable conditions. In the first case, the conclusion is a result of stability properties of π-invariant Markov chains at almost all points, while in the second case regularity is imposed by ensuring that µ cannot place too much probability in possibly problematic regions; for example it guarantees that φ ∈ L 1 (µ).
Proof. By assumption, P is π-irreducible and from Proposition 12 it is aperiodic. Since π(φ) < ∞, we which concludes the proof.
The following theorem provides sufficient conditions for the estimator G y (x) to be unbiased and have finite pth moments.

Unbiased approximation of π(h)
We next demonstrate that the approximation in Definition 7 is indeed unbiased and has finite pth moments under suitable conditions. The proof of Theorem 24 is an application of Theorem 21 and its statement can be compared with Middleton et al. (2020, Theorem 1), which treats the case p = 2, for which we essentially arrive at the same condition for κ and m. The lack-of-bias condition here is less demanding, and we deduce finiteness of higher moments of H for m sufficiently large.
In Propositions 28-29 we establish that the properties obtained for H may be deduced also for averages of lagged and offset estimators that are used in practice following Jacob et al. (2020). We use such unbiased estimators as part of the proposed asymptotic variance estimator SUAVE of Section 3.2. Finally, in Proposition 30 we show that subsampled estimators are also unbiased and have finite pth moments under the same conditions, because subsampling is an important aspect of SUAVE that controls its computational cost.
The following lemma guarantees that if an "independent initialization" is used, and dµ/dπ ≤ M , then E γ [τ κ ] < ∞ is guaranteed by Assumption 1, and similarly lack-of-bias results for γ-almost all (x, y) may be deduced from lack-of-bias results for π ⊗ π-almost all (x, y) from which we may deduce that dµP/dπ ≤ M . It then follows that for A, B ∈ X , from which we may conclude.
Remark 25. By Theorem 24 it is sufficient that κ > 2 and m > 2κ/(κ − 2) for H to have a finite variance. On the other hand, Theorem 14 implies that a CLT holds for h if κ > 1 and m > 2κ/(κ − 1), which is weaker. In particular, the stronger condition for Theorem 24 is because finite second moment of the unbiased estimator is shown via finite second moment of the approximation of g , and this requires that g ∈ L 2 0 (π). For variance reduction purposes, it might be preferable not to use H but instead an unbiased estimator with k-offset and a lag of L between the chains, as described in Appendix A. We show here that the previous results continue to hold for these estimators.
Definition 26. The L-lagged and k-offset approximation of π(h) is is a time-homogeneous Markov chain with Markov kernelP , and t=0 is a Markov chain with transition kernel P and initial distribution π 0 and Y 0 ∼ π 0 independent of (X 0:L ). In particular, ( Remark 27. The approximation may be viewed as the sum of h(X k ) and an unbiased approximation If g (L) ∈ L p (π) then, noting that , since P is a bounded linear operator in L p (π).
The following shows that Theorem 24 holds for general L ≥ 1 and k ≥ 0, and in fact increasing either of these decreases the upper bound on the moments of the estimator.

Proposition 28. Under Assumption 1, let
Proof. Recall that a ∨ b stands for the maximum and a ∧ b for the minimum of a and b. SinceP L is a coupling of P L with itself and P L is π-invariant, we seek to apply Theorem 24 to the approximation where the second term is obtained by considering a chain (X (L,k) , Y (L,k) ) with Markov transition = γ (L)P k = (π 0 P L ⊗ π 0 )P k . This is analogous to the coupled chain described in Section 2.1 which is used in Definition 5. We define and we seek to verify that Assumption 1 holds for τ L,k and obtain a finite bound for E to the law of (X (L,k) , Y (L,k) ). In Theorem 24, we take µ = π 0 P k . If we define (X, Y ) with Markov transition kernelP and (X 0 , Y 0 ) ∼ γ (L) , and define τ = inf{t ≥ 0 : X t = Y t }, we observe that (X (L,k) , Y (L,k) ) may be taken as a skeleton of this chain, i.e. (X Therefore, we may deduce that It follows from (4.5) that for p ≥ 1 such that 1 The following shows that the unbiased signed measure used in Section 3 and described in Appendix A is indeed unbiased for functions with suitably large moments when κ is large enough, and that moments of averaged unbiased estimators are finite under the same conditions as for H. ) k: ] = π(h) and for p ≥ 1 such that from which we may conclude using the fact that the upper bound in (4.6) is non-increasing in k.
We now demonstrate that estimators associated with subsampling the unbiased signed measureπ described in Appendix A are unbiased and have finite pth moments under the same conditions as the standard estimator.

Unbiased approximation of π(h 0 · g)
We now look at combinations of unbiased fishy function estimation and unbiased estimation of π(h).
This will involve estimators G y (x) at several points x drawn randomly. We first define an alternative representation that will avoid some ambiguity in the following developments. We define a probability measure Q such that, with U ∼ Q, and we define an extended distributionπ(dx, du) = π(dx)Q(du) withπ-invariant Markov kernel T (x, u; dy, dv) = P (x, dy)Q(dv), and coupled Markov kernel We useǦ f to indicate thatǦ f y,v (x, u) is an unbiased approximation of the fishy function for T and f as opposed to P and h, i.e. such that is the fishy function for T and f in L 1 0 (π).
The Markov chain (X, U, Y, V ) has the same meeting time as the Markov chain (X, Y ) by construction, so Assumption 1 holds for this chain with π replaced byπ. Similarly if dµ/dπ ≤ M then withμ = µ ⊗ Q andπ = π ⊗ Q we have dμ/dπ ≤ M . Hence, we may apply Theorem 24 to deduce that an unbiased estimator ofπ(φ) satisfies E[Φ] =π(φ) and that E [|Φ| p ] < ∞ if φ ∈ L m (π) and 1 ≥ 1 p > 1 κ + 1 m . The following two lemmas provide conditions for g y ∈ L q (π) and h · g y ∈ L s (π), which are used to analyze both the MCMC estimator of v(P, h) of Section 3.1 and the unbiased estimators of Section 3.2.
Definition 33 (Basic estimator of π(h 0 · g y )). Let H = H (L) k: be an unbiased estimator of π(h) as defined in Proposition 29, with dπ 0 /dπ ≤ M . Then, with random variables independent to those used to define H, let The choice of γ in Definition 33 ensures that E γ [τ κ ] < ∞ under Assumption 1. The following proposition shows that Φ is indeed an unbiased estimator of π(h 0 · g y ), and that it has a pth moment if κ and m are sufficiently large.

For
Proof. To apply Theorem 24 we determine s ≥ 1 such that φ H ∈ L s (π). Let φ 1 (x, u) = h(x)g y (x, u) and φ 2,c (x, u) = cg y (x, u), so that φ H = φ 1 − φ 2,H . We find that for s ≥ 1, We may then deduce that for which it remains to show that the final expectation on the right-hand side is finite. We have and since E[H] = π(h), we may conclude.

Theorem 35. Under Assumption 1 with
be an unbiased estimator of π(h) as defined in Proposition 29, with dπ 0 /dπ ≤ M . Let u), and S R be as in Proposition 30 with the unbiased signed measure approximatingπ involving random variables independent to those used to define H. Then for π-almost all y: 1. E [S R ] = π(h 0 · g y ).

For
Proof. To apply Proposition 30 we determine s ≥ 1 such that φ H ∈ L s (π). By the same argument as in the proof of Proposition 34, we find that φ H ∈ L s (π) for 1 s > We may deduce that E φ H p L s (π) < ∞ from the same argument as in the proof of Proposition 34. Hence, we may conclude that for 1 For the lack-of-bias property, we consider p = 1 and, if m > 2κ/(κ − 2) then, by Proposition 29, and since E[H] = π(h), we may conclude.

Subsampled unbiased asymptotic variance estimator
We show that the basic and subsampled unbiased asymptotic variance estimators are indeed unbiased, and have finite pth moments under the same conditions. Definition 36. Let H 1 and H 2 be two independent unbiased estimators of π(h) and S be a unbiased estimator of π(h 2 ), all of the type described in Proposition 29. Let Φ be as in Definition 33 and define or similarly for any integer R ≥ 1 let S R be as in Theorem 35 and define Lemma 37. Under Assumption 1, let h ∈ L m (π) for some m > 2κ/(κ − 1). Let H 1 , H 2 and S be as in Definition 36. Then

For
Proof. We note that if h ∈ L m (π) then h 2 ∈ L m/2 (π). Hence, we deduce by Proposition 29 that if By Minkowski's inequality, for p ≥ 1: Remark 39. By Minkowski's inequality we may similarly conclude that any average of estimators of the form given in Definition 36 also has lack-of-bias and moments implied by Theorem 38, and hence that the results apply to (3.5) when ξ . . . , N (j) }. Remark 40. If Assumption 1 holds for all κ ≥ 1 (respectively m ≥ 1), one requires only slightly more than 2 moments of h (respectively τ ) to estimate the asymptotic variance consistently and slightly more than 4 moments of h (respectively τ ) to approximate the asymptotic variance with a variance in if M is the number of independent unbiased estimators averaged. This seems close to tight, since 4 moments of h are required for the sample variance to have a finite variance in the simpler setting of independent and identically distributed random variables.

Ergodic Poisson asymptotic variance estimator
The asymptotic variance estimator in (3.2) can be analyzed using somewhat standard convergence theorems. For the CLT in particular, it is helpful to view the Markov chain on the extended space introduced in Section 4.4.

Proposition 41. Under Assumption 1, let X be a Markov chain with Markov kernel P , and h ∈
L m (π) with m > 2κ/(κ − 1). For π-almost all X 0 , the CLT holds for h and for π-almost all y, v(P, h) = −v(π, h) + 2π(h 0 · g y ). The estimator (3.2) and these terms converge almost surely to π(h 2 ) and π(h) 2 , respectively, as t → ∞ by the Markov chain law of large numbers (see, e.g., Douc et al. 2018, Theorem 5.2.9) and continuous mapping, and The assumptions guarantee by Theorem 14 that for π-almost all y, h·g y ∈ L 1 (π) and since E [G y (x)] = g y (x) for π-almost all x by Theorem 21, the first term on the right-hand side converges almost surely to π(h · g y ) by the Markov chain law of large numbers, while similarly for the second term we have π MC (h) → a.s. π(h) and 1 t t−1 s=0 G y (X s ) → a.s. π(g y ), so the second term converges almost surely to π(h)π(g y ). Hence, the left-hand side converges almost surely to π(h 0 · g y ), and we conclude. 2 , f 3 = f 1 · g y and f 4 = g y . We observe thatv Ui) is the empirical measure of the Markov chain introduced in Section 4.4. Hence, we define

Numerical experiments
Through numerical experiments we investigate some distinctive features of the proposed estimators.
The initial distribution is π 0 = Normal(0, 4 2 ). The target distribution is Normal(0, (1 − φ 2 ) −1 ), and We use a reflection-maximal coupling to define the coupled kernelP , see Appendix B. Using coupling-based diagnostics (Biswas et al. 2019) we find evidence that the chain converges to stationarity in less than 500 iterations, and we choose k = 500, L = 250, m = 2500 for unbiased MCMC approximations. The state y ∈ X = R, used to define g = g y in (2.4), is set as y = 0.
We start by estimating the fishy function g(x) using G in ( (1) (1) example.
stat_smooth function of ggplot2 (Wickham 2016) that fits a generalized additive model using the mgcv package (Wood 2017) and shows the fitted values in full line and a sequence of 95% confidence intervals in a gray ribbon.
The performance of the proposed estimator of v(P, h) is shown in Table 5 shows the variance of eachv(P, h), which decreases with R but with diminishing returns, and the last column reports the variance multiplied by the average cost, which is known as the inefficiency (Glynn & Whitt 1992). From the last column, increasing R improves the efficiency with diminishing returns.
Overall we find that we can get accurate estimates of v(P, h) with parallel runs that each costs about 5000 − 20000 units, depending on the choice of R, and we emphasize the strong efficiency gains in increasing R from 1 to 10 or more.
We compare the proposed method to two families of estimators: "batch means" and "spectral variance" estimators from the mcmcse package (Flegal et al. 2020). These estimators are based on batch sizes determined with the method of Liu et al. (2022). For each class of methods, we try three values of the lugsail parameter, r ∈ {1, 2, 3}. We compute these estimators of v(P, h) using chains of lengths in {10 4 , 10 5 , 10 6 }, and different numbers of parallel chains in {4, 8}. We base all the results shown here on 400 independent trajectories, thus we obtain 100 independent estimators each based on four chains, and 50 independent estimators each based on eight chains. For each method, we report the estimate, the bias, the mean squared error (MSE), and the total cost equal to number of chains multiplied by the time horizon. Finally, we report the MSE that would be achieved by our proposed Table 5.2: Batch means estimation of the asymptotic variance v(P, h), here equal to 10 4 , in the AR (1) example.
method, using R = 10 (second row of Table 5.1), when generating a number of independent replicates such that the total costs are matched; in Tables 5.2 and 5.3 we refer to this as "MSE proposed", and it can be compared to "MSE". Table 5.2 presents the results for batch means estimators. We see that r = 1 leads to negativelybiased estimators, while r = 3 leads to positively biased estimators, as anticipated. From the comparison of the MSEs in the last two columns of Table 5.2, we see that the proposed estimator is at least competitive with batch means estimators in this example.
We repeat the same experiments with spectral variance estimation, with the same high-level conclusions; see Table 5.3. We recall that batch means and spectral variance estimators do not in general converge at the Monte Carlo √ t rate, where t is the length of the trajectory; see Chakraborty et al. (2019) for batch means estimators. For example, in the tables we see that using 8 chains instead of 4 systematically doubles the cost but does not always halve the MSE. Likewise, multiplying the time horizon by 10 incurs a deterministic ten-fold increase in cost but the MSE is not systematically divided by 10. This contrasts with the proposed method, where averaging c ≥ 1 times more independent estimators will precisely decrease the variance by c, while the expected cost will be essentially multiplied by c, with the caveat that, when running stochastic simulations on parallel processors, the user might need to wait for the slowest runs to complete in order to avoid biases; see Section 3.3 in Jacob et al. (2020).

Comparing two MCMC algorithms
Here we consider two MCMC algorithms targeting the same distribution. We wish to determine which algorithm leads to the smallest asymptotic variance. We also illustrate the effect of the selection probabilities ξ following the considerations of Section 3.3.
The target is the posterior distribution in a Cauchy location model, with a Normal prior on the location parameter denoted by θ. We observe n = 3 measurements z 1 = −8, z 2 = +8, z 3 = +17, Table 5.3: Spectral variance estimation of the asymptotic variance v(P, h), here equal to 10 4 , in the AR(1) example.
assumed to be realizations of Cauchy(θ, 1). The prior on θ is Normal(0, 100). We consider the test function h : x → x. The first algorithm considered here is the Gibbs sampler described in Robert (1995), alternating between Exponential and Normal draws: where σ 2 = 100 is the prior variance. The coupling of this Gibbs sampler is done using common random numbers for the η-variables, and a maximal coupling for the update of θ; see Appendix B. The second algorithm is a Metropolis-Rosenbluth-Teller-Hastings (MRTH) algorithm with Normal proposal, with standard deviation 10. Its coupling employs a reflection-maximal coupling of the proposals. The state y used to define g = g y in (2.4) is set to zero. The initial distribution, for both chains, is set to π 0 = Normal(0, 1). Using coupling-based diagnostics, both MCMC algorithms appear to converge in less than 100 iterations and we set k = L = 100 and m = 500.   Table 5.4: Estimators of v(P, h) for two MCMC algorithms targeting the posterior distribution in the Cauchy-Normal example, using either optimal or uniform selection probabilities ξ.
probabilities. The results are shown in Table 5. 4. We see that the use of optimal ξ provides a sizable reduction of variance. We note that the estimates, along with the reported standard errors, are precise enough for us to rank the two asymptotic variances with confidence.

High-dimensional Bayesian linear regression
We consider a more challenging example, with a linear regression of n = 71 responses on p = 4088 predictors of the riboflavin data set (Bühlmann et al. 2014). The model, MCMC algorithm and its coupling are taken from Biswas et al. (2022) and a self-contained description is provided in Appendix C; essentially we use a shrinkage prior on the coefficients (e.g. Bhadra et al. 2019) and the target distribution is the posterior distribution of the regression coefficients along with the global precision, the local precisions, and the variance of the observation noise. The state space is of dimension 2p+2 = 8178.
This example motivates the implementation of reservoir sampling to select the states at which to estimate the fishy function, as described in Section 3.3.
Based on preliminary runs, we choose the test function h : x → β 2564 , which is a coordinate of the regression coefficients with a clearly bimodal marginal posterior distribution. Figure 5.3a shows three independent traces of h(X t ) over the first 1000 iterations. Figure 5.3b presents a histogram of β 2564 , obtained from 10 independent chains run for 50,000 iterations each and discarding 1000 iterations as burn-in. Figure 5.3c shows upper bounds on |π t − π| TV obtained with the method of Biswas et al.
(2019), using a lag L = 1000 and 10 3 independent meeting times. From this we choose k = L = 1000 and m = 5k.
To define g = g y , we draw y once from the prior, and keep it fixed. We generate M = 10 3 independent estimates of v(P, h), for R ∈ {1, 2, 5, 10}. The results are summarized in Table 5.5. The estimator v (P, h) has a relatively small standard error, and we observe tangible gains in efficiency when increasing R.    , which is about 2.5 times larger than the inefficiency of MCMC, found to be around 82 according to Table 5.5. As argued in Jacob et al. (2020), that relative inefficiency could be reduced by increasing k and m. Thanks to the proposed SUAVE, we can anticipate that the inefficiency could be reduced by a factor of 2.5 with increased k and m, without running any longer chains.

Discussion
When a Markov chain admits an accessible atom α, a solution of the Poisson equation (e.g. Glynn & Meyn 1996) , where σ α = inf{n ≥ 0 : X n ∈ α}. This allows approximation of g pointwise by simulation if one can identify entries into α and one can approximate h 0 = h − π(h) pointwise. The approach we took in this article is different in that it does not rely on identifying atoms, and hence may be applicable in scenarios where this regenerative approximation is not implementable.
For example, proper atoms may not exist for a given general state space Markov chain, and identification of hitting times of an atom for a suitable split chain as in Mykland et al. (1995) is not always feasible. While it is often possible to define a modified Markov chain that admits an easily identified, accessible atom (Brockwell & Kadane 2005, Lee et al. 2014, the corresponding solution of the Poisson equation may not be similar to that of the original chain. Note also that, when atoms can be identified, one would often use regenerative simulation to approximate the asymptotic variance (Hobert et al. 2002), which can be expressed as v (P, h) 2 (see, e.g., Bednorz et al. 2008).
Coupling techniques have been extensively used to study distances between a Markov chain at time t and its limiting distribution. Couplings yield upper bounds on these distances and on the resulting mixing time, and these bounds are rarely sharp. Therefore, without matching lower bounds, coupling techniques cannot determine whether an algorithm converges faster than another. On the other hand, the proposed coupling-based asymptotic variance estimators are unbiased, even with sub-optimal couplings, under conditions such as those presented in Section 4. Thus, with sufficient computing resources we can determine which algorithm leads to the smallest asymptotic variance.
Our experiments suggest that the proposed estimators of v(P, h) are practical and competitive.
They only apply to settings where the requisite couplings are available, so that the proposed methods demand more from the user, compared to batch means and spectral methods. On the other hand, the proposed estimators are unbiased, with the important consequence that averages of independent copies achieves the Monte Carlo rate of convergence. This compares favourably to batch means estimators, for which, under general assumptions, convergence occurs at rate √ a n , where a n is the number of batches, and a n is often chosen as a n = n 1/2 or a n = n 2/3 where n is the length of the chain (see Instead of the above expressions that involve a test function h, we can consider the signed measurê 4) as an unbiased approximation of π.
If we count the cost of sampling from the kernel P as one unit, and the cost of sampling fromP as one unit if the chains have met already, and two units otherwise, then the random cost of obtaining (A. 2), or (A.4), equals max(L, m + L − τ (L) ) + 2(τ (L) − L) units.
The above unbiased estimators are to be generated independently in parallel, and averaged to obtain a final approximation of π. Since they are unbiased, the mean squared error is equal to the variance. To compare unbiased estimators with different cost, e.g. to compare different configurations of the tuning parameters k, m, L, we can compute the asymptotic inefficiency defined as the expected cost multiplied by the variance, as described in Glynn & Whitt (1992); the lower value, the better.
A by-product of the above estimator is an upper bound on |π t − π| TV (Biswas et al. 2019), given The upper bound can be estimated using independent replications of the meeting time τ (L) . We can use this bound to choose "burn-in" values, by choosing t such that the estimated upper bound is small.

B Coupling of MCMC algorithms
For completeness we describe how Metropolis-Rosenbluth-Teller-Hastings (MRTH) algorithms (Metropolis et al. 1953, Hastings 1970 can be coupled in such a way that exact meetings occur. Various relevant considerations are presented in Wang et al. (2021).
Algorithm B.1 describes a simple way of coupling P (X, ·) and P (Y, ·) where P is the transition associated with MRTH, with proposal transition q and acceptance rate α MRTH . A maximal coupling of two distributions µ(dx) and ν(dy) refers to a joint distribution γ(dx, dy) with marginals µ and ν, and such that P(X = Y ) is maximized over all such joint distributions. Algorithm B.1 describes a Markov kernelP ((X, Y ), ·), that can be employed to generate a pair of chains (X t ) and (Y t ) such that the chains meet at a random time.

Algorithm B.1 A coupled Metropolis-Rosenbluth-Teller-Hastings kernel
We provide here a short and self-contained description of one particular version of the Gibbs sampler and its coupling; many more algorithmic considerations can be found in Biswas et al. (2022).

C.1 Model
The context is that of linear regression, with n individuals and p covariates, with p n. The generative model is described below, where Y is the outcome, X the vector of explanatory variables, β ∈ R p the regression coefficients, σ 2 ∈ R + the observation noise, ξ is called the global precision and η j is the local precision associated with β j for j ∈ {1, . . . , p}, The distribution t(ν) + refers to the Student t-distribution with ν degrees of freedom, truncated on (0, ∞), with density x → (1 + x 2 /ν) −(ν+1)/2 1(x ∈ (0, ∞)) up to a multiplicative constant. The hyperparameters are set as a 0 = 1, b 0 = 1, ν = 2. In our experiments we initialize Markov chains from the prior distribution.

C.2 Gibbs sampler
The main steps of the Gibbs sampler under consideration are as follows.
Overall the computational complexity is of the order of n 2 p operations per iteration, therefore it can be used with large p and moderate values of n. Details on each step can be found below.
η-update. The conditional distribution of η given the rest has density which we can target with the slice sampler described in Algorithm C.1, applied independently componentwise.

C.3 Coupled Gibbs sampler
We consider only one of the variants in Biswas et al. (2022), which is not necessarily the most efficient but achieves good performance in the experiments of Section 5.3 and is simpler than the "two-scale" coupling described in Biswas et al. (2022). We describe how to couple each update, with the first chain in state η, ξ, σ 2 , β and the second in stateη,ξ,σ 2 ,β.
η-update. We consider two strategies to couple the slice sampling updates of η j , as described in Algorithm C.1.
1. We can use a common uniform V in the first step of Algorithm C.1, to define U j for the first chain andŨ j for the second. Then we can sample from a maximal coupling of the distributions of η j |U j andη j |Ũ j , using Algorithm B.2. This strategy results in a non-zero probability for the event {η j =η j }.
2. We can use a common uniform V in the first step of Algorithm C.1, and then a common uniform U * in the third step. This is a pure "common random numbers" (CRN) strategy.
We adopt a "switch-to-CRN" strategy: we scan the components j ∈ {1, . . . , p}, and sample η j ,η j using the maximal coupling strategy above. If any component fails to meet, we switch to the CRN strategy for the remaining components.
ξ-update. To update ξ,ξ, we draw the proposals in the MRTH step using a maximal coupling as in Algorithm B.2. We then employ a common uniform variable for the two acceptance steps. β-update. We use a CRN strategy, which amounts to using the same draws r, δ in the first step of Algorithm C.2 to sample both β andβ.

D Proof of Theorem 14
Lemma 43 ([Rio, 1993, Theorem 1.1]). Let X and Y be integrable random variables such that where the supremum is over all measurable sets. Then where for a random variable Z, Q Z is the tail quantile function of |Z|, i.e. Q Z (u) = inf{t : P(|Z| > t) ≤ u}.
Lemma 44. Assume that, for all k ∈ N, Then for g, h measurable functions, and A, B measurable sets, Lemma 45. Assume that Then with h such that π(h) = 0, and k, ∈ N, where Q 0 is the tail quantile function of |h(X 0 )|.
Proof. We follow the same strategy as in Douc et al. (2018, Lemma 21.4.3). By Lemma 44 with g = P h, we obtain for use in Lemma 43. It follows that where Q Z is the tail quantile function of |Z|. Since h(X k ) has the same distribution as h(X 0 ), Q h(X k ) = Q 0 . On the other hand by Douc et al. (2018, Lemma 21 Hence, by Cauchy-Schwarz, we have By interchanging the use of k and , we obtain the final bound. The following Lemma is similar to Douc et al. (2018, Theorem 21.4.4).
Lemma 46. Assume that and let h ∈ L m 0 (π). Then Proof. We have by Markov's inequality from which we may deduce that Q 0 (u) ≤ h L m (π) u −1/m and so By Lemma 45, it follows that from which we may conclude.

E Verifying the assumptions in the AR(1) case
Assumption 1 may be verified abstractly using polynomial drift and minorization via the argument in Middleton et al. (2020, Theorem 2) by strengthening their Assumption 5 to include π(V ) < ∞. We demonstrate here how one can verify Assumption 1 directly in the AR(1) example, for which we find that Assumption 1 holds with any κ > 1 since the survival function decays geometrically.
Lemma 47. For Algorithm E.1, the meeting time τ = inf{n ≥ 1 : X n = Y n } satisfies where D = (D n ) is an AR ( Proof. We consider the following equivalent construction of X and Y as in Algorithm E.1. Let (W n ) be a sequence of independent Normal(0, 1) random variables. Set X 0 = x 0 and define X n = φX n−1 + σW n , n ≥ 1, and then let Y 0 = y 0 and define Y n = φY n−1 − σW n , n ≥ 1.
We can then define (B n ) to be a sequence of conditionally independent Bernoulli random variables with B n ∼ Bernoulli (1 ∧ g(X n−1 , Y n−1 , W n )) , where g(x, y, w) = ϕ φ (x−y) σ + w /ϕ(w). X is the Markov chain described by the algorithm, while Y is the Markov chain associated with the Y chain when the Bernoulli random variables B n take the value 0 for all n. In particular, if τ = inf{n ≥ 1 : B n = 1) then Y n = Y n for n ∈ {0, . . . , τ − 1} and Y n = X n for n ≥ τ . Hence, the meeting time τ can be determined by analyzing the X and Y chains, and in particular P x0,y0 (τ > n) = E x0,y0 Now define D n = (X n − Y n )/(2σ). We notice that this is an AR(1, φ, 1) chain with D n = φD n−1 + W n , and that we may re-express the conditional distribution of B n as It follows that, with d 0 = (x 0 − y 0 )/(2σ), This equivalence suggests the relevance of that expectation with respect to an AR(1, φ, 1) chain, which we bound below. The proof is delayed to the end of this appendix as it requires several intermediate results.
Proposition 48. Let X be AR(1, φ, 1) and P its corresponding Markov kernel. Then Combining Lemma 47 with Proposition 48, we obtain the following.
We see that the dependence on σ and |x − y| is fairly mild. On the other hand, if one calculates the dependence ofβ on φ, one finds that it deteriorates quickly as φ 1, even though it remains less than 1. In contrast, numerical experiments suggest that the true geometric rate is in fact φ, but we are not aware of a proof technique that is able to capture such a rate. Indeed, the calculations we have used to provide a rigorous bound are similar to those used to provide quantitative convergence rates for Markov chains more generally and these are often loose in practice.
Lemma 50. Let X be a Markov chain with Markov kernel P . Assume there exists V ≥ 1, (β, b) ∈ (0, 1) × [1, ∞) such that for some set C ⊂ X, where C is the complement X \ C. Then for G : X × X → [0, 1], where we may take δ = log h/(log h + log β − log b) ∈ (0, 1) for any (0, 1) Proof. By Hölder's inequality and the assumptions we have for any δ ∈ (0, 1), where the equality is due to the specific choice of δ. Now, since 0 ≤ G ≤ 1 and V ≥ 1, we have It follows that and hence that A n ≤ β nδ V (x) δ .
In combination with Lemma 47, we obtain the following.
In the above, the dependence ofC on (x 0 − y 0 ) 2 is suboptimal, and we can improve this via the the following result.