Limit theorems for cloning algorithms

TL;DR

This paper provides a rigorous theoretical framework for analyzing the convergence and errors of cloning algorithms used in estimating large deviations in continuous-time stochastic processes, applicable to jump processes on compact spaces.

Contribution

It establishes the first comprehensive, fully rigorous bounds on systematic and random errors of cloning algorithms in continuous time, avoiding time discretization.

Findings

Bounds systematic and random errors of cloning algorithms

Applicable to a large class of jump processes on compact spaces

Provides a framework for evaluating and improving algorithm efficiency

Abstract

Large deviations for additive path functionals of stochastic processes have attracted significant research interest, in particular in the context of stochastic particle systems and statistical physics. Efficient numerical `cloning' algorithms have been developed to estimate the scaled cumulant generating function, based on importance sampling via cloning of rare event trajectories. So far, attempts to study the convergence properties of these algorithms in continuous time have only led to partial results for particular cases. Adapting previous results from the literature of particle filters and sequential Monte Carlo methods, we establish a first comprehensive and fully rigorous approach to bound systematic and random errors of cloning algorithms in continuous time. To this end we develop a method to compare different algorithms for particular classes of observables, based on the martingale characterization of stochastic processes. Our results apply to a large class of jump processes on compact state space, and do not involve any time discretization in contrast to previous approaches. This provides a robust and rigorous framework that can also be used to evaluate and improve the efficiency of algorithms.

Equations

$$\Vert f\Vert =\underset{x\in E}{\sup }|f(x)|.$$

$$\Omega:=\big{\{}\omega:[0,\infty)\to E\mathrm{\,right\,continuous\,with\,left\,limits}\big{\}}.$$

$$P(t)\,f(x)=\mathbb{E}_{x}\big{[}f(X_{t})\big{]}\ ,\quad\mu P(t)\,(f):=\int_{E}P(t)f(x)\mu(dx)=\mathbb{E}_{\mu}\big{[}f(X_{t})\big{]}\ ,$$

$$\frac{d}{dt}P(t)\,f=\,\mathcal{L}\big{(}P(t)\,f\big{)}\,=\,P(t)\,\mathcal{L}(f),$$

$$W(x,dy):=\lambda (x)\cdot p(x,dy)$$

$$\mathcal{L}(f)(x)={\int }_{E}W(x,dy)[f(y)-f(x)],\forall f\in {\mathcal{C}}_b(E),x\in E.$$

$${\mathcal{L}}^{\mathcal{V}}(f)(x):=\mathcal{L}(f)(x)+\mathcal{V}(x)f(x)\text{mbox}definedforallf\in C_b(E).$$

$$\frac{d}{dt}P^{\mathcal{V}}(t)f=\mathcal{L}^{\mathcal{V}}\big{(}P^{\mathcal{V}}(t)f\big{)}\quad\mbox{with}\quad P^{\mathcal{V}}(0)f=f$$

$$\big{\|}e^{-t\lambda_{0}}P^{\mathcal{V}}(t)f(\cdot)-\mu_{\infty}(f)\big{\|}\,\leq\,\|f\|\cdot\alpha\rho^{t}\ ,$$

$$\nu_{t,\mu_{0}}(f):=\mu_{0}\big{(}P^{\mathcal{V}}(t)f\big{)},$$

$$\frac{d}{dt}\nu_{t,\mu_{0}}(f)=\nu_{t,\mu_{0}}\big{(}\mathcal{L}^{\mathcal{V}}(f)\big{)}=\nu_{t,\mu_{0}}\big{(}\mathcal{L}(f)+\mathcal{V}\cdot f\big{)},$$

$${\mu }_{t,{\mu }_0}(f):=\frac{{\nu }_{t,{\mu }_0}(f)}{{\nu }_{t,{\mu }_0}(1)},$$

$$\big{|}\mu_{t,\mu_{0}}(f)-\mu_{\infty}(f)\big{|}\leq\|f\|\cdot{\tilde{\alpha}\rho^{t}}\ ,$$

$$\frac{\mu_{\infty}(f)\,-\,\|f\|\,\alpha\cdot\rho^{t}}{1+\alpha\cdot\rho^{t}}\leq\mu_{t,\mu_{0}}(f)\,=\,\frac{\mu_{0}\big{(}e^{-t\lambda_{0}}P^{\mathcal{V}}(t)f\big{)}}{\mu_{0}\big{(}e^{-t\lambda_{0}}P^{\mathcal{V}}(t)1\big{)}}\,\leq\,\frac{\mu_{\infty}(f)\,+\,\|f\|\,\alpha\cdot\rho^{t}}{1-\alpha\cdot\rho^{t}}\ ,$$

$$\frac{d}{dt}\mu_{t}(f)=\mu_{t}\big{(}\mathcal{L}(f)\big{)}+\mu_{t}(\mathcal{V}f)-\mu_{t}(f)\cdot\mu_{t}(\mathcal{V}).$$

$$\frac{d}{dt}{\mu }_t(f)$$

$$\displaystyle=\frac{1}{\nu_{t}(1)}\cdot\nu_{t}\big{(}\mathcal{L}(f)\,+\,\mathcal{V}\cdot f\big{)}-\frac{\nu_{t}(f)}{\nu_{t}(1)^{2}}\,\nu_{t}\big{(}\mathcal{L}(1)\,+\,\mathcal{V}\big{)}$$

$$\displaystyle=\mu_{t}\big{(}\mathcal{L}(f)\big{)}+\mu_{t}(\mathcal{V}f)-\mu_{t}(f)\cdot\mu_{t}(\mathcal{V})\ .$$

$$\mu\big{(}\widetilde{\mathcal{L}}_{\mu}(f)\big{)}\,=\,\mu(\mathcal{V}f)\,-\,\mu(f)\cdot\mu(\mathcal{V})\ ,$$

$$\widetilde{\mathcal{L}}_{\mu}(f)(x)=\int_{E}\widetilde{W}(x,y)\big{(}f(y)-f(x)\big{)}\mu(dy)\ ,$$

$$\mu\big{(}\widetilde{W}(\cdot,x)-\widetilde{W}(x,\cdot)\big{)}\,=\,\mathcal{V}(x)-\mu(\mathcal{V})\ ,$$

$$\tilde{W}(y,x)-\tilde{W}(x,y)=\mathcal{V}(x)-\mathcal{V}(y),$$

$$\displaystyle\mu\big{(}\widetilde{\mathcal{L}}_{\mu}(f)\big{)}$$

$$\displaystyle=\int_{E^{2}}\big{(}\widetilde{W}(y,x)-\widetilde{W}(x,y)\big{)}\,f(x)\mu(dy)\mu(dx)\ .$$

$${\overline{\mathcal{L}}}_{\mu }:=\mathcal{L}+{\tilde{\mathcal{L}}}_{\mu }\text{mbox}forall\mu \in \mathcal{P}(E),$$

$$\frac{d}{dt}\mu_{t}(f)\,=\,\mu_{t}\big{(}\overline{\mathcal{L}}_{\mu_{t}}(f)\big{)}\ ,$$

$$\Theta_{t,T}f(x):=\frac{P^{\mathcal{V}}(T-t)f(x)}{\mu_{t}\big{(}P^{\mathcal{V}}(T-t)1\big{)}}\quad\mbox{such that}\quad\mu_{T}(f)=\mu_{t}(\Theta_{t,T}f)$$

$$\widetilde{W}_{c}(x,y)\,=\,\big{(}\mathcal{V}(x)-c\big{)}^{-}\,+\,\big{(}\mathcal{V}(y)-c\big{)}^{+}\ ,$$

$$\widetilde{W}(x,y)\,=\,\big{(}\mathcal{V}(y)-\mathcal{V}(x)\big{)}^{+}\ .$$

$$m(\underline{x})(dy):=\frac{1}{N}\sum _{i=1}^N{\delta }_{x_i}(dy)\in \mathcal{P}(E).$$

$${\mu }_t^N(\cdot ):=m({\xi }_t)(\cdot ).$$

$${\Gamma }_{{\overline{L}}^N}(\gamma ,\phi ):={\overline{L}}^N(\gamma \cdot \phi )-\gamma \cdot {\overline{L}}^N(\phi )-\phi \cdot {\overline{L}}^N(\gamma ),\gamma ,\phi \in {\mathcal{C}}_b(E^N),$$

Full text

65C35,

60F25,

62L20,

60F10,

60J75,

60K35,

cloning algorithm,

dynamic large deviations,

interacting particle systems,

$L^{p}$ convergence,

Feynman-Kac formulae,

jump processes,

,

and

1 Introduction

Cloning algorithms have been introduced to the theoretical physics literature [1, 2] as numerical methods to study large deviations of particle currents and other dynamic observables in stochastic particle systems. They combine importance sampling with a stochastic selection mechanism which is used to evaluate numerically the scaled cumulant generating function for time-additive path functionals of stochastic processes. Based on classical ideas of evolutionary algorithms [3, 4], a fixed size population of copies of the original system evolves in parallel, subject to cloning or killing in such a way as to favor the realization of atypical trajectories contributing to rare events. Various variants of the approach are now applied on a regular basis to different systems and large deviation phenomena of interest [5, 6, 7], including also current fluctuations of non-equilibrium lattice gas models [8, 6, 9, 10], turbulent flows [11], glassy dynamics [12], heat waves in climate models [13] and pressure of the edge-triangle model [14]. Due to its widespread applications, the mathematical justification and convergence properties of the algorithm have recently become a subject of research interest with only partial progress. Formal approaches so far are based on a branching process interpretation of the algorithm in discrete time [15], with limited and mostly numerical results in continuous time [16, 17, 18, 19].

In this paper, we provide a novel interpretation of cloning algorithms through Feynman-Kac models and their particle approximations (see [20, 21, 22, 23] for comprehensive reviews), which is itself an established approach to understanding sequential Monte Carlo methods and particle filtering. Previous results provide rigorous control on convergence properties and error bounds of particle filters and related algorithms, mostly for models in discrete time, beginning with the chain of research initiated by [24] with a recent survey provided in [22]. Fewer results address continuous-time dynamics, dating back to [25] in the filtering context, with a Feynman-Kac-based treatment provided by [20] and references therein; a survey of the filtering literature is provided by [26, Chapter 9]. In the current context, particularly relevant recent works include [27, 28, 29, 30, 31]. This literature generally considers diffusive dynamics and relies upon approximative time-discretisations of those dynamics. Adapting those results to the context of jump processes on locally compact state spaces, for which exact simulation from the dynamics is possible, we can establish first rigorous convergence results for the cloning algorithm in continuous time including $L^{p}$ bounds on the random error and bounds on the systematic error. These bounds include the explicit dependence on the clone size distribution, which is a key parameter of the cloning algorithm. The setting of finite activity pure jump processes in which cloning algorithms are primarily employed allows these algorithms to avoid time discretisation by simulating exactly from the law of the underlying process and allows the use of different approximating particle systems. Similar methods have been previously employed in the probabilistic rare event analysis literature in both discrete and continuous time, via explicit Feynman-Kac approximations, e.g. [32], and splitting algorithms (see [33] and references therein); however, both the underlying processes and approximations considered are quite different to those for which cloning algorithms are usually employed. Practically, an important contribution of our approach is a systematic method to compare different cloning algorithms and particle approximations for particular classes of observables of interest, based on the martingale characterization of continuous-time stochastic processes.

This framework provides a novel perspective on the underlying structure of cloning algorithms in terms of McKean representations [22, Section 1.2.2], and can be used to systematically explore several degrees of freedom in the design of algorithms that can be used to improve performance, as illustrated in [34] for current large deviations of the inclusion process [10]. Here we focus on presenting full rigorous results obtained by applying this approach to a version of the classical cloning algorithm in continuous time [2]. In contrast to previous work in the context of cloning algorithms [15, 16], our mathematical approach does not require a time discretization and works in the very general setting of a pure jump Markov process on a locally compact state space. This covers in particular any finite-state Markov chain or stochastic particle systems on finite lattices.

The paper is organized as follows. In Section 2 we introduce general Feynman-Kac models associated to pure jump Markov processes and show that they can be interpreted as the law of a non-linear Markov process, known as a McKean interpretation [21]. In Section 3 we introduce particle approximations for Feynman-Kac models, including classical mean-field versions and cloning algorithms. We provide generalized conditions for convergence as our main result (proved in Section 4), and use this to establish rigorous convergence bounds for cloning algorithms. In Section 5 we introduce large deviations and scaled cumulant generating functions (SCGF) of additive observables for pure jump Markov processes and discuss how the results presented in Section 3 can be applied to estimate the SCGF. We conclude with a short discussion in Section 6.

2 Mathematical Setting

2.1 Dynamics and Feynman-Kac models

We consider a continuous-time homogeneous Feller process \big{(}X_{t}:t\geq 0\big{)} taking values on a locally compact Polish state space $(E,\mathcal{B}(E))$, where $\mathcal{B}(E)$ is the Borel field on $E$. We denote by $\mathcal{M}(E)$ and $\mathcal{P}(E)$ the sets of measures and probability measures, respectively, on $(E,\mathcal{B}(E))$. \big{(}P(t):t\geq 0\big{)} describes the semigroup associated with $X_{t}$, which is considered as acting on the Banach space $\mathcal{C}_{b}(E)$ of bounded continuous functions $f:E\to\mathbb{R}$, endowed with the supremum norm

[TABLE]

We use the standard notation $\mathbb{P}$ and $\mathbb{E}$ for the distribution and the corresponding expectation on the usual path space

[TABLE]

The measurable structure on $\Omega$ is given by the Borel $\sigma$-algebra induced by the Skorokhod topology (see [35], Chapter 3). If we want to emphasize a particular initial condition $x\in E$ or distribution $\mu\in\mathcal{P}(E)$ of the process we write $\mathbb{P}_{x}$ and $\mathbb{E}_{x}$, or $\mathbb{P}_{\mu}$ and $\mathbb{E}_{\mu}$, respectively. The semigroup $P(t)$ acts on bounded continuous functions $f$ and probability measures $\mu\in\mathcal{P}(E)$ via

[TABLE]

where the latter provides a weak characterization of the distribution $\mu P(t)$ at time $t\geq 0$. Here and in the following we use the common notation $\mu(f)$ for expectations of $f\in\mathcal{C}_{b}(E)$ w.r.t. measures $\mu$ on $E$.

Using the Hille-Yosida Theorem (see e.g. [36], Chapter 3), it is possible to associate to the above Feller process an infinitesimal generator $\mathcal{L}$ acting on a dense subset $\mathcal{D}\subset\mathcal{C}_{b}(E)$ so that

[TABLE]

for all $f\in\mathcal{D}$ and $t\geq 0$.

In this work, we restrict ourselves to nonexplosive pure jump Feller processes. We denote by $\lambda(x)$ the escape rate from state $x\in E$ and the target state is chosen with the probability kernel $p(x,dy)$, so that the overall transition rate is

[TABLE]

for $(x,y)\in E^{2}$. We assume $\lambda:E\to[0,\infty)$ to be a strictly positive, bounded and continuous function and $x\mapsto p(x,A)$ to be a continuous function for every $A\in\mathcal{B}(E)$. Under these assumptions, the pure jump process possesses an infinitesimal generator [37, p. 162] with full domain $\mathcal{D}=\mathcal{C}_{b}(E)$ given by

[TABLE]

Along with jump processes on continuous spaces such as continuous-time random walks on $\mathbb{R}^{d}$ (see e.g. [38]), this setting includes in particular any finite-state continuous-time Markov chain. Typical compact examples we have in mind are given by stochastic particle systems on $E=S^{\Lambda}$, with finite local state space $S$ and lattice $\Lambda$ which can be finite or countably infinite. These include spin systems with $S=\{-1,1\}$ or exclusion processes with $S=\{0,\,1\}$, in which particles can jump only onto empty sites. Stochastic particle systems such as zero-range processes with $S=\mathbb{N}_{0}$ are locally compact as long as the lattice $\Lambda$ is finite (see e.g. [39] for details).

We will study Feynman-Kac models associated to the jump process by tilting its generator with a diagonal part or potential, which arise in many applications including dynamic large deviations, as explained in detail in Section 5.

Lemma 2.1**.**

Consider a potential function $\mathcal{V}\in C_{b}(E)$ and the tilted generator

[TABLE]

Then the family of operators \big{(}P^{\mathcal{V}}(t):t\geq 0\big{)} with $P^{\mathcal{V}}:C_{b}(E)\to C_{b}(E)$, defined as the solution to the backward equation

[TABLE]

*for all $f\in C_{b}(E)$, forms a non-conservative semigroup, the so-called Feynman-Kac semigroup, and $\mathcal{L}^{\mathcal{V}}$ is its infinitesimal generator in the sense of the Hille-Yosida Theorem. *

Proof.

See [36], Theorem 3.47.

∎

In order to control the asymptotic behaviour of $P^{\mathcal{V}}(t)$, we make the following assumption, which closely resembles [28, Assumption 1], on asymptotic stability.

Assumption 2.2** (Asymptotic Stability).**

The spectrum of $\mathcal{L}^{\mathcal{V}}=\mathcal{L}+\mathcal{V}$ (2) is bounded by a principal eigenvalue $\lambda_{0}$. Moreover, $\lambda_{0}$ is associated to a positive eigenfunction $r\in\mathcal{C}_{b}(E)$ and an eigenmeasure $\mu_{\infty}\in\mathcal{P}(E)$. Finally, there exist constants $\alpha>0$ and $\rho\in(0,1)$ such that

[TABLE]

for every $t\geq 0$ and $f\in\mathcal{C}_{b}(E)$.

Asymptotic stability is for example guaranteed for all irreducible, finite-state continuous-time Markov chains which necessarily have a spectral gap. For alternative sufficient conditions implying asymptotic stability in a more general context including continuous state spaces, see Appendix A.

We introduce the measures $\nu_{t,\mu_{0}}$ for any general initial distribution $\mu_{0}\in\mathcal{P}(E)$ and $t\geq 0$, defined by

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$. In the literature [21], $\nu_{t}$ is known as the unnormalised t-marginal Feynman-Kac measure. Applying Lemma 2.1, we can see that $\nu_{t}$ solves the evolution equation

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$, $t\geq 0$ and $\mu_{0}\in\mathcal{P}(E)$. The measures with which one can most naturally associate a process are the corresponding normalised t-marginal Feynman-Kac measures in $\mathcal{P}(E)$,

[TABLE]

defined for any $t\geq 0$ and $f\in\mathcal{C}_{b}(E)$.

Observe that, as a direct consequence of asymptotic stability (Assumption 2.2), there exist constants $\tilde{\alpha}\geq 0$ and $0<\rho<1$ such that for any $f\in\mathcal{C}_{b}(E)$,

[TABLE]

for any $t\geq 0$ and initial distribution $\mu_{0}\in\mathcal{P}(E)$. In particular $\mu_{t,\mu_{0}}$ converges weakly to $\mu_{\infty}$, as $t\to\infty$. Indeed, by definition of $\mu_{t,\mu_{0}}$ (7) and then by asymptotic stability (Assumption 2.2),

[TABLE]

for any $t>-\log\alpha/\log\rho$ and for some constant $\alpha>0$. This gives the bound (8) for any $t$ large enough. Increasing $\tilde{\alpha}$ accordingly to ensure that the bound holds also for small $t$, we obtain (8) for any $t\geq 0$.

For simplicity, in the rest of this article the initial distribution $\mu_{0}$ is fixed and we write $\mu_{t}$ (resp. $\nu_{t}$) instead of $\mu_{t,\mu_{0}}$ (resp. $\nu_{t,\mu_{0}}$).

2.2 McKean Interpretations

Now, we want to outline the evolution of the time-marginal distribution $\mu_{t}$ in terms of interacting jump-type infinitesimal generators. The content presented in the rest of this section is based on the works of Del Moral and Miclo [21, 22, 20]. In this established framework it is possible to define generic Markov processes with time marginals $\mu_{t}$ and then use Monte Carlo sampling techniques to approximate those marginals.

Lemma 2.3**.**

For every $f\in\mathcal{C}_{b}(E)$ and $t\geq 0$, the normalised $t$-marginal $\mu_{t}$ (7) solves the non-linear evolution equation

[TABLE]

Proof.

Using the evolution equation (6) of $\nu_{t}$, we see that

[TABLE]

∎

The evolution equation (10) results from the unique decomposition of the non-conservative generator $\mathcal{L}+\mathcal{V}$ into a conservative and a diagonal part given by the potential $\mathcal{V}$. The latter, together with the normalization of $\nu_{t}$, leads to the nonlinear second part in (10) which we want to rewrite to be in the form of another infinitesimal generator, that we denote by $\widetilde{\mathcal{L}}_{\mu_{t}}$. Since (10) is non-linear in $\mu_{t}$, this depends itself on the current distribution such that

[TABLE]

for every $\mu\in\mathcal{P}(E)$ and $f\in\mathcal{C}_{b}(E)$. The choice of the non-linear generator $\widetilde{\mathcal{L}}_{\mu}$ is not unique, leading to various representations of the form

[TABLE]

where $\widetilde{W}(x,y)\mu(dy)$ is the overall transition kernel of $\widetilde{\mathcal{L}}_{\mu}$ and depends on the current distribution $\mu$.

Lemma 2.4** (Sufficient conditions).**

An infinitesimal generator in the form (12) satisfies condition (11) if and only if

[TABLE]

for all $\mu\in\mathcal{P}(E)$ and $x\in E$. In particular, a sufficient condition on $\widetilde{\mathcal{L}}_{\mu}$ (12) for (11) to hold is

[TABLE]

for all $x,y\in E$.

Proof.

It is enough to observe that

[TABLE]

∎

Combining $\widetilde{\mathcal{L}}_{\mu}$ with the linear part $\mathcal{L}$ of (10) into a so-called McKean generator on $\mathcal{C}_{b}(E)$,

[TABLE]

the evolution equation (10) can be written as

[TABLE]

for every $f\in\mathcal{C}_{b}(E)$ and $t\geq 0$. Therefore, the normalized Feynman-Kac marginal $\mu_{t}$ can be interpreted as the law of a Markov process \big{(}\overline{X}_{t}\,:\,t\geq 0\big{)} on $E$, associated to the family of generators \big{(}\overline{\mathcal{L}}_{\mu_{t}}:t\geq 0\big{)}. This process is also known as a McKean representation of the process associated to the Feynman-Kac measure $\mu_{t}$, and it is non-linear and in particular time-inhomogeneous. This can be formulated using the propagator

[TABLE]

for all $0\leq t\leq T$, which follows directly from the definition of $\mu_{t}$ (7) and the semigroup characterizing the time evolution for $\nu_{t}$ (3).

While the time evolution of $\mu_{t}$ is uniquely determined by (10) and therefore independent of the choice of (13), Lemma 2.4 leads to various different McKean representations of the form (12) (see e.g. [34, 28]), that can be characterized by the operator $\widetilde{W}$. One common choice related to algorithms in [1, 2] is

[TABLE]

where $c\in\mathbb{R}$ is an arbitrary constant, and we use the standard notation $a^{+}=\max\{0,a\}$ and $a^{-}=\max\{0,-a\}$ for positive and negative part of $a\in\mathbb{R}$.

One other possible representation of (12) we want to mention explicitly here is given by

[TABLE]

This corresponds to a pure jump process on $E$ in which every jump strictly increases the value of the potential $\mathcal{V}$ in contrast to the previous representation (15). We will see in the next section that $\mathcal{V}$ can be interpreted as a fitness potential for the overall process. Further McKean representations of (10) are discussed in [34], here we focus on cloning algorithms which are based on (15).

3 Interacting Particle Approximations

Independent of the particular representation, the rates of the McKean process $(\overline{X}_{t}:t\geq 0)$ depend on the distribution $\mu_{t}$ itself, which is in general not known. A standard approach is to sample such processes through particle approximations [23], which involve running, in parallel, $N$ copies or clones $\xi_{t}:=(\xi_{t}^{1},\dots,\xi_{t}^{N})\in E^{N}$ of the process (called particles), and then approximating $\mu_{t}$ by the empirical distribution $m(\xi_{t})$ of the realizations. For any $\underline{x}\in E^{N}$ the latter is defined as

[TABLE]

We write $\overline{L}^{N}$ for the infinitesimal generator of an $N$-particle system $\xi_{t}$ and also call this an IPS generator, and denote the associated empirical distribution as

[TABLE]

We denote by

[TABLE]

the standard carré-du-champ operator associated to the generator $\overline{L}^{N}$.

3.1 A general convergence result

The full dynamics can be set up in various different ways such that $\mu_{t}^{N}\to\mu_{t}$ converges in an appropriate sense as $N\to\infty$ for any $t\geq 0$. Theoretical convergence results can be obtained under the following assumptions, which are fulfilled by standard mean field particle approximations (as shown in Section 3.2) and cloning algorithms (Section 3.3).

Assumption 3.1**.**

[TABLE]

Remark*.*

Test functions of the form

[TABLE]

describe mean-field observables averaged over the particle ensemble which are generally of most interest, e.g. for the estimator (79) of the SCGF it is sufficient to consider such functions, as shown in Section 5.2. In general the goal is to approximate $\mu_{t}(f)$ for a given $f\in\mathcal{C}(E)$, so it is natural to set up the auxiliary particle approximation in a permutation invariant way and use mean-field observables.

To better understand the above assumptions, recall that the carré du champ of an interacting particle system is a quadratic operator associated to the fluctuations of the process, whereas the generator determines the expected behaviour of the observables $F({\xi}_{t})$. Thus, Assumption 3.1 implies that trajectories of mean-field observables in a particle approximation coincide in expectation with average trajectories of the McKean representation they are based on (19a), and concentrate on their expectation with diverging $N$ (19b). We include the operators $G_{\mu}$ explicitly in (19b), because it allows the condition to be stated in a convenient form and we anticipate it being useful in further analysis. Condition (19c) assures that at any given time only a bounded number of particles can change their state, which is a mild technical assumption, necessary to allow the application of Lemma 4.1 in the proof of the $L^{p}$ error estimates.

Theorem 3.2**.**

Consider a sequence of particle approximations satisfying Assumption 3.1 with empirical distributions $\mu_{t}^{N}$ (18). Under Assumption 2.2, for every $p\geq 2$ there exists a constant $c_{p}>0$ independent of $N$ and $T$ such that

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$. Furthermore, there exists a constant $c^{\prime}>0$ independent of $N$ and $T$ such that

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$ and $N\in\mathbb{N}$ large enough.

Remark*.*

The constants $c_{p}$ and $c^{\prime}$ depend on the Feynman-Kac model of interest, on the choice of the McKean model and on the considered interacting particle approximation.

The proof, presented in Section 4, is an adaptation of the results in [28] and makes use of the propagator (14) of $\mu_{t}$ and the martingale characterization of $(\xi_{t}:t\geq 0)$.

Remark*.*

Observe that, by Markov’s inequality, Theorem 3.2 implies

[TABLE]

for every $\varepsilon,\,t>0$, $f\in\mathcal{C}_{b}(E)$, $N\geq K$ and $p\geq 2$, where $c_{p}>0$ does not depend on $N$. In particular, considering $p>2$, we can see that

[TABLE]

as $N\to\infty$, for any $f\in\mathcal{C}_{b}(E)$, by a Borel-Cantelli argument. The existence of a countable determining class allows this to be further strengthened to the almost sure convergence of $\mu_{t}^{N}$ to $\mu_{t}$ in the weak topology (see, for example, [40, Theorem 4]).

It is important to clarify that the estimators of the Feynman-Kac distribution $\mu_{t}$ given by the empirical measures $\mu^{N}_{t}$ usually have a bias, i.e. $\mathbb{E}[\mu^{N}_{t}(f)]\neq\mu_{t}(f)$ for $f\in\mathcal{C}_{b}(E)$, which vanishes only asymptotically, as illustrated in Theorem 3.2. This arises from the non-linear time evolution of $\mu_{t}$. However, it is straightforward to derive unbiased estimators of the unnormalized measures $\nu_{t}$ (5), as shown by the following result.

Proposition 3.3** (Unbiased Estimators).**

Consider a sequence of particle approximations satisfying (19a) and initial condition (19d), with empirical distributions $\mu_{t}^{N}$ (18). Then, the unnormalized empirical measure

[TABLE]

is an unbiased estimator of the unnormalized $t$-marginal $\nu_{t}$ (5), i.e.

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$.

Proof.

First observe that \mathbb{E}\big{[}\nu_{0}^{N}(f)\big{]}=\nu_{0}(f). Indeed, $\nu^{N}_{0}(f)=\mu^{N}_{0}(f)$ is the average of $N$ i.i.d. random variables with law $f_{\#}\mu_{0}$, and $\mu_{0}$ corresponds to the initial distribution of $\nu_{t}=\nu_{t,\mu_{0}}$ (5).

Note that \mathbb{E}\big{[}\nu_{t}^{N}(f)\big{]} satisfies the evolution equation

[TABLE]

Moreover, by assumption (19a) and using the characterization of $\overline{\mathcal{L}}_{\mu}$ (11)-(13), we have

[TABLE]

Inserting into (24), this simplifies to

[TABLE]

Since $\mathcal{L}+\mathcal{V}$ also generates the time evolution of $\nu_{t}(f)$ (6), a simple Gronwall argument with \mathbb{E}\big{[}\nu_{0}^{N}(f)\big{]}=\nu_{0}(f) gives (23).

∎

A generic version of interacting particle systems, directly related to the above McKean representations has been studied in the applied probability literature in great detail [23, 28], providing quantitative control on error bounds for convergence. After reviewing those results in the next subsection, we present a different approach taken in the theoretical physics literature under the name of cloning algorithms [1, 5], which provides some computational advantages but lacks general rigorous error control so far [15, 16].

3.2 Mean Field Particle Approximation

The most basic particle approximation is simply to run the McKean dynamics in parallel on each of the particles, replacing the distribution $\mu_{t}$ by the empirical measure. Formally, the mean field particle model $(\xi_{t}:t\geq 0)$ with $\xi_{t}=(\xi^{i}_{t}:i=1,\dots,N)$ associated to a McKean generator $\overline{\mathcal{L}}_{\mu_{t}}$ (13), is a Markov process on $E^{N}$ with homogeneous infinitesimal generator $\overline{L}^{N}$ defined by

[TABLE]

for any $F\in\mathcal{C}_{b}(E^{N})$. Here $\overline{\mathcal{L}}_{m(\underline{x})}^{(i)}$ denotes the McKean generator $\overline{\mathcal{L}}_{m(\underline{x})}$ (13) acting on the function $x^{i}\mapsto F(x^{1},\dots,x^{i},\dots,x^{N})$, where the dependence on $\mu$ has been replaced by the empirical distribution $m(\underline{x})$.

In analogy to the decomposition $\overline{\mathcal{L}}_{\mu}=\mathcal{L}+\widetilde{L}_{\mu}$ in (13), the generator (25) can be decomposed as $\overline{L}^{N}=L^{N}+\widetilde{L}^{N}$ with

[TABLE]

where $\mathcal{L}^{(i)}$ and $\widetilde{\mathcal{L}}_{m(\underline{x})}^{(i)}$ stand respectively for the operators $\mathcal{L}$ and $\widetilde{\mathcal{L}}_{m(\underline{x})}$ acting on the function $x^{i}\mapsto F(\underline{x})$, i.e. only on particle $i$.

Moreover, using representation (12) for $\widetilde{\mathcal{L}}_{\mu}$, we can write

[TABLE]

with $\underline{x}^{i,y}:=(x_{1},\dots,x_{i-1},y,x_{i+1},\dots,x_{N})$, which introduces an interaction between the particles. In this decomposition, (26) generates the so-called mutation dynamics, where the particles evolve independently under the dynamics given by the infinitesimal generator $\mathcal{L}$ of the original process, whereas (27) generates the selection dynamics, which leads to mean-field interactions between particles. With (28) the state of particle $i$ gets replaced by that of particle $j$ with rate $\frac{1}{N}\widetilde{W}(x_{i},\,x_{j})$. The total selection rate in the particle approximation is $\frac{1}{N}\sum_{i,j=1}^{N}\widetilde{W}(x_{i},\,x_{j})$, and depends on the McKean representation, in particular the choice of $\widetilde{\mathcal{L}}_{\mu}$ in (12).

From general practical experience it is favourable to minimize the total selection rate in order to improve the estimator’s asymptotic variance; it’s widely understood in the SMC literature that eliminating unnecessary selection events can significantly improve estimator variances, see, for example, [21, Section 7.2.1, 7.4.2] and [41]. For mean-field particle approximations this suggests that (16) is preferable to (15) since

[TABLE]

for all $x,y\in E$ and $c\in\mathbb{R}$. In view of Lemma 2.4, minimizing the total selection rate pertains to maximizing $\sum_{i=1}^{N}\mathcal{V}(x_{i})$, and $\mathcal{V}$ can be interpreted as a fitness function. With (16) every selection event therefore increases the fitness of the particle ensemble, which is not necessarily the case with (15), and there are even more optimal choices than (16) in that sense as discussed in [34]111As a side remark, the mutation part of the McKean dynamics (which is fixed for mean-field particle particle approximations by (26)), can naturally also decrease the fitness of the ensemble.. On the other hand, depending on the particular application, implementing particle approximations with lower total selection rate could be computationally more expensive, leading to a trade-off in lower values for $N$ to be accessible in practice. This is discussed in [34] for a particular example, and is not the subject of this paper.

In order to motivate the choice of the cloning algorithm in the next subsection which is based on the selection rates (15), we note that one can write (27) as

[TABLE]

using a change of summation indices in the second term. With the above discussion this can be interpreted as follows: If particle $i$ is less fit than level $c$ it is killed and replaced by a uniformly chosen particle $j$, and if it is fitter than $c$ it cloned, replacing a uniformly chosen particle $j$.

Observe that, by definition of $\overline{L}^{N}$ (25), for any function $F$ on $E^{N}$ of the form $F(\underline{x})=m(\underline{x})(f)$, with $f\in\mathcal{C}_{b}(E)$, we have that

[TABLE]

thus conditions (19a)-(19b) are satisfied.

Analogous relations hold also for the individual mutation and cloning parts of the generator. Since generators are linear, the identity (30) is immediate. The carré du champ (31) is quadratic in $F$, but off-diagonal terms in the corresponding double sum turn out to vanish in a straightforward computation, leading to the additional factor $1/N$. Furthermore, by construction, for almost every realization $\xi_{t}$, $t>0$, of the mean field particle approximation, there exists at most one particle $i$ such that $\xi^{i}_{t}\neq\xi^{i}_{t-}$, thus condition (19c) is satisfied with $K=1$. Therefore, Theorem 3.2 holds and provides $L^{p}$-error and bias estimates of order $1/\sqrt{N}$ and $1/N$ respectively, in accordance with already established results, e.g. in [28, 23, 22].

3.3 The Cloning Algorithm

Cloning algorithms have been proposed in the theoretical physics literature [1, 2] for evaluating large deviation functions associated to Markov processes similar to the mean field system (25), using the same mutation dynamics. While selection and mutation events are independent in the latter due to the additive structure of $\overline{L}^{N}$ in (26) and (27), in cloning algorithms both are combined to reduce computational cost. We focus the exposition on a variant of the algorithm proposed in [2], but other continuous-time versions can be analysed analogously. This cloning algorithm is constructed from the McKean model $\overline{\mathcal{L}}_{\mu}$ (13) with selection rates \widetilde{W}_{c}(x,y)=\big{(}\mathcal{V}(x)-c\big{)}^{-}+\big{(}\mathcal{V}(y)-c\big{)}^{+} as in (15), and we denote the associated McKean generator by

[TABLE]

We will use in particular the killing/cloning interpretation introduced in (29). We recall that the overall escape rate and probability kernel of the original dynamics $\mathcal{L}$ are denoted respectively by $\lambda(x)$ and $p(x,dy)$.

The infinitesimal description of the cloning algorithm as a continuous-time Markov process on the state space $E^{N}$ is given by the generator

[TABLE]

for any $F\in\mathcal{C}_{b}(E^{N})$ and $\underline{x}\in E$. Here $\mathcal{N}$ is the set of all subsets of $N$ particle indices, $\underline{x}^{A,w}$ denotes the vector $(z_{1},\dots,z_{N})\in E^{N}$, with

[TABLE]

and, similarly, $\underline{x}^{A,w;\,i,y}$ denotes the vector $(z_{1},\dots,z_{N})\in E^{N}$ with

[TABLE]

for any $w,y\in E$. Cloning events are now coupled with mutation, and if $\mathcal{V}(x_{i})>c$, a non-empty set $A$ of particles is chosen at random from the ensemble with probability $\pi_{x_{i}}(A)$ and every particle $j\in A$ is replaced by a clone of $i$, before particle $i$ mutates to a new state $y\in E$. If $\mathcal{V}(x_{i})\leq c$ we set $\pi_{x_{i}}(A)=\delta_{A,\emptyset}$, so that no cloning occurs. Further properties of the cloning distribution $\pi_{x}(\cdot)$, which is the main distinctive feature of this algorithm, are discussed below. The killing part in the second line runs independently and remains unchanged from (29). The algorithm is often applied in situations with $\mathcal{V}(x)\geq c$ for all $x\in E$ (in particular also with $c=0$), leaving cloning coupled with mutation as the only selection events.

In order to simplify the presentation, we make some further assumptions on $\pi_{x}(A)$, which are all satisfied by common choices in the theoretical physics literature. The probability of choosing a set $A$ depends only on its size $|A|$ and not on its elements, i.e. for any $x\in E$

[TABLE]

Denote the mean and second moment of this distribution by

[TABLE]

Of course, $\pi_{x,.}$ and its moments also depend on $N$ and $c$, which we omit in the notation for simplicity. In order to ensure that the third condition in Assumption 3.1, namely (19c), is satisfied, we assume that the support of $\pi_{x,.}$ is uniformly bounded in $N$, i.e.

[TABLE]

Note that this implies that also $M(x)$ and $Q(x)$ are uniformly bounded, i.e. $M,Q\in\mathcal{C}_{b}(E)$. We further assume $N\geq K$, i.e. $N$ is large enough so that the process (33) is well defined.

The most common choice in the physics literature (see, e.g., the recent summary in [7]) for the distribution $\pi_{x,.}$ is

[TABLE]

This corresponds to a binary distribution on the two integers nearest to the prescribed mean, and minimizes the second moment $Q$ of the distribution for a given mean. Note that if $M(x)$ is an integer, $\pi_{x,n}=\delta_{n,M(x)}$ concentrates, which includes the case $M(x)=0$.

The next two results assure respectively that condition (19a) and condition (19b) in Assumption 3.1 are satisfied for the cloning algorithm, so we can apply Theorem 3.2. The only condition is to choose $M(x)$ such that each particle $i$ produces on average \big{(}\mathcal{V}(x_{i})-c\big{)}^{+} clones per unit time, in accordance with the second term in (29).

Proposition 3.4**.**

Consider the cloning generator $\overline{L}^{N}_{c}$ (33) with $\pi_{x}(A)$ as in (3.3) and (36), such that the mean of the cloning size (35) is

[TABLE]

and $\sup_{x\in E}M(x)<\infty$. Then, for any test function of the form $F(\underline{x})=m(\underline{x})(f)$, with $f\in\mathcal{C}_{b}(E)$ and $N$ large enough, we get

[TABLE]

where $\overline{\mathcal{L}}_{m(\underline{x}),c}$ is the McKean generator given in (32).

Remark*.*

Note that $\sup_{x\in E}M(x)<\infty$ is essential for (36) and (19c), and a simple sufficient condition is for the escape rates to be uniformly bounded below, i.e. $\inf_{x\in E}\lambda(x)>0$.

Proof.

We start by considering the first term in the expression of $\overline{L}^{N}_{c}$ (33). Observe that with $F(\underline{x})=m(\underline{x})(f)$,

[TABLE]

Thus, we can write

[TABLE]

Moreover, by (3.3), we have that, for any $j\in\{1,\dots,N\}$,

[TABLE]

Therefore,

[TABLE]

Thus, (33) can be rewritten as

[TABLE]

by changing summation variables in the cloning term and using (32).

∎

Proposition 3.5**.**

Let $\overline{L}^{N}_{c}$ be a cloning generator satisfying the conditions in Proposition 3.4. Then, for any test function of the form $F(\underline{x})=m(\underline{x})(f)$, with $f\in\mathcal{C}_{b}(E)$,

[TABLE]

as $N\to\infty$, where $\|\Delta_{N}(.,\,f)\|\,\leq\,C\,\frac{\|f\|^{2}}{N^{2}}$ for some constant $C>0$ independent of $f$ and $N$, and

[TABLE]

with

[TABLE]

and

[TABLE]

Remark*.*

Due to the linearity of the generator, the combined mutation/cloning events in the cloning algorithm can be decomposed easily, which leads to extra terms only in the quadratic carré du champ. In the expression of the operator $G_{\mu}$ (42), the term

[TABLE]

is due to the dependence between mutation and cloning dynamics and its sign is not known a priori. Whereas, the term \lambda(x)\,\big{(}Q(x)-M(x)\big{)}\cdot\big{(}\ell_{\mu}f(x)\big{)}^{2} arises from the dependence between clones (since multiple cloning events are allowed at the same time) and is always non-negative. In particular, in any setting in which there is at most one clone per event, i.e. when $Q(x)=M(x)$, the term vanishes. Furthermore, minimizing $Q(x)$ as in (37) for given $M(x)$ (38) leads to the best bound on the carré du champ and convergence properties of the algorithm.

Proof.

Consider the carré du champ of $\overline{L}_{c}^{N}$,

[TABLE]

Using (39), the first term can be decomposed as

[TABLE]

where with (40) and (41) the last line can be rewritten as

[TABLE]

Substituting in the expression of the carré du champ $\Gamma_{\overline{L}_{c}^{N}}$, we obtain

[TABLE]

The first line in (43) is simply

[TABLE]

Now, considering the second line of (43), we can write

[TABLE]

since

[TABLE]

for every $j,k\in\{1,\dots,N\}$ such that $j\neq k$.

Recalling that \lambda(x)M(x)=\big{(}\mathcal{V}(x)-c\big{)}^{+}, exchanging summation indices and combining with the third line of (43), we see that

[TABLE]

Moreover,

[TABLE]

with

[TABLE]

for all $\underline{x}\in E$, for some constant $C>0$, since $M(x)$ and $Q(x)$ are bounded by condition (36) and $\lambda(x)$ is bounded by assumption. Combining all together, we obtain the statement.

∎

Proposition 3.4 and Proposition 3.5 show in particular that Assumption 3.1 is satisfied for cloning algorithms, hence Theorem 3.2 holds and provides bias and $L^{p}$ error bounds.

4 Proof of Theorem 3.2

This section is devoted to the proof of Theorem 3.2, which is an adaptation of the results presented by M. Rousset in [28]. Throughout this section we consider a generic sequence of IPS generators $(\overline{L}^{N})_{N\in\mathbb{N}}$ satisfying Assumption 3.1 for some McKean generator $\overline{\mathcal{L}}_{\mu}$ (13). Furthermore, we assume that the normalized Feynman-Kac measure $\mu_{t}$ is asymptotically stable, i.e. Assumption 2.2 holds.

The proof makes use of the propagator $\Theta_{t,T}$ of $\mu_{t}$ defined in (14), and the martingale characterization of $\overline{L}^{N}$. We denote by $\mathcal{C}_{b}^{0,1}(E\times\mathbb{R}^{+})$ the set of bounded functions $\varphi_{\cdot}$ such that $\varphi_{t}(\cdot)$ is continuous on $E$ for every $t\in\mathbb{R}^{+}$ and $\varphi_{\cdot}(x)$ has continuous time derivative for every $x\in E$. Following the standard martingale characterization of Feller-type Markov processes, using Itô’s formula and (19a) one can show that (see also [28], Proposition 3.3), for every $\varphi_{\cdot}\in\mathcal{C}_{b}^{0,1}(E\times\mathbb{R}^{+})$, the process

[TABLE]

is a local martingale. With (19b) its predictable quadratic variation is bounded by

[TABLE]

for some constant $C\geq 0$ independent of $\varphi$ and $N$, and with (19c) jumps are bounded by

[TABLE]

The following technical Lemma for martingales will play a central role in the proof of Theorem 3.2.

Lemma 4.1**.**

Let $\mathcal{M}$ be a locally square-integrable martingale with continuous predictable quadratic variation $\langle\mathcal{M}\rangle$, $\mathcal{M}_{0}=0$ and uniformly bounded jumps $\sup_{t}|\Delta\mathcal{M}_{t}|\leq a<\infty$. Then, for every $q\in\mathbb{N}_{0}$ and $T\geq 0$, there exists a constant $C_{q}>0$ such that

[TABLE]

Proof.

See [28], Lemma 6.2.

∎

4.1 Properties of the normalized propagator

Lemma 4.2**.**

For any test function $f\in\mathcal{C}_{b}(E)$ and $0\leq t\leq T$, we have for the normalized propagator (14)

[TABLE]

Proof.

See [28], p. 836. The idea of the proof is to substitute (3) into the time derivative of $\Theta_{t,T}f$ (14).

∎

Lemma 4.3**.**

Under Assumption 2.2 on asymptotic stability, for any $0\leq t\leq T$ and $n\in\mathbb{N}$ and $f\in\mathcal{C}_{b}(E)$, there exists a constant $\beta>0$ such that

[TABLE]

Moreover, for any $\overline{f}:=f-\mu_{T}(f)$, there exists some $0<\rho<1$, such that

[TABLE]

Proof.

The proof can be found in [28, Lemma 5.1] and the result is due to the asymptotic stability of the Feynman-Kac model.

∎

Observe that, applying Lemma 4.2 to the martingale characterization (44) of $\overline{L}^{N}$, we obtain

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$, where the last equality follows by the characterization (11) of McKean models. By (47), we obtain the stochastic differential equation

[TABLE]

Moreover, applying Lemma 4.3 to the predictable quadratic variation (45), we obtain that almost surely,

[TABLE]

where $\overline{G}(f,f)=\sup_{\mu\in\mathcal{P}(E)}G_{\mu}(f,f)$ .

Note that Equation (47) for centered test functions $\overline{f}=f-\mu_{T}(f)$ can be rewritten as

[TABLE]

The martingale characterization (47)-(50) will be the key element in the proof of Theorem 3.2.

4.2 $L^{p}$ and bias estimates

Define

[TABLE]

with $\mu\in\mathcal{P}(E)$ and $0\leq t\leq T$. Observe that the measure $\Phi_{t,T}(\mu)$ can be also rewritten in terms of $\Theta_{t,T}$ (14) as

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$. To prove Theorem 3.2, we consider the decomposition

[TABLE]

for any $0\leq t\leq T$. The proof is structured as follows:

• •

In Lemma 4.4, we bound the first term of the decomposition under Assumptions 2.2 and 3.1;

• •

In Lemma 4.5, we bound the second term under Assumption 2.2;

• •

In Lemma 4.6, we combine Lemma 4.4 and Lemma 4.5 to obtain $L^{p}$-error estimates of order $1/N^{\delta/2}$, for some $\delta\in(0,1)$;

• •

Finally, from Lemma 4.6 we derive, by iteration, $L^{p}$ estimates of order $1/\sqrt{N}$, as presented in Theorem 3.2.

Lemma 4.4**.**

Consider a sequence of particle approximations satisfying Assumption 3.1 with empirical distributions $\mu_{t}^{N}$ (18). Under Assumption 2.2 on asymptotic stability, for any $p\geq 2$ there exists a constant $c_{p}>0$ such that

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$ and $0\leq t\leq T$.

Proof.

This is an adaptation of the first part of the proof of Lemma 5.3 in [28]. First, consider

[TABLE]

with $0\leq t_{1}\leq t_{2}$. Observe that, by the stochastic differential equation (48), we can write

[TABLE]

for any $t\leq s\leq T$. Therefore,

[TABLE]

Fixing $0\leq t\leq T$, the process

[TABLE]

with $t\leq\tau\leq T$, as the integral of a progressively measurable process with respect to a local martingale, is itself a local martingale with predictable quadratic variation given by

[TABLE]

and jumps bounded by

[TABLE]

by Assumption (19c) on bounded jumps, (46) and Lemma 4.3.

Moreover, with (52), we can write

[TABLE]

where the last equality follows by (55). Noting that (A_{t}^{T})^{-1}\leq\exp\big{(}2(T-t)\cdot\|\mathcal{V}\|\big{)} by definition (54), we get

[TABLE]

By Lemma 4.1, we have that, for any $q\in\mathbb{N}_{0}$,

[TABLE]

where the last inequality follows by (49). Therefore, for $p=2^{q+1}$, $q\in\mathbb{N}_{0}$, we get

[TABLE]

By Jensen’s inequality, this bound holds for any $p\geq 2$. Applying this to inequality (56), we obtain the result.

∎

Lemma 4.5**.**

Under Assumption 2.2 on asymptotic stability with constants $\alpha>0$ and $\rho\in(0,1)$, we have that for any $p\geq 2$ and any $0\leq t\leq T$ such that $T-t\geq(\log\varepsilon\,-\,\log\alpha)/\log\rho$ for some $\varepsilon\in(0,1)$, the following bound holds

[TABLE]

Furthermore, when $t=0$, there exists a constant $C_{p}>0$ depending on $p$ such that

[TABLE]

Proof.

By definition (51) of $\Phi_{t,T}$, for any $\eta\in\mathcal{P}(E)$ and $\lambda\in\mathbb{R}$ we have

[TABLE]

Taking $\lambda$ to be the principal eigenvalue of $\mathcal{L}+\mathcal{V}$, using Assumption 2.2 on asymptotic stability and the basic fact $\eta(1)=1$, we can write

[TABLE]

Therefore, for $T-t\geq(\log\varepsilon\,-\,\log\alpha)/\log\rho$, for some $\varepsilon\in(0,1)$, we have

[TABLE]

and similarly

[TABLE]

Therefore,

[TABLE]

Now, for $t=0$, observe that

[TABLE]

Using the basic fact $1=\mu_{0}(\Theta_{0,T}(1))$, to conclude it is enough to observe that, for any $f\in\mathcal{C}_{b}(E)$,

[TABLE]

with $C_{p}>0$ constant depending on $p$. Indeed, with (19d) at time $t=0$, $\mu^{N}_{0}(f)$ is the sum of $N$ i.i.d. random variables with law $f_{\#}\mu_{0}$. Inequality (57) is then a direct application of Marcinkiewicz-Zygmund/BDG inequalities for i.i.d. variables.

∎

Lemma 4.6**.**

Consider a sequence of particle approximations satisfying Assumption 3.1 with empirical distributions $\mu_{t}^{N}$ (18). Under Assumption 2.2, there exists $\delta\in(0,1)$ such that for any $p\geq 2$ there exist $c_{p}>0$ such that

[TABLE]

for any $N\in\mathbb{N}$ large enough.

Proof.

Recalling decomposition (53), where the first term is estimated in Lemma 4.4 and the second in Lemma 4.5, and using the basic fact $T-t\leq e^{T-t}$, we obtain

[TABLE]

taking $t=0$, and

[TABLE]

taking $0\leq t\leq T$ such that $T-t$ is large enough.

The idea is to find $t\geq 0$ and $\varepsilon\in(0,1)$ such that

[TABLE]

Recalling that $\log\rho<0$, the solution is given by

[TABLE]

provided $T\geq{\log N}/({8\|\mathcal{V}\|+1-2\log\rho})$ to ensure that $t\geq 0$. Also observe that for $N$ large enough, $T-t$ satisfies the conditions in Lemma 4.5.

Otherwise, in case $T<{\log N}/({8\|\mathcal{V}\|+1-2\log\rho})$, we consider the bound (58) instead, and we obtain

[TABLE]

with

[TABLE]

Taking $\delta=\min\{\varepsilon,\,\overline{\varepsilon}\}$ the result follows from observing that

[TABLE]

for $t=0$ and $T$ at most of order $\log N$ as above, or for $t\geq 0$ given by (60).

∎

Proof of Theorem 3.2.

We denote

[TABLE]

in accordance with Rousset [28], Section 5.2. Using (50), we have

[TABLE]

with $\overline{f}=f-\mu_{T}(f)$ for any $f\in\mathcal{C}_{b}(E)$.

First, observe that, similarly to (57), we have

[TABLE]

for some constant $C_{p}>0$ depending on $p$. Moreover, by Lemma 4.1 and bound (49), we get with another $p$-dependent constant

[TABLE]

Finally, writing

[TABLE]

and using Hölder’s inequality, we get

[TABLE]

by Lemma 4.3. Using the fact that

[TABLE]

for centered test functions, and applying the Cauchy-Schwarz inequality, we get

[TABLE]

Combining all together, we obtain

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$ and $T\geq 0$. In particular,

[TABLE]

for any $p\geq 2$. Applying Lemma 4.6, we get

[TABLE]

for any $k\in\mathbb{N}$, by iteration of (62). Thus, we can conclude

[TABLE]

This proves the $L^{p}$-error estimate (20).

We conclude by proving the bias estimate (21). By Equation (50), we have

[TABLE]

By (61) for $p=1$, we obtain

[TABLE]

∎

5 Interacting particle approximations for dynamic large deviations

5.1 Large deviations and Feynman-Kac models

Dynamic large deviations of continuous-time jump processes are a common application area of cloning algorithms [1, 2]. For a given process $(X_{t}:t\geq 0)$ with bounded rates $W(x,dy)=\lambda(x)p(x,dy)$ (1) and path space $\Omega$ as outlined in Section 2, we consider a time-additive observable $A_{T}:\Omega\to\mathbb{R}$, taken to be a real measurable function of the paths of $X_{t}$ over the time interval $[0,T]$ of the form [42]

[TABLE]

Here $g\in\mathcal{C}_{b}(E^{2})$ is such that $g(x,x)=0$, for any $x\in E$, and $h\in\mathcal{C}_{b}(E)$, with $\omega\in\Omega$ a realization of $(X_{t}:t\geq 0)$. Note that $A_{T}$ is well defined since the bound on $\lambda(x)$ implies that the process does not explode and the first sum contains almost surely only finitely many non-zero terms for any $T\geq 0$.

More precisely, we are interested in studying the limiting behaviour, as $T\to\infty$, of the family of probability measures $\mathbb{P}_{\mu_{0}}(A_{T}\in\,\cdot\,)=\mathbb{P}_{\mu_{0}}\circ A_{T}^{-1}$ on $(\mathbb{R},\,\mathcal{B}(\mathbb{R}))$, where $\mu_{0}$ represents the initial distribution of the underlying process. This can be characterized by the large deviation principle (LDP) [43, 44], in terms of a rate function. We assume that an LDP with convex rate function $I$ holds, which can be written as

[TABLE]

for every $C\subseteq\mathbb{R}$ closed and $O\subseteq\mathbb{R}$ open. For the study of large deviations, a key role is played by the scaled cumulant generating function (SCGF)

[TABLE]

Indeed, if the rate function $I$ is convex and the limit $\Lambda_{k}$ in (64) exists and is finite for every $k\in\mathbb{R}$, then $I$ is fully characterized by the SCGF via Legendre duality (see [43], Theorem 4.5.10), i.e.

[TABLE]

The SCGF is also the object that can be numerically approximated by cloning algorithms [1, 2] and related approaches and our main aim in this section is to illustrate how our results on Feynman-Kac models can be applied here. Possible subtleties regarding the LDP are not our focus and we restrict ourselves to settings where $\Lambda_{k}$ exists and is finite. In the following we introduce the associated Feynman-Kac models in the notation that is established in this context.

Lemma 5.1**.**

For any $k\in\mathbb{R}$ the family of operators \big{(}P_{k}(t):t\geq 0\big{)} on $\mathcal{C}_{b}(E)$ defined by

[TABLE]

with $f\in\mathcal{C}_{b}(E)$, is well defined and it is a non-conservative semigroup, the so-called tilted semigroup.

Moreover, the infinitesimal generator associated with \big{(}P_{k}(t):t\geq 0\big{)}, in the sense of the Hille-Yosida Theorem, can be written in the form

[TABLE]

for $f\in\mathcal{C}_{b}(E)$ and all $x\in E$, with $g$ and $h$ the bounded continuous functions which characterize $A_{T}$ via (63). In particular, the semigroup $P_{k}(t)$ satisfies the differential equations

[TABLE]

for all $f\in\mathcal{C}_{b}(E)$ and $t\geq 0$.

Proof.

See [42], Appendix A.1.

∎

Observe that, if the SCGF (64) is independent of the choice of the initial distribution $\mu_{0}$, it can be written in terms of the tilted semigroup as

[TABLE]

for all $x\in E$, moreover $\Lambda_{k}$ is the spectral radius of the generator $\mathcal{L}_{k}$ (see also (70) below). With Assumption 2.2 on asymptotic stability, $\Lambda_{k}$ is also the principal eigenvalue of $\mathcal{L}_{k}$ and there exists a probability measure $\mu_{\infty}=\mu_{\infty,k}\in\mathcal{P}(E)$222To avoid notation overload, we omit writing explicitly the dependence of certain quantities on the fixed parameter $k$ in the rest of this section. and constants $\alpha>0$ and $\rho\in(0,1)$ such that

[TABLE]

for every $t\geq 0$ and $f\in\mathcal{C}_{b}(E)$. Note that this implies the independence of the SCGF from the initial distribution, $\mu_{0}$, and thus (68) holds for every initial state $x\in E$. Note that (69) implies in particular that $\mu_{0}e^{-t\Lambda_{k}}P_{k}(t)$ converges weakly to $\mu_{\infty}$ for all initial distributions $\mu_{0}$, and that $\mu_{\infty}$ is the unique invariant probability measure for the modified semigroup $t\mapsto e^{-t\Lambda_{k}}P_{k}(t)$. Therefore we have from the generator $\mathcal{L}_{k}-\Lambda_{k}$ of this semigroup that

[TABLE]

Neither the semigroup $P_{k}(t)$ nor the modified one $e^{-t\Lambda_{k}}P_{k}(t)$ conserve probability, and therefore they do not provide a corresponding process to sample from and use standard MCMC methods to estimate the SCGF $\Lambda_{k}$. This can be achieved by interpreting the tilted generator $\mathcal{L}_{k}$ through Feynman-Kac models analogous to Lemma 2.1, so that we can apply our results from Section 3.

Lemma 5.2**.**

The infinitesimal generator $\mathcal{L}_{k}$ (66) can be written as

[TABLE]

for all $f\in\mathcal{C}_{b}(E)$ and $x\in E$. Here

[TABLE]

is the generator of a pure jump process with modified rates $W(x,dy)\,e^{kg(x,y)}$, and

[TABLE]

is a diagonal potential term where $\widehat{\lambda}_{k}(x):=\int_{E}W(x,dy)e^{kg(x,y)}$ is the escape rate of $\widehat{\mathcal{L}}_{k}$.

Proof.

Follows directly from the definition of $\mathcal{L}_{k}$ in (66).

∎

In analogy with (1), in the following we also use the notation with a probability kernel

[TABLE]

Observe that

[TABLE]

thus, we get with (70) another representation of the SCGF,

[TABLE]

Recall the unnormalized and normalized versions of the Feynman-Kac measures defined in (5) and (7) for a given initial distribution $\mu_{0}\in\mathcal{P}(E)$,

[TABLE]

and that asymptotic stability (69) implies that $\mu_{t}\to\mu_{\infty}$ weakly as $t\to\infty$. This suggests the following finite-time approximations for $\Lambda_{k}$.

Proposition 5.3**.**

For any $k\in\mathbb{R}$ and every $t\geq 0$, we have that

[TABLE]

where $\mathcal{V}_{k}$ is defined in (73). In particular, if asymptotic stability (69) is satisfied,

[TABLE]

Proof.

Recalling the evolution equation (6) of $\nu_{t}$, we have

[TABLE]

And, thus,

[TABLE]

since $\nu_{0}(1)=1$. We can conclude by observing that $\mathcal{L}_{k}(1)(x)=\mathcal{V}_{k}(x)$ and

[TABLE]

using that the SCGF is well defined under asymptotic stability (69).

∎

For any $t<T$, we define

[TABLE]

as a finite-time approximation for $\Lambda_{k}$.

Lemma 5.4**.**

For any $k\in\mathbb{R}$, under asymptotic stability (69) with $\rho\in(0,1)$, there exists a constant $\alpha^{\prime}>0$ such that

[TABLE]

for any given $a\in[0,1)$ and $T>0$.

Proof.

By (8), we have

[TABLE]

where $\alpha^{\prime}:=\tilde{\alpha}/(-\log\rho)>0$, using the basic fact $0\leq\rho^{aT}-\rho^{T}\leq\rho^{aT}$. In particular, $\lim_{T\to\infty}\Lambda^{aT,T}(k)=\mu_{\infty}(\mathcal{V}_{k})=\Lambda(k)$, by (76).

∎

Note that for $a=0$ the above result only implies a convergence rate of order $1/T$, since errors from the arbitrary initial condition have to be averaged out over time. In contrast for $a>0$ (corresponding to the usual idea of burn-in in conventional Markov chain Monte Carlo approximations – see [45], for example), we get a much better exponential rate of convergence dominated by the asymptotic stability parameter $\rho\in(0,1)$.

5.2 Estimation of the SCGF

In this section we establish the convergence of estimators of the SCGF, $\Lambda_{k}$ (64), provided by interacting particle approximations. Approximating $\mu_{t}$ by the empirical distribution $\mu^{N}_{t}$ (18) associated to an interacting particle system, we can estimate $\Lambda^{t,T}_{k}$ with

[TABLE]

Note that, choosing $f\equiv 1$ in Proposition 3.3 and (77) implies that \exp\big{(}{t\cdot\Lambda_{k}^{0,t,N}}\big{)} is an unbiased estimator of \exp\big{(}{t\cdot\Lambda_{k}^{0,t}}\big{)}. Recall that particle approximations are characterized by a sequence of IPS generators $(\overline{L}^{N})_{N\in\mathbb{N}}$ on $\mathcal{C}_{b}(E^{N})$, based on the McKean generators (13)

[TABLE]

where $\widetilde{\mathcal{L}}_{\mu,k}$ describes the selection dynamics of the McKean model as in Lemma 2.4, with examples in (15) or (16). Due to tilted dynamics explained in Lemma 5.2 we have an additional dependence on the parameter $k$.

Proposition 5.5**.**

Given $k\in\mathbb{R}$, let $(\overline{L}^{N}_{k})_{N\in\mathbb{N}}$ be a sequence of IPS generators satisfying Assumption 3.1 with McKean generators $\overline{\mathcal{L}}_{\mu,k}$. Under asymptotic stability (69) with $\rho\in(0,1)$, for every $p\geq 2$ and $a\in[0,1)$ there exist constants $c_{p},\,c^{\prime},\,\,\alpha^{\prime}>0$ independent of $N$ and $T$, such that

[TABLE]

and

[TABLE]

for any $N\in\mathbb{N}$ large enough and $T>0$.

Proof.

First, note that

[TABLE]

The bound for the second term is given in Lemma 5.4, whereas we can bound the first term by observing that

[TABLE]

and applying Theorem 3.2. The second claim can be established similarly.

∎

Proposition 5.5 provides the $L^{p}$ and bias estimates of the approximation error with order of convergence respectively given by $1/\sqrt{N}$ and $1/N$. The necessarily finite simulation time $T$ leads to an additional error of order ${\rho^{aT}}/{T}$, with $\rho\in(0,1)$, which is controlled by asymptotic stability properties of the process as summarized in Lemma 5.4. Ideally, during simulations we want to choose the final time $T=T(N)$ with respect to the population size $N$ in order to balance both terms in (80), resp. (81). The details depend on asymptotic stability properties of the process and values of constants, but it is clear in general that choosing any $T(N)\gg N$ would only give the same order of convergence as $T(N)\approx N$, which is computationally cheaper. Proposition 5.5 also implies that $\Lambda_{k}^{aT,T,N}$ converges almost surely to $\Lambda_{k}^{aT,T}$ as $N\to\infty$.

5.3 The Cloning Factor

Most results in the physics literature do not use the estimator $\Lambda_{k}^{aT,T,N}$ (79) based on the ergodic average of the mean fitness of the clone ensemble, but an estimator based on a so-called ‘cloning factor’ (see, e.g., [1, 5, 7]). This is essentially a continuous-time jump process $(C^{N}_{t}:t\geq 0)$ on $(0,\infty)$ with $C^{N}_{0}=1$, where at each cloning event of size $n\in\mathbb{N}_{0}\cup\{-1\}$ at a given time $\tau$, the value is updated as

[TABLE]

where $n=-1$ occurs when there is a ’killing’ event. In our context, we can define the dynamics of $C^{N}_{t}$ jointly with the cloning algorithm via an extension of the cloning generator $\overline{L}^{N}_{c,k}$ (33) as introduced in Section 3.3, with exit rate $\lambda(x)$ and probability kernel $p(x,dy)$ replaced by $\widehat{\lambda}_{k}$ and $\widehat{p}_{k}$, respectively. On the state space $E^{N}\times(0,\infty)$ define

[TABLE]

where the test function $F^{\star}:E^{N}\times(0,\infty)\to\mathbb{R}$ now has a second counting coordinate, and we denote \varsigma_{n}:=\varsigma\cdot\big{(}1{+}\tfrac{n}{N}\big{)}, with $n\in\mathbb{N}_{0}\cup\{-1\}$. Also recall that the cloning algorithm is based on a McKean model with parameter $c\in\mathbb{R}$ as given in (15).

We introduce the coordinate projection $G(\underline{x},\,\varsigma):=\varsigma$ in order to observe only the cloning factor, $G(\zeta_{t}^{N},\,C^{N}_{t})=C^{N}_{t}$. Note that $E^{N}\times(0,\infty)$ is not compact, and $G$ is an unbounded test function. However, since the range of the clone size distribution is uniformly bounded (condition 19c), $t\mapsto\log C^{N}_{t}$ is a birth-death process on $[0,\infty)$ with bounded jump length, and the generator (82) and associated semigroup is therefore well defined for the test function $G$ (see e.g. [46]) and all $t\geq 0$.

The following result provides an unbiased estimator for the unnormalized quantity $\nu_{t}(1)$ based on the cloning factor.

Proposition 5.6**.**

Let $\overline{L}_{c,k}^{(N,\star)}$ be the extension (82) of the cloning generator $\overline{L}_{c,k}^{N}$ (33). Then, the quantity $e^{tc}C^{N}_{t}$ is an unbiased estimator for $\nu_{t}(1)$ (5), i.e.

[TABLE]

for every $t\geq 0$ and $N\geq 1$, and all choices of the parameter $c\in\mathbb{R}$ (cf. (15)).

Proof.

First, observe that following (40) and (41)

[TABLE]

using the mean $M(x_{i})$ of the distribution $\pi_{x_{i},n}$ as given in Proposition 3.4. Therefore,

[TABLE]

and analogously to (24), the expected time evolution of $C^{N}_{t}$ is then given by

[TABLE]

This is also the evolution of $\nu_{t}^{N}(e^{-tc})=e^{-tc}\nu_{t}^{N}(1)$, since

[TABLE]

With initial conditions $C_{k}^{N}(0)=1=\nu_{t}^{N}(1)$, the statement follows by a Gronwall argument analogous to (23) and by Proposition 3.3.

∎

Proposition 5.6 leads to an alternative estimator for $\Lambda_{k}^{t,T}$ (78) given by

[TABLE]

Note that this is not itself unbiased as a consequence of the nonlinear transformation involving the logarithm.

In order to study the convergence of the new estimator to the SCGF, it is convenient to use the martingale characterization of the process, which is given by the following result.

Proposition 5.7**.**

Let $\overline{L}_{c,k}^{(N,\star)}$ be the extension (82) of the cloning generator $\overline{L}_{c,k}^{N}$. Then, the process

[TABLE]

with $H(\underline{x},\varsigma)=\log\varsigma$, is a local martingale satisfying

[TABLE]

and with predictable quadratic variation

[TABLE]

where $Q(x_{i})$ is the second moment of the distribution $\pi_{x_{i},n}$ (35).

Remark*.*

Note that, in the case in which there is at most one clone per transition event, i.e. if $Q(x_{i})=M(x_{i})=(\mathcal{V}_{k}(x_{i})-c)^{+}/\widehat{\lambda}_{k}(x_{i})$, then

[TABLE]

Proof.

Observe that we can rewrite (82) as

[TABLE]

using the expansion $\log(1+x)=x+O(x^{2})$ as $x\to 0$. Similarly,

[TABLE]

The statement corresponds to the martingale problem associated to $\overline{L}_{c,k}^{(N,\star)}(H)$.

∎

By Proposition 5.7 and recalling the definition of the SCGF estimators $\Lambda_{k}^{t,T,N}$ (79) and $\overline{\Lambda}_{k}^{t,T,N}$ (84) we immediately get

[TABLE]

In what follows, we discuss the convergence of the estimator $\overline{\Lambda}_{k}^{aT,T,N}$ to the SCGF $\Lambda_{k}$, which is based on the cloning factor.

Theorem 5.8**.**

Let $\overline{L}_{c,k}^{(N,\star)}$ be the extension (82) of the cloning generator $\overline{L}_{c,k}^{N}$. Then, for every $p\geq 2$ and $a\in[0,1)$, there exists a constant $C_{p}^{\star}>0$ such that for all $N$ large enough

[TABLE]

If in addition Assumption (69) on asymptotic stability holds, there exist constants $\gamma_{p}^{\star},\,c_{p}^{\star},\alpha^{\prime}>0$ and $0<\rho<1$ (dependent on $a,\,p,\,\widehat{\lambda}_{k},\,{Q}$ and $\mathcal{V}_{k}$) such that

[TABLE]

for every $T\geq 1$.

Proof.

Thanks to Jensen’s inequality, it is enough to prove the inequality for all $p=2^{q}$, $q\in\mathbb{N}$. First, we can write

[TABLE]

Observe that \sup_{t\leq T}\big{|}\mathcal{M}^{\star}_{t}\big{|}<\infty, so the assumptions of Lemma 4.1 are satisfied. Thus, using Lemma 4.1, we obtain

[TABLE]

The second part of the Theorem follows directly by Proposition 5.5.

∎

Therefore, the $L^{p}$-error for estimator $\overline{\Lambda}_{k}^{aT,T,N}$ has the same rate of convergence $1/\sqrt{N}$ as $\Lambda^{aT,T,N}_{k}$. Analogous results hold for the bias estimates, which have order of convergence $1/N$ as for the estimator $\Lambda_{k}^{aT,T,N}$ (Proposition 5.5), since with (85) the difference of both estimators is only of order $N^{-1-1/p}$.

6 Discussion

In this work we have established a framework to compare variants of cloning algorithms and understand their connections with mean field particle approximations. This allowed us to obtain first rigorous results on the convergence properties of cloning algorithms in continuous time. Our results apply in the general setting of jump Markov processes on locally compact state spaces. Essential conditions for our approach are summarized in Assumptions 2.2 on asymptotic stability of the process and 3.1 on the particle approximation, which are usually straightforward to check for practical applications. We summarize further sufficient conditions for asymptotic stability in the Appendix A.

In certain situations the cloning algorithm is computationally cheaper and simpler to implement than mean field particle systems, since only the mutation process has to be sampled independently for all particles and cloning events happen simultaneously. However, as discussed in [34], this choice reduces in general the accuracy of the estimator since it does not consider the fitness potential of the replaced particles during the cloning events. Adjusting the algorithm by allowing only substitutions of particles with lower fitness based on different McKean models could improve the accuracy. The approach developed in this paper can be used to conduct a systematic study of this question, which is current work in progress.

Appendix A Asymptotic Stability

We present sufficient conditions for asymptotic stability as presented in Assumption 2.2. The discussion is based on the work of Tweedie et al. [47, 48], which we briefly recall in Lemma A.3 below.

Definition A.1**.**

A Feller process $Y_{t}$ is said to be $\phi$-irreducible for a non-trivial measure $\phi$ (i.e. $\phi(E)>0$) on $(E,\mathcal{B}(E))$, if \mathbb{E}_{x}\big{[}\int_{0}^{\infty}\mathbbm{1}_{Y_{t}\in A}dt\big{]}>0 for every $x\in E$ and every set $A\in\mathcal{B}(E)$ such that $\phi(A)>0$. We simply say that $Y_{t}$ is irreducible if it is $\phi$-irreducible for some $\phi$.

Definition A.2**.**

A $\phi$-irreducible Feller process $Y_{t}$ is called aperiodic if there exists a small set $C\in\mathcal{B}(E)$, $\phi(C)>0$, such that the associated Markov semigroup $P(t)$ satisfies the following conditions:

• •

there exists a non-trivial measure $\eta$ and $t>0$ such that $P(t)\,(x,B)\geq\eta(B)$, for all $x\in C$ and $B\in\mathcal{B}(E)$;

• •

there exists a time $\tau\geq 0$ such that $P(t)\,(x,C)>0$, for all $t\geq\tau$ and $x\in C$.

Lemma A.3**.**

Let $Y_{t}$ be a $\phi$-irreducible and aperiodic Feller process on a locally compact state space $E$ such that $\mathrm{supp}\,\phi$ has non-empty interior. Denote by $\mathcal{L}$ and $P(t)$ the associated infinitesimal generator and the semigroup, respectively. Assume that for a given function $h\in\mathcal{C}_{b}(E)$ such that $h\geq 1$, there exist constants $b,c>0$ and a compact set $S\in\mathcal{B}(E)$ such that for all $x\in E$

[TABLE]

Then there exist constants $\alpha\geq 0$ and $\rho\in(0,1)$ such that for any test function $f\in\mathcal{C}_{b}(E)$ and $t\geq 0$,

[TABLE]

for any $x\in E$, where $\pi$ is the (unique) invariant measure of $Y_{t}$.

Proof.

See [47], Theorem 5.2(c), using the fact that if a Feller process $Y_{t}$ is $\phi$-irreducible and $\mathrm{supp}\,\phi$ has non-empty interior, then every compact set is petite (See [48], Theorem 7.1 and Theorem 5.1).

∎

In the following we discuss how the spectral properties of the tilted generator $\mathcal{L}^{\mathcal{V}}$ in Assumption 2.2 can imply asymptotic stability in the sense of (4).

Assumption A.4**.**

We assume that the spectrum of $\mathcal{L}^{\mathcal{V}}$ (2) is bounded by a greatest eigenvalue $\lambda_{0}$. Moreover, there exist a positive function $r\in\mathcal{C}_{b}(E)$, unique up to multiplicative constants, and a probability measure $\mu_{\infty}\in\mathcal{P}(E)$ satisfying respectively

[TABLE]

and

[TABLE]

Without loss of generality, we can assume $\mu_{\infty}(r)=1$.

Remark*.*

Sufficient conditions for Assumption A.4 to hold can be found, for instance, in [49, 50]. These are of course satisfied if the original process with generator $\mathcal{L}$ is an irreducible, finite-state Markov chain, including for example stochastic particle systems on finite lattices with a fixed number of particles.

Under Assumption A.4, we define the generator

[TABLE]

which is known in the literature as Doob’s $h$-transform of $\mathcal{L}^{\mathcal{V}}$ [42] or twisted Markov kernel [51]. Observe that $\mathcal{L}^{\mathcal{V}}_{r}1=0$, so that it is a probability generator associated to a Markov process with probability semigroup defined for any $f\in\mathcal{C}_{b}(E)$ by

[TABLE]

Proposition A.5** (Asymptotic stability).**

Assume that there exists $\varepsilon>0$ such that the set

[TABLE]

is compact. Under Assumption A.4, if the initial pure jump process $(X_{t}:t\geq 0)$ with generator $\mathcal{L}$ is $\phi$-irreducible for some $\phi$ for which $\textrm{supp}\,\phi$ has non-empty interior, and aperiodic as defined above then (4) holds, i.e. there exists $\alpha>0$ and $\rho\in(0,1)$ such that

[TABLE]

for every $t\geq 0$ and $f\in\mathcal{C}_{b}(E)$.

Proof.

First, note that if the initial process $X_{t}$ is irreducible and aperiodic, then also the process associated to $\mathcal{L}^{\mathcal{V}}_{r}$ is irreducible and aperiodic. Moreover, $\mathcal{L}^{\mathcal{V}}_{r}$ is bounded in ${K_{\varepsilon}}$ and $\mathcal{L}^{\mathcal{V}}_{r}(r^{-1})\leq-\varepsilon\,r^{-1}$ for every $x\not\in K_{\varepsilon}$. Therefore, the hypotheses of Lemma A.3 are satisfied for the generator $\mathcal{L}^{\mathcal{V}}_{r}$ acting on the function $h=r^{-1}$. Thus, applying the lemma we obtain

[TABLE]

for any $f\in\mathcal{C}_{b}(E)$ and $x\in E$, where $\pi(\cdot)=\mu_{\infty}(r\,\cdot)\in\mathcal{P}(E)$ is the invariant measure for $\mathcal{L}^{\mathcal{V}}_{r}$. Dividing by $r^{-1}(x)$ and substituting $f$ with $r^{-1}f\in\mathcal{C}_{b}(E)$, we obtain the statement ($\|r^{-1}\|<\infty$ and can be included in the constant $\alpha$).

∎

Acknowledgments

This work was supported by The Alan Turing Institute under the EPSRC grant EP/N510129/1 and the Lloyd’s Register Foundation–Alan Turing Institute Programme on Data-Centric Engineering; AMJ was partially supported by EPSRC grants EP/R034710/1 and EP/T004134/1.