# A non-asymptotic error analysis for parallel Monte Carlo estimation from
many short Markov chains

Austin Brown ,
Department of Statistical Sciences, University of Toronto,
<EMAIL_ADDRESS>0000-0003-1576-8381

###### Abstract

Single-chain Markov chain Monte Carlo simulates realizations from a Markov
chain to estimate expectations with the empirical average. The single-chain
simulation is generally of considerable length and restricts many advantages
of modern parallel computation. This paper constructs a novel many-short-
chains Monte Carlo (MSC) estimator by averaging over multiple independent sums
from Markov chains of a guaranteed short length. The computational advantage
is the independent Markov chain simulations can be fast and may be run in
parallel. The MSC estimator requires an importance sampling proposal and a
drift condition on the Markov chain without requiring convergence analysis on
the Markov chain. A non-asymptotic error analysis is developed for the MSC
estimator under both geometric and multiplicative drift conditions. Empirical
performance is illustrated on an autoregressive process and the Pólya-Gamma
Gibbs sampler for Bayesian logistic regression to predict cardiovascular
disease.

Keywords: concentration inequalities; parallel Gibbs sampling; parallel Markov
chain estimation;

MSC: 60J27; 60J20

##  1 Introduction

Markov chain Monte Carlo (MCMC) is a widely applicable simulation technique to
estimate a integrals with respect to a complex target probability distribution
$\Pi$. Accurate estimation of integrals is central to modern Bayesian
inference where the target probability distribution may not be simulated from
directly. In the single-chain regime, MCMC simulates realizations from one
Markov chain $X_{0},\ldots,X_{t}$ until the marginal distribution is
approximately from $\Pi$. After discarding realizations before this ”burn-in”
time, the Markov chain is simulated further and the empirical average is used
to approximate integrals.

A persistent drawback of the single-chain regime is a necessarily long
simulation length required to reduce the bias in the estimation. It remains an
active research area to develop explicit and practically useful convergence
analyses on the mixing times of Markov chains. This computational drawback has
sparked research and debate into the many-short-chains regime by instead
combining simulations from multiple independent Markov chains of substantially
shorter length [Gelman and Rubin, 1992, Geyer, 1992, Margossian et al., 2023].
An advantage of the many-short-chains regime is modern parallel computation
allows simulating thousands, even millions, of independent Markov chains at
the same time. In applications, the many-short-chains approach appears to
empirically perform well for estimation in certain scenarios, but the
shortened simulation length of the Markov chains may theoretically result in a
large or even unknown bias.

Our main contribution constructs a novel many-short-chains Monte Carlo (MSC)
estimator and develops a theoretical non-asymptotic error analysis for
estimation. The MSC estimator is the empirical average over multiple sums from
independent Markov chains until their first return time to a particular set. A
significantly unique property of the MSC estimator is that it does not require
a convergence analysis on the Markov chain. To the best of our knowledge, the
non-asymptotic analysis developed is the first for a general purpose parallel
MCMC estimator that does not require a mixing time analysis of the Markov
chain. From a solely theoretical viewpoint, the MSC estimator can provide a
new perspective into the empirical performance of MCMC estimation in the many-
short-chains regime. The estimator may also be useful in non-trivial
applications as it may be simulated with general Markov chains such as Gibbs
samplers and is described in detail in Section 2 and applied in Section 6.

The non-asymptotic estimation error analysis ensures reliable performance of
the MSC estimator but requires a proper and careful construction. The
estimator is built upon essential theory in Markov chains for the
representation of invariant measures [Nummelin and Arjas, 1976, lemma 2] and
[Meyn and Tweedie, 2009, Theorem 10.4.9]. The independently generated Markov
chains require a random initial distribution combined with a geometric or
multiplicative drift condition. The random initial distribution uses an
importance sampling proposal for $\Pi$ that is constructed without the need to
compute the normalizing constant of $\Pi$. The length of the Markov chain
simulation is determined by the time to return to a set $C$ and controlled by
the drift condition. Since the drift condition ensures the return time to the
set $C$ is sufficiently fast, then the MSC estimator is guaranteed to remain
in the many-short-chains regime. Section 3 develops a non-asymptotic upper
bound on the mean squared error under a geometric drift condition and when the
importance weights have finite variance. Section 4 develops a concentration
inequality under a stronger multiplicative drift condition and uniformly
bounded importance weights. In Section 5 and Section 6, practical applications
are illustrated on a simple autoregressive process and a non-trivial example
with the Pólya-Gamma Gibbs sampler to predict cardiovascular disease.

The novel MSC estimator contributes to the recent research in parallel
implementations of Markov chain Monte Carlo ([Craiu and Meng, 2005, Craiu et
al., 2009, Lao et al., 2020, Jacob et al., 2020a] among others). Unbiased
Markov chain Monte Carlo is a recent approach to parallel estimation which
constructs an unbiased estimator from multiple Markov chains [Glynn and Rhee,
2014, Jacob et al., 2020b]. This requires constructing and simulating from a
joint Markov kernel between a Markov chain and lagged version of itself until
the first stopping time when the two chains exactly meet. Constructing optimal
joint Markov kernels is an active area of research [Wang et al., 2021] and
obtaining explicit convergence bounds on the time when two Markov chains meet,
even in moderate dimensions, can be challenging. In particular, the MSC
estimator may be beneficial in certain applications where Unbiased Markov
chain Monte Carlo may not readily accessible. For example, the length of the
Markov chain simulations in the MSC estimator depend only on the drift
condition which can be simpler to establish and scale efficiently to high
dimensions. However, this may not always be the case as explicit drift
conditions are often available for Markov chains such as Gibbs samplers but
can also be rare for others such as Metropolis-Hastings. Section 7 discusses
benefits and drawbacks of our results and future research possibilities.

##  2 The MSC estimator

Let X be a Borel measurable metric space with Borel sigma field
$\mathcal{B}(\textbf{X})$. We will inherently assume all functions and sets
are Borel measurable unless otherwise stated. Let $\Pi$ be a target Borel
probability measure on X and our goal is to estimate integrals through
simulation of a suitable function $\varphi:\textbf{X}\to\mathbb{R}$, that is,
to estimate

$\int_{\textbf{X}}\varphi d\Pi.$

A major motivation is estimating integrals with modern computation to perform
Bayesian inference.

The MSC estimator uses independent sample paths from a time-homogeneous Markov
chain $(X_{t})_{t=0}^{\infty}$ initialized from a random initial distribution
constructed using importance sampling. To construct the distribution
initializing the Markov chains, let $Q$ be a Borel probability measure on X,
which must be chosen beforehand, used as an importance sampling proposal.
Define the importance weight $w=d\Pi/dQ$ which is the Radon-Nikodym derivative
between $\Pi$ and $Q$. Denote $\mathbb{Z}_{+}$ as the strictly positive
integers. For $N\in\mathbb{Z}_{+}$, let $Y=(Y_{1},\ldots,Y_{N})^{T}$ where
$Y_{j}\sim Q$ are independent for $j=1,\ldots,N$ and define the self-
normalized importance weights

$w^{N}_{Y}(\cdot)=\frac{w(\cdot)}{\sum_{j=1}^{N}w(Y_{j})}.$

For a point $x\in\textbf{X}$, denote the Dirac probability measure by
$\delta_{x}$. Define the random initial probability distribution

$\Pi^{N}_{Y}=\sum_{i=1}^{N}w^{N}_{Y}(Y_{i})\delta_{Y_{i}}.$

Since $\sum_{i=1}^{N}w^{N}_{Y}(Y_{i})=1$, this defines a valid probability
measure.

For each $x\in\textbf{X}$, the conditional distribution of $X_{t}|X_{t-1}=x$
is defined by a Markov transition kernel $\mathcal{P}$ on
$\textbf{X}\times\mathcal{B}(\textbf{X})$ for $t\in\mathbb{Z}_{+}$ with
$\mathcal{P}^{1}=\mathcal{P}$ and

$\mathcal{P}^{t}(x,dy)=\int_{\textbf{X}}\mathcal{P}(z,dy)\mathcal{P}^{t-1}(x,dz).$

By this construction, we will assume $X_{t}|X_{0}=x$ is independent of $Y$
since $\mathcal{P}$ is independent of $Y$. We will also assume throughout that
the Markov kernel $\mathcal{P}$ has a unique invariant measure $\Pi$, that is,

$\displaystyle\Pi\mathcal{P}=\Pi$ (1)

is uniquely satisfied by $\Pi$. In particular, if $X_{0}\sim\Pi$, then
$X_{t}\sim\Pi$ for all $t$ thereafter. This assumption is satisfied for many
useful Markov chains and verifiable conditions to ensure uniqueness for
irreducible chains exist [Meyn and Tweedie, 2009, Theorem 10.4.4]. For
example, Metropolis-Hastings and Gibbs samplers often satisfy this assumption
and is not restrictive for many practical applications.

We can then define a Markov chain $(X_{t})_{t=0}^{\infty}$ with a random
distribution, in the sense of $Y$ being random, using initialization at
$X_{0}\sim\Pi^{N}_{Y}$ with conditional
$X_{t}|X_{t-1}=x\sim\mathcal{P}(x,\cdot)$ and marginal
$X_{t}\sim\Pi^{N}_{Y}\mathcal{P}^{t}$ defined by the random probability
measure

$\Pi^{N}_{Y}\mathcal{P}^{t}(dx^{\prime})=\int_{\textbf{X}}\mathcal{P}^{t}(x,dx^{\prime})\Pi^{N}_{Y}(dx).$

With this Markov chain, we sum the sample path of the Markov chain until its
first return to a set $C$. For a set $C\subseteq\textbf{X}$, define the first
return time by $\tau_{C}=\inf\\{t\geq 1:X_{t}\in C\\}$. For functions
$\varphi:\textbf{X}\to\mathbb{R}$, define the sum

$\displaystyle
S_{\tau_{C}}(\varphi)=\begin{cases}\sum_{k=1}^{\tau_{C}}\varphi(X_{k})&X_{0}\in
C\\\ 0&X_{0}\not\in C.\end{cases}$

Let $M$ denote the the number of independent Markov chains and conditioned on
$Y$, let $S^{m}_{\tau_{C}}(\varphi)|Y$ be independent for $m=1,\ldots,M$. The
MSC estimator is the empirical average over these independent sums of Markov
chains, that is,

$\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\varphi).$

The algorithm used to simulate the MSC estimator is outlined below in
Algorithm (1).

Simulate independent $Y_{i}\sim Q$ for $i=1,\ldots,N$ to construct
$\Pi_{Y}^{N}$

for  _$m=1,\ldots,M$_ do

Simulate an independent Markov chain path $X^{m}_{0},X_{1}^{m},\ldots
X_{\tau_{C}}^{m}$ and compute $S^{m}_{\tau_{C}}(\varphi)$

end for

Return $\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\varphi)$.

Algorithm 1 MSC estimator

##  3 Mean squared error analysis under geometric drift conditions

The main result of this section is developing a non-asymptotic analysis on the
mean squared error in Theorem 1 if the importance weights have finite variance
and under a geometric drift condition on the Markov chain. However, some
results in this section do not necessarily require the drift to be geometric
and we define the drift condition more generally. We will assume the following
general drift condition (see [Meyn and Tweedie, 2009, Chapter 11]) to control
return times of the Markov chain to a specific set.

