Uniform Ergodicity of the Iterated Conditional SMC and Geometric Ergodicity of Particle Gibbs samplers

Introduction

Particle Markov chain Monte Carlo (P-MCMC) methods are a set of recently proposed sampling techniques particularly well suited to the Bayesian estimation of static parameters in general state-space models , although their scope extends beyond this class of models. At an abstract level, once the likelihood function and prior are defined, inference for this class of models relies on a probability distribution \pi\big{(}{\rm d}\theta\times{\rm d}x\big{)}, defined on some measurable space $\left(\Theta\times\mathsf{X},\mathcal{B}(\Theta)\times\mathcal{B}\mathsf{(X)}\right)$, where $\theta$ is generally a low dimensional static parameter, the static parameter, while $x$, the hidden state of the system, is a large vector with a non-trivial dependence structure. Here, $\mathcal{B}(\,\cdot\,)$ denotes the $\sigma$-algebra related to the corresponding space. In practice the complexity of such probability distributions requires the use of sampling techniques to effectively carry out inference. When $\theta$ is known sequential Monte Carlo methods (SMC), or particle filters, are particularly suitable to carry out inference about $x$ by approximately sampling from the conditional distribution \pi_{\theta}\big{(}{\rm d}x\big{)}. These algorithms rely on interacting particle systems and their performance and accuracy can be improved by increasing the number $N$ of such particles. P-MCMC realises the synthesis between SMC methods and classical Markov chain Monte Carlo (MCMC) methods, that is it allows the construction of Markov transition probabilities leaving \pi\big{(}{\rm d}\theta\times{\rm d}x\big{)} at least marginally invariant and from which it is possible to sample realisations $\{(\theta_{i},X_{i}),i\geq 0\}$ with attractive efficiency properties.

The particle marginal Metropolis–Hastings (PMMH) method is one such algorithm, which takes advantage of the availability of unbiased estimators of the likelihood function to provide an exact approximation of an idealized algorithm which computes the likelihood function exactly. The algorithm simply consists of replacing the true value of the likelihood function required to implement the standard Metropolis–Hastings (MH) algorithm with estimators, but is nevertheless guaranteed to be correct in that it leaves the required distribution of interest marginally invariant. In PMMH, the estimator of the likelihood is a byproduct of a sequential Monte Carlo (SMC) algorithm, whose accuracy can be improved by increasing $N$.

In contrast, the particle Gibbs (PGibbs) sampler involves approximating a Gibbs sampler which consists of constructing a Markov chain $\{(\theta_{i},X_{i}),i\geq 0\}$, by repeatedly sampling from \pi_{\theta}\big{(}{\rm d}x\big{)} and \pi_{x}\big{(}{\rm d}\theta\big{)} in turn. In practice sampling from \pi_{\theta}\big{(}{\rm d}x\big{)} may be particularly difficult and the conditional SMC (cSMC) update is a Markov transition probability $P_{N,\theta}$ which leaves \pi_{\theta}\big{(}{\rm d}x\big{)} invariant, therefore allowing the implementation of a Metropolis-within-Gibbs algorithm, that is a Markov transition probability leaving \pi\big{(}{\rm d}\theta\times{\rm d}x\big{)} invariant. The cSMC relies for its construction, as suggested by its name, on an SMC-like procedure and it is expected that as $N$ increases $P_{N,\theta}$ approaches \pi_{\theta}\big{(}{\rm d}x\big{)}.

While PMMH methods have been studied in a series of papers , a theoretical study of the PGibbs is still missing. Indeed it has been shown that as $N$ increases, performance of the PMMH approaches that of the exact MH algorithm but the question of the approximation of the Gibbs sampler by a PGibbs has not been addressed to date. We note however that a study of one of its components, the cSMC update, has recently been undertaken in , in which a coupling argument is central to their analysis. We refer to the Markov chain obtained by iterating the cSMC algorithm for a fixed target distribution as iterated i-cSMC here in order to distinguish it from that of the PGibbs. The present manuscript addresses questions concerning the i-cSMC similar to those of , but our results differ in many respects and complement their findings in several directions. At a technical level our approach seems to be more straightforward in the scenario considered, relies on weaker assumptions for uniform convergence which we prove are necessary and sufficient and lead to quantitative bounds on performance measures in terms of the number $N$ of particles involved. We additionally transfer sufficient conditions for uniform ergodicity of the i-cSMC Markov chain into sufficient conditions for geometric ergodicity of the associated PGibbs Markov chain, the main motivation behind our work. This allows us in particular to show that under some conditions PGibbs is asymptotically as efficient as the Gibbs sampler as the number $N$ of particles increases.

Contemporary to the first version of the present manuscript , have also provided essentially the same sufficient conditions for the uniform convergence of the i-cSMC Markov chain (Theorem 1, Section 3) using a different proof technique. Here we have further established that the aforementioned conditions are also necessary for uniform convergence in general, but also geometric ergodicity in many realistic scenarios (Section 6). Similarly to us also provide quantitative bounds and associated scaling properties of the i-cSMC, albeit for a different set of specialised conditions (a detailed comparison of the assumptions is provided after Theorem 5 at the end of Section 3). We have also very recently become aware of the contribution to the analysis of the properties of the cSMC, established using the formalism of , but their practical implications are unclear. Similarly to , do not attempt to address the practically important question of how uniform ergodicity of the i-cSMC can be translated into geometric ergodicity of the PGibbs sampler, an issue we address in Section 7. In Section 8 we contrast the results obtained in this paper concerning the i-cSMC and PGibbs algorithm with known results concerned with other particle MCMC methods and draw final conclusions.

Similarly to SMC methods, the cSMC and associated algorithms are complex mathematical objects which require the introduction of sometimes overwhelming notation which may obscure the main ideas. In the next section we attempt to remedy this by presenting our results in a simplified scenario, which captures our main ideas, before moving on to the general scenario.

Statement of our results in a simplified scenario

We first explain our results on a particularly simple instance of the i-cSMC algorithm. This should provide the reader with the essence of the results proved later on in the general scenario, while its simple structure will allow us to outline the main idea behind our proof in the general set-up (in Section 4).

with $G(x):=\pi({\rm d}x)/M({\rm d}x)$ and the convention $z^{1}=x$. Our first results are concerned with properties of the homogeneous Markov chain with transition probability $P_{N}$, in terms of $\bar{G}:=\pi-{\rm ess}\sup_{x}G(x)$ and $N$. We refer to the resulting algorithm as iterated SIR (i-SIR).

We briefly introduce notions that allow us to make quantitative statements about the Markov chains under study. We use classical Hilbert space techniques for the analysis of reversible Markov chains. Letting \mu\bigl{(}\cdot\bigr{)} be a probability distribution defined on some measurable space \bigl{(}\mathsf{E},\mathcal{B}\bigl{(}\mathsf{E}\bigr{)}\bigr{)}, we define the function space

We denote \nu\Pi^{k}f:=\nu\bigl{(}\Pi^{k}f\bigr{)} and refer to $\nu\Pi^{k}$ as either a probability measure or its corresponding operator on $L^{2}(\mathsf{E},\mu)$. For f\in L^{2}\bigl{(}\mathsf{E},\mu\bigr{)}, we define the variance of $f$ under $\mu$ as ${\rm var}_{\mu}(f):=\mu(f^{2})-\mu(f)^{2}$ and the “asymptotic variance” of M^{-1}\sum_{i=1}^{M}f\big{(}\xi_{i}\big{)} for stationary realizations $\{\xi_{i},i\geq 0\}$ associated to the homogeneous Markov chain with transition $\Pi$ as

Some of our results involve norms of signed measures. As in, e.g., , for any signed measure $\nu$ on \bigl{(}\mathsf{E},\mathcal{B}\bigl{(}\mathsf{E}\bigr{)}\bigr{)} we let

denote the total variation distance and for $\nu\ll\mu$,

$P_{N}$ is reversible with respect to $\pi$ and positive, that is the i-SIR Markov chain has non-negative stationary autocorrelations.

If $\bar{G}<\infty$, and $N\geq 2$, the i-SIR Markov chain is uniformly ergodic with for any $x\in\mathsf{X}$,

