On backward smoothing algorithms

A state-space model is composed of an unobserved Markov process $X_{0},\ldots,X_{T}$ and observed data $Y_{0},\ldots,Y_{T}$. Given $X_{0},\ldots,X_{T}$, the data $Y_{0},\ldots,Y_{T}$ are independent and generated through some specified emission distribution $Y_{t}|X_{t}\sim\boldsymbol{f}_{t}(\mathrm{d}y_{t}|x_{t})$. These models have wide-ranging applications (e.g. in biology, economics and engineering). Two important inference tasks related to state-space models are filtering (computing the distribution of $X_{t}$ given $Y_{0},\ldots,Y_{t}$) and smoothing (computing the distribution of the whole trajectory $(X_{0},\ldots,X_{t})$, again given all data until time $t$). Filtering is usually carried out through a particle filter, that is, a sequential Monte Carlo algorithm that propagates $N$ weighted particles (realisations) through Markov and importance sampling steps; see Chopin and Papaspiliopoulos, (2020) for a general introduction to state-space models (Chapter 2) and particle filters (Chapter 10).

This paper is concerned with skeleton-based smoothing algorithms, i.e. algorithms that approximate the smoothing distributions with empirical distributions based on the output of a particle filter (i.e. the locations and weights of the $N$ particles at each time step). A simple example is genealogy tracking (empirically proposed in Kitagawa,, 1996 and theoretically analysed in Del Moral and Miclo,, 2001) which keeps track of the ancestry (past states) of each particles. This smoother suffers from degeneracy: for $t$ large enough, all the particles have the same ancestor at time $0$.

The forward filtering backward smoothing (FFBS) algorithm (Godsill et al.,, 2004) has been proposed as a solution to this problem. Starting from the filtering approximation at time $t$, the algorithm samples successively particles at times $t-1$, $t-2$, etc. using backward kernels. Its theoretical properties, in particular the stability as $t\to\infty$, have been studied by Del Moral et al., (2010); Douc et al., (2011).

In many applications, one is mainly interested in approximating smoothing expectations of additive functions of the form

for some functions $\psi_{0},\ldots,\psi_{t}$. Such expectations can be approximated in an online fashion by a procedure described in Del Moral et al., (2010). Inspired by this, the particle-based, rapid incremental smoother (PaRIS) algorithm of Olsson and Westerborn, (2017) replaces some of the calculations with an additional layer of Monte Carlo approximation.

The backward sampling operation is central to both the FFBS and the PaRIS algorithms. The naive implementation has an $\mathcal{O}(N^{2})$ cost. There have been numerous attempts at alleviating this problem in the literature, but, to our knowledge, they all lack formal support in terms of either computation complexity or stability as $T\to\infty$.

In the following five paragraphs, we elaborate on this limitation for each of the three major contenders, and we point out two related challenges with current backward sampling algorithms that we try to resolve in this article.

First, Douc et al., (2011) proposed to use rejection sampling for the generation of backward indices in FFBS, and Olsson and Westerborn, (2017) extended this technique to PaRIS. If the model has upper-and-lower bounded transition densities, this sampler has an $\mathcal{O}(N)$ expected execution time (Douc et al., 2011, Proposition 2 and Olsson and Westerborn, 2017, Theorem 10). Unfortunately, most practical state space models (including linear Gaussian ones) violate this assumption, and the behaviour of the algorithm in this case is unclear. Empirically, it has been observed (Taghavi et al., 2013; Bunch and Godsill, 2013; Olsson and Westerborn, 2017, Section 4.3) that in real examples, FFBS-reject and PaRIS-reject frequently suffer from low acceptance rates, in contrary to what users would expect from an algorithm with linear complexity. To cite Bunch and Godsill, (2013), “[a]lthough theoretically elegant, the […] algorithm has been found to suffer from such high rejection rates as to render it consistently slower than direct sampling implementation on problems with more than one state dimension”. To the best of our knowledge, no theoretical result has been put forward to formalise or quantify this bad behaviour.