######  Assumption 1.

(Drift condition) Suppose there are functions $V:\textbf{X}\to[1,\infty)$,
$f:\textbf{X}\to[1,\infty)$, and constants $\gamma\in(0,1)$ and
$K\in(0,\infty)$ such that for each $x\in\textbf{X}$

$(\mathcal{P}V)(x)-V(x)\leq-(1-\gamma)f(x)+K.$

Under the drift condition (1), define the set $C$ for any $R>K/(1-\gamma)$ by

$\displaystyle C=\\{x\in\textbf{X}:f(x)\leq R\\}.$ (2)

The drift condition (1) controls the return times to $C$ and in particular,
the drift condition (1) is geometric if $f=V$. In this case, the return time
to a large enough sublevel set of the function $V$ can be rapid resulting in a
short simulation of the Markov chain paths in Algorithm (1). We first
introduce a conditional measure $\mathcal{K}_{C}$ for every set
$B\subseteq\textbf{X}$ and every $x\in\textbf{X}$ by

$\displaystyle\mathcal{K}_{C}(x,B)=\mathbb{E}\left[S_{\tau_{C}}(I_{B})I_{C}(x)\bigm{|}X_{0}=x\right].$

Under Theorem 3, we bound the mean squared bias for possibly unbounded
functions dominated by $f$ defined in the drift condition (1). The result is
inspired by an importance sampling result [Agapiou et al., 2017, Theorem 2.1].

######  Proposition 1.

Assume the drift condition (1) holds, and the set $C$ is defined by (2). If
$Y_{i}\sim Q$ for $i=1,\ldots,N$ are independent and $\sup_{x\in
C}V(x)<\infty$,

$\sup_{|\varphi|\leq f}\mathbb{E}\left|\int_{\textbf{X}}\mathcal{K}_{C}\varphi
d\Pi_{Y}^{N}-\int_{\textbf{X}}\varphi d\Pi\right|^{2}\leq\frac{4A_{R}^{2}\int
wd\Pi}{N(1-\gamma_{R})^{2}}-\frac{2A_{R}^{2}}{N(1-\gamma_{R})^{2}}$

where $\gamma_{R}=\gamma+K/R$ and $A_{R}=\sup_{x\in
C}\left[V(x)+(\gamma-1)f(x)\right]+2K-1.$

An interesting property of the mean squared bias in Proposition 1 is that it
holds under relatively mild conditions such as subgeometric drift conditions.
We can also see limitations of the upper bound if $V$ is unbounded on $C$. We
now bound the variance conditioned on the importance samples $Y$. Here we will
require a geometric drift condition (1) where $f=V$ and a slightly more
restrictive class of functions $|\varphi|\leq\sqrt{V}$.

######  Proposition 2.

Assume the geometric drift condition (1) with $f=V$ holds and the set $C$ is
defined by (2). If $S^{m}_{\tau_{C}}(\varphi)|Y$ are independent for
$m=1,\ldots,M$, then

$\sup_{|\varphi|\leq\sqrt{V}}\mathbb{E}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\varphi)-\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right|^{2}\right]\leq\frac{[R+K]}{M}\left[\frac{\gamma_{R}/2}{1-\gamma_{R}/2}\right]^{2}.$

Intuitively, the upper bound on the variance will increase if the radius $R$
of the set $C$ increases, which is a sublevel set of $V$. We may now combine
the results on the mean squared bias and variance to prove an upper bound on
the mean squared error.

######  Theorem 1.

Assume the a geometric drift condition (1) with $f=V$ holds, and the set $C$
is defined by (2). If $Y_{i}$ are independent for $i=1,\ldots,N$ and
$S^{m}_{\tau_{C}}(\varphi)|Y$ are independent for $m=1,\ldots,M$, then

$\displaystyle\sup_{|\varphi|\leq\sqrt{V}}\mathbb{E}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\varphi)-\int_{\textbf{X}}\varphi
d\Pi\right|^{2}\right]$
$\displaystyle\leq\left(\frac{\gamma_{R}\sqrt{R+K}}{\sqrt{M}(2-\gamma_{R})}+\frac{2A_{R}\sqrt{\int
wd\Pi}}{\sqrt{N}(1-\gamma_{R})}\right)^{2}.$

The upper bound in Theorem 1 is a simple combination of Proposition 1 and 2.
Interestingly, the mean squared error bound holds for any sublevel set of $V$
defined in (2) so long as $R$ is large enough. To enforce shorter Markov chain
path lengths, and hence short simulation lengths, one may choose a large set
$C$. However, this upper bound illustrates the effect of choosing a larger set
$C$ as the number of independent chains will then need to increase.

An application of Theorem 1 and Markov’s inequality yields an error analysis
that can be useful for assessing the reliability of estimation in Bayesian
inference and is demonstrated in applications in Section 6. Let
$\left\lVert\cdot\right\rVert_{p}$ for $p\geq 1$ denote the standard p-norms.
Roughly speaking, the main application in mind is in Euclidean space
$\mathbb{R}^{d}$ where $\left\lVert x\right\rVert_{2}\leq\sqrt{V(x)}$. Then if
the variance of the importance weights is well-behaved, we have an explicit
requirement on $N$ and $M$ so it is guaranteed with high probability, for any
$|\varphi(x)|\leq\left\lVert x\right\rVert_{2}$, the MSC estimator

$\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\varphi)\approx\int_{\mathbb{R}^{d}}\varphi
d\Pi.$

This is of practical importance in Bayesian estimation where
$\int_{\mathbb{R}^{d}}x\Pi(dx)$ corresponds to the posterior mean and the
error analysis developed here may be applied coordinate-wise.

Alternative upper bounds to Theorem 1 can be developed as follows. If the
ratio satisfies $N/M\geq[\gamma R+2K]\int wd\Pi$, then we can obtain the
following bound by optimizing $\alpha$ below to get

$\displaystyle\mathbb{E}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S_{\tau_{C}}^{m}(\varphi)-\int_{\textbf{X}}\varphi
d\Pi\right|^{2}\right]$
$\displaystyle\leq(1+\alpha)\frac{[R+K](\gamma_{R}/2)^{2}}{M(1-\gamma_{R}/2)^{2}}+\frac{(1+\alpha^{-1})4[\gamma
R+2K]}{M(1-\gamma_{R})^{2}}$
$\displaystyle\leq\frac{\gamma_{R}R+2K}{M(1-\gamma_{R})^{2}}\left\\{(1+\alpha)\frac{1}{4}+(1+\alpha^{-1})4\right\\}$
$\displaystyle\leq\frac{6.25\gamma_{R}R+12.5K}{M(1-\gamma_{R})^{2}}.$

Here the scaling of $M$ is more clear with respect to $R,K,\gamma$, but the
ratio $N/M$ plays a part in this simplified estimate. However, this can be
estimated with importance sampling in practice to ensure such a requirement
holds.

##  4 Concentration results under multiplicative drift conditions

The goal of this section is to improve the error analysis in the previous
section but introducing stronger conditions. We look at a stronger
multiplicative drift condition (1) that is used for large deviation theory in
Markov chains [Kontoyiannis and Meyn, 2005] and has received attention
recently [Devraj et al., 2020].

######  Assumption 2.

(Multiplicative drift condition) Let $V:\textbf{X}\to[0,\infty)$ and
$f:\textbf{X}\to[0,\infty)$. Suppose for some $\gamma\in(0,1)$ and
$K\in[0,\infty)$ that for each $x\in\textbf{X}$,

$\log\\{[\mathcal{P}(\exp(V)](x)\\}\leq V(x)-(1-\gamma)f(x)+K.$

We first have a sub-Gaussian [Hoeffding, 1963] concentration bound for the
bias when the importance weights are uniformly bounded. In many applications,
assuming the importance weights are bounded is unreasonable, but if this is
condition holds, the error analysis provided in this section can produce
sharper results.

######  Proposition 3.

Assume the multiplicative drift condition (2) holds, and the set $C$ is
defined by (2). Suppose $\sup_{x\in C}V(x)<\infty$ and suppose
$\sup_{x\in\textbf{X}}w(x)=w_{*}<\infty$. If $Y_{i}\sim Q$ for $i=1,\ldots,N$
are independent, then for any $\epsilon\in(0,1)$,

$\sup_{|\varphi|\leq
f}\mathbb{P}\left(\left|\int_{\textbf{X}}\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}-\int_{\textbf{X}}\varphi d\Pi\right|\geq\epsilon\right)\leq
4\exp\left(-\frac{N\epsilon^{2}(1-\gamma_{R})^{2}}{2\exp[2B_{R}]w_{*}^{2}}\right)$

where $\gamma_{R}=\gamma+K/R$ and $B_{R}=\sup_{x\in
C}\left[V(x)-(1-\gamma)f(x)\right]+2K$.

Next, we obtain a sub-Gaussian concentration inequality under the
multiplicative drift condition.

######  Proposition 4.

Assume the multiplicative drift condition (2) holds and the set $C$ is defined
by (2). Assume $\sup_{x\in C}V(x)<\infty$. If $S^{m}_{\tau_{C}}(\varphi)|Y$
are independent for $m=1,\ldots,M$, then for any $\epsilon\in(0,1)$,

$\sup_{|\varphi|\leq
f}\mathbb{P}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\varphi)-\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right|\geq\epsilon\right]\leq
2\exp\left[-\frac{M\epsilon^{2}(1-\gamma_{R})^{2}}{9\exp(2B_{R}))}\right].$

We now have the main result of this section developing a concentration
inequality in Theorem 2 under the stronger multiplicative drift condition and
uniformly bounded importance weights.

######  Theorem 2.

Assume the drift condition (2) holds, and the set $C$ is defined by (2).
Suppose $\sup_{x\in C}V(x)<\infty$ and
$\sup_{x\in\textbf{X}}w(x)=w_{*}<\infty$. If $Y_{i}$ are independent for
$i=1,\ldots,N$ and $S^{m}_{\tau_{C}}(\varphi)|Y$ are independent for
$m=1,\ldots,M$, then for any $\epsilon\in(0,1)$,

$\displaystyle\sup_{|\varphi|\leq
f}\mathbb{P}\left(\left|\frac{1}{M}\sum_{m=1}^{M}S_{\tau_{C}}^{m}(\varphi)-\int_{\textbf{X}}\varphi
d\Pi\right|\geq\epsilon\right)\leq
6\exp\left(-\frac{\epsilon^{2}(1-\gamma_{R})^{2}\min\left\\{2M/9,N/w_{*}^{2}\right\\}}{8\exp(2B_{R})}\right).$

##  5 Toy example: Autoregressive process

Consider the simple Markov chain on $\mathbb{R}^{d}$ defined for
$\rho\in(0,1)$ and independent $\xi_{t}\sim N(0,I_{d})$ for
$t\in\mathbb{Z}_{+}$ by

$\displaystyle X_{t}=\rho X_{t-1}+\sqrt{1-\rho^{2}}\xi_{t}.$ (3)