If $\bar{G}<\infty$, then for any $f\in L^{2}(\mathsf{X},\pi)$,

If $\bar{G}=\infty$ then the i-SIR Markov chain cannot be geometrically ergodic for any finite $N$.

The second and third points provide quantitative bounds on standard measures of performance for MCMC algorithms, where the second provides a bound on the uniform (or equivalently uniformly geometric) rate of convergence of the Markov chain. Interest in algorithms such as i-SIR is motivated empirically from observed behaviour in line with the above bounds, as performance improves as $N$ increases, and part of our purpose here is to confirm and quantify theoretically such empirical successes. Moreover, this improvement can often be obtained with little extra computational effort, since on a parallel architecture one can sample from $M$ and evaluate $G$ in parallel, a characteristic of SMC algorithms more generally .

While i-SIR can be used alone to sample from fairly general distributions, it can also be used as a constituent element of more elaborate MCMC schemes. Assume now that we wish to sample from a distribution $\pi$ defined on some measurable space \left(\Theta\times\mathsf{X},\mathcal{B}(\Theta)\times\mathcal{B}\bigl{(}\mathsf{X}\bigr{)}\right), often defined for some $S\in\mathcal{B}(\Theta)\times\mathcal{B}(\mathsf{X})$ via (note the different nature of $\pi$ as compared to earlier)

where $\{G_{\theta},\theta\in\Theta\}$ is a collection of non-negative potential functions and $\{M_{\theta},\theta\in\Theta\}$ a collection of probability measures which define for each $\theta\in\Theta$ the conditional distributions \pi_{\theta}\bigl{(}{\rm d}x\bigr{)}:=M_{\theta}\bigl{(}{\rm d}x\bigr{)}G_{\theta}(x)/\gamma_{\theta} with

The interpretation in a statistical context is that $\varpi$ is the prior distribution for some parameter $\theta$ of interest, whilst $\gamma_{\theta}$ is the likelihood function associated with some observed data and $x$ corresponds to the so-called latent variable(s). The form of $\gamma_{\theta}$ is often derived from the data being explained by the latent variable $x$ whose a priori distribution conditional upon $\theta$ is $M_{\theta}$ and the likelihood function given the data and $x$ is $G_{\theta}(x)$. Assume here that we are able to sample from $\pi_{x}$, the conditional distribution of $\theta$ given $X=x$. For any $\theta\in\Theta$ one can define the i-SIR kernel for any (x,S)\in\mathsf{X}\times\mathcal{B}\bigl{(}\mathsf{X}\bigr{)} via

with $z^{1}=x$, so that the invariant distribution associated with $P_{N,\theta}$ is $\pi_{\theta}$, the conditional distribution of $X$ given $\theta$. One can sample from $\pi({\rm d}\theta\times{\rm d}x)$ with the following Markov transition, defined for any \bigl{(}\theta_{0},x,S\bigr{)}\in\Theta\times\mathsf{X}\times\big{(}\mathcal{B}(\Theta)\times\mathcal{B}(\mathsf{X})\big{)} via

which can be viewed as an exact approximation of the Gibbs sampler defined via

Assume the $\Gamma$ Markov chain is such that there exists $\beta\in(0,1]$ such that for any $f:\mathsf{X}\rightarrow$ and $\nu\ll\pi$

If $\bar{G}<\infty$, and $N\geq 2$, then for any $f:\mathsf{X}\rightarrow$ and $\nu\ll\pi$

For any $f\in L^{2}(\mathsf{X},\pi)$ and $N\geq 2$, the asymptotic variance ${\rm var}(f,\Phi_{N})$ satisfies

For any f\in L^{2}\bigl{(}\Theta,\pi\bigr{)} and $N\geq 2$, the asymptotic variance ${\rm var}(f,\Phi_{N})$ satisfies

In the sequel, we prove similar results in the more general (and complex) scenario where $P_{N,\theta}$ is defined by a general cSMC algorithm with multinomial resampling, but the key ideas and results are similar (Section 3). The results concerning the general form of the PGibbs sampler, from which its convergence in the sense of points 1–3 above follows, can be found in Section 7.

The i-cSMC and its properties

and can define for any $S\in\mathcal{B}(\mathsf{X})$ the probability distribution $\pi$ (which will be the target distribution of interest)

Note in particular that with the convention above, for any $l\geq 2$ and $z_{0}\in\mathsf{Z}$, M_{0,l}\bigl{(}z_{0},{\rm d}z_{1:l}\bigr{)}:=M_{1}({\rm d}z_{1})\times M_{1,l}\bigl{(}z_{1},{\rm d}z_{2:l}\bigr{)}.

where we keep $k$ to emphasize that we are sampling from that mixture. For the last iteration we only require one index and point out that whereas $A_{t}\in[N]^{N}$ for $t=1,\ldots,T-1$, we have $A_{T}\in[N]$ following

For any \mathbf{i}:=\bigl{(}i_{1},i_{2},\ldots,i_{T}\bigr{)}\in[N]^{T}, z_{1:T}\in\bigl{(}\mathsf{Z}^{N}\bigr{)}^{T}, a_{1:T}:=(a_{1},\ldots,a_{T})\in\bigl{(}[N]^{N}\bigr{)}^{T-1}\times[N] and S\in\mathcal{B}\bigl{(}\mathsf{X}\bigr{)} define

Then the transition kernel of the iterated conditional SMC (i-cSMC), in the multinomial sampling scenario, is given for any $x\in\mathsf{X}$ and $S\in\mathcal{B}\left(\mathsf{X}\right)$ by

that is, conditional upon $x$ we consider the probability distribution of those trajectories $Z_{1:T}^{\mathbf{i}}$ generated by the cSMC which form a lineage compatible with the lineages defined by the random variables $A_{1:T}$. Our main results concerning the i-cSMC algorithm are the following (our results concerning the particle Gibbs sampler are provided in Section 7). We will denote by $\pi_{t}$ the corresponding marginal distribution of $\pi$ (see (8) for a precise definition).

For $N\geq 2$ the i-cSMC algorithm with kernel $P_{N}$

is reversible with respect to $\pi$ and defines a positive operator,

if for all $t\in\{1,\ldots,T\}$ $\pi_{t}-{\rm ess}\sup_{z_{t}}G_{t}(z_{t})<\infty$ then there exists $\epsilon_{N}>0$ such that

for any (x,S)\in\mathsf{X}\times\mathcal{B}\bigl{(}\mathsf{X}\bigr{)},

for any probability distribution $\nu\ll\pi$ on \bigl{(}\mathsf{X},\mathcal{B}(\mathsf{X})\bigr{)} and $k\geq 1$

for any f\in L^{2}\bigl{(}\mathsf{X},\pi\bigr{)}

From Lemma 24, statement (d) holds under a more abstract assumption, but we have chosen this explicit simplified statement for clarity at this point. In fact we suspect that (d) holds under the assumption $\pi_{t}-{\rm ess}\sup_{z_{t}}G_{t}(z_{t})=\infty$ for some $t\in[T]$ only, that is essential boundedness is a necessary condition for geometric ergodicity; see Conjecture 26.

and \bar{M}_{p,p+1}\bigl{(}z,\cdot\bigr{)}=M_{p+1}\bigl{(}z,\cdot\bigr{)} and for $q>p\geq 0$ we have the recursive definition, for any $z_{p}\in\mathsf{Z}$,

The first condition is rather abstract, and can be viewed as a condition on the $h$-functions investigated in in the context of stability properties of standard SMC algorithms.

One can however show that (A3) is implied by the following stronger assumption (see Lemma 18).

There exists a constant $1\leq\beta<\infty$ such that for any $p\geq 1$ and any $(z,z^{\prime})\in\mathsf{Z}$ and $S\in\mathcal{B}(\mathsf{Z})$,

The potential functions $G_{p}$ satisfy, for some $\delta<\infty$,

Similar results for the PGibbs sampler are provided in Section 7.

The proofs of the various results are the subject of the following sections. More specifically, statement

follows from Lemma 10 (the latter property was established in and the former noted/proved in ),

all parts follow from Corollary 14 and [2, Proposition 33], which gathers generic results on $\pi-$invariant Markov chains satisfying (b)(b)(i),