Second, Taghavi et al., (2013) and Olsson and Westerborn, (2017, Section 4.3) suggest putting a threshold on the number of rejection sampling trials to get more stable execution times. The thresholds are either chosen adaptively using a Kalman filter in Taghavi et al., (2013) or fixed at $\sqrt{N}$ in Olsson and Westerborn, (2017, Section 4.3). Although improvements are empirically observed, to the best of our knowledge, no theoretical analysis of the complexity of the resulting algorithm or formal justification of the proposed threshold is available.

Third, Bunch and Godsill, (2013) use MCMC moves starting from the filtering ancestor instead of a full backward sampling step. Empirically, this algorithm seems to prevent the degeneracy associated with the genealogy tracking smoother using a very small number of MCMC steps (e.g. less than five). Unfortunately, as far as we know, this stability property is not proved anywhere in the literature, which deters the adoption of the algorithm. Using MCMC moves provides a procedure with truly linear and deterministic run time, and a stability result is the only missing piece of the puzzle to completely resolve the $\mathcal{O}(N^{2})$ problem. We believe one reason for the current state of affair is that the stability proof techniques employed by Douc et al., (2011) and Olsson and Westerborn, (2017) are difficult to extend to the MCMC case.

Fourth, and this is related to the third point, the stability of the PaRIS algorithm has only been proved in the asymptotic regime. More specifically, Olsson and Westerborn, (2017) established a central limit theorem as $N\to\infty$ in Corollary 5, then showed that the corresponding asymptotic error remains controlled as $T\to\infty$ in Theorem 8 and Theorem 9. While non-asymptotic stability bounds for the FFBS algorithm are already available in Del Moral et al., (2010); Douc et al., (2011); Dubarry and Le Corff, (2013), we do not think that they can be extended straightforwardly to PaRIS and we are not aware of any such attempt in the literature.

Fifth, all backward samplers mentioned thus far require access to the transition density. Many models have dynamics that can be simulated from but transition densities that are not explicitly calculable. Enabling backward sampling in this scenario is challenging and will certainly require some kind of problem-specific knowledge to extract information from the transition densities, despite not being able to evaluate them exactly.

Section 2 presents a general framework which admits as particular cases a wide variety of existing algorithms (e.g. FFBS, forward-additive smoothing, PaRIS) as well as the novel ones considered later in the paper. It allows to simultaneously prove the consistency as $N\to\infty$ and the stability as $T\to\infty$ for all of them. The main ingredient is the discrete backward kernels, which are essentially random $N\times N$ matrices employed differently in the offline and the online settings. On the technical side, the stability result is proved using a new technique, yielding a non-asymptotic bound that addresses the fourth point in subsection 1.2.

Next, we closely look at the use of rejection sampling and realise that in many models, the resulting execution time may be significantly heavy-tailed; see Section 3. For instance, the run time of PaRIS may have infinite expectation, whereas the run time of FFBS may have infinite variance. (Since it is technically more involved, the material for the FFBS algorithm is delegated to Supplement B.) These results address the first point in subsection 1.2 and we discuss their severe practical implications.

We then derive and analyse hybrid rejection sampling schemes (i.e. schemes that use rejection sampling only up to a certain number of attempts, and then switch to the standard method). We show that they lead to a nearly $\mathcal{O}(N)$ algorithm (up to some $\log$ factor) in Gaussian models; again see Section 3. This stems from the subtle interaction between the tail of Gaussian densities and the finite Feynman-Kac approximation. Outside this class of model, the hybrid algorithm can still escape the $\mathcal{O}(N^{2})$ complexity, although it might not reach the ideal linear run time target. These results shed some light on the second issue mentioned in subsection 1.2.

In Section 4, we look at backward kernels that are more efficient to simulate than the FFBS and the PaRIS ones. Section 4.1 describes backward kernels based on MCMC (Markov chain Monte Carlo) following Bunch and Godsill, (2013) and extends them to the online scenario. We cast this family of algorithms as a particular case of the general framework developed in Section 2, which allows convergence and stability to follow immediately. This solves the long-standing problem described in the third point of subsection 1.2.