Denote the normal distribution by $N(\mu,\Sigma)$ with mean $\mu$ and
covariance matrix $\Sigma$. With $\Gamma\equiv N(0,I_{d})$, it is readily seen
that if $X_{0}\sim\Gamma$, then $X_{t}\sim\Gamma$. We are interested in
estimating the mean $\int x\Gamma(dx)\equiv 0$ through simulation.

We will choose importance sampling proposal $\Gamma_{h}\equiv
N(0,(1/2+h)I_{d})$ with $h>0$ to sample independently $Y_{j}\sim\Gamma_{h}$
for $j=1,\ldots,N$ and construct the random initial distribution $\Pi^{N}_{Y}$
for the autoregressive process $X_{t}$ defined by (3). Using this chosen
proposal, we can directly compute

$\displaystyle\int wd\Gamma$
$\displaystyle=\frac{[1/2+h]^{d/2}}{(2\pi)^{d/2}}\int\exp\left[-\frac{1}{2}\left(\frac{2h}{1/2+h}\right)\left\lVert
x\right\rVert_{2}^{2}\right]dx=\frac{[1/2+h]^{d}}{(2h)^{d/2}}$
$\displaystyle=\left(\frac{1}{2\sqrt{2h}}+\frac{\sqrt{h}}{\sqrt{2}}\right)^{d}.$

The optimal $h=1/2$ results in the invariant distribution for the proposal.

Using the identity for the variance, we have the geometric drift condition
(1). Indeed for any $x\in\mathbb{R}^{d}$,

$\displaystyle\mathbb{E}\left[1+\left\lVert
X_{1}\right\rVert_{2}^{2}\bigm{|}X_{0}=x\right]$
$\displaystyle=\rho^{2}(1+\left\lVert
x\right\rVert_{2}^{2})+(1-\rho^{2})(1+d).$

Since the level sets of the drift function are compact, it can be shown the
Markov chain converges geometrically fast and the invariant measure is unique
[Hairer and Mattingly, 2011]. For any $r>1$, define the set
$C_{r}=\\{x\in\mathbb{R}^{d}:\left\lVert x\right\rVert_{2}^{2}\leq rd\\}$ and
let $S^{m}_{\tau_{C_{r}}}|Y$ independent $m=1,\ldots,N$ using the
autoregressive process (3).

The drift controls the return times to the set $C_{r}$ and combining these two
results, we can apply Theorem 1 to estimate functions satisfying
$|\varphi(x)|\leq\sqrt{1+\left\lVert x\right\rVert_{2}^{2}}$. Indeed, Theorem
1 says that if

$\displaystyle
N\propto\frac{r^{2}d^{2}}{\delta\epsilon^{2}(1-\rho_{r})^{2}}\left(\frac{1}{2\sqrt{2h}}+\sqrt{\frac{h}{2}}\right)^{d}$
$\displaystyle M\propto\frac{rd}{\delta\epsilon^{2}(2-\rho_{r})^{2}}$ (4)

where $\rho_{r}=\rho^{2}+(1-\rho^{2})(d+1)/(rd+1)$, then with probability at
least $1-\delta$,

$\left|\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C_{r}}}(\varphi)-\int_{\mathbb{R}^{d}}\varphi
d\Gamma\right|\leq\epsilon.$

Figure 1 plots the precise theoretical lower bounds (4) on the number of
samples for $N$ and $M$ to increasing dimensions $\\{1,5,10,15,20,25,30\\}$
with tuning parameters $\epsilon=\delta=.1$, $\rho=.9$, $h=.49$, and $r=1.5$.
We can see a linear scaling with the dimension for the number of Markov chains
$M$. We also see roughly quadratic scaling with respect to the dimension for
$N$ although the lower bound on the required number of samples $N$ in the
initialization is much larger. However, alternative tuning parameter choices
on the proposal may lead to approximately exponential scaling with respect to
the dimension.

Figure 1:  Computation of the lower bound on the number of Markov chains $M$
and initialization samples $N$ for the MSC estimator using the autoregressive
process with respect to increasing dimensions.

We look to investigate the empirical estimation performance of the MSC
estimator using this autoregressive process. We attempt to estimate
coordinates of the mean so that $\varphi(x)=x_{k}$ for $k\leq d$ via an MSC
simulation. Use the tuning parameters introduced previously, we simulate the
MSC estimate with $N=10^{6}$ and $M=10^{5}$ independent samples in dimension
$d=2$. In this case, we can directly sample the target distribution and
compare to standard Monte Carlo estimation.

A comparison is made to the MSC simulation estimate within 2 estimated
standard errors to standard Monte Carlo realizations and its estimated
standard error. Figures 2(a) and 2(b) illustrate the results of the simulation
by plotting both coordinates individually. Compared to standard Monte Carlo,
the MSC estimate shows increased variability but can be utlized when standard
Monte Carlo is not available. Figure 2(c) shows the run length of the
independent Markov chain simulations until each Markov chain returns to the
set $C_{r}$ defined in this section. Here we observe an alignment with the
theoretical analysis that the Markov chain has a rapid return time and remains
in the many-short-chains regime. In simulations not presented here, we
observed it is possible to increase the radius of the set $C_{r}$ to further
reduce the run length while also observing stable estimation performance.

(a)

(b)

(c)

Figure 2:  Figures 2(a) and Figure 2(b) plot the coordinates of the MSC mean
within 2 estimated standard errors compared to standard MC. Figure 2(b) plots
the length of the Markov sample paths over the independent runs of the Markov
chains.

##  6 Cardiovascular disease prediction with the Pólya-Gamma Gibbs sampler

We look now to a non-trivial example on a real data set with Bayesian logistic
regression and the Pólya-Gamma Gibbs sampler [Nicholas G. Polson and Windle,
2013]. A convergence analysis in the single-chain regime has been developed
[Choi and Hobert, 2013], but the scaling of the constants in the convergence
rate may lead to a long simulation length in moderate to high dimensions. It
may then be an advantage to incorporate parallel computation using many-short-
chains for this Markov chain that does not require such a convergence
analysis.

Here we use standard notation in Bayesian inference although it overlaps
slightly with our previous notation. Let $(Y_{i},X_{i})_{i=1}^{n}$ with
$Y_{i}\in\\{0,1\\}$ and $X_{i}\in\mathbb{R}^{d}$ and define
$Y=(Y_{1},\ldots,Y_{n})^{T}\in\\{0,1\\}^{n}$ and
$X=(X_{1},\ldots,X_{n})^{T}\in\mathbb{R}^{n\times d}$. For $x\in\mathbb{R}$,
define the sigmoid function by $\sigma(x)=(1+\exp(-x))$. Consider the Bayesian
logistic regression model with a Gaussian prior satisfying

$\displaystyle Y_{i}|X_{i},\beta\sim\text{Bern}(\sigma(\beta^{T}X_{i}))$
$\displaystyle\beta\sim N(0,\Sigma)$

where $\Sigma\in\mathbb{R}^{d\times d}$ is a symmetric, positive-definite
(SPD) covariance matrix. Define the negative log-likelihood by
$\ell_{n}(\beta)=\sum_{i=1}^{n}\left[\log(1+\exp(X_{i}^{T}\beta))-Y_{i}X_{i}^{T}\beta\right].$
The posterior for this model $\Pi_{n}$ has a Lebesgue density defined by

$\pi_{n}(\beta)\propto\exp\left[-\ell_{n}(\beta)-\frac{1}{2}\beta^{T}\Sigma^{-1}\beta\right].$

We will be interested in estimating the posterior mean of the coefficients
$\int\beta\Pi_{n}(d\beta)$ to perform Bayesian inference.

The Pólya-Gamma Gibbs sampler uses Pólya-Gamma random variables to construct
an augmented two-variable deterministic scan Gibbs sampler. This results in a
Markov chain $(\omega_{t},\beta_{t})_{t=0}^{\infty}$ and a marginal chain
$(\beta_{t})_{t=0}^{\infty}$ that has the posterior as its invariant
distribution. Let $\text{PG}(1,b)$ denote the Pólya-Gamma distribution with
parameter $b\geq 0$ [Nicholas G. Polson and Windle, 2013]. For
$\omega\in\mathbb{R}^{n}$, define

$\displaystyle\Omega=\text{diag}(\omega)$
$\displaystyle\Sigma(\omega)=(X^{T}\Omega X+\Sigma^{-1})^{-1}$
$\displaystyle\mu(\omega)=\Sigma(\omega)X^{T}(Y-(1/2)1_{d}).$

For $t=1,\ldots$, this Pólya-Gamma Markov chain alternates

$\displaystyle\omega_{t}=(\omega^{1}_{t},\ldots,\omega^{n}_{t})|\beta_{t-1}$
with independent
$\displaystyle\omega^{i}_{t}\sim\text{PG}(1,|X_{i}^{T}\beta_{t-1}|)$
$\displaystyle\beta_{t}|\omega_{t}\sim
N\left[\mu(\omega_{t}),\Sigma(\omega_{t})\right].$

The posterior density is log-concave and we prove a general result to choose
the importance sampling proposal in this scenario. For functions
$g:\mathbb{R}^{d}\to\mathbb{R}$, denote the gradient and Hessian at
$x\in\mathbb{R}^{d}$ by $\nabla g(x)$ and $\nabla^{2}g(x)$ respectively. For
an arbitrary matrix $A$, let $\left\lVert A\right\rVert_{2}$ denote the square
root of the largest eigenvalue of $A^{T}A$. The following Proposition allows
choosing a proposal for many log-concave distributions and we will apply it to
this example.

######  Proposition 5.

Let $\ell:\mathbb{R}^{d}\to\mathbb{R}$ and suppose $\ell$ is a twice
continuously differentiable convex function such that for some $\lambda^{*}>0$

$\displaystyle\sup_{x\in\mathbb{R}^{d}}\left\lVert\nabla^{2}\ell(x)\right\rVert_{2}\leq\lambda^{*}.$
(5)

Let $\Sigma$ be a SPD matrix and define a function $f$ by

$\displaystyle f(x)=\ell(x)+\frac{1}{2}x^{T}\Sigma^{-1}x$

and define a probability density by $\pi(x)=\exp(-f(x))/Z$ where the
normalizing constant $Z=\int\exp(-f(x))dx$. Let
$x^{*}=\text{argmin}\\{f(x):x\in\mathbb{R}^{d}\\}$ and for $h\in(0,1/2]$, let
$Q\sim N(x^{*},(1/2+h)\Sigma)$ with density $q$. Then the upper bound holds

$\int_{\mathbb{R}^{d}}\frac{\pi(x)^{2}}{q(x)}dx\leq\det(\lambda^{*}\Sigma+I_{d})\left(\frac{1}{2\sqrt{2h}}+\sqrt{\frac{h}{2}}\right)^{d}.$

Let $\beta^{*}_{n}$ be the maximum of the posterior density $\pi_{n}$, which
always exists. With Proposition 5 in mind, we will choose the importance
sampling proposal $Q_{h}\equiv N(\beta_{n}^{*},(1/2+h)\Sigma)$ with
$h\in(0,1/2]$. Let $\theta=(\theta_{1},\ldots,\theta_{N})$ where
$\theta_{j}\sim Q_{h}$ are independent to construct the random initial
distribution $\Pi_{\theta}^{N}$. Define
$L=\left\lVert\Sigma\right\rVert_{2}^{2}\left\lVert
X^{T}(Y-1/21_{d})\right\rVert_{2}^{2}$ and for any $r>1$, define the set