follows from Proposition 22 and Lemma 24; Remark 25.

Follows from Proposition 15, Corollary 16 and Lemma 18. ∎

As pointed out in the introduction, soon after completing this work we have become aware of , where a subset of our results have also been independently discovered. This motivates the following comparison. Result (b)(b)(i) of Theorem 1 is identical to Theorem 1 of , but relies on a different proof. Results (b)(b)(ii)–(b)(iv) rely on standard arguments, although (b)(iv) does not seem to be well known and establishes informative quantitative bounds. The study of the necessity of our conditions to imply uniform or geometric ergodicity is not addressed in . The result of Theorem 5 corresponds to Proposition 5 of . The conditions under which Theorem 5 holds are rather stringent for some applications, in particular in the state-space model scenario. As discussed by in that scenario (A4) will essentially only hold in the case where $\mathsf{X}$ is compact. The condition (A3) is weaker and more natural in our analysis, but is not currently easy to verify in applications except through (A4).

In an attempt to relax (A4), the authors of investigate another set of specialised assumptions guaranteeing that the result of Theorem 5 holds even in some non-compact scenarios provided the number of particles $N$ grows at a rate $T^{1/\gamma}$ for any $\gamma\in(0,1)$, a result in line with what is obtained with the stronger assumption (A4), for which $\gamma=1$ is permissible. This requires the specification of a “moment assumption” which aims at controlling the variations of the various quantities involved under the law of the observation process $\{Y_{t},t\geq 0\}$. Their approach, however, does not seem to allow one to consider the scaling properties of the PGibbs sampler (i.e. not just the i-cSMC); see their Theorem 6 and Remark 7. More importantly we note that their results require the law of the data to coincide with that of the specified model for some $\theta^{\star}\in\Theta$ which, although suggestive of what may happen in practice, is always an idealization. This delicate work is the main focus of the remainder of their investigation while here, in addition to establishing the necessity of some of the conditions, we have focused on the transference of the results obtained for the i-cSMC to the PGibbs sampler (Section 7) with the aim of showing that the PGibbs has performance inferior to that of the Gibbs sampler, but arbitrarily close if we increase $N$.

Establishing the uniform minorization condition

Before proceeding we turn to the i-SIR which is particularly simple to analyze. The reason for detailing the short analysis of this simple scenario is to provide the reader with an overview of the developments which are to follow – the remainder of the paper essentially replicates the key steps of the argument below, albeit in the more complex SMC framework. Notice that in this scenario $\mathsf{X}=\mathsf{Z}$ since $T=1$. We let $G(x):=\pi({\rm d}x)/M({\rm d}x)$ for any $x\in\mathsf{X}$ and assume that $\bar{G}:=\sup_{x\in\mathsf{X}}G(x)<\infty$. Then for (x,S)\in\mathsf{X}\times\mathcal{B}\bigl{(}\mathsf{X}\bigr{)} we can rewrite

This is a uniform minorization condition which immediately implies uniform geometric convergence (see the outline of our results in Section 1), but in the present situation the result is even stronger in that, in particular, it provides us with quantitative bounds on the dependence of the performance of the algorithm on $N$. Indeed it is a standard result that the minorization constant

provides the upper bound $1-\epsilon_{N}$ on the (geometric) rate of convergence of the algorithm, which here vanishes at an asymptotic rate $N^{-1}$ as $N$ increases. As we shall see the fact that the minorization measure is the invariant distribution leads to a direct lower bound on associated Dirichlet forms associated to $P_{N}$ which in turn provide quantitative bounds on the spectral gap and the associated asymptotic variance. In the remainder of the section we generalize the representation of $P_{N}$ in terms of the c2SMC algorithm and “the estimator of the normalizing constant” which suggests applying Jensen’s inequality as above. This requires us to consider estimates of the resulting expectation in Section 5.

In order to proceed further it is required to define the c2SMC process, which is essentially similar to the cSMC process but where conditioning is now upon two trajectories $x,y\in\mathsf{X}$. The definition is therefore similar, but for reasons which will become clearer below the second fixed trajectory is set to have a lineage of the general form $\mathbf{k}:=k_{1:T}\in[N]^{T}$. We will use below the convention that \delta_{a,b}\bigl{(}{\rm d}z^{1}\times{\rm d}z^{k}\bigr{)} reduces to $\delta_{a}({\rm d}z^{1})$ whenever $k=1$. The definition of this process is similar to that of the cSMC algorithm and the distributions involved are defined for $x,y\in\mathsf{X}$ and $\mathbf{k}\in[N]^{T}$ as follows

and for $t=2,\ldots,T-1$ (with the convention $a_{t-1}^{k,l}:=(a_{t-1}^{k},a_{t-1}^{l})$)

For $\mathbf{i}\in\{2,\ldots,N\}^{T}$ and $x\in\mathsf{X}$,

As we shall see the concentration properties of the “estimator of the normalizing constant” plays a central role for any z_{1:T}\in\bigl{(}\mathsf{Z}^{N}\bigr{)}^{T}

We first obtain a uniform minorization condition for the cSMC transition probability. This simple result establishes the expectation of \hat{\gamma}_{T}^{N}\bigl{(}Z_{1:T}\bigr{)} with respect to a c2SMC algorithm as a key quantity of interest, and motivates the non-asymptotic analysis and bounds of Section 5.

For any (x,S)\in\mathsf{X}\times\mathcal{B}\bigl{(}\mathsf{X}\bigr{)} and $N\geq 2$ we have

Using (7), we only keep the trajectories for which there is no coalescence with the first trajectory, i.e., we exclude terms such that $i_{t}=1$ for some $t\in[T]$ and obtain

then for any (x,S)\in\mathsf{X}\times\mathcal{B}\bigl{(}\mathsf{X}\bigr{)}, $P_{N}(x,S)\geq\epsilon_{N}\pi(S)$ and from Proposition 8 all the properties of [2, Proposition 33] apply to the i-cSMC with $\varepsilon=\epsilon_{N}$.

Before proceeding to novel analysis, for completeness we gather two known properties of the i-cSMC (in the general set-up) in the following lemma which will be exploited throughout the remainder of the paper. Both results are immediate upon noticing that the i-cSMC is a two stage Gibbs sampler on an artificial joint distribution (see (14) in [2, Appendix B], which is a generalization of (1)). The results have also been shown in detail in . A proof is included in [2, Appendix B] for completeness.

$P_{N}$, viewed as an operator on L^{2}\bigl{(}\mathsf{X},\pi\bigr{)}, is self-adjoint and positive.

Quantitative bounds for the doubly conditional i-cSMC expectation

Let $x,y\in\mathsf{X}$ and $N\geq 2$. Then,

While the expectation of interest here has been hitherto uninvestigated, the form of Proposition 11 is reminiscent of non-asymptotic results in , in which second moments of \hat{\gamma}_{T}^{N}\bigl{(}Z_{1:T}\bigr{)} are analyzed with respect to the law of a standard SMC algorithm.

We now turn to estimates of the expectation above, starting with very minimal assumptions which allow us to establish the minorization condition required to apply [2, Proposition 33] and deduce most of our results, without the need for assumptions on the dynamic of the system—the number of particles is however required to grow exponentially in order to maintain a set level of performance. We show subsequently that with stronger assumptions on $\{M_{t},G_{t}\}_{t=1}^{T}$ it is possible to show that $N$ should grow linearly with $T$ to ensure that a set level of performance is maintained.

Assume that for all $t\in\{1,\ldots,T\}$, $\bar{G}_{t}:=\sup_{z\in\mathsf{Z}}G_{t}(z)<\infty$, then for any $N\geq 2$

Propositions 8 and 13 together imply that for any $x,S\in\mathsf{X}\times\mathcal{B}(\mathsf{X})$,

and $\lim_{N\rightarrow\infty}\epsilon_{N}=1$.