MCMC methods require evaluation of the likelihood and thus cannot be applied to models with intractable transition densities. In Section 4.2, we show how the use of forward coupling can replace the role of backward MCMC steps in these scenarios. This makes it possible to obtain stable performance in both on-line and off-line scenarios (with intractable transition densities) and provides a possible solution to the fifth challenge describe in subsection 1.2.

Section 5 illustrates the aforementioned algorithms in both on-line and off-line uses. We highlight how hybrid and MCMC samplers lead to a more user-friendly (i.e. smaller, less random and less model-dependent) execution time than the pure rejection sampler. We also apply our smoother for intractable densities to a continuous-time diffusion process with discretization. We observe that our procedure can indeed prevent degeneracy as $T\to\infty$, provided that some care is taken to build couplings with good performance. Section 6 concludes the paper with final practical recommendations and further research directions. Table 1 gives an overview of existing and novel algorithms as well as our contributions for each.

Proposition 1 of Douc et al., (2011) states that under certain assumptions, the FFBS-reject algorithm has an asymptotic $\mathcal{O}_{\mathbb{P}}(N)$ complexity. This does not contradict our results, which point out the undesirable properties of the non-asymptotic execution time. Clearly, non-asymptotic behaviours are what users really observe in practice. From a technical point of view, the proof of Douc et al., (2011, Prop. 1) is a simple application of Theorem 5 of the same paper. In contrast, non-asymptotic results such as Theorem B.1 and Theorem B.2 require more delicate finite sample analyses.

Figure 1 of Olsson and Westerborn, (2017) and the accompanying discussion provide an excellent intuition on the stability of smoothing algorithms based on the support size of the backward kernels. We formalise this support size condition for the first time by the inequality  (13) and construct a novel non-asymptotic stability result based on it. In contrast, Olsson and Westerborn, (2017) depart from their initial intuition and use an entirely different technique to establish stability. Their result is asymptotic in nature.

Gloaguen et al., (2022) briefly mention the use of MCMC in PaRIS algorithm, but their algorithm is fundamentally different to and less efficient than Bunch and Godsill, (2013). Indeed, they do not start the MCMC chains at the ancestors previously obtained during the filtering step. They are thus obliged to perform a large number of MCMC iterations for decorrelation, whereas the algorithms described in our Proposition 4, built upon the ideas of Bunch and Godsill, (2013), only require a single MCMC step to guarantee stability. However, we would like to stress again that Bunch and Godsill, (2013) did not prove this important fact.

Another way to reduce the computation time is to perform the expensive backward sampling steps at certain times $t$ only. For other values of $t$, the naive genealogy tracking smoother is used instead. This idea has been recently proposed by Mastrototaro et al., (2021), who also provided a practical recipe for deciding at which values of $t$ the backward sampling should take place and derived corresponding theoretical results.

Smoothing in models with intractable transition densities is very challenging. If these densities can be estimated accurately, the algorithms proposed by Gloaguen et al., (2022) permit to attack this problem. A case of particular interest is diffusion models, where unbiased transition density estimators are provided in Beskos et al., (2006); Fearnhead et al., (2008). More recently, Yonekura and Beskos, (2022) use a special bridge path-space construction to overcome the unavailability of transition densities when the diffusion (possibly with jumps) must be discretised.

Our smoother for intractable models are based on a general coupling principle that is not specific to diffusions. We only require users to be able to simulate their dynamics (e.g. using discretisation in the case of diffusions) and to manipulate random numbers in their simulations so that dynamics starting from two different points can meet with some probability. Our method does not directly provide an estimator for the gradient of the transition density with respect to model parameters and thus cannot be used in its current form to perform maximum likelihood estimation (MLE) in intractable models; whereas the aforementioned work have been able to do so in the case of diffusions. However, the main advantage of our approach lies in its generality beyond the diffusion case. Furthermore, modifications allowing to perform MLE are possible and might be explored in further work specifically dedicated to the parameter estimation problem.