$\displaystyle
C_{r}=\\{\beta\in\mathbb{R}^{d}:\left\lVert\beta\right\rVert_{2}^{2}\leq
rL\\}.$ (6)

A drift condition is shown in Proposition 6 to this set $C_{r}$ and yields the
main result for the error analysis of the MSC estimator using this Pólya-Gamma
sampler.

######  Proposition 6.

With the set $C_{r}$ defined by (6), let $\theta_{j}\sim Q_{h}$ be independent
for $j=1,\ldots,N$ and $S^{m}_{\tau_{C_{r}}}(\varphi)|\theta$ be independent
$m=1,\ldots,M$ using marginal Pólya-Gamma Markov chains
$(\beta^{m}_{t})_{t=1}^{\infty}$ initialized at
$\beta^{m}_{0}\sim\Pi_{\theta}^{N}$. Then for all functions with
$|\varphi(\beta)|\leq\sqrt{\left\lVert\beta\right\rVert_{2}^{2}+1}$,

$\displaystyle\mathbb{E}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C_{r}}}(\varphi)-\int_{\textbf{X}}\varphi
d\Pi\right|^{2}\right]\leq\left(\frac{\gamma_{r}\sqrt{(r+1)L+2}}{\sqrt{M}(2-\gamma_{r})}+\frac{[2L+1]W_{d}}{\sqrt{N}(1-\gamma_{r})}\right)^{2}.$

where $\gamma_{r}=(1+L)/(1+rL)$ and

$W_{d}=\det(\left\lVert
X\right\rVert_{2}^{2}\Sigma/4+I_{d})\left(\frac{1}{2\sqrt{2h}}+\sqrt{\frac{h}{2}}\right)^{d}.$

We now investigate the empirical performance on a real data set of the MSC
estimator from Proposition 6. The aim is to predict cardiovascular disease
from a well-known data set provided by the Cleveland Clinic [Detrano et al.,
1989]. The data consists of $303$ patients with binary responses determining
if cardiovascular disease is present and $14$ predictor variables based on
patient characteristics. After converting categorical predictors, a total of
$21$ covariates are used. We will investigate the coefficients $\beta_{1}$
indicating male versus female patients, $\beta_{2}$ for the number of major
blood vessels colored by flourosopy, and $\beta_{3}$ for resting blood
pressure in mm/Hg on admission to the hospital. For the prior, a covariance
$\Sigma=10I$ is chosen.

We compare the MSC estimation performance to a single-chain Pólya-Gamma Gibbs
simulation to estimate the posterior mean of the covariates of interest. For
an additional comparison, a tuned random-walk Metropolis-Hastings (RWM) is
also simulated using scaling parameter according to optimal scaling [Roberts
et al., 1997]. We generate the Pólya-Gamma Gibbs sampler and RWM for $10^{6}$
iterations starting from the random initial distribution used the MSC
estimation. For the MSC estimate, $N=10^{7}$ independent samples are used for
the initial distribution and $M=10^{6}$ independent Markov chains are used
with tuning parameters $h=.49$ and $r=1.001$.

Figure 3 shows the means of the MSC estimates within 2 estimated standard
errors of the simulations compared to the means for the single-chain
simulations. Similar performance is observed of all algorithms up to a small
Monte Carlo error. The drift condition here results in a large radius for the
set $C_{r}$ and the simulation results in Markov chain paths of essentially
length $1$. In particular, the estimate in the simulation is an empirical mean
of a combination of only $1$ step from the Gibbs sampler. However, the initial
distribution must be carefully constructed for the MSC estimate when compared
to the single-chain MCMC algorithms.

(a)

(b)

(c)

Figure 3:  The mean within 2 estimated standard deviations of MSC the Pólya-
Gamma sampler using compared to the single-chain Pólya-Gamma sampler and RWM
algorithm.

##  7 Discussion and conclusions

The MSC estimator based on the representation [Meyn and Tweedie, 2009, Theorem
10.4.9] provides interesting insight into the many-short-chains regime. It is
remarkable that the MSC estimator relies only on the drift condition for the
Markov chain and not its convergence rate where the latter often scales poorly
with the dimension. Another advantage of the MSC estimator is a simple
uncertainty quantification by estimating the standard error since it is based
on independent Markov chains. In comparison, determining the termination time
of an MCMC simulation in the single-chain regime is a challenging problem and
actively researched.

The parallel MSC estimator developed here may be practically useful when a
reasonable importance sampling proposal is available along with an explicit
geometric or multiplicative drift condition. However, this may not always be
the case and other parallel MCMC estimation techniques may be suitable instead
[Jacob et al., 2020a]. For moderately sized problems, the ability to utilize
parallel computation is an important property. However, for some high
dimensional problems, the non-asymptotic bounds may not be reasonable unless
strong conditions are satisfied so that the non-asymptotic bounds scale
appropriately. There are also some criticisms and open questions on the
benefit of using the MSC estimator over self-normalized importance sampling.
The current importance sampling bounds [Agapiou et al., 2017, Theorem 2.1]
hold for bounded functions and require stronger conditions for unbounded
functions [Agapiou et al., 2017, Theorem 2.3]. An alternative MSC estimator
can also be defined based on a second representation of the invariant measure
and this appears to bear some further relations to importance sampling [Meyn
and Tweedie, 2009, Theorem 10.4.9].

There appear to be many future research directions. For example, the optimal
set $C$ and drift condition is not understood, but this is also the case in
convergence analysis using drift and minorization conditions for Markov
chains. Empirically, the MSC estimator has different run lengths based on the
size of the set $C$ and benefits and drawbacks of this still still remain
somewhat elusive. Another possible research direction is to develop error
bounds for not only estimation but not for stronger discrepancies such as
Wasserstein distances and other probability discrepancies. However, the lower
bounds for such distances are poor in high dimensions [Dudley, 1969].

## Acknowledgements

Thanks to Radu Craiu, Jeffrey Rosenthal, Qiang Sun, Galin Jones, and many
others at the University of Toronto for insightful conversations while writing
this manuscript.

## Supplementary material

The Pypolyagamma library https://github.com/slinderman/pypolyagamma is used by
Scott Linderman to generate Pólya random variables. Simulation and all code is
available https://github.com/austindavidbrown/msc-estimator. The
cardiovascular data is available at the UCI database
https://archive.ics.uci.edu/dataset/45/heart+disease.

##  Appendix A Supporting technical results for Section 3

The following bound is well-known.

######  Lemma 1.

Assume the drift condition (1) holds. Then

$\int_{\textbf{X}}fd\Pi\leq\frac{K}{1-\gamma}.$

###### Proof.

The drift condition implies

$\displaystyle P^{n}V(x)-V(x)$
$\displaystyle=\sum_{k=1}^{n-1}\left[P^{k+1}V(x)-P^{k}V(x)\right]\leq-(1-\gamma)\sum_{k=0}^{n-1}P^{k}f(x)+nK.$

By Fatou’s lemma [Bogachev and Ruas, 2007, Theorem 2.8.3]

$\displaystyle\int_{\textbf{X}}fd\Pi=\limsup_{n}\int_{\textbf{X}}\frac{1}{n}\sum_{k=0}^{n-1}P^{k}f(x)d\Pi\leq\int_{\textbf{X}}\limsup_{n}\frac{1}{n}\sum_{k=0}^{n-1}P^{k}f(x)d\Pi\leq\frac{K}{1-\gamma}.$

∎

The following bound in Lemma 2 is implied by the drift condition. The result
is similar to [Meyn and Tweedie, 2009, Proposition 11.3.2]. Denote $a\wedge
b,a\vee b$ for the minimum and maximum of two numbers $a,b$.

######  Lemma 2.

If the drift condition (1) holds, then for every $x\in C$,

$\sup_{|\varphi|\leq
f}(\mathcal{K}_{C}|\varphi|)(x)\leq\frac{V(x)-1-(1-\gamma)f(x)+2K}{1-\gamma_{R}}$

where $\gamma_{R}=\gamma+K/R$.

###### Proof of Lemma 2.

The drift condition (1) implies

$\mathcal{P}V(x)-V(x)\leq-(1-\gamma_{R})f(x)+K(1-\frac{f(x)}{R})\leq-(1-\gamma_{R})f(x)+KI_{C}(x).$

For $n\in\mathbb{Z}_{+}$, define

$\displaystyle
M_{n+1}=V(X_{n+1})+(1-\gamma_{R})\sum_{k=1}^{n}f(X_{k})-K\sum_{k=1}^{n}I_{C}(X_{k}).$
(7)

Using the drift condition, $M_{n+1}$ is a supermartingale with respect to the
standard filtration $\mathcal{F}_{n+1}$ for $X_{0},X_{1},\ldots,X_{n+1}$. By
the optional sampling theorem [Kallenberg, 2021, Theorem 9.12],
$M_{n\wedge\tau_{C}+1}$ is also a supermartingale with its respective
filtration and

$\displaystyle(1-\gamma_{R})\mathbb{E}\left[\sum_{k=1}^{\tau_{C}\wedge
n}f(X_{n})\bigm{|}X_{0}=x\right]+1-K\leq\mathbb{E}[M_{\tau_{c}\wedge
n+1}\bigm{|}X_{0}=x]$
$\displaystyle\leq\mathbb{E}[M_{2}\bigm{|}X_{0}=x]\leq\mathbb{E}[V(X_{1})\bigm{|}X_{0}=x].$

Since $f\geq 1$, using monotone convergence [Bogachev and Ruas, 2007, Theorem
2.8.2],

$\displaystyle\mathbb{E}\left[\sum_{k=1}^{\tau_{C}}f(X_{n})\bigm{|}X_{0}=x\right]\leq\frac{\mathbb{E}[V(X_{1})\bigm{|}X_{0}=x]-1+K}{1-\gamma_{R}}$
$\displaystyle\leq\frac{V(x)-1-(1-\gamma)f(x)+2K}{1-\gamma_{R}}.$

∎

The following result shows that $\Pi$ is invariant for $\mathcal{K}_{C}$ and
we have an identity for many integrals with respect to the invariant measure
under the drift condition (1).

######  Theorem 3.

Assume the Markov kernel $\mathcal{P}$ has unique invariant measure $\Pi$ and
the drift condition (1) holds. If $\int_{C}Vd\Pi<\infty$, then for every
$\varphi:\textbf{X}\to\mathbb{R}$ such that $|\varphi|\leq f$, we have
$\int_{\textbf{X}}\varphi d\Pi<\infty$ and

$\displaystyle\int_{\textbf{X}}\mathcal{K}_{C}\varphi
d\Pi=\int_{\textbf{X}}\varphi d\Pi.$

###### Proof for Theorem 3.