It should be clear that despite Corollary 14, the term $\prod_{t=1}^{T}\bar{G}_{t}/\gamma_{T}$ typically grows exponentially fast with $T$ whenever the potentials are not constant functions. Therefore, Proposition 13 suggests that the number of particles $N$ should grow exponentially with $T$ in general. However, stronger assumptions on the system under consideration will allow us to maintain a given lower bound on $\epsilon_{N}$ by increasing $N$ only linearly with $T$. We first state our main result using the abstract condition (A3) and then show that classical strong mixing conditions (A4) imply (A3).

First notice that for any $1\leq k\leq n$

and therefore for any $s\in\{1,\ldots,T\}$ and $0<i_{1}<\cdots<i_{s-1}<i_{s}=T+1$ with the notation defined earlier,

and we conclude by an application of the binomial theorem.∎

Propositions 8 and 13 together imply that for any $(x,S)\in\mathsf{X}\times\mathcal{B}(\mathsf{X})$,

Now, let $N-1\geq CT$ for some $C>0$. Then $\epsilon_{N}\geq\exp\left(-\frac{2\alpha-1}{C}\right)$.

Propositions 8 and 15 together imply that

Since $(N-1)\geq CT$ for some $C>0$, and $\log(1+x)\leq x$ for all $x\geq 0$,

The combination of the upper bound of ${\rm var}(f,P_{N})$ in Theorem 1 with Corollary 16 suggests a rough rule of thumb to select $N$ for the i-cSMC Markov kernel. In particular, there is generally a tradeoff between iterating a less computationally intensive Markov kernel more times and iterating a more computationally intensive expensive fewer times. This suggests that one should minimize the function $f(N):=N{\rm var}(f,P_{N})$. While an analytic expression for ${\rm var}(f,P_{N})$ is not available we can minimize its upper bound

with respect to $C$. Assuming that we are in the scenario where $N\gg 1$ and therefore $CT+1\approx CT$ one then finds the unique minimum

(where ${\rm Lambert_{W}}$ is the principal branch of the Lambert W function) or correspondingly

Hence, under (A3) it is only required for $N$ to scale linearly with $T$ in order to maintain a non-vanishing ergodicity rate. Following, e.g., we make the following assumptions on $\{M_{t}\}$ and the potentials $\{G_{t}\}$ which combined define an $m$-step “strong mixing” condition which automatically implies (A3). The following result relies on classical arguments [9, 7, Lemma 4.3]

Necessity of the boundedness assumption and a conjecture

Proposition 13 showed that the i-cSMC kernel is uniformly ergodic if the potentials are bounded. We study here the opposite case, where at least one of the potentials is unbounded. We discover that then the algorithm cannot be uniformly ergodic (Proposition 19), and in many cases the algorithm cannot be geometrically ergodic (Proposition 22 and Lemma 24; Remark 25). We believe that the latter holds in general (Conjecture 26), but a proof has remained elusive. This dichotomy of algorithms which are uniformly ergodic and sub-geometrically ergodic would be in perfect analogy with the behaviour of the independent Metropolis–Hastings [20, Theorem 2.1].

We will denote hereafter the marginal densities of $\pi$ by

where $1\leq t\leq u\leq T$ and we use the shorthand $\pi_{t}(A):=\pi_{t:t}(A)$.

In this section, we will assume that $\mathsf{S}\in\mathcal{B}(\mathsf{Z})^{T}$ is a fixed set such that for all $x\in\mathsf{S}$, $\prod_{t=1}^{T}G_{t}(x_{t})>0$ and $\pi(\mathsf{S})=1$. Further, $\mathsf{S}$ contains all possible starting points of the algorithm, that is, we assume that the state space of the i-cSMC is $\mathsf{S}$. In the discrete case, the minimal $\mathsf{S}$ consists of the points of positive $\pi$-measure, and in the continuous case where $\pi$ admits a density, the set $\mathsf{S}$ can be taken as the set where the density is positive.

Further, we will assume that $\pi_{1}$ is not concentrated on a single point. We can do this without loss of generality, because if $\pi_{1},\ldots,\pi_{t}$ were concentrated on single points of the state space, the algorithm would be deterministic until $\pi_{t+1}$ and we could consider the i-cSMC for $\pi^{\prime}=\pi_{t+1:T}$.

If the i-cSMC kernel is uniformly ergodic, then there exist $K<\infty$ and $\rho\in(0,1)$ such that

Denote the level set $L_{t}(\underline{G}):=\{x_{t}\in\mathsf{Z}\,:\,G_{t}(x_{t})\leq\underline{G}\}$. Lemma 20 shows that there exists $c_{2}=c_{2}(N)\in[1,\infty)$ such that for $G_{t}(x_{t})\geq\underline{G}$

We may estimate for any $i\in[n]$ and all $x\in\mathsf{S}$ such that $G_{t}(x_{t})\geq\delta^{i-1}G_{*}$,

We conclude that for $x\in\mathsf{S}$ such that $G_{t}(x_{t})\geq G_{*}$,

This proves the claim, as $\epsilon>0$ was arbitrary.∎

$P_{N}(x,\{x_{1}\}^{\complement}\times\mathsf{Z}^{T-1})\leq\phi(G(x_{t}))$,

P_{N}(x,\mathsf{Z}^{t-1}\times L_{t}(\underline{G})\times\mathsf{Z}^{T-t})\leq(N-1)^{2}\underline{G}/G_{t}(x_{t})\quad\text{wheneverG_{t}(x_{t})\geq\underline{G}}.

In both cases, we consider the case $t<T$; the special case $t=T$ can be treated similarly. In order to facilitate the theoretical analysis, we introduce a non-standard implementation of the cSMC which relies on the remark that at any time instant a given particle can only have a maximum number $N$ of children. Hence when implementing the cSMC it is always possible to draw $N$ children first and then decide who is carried forward according to the standard selection mechanism. It is in fact possible to push this idea further and, given a fixed $x\in\mathsf{S}$, to sample the following $N$-ary tree of random variables first

and then prune the tree using the selection mechanism of the cSMC algorithm with fixed path $x\in\mathsf{S}$. As a result, each $Z_{t}^{j}$ in the cSMC is associated with some $\hat{Z}_{t}^{i}$. The construction above permits the bound

where $U$ corresponds to the sum of potentials associated with those $Z_{t}^{j}$ whose ancestral lineage does not contain the value $1$. It therefore follows that

because $u\mapsto u/(g+u)$ is increasing. Now, $V$ is a finite non-negative random variable independent of $x$. We may define

which satisfies $\lim_{g\to\infty}\phi(g)=0$ by the monotone convergence theorem.

For the second inequality, we can show similarly that for $G_{t}(x_{t})\geq\underline{G}$

To establish that $P_{N}$ cannot be even geometrically ergodic whenever $\pi_{t}$-${\rm ess}\sup_{x_{t}}G_{t}(x_{t})=\infty$ for some $t\in[T]$ in many settings, we use Proposition 21. This allows for the developments of Proposition 22 and Lemma 24, leading to the desired result under assumptions satisfied in many applications; see Remark 25.

Suppose $P$ is an ergodic Markov kernel on a state space \big{(}\mathsf{X},\mathcal{B}(\mathsf{X})\big{)} with invariant distribution $\pi$. Suppose that for any $\epsilon,\delta>0$ there exists a set $A\in\mathcal{B}(\mathsf{X})$ such that $\pi(A)\in(0,\delta)$ and $\inf_{x\in A}P(x,A)\geq 1-\epsilon$. Then $P$ is not geometrically ergodic.

The result follows directly by following the proof of [26, Theorem 3.1], or by a conductance argument [16, Theorem 1].∎

Then $P_{N}$ cannot be geometrically ergodic.

Because of Proposition 21 it suffices to establish that

because $A_{t}^{1}=1$ by construction. We emphasize that $A_{t}^{i}$ are independent of $x_{t+1:T}$. Now (10) follows directly from (9) because for $i\in\{2,\ldots,N\}$,

For any $\epsilon,\delta>0$ there exists $A_{\epsilon,\delta}$ such that $\pi(A_{\epsilon,\delta})>0$ and for $x\in A_{\epsilon,\delta}$

Because of exchangeability, for any $x$ and $2\leq k\leq N$,

Denote B=\big{\{}\frac{\sum_{k=2}^{N}G_{t}(Z_{t}^{k})}{G_{t}(x_{t})}\geq(N-1)\epsilon\big{\}}, then for $x\in A_{\epsilon,\delta}$ also