The idea of coupling has been incorporated in the smoothing problem in a different manner by Jacob et al., (2019). There, the goal is to provide offline unbiased estimates of the expectation under the smoothing distribution. Coupling and more generally ideas based on correlated random numbers are also useful in the context of partially observed diffusions via the multilevel approach (Jasra et al.,, 2017).

In this work, we consider smoothing algorithms that are based on a unique pass of the particle filter. Offline smoothing can be done using repeated iterations of the conditional particle filter (Andrieu et al.,, 2010). Full trajectories can also be constructed in an online manner if one is willing to accept some lag approximations (Duffield and Singh,, 2022). Another approach to smoothing consists of using an additional information filter (Fearnhead et al.,, 2010), but it is limited to functions depending on one state only. Each of these algorithmic families has their own advantages and disadvantages, of which a detailed discussion is out of the scope of this article (see however Nordh and Antonsson,, 2015).

Measure-kernel-function notations. Let $\mathcal{X}$ and $\mathcal{Y}$ be two measurable spaces with respective $\sigma$-algebras $\mathcal{B}(\mathcal{X})$ and $\mathcal{B}(\mathcal{Y})$. The following definitions involve integrals and only make sense when they are well-defined. For a measure $\mu$ on $\mathcal{X}$ and a function $f:\mathcal{X}\to\mathbb{R}$, the notations $\mu f$ and $\mu(f)$ refer to $\int f(x)\mu(\mathrm{d}x)$. A kernel (resp. Markov kernel) $K$ is a mapping from $\mathcal{X}\times\mathcal{B}(\mathcal{Y})$ to $\mathbb{R}$ (resp. $[0,1]$) such that, for $B\in\mathcal{B}(\mathcal{Y})$ fixed, $x\mapsto K(x,B)$ is a measurable function on $\mathcal{X}$; and for $x$ fixed, $B\mapsto K(x,B)$ is a measure (resp. probability measure) on $\mathcal{Y}$. For a real-valued function $g$ defined on $\mathcal{Y}$, let $Kg:\mathcal{X}\to\mathbb{R}$ be the function $Kg(x):=\int g(y)K(x,\mathrm{d}y)$. We sometimes write $K(x,g)$ for the same expression. The product of the measure $\mu$ on $\mathcal{X}$ and the kernel $K$ is a measure on $\mathcal{Y}$, defined by $\mu K(B):=\int K(x,B)\mu(\mathrm{d}x)$. Other notations. • The notation $X_{0:t}$ is a shorthand for $(X_{0},\ldots,X_{t})$ • We denote by $\mathcal{M}(W^{1:N})$ the multinomial distribution supported on $\left\{1,2,\ldots,N\right\}$. The respective probabilities are $W_{1},\ldots,W_{N}$. If they do not sum to $1$, we implicitly refer to the normalised version obtained by multiplication of the weights with the appropriate constant • The symbol $\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}$ means convergence in probability and $\Rightarrow$ means convergence in distribution • The geometric distribution with parameter $\lambda$ is supported on $\mathbb{Z}_{\geq 1}$, has probability mass function $f(n)=\lambda(1-\lambda)^{n-1}$ and is noted by $\operatorname{Geo}(\lambda)$ • Let $\mathcal{X}$ and $\mathcal{Y}$ be two measurable spaces. Let $\mu$ and $\nu$ be two probability measures on $\mathcal{X}$ and $\mathcal{Y}$ respectively. The o-times product measure $\mu\otimes\nu$ is defined via $(\mu\otimes\nu)(h):=\iint h(x,y)\mu(\mathrm{d}x)\nu(\mathrm{d}y)$ for bounded functions $h:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}$. If $X\sim\mu$ and $Y\sim\nu$, we sometimes note $\mu\otimes\nu$ by $X\otimes Y$.