Under the drift condition, then irreducibile Markov chains with unique
invariant distribution have a representation of the invariant measure $\Pi$
for non-negative simple functions [Meyn and Tweedie, 2009, Theorem 10.4.9].
The assumption of an irreducible Markov chain is not required as seen in the
proof [Meyn and Tweedie, 2009, Theorem 10.4.9] and we include the proof only
for completeness. Let $\varphi:\textbf{X}\to[0,\infty]$ be simple so that
$\varphi=\sum_{k=1}^{J}a_{k}I_{B_{k}}$ for some constants $a_{k}\geq 0$ and
sets $B_{k}$. By invariance and induction it is readily shown that

$\int\varphi
d\Pi\geq\sum_{k=1}^{n}\int_{C}\mathbb{E}\left[\varphi(X_{k})I_{\tau_{C}\geq
k}\bigm{|}X_{0}=x\right]\Pi(dx).$

For sets $B\subseteq\textbf{X}$, this implies $\Pi(B\cap
C)\geq\int_{C}\mathbb{P}\left[X_{\tau_{C}}\in
B,\tau_{C}<\infty\bigm{|}X_{0}=x\right]\Pi(dx).$ Suppose by contradiction that
there is some set $B_{0}\subseteq C$ with

$\Pi(B_{0})>\int_{C}\mathbb{P}\left[X_{\tau_{C}}\in
B_{0},\tau_{C}<\infty\bigm{|}X_{0}=x\right]\Pi(dx).$

We must have

$\displaystyle\Pi(C)=\Pi(B_{0})+\Pi(C\cap
B_{0}^{c})>\int_{C}\mathbb{P}\left[\tau_{C}<\infty\bigm{|}X_{0}=x\right]\Pi(dx).$

Since $f\geq 1$, the drift condition implies
$\mathbb{P}\left[\tau_{C}<\infty\bigm{|}X_{0}=x\right]=1$ for all $x\in C$ and
so this is a contradiction and for every set $B\subseteq\textbf{X}$

$\Pi(B\cap C)\leq\int_{C}\mathbb{P}\left[X_{\tau_{C}}\in
B,\tau_{C}<\infty\bigm{|}X_{0}=x\right]\Pi(dx)\leq\Pi(B\cap C).$

Combining these inequalities implies

$\displaystyle\int_{C}\mathbb{E}\left[\varphi(X_{\tau_{C}})I_{\tau_{C}<\infty}\bigm{|}X_{0}=x\right]\Pi(dx)=\int_{C}\varphi
d\Pi.$ (8)

Combining (8) and the Markov property

$\displaystyle\int_{C}\mathbb{E}\left[\sum_{k=1}^{\tau_{C}}\mathcal{P}\varphi(X_{k})\bigm{|}X_{0}=x\right]\Pi(dx)$
$\displaystyle=\sum_{k=1}^{\infty}\int_{C}\mathbb{E}\left[\mathcal{P}\varphi(X_{k+1})I_{\tau_{C}\geq
k+1}\bigm{|}X_{0}=x\right]\Pi(dx)+\sum_{k=1}^{\infty}\int_{C}\mathbb{E}\left[\mathcal{P}\varphi(X_{k})I_{\tau_{C}=k}\bigm{|}X_{0}=x\right]\Pi(dx)$
$\displaystyle=\int_{C}\mathbb{E}\left[\sum_{k=2}^{\tau_{C}}\varphi(X_{k})\bigm{|}X_{0}=x\right]\Pi(dx)+\int_{C}\mathbb{E}\left[\mathcal{P}\varphi(X_{\tau_{C}})I_{\tau_{C}<\infty}\bigm{|}X_{0}=x\right]\Pi(dx)$
$\displaystyle=\int_{C}\mathbb{E}\left[\sum_{k=2}^{\tau_{C}}\varphi(X_{k})\bigm{|}X_{0}=x\right]\Pi(dx)+\int_{C}\mathcal{P}\varphi
d\Pi$
$\displaystyle=\int_{C}\mathbb{E}\left[\sum_{k=1}^{\tau_{C}}\varphi(X_{k})\bigm{|}X_{0}=x\right]\Pi(dx).$

Since $\Pi$ is the unique invariant measure and under the drift condition, we
have

$\int_{\textbf{X}}\varphi
d\Pi=\int_{C}\mathbb{E}\left[\sum_{k=1}^{\tau_{C}}\varphi(X_{k})\bigm{|}X_{0}=x\right]\Pi(dx).$

This representation extends to simple functions
$\varphi:\textbf{X}\to\mathbb{R}$ taking the positive and negative parts.

It remains to extend this representation to more general functions. Now let
$\varphi:\textbf{X}\to\mathbb{R}$ such that $|\varphi|\leq f$ and Lemma 1
implies $\int_{\textbf{X}}|\varphi|d\Pi\leq K/(1-\gamma).$ Since $\varphi$ is
measurable, we can choose a sequence of simple functions
$\varphi_{n}\to\varphi$ pointwise with $|\varphi_{n}|\leq|\varphi_{n+1}|$.
Using Lemma 2 and dominated convergence [Bogachev and Ruas, 2007, Theorem
2.8.1], for each $x\in C$
$\lim_{n\to\infty}(\mathcal{K}_{C}\varphi_{n})(x)=(\mathcal{K}_{C}\varphi)(x).$
Using Lemma 2 and since we have assumed $\int_{C}Vd\Pi<\infty$, we have by
dominated convergence [Bogachev and Ruas, 2007, Theorem 2.8.1],

$\displaystyle\int_{\textbf{X}}\mathcal{K}_{C}\varphi d\Pi$
$\displaystyle=\int_{C}\lim_{n\to\infty}(\mathcal{K}_{C}\varphi_{n})d\Pi=\lim_{n\to\infty}\int_{C}(\mathcal{K}_{C}\varphi_{n})\Pi(dx)$
$\displaystyle=\lim_{n\to\infty}\int_{\textbf{X}}\varphi_{n}d\Pi=\int_{\textbf{X}}\lim_{n\to\infty}\varphi_{n}d\Pi=\int_{\textbf{X}}\varphi
d\Pi.$

∎

###### Proof of Proposition 1.

Lemma 2 implies

$\displaystyle\sup_{x\in\textbf{X}}|(\mathcal{K}_{C}\varphi)(x)|\leq\frac{A_{R}}{1-\gamma_{R}}.$
(9)

Since we assumed $\sup_{x\in C}V(x)<\infty$, then $A_{R}<\infty$ and we have
shown that $\mathcal{K}_{C}\varphi$ is bounded. Here we used that $x$ is
restricted to $C$ due to the definition of $\mathcal{K}_{C}$. Since we have
assumed the drift condition, the invariant distribution is unique, and
$\sup_{x\in C}V(x)<\infty$, then Theorem 3 implies the identity

$\mathbb{E}[(\mathcal{K}_{C}\varphi)(Y_{1})w(Y_{1})]=\int_{C}\mathcal{K}_{C}\varphi
d\Pi=\int_{\textbf{X}}\varphi d\Pi.$

By [Agapiou et al., 2017, Theorem 2.1], we have the error bound for self-
normalized importance sampling with the bounded function
$\mathcal{K}_{C}\varphi$. We can slightly improve the bound, but the proof
technique is essentially the same. Using Cauchy-Schwarz, for any $\alpha>0$,
we have the upper bound

$\displaystyle\mathbb{E}\left|\sum_{i=1}^{N}w^{N}_{Y}(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\mathbb{E}[(\mathcal{K}_{C}\varphi)w]\right|^{2}$
(10)
$\displaystyle\leq(1+\alpha)\mathbb{E}\left|\sum_{i=1}^{N}w^{N}_{Y}(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\sum_{i=1}^{N}\frac{w(Y_{i})}{N}(\mathcal{K}_{C}\varphi)(Y_{i})\right|^{2}$
(11) $\displaystyle\hskip
10.00002pt+(1+\alpha^{-1})\mathbb{E}\left|\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\mathbb{E}[(\mathcal{K}_{C}\varphi)w]\right|^{2}.$
(12)

This second term in (12) is the variance of independent random variables:

$\displaystyle\text{Var}\left[\frac{1}{N}\sum_{i=1}^{N}(\mathcal{K}_{C}\varphi)(Y_{i})w(Y_{i})\right]$
$\displaystyle=\frac{1}{N}\text{Var}\left[(\mathcal{K}_{C}\varphi)(Y)w(Y)\right]\leq[\sup_{x\in\textbf{X}}|(\mathcal{K}_{C}\varphi)(x)|]^{2}\frac{\mathbb{E}\left[w(Y_{1})^{2}\right]}{N}.$

For the first term in (12) , we have

$\displaystyle\left|\sum_{i=1}^{N}\frac{w(Y_{i})}{\sum_{j}w(Y_{j})}\mathcal{K}_{C}\varphi(Y_{i})-\sum_{i=1}^{N}\frac{w(Y_{i})}{N}\mathcal{K}_{C}\varphi(Y_{i})\right|$
$\displaystyle=\left|\sum_{i=1}^{N}\mathcal{K}_{C}\varphi(Y_{i})\frac{w(Y_{i})}{\sum_{j}w(Y_{j})}\left[\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})-1\right]\right|$
$\displaystyle\leq\max_{i\in
1,\ldots,N}|(\mathcal{K}_{C}\varphi)(X_{i})|\left|\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})-1\right|.$

By a standard Monte Carlo argument,
$\text{Var}[\overline{w}]=\text{Var}[w(Y_{1})]/N,$ it follows then that

$\displaystyle\mathbb{E}\left|\sum_{i=1}^{N}\frac{w(Y_{i})}{\sum_{j}w(Y_{j})}(\mathcal{K}_{C}\varphi)(Y_{i})-\sum_{i=1}^{N}\frac{w(Y_{i})}{N}(\mathcal{K}_{C}\varphi)(Y_{i})\right|^{2}\leq[\sup_{x\in\textbf{X}}|(\mathcal{K}_{C}f)(x)|]^{2}\frac{\text{Var}[w]}{N}.$

Combining these bounds choosing $\alpha=1$, we have

$\displaystyle\mathbb{E}\left|\sum_{i=1}^{N}w^{N}_{Y}(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\mathbb{E}[(\mathcal{K}_{C}\varphi)w]\right|^{2}$
$\displaystyle\leq\frac{2\int
wd\Pi}{N}\frac{A_{R}^{2}}{(1-\gamma_{R})^{2}}+\frac{2[\int
wd\Pi-1]}{N}\frac{A_{R}^{2}}{(1-\gamma_{R})^{2}}$
$\displaystyle\leq\frac{4A_{R}^{2}\int
wd\Pi}{N(1-\gamma_{R})^{2}}-\frac{2A_{R}^{2}}{N(1-\gamma_{R})^{2}}.$

∎

###### Proof of Proposition 2.

Since $S_{\tau_{C}}^{m}(\varphi)|Y$ are assumed independent, we have
$\mathbb{E}[S_{\tau_{C}}^{m}(\varphi)|Y]=\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}$ and the variance bound

$\mathbb{E}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\varphi)-\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right|^{2}\Bigm{|}Y\right]\leq\frac{\mathbb{E}\left[|S^{m}_{\tau_{C}}(\varphi)|^{2}\Bigm{|}Y\right]}{M}.$