We may bound for any $x\in A_{\epsilon,\delta}$,

Letting $\epsilon,\delta\to 0$ completes the proof.∎

Assume that there exists $t\in[T]$ such that $\pi_{t}$-${\rm ess}\sup_{x_{t}}G_{t}(x_{t})=\infty$, and if $t\geq 2$, suppose also that for any $A\in\mathcal{B}(\mathsf{Z}^{1:t-1})$ and $B\in\mathcal{B}(\mathsf{Z})$,

Then, the assumption of Lemma 23 and consequently (9) holds for $t$.

An immediate implication of Propositions 22 and 13 and Lemma 24 is that if $\pi$ is equivalent to a Lebesgue or counting measure on $\mathsf{X}$ then $P_{N}$ is geometrically ergodic for any $N\geq 2$ if and only if $\pi_{t}$-${\rm ess}\sup_{x_{t}}G_{t}(x_{t})<\infty$ for all $t\in[T]$. This covers many applications in statistics, where often the potentials $G_{t}$ are strictly positive and for any $x_{t}\in\mathsf{Z}$, the Markov kernel $M_{t}(x_{t},\cdot)$ is equivalent to a Lebesgue or counting measure on $\mathsf{Z}$.

Proposition 22 does not characterize all situations in which $P_{N}$ fails to be geometrically ergodic. Indeed, in the following example (9) does not hold, and $P_{N}$ still fails to be geometrically ergodic.

Our findings above suggest that the essential boundedness of the potentials could in fact be a necessary condition for geometric ergodicity. We have considered also various other examples, and it seems that in any specific scenario it is easy to identify “sticky” sets and conclude by Lemma 21. However, we have yet to identify such sets in general, and so have resorted to stating the following.

The particle Gibbs sampler

In numerous situations of practical interest one is interested in sampling from a probability distribution \pi\bigl{(}{\rm d}\theta\times{\rm d}x\bigr{)} defined on some measurable space \bigl{(}\Theta\times\mathsf{X},\mathcal{B}(\Theta)\times\mathsf{\mathcal{B}(X)}\bigr{)} for which direct sampling is difficult, but sampling from the associated conditional probability distributions $\pi_{\theta}({\rm d}x)$ and $\pi_{x}({\rm d}\theta)$ for any $(\theta,x)\in\Theta\times\mathsf{X}$ turns out to be easier. In fact when sampling exactly from these conditionals is possible one can define the two stage Gibbs sampler which alternately samples from these conditional distributions. More precisely, let us define, for any $(\theta,x)\in\Theta\times\mathsf{X}$ and $S\in\mathcal{B}(\Theta)\times\mathsf{\mathcal{B}(X)}$,

This can be interpreted as a Markov transition probability, and is precisely the Markov kernel underpinning the standard two stage Gibbs sampler. The corresponding Markov chain $\{(\theta_{i},X_{i}),i\geq 0\}$ on $\Theta\times\mathsf{X}$ leaves $\pi$ invariant and is ergodic under fairly general and natural conditions. In fact it can be shown that $\{X_{i},i\geq 0\}$ and $\{\theta_{i},i\geq 0\}$ are themselves Markov chains leaving the marginals \pi\bigl{(}{\rm d}x\bigr{)} and $\pi({\rm d}\theta)$ invariant respectively. For reasons which will appear clearer below, we define for any (x_{0},S)\in\mathsf{X}\times\mathcal{B}\bigl{(}\mathsf{X}\bigr{)} the Markov transition probability \Gamma_{x}\bigl{(}x_{0},S\bigr{)}:=\Gamma\bigl{(}x_{0},\Theta\times S\bigr{)} corresponding to the Markov chain $\{X_{i},i\geq 0\}$ (we point out that the index $x$ in this notation is a name, not a variable). In some situations, however, while sampling from the conditional distribution \pi_{x}\bigl{(}{\rm d}\theta\bigr{)} may be routine, sampling from $\pi_{\theta}({\rm d}x)$ may be difficult and this step is instead replaced by a Markov transition probability $\Pi_{\theta}(x,{\rm d}y)$ leaving $\pi_{\theta}({\rm d}x)$ invariant for any $\theta\in\Theta$. The resulting algorithm, whose transition kernel $\Phi$ is given below, is often referred to as “Metropolis-within-Gibbs” in the common situation where $\Pi_{\theta}$ is a Metropolis–Hastings transition kernel—we will however use this name in order to refer to the general scenario. In the particular situation where $\Pi_{\theta}$ is a cSMC transition kernel the resulting algorithm is known as the particle Gibbs (PGibbs) sampler . We note that in the general scenario, for any (\theta_{0},x,S)\in\Theta\times\mathsf{X}\times\big{(}\mathcal{B}(\Theta)\times\mathsf{\mathcal{B}(X)}\big{)}

Similarly to above one can show that $\{X_{i},i\geq 1\}$ defines a Markov chain, with transition kernel, for (x_{0},S)\in\mathsf{X}\times\mathcal{B}\bigl{(}\mathsf{X}\bigr{)}, $\Phi_{x}(x_{0},S):=\Phi(x_{0},\Theta\times S)$ which is \pi\bigl{(}{\rm d}x\bigr{)}-reversible, and positive as soon as $\Pi_{\theta}$ defines a positive operator for any $\theta\in\Theta$. Indeed since for any f,g\in L^{2}\bigl{(}\mathsf{X},\pi\bigr{)},

Let $\pi$ be a probability distribution defined on \bigl{(}\Theta\times\mathsf{X},\mathcal{B}(\Theta)\times\mathcal{B}(\mathsf{X)}\bigr{)} and let $\left\{\Pi_{\theta},\theta\in\Theta\right\}$ be a family of Markov transition probabilities $\left\{\Pi_{\theta},\theta\in\Theta\right\}$ such that for any $\theta\in\Theta$ the Markov kernel $\Pi_{\theta}$ is reversible with respect to $\pi_{\theta}$, and let $\Gamma$ and $\Phi$ be as in (11) and (12). Define

Then, for any $f\in L^{2}(\mathsf{X},\pi)$ we have the following inequalities,

where the latter inequality holds for $\varrho>0.$

there exist $\epsilon>0$ such that for all $\theta\in\Theta$ and all (x,B)\in\mathsf{X}\times\mathcal{B}\big{(}\mathsf{X}\big{)}, the minorisation inequality \Pi_{\theta}\big{(}x,B\big{)}\geq\epsilon\pi_{\theta}\big{(}B\big{)} holds, then for any $f\in L^{2}(\mathsf{X},\pi)$

for all $\theta\in\Theta$, $\Pi_{\theta}$ is a positive operator then for any $f\in L^{2}(\mathsf{X},\pi)$

We prove the first point. Without loss of generality we consider any f\in L_{0}^{2}\bigl{(}\mathsf{X},\pi\bigr{)} and notice that

Now using that $\mathcal{E}_{\Gamma_{x}}(f)=\frac{1}{2}\int_{\Theta\times\mathsf{X}^{2}}\pi({\rm d}x)\pi_{x}({\rm d}\theta)\pi_{\theta}({\rm d}y)\left[f(x)-f(y)\right]^{2}=\int_{\Theta}\pi({\rm d}\theta){\rm var}_{\pi_{\theta}}\left(f\right)$ and letting $\bar{f}_{\theta}:=f-\pi_{\theta}(f)$ for any $\theta\in\Theta$, we obtain

where we have used that for any g\in L_{0}^{2}\bigl{(}\mathsf{X},\pi\bigr{)}, \mathcal{E}_{\Pi_{\theta}}\bigl{(}g\bigr{)}\leq 2{\rm var}_{\pi_{\theta}}\left(g\right) and that the set A:=\bigl{\{}\theta\in\Theta:{\rm var}_{\pi_{\theta}}(\bar{f}_{\theta})=\infty\bigr{\}} satisfies \pi\bigl{(}A\times\mathsf{X}\bigr{)}=0. The latter result follows from ${\rm var}_{\pi}(f)<\infty$ and the variance decomposition identity: $\|f\|_{\pi}^{2}=\|f-\bar{f}_{\theta}\|_{\pi}^{2}+\|\bar{f}_{\theta}\|_{\pi}^{2}$. We deduce (a) from the last inequality. Points (b) and (c) then follow from [2, Lemma 34].