Let $\mathcal{X}_{0:T}$ be a sequence of measurable spaces and $M_{1:T}$ be a sequence of Markov kernels such that $M_{t}$ is a kernel from $\mathcal{X}_{t-1}$ to $\mathcal{X}_{t}$. Let $X_{0:T}$ be an unobserved inhomogeneous Markov chain with starting distribution $X_{0}\sim\mathbb{M}_{0}(\mathrm{d}x_{0})$ and Markov kernels $M_{1:T}$; i.e. $X_{t}|X_{t-1}\sim M_{t}(X_{t-1},\mathrm{d}x_{t})$ for $t\geq 1$. We aim to study the distribution of $X_{0:T}$ given observed data $Y_{0:T}$. Conditioned on $X_{0:T}$, the data $Y_{0},\ldots,Y_{T}$ are independent and $Y_{t}|X_{0:T}\equiv Y_{t}|X_{t}\sim\boldsymbol{f}_{t}(\cdot|X_{t})$ for a certain emission distribution $\boldsymbol{f}_{t}(\mathrm{d}y_{t}|x_{t})$. Assume that there exists dominating measures $\tilde{\lambda}_{t}$ not depending on $x_{t}$ such that

The distribution of $X_{0:t}|Y_{0:t}$ is then given by

where $G_{s}(x_{s}):=f(y_{s}|x_{s})$ and $L_{t}>0$ is the normalising constant. Moreover, $\mathbb{Q}_{-1}:=\mathbb{M}_{0}$ and $L_{-1}:=1$ by convention. Equation (1) defines a Feynman-Kac model (Del Moral,, 2004). It does not require $M_{t}$ to admit a transition density, although herein we only consider models where this assumption holds. Let $\lambda_{t}$ be a dominating measure on $\mathcal{X}_{t}$ in the sense that there exists a function $m_{t}$ (not necessarily tractable) such that

A special case of the current framework are linear Gaussian state space models. They will serve as a running example for the article, and some of the results will be specifically demonstrated for models of this class. The rationale is that many real-world dynamics are partly, or close to, Gaussian. The notations for linear Gaussian models are given in Supplement A.1 and we will refer to them whenever this model class is discussed.

Particle filters are algorithms that sample from $\mathbb{Q}_{t}(\mathrm{d}x_{t})$ in an on-line manner. In this article, we only consider the bootstrap particle filter (Gordon et al.,, 1993) and we detail its notations in Algorithm 1. Many results in the following do apply to the auxiliary filter (Pitt and Shephard,, 1999) as well, and we shall as a rule indicate explicitly when it is not the case.

We end this subsection with the definition of two sigma-algebras that will be referred to throughout the paper. Using the notations of Algorithm 1, let

In this subsection, we first describe three examples of backward kernels, in which we emphasise both the random measure and the random matrix viewpoints. We then formalise their use by stating a generic off-line smoothing algorithm.

We formalise how off-line smoothing can be done given random matrices $\hat{B}_{1:T}^{N}$; see Algorithm 2. Note that in the above examples, our matrices $\hat{B}_{t}^{N}$ are $\mathcal{F}_{t}$-measurable (i.e. they depend on particles and indices up to time $t$), but this is not necessarily the case in general (i.e. they may also depend on additional random variables, see Section 2.5). Furthermore, Algorithm 2 describes how to perform smoothing using the matrices $\hat{B}_{1:T}^{N}$, but does not say where they come from. At this point, it is useful to keep in mind the above three examples. In Section 2.4, we will give a general recipe for constructing valid matrices $\hat{B}_{t}^{N}$ (i.e. those that give a consistent algorithm).

Algorithm 2 simulates, given $\mathcal{F}_{T}$ and $\hat{B}_{1:T}^{N}$, $N$ i.i.d. index sequences $\mathcal{I}_{0:T}^{n}$, each distributed according to