Let $\mathcal{F}_{n}$ denote the standard filtration with respect to
$X_{0},\ldots,X_{n}$. Repeated application of the geometric drift condition
for $x\in\textbf{X}$ gives

$\displaystyle\mathbb{E}[V(X_{k})I_{\tau_{C}\geq k}\bigm{|}X_{0}=x,Y]$
$\displaystyle=\mathbb{E}\left[I_{\tau_{C}\geq
k}\mathbb{E}(V(X_{k})\bigm{|}\mathcal{F}_{k-1},X_{0}=x,Y)|Y\right]$
$\displaystyle\leq\mathbb{E}\left[\gamma_{R}V(X_{k-1})I_{\tau_{C}\geq
k}+KI_{C}(X_{k-1})I_{\tau_{C}\geq k}\bigm{|}X_{0}=x,Y\right]$
$\displaystyle\leq\gamma_{R}\mathbb{E}\left[V(X_{k-1})I_{\tau_{C}\geq
k-1}\bigm{|}X_{0}=x,Y\right]$
$\displaystyle\leq\gamma_{R}^{k}[V(x)+KI_{C}(x)].$

Here we used that for any $k\in\mathbb{Z}_{+}$, $I_{\tau_{C}\geq k}\leq
I_{\tau_{C}\geq k-1}$ and $\\{\tau_{C}\geq k\\}$ is $\mathcal{F}_{k-1}$
measurable.

Using the Cauchy-Schwarz inequality,

$\displaystyle\mathbb{E}\left[\left(\sum_{k=1}^{n}\varphi(X_{k})I_{\tau_{C}\geq
k}\right)^{2}\bigm{|}X_{0},Y\right]$
$\displaystyle=\sum_{k=1}^{n}\sum_{\ell=1}^{n}\mathbb{E}\left[\varphi(X_{k})I_{\tau_{C}\geq
k}\varphi(X_{\ell})I_{\tau_{C}\geq\ell}\bigm{|}X_{0},Y\right]$
$\displaystyle\leq\sum_{k=1}^{n}\sum_{\ell=1}^{n}\sqrt{\mathbb{E}\left[\varphi(X_{k})^{2}I_{\tau_{C}\geq
k}\bigm{|}X_{0},Y\right]}\sqrt{\mathbb{E}\left[\varphi(X_{\ell})^{2}I_{\tau_{C}\geq\ell}\bigm{|}X_{0},Y\right]}.$

Now applying the drift condition,

$\displaystyle\mathbb{E}\left[\left(\sum_{k=1}^{\tau_{C}\wedge
n}\varphi(X_{k})\right)^{2}\bigm{|}X_{0},Y\right]$
$\displaystyle\leq\sum_{k=1}^{n}\sum_{\ell=1}^{n}\sqrt{\mathbb{E}\left[V(X_{k})I_{\tau_{C}\geq
k}\bigm{|}X_{0},Y\right]}\sqrt{\mathbb{E}\left[V(X_{\ell})I_{\tau_{C}\geq\ell}\bigm{|}X_{0},Y\right]}$
$\displaystyle\leq[V(X_{0})+K]\sum_{k=1}^{n}\sum_{\ell=1}^{n}\gamma_{R}^{k/2}\gamma_{R}^{\ell/2}$
$\displaystyle\leq[V(X_{0})+K]\sum_{k=1}^{n}\sum_{\ell=1}^{n}\gamma_{R}^{k/2}\gamma_{R}^{\ell/2}$
$\displaystyle\leq[V(X_{0})+K]\left[\frac{\gamma_{R}/2}{1-\gamma_{R}/2}\right]^{2}.$

Using Fatou’s lemma [Bogachev and Ruas, 2007, Theorem 2.8.3] and since
$\sup_{x\in C}V(x)\leq R$, we conclude

$\mathbb{E}\left[S^{m}_{\tau_{C}}(\varphi)^{2}\bigm{|}X_{0},Y\right]\leq[R+K]\left[\frac{\gamma_{R}/2}{1-\gamma_{R}/2}\right]^{2}.$

Taking the iterated expectation with respect to $X_{0}\sim\Pi^{N}_{Y}$ and
then $Y_{j}\sim Q$ for $j=1,\ldots,N$, completes the proof. ∎

###### Proof of Theorem 1.

Since $V=f\geq 1$, then $|\varphi|\leq\sqrt{V}\leq V$ and also $\sup_{x\in
C}V(x)\leq R$, so we are able to apply both upper bounds in Proposition 1 and
Proposition 2. By Lemma 1,
$\int_{\textbf{X}}|\varphi|d\Pi\leq\int_{\textbf{X}}Vd\Pi\leq K/(1-\gamma).$
Using Cauchy-Schwarz for any $\alpha>0$, we have the upper bound

$\displaystyle\mathbb{E}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S_{\tau_{C}}^{m}(\varphi)-\int_{\textbf{X}}\varphi
d\Pi\right|^{2}\right]$
$\displaystyle\leq(1+\alpha)\mathbb{E}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S_{\tau_{C}}^{m}(\varphi)-\int_{\textbf{X}}\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right|^{2}\right]$ $\displaystyle\hskip
10.00002pt+(1+\alpha^{-1})\mathbb{E}\left[\left|\int_{\textbf{X}}\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}-\int_{\textbf{X}}\varphi d\Pi\right|^{2}\right].$

Using Proposition 1, Proposition 2 and optimizing over $\alpha$,

$\displaystyle\mathbb{E}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S_{\tau_{C}}^{m}(\varphi)-\int_{\textbf{X}}\varphi
d\Pi\right|^{2}\right]$
$\displaystyle\leq(1+\alpha)\frac{[R+K]}{M}\left[\frac{\gamma_{R}/2}{1-\gamma_{R}/2}\right]^{2}+(1+\alpha^{-1})\frac{4A_{R}^{2}\int
wd\Pi}{N(1-\gamma_{R})^{2}}$
$\displaystyle\leq\left(\frac{\gamma_{R}\sqrt{R+K}}{\sqrt{M}(2-\gamma_{R})}+\frac{2A_{R}\sqrt{\int
wd\Pi}}{\sqrt{N}(1-\gamma_{R})}\right)^{2}.$

∎

##  Appendix B Supporting technical results for Section 4

Similar to Lemma 2, the following upper bound is implied by the drift
condition.

######  Lemma 3.

If the multiplicative drift condition (2) holds, then for every $x\in C$

$\displaystyle\mathbb{E}\left(\exp\left[(1-\gamma_{R})\sum_{k=1}^{\tau_{C}}f(X_{k})\bigm{|}X_{0}=x\right]\right)$
$\displaystyle\leq\exp\left[V(x)-(1-\gamma)f(x)+2K\right]$

where $\gamma_{R}=\gamma+K/R$.

###### Proof.

For any $x\in\textbf{X}$, the drift condition implies

$\displaystyle\log[\exp(-V(x))(\mathcal{P}\exp(V))(x)]$
$\displaystyle\leq-(1-\gamma-K/R)f(x)+K(1-f(x)/R)$
$\displaystyle\leq-(1-\gamma_{R})f(x)+KI_{C}(x).$

Define

$M_{n+1}=\exp\left[V(X_{n+1})+(1-\gamma_{R})\sum_{k=1}^{n}f(X_{k})-K\sum_{k=1}^{n}I_{C}(X_{k})\right].$

For every $n\in\mathbb{Z}_{+}$, let $\mathcal{F}_{n}$ be the standard
filtration for $X_{0},X_{1},\ldots,X_{n}$ and then using the drift condition,
$\mathbb{E}_{x}(M_{n+1}|\mathcal{F}_{n})\leq M_{n}$ and $M_{n}$ is a
supermartingale. By the optional sampling theorem [Kallenberg, 2021, Theorem
9.12],

$\displaystyle\mathbb{E}\left(\exp\left[-K+(1-\gamma_{R})\sum_{k=1}^{\tau_{C}\wedge
n}f(X_{k})\bigm{|}X_{0}=x\right]\right)\leq\mathbb{E}(M_{\tau_{C}\wedge
n}\bigm{|}X_{0}=x)$ $\displaystyle\leq\mathbb{E}(M_{2}\bigm{|}X_{0}=x)$
$\displaystyle\leq\mathbb{E}(\exp(V(X_{1}))\bigm{|}X_{0}=x).$

Taking the limit and using Fatou’s lemma [Bogachev and Ruas, 2007, Theorem
2.8.3],

$\displaystyle\mathbb{E}\left(\exp\left[(1-\gamma_{R})\sum_{k=1}^{\tau_{C}}f(X_{k})\bigm{|}X_{0}=x\right]\right)$
$\displaystyle\leq\mathbb{E}\left[\exp(V(X_{1}))+K)\bigm{|}X_{0}=x\right]$
$\displaystyle\leq\exp\left[V(x)-(1-\gamma)f(x)+2K\right].$

∎

We have the identity under the multiplicative drift. The following proof is
similar to Theorem 3.

######  Theorem 4.

Assume the Markov kernel $\mathcal{P}$ has unique invariant measure $\Pi$ and
the multiplicative drift condition (2) holds. If $sup_{x\in C}V(x)<\infty$,
then for every $\varphi:\textbf{X}\to\mathbb{R}$ such that $|\varphi|\leq f$,
then $\int_{\textbf{X}}\varphi d\Pi<\infty$ and

$\displaystyle\int_{\textbf{X}}\mathcal{K}_{C}\varphi
d\Pi=\int_{\textbf{X}}\varphi d\Pi.$

###### Proof of Theorem 4.

Using the multiplicative drift and since $\Pi$ is the unique invariant measure
and under the drift condition, we can show as in Theorem 3 for every simple
function $\varphi:\textbf{X}\to\mathbb{R}$,

$\int_{\textbf{X}}\varphi d\Pi=\int_{C}\mathcal{K}_{C}\varphi d\Pi.$

Now let $\varphi:\textbf{X}\to\mathbb{R}$ such that $|\varphi|\leq f$. Since
$\varphi$ is measurable, we can choose a sequence of simple functions
$\varphi_{n}\to\varphi$ pointwise with $|\varphi_{n}|\leq|\varphi_{n+1}|$.
Using Lemma 3 and $1+x\leq e^{x}$ for $x\geq 0$,

$\displaystyle(1-\gamma_{R})(\mathcal{K}_{C}f)(x)$
$\displaystyle\leq\mathbb{E}\left\\{\exp\left[(1-\gamma_{R})\sum_{k=1}^{\tau_{C}}f(X_{k})\right]\bigm{|}X_{0}=x\right\\}$
$\displaystyle\leq\exp(V(x)+2K).$

Since we have assumed $\sup_{x\in C}V(x)<\infty$, then $\sup_{x\in
C}(\mathcal{K}_{C}f)(x)<\infty$. By dominated convergence [Bogachev and Ruas,
2007, Theorem 2.8.1],

$\displaystyle\int_{\textbf{X}}\mathcal{K}_{C}\varphi
d\Pi=\int_{\textbf{X}}\lim_{n\to\infty}\mathcal{K}_{C}\varphi_{n}(x)\Pi(dx)=\lim_{n\to\infty}\int_{\textbf{X}}\mathcal{K}_{C}\varphi_{n}(x)\Pi(dx)=\int_{\textbf{X}}\varphi
d\Pi.$