We next turn into (d). As above, we find that

it may be easier in practice to use the lower bound \underline{\varrho}:=\inf_{\theta\in\Theta}{\rm Gap}\bigl{(}\Pi_{\theta}\bigr{)}\leq\varrho which leads to {\rm Gap}\bigl{(}\Phi_{x}\bigr{)}\geq\underline{\varrho}\times{\rm Gap}\left(\Gamma_{x}\right) and {\rm var}\bigl{(}f,\Phi_{x}\bigr{)}\leq(\underline{\varrho}^{-1}-1){\rm var}_{\pi}(f)+{\rm\underline{\varrho}^{-1}}{\rm var}\left(f,\Gamma_{x}\right) when $\underline{\varrho}>0$,

one could suggest iterating $\Pi_{\theta}$ sufficiently many times, say $k_{\theta}$ times, in order to ensure that $\Pi_{\theta}^{k_{\theta}}$ satisfies the uniform in $\theta$ properties of the type suggested above. This would require however a computable quantitative bound on the spectral gap of $\Pi_{\theta}$ ,

the lower bound in (c) is motivated by the fact that $\{\Pi_{\theta},\theta\in\Theta\}$ may be a family with non-positive elements, which may introduce negative correlations. On the contrary in the situation where $\{\Pi_{\theta},\theta\in\Theta\}$ is a collection of positive operators (e.g. cSMC kernels) then (b) implies that $\Phi_{x}$ is geometrically ergodic as soon as $\Gamma_{x}$ is geometrically ergodic and $\varrho>0$ (and of course $\Gamma_{x}$ is always positive) and (d)(ii) that $\Phi_{x}$ is always inferior to $\Gamma_{x}$ in terms of asymptotic variance. In the context of the PGibbs sampler the latter result parallels what is known for pseudo-marginal algorithms ,

we note that from [25, Theorem 1; Proposition 1] $\Phi$ is geometrically ergodic as soon as $\Phi_{x}$ is geometrically ergodic.

Now we show how these results can be transferred to the $\{\theta_{i}\}$ chain.

Let the notation be as in Theorem 27. Then,

for any f\in L^{2}\bigl{(}\Theta,\pi\bigr{)}, letting for any $x\in\mathsf{X}$ \bar{f}(x):=\pi_{x}\bigl{(}f\bigr{)}\in L^{2}\bigl{(}\mathsf{X},\pi\bigr{)}, we have for any $k\geq 1$

if $\varrho>0$ defined in 13, then for f\in L^{2}\bigl{(}\Theta,\pi\bigr{)}

if for all $\theta\in\Theta$, $\Pi_{\theta}$ is a positive operator, then for f\in L^{2}\bigl{(}\Theta,\pi\bigr{)} ${\rm var}(f,\Phi)\geq{\rm var}(f,\Gamma)$.

We remark that without loss of generality we can let $f\in L_{0}^{2}(\Theta,\pi)$ throughout. First note that for f\in L_{0}^{2}\bigl{(}\Theta,\pi\bigr{)} and any \bigl{(}\theta,x_{0}\bigr{)}\in\Theta\times\mathsf{X}

and for g\in L^{2}\bigl{(}\mathsf{X},\pi\bigr{)} and any $p\geq 1$, $\Phi^{p}(\theta,x_{0};g)=\Phi_{x}^{p}(x_{0};g)$. The first result is straightforward upon remarking that for $k\geq 1$

For the second and third point, using the remarks above, for f\in L_{0}^{2}\bigl{(}\Theta,\pi\bigr{)} and $k\geq 1$

Now $\|f\|_{\pi}^{2}=\left\langle f-\bar{f}+\bar{f},f-\bar{f}+\bar{f}\right\rangle_{\pi}=\|\bar{f}\|_{\pi}^{2}+\|f-\bar{f}\|_{\pi}^{2}$, which is the variance decomposition identity and by noting that \pi\bigl{(}\bar{f}\bigr{)}=0 lets us deduce that f\in L_{0}^{2}\bigl{(}\Theta,\pi\bigr{)} implies that \bar{f}\in L_{0}^{2}\bigl{(}\mathsf{X},\pi\bigr{)}. Now,

We conclude by noting that for f\in L^{2}\bigl{(}\mathsf{X},\pi\bigr{)} then {\rm var}_{\pi}\bigl{(}f\bigr{)}=\|f-\pi(f)\|_{\pi}^{2} and {\rm var}_{\pi}\bigl{(}\bar{f}\bigr{)}=\|\bar{f}-\pi(f)\|_{\pi}^{2}=\|\overline{f-\pi(f)}\|_{\pi}^{2}. We will also use the equality above for $\Gamma$ and $\Gamma_{x}$, since again the latter corresponds to a particular instance of the above. We can now use the bound from Theorem 27, which leads, for f\in L_{0}^{2}\bigl{(}\Theta,\pi\bigr{)}, to

We conclude as above. The final statement follows from {\rm var}\bigl{(}\bar{f},\Phi_{x}\bigr{)}\geq{\rm var}\bigl{(}\bar{f},\Gamma_{x}\bigr{)} (see Theorem 27) and the equality established above for $\Phi$ and $\Phi_{x}$ and $\Gamma$ and $\Gamma_{x}$.∎

Consider the PGibbs sampler with $N\geq 2$ particles with kernel $\Phi_{N}$ defined as in (12) such that for any $\theta\in\Theta$, $\Pi_{\theta}=P_{\theta,N}$ is the i-cSMC kernel as defined in Section 3 for the families $\{M_{\theta,t}\}$and $\{G_{\theta,t}\}$ of kernels and potentials on \mathsf{Z}\times\mathcal{B}\bigl{(}\mathsf{Z}\bigr{)} and $\mathsf{Z}$ respectively. For any $\theta\in\Theta$ we let $\gamma_{\theta,T}$ be the corresponding normalizing constant as defined below (3). Then, the results of Theorems 27 and 29 hold as follows:

then $\varrho\geq\epsilon_{N}$ as defined in Corollary 14,

or we have the uniform mixing condition, for some $0\leq\alpha<\infty$,

then $\varrho\geq\epsilon_{N}$ as defined in Corollary 16.

In particular, in both cases $\varrho$ convergences to one as $N\rightarrow\infty$, implying that the spectral gaps and the asymptotic variances associated with the PGibbs sampler converge to those of the related Gibbs sampler.

It is worth noting that terms related to $\gamma_{\theta,T}$ appear in all these bounds. So, for example in the first part it is not sufficient that our potentials $\{G_{\theta,t}\}$ are essentially bounded, but it is sufficient if, for all $t\in[T]$, $\pi_{t}-{\rm ess}\sup_{\theta,x_{t}}G_{\theta,t}(x_{t})/\eta_{\theta,t}(G_{t})$ is bounded.

Discussion

The developments above go some way in characterizing the behaviour of i-cSMC and associated PGibbs Markov chains, and raise a number of possible future directions for research. We have already embarked upon investigating some potentially practical uses of the minorization conditions and spectral properties for these chains. Of particular interest in practice is how to choose $N$ in the i-cSMC algorithm so as to balance the trade off between mixing properties of $P_{N}$ and the total number of iterations that can be performed with limited computational resources. Remark 17, for example, can be used to find approximately good values of $N$ in this spirit, but can only serve as a heuristic. In particular, while Proposition 8 may provide a fairly accurate bound in the large $N$ regime, it is unclear how much is lost in applying Jensen’s inequality, and consequently how accurate estimates such as those in Remark 17 can be. It is possible that results such as those in may provide a way to exploit additional structure often found in statistical applications.