Once the indices $\mathcal{I}_{0:T}^{1:N}$ are simulated, the $N$ smoothed trajectories are returned as $(X_{0}^{\mathcal{I}_{0}^{n}},\ldots,X_{T}^{\mathcal{I}_{T}^{n}})$. Given $\mathcal{F}_{T}$ and $\hat{B}_{1:T}^{N}$, they are thus conditionally i.i.d. and their conditional distribution is described by the $x_{0:T}$ component of the joint distribution

Throughout the paper, the symbol $\bar{\mathbb{Q}}_{T}^{N}$ will refer to this joint distribution, while the symbol $\mathbb{Q}_{T}^{N}$ will refer to the $x_{0:T}$-marginal of $\bar{\mathbb{Q}}_{T}^{N}$ only. This allows the notation $\mathbb{Q}_{T}^{N}\varphi$ to make sense, where $\varphi=\varphi(x_{0},\ldots,x_{T})$ is a real-valued function defined on the hidden states.

The kernels $B_{t}^{N,\text{FFBS}}$ and $B_{t}^{N,\text{GT}}$ are both valid backward kernels to generate convergent approximation of the smoothing distribution (Del Moral,, 2004; Douc et al.,, 2011). This subsection shows that they are not the only ones and gives a sufficient condition for a backward kernel to be valid. It will prove a necessary tool to build more efficient $B_{t}^{N}$ later in the paper.

Recall that Algorithm 1 outputs particles $X_{0:T}^{1:N}$, weights $W_{0:T}^{1:N}$ and ancestor variables $A_{1:T}^{1:N}$. Imagine that the $A_{1:T}^{1:N}$ were discarded after filtering has been done and we wish to simulate them back. We note that, since the $X_{0:T}^{1:N}$ are given, the $T\times N$ variables $A_{1:T}^{1:N}$ are conditionally i.i.d. We can thus simulate them back from

This is precisely the distribution of $B_{t}^{N,\text{FFBS}}(n,\cdot)$. It turns out that any other invariant kernel that can be used for simulating back the discarded $A_{1:T}^{1:N}$ will lead to a convergent algorithm as well. For instance, $B_{t}^{N,\mathrm{GT}}(n,\cdot)$ (Example 2) simply returns back the old $A_{t}^{n}$, unlike $B_{t}^{N,\mathrm{FFBS}}(n,\cdot)$ which creates a new version. The kernel $B_{t}^{N,\mathrm{IMH}}(n,\cdot)$ (Example 3) is somewhat an intermediate between the two. We formalise these intuitions in the following theorem. It is stated for the bootstrap particle filter, but as a matter of fact, the proof can be extended straightforwardly to auxiliary particle filters as well.

A typical relation between variables defined in the statement of the theorem is illustrated by a graphical model in Figure 1. (See Bishop, 2006, Chapter 8 for the formal definition of graphical models and how to use them.) By “typical”, we mean that Theorem 2.1 technically allows for more complicated relations, but the aforementioned figure captures the most essential cases.

Theorem 2.1 is a generalisation of Douc et al., (2011, Theorem 5). Its proof thus follows the same lines (Supplement E.1). However, in our case the measure $\mathbb{Q}_{T}^{N}(\mathrm{d}x_{0:T})$ is no longer Markovian. This is because the backward kernel $B_{t}^{N}(i_{t},\mathrm{d}i_{t-1})$ does not depend on $X_{t}^{i_{t}}$ alone, but also possibly on its ancestor and extra random variables. This small difference has a big consequence: compared to Douc et al., (2011, Theorem 5), Theorem 2.1 has a much broader applicability and encompasses, for instance, the MCMC-based algorithms presented in Section 4.1 and novel kernels presented in Section 4.2 for intractable densities.

As we have seen in (7), $\mathbb{Q}_{T}^{N}$ is fundamentally a discrete measure of which the support contains $N^{T+1}$ elements. As such, $\mathbb{Q}_{T}^{N}\varphi$ cannot be computed exactly in general and must be approximated using $N$ trajectories $(X_{0}^{\mathcal{I}_{0}^{n}},\ldots,X_{T}^{\mathcal{I}_{T}^{n}})$ simulated via Algorithm 2. Theorem 2.1 is thus completed by the following corollary, which is an immediate consequence of Hoeffding inequality (Supplement E.13).