∎

###### Proof of Proposition 3.

We have the upper bound

$\displaystyle\left|\sum_{i=1}^{N}\frac{w(Y_{i})}{\sum_{j}w(Y_{j})}(\mathcal{K}_{C}\varphi)(Y_{i})-\sum_{i=1}^{N}\frac{w(Y_{i})}{N}(\mathcal{K}_{C}\varphi)(Y_{i})\right|$
$\displaystyle\leq\max_{i\in
1,\ldots,N}(\mathcal{K}_{C}f)(X_{i})\left|\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})-1\right|.$

This implies the upper bound

$\displaystyle\left|\sum_{i=1}^{N}w^{N}_{Y}(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\mathbb{E}[(\mathcal{K}_{C}\varphi)(Y_{1})w(Y_{1})]\right|$
$\displaystyle\leq\left|\sum_{i=1}^{N}w^{N}_{Y}(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\sum_{i=1}^{N}\frac{w(Y_{i})}{N}(\mathcal{K}_{C}\varphi)(Y_{i})\right|$
$\displaystyle\hskip
10.00002pt+\left|\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\mathbb{E}[(\mathcal{K}_{C}\varphi)(Y_{1})w(Y_{1})]\right|$
$\displaystyle\leq\sup_{x}(\mathcal{K}_{C}f)(x)\left|\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})-1\right|$
$\displaystyle\hskip
10.00002pt+\left|\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\mathbb{E}[(\mathcal{K}_{C}\varphi)(Y_{1})w(Y_{1})]\right|.$

Since for $a,b\geq 0$, $a,b<\epsilon/2$ implies $a+b<\epsilon$, a union
probability bound and Hoeffding’s inequality [Hoeffding, 1963],

$\displaystyle\mathbb{P}\left(\left|\sum_{i=1}^{N}w^{N}_{Y}(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\mathbb{E}[(\mathcal{K}_{C}\varphi)(Y_{1})w(Y_{1})]\right|\geq\epsilon\right)$
$\displaystyle\leq\mathbb{P}\left(\sup_{x}(\mathcal{K}_{C}f)(x)\left|\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})-1\right|\geq\epsilon/2\right)$
$\displaystyle\hskip
10.00002pt+\mathbb{P}\left(\left|\frac{1}{N}\sum_{i=1}^{N}w(Y_{i})(\mathcal{K}_{C}\varphi)(Y_{i})-\mathbb{E}[(\mathcal{K}_{C}\varphi)(Y_{1})w(Y_{1})]\right|\geq\epsilon/2\right)$
$\displaystyle\leq
4\exp\left(-\frac{\epsilon^{2}N}{2\left[\sup_{x}(\mathcal{K}_{C}f)(x)\right]^{2}w_{*}^{2}}\right).$

Since the multiplicative drift condition holds, $\Pi$ is the unique invariant
distribution and $\sup_{x\in C}V(x)<\infty$, then Theorem 4 implies

$\mathbb{E}[(\mathcal{K}_{C}\varphi)(Y_{1})w(Y_{1})]=\int_{C}(\mathcal{K}_{C}\varphi)d\Pi=\int_{\textbf{X}}\varphi
d\Pi.$

By Lemma 3 and using $1+x\leq e^{x}$, then for $x\in C$,

$\displaystyle(1-\gamma_{R})(\mathcal{K}_{C}|\varphi|)(x)$
$\displaystyle\leq(1-\gamma_{R})(\mathcal{K}_{C}f)(x)$
$\displaystyle\leq\exp\left[V(x)-(1-\gamma)f(x)+2K\right].$

∎

###### Proof of Proposition 4.

Let $\lambda>0$ and let $\varphi:\textbf{X}\to\mathbb{R}$ with
$|\varphi|\leq(1-\gamma_{R})f/(1+\lambda)$. Since $x^{2}/2+1\leq e^{x}$ for
$x\geq 0$ and using a second order Taylor expansion of $\exp(\lambda\cdot)$,

$\displaystyle\mathbb{E}\left[\exp\left(\lambda
S^{m}_{\tau_{C}}(\varphi)-\lambda\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right)\bigm{|}Y\right]$
$\displaystyle=1+\lambda^{2}\mathbb{E}\left[\int_{0}^{1}\exp\left[t\lambda\left(S^{m}_{\tau_{C}}(\varphi)-\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right)\right]\left(S^{m}_{\tau_{C}}(\varphi)-\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right)^{2}(1-t)dt\bigm{|}Y\right]$ $\displaystyle\leq
1+\lambda^{2}\mathbb{E}\left[\exp\left((1+\lambda)S^{m}_{\tau_{C}}(|\varphi|)+(1+\lambda)\int\mathcal{K}_{C}|\varphi|d\Pi^{N}_{Y}\right)\bigm{|}Y\right].$
$\displaystyle\leq
1+\lambda^{2}\mathbb{E}\left[\exp\left((1-\gamma_{R})S^{m}_{\tau_{C}}(f)+(1-\gamma_{R})\int\mathcal{K}_{C}fd\Pi^{N}_{Y}\right)\bigm{|}Y\right].$

By Lemma 3, we have the upper bound

$\sup_{x\in
C}\mathbb{E}\left[(1-\gamma_{R})\exp(S^{m}_{\tau_{C}}(f))\Bigm{|}X_{0}=x,Y\right]\leq\exp[B_{R}].$

Since $\exp$ is convex, we have by Jensen’s inequality,

$\exp\left[(1-\gamma_{R})\int\mathcal{K}_{C}fd\Pi^{N}_{Y}\right]\leq\int\mathbb{E}\left\\{\exp\left[(1-\gamma_{R})S_{\tau_{C}}(f)\bigm{|}X_{0}=x,Y\right]\right\\}\Pi^{N}_{Y}(dx)\leq\exp[B_{R}].$

Using this upper bound and $x+1\leq e^{x}$ for $x\geq 0$,

$\displaystyle\mathbb{E}\left[\exp\left(\lambda
S_{\tau_{C}}(\varphi)-\lambda\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right)\bigm{|}Y\right]$ $\displaystyle\leq
1+\lambda^{2}\exp(2B_{R})$
$\displaystyle\leq\exp\left[\lambda^{2}\exp(2B_{R})\right].$

We have shown this is a sub-Gaussian random variable. Since
$S_{\tau_{C}}^{m}(\varphi)|Y$ are assumed independent and
$\mathbb{E}[S_{\tau_{C}}^{m}(\varphi)|Y]=\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}$, taking the Chernoff bound with the optimal
$\lambda=\epsilon\exp(-2B_{R})/2$,

$\displaystyle\mathbb{P}\left[\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\varphi)-\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\geq\epsilon\Bigm{|}Y\right]$
$\displaystyle\leq\exp\left[M\epsilon\left\\{\lambda^{2}\frac{\exp(2B_{R})}{\epsilon}-\lambda\right\\}\right]$
$\displaystyle\leq\exp\left[-\frac{M\epsilon^{2}}{4\exp(2B_{R})}\right].$

Taking the union probability bound and substituting for the function
$\varphi=(1-\gamma_{R})\psi/(1+\lambda)$ with $|\psi|\leq f$,

$\displaystyle\mathbb{P}\left[\left|\frac{1}{M}\sum_{m=1}^{M}S^{m}_{\tau_{C}}(\psi)-\int\mathcal{K}_{C}\psi
d\Pi^{N}_{Y}\right|\geq\epsilon\Bigm{|}Y\right]\leq
2\exp\left[-\frac{M\epsilon^{2}(1-\gamma_{R})^{2}}{4(1+\lambda)^{2}\exp[2B_{R}]}\right].$

Using that $\lambda\leq 1/2$ and taking the iterated expectation completes the
proof. ∎

###### Proof of Theorem 2.

A probability union bound implies

$\displaystyle\mathbb{P}\left(\left|\frac{1}{N}\sum_{m=1}^{M}S_{\tau_{C}}^{m}(\varphi)-\int_{\textbf{X}}\varphi
d\Pi\right|\geq\epsilon\right)$
$\displaystyle\leq\mathbb{P}\left(\left|\frac{1}{N}\sum_{m=1}^{M}S_{\tau_{C}}^{m}(\varphi)-\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}\right|\geq\epsilon/2\right)+\mathbb{P}\left(\left|\int\mathcal{K}_{C}\varphi
d\Pi^{N}_{Y}-\int\varphi d\Pi\right|\geq\epsilon/2\right).$

The proof is completed using Proposition 4, and Proposition 3. ∎

##  Appendix C Supporting technical results for Section 6

###### Proof of Proposition 5.

By (5) the largest eigenvalue of $\nabla^{2}\ell(x)$ uniformly over
$x\in\mathbb{R}^{d}$ is bounded by $\lambda^{*}$. Since $\ell$ is twice
continuously differentiable, using Taylor expansion, we have for all
$x\in\mathbb{R}^{d}$

$\displaystyle f(x)-f(x^{*})$
$\displaystyle=(x-x^{*})^{T}\int_{0}^{1}[\nabla^{2}\ell(x^{*}+t(x-x^{*}))+\Sigma^{-1}](x-x^{*})(1-t)dt$
$\displaystyle\leq\frac{1}{2}(x-x^{*})^{T}[\lambda^{*}+\Sigma^{-1}](x-x^{*}).$

This then implies a lower bound on the normalizing constant using properties
of the determinant