The results for the i-cSMC and PGibbs Markov chains developed here can be compared and contrasted with similar results for the Particle Independent Metropolis–Hastings (PIMH) and PMMH Markov chains . We summarize here the detailed comparison provided in [2, Appendix F]. Like i-cSMC, PIMH is an exact approximation of an independent sampler but PMMH is an exact approximation of an idealized Metropolis–Hastings kernel, rather than a Gibbs sampler. Just as i-cSMC can be viewed as a constituent element of PGibbs, PIMH can be viewed as playing the same role within PMMH. Central to the analysis of PIMH is the essential supremum of the normalizing constant estimate \hat{\gamma}_{T}^{N}\bigl{(}Z_{1:T}\bigr{)} introduced in Section 4 with respect to the law of a standard SMC algorithm and indeed the PIMH Markov chain is (uniformly) geometrically ergodic if and only if this supremum is finite as a consequence of the characterisation of independent Metropolis–Hastings chains in . However, it can also be seen that the rate of convergence of PIMH will typically not improve as $N$ increases, in contrast with the convergence for the i-cSMC (see Propositions 13 and 15).

For PMMH, show that if the essential supremum of the relative normalizing constant estimate \hat{\gamma}_{\theta,T}^{N}\bigl{(}Z_{1:T}\bigr{)}/\gamma_{\theta,T} is moreover bounded essentially uniformly in $\theta$ then the existence of a spectral gap of the idealized Metropolis–Hastings Markov kernel it approximates is inherited by PMMH. However, the rate of convergence of the PMMH Markov chain when this occurs does not improve in general as $N$ increases, in contrast to our results for PGibbs Markov chains. In this context, weak convergence in $N$ of the asymptotic variance of estimates of $\pi(f)$ to the corresponding asymptotic variance of the Metropolis–Hasting kernel is nevertheless provided by [5, Proposition 19] for all $f\in L^{2}(\Theta,\pi)$ but this can be contrasted with quantitative bounds obtained in Theorem 29.

The one step uniform minorization condition in Corollary 9, where the minorization measure is the invariant distribution of the Markov chain, suggests that it may be possible to apply coupling from the past techniques (see, e.g., ) in order to produce samples from exactly this distribution. It is, however, not clear how to implement such an algorithm in general, although provides a perfect simulation algorithm motivated by Theorem 1. Finally, our analysis has focused mainly on the case where the essential boundedness condition holds. However, a refined analysis may permit characterization of the i-cSMC and hence the PGibbs Markov chains even in the absence of this condition, with parallels to .

CA’s research was supported by EPSRC EP/K009575/1 Bayesian Inference for Big Data with Stochastic Gradient Markov Chain Monte Carlo and EP/K014463/1 Intractable Likelihood: New Challenges from Modern Applications (ILike). MV was supported by Academy of Finland grant 250575.

Supplementary material

The proof of Lemma 7 is a simple consequence of Lemma 32 (b). We introduce the set of indices $\mathcal{J}_{T}:=\bigcup_{m=0}^{T}\{1\}^{m}\times\{2,\ldots,N\}^{T-m}$, which will allow us to define the lineages coalescing with $\mathbf{1}\in\{1\}^{T}$ at some point in the past, and $m_{\mathbf{i}}:=\max\{k:i_{k}=1\}$ (with the convention that $\max\emptyset=0$) the time at which coalescence occurs.

For any $x\in\mathsf{X}$, $z_{1:T}\in\mathsf{X}^{T}$ and $a_{1:T}\in[N]^{N(T-1)}\times[N]$,

for any $y_{2:T}\in\mathsf{Z}^{T-1}$ and $\mathbf{k}=k_{1:T}\in[N]^{T}$ such that $k^{1}\neq 1$

and for $t\in\{2,\ldots,T\}$, any $(y_{1},\ldots,y_{t-1},y_{t+1},\ldots,y_{T})\in\mathsf{Z}^{T-1}$, $\mathbf{k}\in[N]^{T}$ such that $k_{t}\neq 1$ and $a_{t-1}^{k_{t}}=k_{t-1}$

for $\mathbf{i}\in\mathcal{J}_{T}$ and $y_{1:m_{\mathbf{i}}}=x_{1:m_{\mathbf{i}}}$ we have

We note that the above is well defined for $m_{\mathbf{i}}=0$ from the definition of $M_{p,l}$ in Section 3 and associated remark, and the convention that $x_{1:0}=y_{1:0}$ should be ignored in this case.

In order to alleviate notation we omit $Z_{t}\in\cdot,Z_{t-1}=\cdot$ and $A_{t-1}=\cdot$ and set G_{t}^{k}:=G_{t}\bigl{(}z_{t}^{k}\bigr{)}. For the first point we note the independence on $(y_{2},\ldots,y_{T})\in\mathsf{Z}^{T-1}$ of

and conclude from (4). Similarly we note the independence on $(y_{1},\ldots,y_{t-1},y_{t+1},\ldots,y_{T})\in\mathsf{Z}^{T-1}$ of

and we conclude with (5). For the second point, let $\mathbf{i}\in\mathcal{J}_{T}$, a_{1:T}\in\bigl{(}[N]^{N}\bigr{)}^{T-1}\times[N] such that $a_{t-1}^{i_{t}}=i_{t-1}$ for $t=m_{\mathbf{i}}+1,\ldots,T$, $a_{T}=i_{T}$ and $y_{1:m_{\mathbf{i}}}=x_{1:m_{\mathbf{i}}}$ then, with an obvious convention when $m_{\mathbf{i}}=1$ (i.e. $a_{0}$ does not exist and should be ignored), we have

Appendix B Proof of Lemma 10

We can define the artificial joint distribution

It is straightforward to check that the conditional distribution of ${\bf K}$ given $(z_{1:T},a_{1:T-1})$ can be written

Appendix C Supplementary material for Section 4

In the next proposition we gather general properties for generic reversible Markov chains satisfying a uniform minorization condition for which the minorization probability is precisely the invariant distribution of the Markov chain. We suspect these results to be widely known, but could not find a relevant reference. Let $L^{2}(\mathsf{E},\mu)$ and L_{0}^{2}\bigl{(}\mathsf{E},\mu\bigr{)}:=\bigl{\{}f\in L^{2}(\mathsf{E},\mu):\mu\bigl{(}f\bigr{)}=0\bigr{\}} both endowed with the inner product defined for any $f,g\in L^{2}(\mathsf{E},\mu)$ as $\left\langle f,g\right\rangle_{\mu}:=\int_{\mathsf{E}}f(x)g(x)\mu({\rm d}x)$, which yields the associated norm $\|f\|_{\mu}:=\sqrt{\left\langle f,f\right\rangle_{\mu}}$. For any $f\in L^{2}(\mathsf{E},\mu)$ we define the Dirichlet forms

where $I$ is the identity operator. The right and left spectral gaps of a generic reversible Markov transition kernel have the following variational representation

The condition ${\rm Gap}\left(\Pi\right)>0$ and ${\rm Gap}_{L}\left(\Pi\right)>0$ implies geometric ergodicity of the Markov chain. It turns out that convergence is in fact uniformly geometric in the following scenario.

Let $\mu$ be a probability distribution on some measurable space \bigl{(}\mathsf{E},\mathcal{B}\bigl{(}\mathsf{E}\bigr{)}\bigr{)} and let \Pi:\mathsf{E}\times\mathcal{B}\bigl{(}\mathsf{E}\bigr{)}\rightarrow be a Markov transition kernel reversible with respect to $\mu$. Assume that there exists $\varepsilon>0$ such that for any (x,A)\in\mathsf{E}\times\mathcal{B}\bigl{(}\mathsf{E}\bigr{)},

the Dirichlet forms satisfy for any f\in L^{2}\bigl{(}\mathsf{E},\mu\bigr{)}

for any probability distribution $\nu\ll\mu$ we have

and for any f\in L^{2}\bigl{(}\mathsf{E},\mu\bigr{)}

and if $\Pi$ is a positive operator then naturally {\rm var}\bigl{(}f,\Pi\bigr{)}\geq{\rm var}_{\mu}\bigl{(}f\bigr{)}.

First, from the minorization condition one can write $\Pi(x,{\rm d}y)=\varepsilon\mu({\rm d}y)+(1-\varepsilon)R_{\Pi,\varepsilon}(x,{\rm d}y)$, where $R_{\Pi,\varepsilon}(x,A):=\frac{\Pi(x,A)-\varepsilon\mu(A)}{1-\varepsilon}$ is $\mu-$invariant. Now for f\in L_{0}^{2}\big{(}\mathsf{E},\mu\big{)}