As we have seen in Section 2.3 and Section 2.4, in general, the expectation $\mathbb{Q}_{T}^{N}\varphi$, for a real-valued function $\varphi=\varphi(x_{0},\ldots,x_{T})$ of the hidden states, cannot be computed exactly due to the large support ($N^{T+1}$ elements) of $\mathbb{Q}_{T}^{N}$. Moreover, in certain settings we are interested in the quantities $\mathbb{Q}_{t}^{N}\varphi_{t}$ for different functions $\varphi_{t}$. They cannot be approximated in an on-line manner without more assumptions on the connection between $\varphi_{t-1}$ and $\varphi_{t}$. If the family $(\varphi_{t})$ is additive, i.e. there exists functions $\psi_{t}$ such that

then we can calculate $\mathbb{Q}_{t}^{N}\varphi_{t}$ both exactly and on-line. The procedure was first described in Del Moral et al., (2010) for the kernel $\mathbb{Q}_{t}^{N,\text{FFBS}}$ (i.e. the measure defined by (7) and the random kernels $B_{t}^{N,\text{FFBS}}$), but we will use the idea for other kernels as well. In this subsection, we first explain the principle of the method, then discuss its computational complexity and the link to the PaRIS algorithm (Olsson and Westerborn,, 2017).

When $\hat{B}_{t}^{N}\equiv\hat{B}_{t}^{N,\text{GT}}$, Algorithms 2 and 3 reduce to the genealogy tracking smoother (Kitagawa,, 1996). The matrix $\hat{B}_{t}^{N,\text{GT}}$ is indeed sparse, leading to the well-known $\mathcal{O}(N)$ complexity of this on-line procedure. As per Theorem 2.1, smoothing via genealogy tracking is convergent at rate $\mathcal{O}(N^{-1/2})$ if $T$ is fixed. When $T\to\infty$ however, all particles will eventually share the same ancestor at time $0$ (or any fixed time $t$). Mathematically, this phenomenon is manifested in two ways: (a) for fixed $t$ and function $\phi_{t}:\mathcal{X}_{t}\to\mathbb{R}$, the error of estimating $\mathbb{E}[\phi_{t}(X_{t})|Y_{0:T}]$ grows linearly with $T$; and (b) the error of estimating $\mathbb{E}\left[\left.{\sum_{t=0}^{T}\psi_{t}(x_{t-1},x_{t})}\right|{Y_{0:T}}\right]$ grows quadratically with $T$. These correspond respectively to the degeneracy for the fixed marginal smoothing and the additive smoothing problems; see also the introductory section of Olsson and Westerborn, (2017) for a discussion. The random matrices $\hat{B}_{t}^{N,\text{GT}}$ are therefore said to be unstable as $T\to\infty$, which is not the case for $\hat{B}_{t}^{N,\text{FFBS}}$ or $\hat{B}_{t}^{N,\text{PaRIS}}$. This subsection gives sufficient conditions to ensure the stability of a general $\hat{B}_{t}^{N}$.

The essential point behind smoothing stability is simple: the support of $B_{t}^{N,\mathrm{FFBS}}(n,\cdot)$ or $B_{t}^{N,\mathrm{PaRIS}}(n,\cdot)$ for $\tilde{N}\geq 2$ contains more than one element, contrary to that of $B_{t}^{N,\mathrm{GT}}(n,\cdot)$. This property is formalised by (13). To explain the intuitions, we use the notations of Algorithm 2 and consider the estimate