$\displaystyle\int_{\mathbb{R}^{d}}\exp[-(f(x)-f(x^{*})]dx$
$\displaystyle\geq\int_{\mathbb{R}^{d}}\exp\left[-\frac{1}{2}(x-x^{*})^{T}[\lambda^{*}+\Sigma^{-1}](x-x^{*})\right]dx$
$\displaystyle\geq(2\pi)^{d/2}\det(\lambda^{*}I+\Sigma^{-1})^{-1/2}.$

By convexity of $\ell$, Taylor expansion implies for each $x\in\mathbb{R}^{d}$

$f(x)-f(x^{*})\geq\frac{1}{2}(x-x^{*})^{T}\Sigma^{-1}(x-x^{*}).$

Combining these estimates and using properties of the determinant, we have the
upper bound

$\displaystyle\int_{\mathbb{R}^{d}}\frac{\pi(x)}{q(x)}\pi(x)dx$
$\displaystyle=(2\pi)^{d/2}(1/2+h)^{d/2}\det(\Sigma)^{1/2}\frac{1}{Z^{2}}\int_{\mathbb{R}^{d}}\exp\left[\frac{1}{2(1/2+h)}(x-x^{*})^{T}\Sigma^{-1}(x-x^{*})-2f(x)\right]dx$
$\displaystyle\leq\frac{(1/2+h)^{d/2}\det(\lambda^{*}I+\Sigma^{-1})\det(\Sigma)^{1/2}}{(2\pi)^{d/2}}\int_{\mathbb{R}^{d}}\exp\left[-\frac{1}{2}\frac{2h}{1/2+h}(x-x^{*})^{T}\Sigma^{-1}(x-x^{*})\right]dx$
$\displaystyle\leq\det(\lambda^{*}I+\Sigma^{-1})\det(\Sigma)\frac{(1/2+h)^{d}}{(2h)^{d/2}}$
$\displaystyle\leq\det(\lambda^{*}\Sigma+I)\left(\frac{1}{2\sqrt{2h}}+\sqrt{\frac{h}{2}}\right)^{d}.$

∎

###### Proof of Proposition 6.

We look to apply Theorem 1. It is readily seen that the negative log-
likelihood is infinitely differentiable, and for every $v\in\mathbb{R}^{d}$,

$v^{T}\nabla^{2}\ell_{n}(\beta)v=v^{T}X^{T}\text{diag}[\sigma(X\beta)(1-\sigma(X\beta))]Xv\leq\frac{\left\lVert
X^{T}X\right\rVert_{2}}{4}\left\lVert v\right\rVert_{2}^{2}.$

Since the Hessian is symmetric, we have shown that

$\left\lVert\nabla^{2}\ell_{n}(\beta)\right\rVert_{2}\leq\frac{\left\lVert
X^{T}X\right\rVert_{2}}{4}.$

We may then apply the bound in Proposition 5 to the chosen importance sampling
proposal $Q_{h}$.

The proof will be complete if we obtain a drift condition with $V$ defined by
$V(\beta)=1+\left\lVert\beta\right\rVert_{2}^{2}$. Fix
$\beta\in\mathbb{R}^{d}$ and we will take the conditional expectation of the
first step of the Gibbs sampler with $\beta_{1}|\omega_{1}$ and then
$\omega_{1}|\beta$. We have the variance identity

$\displaystyle\mathbb{E}\left(\left\lVert\beta_{1}\right\rVert_{2}^{2}|\omega_{1},\beta\right)$
$\displaystyle=\left\lVert\mu(\omega_{1})\right\rVert_{2}^{2}+\text{tr}(\Sigma(\omega_{1})).$

Let $\Omega_{1}^{1/2}=\text{diag}(\omega^{1/2})$. Since $\Sigma$ is SPD, there
exists a decomposition with a SPD matrix $\Sigma^{-1/2}$ so that the
$\Sigma^{-1}=\Sigma^{-1/2}\Sigma^{-1/2}$. By singular value decomposition
[Horn and Johnson, 2012, Theorem 2.6.3], we may write
$\Omega_{1}^{1/2}X\Sigma^{1/2}=U_{\omega_{1}}D_{\omega_{1}}V_{\omega_{1}}^{T}$
where $V_{\omega_{1}}\in\mathbb{R}^{n\times n}$ is an orthogonal matrix so
that $V_{\omega_{1}}^{T}V_{\omega_{1}}=V_{\omega_{1}}V_{\omega_{1}}^{T}=I_{n}$
and $U_{\omega_{1}}\in\mathbb{R}^{d\times d}$ is also orthogonal and
$D_{\omega_{1}}\in\mathbb{R}^{n\times d}$ is a block diagonal matrix. Then we
have the decomposition

$\displaystyle\Sigma(\omega_{1})$
$\displaystyle=\Sigma^{1/2}([\Omega_{1}^{1/2}X\Sigma^{1/2}]^{T}\Omega_{1}^{1/2}X\Sigma^{1/2}+I)^{-1}\Sigma^{1/2}$
$\displaystyle=\Sigma^{1/2}V_{\omega_{1}}(D_{\omega_{1}}^{T}D_{\omega_{1}}+I)^{-1}V_{\omega_{1}}^{T}\Sigma^{1/2}.$

This decomposiiton implies by the submultiplicative property of the matrix
norm,

$\displaystyle\left\lVert\mu(\omega_{1})\right\rVert_{2}^{2}\leq\left\lVert\Sigma(\omega_{1})\right\rVert_{2}^{2}\left\lVert
X^{T}(Y-(1/2)1_{d})\right\rVert_{2}^{2}\leq\left\lVert\Sigma^{1/2}\right\rVert_{2}^{2}\left\lVert
X^{T}(Y-(1/2)1_{d})\right\rVert_{2}^{2}.$

Using cyclic property of the trace, we also have

$\displaystyle\text{tr}\left[\Sigma(\omega_{1})\right]\leq\text{tr}\left[V_{\omega_{1}}^{T}\Sigma
V_{\omega_{1}}(D_{\omega_{1}}^{T}D_{\omega_{1}}+I)^{-1}\right]\leq\text{tr}\left[V_{\omega_{1}}^{T}\Sigma
V_{\omega_{1}}\right]\leq\text{tr}\left[V_{\omega_{1}}V_{\omega_{1}}^{T}\Sigma\right]\leq\text{tr}\left[\Sigma\right].$

Combining these estimates and taking the iterated expectation

$\displaystyle\mathbb{E}\left(\left\lVert\beta_{1}\right\rVert_{2}^{2}+1|\beta\right)$
$\displaystyle\leq 1+\left\lVert\Sigma\right\rVert_{2}^{2}\left\lVert
X^{T}(Y-1/21_{d})\right\rVert_{2}^{2}+\text{tr}(\Sigma).$

.

This drift condition combined with a minorization shows the posterior is the
unique invariant measure [Hairer, 2021] or Harris recurrence would work as
well [Meyn and Tweedie, 2009, Theorem 10.4.4]. ∎

## References

  * Agapiou et al. [2017] S. Agapiou, O. Papaspiliopoulos, D. Sanz-Alonso, and A. M. Stuart.  Importance sampling: intrinsic dimension and computational cost.  _Statistical Science_ , 32(3):405–431, 2017.
  * Bogachev and Ruas [2007] Vladimir Igorevich Bogachev and Maria Aparecida Soares Ruas.  _Measure theory_ , volume 1.  Springer, 2007.
  * Choi and Hobert [2013] Hee Min Choi and James P. Hobert.  The Polya-Gamma Gibbs sampler for Bayesian logistic regression is uniformly ergodic.  _Electronic Journal of Statistics_ , 7(none):2054–2064, 2013.
  * Craiu and Meng [2005] Radu V. Craiu and Xiao-Li Meng.  Multiprocess parallel antithetic coupling for backward and forward markov chain monte carlo.  _The Annals of Statistics_ , 33(2):661–697, 2005\.
  * Craiu et al. [2009] Radu V. Craiu, Jeffrey Rosenthal, and Chao Yang.  Learn from thy neighbor: Parallel-chain and regional adaptive mcmc.  _Journal of the American Statistical Association_ , 104(488):1454–1466, 2009.
  * Detrano et al. [1989] Robert Detrano, Andras Janosi, Walter Steinbrunn, Matthias Pfisterer, Johann-Jakob Schmid, Sarbjit Sandhu, Kern H. Guppy, Stella Lee, and Victor Froelicher.  International application of a new probability algorithm for the diagnosis of coronary artery disease.  _The American Journal of Cardiology_ , 64(5):304–310, 1989.  ISSN 0002-9149.
  * Devraj et al. [2020] Adithya Devraj, Ioannis Kontoyiannis, and Sean Meyn.  Geometric ergodicity in a weighted Sobolev space.  _The Annals of Probability_ , 48(1):380–403, 2020\.
  * Dudley [1969] R. M. Dudley.  The speed of mean glivenko-cantelli convergence.  _The Annals of Mathematical Statistics_ , 40(1):40–50, 1969.
  * Gelman and Rubin [1992] Andrew Gelman and Donald B. Rubin.  Inference from Iterative Simulation Using Multiple Sequences.  _Statistical Science_ , 7(4):457–472, 1992.
  * Geyer [1992] Charles J. Geyer.  Practical Markov Chain Monte Carlo.  _Statistical Science_ , 7(4):473–483, 1992.
  * Glynn and Rhee [2014] Peter W. Glynn and Chang-Han Rhee.  Exact estimation for markov chain equilibrium expectations.  _Journal of Applied Probability_ , 51a:377–389, 2014.
  * Hairer [2021] Martin Hairer.  Convergence of markov processes, 2021.
  * Hairer and Mattingly [2011] Martin Hairer and Jonathan C. Mattingly.  Yet another look at Harris’ ergodic theorem for Markov chains.  _Sem. on Stoch. Anal., Random Fields and Appl. VI_ , 63, 2011.
  * Hoeffding [1963] Wassily Hoeffding.  Probability inequalities for sums of bounded random variables.  _Journal of the American Statistical Association_ , 58(301):13–30, 1963.
  * Horn and Johnson [2012] Roger A. Horn and Charles R. Johnson.  _Matrix Analysis_.  Cambridge University Press, 2012.
  * Jacob et al. [2020a] Pierre E. Jacob, John O’Leary, and Yves F. Atchadé.  Unbiased Markov Chain Monte Carlo Methods with Couplings.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(3):543–600, 05 2020a.
  * Jacob et al. [2020b] Pierre E. Jacob, John O’Leary, and Yves F. Atchadé.  Unbiased Markov Chain Monte Carlo Methods with Couplings.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(3):543–600, 05 2020b.
  * Kallenberg [2021] Olav Kallenberg.  _Foundations of Modern Probability_.  Springer, Cham, 3 edition, 2021.
  * Kontoyiannis and Meyn [2005] Ioannis Kontoyiannis and Sean Meyn.  Large Deviations Asymptotics and the Spectral Theory of Multiplicatively Regular Markov Processes.  _Electronic Journal of Probability_ , 10(none):61–123, 2005.
  * Lao et al. [2020] Junpeng Lao, Christopher Suter, Ian Langmore, Cyril Chimisov, Ashish Saxena, Pavel Sountsov, Dave Moore, Rif A. Saurous, Matthew D. Hoffman, and Joshua V. Dillon.  tfp.mcmc: Modern markov chain monte carlo tools built for modern hardware.  _arXiv preprint_ , 2020.
  * Margossian et al. [2023] Charles C. Margossian, Matthew D. Hoffman, Pavel Sountsov, Lionel Riou-Durand, Aki Vehtari, and Andrew Gelman.  Nested $\hat{r}$: Assessing the convergence of markov chain monte carlo when running many short chains.  _arXiv preprint_ , 2023.
  * Meyn and Tweedie [2009] Sean P. Meyn and Richard L. Tweedie.  _Markov Chains and Stochastic Stability_.  Cambridge University Press, Cambridge, 2 edition, 2009.  ISBN 0521731828.
  * Nicholas G. Polson and Windle [2013] James G. Scott Nicholas G. Polson and Jesse Windle.  Bayesian inference for logistic models using pólya–gamma latent variables.  _Journal of the American Statistical Association_ , 108(504):1339–1349, 2013.
  * Nummelin and Arjas [1976] Esa Nummelin and Elja Arjas.  A direct construction of the $r$-invariant measure for a markov chain on a general state space.  _The Annals of Probability_ , 4(4):674–679, 1976\.
  * Roberts et al. [1997] Gareth O. Roberts, A. Gelman, and W. R. Gilks.  Weak convergence and optimal scaling of random walk metropolis algorithms.  _Ann. Appl. Probab._ , 7(1):110–120, 1997.
  * Wang et al. [2021] Guanyang Wang, John O’Leary, and Pierre Jacob.  Maximal couplings of the metropolis-hastings algorithm.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1225–1233. PMLR, 2021.