and therefore with $\mathcal{E}_{\mu}(f)=\left\langle f,f\right\rangle_{\mu}$ the Dirichlet form of the (reversible) “independent samples” Markov chain we deduce

which implies (a). The bounds on the spectral gaps (b) follow immediately and the results in points (c) and (d) are now a consequence of the resulting property of the spectrum and e.g. [27, Proposition 3.12, p. 44] and [14, Proposition 1.5]. Result (e) is due to Doeblin , while the two bounds on the asymptotic variance are direct consequences of Lemma 34 and coincide in this case with the “Kipnis–Varadhan” upper bound . ∎

Let $\Pi_{1},\Pi_{2}$ be reversible with respect to $\mu$ and assume that there exists $\varrho\geq 0$ such that for any f\in L_{0}^{2}\bigl{(}\mathsf{E},\mu\bigr{)}

The first result is straightforward. For the second result, first notice that

and since {\rm var}\bigl{(}f,\Pi\bigr{)}=2\bigl{[}\sup_{g\in L_{0}^{2}\bigl{(}\mathsf{E},\mu\bigr{)}}2\bigl{\langle}f,g\bigr{\rangle}_{\mu}-\mathcal{E}_{\Pi}(g)\bigr{]}-\|f\|_{\mu}^{2} we conclude that

Appendix D Supplementary material for Section 5

The proof of Proposition 11 relies on the following technical lemma, and is given after this intermediate result.

for any $k=1,\dots,T-1$, any $z_{1:T-k}\in\mathsf{Z}^{N(T-k)}$ such that $(z_{1:T-k}^{1},z_{1:T-\text{k}}^{2})=(x_{1:T-k},y_{1:T-k})$

and $\mathcal{I}_{k,s}$ and $C_{k,s}$ are as in Proposition 11.

The property in (a) is immediate from the linearity of the expectation and the definition of the process. We now prove property (b) by induction on $k=1,\ldots,T-1$. In order to alleviate notation we let G_{p,q}^{i}:=G_{p,q}\bigl{(}Z_{p}^{i}\bigr{)} when found inside an expectation and G_{p,q}^{i}:=G_{p,q}\bigl{(}z_{p}^{i}\bigr{)} otherwise, G_{p,q}^{1+2}:=G_{p,q}\bigl{(}x_{p}\bigr{)}+G_{p,q}\bigl{(}y_{p}\bigr{)} and $C_{k,s}(\mathbf{i}):=C_{k,s}(\mathbf{i},x,y)$. The case $k=1$ follows from (a) with $t=T$ by observing that \mathcal{I}_{1,1}=\bigl{\{}T+1\bigr{\}}, C_{1,1}\bigl{(}\mathbf{i},x,y\bigr{)}=1 and that G_{T-1,T+1}^{r}=Q_{T-1,T}\bigl{(}G_{T}\bigr{)}\bigl{(}z_{T-1}^{r}\bigr{)} :

Now we assume the property true for some $k\in\{1,\ldots,T-2\}$ and establish it for $k+1$. We have

and we deal with the two terms separately. Observe that $A_{T-k}$ only depends on $x_{T-k+1:T}$ and $y_{T-k+1:T}$, then by application of the first result of the lemma we obtain

and, noting that $C_{k,s}(\mathbf{i})$ depends on $x_{T-k+2:T},y_{T-k+2:T}$ only

where we have again applied the first result of the lemma. Consequently we can group the terms as follows

Now we first focus on the first term on on the RHS on the first line (with the sum now written in extension in order to help and we note that we do not use the double indexing $i_{j}^{s}$ in order to keep notation simple),

where we have used the following changes of variables: $j_{m}=i_{m-1}$ for $m=1,\ldots,s+1$ followed by $s=s^{\prime}-1$. Note that we can extend the sum in order to include the term $s^{\prime}=1$, since we cannot have $j_{1}=T+1\neq T-k+1=j_{1}$. We examine the second term on the RHS of the first line of (15)

and we notice that we can extend the sum in order to include the term $s=k+1$ because $\sharp\left\{T-k+2,\ldots,T+1\right\}=k$, which implies that $\mathcal{I}_{k,k+1}=\emptyset$. Consequently we deduce that

We now turn to the second line of (15) and examine the two terms within the brackets and use similar ideas. First we have

We start with the second result of Lemma 35 for $k=T-1$ and we proceed as in the beginning of the proof of that lemma, using similar notation and arguments. Here we have however

and using arguments similar to those of the proof of Lemma 35,

which can be rewritten as (again we use $s^{\prime}-1=s$ and the fact that $G_{0,1}=G_{0}=1$ by convention)

Consider first the case where $k\geq m$, for $z_{p},z^{\prime}_{p}\in\mathsf{Z}^{2}$,

by using (A4)(b) and a straightforward induction. Now we can conclude by using (A4)(a). When $k<m$ we simply note that, proceeding as above, for any $z_{p},z^{\prime}_{p}\in\mathsf{X}^{2}$,

Appendix E Supplementary material for Section 6

Denote by $B_{r}$ the closed ball of radius $r$ centred at the origin and define the sets

We may define $A:=\cap_{k\geq 1}A_{k,r_{\epsilon,k}}$ which satisfies $\xi(A^{\complement})\leq\sum_{k=1}^{\infty}\epsilon 2^{-k}=\epsilon$. ∎

Appendix F Detailed comparisons with the PIMH and PMMH

In this section we contrast the performance properties of the i-cSMC (resp. PGibbs sampler), as established in Section 5 (resp. Section 7), with those of the Particle Independent Metropolis–Hastings kernel (PIMH) (resp. particle Marginal Metropolis–Hastings (PMMH)) also proposed in , which also aims to (indirectly) sample from $\pi$ as defined in Section 3 (resp. Section 7). We use notation similar to that used in Section 3 for the i-cSMC algorithm. The Markov kernel of the PIMH can be defined for $(a,z)\in\mathsf{W}$ (with an obvious abuse of notation in order to alleviate notation), W\in\mathcal{B}\bigr{(}\mathsf{W}\bigr{)} and $N\geq 1$ as

$A=A_{1:T}$ and $Z:=Z_{1:T}$. We note that $\hat{\gamma}_{T}^{N}(z)$ is not a random quantity in the definition of $\check{P}_{N}(a,z;W)$. The invariant distribution of the Markov chain, which evolves on $\mathsf{W}$, is given for any W\in\mathcal{B}\bigr{(}\mathsf{W}\bigr{)} by

Clearly $\check{\epsilon}_{N}>0$ whenever $\check{\pi}_{N}-{\rm ess}\sup_{z}\hat{\gamma}_{T}^{N}(z)<\infty$, which is similar to what we have obtained in Propositions 13 and 15 for the i-cSMC. An important difference, which may explain the widely perceived superiority of the i-cSMC, is that the rate of convergence of PIMH will typically not improve (and in particular converge to $1$) as $N$ increases, even for bounded potentials, which is in contrast with the corresponding convergence rate of the i-cSMC (see Propositions 13 and 15).

We can also compare the results of Section 7 for the PGibbs sampler with the corresponding results for the PMMH algorithm . This latter algorithm evolves on $\Theta\times\mathsf{W}$ with transition probability

${\rm Gap}(\check{\Phi}_{N})>0$ whenever ${\rm Gap}(\Phi^{*})>0$, i.e. the existence of a spectral gap of $\Phi^{*}$ is “inherited” by $\check{\Phi}_{N}$. This coincides in many cases with inheritance of geometric ergodicity, for example when $\check{\Phi}_{N}$ is positive.

The rate of convergence of a geometrically ergodic PMMH Markov chain does not improve in general as $N$ increases, in contrast to our results for PGibbs Markov chains. In this context, weak convergence in $N$ of the asymptotic variance of estimates of $\pi(f)$ using $\check{\Phi}_{N}$ to that of $\Phi^{*}$ is nevertheless provided by [5, Proposition 19] for all $f\in L^{2}(\Theta,\pi)$. This can be contrasted with quantitative bounds obtained in Theorem 29.

References