of $\mathbb{E}\left[\left.{\psi_{0}(X_{0})}\right|{Y_{0:T}}\right]$ when $T\to\infty$. The variance of the quantity above is a sum of $\mathrm{Cov}(\psi_{0}(X_{0}^{\mathcal{I}_{0}^{i}}),\psi_{0}(X_{0}^{\mathcal{I}_{0}^{j}}))$ terms. It can therefore be understood by looking at a pair of trajectories simulated using Algorithm 2.

At final time $t=T$, $\mathcal{I}_{T}^{1}$ and $\mathcal{I}_{T}^{2}$ both follow the $\mathcal{M}(W_{T}^{1:N})$ distribution. Under regularity conditions (e.g. no extreme weights), they are likely to be different, i.e., $\mathbb{P}(\mathcal{I}_{T}^{1}=\mathcal{I}_{T}^{2})=\mathcal{O}(1/N)$. This property can be propagated backward: as long as $\mathcal{I}_{t}^{1}\neq\mathcal{I}_{t}^{2}$, the two variables $\mathcal{I}_{t-1}^{1}$ and $\mathcal{I}_{t-1}^{2}$ are also likely to be different, with however a small $\mathcal{O}(1/N)$ chance of being equal. Moreover, as long as the two trajectories have not met, they can be simulated independently given $\mathcal{F}_{T}^{-}$ (the sigma algebra defined in (3)). In mathematical terms, under the two hypotheses of Theorem 2.1, given $\mathcal{F}_{T}^{-}$ and $\mathcal{I}_{t:T}^{1,2}$, it can be proved that the two variables $\mathcal{I}_{t-1}^{1}$ and $\mathcal{I}_{t-1}^{2}$ are independent if $\mathcal{I}_{t}^{1}\neq\mathcal{I}_{t}^{2}$ (Lemma E.1, Supplement E.3).

Since there is an $\mathcal{O}(1/N)$ chance of meeting at each time step, if $T\gg N$, it is likely that the two paths will meet at some point $t\gg 0$. When $\mathcal{I}_{t}^{1}=\mathcal{I}_{t}^{2}$, the two indices $\mathcal{I}_{t-1}$ and $\mathcal{I}_{t-2}$ are both simulated according to $B_{t}^{N}(\mathcal{I}_{t}^{1},\cdot)$. In the genealogy tracking algorithm, $B_{t}^{N,\text{GT}}(i,\cdot)$ is a Dirac measure, leading to $\mathcal{I}_{t-1}^{1}=\mathcal{I}_{t-1}^{2}$ almost surely. This spreads until time $0$, so $\operatorname{Corr}(\psi_{0}(X_{0}^{\mathcal{I}_{0}^{1}}),\psi_{0}(X_{0}^{\mathcal{I}_{0}^{2}}))$ is almost $1$ if $T\gg N$.

Other kernels like $B_{t}^{N,\text{FFBS}}$ or $B_{t}^{N,\text{PaRIS}}$ do not suffer from the same problem. For these, the support size of $B_{t}^{N}(\mathcal{I}_{t}^{1},\cdot)$ is greater than one and thus there is some real chance that $\mathcal{I}_{t-1}^{1}\neq\mathcal{I}_{t-1}^{2}$. If that does happen, we are again back to the regime where the next states of the two paths can be simulated independently. Note also that the support of $B_{t}^{N}(\mathcal{I}_{t}^{1},\cdot)$ does not need to be large and can contain as few as $2$ elements. Even if $\mathcal{I}_{t-1}^{1}$ might still be equal to $\mathcal{I}_{t-1}^{2}$ with some probability, the two paths will have new chances to diverge at times $t-2$, $t-3$ and so on. Overall, this makes $\operatorname{Corr}(\psi_{0}(X_{0}^{\mathcal{I}_{0}^{1}}),\psi_{0}(X_{0}^{\mathcal{I}_{0}^{2}}))$ quite small (Lemma E.3, Supplement E.3).

We formalise these arguments in the following theorem, whose proof (Supplement E.3) follows them very closely. The price for proof intuitiveness is that the theorem is specific to the bootstrap filter, although numerical evidence (Section 5) suggests that other filters are stable as well.