1 Introduction↩︎

1.1 Resampling↩︎

SMC algorithms (also termed particle filters) are a class of Monte Carlo methods for approximating (sequences of) probability measures and their normalising constants such as marginal likelihoods [1]. The closely related CSMC algorithm induces a Markov chain which can be used to efficiently approximate the joint posterior distribution of all the latent states in state-space models (a.k.a., the joint smoothing distribution), and more generally, due to its ability to exploit the ‘decorrelation-over-time’ structure of the target distribution [2]–[4]. CSMC methods are also at the core of more complex schemes, e.g., for unbiased estimation [5]–[7] or diffusion models [8], [9].

A key ingredient in SMC algorithms is resampling and CSMC algorithms rely on a related conditional resampling mechanism. Numerous resampling schemes have been proposed, e.g., multinomial resampling [10], stratified resampling [11], systematic resampling [12], residual resampling [13], [14], adaptive resampling [15], killing resampling [16], [17], Metropolis resampling [18], Hilbert-curve resampling [19], chopthin resampling [20], as well as further sophisticated resampling schemes in [17], [21]. For practitioners, it is therefore important to have verifiable conditions which guarantee that their chosen resampling scheme leads to SMC or CSMC algorithms that are valid. In this work, we introduce a perspective which allows us to view most of these algorithms within a common framework which provides such guarantees.

1.2 Existing frameworks↩︎

We will consider ‘validity’ of SMC and CSMC algorithms only in the sense that the former yield unbiased estimates of normalising constants (e.g., marginal likelihoods) and that the latter induce Markov kernels which leave the target distribution (e.g., joint smoothing distribution) invariant. We do not analyse the efficiency (e.g.convergence rates) associated with different resampling schemes. For the latter, we refer the reader to [14], [17], [21]–[24].

When homogeneous (i.e., equal) post-resampling weights are sought, a common requirement used to ensure validity (in the above-mentioned sense) is the conventional ‘unbiasedness’ condition from [25] which says that the expected number of descendants of each particle should be proportional to its weight.

For CSMC algorithms, the seminal work of [2] additionally required that resampling schemes should be exchangeable. They also noted that exchangeability can always be ensured by randomly permuting the order of the ancestor indices after resampling. This idea was used in [26] to derive a conditional version of residual-multinomial resampling. However, such random permutations can incur a non-negligible computational cost, especially on parallel architectures [18], and the framework from [2] does not cover adaptive resampling.

Alternatively, [26], [27] imposed a so-called marginal ‘unbiasedness’ condition which implies the conventional ‘unbiasedness’ condition but does not require exchangeability. For some non-exchangeable resampling schemes (i.e., systematic resampling and killing resampling) the authors were thus able to derive a valid CSMC algorithm by applying a random shift to the order of the ancestor indices after resampling. However, the random shift again incurs additional computational cost and their framework does not cover adaptive resampling.

It has long been known that the conventional ‘unbiasedness’ condition is not necessary and can be dropped at the expense of non-uniform post-resampling weights [28]. [29] justified this strategy as changing the sequence of intermediate distributions of the SMC or CSMC algorithm. However, this justification does not apply to non-exchangeable resampling schemes such as the partial resampling scheme from [30] or to the optimal finite-state resampling scheme from [31].

Finally, [22], [32] proved the validity of SMC and CSMC algorithms for (adaptive) resampling schemes that fit into their \(\alpha\)SMC framework. However, this framework’s resampling-matrix representation only exists for resampling schemes in which the ancestor indices are conditionally independent and thus does not apply to the widely used systematic resampling scheme nor to Metropolis resampling or Hilbert-curve resampling.

1.3 Contributions↩︎

We show how to construct valid SMC and CSMC algorithms based on a wide range of resampling schemes and in particular:

1. In Section 2, we present a simple framework for resampling which covers most known resampling schemes and their conditional versions as special cases, including those with a complicated dependence structure like (conditional) systematic resampling and those which do not satisfy the conventional ‘unbiasedness’ condition from [25] such as partial resampling [30], optimal finite-state resampling scheme [31] or a version of chopthin resampling [20].

2. Key to our framework is the introduction of an index distribution \(\lambda(\,\cdot\,|n)\) which governs the law of the particle index of the descendant of particle \(n\) on a suitably extended space. By making the rôle and the (potentially non-uniform) choice of this index distribution explicit, our work automatically provides conditional versions of any resampling scheme within our framework. For example, our work yields a novel conditional version of chopthin resampling. Further, our work makes it possible to implement conditional versions of non-exchangeable resampling schemes without having to randomly permute or shift the order of the ancestor indices.

3. The assumptions of our framework are weaker than those imposed in the existing frameworks from [2], [22], [26], [27], [32], requiring only a mild persistence condition (Definition 8) which ensures that the post-resampling weights cannot all be simultaneously zero and a proper weighting assumption (Definition 5). Additionally, we specify balanced index distribution (Definition 6) which automatically ensures proper weighting. Section 3 shows that most known resampling schemes implicitly use a balanced index distribution.

4. In Section 4, we prove that our framework makes it simple to guarantee (a) the unbiasedness of normalising-constant estimates generated by SMC algorithms; (b) the validity of the corresponding CSMC algorithm (with backward/ancestor sampling), for a wide range of resampling schemes—including adaptive/systematic resampling. In Appendix 6, we further demonstrate that naïve implementations of conditional resampling sometimes used by practitioners can bias the Markov chain induced by the CSMC algorithm and may even worsen its autocorrelation.

1.4 Notation↩︎

We let \([n] \mathrel{\vcenter{:}}= \{1,\dotsc,n\}\), for any \(n \in \mathbb{N}\mathrel{\vcenter{:}}= \{1, 2, \dotsc\}\) and use the convention \(\tfrac{0}{0} \mathrel{\vcenter{:}}= 0\). Throughout this work, we assume some background probability space \((\varOmega, \mathcal{A}, \mathop{\mathrm{\mathbb{P}}})\); \(\mathcal{B}(\mathsf{X})\) denotes the Borel \(\sigma\)-algebra on some topological space \(\mathsf{X}\). For any measure space \((\mathsf{X}, \mathcal{X}, \mu)\), where \(\mu\) may be random, we write \(\smash{\mu(f) \mathrel{\vcenter{:}}= \int_\mathsf{X}f(x) \mu(\mathrm{d}x)}\), for \(\mu\)-integrable functions \(f \colon \mathsf{X}\to \mathbb{R}\). Abusing notation, we frequently use the same symbol for measures and their densities (w.r.t.a suitable dominating measure, e.g., Lebesgue or counting measure); \(\mathop{\mathrm{vol}}(A)\) is the volume (i.e., Lebesgue measure) of measurable \(A \in \mathbb{R}^d\); \(\smash{\mathop{\mathrm{Cat}}(p_{1:K})}\) denotes the categorical distribution on \([K]\) with probability vector \(\smash{p_{1:K} \in [0,1]^K}\) such that \(\smash{\sum_{k=1}^K p_k = 1}\). For a suitable set \(B\), \(\mathop{\mathrm{Unif}}_B\) is the (discrete or continuous) uniform distribution on \(B\); \(\delta_x(B) = \mathop{\mathrm{\mathbf{1}}}\{x \in B\}\) is the point mass/Dirac measure at \(x\); if \(B = \{b\}\) is a singleton, we write \(\smash{\delta_x(\{b\}) =: \delta_x(b)}\).

2 Framework for the validity of resampling schemes↩︎

2.1 Resampling↩︎

2.1.1 Goal of resampling↩︎

In this section, we present a generic framework for (conditional) resampling. We acknowledge that this necessitates a level of abstraction whose utility will become fully apparent in Section 3, where we illustrate how most known resampling schemes fit into our framework. Due to the notational burden of SMC and CSMC algorithms, we defer these methods to Section 4 and treat (conditional) resampling in a simpler setting until then.

Hereafter, assume some measurable space \((\mathsf{X}, \mathcal{B}(\mathsf{X}))\) where \(\mathsf{X}\) is Polish, e.g., \(\smash{\mathsf{X}= \mathbb{R}^d}\). Consider a (possibly random) weighted sample \((\mathbf{X}, \mathbf{W}) \mathrel{\vcenter{:}}= (X^{1:N}, W^{1:N})\), for some \(N \in \mathbb{N}\). Here, for all \(n \in [N]\), \(X^n\) is a sample point (‘particle’) that takes values in \(\mathsf{X}\) and \(W^n\) is a non-negative weight, such that \(\sum_{n=1}^N W^n = 1\). Such a weighted sample can be associated with a probability measure on \((\mathsf{X}, \mathcal{B}(\mathsf{X}))\): \[\begin{align} \hat{\pi}\mathrel{\vcenter{:}}= \sum_{n=1}^N W^n \delta_{x^n}. \end{align}\] Resampling is a procedure that replaces \((\mathbf{X}, \mathbf{W})\) with random variables \((\widetilde{\mathbf{X}}, \widetilde{\mathbf{w}}) \mathrel{\vcenter{:}}= (\widetilde{X}^{1:M}, \widetilde{w}^{1:M})\), for some \(M \in \mathbb{N}\), which can be associated with the random measure: \[\begin{align} \tilde{\pi}\mathrel{\vcenter{:}}= \sum_{m=1}^M \widetilde{w}^m \delta_{\widetilde{X}^m}. \end{align}\] A typical goal of resampling is to ensure that \(\tilde{\pi}\) is close to \(\smash{\hat{\pi}}\) with \(\smash{\widetilde{w}^{1:M}}\) more homogeneous than \(\smash{W^{1:N}}\), in a suitable sense, whilst maintaining unbiasedness i.e., that \(\smash{\mathop{\mathrm{\mathbb{E}}}[\tilde{\pi}(f)|\mathcal{F}] = \hat{\pi}(f)}\), for any suitably integrable test function \(f\colon \mathsf{X}\to \mathbb{R}\) where \(\mathcal{F}\mathrel{\vcenter{:}}= \sigma(\mathbf{X}, \mathbf{W})\) is the \(\sigma\)-algebra generated by \((\mathbf{X}, \mathbf{W})\).

We will refer to the \(\smash{\widetilde{w}^m} \geq 0\) as post-resampling ‘weights’, though it is possible that \(\sum_{m=1}^M \widetilde{w}^m \neq 1\) (equality holds for most of the commonly used resampling schemes and must hold at least in expectation, i.e., \(\smash{\mathop{\mathrm{\mathbb{E}}}[\sum_{m=1}^M \widetilde{w}^m|\mathcal{F}]} = 1\), almost surely, for resampling schemes that are unbiased).

2.1.2 Specific class of resampling schemes treated in this work↩︎

Although resampling is very general and imposes almost no constraints on how \(\smash{(\widetilde{\mathbf{X}}, \widetilde{\mathbf{w}})}\) is produced, we focus in the sequel on the a class of resampling schemes which specify the resampled particles as \(\smash{\widetilde{X}^m \mathrel{\vcenter{:}}= X^{A^m}}\), where \(A^m\) is an ancestor index taking a value in \([N]\). Definition 1 gives a sampling distribution for \(\mathbf{A}\mathrel{\vcenter{:}}= A^{1:M}\) and for some associated optional collection of auxiliary variables \(U\) used in the simulation of \(\mathbf{A}\). The latter take values in some set \(\mathsf{U}\); if no such auxiliary variables are used, we set \(\mathsf{U}\mathrel{\vcenter{:}}= \emptyset\) and suppress \(U\) and \(u\) from the notation.

For any \(m \in [M]\), we define the following notation: \(\rho^m(a, \mathrm{d}u) \mathrel{\vcenter{:}}= \mathop{\mathrm{\mathbb{P}}}(A^m = a, U\in \mathrm{d}u)\) is the marginal distribution of \((A^m\!, U)\) under \(\rho\); \(\rho^m(a) \mathrel{\vcenter{:}}= \mathop{\mathrm{\mathbb{P}}}(A^m = a)\) is the marginal distribution of \(A^m\) under \(\rho\); \(\rho^{-m}(\mathbf{a}^{-m}|a, u) \mathrel{\vcenter{:}}= \mathop{\mathrm{\mathbb{P}}}(\mathbf{A}^{-m} = \mathbf{a}^{-m}| A^m = a, U= u)\) (assuming \(M > 1\)) is the full conditional distribution of \(\mathbf{A}^{-m} \mathrel{\vcenter{:}}= (A^1, \dotsc, A^{m-1}, A^{m+1}\!, \dotsc, A^M)\) under \(\rho\).

The key rôle of index distributions for validity of (conditional) resampling is the focus of this work.

When necessary, we will make the dependence of \(\rho\), \(\rho^m\), \(\rho^{-m}(\,\cdot\,;n)\), \(\lambda(\,\cdot\,|n)\) and \(\phi\) on the weights \(\mathbf{W}= W^{1:N}\) explicit in the notation by writing \(\rho= \rho(\,\cdot\,; \mathbf{W})\) and \(\rho^m = \rho^m(\,\cdot\,; \mathbf{W})\), \(\rho^{-m}(\,\cdot\,|n; \mathbf{W})\), \(\lambda(\,\cdot\,|n; \mathbf{W})\) and \(\phi(\,\cdot\,; \mathbf{W})\); and, in principle, these could also depend on \(\mathbf{X}\). The dependence of these quantities (and also of \(\mathsf{U}\)) on \(N\) and \(M\) will always be clear from the context.

2.1.3 Sufficient condition for unbiasedness↩︎

A frequently-imposed requirement on resampling schemes is that \(\smash{\tilde{\pi}}\) is an unbiased estimate of \(\hat{\pi}\) in the following sense, where we remember that \(\mathcal{F}= \sigma(\mathbf{X}, \mathbf{W})\) is the \(\sigma\)-algebra generated by \((\mathbf{X}, \mathbf{W})\).

The following is the main result of this section. It shows that proper weighting is sufficient for unbiasedness. Motivated by the fact that \(\smash{\hat{\pi}(f) = \mathop{\mathrm{\mathbb{E}}}_{A \sim \mathop{\mathrm{Cat}}(\mathbf{W})}[f(X^A)|\mathcal{F}]}\), the proof interprets resampling as importance sampling (with one sample point) targetting \(\mathop{\mathrm{Cat}}(\mathbf{W})\). More precisely, we perform importance sampling on the extended space that includes \((K, A, \mathbf{A}, U)\), where \(K\in [M]\) is the index of a single ancestor index and \(A \in [N]\) is its value which is almost surely equal to \(A^K\). Proper weighting then guarantees that the target distribution of this importance-sampling construction on \([M] \times [N] \times [N]^M \times \mathsf{U}\): \[\begin{align} \mu(k, a, \mathbf{a}, \mathrm{d}u) & \mathrel{\vcenter{:}}= \smash{\mathop{\mathrm{Cat}}(a; \mathbf{W}) \lambda(k, \mathrm{d}u|a) \delta_a(a^k) \rho^{-k}(\mathbf{a}^{-k}|a^k\!, u)}, \label{eq:prop:proper95weighting95implies95unbiasedness:proof:1} \end{align}\tag{2}\] (which admits \(\mathop{\mathrm{Cat}}(\mathbf{W})\) as a marginal) is absolutely continuous w.r.t.the sampling distribution \[\begin{align} q(k, a, \mathbf{a}, \mathrm{d}u) & \mathrel{\vcenter{:}}= \rho(\mathbf{a}, \mathrm{d}u) \mathop{\mathrm{Unif}}_{[M]}(k) \delta_{a^k}(a). \label{eq:prop:proper95weighting95implies95unbiasedness:proof:2} \end{align}\tag{3}\]

2.1.4 Standard choice of index distribution↩︎

The proper weighting definition leaves a degree of freedom in choosing the family of index distributions. The following class of resampling schemes provides a standard choice which automatically guarantees proper weighting by taking \(\lambda(m, \mathrm{d}u|n) \propto \rho^m(n, \mathrm{d}u)\).

The name balanced is motivated by the fact that 4 is essentially the balance heuristic for multiple importance sampling from [33] (see also [34] for further discussions on the balance heuristic and its relationship with extended state-space justifications of importance sampling similar to 2 –3 ).

Within the importance-sampling construction used in the proof of Proposition 1, 4 implies that \(\lambda(m, \mathrm{d}u|n) = q(m, \mathrm{d}u|n)\) (whenever this conditional distribution is well defined). Balanced resampling schemes have that \(\phi(m, n, u) = M W^n / \sum_{l=1}^M \rho^l(n)\) is almost-surely constant in \((m, u)\), i.e., in the proof of Proposition 1, they perform importance sampling on a smaller space which does not include \((K, U)\).

As Section 3 will demonstrate, most known resampling schemes are balanced, i.e., they implicitly use the index distribution from 4 . The definition of proper weighting implies the following result.

The following example illustrates that using a balanced resampling scheme not only guarantees proper weighting but also ensures that the conditional variance of resampling is robust to single small marginal probabilities \(\rho^m(A^m)\) under the resampling distribution.

As noted above, a frequent goal of resampling is to make the resampled weights \(\widetilde{w}^{1:M}\) more homogeneous than \(W^{1:N}\). The following class of resampling schemes achieves \(\widetilde{w}^m = 1/M\), for all \(m \in [M]\).

Since achieving \(\widetilde{w}^m = 1/M\) is frequently the goal of resampling, 5 is itself considered the unbiasedness condition [25]. We will thus call 5 the conventional ‘unbiasedness’ condition.

Persistence guarantees that the post-resampling weights \(\smash{\widetilde{w}^m}\) can (almost surely) be self-normalised so that another round of resampling could be applied to \(\smash{(\widetilde{X}^{1:M}, \widetilde{W}^{1:M})}\), where \(\smash{\widetilde{W}^m \mathrel{\vcenter{:}}= \widetilde{w}^m / \sum_{l=1}^M \widetilde{w}^l}\). The following sufficient condition for persistence is immediate from the definition.

2.2 Conditional resampling↩︎

We will now turn to conditional resampling. This is a key ingredient in CSMC methods but we present it here without reference to the latter (the extension to CSMC algorithms will be discussed in Section 4.2). Admittedly, this makes our discussion below somewhat abstract. However, we believe that abstraction is justified by the fact that doing so avoids the notational burden incurred by CSMC methods.

Throughout this section, we assume that \(M \geq 2\) and that the resampling scheme \(\langle \rho, \lambda\rangle\) is persistent (so that the post-resampling ‘weights’ can be normalised). Given some \(a \in [N]\) with \(\smash{W^a > 0}\), conditional resampling (Algorithm 2) forms a (possibly random) probability measure \[\begin{align} \tilde{\pi}^a \mathrel{\vcenter{:}}= \sum_{m=1}^M \frac{\widetilde{w}^m}{\sum_{l=1}^M \widetilde{w}^l} \delta_{\widetilde{X}^m}, \end{align}\] in such a way that the support of \(\tilde{\pi}^a\) includes \(X^a\).

The following requirement says that if we first sample \(A = a\) with probability \(W^a\) and then perform conditional resampling, \(\smash{\tilde{\pi}^a}\) should again induce a random measure that is unbiased for \(\hat{\pi}\).

The terminology ‘invariance’ is motivated by the fact that Definition 9 holds if executing Algorithm 2 and then returning the ancestor index \(A^L\) with probability proportional to \(\widetilde{w}^L\) induces a \(\mathop{\mathrm{Cat}}(\mathbf{W})\)-invariant Markov kernel. The main result of this section is the following.

2.3 Implementation of conditional resampling with fixed reference index↩︎

In practice, conditional resampling (within CSMC algorithms) is often implemented by placing the ancestor index \(a\) in a fixed location \(K = 1\) or \(K = M\), say, in Line [enum:alg:conditional95resampling:sample95k] of Algorithm 2. While this strategy is generally invalid (invariance is lost), the next proposition (proved in Appendix 8) shows that it is valid when the resampling distribution is symmetric in the following sense.

The reindexability assumption of Proposition 5 is satisfied for exchangeable resampling schemes.

3 Examples↩︎

We now demonstrate that a number of well known resampling schemes \(\langle \rho, \lambda\rangle\) can be viewed as special cases of our framework. In particular, this allows us to present valid conditional versions of popular resampling schemes that are not exchangeable, without the need for randomly permuting or shifting the order of the generated ancestor indices.

We start with a number of elementary resampling schemes (multinomial, stratified, systematic or killing resampling and even ‘trivial’ resampling, i.e., not resampling). Then, we discuss composite resampling schemes which build more sophisticated resampling schemes from elementary ones.

3.1 Elementary resampling schemes↩︎

The popular elementary resampling schemes discussed in this section are persistent and balanced and all (except ‘trivial’ resampling) are even weight-homogeneous.

3.1.1 Multinomial resampling↩︎

Multinomial resampling [10]—also known as stochastic sampling with replacement in the genetic algorithms literature [35]—takes \(\mathsf{U}\mathrel{\vcenter{:}}= \emptyset\) and \(\rho(\mathbf{a}) \mathrel{\vcenter{:}}= \smash{\prod_{m=1}^M \rho^m(a^m)}\), where \(\smash{\rho^m \mathrel{\vcenter{:}}= \mathop{\mathrm{Cat}}(\mathbf{W})}\), and \(\smash{\lambda(m|n) \mathrel{\vcenter{:}}= \mathop{\mathrm{Unif}}_{[M]}(m)}\) and hence \(\phi\equiv 1\). Clearly, multinomial resampling is weight-homogeneous (and hence balanced and persistent) and also exchangeable (so that Proposition 5 would allow us to simply set \(K \mathrel{\vcenter{:}}= 1\), say, in Algorithm 4). Line [alg:multinomial95resampling:2] of Algorithm 3 and Line [alg:conditional95multinomial95resampling:2] of Algorithm 4 can be implemented in \(\mathop{\mathrm{\mathcal{O}}}(\max\{M, N\})\) operations using aliasing methods [36], [37].

3.1.2 Stratified resampling↩︎

Stratified resampling [11] takes \(\mathsf{U}\mathrel{\vcenter{:}}= [0,1]^M\) and, writing \(U= U^{1:M}\): \[\begin{align} \rho(\mathbf{a}, \mathrm{d}u) & \mathrel{\vcenter{:}}= \prod_{m=1}^M \mathop{\mathrm{Unif}}_{[0,1]}(\mathrm{d}u^m) \delta_{F^{-1}([u^m + m - 1] / M )}(a^m), \end{align}\] where \(F(n) \mathrel{\vcenter{:}}= \sum_{l=1}^n W^l\), for \(n \in [N]\) and its generalized invers \(F^{-1}(p) \mathrel{\vcenter{:}}= \min\{n \in [N] \mid p \leq F(n)\}\), for \(p \in [0,1]\). For what follows, it will be convenient to use that for any \(u \in [0, 1]\): \[\begin{align} n = F^{-1}\biggl(\frac{u+m-1}{M}\biggr) & \Longleftrightarrow u \in I^m(n), \end{align}\] for any \((m, n, u) \in [M] \times [N] \times [0,1]\), where \(\{I^m(n)\}_{n \in [N]}\) is the partition of \((0, 1]\) given by \[\begin{align} I^m(n) & \mathrel{\vcenter{:}}= (\max\{0, M F(n-1) -m +1\}, \min\{1, M F(n) - m +1\}]. \end{align}\] with convention \(F(0) \mathrel{\vcenter{:}}= 0\). With this notation, we can express the resampling distribution equivalently as \(\smash{\rho(\mathbf{a}, \mathrm{d}u) \mathrel{\vcenter{:}}= \prod_{m=1}^M \rho^m(a^m) \mathop{\mathrm{Unif}}_{I^m(a^m)}(\mathrm{d}u^m)}\), where \[\begin{align} \rho^m(n) = \mathop{\mathrm{Cat}}(n; (\mathop{\mathrm{vol}}(I^m(l)))_{l \in [N]}). \label{eq:stratified95resampling95marginal} \end{align}\tag{6}\] Stratified resampling is not exchangeable. However, \(\rho\) satisfies the conventional ‘unbiasedness’ condition 5 : \[\begin{align} \sum_{m=1}^M \rho^m(n) = \sum_{m=1}^M \mathop{\mathrm{vol}}((\max\{m-1, M F(n-1)\}, \min\{m, M F(n)\}]) = M W^n. \end{align}\] Consequently, stratified resampling is weight-homogeneous (and hence balanced and thus properly weighted as well as persistent) with \[\begin{align} \lambda(m, u|n) & \mathrel{\vcenter{:}}= \frac{\rho^m(n, u)}{M W^n}\\ & = \mathop{\mathrm{Cat}}(m; (\mathop{\mathrm{vol}}(I^l(n)) / (MW^n))_{l \in [M]}) \mathop{\mathrm{Unif}}_{I^m(n)}(u^m) \smash{ \prod_{\substack{\mathllap{l} = \mathrlap{1}\\ \mathllap{l} \neq \mathrlap{m}}}^{\smash{M}}} \mathop{\mathrm{Unif}}_{[0,1]}(u^l). \label{eq:lambda95stratified95resampling} \end{align}\tag{7}\]

Note that in Algorithm 6, by construction, \(\smash{A^K = a}\) with probability \(1\). It is well known that Line [alg:stratified95resampling:2] in Algorithm 5 and Line [alg:conditional95stratified95resampling:2] in Algorithm 6 can be implemented in \(\mathop{\mathrm{\mathcal{O}}}(\max\{M, N\})\) operations as outlined in Algorithm 7.

3.1.3 Systematic resampling↩︎

Systematic resampling [12]—also known as stochastic universal sampling in genetic algorithms [35]—differs from stratified resampling in that it shares the same uniform random variable \(U\) across all strata: \[\begin{align} \rho(\mathbf{a}, u) \mathrel{\vcenter{:}}= \mathop{\mathrm{Unif}}_{[0,1]}(u) \prod_{m=1}^M \delta_{F^{-1}([u + m - 1] / M )}(a^m), \end{align}\] with \(\mathsf{U}\mathrel{\vcenter{:}}= [0,1]\). Thus, \(\smash{\rho^m(n) = \mathop{\mathrm{vol}}(I^m(n))}\) as in stratified resampling 6 so that systematic resampling again satisfies the conventional ‘unbiasedness’ condition 5 ; and setting \[\begin{align} \lambda(m, u|n) \mathrel{\vcenter{:}}= \frac{\rho^m(n, u)}{M W^n} = \mathop{\mathrm{Cat}}(m; (\mathop{\mathrm{vol}}(I^l(n)) / (MW^n))_{l \in [M]}) \mathop{\mathrm{Unif}}_{I^m(n)}(u), \label{eq:lambda95systematic95resampling} \end{align}\tag{8}\] ensures that systematic resampling is weight-homogeneous (and hence balanced and thus properly weighted as well as persistent).

Previously, [26] derived a version of conditional systematic resampling which coincides with Algorithm 9 up to an additional random shift (see Remark 2).

3.1.4 Killing resampling↩︎

Let \(M = N\). Killing resampling [17] [38] takes \(\mathsf{U}\mathrel{\vcenter{:}}= \emptyset\), \(\smash{\rho(\mathbf{a})} \mathrel{\vcenter{:}}= \smash{\prod_{m=1}^M \rho^m(a^m)}\), where, writing \(\smash{W^\star} \mathrel{\vcenter{:}}= \smash{\max \{W^1, \dotsc, W^N\}}\), \[\begin{align} \rho^m(n) & \mathrel{\vcenter{:}}= \frac{W^m}{W^\star} \delta_m(n) + \Bigl(1 - \frac{W^m}{W^\star}\Bigr) \mathop{\mathrm{Cat}}(n; \mathbf{W}). \end{align}\] It can be verified that \(\rho\) thus satisfies the conventional ‘unbiasedness’ condition 5 . Thus, taking \[\begin{align} \lambda(m|n) & \mathrel{\vcenter{:}}= \frac{\rho^m(n)}{M W^n} = \frac{W^\star - W^n +1}{MW^\star} \mathop{\mathrm{\mathbf{1}}}\{m = n\} + \frac{W^\star - W^m}{MW^\star} \mathop{\mathrm{\mathbf{1}}}\{m \in [M] \setminus \{n\}\}, \label{eq:killing95resampling95lambda} \end{align}\tag{9}\] means that killing resampling is weight-homogeneous (and hence properly weighted and persistent).

Previously, [27] derived a version of conditional killing resampling which coincides with Algorithm 11 up to an additional random shift (see Remark 2).

3.1.5 Trivial resampling (i.e., not resampling)↩︎

We now show that ‘not resampling’ can be viewed as a resampling scheme (termed trivial resampling) within our framework. This will become useful in Section 3.2.3 where justify adaptive resampling as mixing trivial resampling with another, ‘non-trivial’ resampling scheme.

Let \(M \mathrel{\vcenter{:}}= N\). Then trivial resampling can be formalised as \(\mathsf{U}\mathrel{\vcenter{:}}= \emptyset\), \(\rho(\mathbf{a}) \mathrel{\vcenter{:}}= \smash{\prod_{m=1}^M \rho^m(a^m)}\), where \(\smash{\rho^m \mathrel{\vcenter{:}}= \delta_{m}}\). With \(\smash{\lambda(\,\cdot\,|n) \mathrel{\vcenter{:}}= \delta_n}\), trivial resampling is then balanced (and hence properly weighted) and persistent, and has \(\smash{\phi(m, n, u)} = \smash{M W^n \mathop{\mathrm{\mathbf{1}}}\{m = n\}}\).

3.2 Composite resampling schemes↩︎

Let \(\langle \varrho, \kappa\rangle\) denote some properly weighted and persistent resampling scheme (e.g., one of the elementary resampling schemes from the previous section). In this section, we illustrate that \(\langle \varrho, \kappa\rangle\) can be modified to obtain a new properly weighted and persistent resampling scheme \(\langle \rho, \lambda\rangle\) with potentially more favourable properties. This permits a simple justification of many sophisticated resampling schemes that have been proposed in the literature such as residual resampling, resampling with weight re-ordering, adaptive resampling, optimal finite-state resampling and chopthin resampling. In particular, our framework shows how to implement conditional resampling for each of these.

3.2.1 Residual resampling↩︎

Let \(\smash{M_0 \mathrel{\vcenter{:}}= \sum_{n=1}^N \lfloor M W^n \rfloor}\) and \[\begin{align} \smash{\alpha = \alpha^{1:M_0} = (\overbrace{1, \dotsc, 1}^{\mathclap{\text{\lfloor M W^1 \rfloor}\atop\text{times}}}, \overbrace{2, \dotsc, 2}^{\mathclap{\text{\lfloor M W^2 \rfloor}\atop\text{times}}}, \dotsc, \overbrace{N, \dotsc, N}^{\mathclap{\text{\lfloor M W^N \rfloor}\atop\text{times}}})}. \end{align}\] Additionally, let \(\langle \varrho, \kappa\rangle\) be some weight-homogeneous resampling scheme (e.g.any of the elementary resampling schemes from Section 3.1 except for ‘trivial’ resampling) which generates \(M - M_0\) ancestor indices on \([N]\) based on weights \(\mathbf{V} \mathrel{\vcenter{:}}= V^{1:N}\), where \[\begin{align} V^n & \mathrel{\vcenter{:}}= \frac{M W^n - \lfloor M W^n \rfloor}{M - M_0}. \end{align}\] Residual resampling then takes the following resampling distribution (other orderings of \(\smash{\mathbf{a}= a^{1:M}}\) are possible and can facilitate ‘in place’ implementations): \[\begin{align} \rho(\mathbf{a}, u) \mathrel{\vcenter{:}}= \delta_{\alpha}(a^{1:M_0}) \smashoperator{\sum_{\mathbf{b} \in [N]^{M-M_0}}}\varrho(\mathbf{b}, u; \mathbf{V}) \delta_{\mathbf{b}}(a^{(M_0+1):M}). \end{align}\] Consequently, \(\smash{\rho^m(n) = \mathop{\mathrm{\mathbf{1}}}\{m \in [M_0]\} \delta_{\alpha^m}(n) + \mathop{\mathrm{\mathbf{1}}}\{m \in [M] \setminus [M_0]\} \varrho^{m-M_0}(n; \mathbf{V})}\). It can be verified that \(\rho\) satisfies the conventional ‘unbiasedness’ condition 5 . Thus, setting \[\begin{align} \lambda(m, u|n) & = \frac{\rho^m(n, u)}{M W^n} = p^n \mathop{\mathrm{Unif}}_{\{l \in [M_0] \mid \alpha^l = n\}}(m) \varrho(u) + (1 - p^n) \kappa(m-M_0, u|n; \mathbf{V}), \end{align}\] where \(\smash{p^n \mathrel{\vcenter{:}}= \lfloor M W^n \rfloor / (MW^n)}\), makes the residual resampling scheme \(\langle \rho, \lambda\rangle\) is weight-homogenous. If \(\langle \varrho, \kappa\rangle\) is multinomial resampling then \(\langle \rho, \lambda\rangle\) is called residual-multinomial resampling [13]—also called remainder stochastic sampling with replacement in genetic algorithms [35]—but other choices, e.g., residual-systematic resampling [39] are possible.

Previously, [26] derived a version of conditional residual-multinomial resampling which coincides with Algorithm 15 up to an additional random shift (see Remark 2).

3.2.2 Re-ordering weights prior to resampling↩︎

Let \(\langle \varrho, \kappa\rangle\) be a balanced and persistent resampling scheme. If \(\langle \varrho, \kappa\rangle\) is non-exchangeable, then it can sometimes be desirable to re-order the weights before resampling, e.g., to improve convergence rates [11], [14], [21].

More formally, let \(\varpi \in \mathsf{S}_N\) be some \(\mathcal{F}\)-measurable permutation (we only consider deterministic weight reordering for simplicity but random reorderings would be possible, too). Furthermore, write \(\smash{\mathbf{W}^\varpi} = \smash{(W^{\varpi(1)}, \dotsc, W^{\varpi(N)})}\). We also now include the (order of the) weights explicitly into the resampling distribution. Then, sampling \(\smash{(\mathbf{B}, U)} = \smash{(B^{1:M}, U)} \sim \smash{\varrho(\,\cdot\,; \mathbf{W}^\varpi)}\)—and subsequently replacing the \(m\)th generated ancestor index \(\smash{B^m}\) by \(\smash{A^m \mathrel{\vcenter{:}}= \varpi(B^m)}\)—induces another resampling scheme \(\langle \rho, \lambda\rangle\) with \[\begin{align} \rho(\mathbf{a}, u) & \mathrel{\vcenter{:}}= \smash{\sum_{\mathbf{b} \in [N]^M} \varrho(\mathbf{b}, u; \mathbf{W}^\varpi) \delta_{\varpi(\mathbf{b})}(\mathbf{a}) = \varrho(\varpi^{-1}(\mathbf{a}), u; \mathbf{W}^\varpi),}\vphantom{\prod_i}\\ \lambda(m, u|n) & \mathrel{\vcenter{:}}= \frac{\varrho^m(\varpi^{-1}(n); \mathbf{W}^\varpi)}{\sum_{l=1}^M \varrho^l(\varpi^{-1}(n); \mathbf{W}^\varpi)} \varrho^m(u| \varpi^{-1}(n); \mathbf{W}^\varpi) = \kappa(m, u|\varpi^{-1}(n), \mathbf{W}^\varpi), \end{align}\] where, for any \(\smash{\mathbf{b} \mathrel{\vcenter{:}}= b^{1:N} \in [N]^N}\) and any permutation \(\sigma \in \mathsf{S}_N\), we write \(\smash{\sigma(\mathbf{b}) \mathrel{\vcenter{:}}= (\sigma(b^1), \dotsc, \sigma(b^N))}\). By assumption on \(\langle \varrho, \kappa\rangle\), \(\langle \rho, \lambda\rangle\) is then balanced (and hence properly weighted) and persistent.

In particular, note that if \(\langle \varrho, \kappa\rangle\) is weight-homogeneous (which holds for all the elementary resampling schemes discussed in Section 3.1 except for ‘trivial’ resampling) then \(\langle \rho, \lambda\rangle\) is again weight-homogeneous. If \(\langle \varrho, \kappa\rangle\) is the stratified or systematic resampling scheme and if \(\varpi\) is a mean partition order for \(\mathbf{W}\), i.e., \(W^{\varpi(1)}, \dotsc, W^{\varpi(N_0)} \leq \frac{1}{N}\sum_{n=1}^N W^n\), and \(W^{\varpi(N_0)}, \dotsc, W^{\varpi(N)} > \frac{1}{N}\sum_{n=1}^N W^n\), for some \(N_0 \in [N]\), then Algorithm 16 recovers the mean-order versions of stratified/systematic resampling proposed in [17]; [27] recently justified its conditional version by using random shifts.

3.2.3 Mixture resampling and adaptive resampling↩︎

For some \(I \in \mathbb{N}\), let \(\{\langle \varrho_i, \kappa_i\rangle\}_{i \in [I]}\) be family of balanced and persistent resampling schemes each of which may make use of some auxiliary variables \(U_i\) taking values in some space \(\mathsf{U}_i\). Additionally, let \(V \in [I]\) be some random indicator whose distribution \(\rho(v) = \rho(v; \mathbf{W})\) may depend on the weights \(\smash{\mathbf{W}= W^{1:N}}\) although we do not make this dependence explicit in the notation.

We will now specify a mixture resampling scheme \(\langle \rho, \lambda\rangle\) which performs resampling according to \(\langle \varrho_v, \kappa_v\rangle\) with probability \(\rho(v)\): Write \(\smash{\mathsf{U}\mathrel{\vcenter{:}}= [I] \times \bigtimes_{i=1}^I \mathsf{U}_i}\) as well as \(\smash{U\mathrel{\vcenter{:}}= (V, U_{1:I})}\) and define \[\begin{align} \rho(\mathbf{a}, u) & \mathrel{\vcenter{:}}= \rho(v) \varrho_{v}(\mathbf{a}, u_{v}) \smash{ \prod_{\substack{\mathllap{i} = \mathrlap{1}\\ \mathllap{i} \neq \mathrlap{v}}}^I \varrho_i(u_i),}\vphantom{\prod^I} \end{align}\] and (whenever \(\smash{W^n > 0}\)):
\[\begin{align} \lambda(m, u| n) & \mathrel{\vcenter{:}}= \smash{\frac{\smash{\rho(v) \varrho_{v}^m(n, u_v)}}{\sum_{l=1}^M \varrho_{v}^l(n)} \prod_{\substack{\mathllap{i} = \mathrlap{1}\\ \mathllap{i} \neq \mathrlap{v}}}^{\smash{I}} \varrho_i(u_i) = \rho(v) \kappa_{v}(m, u_{v}|n) \prod_{\substack{\mathllap{i} = \mathrlap{1}\\ \mathllap{i} \neq \mathrlap{v}}}^{\smash{I}} \varrho_i(u_i).}\vphantom{\prod_{i_{i_i}}^I} \end{align}\] Then \(\langle \rho, \lambda\rangle\) is properly weighted and persistent but not generally balanced (unless \(\rho(i) = 1\) for some \(i \in [I]\)). To see this, note that \(\phi(m, n, u) = \smash{M W^n / \sum_{l=1}^M \varrho_{v}^l(n)}\) still depends on the component \(v\) of \(u\). We could have made \(\langle \rho, \lambda\rangle\) balanced by setting \(\lambda(m, \mathrm{d}u|n) \propto \rho^m(n, \mathrm{d}u)\). However, then \(\smash{\phi(m, n, u) = M W^n / \sum_{i=1}^I \rho(i) \sum_{l=1}^M \varrho_{i}^l(n) }\), i.e., the cost of evaluating post-resampling weights grows with the number of mixture components, \(I\).

A special case of such a mixture resampling scheme is frequently employed within SMC methods (to be discussed in Section 4). Assume that \(M = N\), \(I = 2\), and that \(\langle \varrho_2, \kappa_2 \rangle\) is the ‘trivial’ resampling scheme from Section 3.1.5. A common choice, first proposed in [15], is then to decide ‘whether or not to resample’ based on some \(\mathcal{F}\)-measurable function such as the ESS of the weights [40]: \[\begin{align} \mathop{\mathrm{ess}}(\mathbf{W}) \mathrel{\vcenter{:}}= \smash{\biggl[\sum_{n=1}^N (W^n)^2\biggr]^{-1},} \vphantom{\biggl[\biggr]} \end{align}\] which corresponds to the degenerate distribution \[\begin{align} \rho(v) = \rho(v; \mathbf{W}) \mathrel{\vcenter{:}}= \mathbf{1}\{\mathop{\mathrm{ess}}(\mathbf{W}) \leq \eta N\} \delta_1(v) + \mathbf{1}\{\mathop{\mathrm{ess}}(\mathbf{W}) > \eta N\} \delta_2(v), \end{align}\] for some user-defined threshold \(\eta \in [0, 1]\). In the remainder of this work, we refer to this special case of mixture resampling as ESS-adaptive resampling.

3.2.4 Partitioned resampling, partial resampling and optimal finite-state resampling↩︎

Let \(I \in \mathbb{N}\) and for \(i \in [I]\), let \(\smash{L_i \mathrel{\vcenter{:}}= \{l_i^1, \dotsc, l_i^{N_i}\} \subseteq [N]}\), be such that \(\smash{\{L_1, \dotsc, L_I\}}\) is a partition of \([N]\). Then we can ‘partition’ the weighted sample \((\mathbf{X}, \mathbf{W})\) into \(I\) disjoint blocks of weighted samples of the form \(\smash{(\mathbf{X}_i, \mathbf{W}_i) \mathrel{\vcenter{:}}= (X_i^{1:N_i}, W_i^{1:{N_i}})}\), where for all \(n \in [N_i]\): \[\begin{align} X_i^n \mathrel{\vcenter{:}}= X^{l_i^n}, \qquad \text{and} \qquad W_i^n \mathrel{\vcenter{:}}= \smash{\frac{W^{l_i^n}}{\sum_{m=1}^{N_i} W^{l_i^m}}.}\vphantom{\prod_i^I} \end{align}\] For \(i \in [I]\), let \(\langle \varrho_i, \kappa_i \rangle\) be some balanced and persistent resampling scheme which generates \(M_i\) ancestor indices \(\smash{\mathbf{A}_i = A_i^{1:M_i} \in [N_i]^{M_i}}\) for the weighted sample \((\mathbf{X}_i, \mathbf{W}_i)\) and may make use of some auxiliary variables \(U_i\) taking values in some space \(\mathsf{U}_i\).

We will now specify a partitioned resampling scheme \(\langle \rho, \lambda\rangle\) which generates ancestor indices \(\smash{\mathbf{A}\in \bigtimes_{i=1}^I L_i^{M_i}}\) from \(\mathbf{A}_1, \dotsc, \mathbf{A}_I\) by setting \[\begin{align} \smash{\mathbf{A}\mathrel{\vcenter{:}}= h(\mathbf{A}_1, \dotsc, \mathbf{A}_I) \mathrel{\vcenter{:}}= (l_1^{A_1^1}, \dotsc, l_1^{A_1^{M_1}}, l_2^{A_2^1}, \dotsc, l_2^{A_2^{M_2}}, \dotsc, l_I^{A_I^1}, \dotsc, l_I^{A_I^{M_I}}).} \end{align}\] More precisely, we set \(\smash{M \mathrel{\vcenter{:}}= \sum_{i=1}^I M_i}\), \(\smash{\mathsf{U}\mathrel{\vcenter{:}}= \bigtimes_{i=1}^I \mathsf{U}_i}\), \(\smash{U\mathrel{\vcenter{:}}= U_{1:I}}\) and \[\begin{align} \smash{\rho(\mathbf{a}, u)} & \mathrel{\vcenter{:}}= \smash{\sum_{(\mathbf{a}_1, \dotsc, \mathbf{a}_I) \in \bigtimes_{i=1}^I [N_i]^{M_i}} \delta_{h(\mathbf{a}_1, \dotsc, \mathbf{a}_I)}(\mathbf{a}) \prod_{i=1}^I \varrho_i(\mathbf{a}_i, u_i),} \vphantom{\prod_i^i}\\ \smash{\lambda(m, u| n)} & \mathrel{\vcenter{:}}= \smash{\frac{\rho^m(n, u)}{\sum_{m=1}^M \rho^m(n)} = \kappa_{i(n)}\Bigl(m - \smashoperator{\sum_{j=1}^{i(n)-1}} M_j, u_{i(n)} \,\Big|\,q(n)\Bigr) \prod_{\substack{\mathllap{j} = \mathrlap{1}\\ \mathllap{j} \neq \mathrlap{i(n)}}}^I \varrho_j(u_j),} \vphantom{\prod_{I_I}^i}\\ \phi(m, n, u) & = \frac{W^n \lambda(m, u| n)}{M^{-1}\rho^m(n, u)} = \frac{M W^n}{\sum_{l=1}^{M_{i(n)}} \varrho_{i(n)}^l(q(n))}, \end{align}\] where \(i(n) \in [I]\) and \(q(n) \in \smash{[N_{i(n)}]}\) are the unique solutions to \(\smash{n \in L_{i(n)}}\) and \(\smash{n = l_{i(n)}^{q(n)}}\), respectively. Then \(\langle \rho, \lambda\rangle\) is again balanced (and thus properly weighted) and persistent.

We end this section with two special cases of partitioned resampling that have been proposed in the literature. Let \(I = 2\) and let \(\langle \varrho_1, \kappa_1 \rangle\) be the trivial resampling scheme from Section 3.1.5 (with \(M_1 = N_1\)) and let \(\langle \varrho_2, \kappa_2 \rangle\) be some other weight-homogeneous (for simplicity) resampling scheme.

Then \(\langle \rho, \lambda\rangle\) reduces to the partial resampling scheme from [30] (a related method was proposed in [41]) which resamples only the particles with indices in \(L_2\) but keeps those with indices in \(L_1\) unchanged. Weight-homogeneity of \(\smash{\langle \varrho_2, \kappa_2 \rangle}\) then implies \(\smash{\sum_{m=1}^{M_2} \varrho_{2}^m(n) = M_2 W_2^n}\), for any \(n \in [N_2]\), so that \[\begin{align} \phi(m, n, u) & = \frac{MW^n}{\sum_{l=1}^{M_{i(n)}} \varrho_{i(n)}^l(q(n))} = \begin{dcases} \smash{MW^n}, & \text{if i(n) = 1,}\\ \frac{\smash{M W^n}}{M_2 W_2^{q(n)}} = \frac{M}{M_2} \sum_{j=1}^{\smash{N_2}} W^{l_2^j}, & \text{if i(n) = 2.} \end{dcases} \label{eq:partial95resampling95post95resampling95weights} \end{align}\tag{10}\]

A sophisticated version of such a partial resampling scheme known as optimal finite-state resampling was proposed in [31], [42]; and its conditional version in [43]. It uses [42] to solve \[\begin{align} \sum_{n=1}^{N} \min\{1, C W^n\} = M, \label{eq:implicit95definition95of95c} \end{align}\tag{11}\] for \(C > 0\) and then sets \(\smash{L_1} \mathrel{\vcenter{:}}= \smash{\{n \in [N] \mid C W^n \geq 1\}}\) and \(L_2 \mathrel{\vcenter{:}}= \smash{\{n \in [N] \mid C W^n < 1\}}\). By 11 , we then have \(\smash{\frac{1}{M_2} \sum_{n=1}^{N_2} W^{l_2^n} = 1 / C}\) so that 10 takes the following simple form: \(\smash{\phi(m, n, u)} = \smash{M W^n / \min\{1, C W^n\}}\). If \(M \leq N\), and if \(\langle \varrho_2, \kappa_2\rangle\) is stratified or systematic resampling, this scheme ensures that the family of ancestor indices \(\mathbf{A}\) contains no duplicates. This minimises the variance induced by the resampling operation when \(\mathsf{X}\) is finite which is why optimal finite-state resampling is at the heart of the discrete particle filters from [31], [42] and related SMC and CSMC algorithms for change-point models [44], [45].

3.2.5 Modified-weight and chopthin resampling↩︎

Consider a weight-homogeneous resampling scheme which generates \(M\) ancestor indices using weights \(W^{1:N}\). Then by the conventional ‘unbiasedness’ condition 5 , the expected number of descendants of the \(n\)th particle is \(M \smash{W^n}\). More generally, we may wish to enforce that the expected number of descendants of the \(n\)th particle is \(\smash{H^n} \neq M\smash{W^n}\), for some \(\smash{H^1, \dotsc, H^N} \geq 0\) such that \[\begin{align} \forall\, n \in [N]: H^n > 0 \text{ if and only if } W^n > 0; \qquad \text{and} \qquad \sum_{n=1}^N H^n = M. \label{eq:modified95weight95resampling95assumption} \end{align}\tag{12}\] This can be achieved as follows. Let \(\langle \varrho, \kappa\rangle\) be some weight-homogeneous resampling scheme which generates \(M\) ancestors based on the weights \(\smash{V^{1:N}}\), where \(\smash{V^n = H^n / M}\). Then, taking \[\begin{align} \rho(\mathbf{a}, u) \mathrel{\vcenter{:}}= \smash{\varrho(\mathbf{a}, u; V^{1:N})}, \qquad \lambda(m, u|n) \mathrel{\vcenter{:}}= \smash{\kappa(m, u|n; V^{1:N})},\label{eq:modified95weight95resampling} \end{align}\tag{13}\] i.e., \(\phi(m,n,u) = \smash{M W^n / H^n}\), gives a generally no longer weight-homogeneous but still balanced and persistent resampling scheme \(\langle \rho, \lambda\rangle\) which satisfies 12 . Resampling according to alternative weights has a long history (e.g., [28], [29], [46]).

As an example of this idea, [20] sought a resampling scheme such that:

1. particles whose weights are below some threshold \(\alpha > 0\) have zero or one descendants (‘thin’);

2. particles whose weights are equal to or above \(\alpha\) have one or more offspring (‘chop’);

3. the variability of the post-resampling ‘weights’ can be controlled, i.e., for some \(\tilde{\beta} \geq 1\) chosen by the user: \(\smash{\max_{(k,l) \in [M]^2}\widetilde{w}^k / \widetilde{w}^l \leq \tilde{\beta}}\).

These goals could be ensured by Algorithm 22 or 23 with systematic resampling as the weight-homogeneous resampling scheme \(\langle \varrho, \kappa\rangle\) and \(\smash{H^n} \mathrel{\vcenter{:}}= \smash{h_{\alpha^{\mathrlap{*}},\beta}(W^n)}\) defined through the continuous mapping \[\begin{align} [0,1] \ni W \mapsto h_{\alpha, \beta}(W) \mathrel{\vcenter{:}}= \begin{cases} W / \alpha, & \text{if W < \alpha,}\\ 1, & \text{if \alpha \leq W < \alpha \beta / 2,}\\ 2 W / (\alpha \beta), & \text{if \alpha \beta / 2 \leq W.} \end{cases} \end{align}\] Above, given \(\beta \geq 4\), \(\alpha^* \in (0,1]\)—found via [20]—solves \[\begin{align} \forall\, n \in [N]: h_{\alpha^{\mathrlap{*}},\beta}(W^n) > 0 \text{ if and only if } W^n > 0; \qquad \text{and} \qquad \sum_{n=1}^N h_{\alpha^{\mathrlap{*}}, \beta}(W^n) = M. \label{eq:chopthin95condition95on95alpha} \end{align}\tag{14}\] Thus, \(\smash{\widetilde{w}^m = W^{A^m} / H^{A^m}}\), i.e., \(\widetilde{w}^m = \alpha^*\), if \(\smash{W^{A^m} < \alpha^*}\), and \(\smash{\widetilde{w}^m = \min\{W^{A^m}, \alpha^* \beta / 2\}}\), if \(\smash{W^{A^m} \geq \alpha^*}\).

However, [20] instead set \(\smash{\widetilde{w}^m \mathrel{\vcenter{:}}= W^{A^m} / \sum_{l=1}^M \mathop{\mathrm{\mathbf{1}}}\{A^l = A^m\}}\), if \(\smash{W^{A^m} \geq \alpha^*}\). We now justify this within our framework using an enlarged set of auxiliary variables and an imbalanced index distribution. Write \(u\mathrel{\vcenter{:}}= (u_0, \mathbf{b})\in \smash{\mathsf{U}_0 \times [N]^M \eqqcolon \mathsf{U}}\) with \(u_0 \in \mathsf{U}_0\) and \(\mathbf{b} \mathrel{\vcenter{:}}= \smash{b^{1:M} \in [N]^M}\) being, respectively, the auxiliary variables and ancestor indices generated by \(\smash{\varrho(\,\cdot\,; V^{1:N})}\) and set: \[\begin{align} \rho(\mathbf{a}, u) & \mathrel{\vcenter{:}}= \smash{\varrho(u; V^{1:N})} \delta_{\mathbf{b}}(\mathbf{a}), \tag{15}\\ \lambda(m, u|n) & \mathrel{\vcenter{:}}= \begin{cases} \frac{\rho^m(n, u)}{\sum_{l=1}^M \rho^l(n)} = \smash{\kappa(m, u_0|n; V^{1:N}) \delta_{n}(b^m) \varrho^{-m}(\mathbf{b}^{-m}|b^m, u_0; V^{1:N})}, & \text{if \smash{W^n < \alpha^*},}\\ \smash{\varrho(u; V^{1:N}) \mathop{\mathrm{Unif}}_{\{k \in [M]\mid b^k = n\}}(m)}, & \text{if \smash{W^n \geq \alpha^*},} \end{cases} \tag{16}\\ \phi(m,n,u) & = \frac{W^n \lambda(m, u|n)}{M^{-1}\rho^m(n,u)} = \begin{cases} \smash{M \alpha^* \mathop{\mathrm{\mathbf{1}}}\{b^m = n\}}, & \text{if \smash{W^n < \alpha^*},}\\ \smash{M W^n \mathop{\mathrm{\mathbf{1}}}\{b^m = n\}} / \textstyle \smash{\sum_{l=1}^M \mathop{\mathrm{\mathbf{1}}}\{b^l = n\}}, & \text{if \smash{W^n \geq \alpha^*}.} \end{cases}\tag{17} \end{align}\] For intuition, note that Particle \(n\) is guaranteed to have at least one descendant under \(\rho\) whenever \(\smash{W^n \geq \alpha}\); thus, we have greater flexibility in choosing \(\lambda(m,u|n)\) whenever \(\smash{W^n \geq \alpha}\).

The resampling scheme \(\langle \rho, \lambda\rangle\), summarised in Algorithm 24, is then (the version of) the chopthin resampling scheme from [20] which was demonstrated by the authors to outperform systematic resampling in some applications. Having cast chopthin resampling within our framework immediately shows that it is properly weighted and hence unbiased (and also persistent); it also allows us to derive the corresponding conditional version, presented in Algorithm 25, which is novel.

A later version of chopthin resampling [47] differs from the one of [20] in that it additionally (a) uses weight reordering with a permutation \(\varpi\) such that \(W^{\varpi(1)}, \dotsc, W^{\varpi(N_0)} < \alpha^*\), and \(W^{\varpi(N_0)}, \dotsc, W^{\varpi(N)} \geq \alpha^*\), for some \(N_0 \in [N]\), and (b) ensures that the post-resampling weights are again normalised, i.e., \(\sum_{m =1}^M \widetilde{w}^m = 1\). The latter is achieved by adjusting the post-resampling weights as \(\smash{\widetilde{w}^m = (W^{A^m} + \xi^{A^m})/ \sum_{l=1}^M \mathop{\mathrm{\mathbf{1}}}\{A^l = A^m\}}\), if \(\smash{W^{A^m} \geq \alpha^*}\), where \(\xi^n\) is a certain random variable satisfying \(\mathop{\mathrm{\mathbb{E}}}[\xi^n|\mathcal{F}] = 0\) and \(\xi^n - W^n \geq 0\), almost surely. Whether (b) can be justified within our framework (e.g., in order to derive a conditional version) remains an open question.

4 Application to sequential Monte Carlo methods↩︎

In this section, we apply our resampling framework to SMC [38], [48] and CSMC [2] methods [1]. Specifically, our framework gives simple proofs that SMC and CSMC algorithms can be implemented in a valid manner even if the resampling scheme (a) violates the (‘unbiasedness’, exchangeability or ‘marginal unbiasedness’) conditions commonly imposed in the literature, e.g., in [2], [26], [27]; (b) has a complicated dependence structure, e.g., as in systematic or chopthin resampling; or (c) is adaptive.

4.1 Sequential Monte Carlo↩︎

Let \(M_t(\,\cdot\,|x_{t-1})\) be conditional distributions on \(\mathsf{X}\) (‘mutation kernels’) and \(\smash{G_t\colon \mathsf{X}^2 \to (0, \infty)}\) (‘potential functions’). Our goal is to estimate expectations w.r.t.the distribution \(\pi_T(\mathrm{d}x_{1:T})\mathrel{\vcenter{:}}= \gamma_T(\mathrm{d}x_{1:T}) / \mathcal{Z}_T\), where \(\smash{\gamma_T(\mathrm{d}x_{1:T}) \mathrel{\vcenter{:}}= \prod_{t=1}^T G_t(x_{t-1:t}) M_t(\mathrm{d}x_t|x_{t-1})}\), with typically unknown normalising constant \(\smash{\mathcal{Z}_T} \mathrel{\vcenter{:}}= \smash{\int_{\mathsf{X}^T} \gamma_T(\mathrm{d}x_{1:T})}\); we use the convention that any quantity with subscript \(t < 1\) is to be ignored, i.e., \(M_1(\mathrm{d}x_1|x_0) = M_1(\mathrm{d}x_1)\) and \(G_1(x_{0:1}) = G_1(x_1)\).

The SMC algorithm is outlined in Algorithm 26. It uses \(N_t \in \mathbb{N}\) particles at time \(t\). For simplicity, we assume that \(N_t\) is deterministic but random numbers of particles could easily be accommodated by further extending the state space to explicitly include \(N_t\). Additionally, we let \(\lambda_1\) be some distribution on \([N_1]\) (usually, \(\lambda_1 \mathrel{\vcenter{:}}= \mathop{\mathrm{Unif}}_{[N_1]}\)). For any \(t > 1\), we let \(\langle \rho_{t-1}, \lambda_{t}\rangle\) be a properly weighted and persistent resampling scheme which draws \(\smash{N_t}\) ancestors based on \(N_{t-1}\) weights \(\mathbf{W}_{t-1} = \smash{W_{t-1}^{1:N_{t-1}}}\) that are defined below and which specifies post-resampling weights via \[\begin{align} \phi_t(m, n, u; \mathbf{W}_{t-1}) \mathrel{\vcenter{:}}= \dfrac{W_{t-1}^{n}\lambda_t(m, u|n; \mathbf{W}_{t-1})}{N_t^{-1}\rho_{t-1}^m(n, u; \mathbf{W}_{t-1})}. \label{eq:smc95post95resampling95weight95function} \end{align}\tag{18}\]

At the end of Algorithm 26, expectations \(\pi_T(f_T)\), integrals \(\gamma_T(f_T) = \mathcal{Z}_T \pi_T(f)\), for some suitably integrable function \(f_T \colon \mathsf{X}^T \to \mathbb{R}\), and the normalising constant \(\mathcal{Z}_T\) can be approximated as \[\begin{align} \hat{\pi}_T(f_T) \mathrel{\vcenter{:}}= \sum_{n=1}^{N_T} W_T^n f_T(X_{1:T}^{(n)}), \qquad \hat{\gamma}_T(f_T) \mathrel{\vcenter{:}}= \widehat{\mathcal{Z}}_T \hat{\pi}_T(f_T), \qquad \widehat{\mathcal{Z}}_T & \mathrel{\vcenter{:}}= \prod_{t=1}^T \sum_{n=1}^{N_t} w_t^n, \end{align}\] where \(\smash{X_{1:t}^{(n)}}\) is recursively defined as \(\smash{X_{1}^{(n)} \mathrel{\vcenter{:}}= X_1^n}\) and \(\smash{X_{1:t}^{(n)} \mathrel{\vcenter{:}}= (X_{1:t-1}^{(A_{t-1}^n)}, X_t^n)}\), for \(t > 1\).

The following result (proved in Appendix 9.3) generalises the unbiasedness of \(\alpha\)SMC from [22] which uses a framework for resampling that encompasses a number of the (adaptive) resampling schemes from Section 3 but not resampling schemes with a complex dependence structure like (adaptive) systematic resampling or chopthin resampling (Appendix 7.2 shows that \(\alpha\)SMC resampling is a strict special case of our framework).

The following simple example illustrates that proper weighting/unbiasedness of the resampling scheme is not sufficient to ensure convergence of estimators.

4.2 Conditional sequential Monte Carlo↩︎

Hereafter, assume that \(N_t \geq 2\), for any \(t \in [T]\). Then sampling \(\smash{X_{1:T}'}\) conditional on \(\smash{x_{1:T}}\) via the CSMC algorithm from Algorithm 27 induces a Markov kernel \(P\) on \(\smash{\mathsf{X}^T}\). The following result is proved in Appendix 9.4.

Proposition 8 generalises [32] which proved the validity of CSMC algorithms based around the \(\alpha\)SMC framework and thus encompasses a number of the resampling schemes from Section 3 but not resampling schemes with a complex dependence structure like (adaptive) systematic resampling or chopthin resampling nor the backward or ancestor-sampling extensions.

The following example illustrates that proper weighting/invariance of the (conditional) resampling scheme is not sufficient for ergodicity of the CSMC kernel.

4.3 A note on (over)simplification↩︎

In the case of multinomial resampling (and other exchangeable schemes), the CSMC algorithm remains valid if we fix \(K_t \mathrel{\vcenter{:}}= 1\), say, instead of sampling \(K_t\) from the index distribution, for \(t \in [T]\). In the case of multinomial resampling, fixing \(K_t \mathrel{\vcenter{:}}= 1\), can be implemented straightforwardly using one of the following two strategies whose outcomes are statistically equivalent:

1. Draw \(\smash{N_t}\) ancestor indices via standard multinomial resampling and then overwrite the first ancestor index with \(1\).

2. Set the first ancestor index to \(1\) and then draw the remaining ancestor indices via standard multinomial resampling (with \(\smash{N_t-1}\) instead of \(\smash{N_t}\) descendants).

These strategies are appealing because they allow for existing implementations of ‘unconditional’ resampling to be used for conditional resampling. One might thus hope that they are valid more generally.

Unfortunately, this is not the case. For instance, as we illustrate in Appendix 6, conditional systematic and residual-multinomial resampling implemented using the Strategies [enum:naive:i] or [enum:naive:ii] can induce a bias in the sense that the CSMC algorithm no longer leaves \(\smash{\pi_T}\) invariant. Further, these naïve implementations can also worsen the autocorrelation of the Markov chain induced by the CSMC algorithm.

5 Discussion↩︎

We have presented a simple framework for guaranteeing the validity of resampling schemes within SMC algorithms and their conditional analogues within CSMC algorithms under minimal assumptions, covering almost all standard resampling mechanisms as well as more exotic approaches such as chopthin resampling, whose conditional version we derive for the first time. Empirical results demonstrate that neglecting careful conditional resampling can lead to substantial bias. For simplicity, we have focused on basic SMC and CSMC algorithms. However, our arguments straightforwardly extend to models with non-Markovian dynamics, potential functions which depend on the entire history of the state sequence, resampling schemes with a random number of particles, and to the related algorithms from [51]–[56], but the details have been omitted here in the interests of brevity.

Acknowledgments↩︎

Large language models were used occasionally to assist with phrasing and with debugging and refining isolated sections of the accompanying code. All original ideas, mathematical derivations, and proofs are entirely the authors’ own work, and the manuscript was written predominantly without AI assistance. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission. Correspondence should be addressed to Axel Finke (axel.finke@ncl.ac.uk).

Funding↩︎

This work was supported by EPSRC grant no EP/K032208/1. AMJ acknowledges financial support from the Engineering and Physical Sciences Research Council (EPSRC; grants EP/R034710/1 and EP/T004134/1) and by United Kingdom Research and Innovation (UKRI) via grant EP/Y014650/1, as part of the ERC Synergy project OCEAN. Research of AL was supported by EPSRC grants EP/R034710/1 and EP/Y028783/1.

Data availability↩︎

No new data were generated or analysed during this study. The GitHub repository https://github.com/AxelFinke/resampling-in-conditional-smc-algorithms provides the code for reproducing the results referenced in Section 4.3 and shown in Appendix 6.

6 On naïve implementations of conditional resampling↩︎

In this section, we illustrate that the two naïve strategies for implementing conditional resampling discussed in Section 4.3 of the main manuscript do not lead to a valid CSMC algorithm. Specifically, we will illustrate that conditional systematic and conditional residual-multinomial resampling implemented using the Strategies [enum:naive:i]/[enum:naive:ii] (hereafter referred to as naïve systematic resampling [enum:naive:i]/[enum:naive:ii] and naïve residual-multinomial resampling [enum:naive:i]/[enum:naive:ii]) can induce a bias in the sense that the CSMC algorithm no longer leaves \(\smash{\pi_T}\) invariant. Further, these naïve implementations can also worsen the autocorrelation of the Markov chain induced by the CSMC algorithm.

6.1 Binary state-space model↩︎

We illustrate the issue on a toy binary state-space model for which it is possible to evaluate the invariant distribution and autocorrelation of the Markov chain induced by a CSMC algorithm for small values of \(N\) and \(T\). The model is as follows. For any \(t \in [T]\), \(\varepsilon, \delta \in [0,1]\) and any \(x_t \in \mathsf{X}\mathrel{\vcenter{:}}= \{0, 1\}\), let \(\smash{M_1} = \smash{(1-\varepsilon) \delta_0 + \varepsilon \delta_1}\), \(\smash{M_{t+1}(\,\cdot\,| x_t)} = \smash{(1-\varepsilon) \delta_{x_t} + \varepsilon \delta_{1 - x_t}}\) and \(\smash{G_t(x_{t-1:t})} = \smash{G_t(x_t)} = \smash{(1-\delta) \mathop{\mathrm{\mathbf{1}}}\{x_t = 0\} + \delta \mathop{\mathrm{\mathbf{1}}}\{x_t = 1\}}\). Additionally, we let \(\smash{N_1 = \dotsc = N_T = N}\) and \(\smash{\lambda_1 = \mathop{\mathrm{Unif}}_{[N]}}\).

As illustrated in Figure 28 (a), naïve conditional residual-multinomial resampling [enum:naive:i] suffers from substantial bias and increasing the number of particles or employing backward sampling reduces but does not remove this bias. Figure 28 (b) demonstrates that the other naïve implementations also induce bias (albeit of a much smaller magnitude) with two exceptions:

1. Whenever the underlying resampling scheme is unbiased and \(N = 2\), the naïve implementation of conditional resampling using Strategy [enum:naive:ii] remains exact (in the sense that the Markov chain leaves \(\smash{\pi_T}\) invariant) because it reduces to conditional multinomial resampling. We note that this holds generally, i.e., not just in this specific model.

2. Additionally, in this specific model, it can be shown that naïve conditional residual-multinomial resampling [enum:naive:ii] reduces to conditional multinomial resampling (and is thus exact) whenever \(\smash{\lvert \delta - \tfrac{1}{2}\rvert < 1 / (4N-6)}\).

As illustrated in Figure 29, the fact that naïve conditional resampling [enum:naive:ii] reduces to conditional multinomial resampling if \(N = 2\) (and can still behave similarly to conditional multinomial resampling if \(N > 2\)) can also negatively impact the autocorrelation of the Markov chain induced by a CSMC algorithm, especially when the potential functions oscillate little.

6.2 Linear-Gaussian state-space model↩︎

One might speculate that the extent of the bias induced by naïve conditional resampling [enum:naive:i] was somehow due to the specific finite-state-space structure of the above example. Here, we illustrate that a similar bias is incurred in a simple univariate linear-Gaussian state-space model, i.e., if \(M_1(x_1) = \mathop{\mathrm{N}}(x_1; 0, 1)\), \(M_t(x_t|x_{t-1}) = \mathop{\mathrm{N}}(x_t; x_{t-1}, 1)\), for \(t > 1\), and \(G_t(x_{t-1:t}) = G_t(x_t) = \mathop{\mathrm{N}}(y_t; x_t, 1)\), for given observations \(y_t \in \mathbb{R}\).

Figure 30 illustrates the bias induced in the approximation of the marginal posterior distribution at time \(t=2\) by naïve conditional resampling [enum:naive:i] in the case that \(T = 4\) and \(y_{1:4} = (0.75, 0.62, 1.91, 0.12)\) (the shown behaviour is typical for other observation sequences sampled from the model).

7 Relationship with other frameworks↩︎

7.1 Relationship with (conditional) resampling based on randomly permuting or shifting the order of the ancestor indices↩︎

7.1.1 Random permutations↩︎

The framework from [2] assumed, for simplicity, that the resampling schemes used within a (conditional) SMC algorithm are exchangeable. And although many popular resampling schemes (like systematic resampling) are not exchangeable, [2] noted that any resampling schemes can be made exchangeable by randomly permuting the order of the ancestor indices after they have been sampled. This idea was used to derive a conditional version of residual-multinomial resampling in [26] in which \(K\) is always set to \(1\). In this section, we show that unless an implementation of conditional resampling with \(K = 1\), say, is required, the use such random permutations is unnecessary (and only adds computational cost).

Formally, let \(\langle \varrho, \kappa\rangle\) be some balanced resampling scheme which may not be exchangeable. We can then derive an exchangeable resampling scheme \(\langle \rho, \lambda\rangle\) by setting \[\begin{align} \rho(\mathbf{a}, u) = \sum_{\mathbf{b} \in [N]^M} \varrho(\mathbf{b}, u) \sum_{\sigma \in \mathsf{S}_M} \mathop{\mathrm{Unif}}_{\mathsf{S}_M}(\sigma) \delta_{\mathbf{b}^\sigma}(\mathbf{a}), \label{eq:exchangeable95resampling95distribution95through95random95permutation} \end{align}\tag{19}\] where for any \(\mathbf{b} \mathrel{\vcenter{:}}= \smash{b^{1:M}}\) and any permutation \(\sigma \in \mathsf{S}_M \mathrel{\vcenter{:}}= \smash{\{\tau\colon [M] \to [M] \mid \text{\tau is bijective}\}}\), we write \(\smash{\mathbf{b}^\sigma} \mathrel{\vcenter{:}}= \smash{(b^{\sigma(1)}, \dotsc, b^{\sigma(M)})}\); as well as by setting \(\lambda(m, u|n) = \smash{\mathop{\mathrm{Unif}}_{[M]}(m)\rho^m(u|n)}\), for any \(n \in [N]\); then \[\begin{align} \phi(m, n, u) & = M W^n\frac{\lambda(m, u|n)}{\rho^m(n, u)}\\ & = W^n \frac{1}{\frac{1}{\# \mathsf{S}_M}\sum_{\sigma \in \mathsf{S}_M} \sum_{l=1}^M \mathop{\mathrm{\mathbf{1}}}\{\sigma(m) = l\} \varrho^l(n)} = \frac{M W^n}{\sum_{l=1}^M \varrho^l(n)}. \end{align}\] Thus, application of the random permutation does not change the law of the resampled measure \(\tilde{\pi}\).

Comparing Algorithm 31 with Algorithm 1, we see that randomly permuting the the order of the ancestor indices only adds computational cost Specifically, the former requires the additional Line [alg:resampling95with95random95permutation:extra95step].

Sampling the index \(K\) within the conditional version of \(\langle \rho, \lambda\rangle\) becomes simpler compared to \(\langle \varrho, \kappa\rangle\) since the distribution of \(K\) is now uniform under \(\lambda(\,\cdot\,|a)\). However, conditional resampling also requires sampling from \(\smash{\rho^k(u|a^k) \rho^{-k}(\mathbf{a}^{-k}|a^k, u) \delta_a(a^k)}\) and the random permutation makes the latter more involved: \[\begin{align} \MoveEqLeft \rho^k(u|a^k) \rho^{-k}(\mathbf{a}^{-k}|a^k, u) \delta_a(a^k)\\ & = \frac{\rho(\mathbf{a}, u) \mathop{\mathrm{\mathbf{1}}}\{a^k = a\}}{\rho^k(a)}\\ & = \frac{\sum_{\mathbf{b} \in [N]^M} \varrho(\mathbf{b}, u) \sum_{\sigma \in \mathsf{S}_M} \mathop{\mathrm{Unif}}_{\mathsf{S}_M}(\sigma) \delta_{\mathbf{b}^\sigma}(\mathbf{a})\mathop{\mathrm{\mathbf{1}}}\{a^k = a\}}{\frac{1}{M}\sum_{m=1}^M \varrho^m(a)} \\ & = \sum_{\mathbf{b} \in [N]^M} \sum_{l=1}^M \frac{\varrho^l(a)}{\sum_{m=1}^M \varrho^m(a)} \varrho^l(u|a) \delta_{a}(b^l) \varrho^{-l}(\mathbf{b}^{-l} |b^l, u) \smashoperator{\sum_{\sigma \in \mathsf{S}_M^k(l)}} \mathop{\mathrm{Unif}}_{\mathsf{S}_M^k(l)}(\sigma) \delta_{\mathbf{b}^\sigma}(\mathbf{a})\\ & = \sum_{\mathbf{b} \in [N]^M} \sum_{l=1}^M \kappa(l, u|a) \delta_{a}(b^l) \varrho^{-l}(\mathbf{b}^{-l} | b^l, u) \smashoperator{\sum_{\sigma \in \mathsf{S}_M^k(l)}} \mathop{\mathrm{Unif}}_{\mathsf{S}_M^k(l)}(\sigma) \delta_{\mathbf{b}^\sigma}(\mathbf{a}), \end{align}\] where \(\smash{\mathsf{S}_M^k(l) \mathrel{\vcenter{:}}= \{\sigma \in \mathsf{S}_M \mid \sigma(k) = l\}}\) which, in particular, means that in the last two lines, \(\smash{\sigma \in \mathsf{S}_M^k(l)}\) implies that \(\smash{b^l} = \smash{b^{\sigma(k)} = a^{\sigma^{-1}(l)}} = \smash{a^k}\).

Again by comparing Algorithm 32 with Algorithm 2, we see that the random permutation incurs the additional Line [alg:conditional95resampling95with95random95permutation:extra95step].

In summary, our work shows that randomly permuting the order of the ancestor indices is not needed (and only adds unnecessary computational cost) unless one wishes to exploit Proposition 5 to implement conditional resampling in a way that always sets \(K \mathrel{\vcenter{:}}= 1\) or \(K \mathrel{\vcenter{:}}= M\), say, in Line [alg:conditional95resampling95with95random95permutation:extra95step].

7.1.2 Random cyclical shifts↩︎

All the developments from the previous section remain valid if we replace the random permutations \(\sigma \in \mathsf{S}_M\) by the subset of \(\mathsf{S}_M\) which correspond to random cyclical shifts \(\sigma \in \mathsf{C}_M \mathrel{\vcenter{:}}= \{\tau\colon [M] \to [M] \mid \exists\, s \in [M]\; \forall\, m \in [M]: \tau(m) = m + s \mod M\}\) (and if we similarly replace \(\mathsf{S}_M^k(l)\) by \(\mathsf{C}_M^k(l) \mathrel{\vcenter{:}}= \{\sigma \in \mathsf{C}_M \mid \sigma(k) = l\}\)). In particular, note that the resampling scheme \(\langle \rho, \lambda\rangle\) induced by applying such random cyclical shifts to the order of the ancestor indices generated by \(\varrho\) still satisfies the assumption of Proposition 5.

Such random cyclical shifts were employed in [26] to derive a conditional version of systematic resampling in which \(K\) can always be set to \(1\). Specifically, [26] can be derived from Algorithm 34 by taking \(\langle \varrho, \kappa\rangle\) to be the systematic resampling scheme and \(a = 1\) as follows:

• in Line [alg:conditional95resampling95with95random95shift:1] of Algorithm 34, sample \[\begin{align} U\sim \frac{\lfloor MW^1\rfloor}{MW^1}\mathop{\mathrm{Unif}}_{[0,1]}(u) + \biggl(1 - \frac{\lfloor MW^1\rfloor}{MW^1}\biggr) \mathop{\mathrm{Unif}}_{(0, MW^1 - \lfloor MW^1\rfloor]}(u), \label{eq:conditional95systematic95resampling95with95random95shift95conditional95distribution95of95u95given95a} \end{align}\tag{20}\] which is justified by the fact that \(\kappa(u|a)\) in \[\begin{align} \kappa(l, u|a) & = \kappa(u|a)\kappa(l|u,a)\\ & = \biggl[\sum_{m=1}^M \frac{\mathop{\mathrm{vol}}(I^m(a))}{MW^a} \mathop{\mathrm{Unif}}_{I^m(a)}(u)\biggr] \mathop{\mathrm{Unif}}_{\{m \in [M] \mid u \in I^m(a)\}}(l), \end{align}\] simplifies to the r.h.s.in 20 due to \(a = 1\);

• also in Line [alg:conditional95resampling95with95random95shift:1] of Algorithm 34, set \(B^m \mathrel{\vcenter{:}}= F^{-1}([U + m - 1]/M)\), for \(m \in [M]\), using Algorithm 7 with \(U^1 = \ldots = U^M = U\) and simply set \(L \mathrel{\vcenter{:}}= 1\) both of which are justified by the fact that \(a = 1\) and that given \(u \sim \kappa(u|a)\), \[\begin{align} \delta_{a}(b^l) \varrho^{-l}(\mathbf{b}^{-l}|b^l, u) & = \prod_{m=1}^M \delta_{F^{-1}([u + m - 1] / M)}(b^m) \end{align}\] does not depend on \(l\);

• in Line [alg:conditional95resampling95with95random95shift:extra95step] of Algorithm 34, appeal to Proposition 5 to set \(K \mathrel{\vcenter{:}}= 1\).

Furthermore, [27] introduced conditional versions of killing resampling and of systematic resampling with weight re-ordering based on such random cyclical shifts which could likewise be justified as instances of Algorithm 34 in which \(\langle \varrho, \kappa\rangle\) is killing resampling or systematic resampling with weight re-ordering, respectively.

Our work shows that such random cyclical shifts are only necessary if one wishes to exploit Proposition 5 to implement conditional resampling in a way that always sets \(K \mathrel{\vcenter{:}}= 1\) or \(K \mathrel{\vcenter{:}}= M\), say, in Line [alg:conditional95resampling95with95random95shift:extra95step].

7.2 Relationship with \(\alpha\)SMC resampling/resamping matrices↩︎

The \(\alpha\)SMC framework from [22] resampling provides a matrix representation of resampling schemes for which \(\mathsf{U}\mathrel{\vcenter{:}}= \emptyset\) and the resampling distribution factorises as \[\begin{align} \rho(\mathbf{a}) = \prod_{m=1}^M \rho^m(a^m). \label{eq:alpha95smc95factorisation} \end{align}\tag{21}\] Let \(\alpha = (\alpha^{m,n}) \in [0, \infty)^{M \times N}\) be some \(\mathcal{F}\)-measurable resampling matrix which is such that each row defines a probability measure on \([N]\), i.e., such that \(\sum_{n=1}^N \alpha^{m,n} = 1\), for any \(m \in [M]\). We will also assume that \[\begin{align} \sum_{m=1}^M \alpha^{m,n} > 0, \label{eq:alpha95smc95support95condition} \end{align}\tag{22}\] for any \(n \in [N]\) with \(W^n > 0\).

Then \(\alpha\)SMC represents \(\rho^m(n)\) as \[\begin{align} \rho^m(n) \mathrel{\vcenter{:}}= \mathop{\mathrm{Cat}}\biggl(n; \biggl(\frac{\alpha^{m,k} W^k}{\sum_{l=1}^N \alpha^{m,l} W^l}\biggr)_{k = 1}^N \biggr), \end{align}\] which is well defined by 22 (and the usual convention \(\tfrac{0}{0} = 0\)).

Use of such resampling matrices goes back to at least [57]

To embed \(\alpha\)SMC resampling within our framework, set \[\begin{align} \lambda(m|n) \mathrel{\vcenter{:}}= \frac{\alpha^{m,n}}{\sum_{l=1}^M \alpha^{l,n}}, \end{align}\] whenever \(W^n > 0\), so that \[\begin{align} \phi(m,n) & = M W^n \frac{\lambda(m|n)}{\rho^m(n)} = M \sum_{l=1}^N \frac{\alpha^{m,l}}{\sum_{k=1}^M \alpha^{k,l}} W^l. \end{align}\] Then the \(\alpha\)SMC resampling scheme \(\langle \rho, \lambda\rangle\) is not necessarily balanced but still properly weighted and persistent.

The factorisation 21 holds for most of the elementary resampling schemes from Section 3.1, including trivial resampling, i.e., not resampling (recall that as explained in Remark 3, whether or not we explicitly include the auxiliary variable \(U\) into the space is immaterial for balanced resampling schemes). And since the resampling matrix \(\alpha\) is allowed to be some (measurable) function of the particle system, \(\alpha\)SMC even certain kinds of ESS-adaptive resampling. Additionally, a conditional version of \(\alpha\)SMC resampling was derived in [32].

However, the \(\alpha\)SMC framework does not cover resampling distributions which violate the factorisation from 21 , such as the widely used (adaptive) systematic resampling scheme or more ‘exotic’ schemes like chopthin resampling. This shows that our framework provides a strict generalisation of \(\alpha\)SMC.

8 Proofs for Section 2↩︎

9 Proofs for Section 4↩︎

9.1 CSMC with backward or ancestor sampling↩︎

We again assume that \(N_t \geq 2\) and make the following assumption on the mutation kernels to ensure that the transition densities \(\smash{m_t(\,\cdot\,|x) \mathrel{\vcenter{:}}= \mathrm{d}M_t(\,\cdot\,|x) / \mathrm{d}\nu}\) are well defined for any \(x \in \mathsf{X}\) and \(t \in [T]\).

Assumption 1 guarantees the existence of the following backward kernel, for any \(k \in [N_t]\), any \(t < T\) and with the notation \(\smash{\mathbf{x}_t \mathrel{\vcenter{:}}= x^{1:N_t}}\): \[\begin{align} B_t(k|\mathbf{x}_t, x_{t+1}; \mathbf{W}_t) & \mathrel{\vcenter{:}}= \dfrac{W_t^k G_{t+1}(x_t^k, x_{t+1}) m_{t+1}(x_{t+1}|x_t^k)}{\sum_{n=1}^{N_t} W_t^n G_{t+1}(x_t^n, x_{t+1}) m_{t+1}(x_{t+1}|x_t^n)}. \label{eq:backward95kernel} \end{align}\tag{23}\]

The following result—proved in Appendix 9.5—extends Proposition 8.

9.2 Extended state-space construction↩︎

Throughout this section, we write \(\smash{\mathbf{X}_t \mathrel{\vcenter{:}}= X_t^{1:N_t}}\) as well as \(\mathbf{Z}_t \mathrel{\vcenter{:}}= (\mathbf{X}_t, \mathbf{A}_{t-1}, U_{t-1})\), for \(t > 1\) and \(\mathbf{Z}_1 \mathrel{\vcenter{:}}= \mathbf{X}_1\). The weights \(\smash{\mathbf{W}_{t} = W_{t}^{1:N_{t}}}\) then depend on \(\mathbf{Z}_{1:T}\) although we do not make this explicit to simplify the notation. The following distributions are well defined if the resampling schemes \(\langle \rho_{t-1}, \lambda_{t}\rangle\) are all properly weighted and persistent.

1. The law of all the random variables generated by the SMC algorithm (Algorithm 26) is then \[\begin{align} \mathrm{smc}_T(\mathrm{d}\mathbf{z}_{1:T}) & \mathrel{\vcenter{:}}= \Biggl[ \prod_{n=1}^{N_1} M_1(\mathrm{d}x_1^n)\biggr] \prod_{t=2}^T \biggl[\rho_{t-1}(\mathrm{d}[\mathbf{a}_{t-1}, u_{t-1}]; \mathbf{W}_{t-1}) \prod_{n=1}^{N_t} M_t(\mathrm{d}x_t^n|x_{t-1}^{a_{t-1}^n})\biggr]. \end{align}\]

2. Similarly, the law of all the random variables generated by the CSMC algorithm with backward sampling (Algorithm 35)—under Assumption 1—is \[\begin{align} \mathrm{csmc}_T(\mathrm{d}[k_{1:T}, \mathbf{z}_{1:T}]|x_{1:T}) \mathrm{back}_T(\mathrm{d}[l_{1:T}, x_{1:T}']| \mathbf{z}_{1:T}), \end{align}\] where the first term is the law of all the random variables generated by Lines [alg:csmc95with95backward95sampling:1a]–[alg:csmc95with95backward95sampling:1b]: \[\begin{align} \MoveEqLeft\mathrm{csmc}_T(\mathrm{d}[k_{1:T}, \mathbf{z}_{1:T}]|x_{1:T})\\ & \mathrel{\vcenter{:}}= \Biggl[ \lambda_1(\mathrm{d}k_1) \delta_{x_1}(\mathrm{d}x_1^{k_1}) \prod_{\substack{\mathllap{n} = \mathrlap{1}\\ \mathllap{n} \neq \mathrlap{k_1}}}^{N_1} M_1(\mathrm{d}x_1^n)\Biggr] \prod_{t=2}^T \Biggl[ \lambda_t(\mathrm{d}[k_t, u_{t-1}] | k_{t-1}, \mathbf{W}_{t-1})\\[-2ex] & \qquad \times \delta_{k_{t-1}}(\mathrm{d}a_{t-1}^{k_t}) \rho_{t-1}^{-k_t}(\mathrm{d}\mathbf{a}_{t-1}^{-k_t} | a_{t-1}^{k_t}, u_{t-1}; \mathbf{W}_{t-1}) \delta_{x_t}(\mathrm{d}x_t^{k_t}) \prod_{\substack{\mathllap{n} = \mathrlap{1}\\ \mathllap{n} \neq \mathrlap{k_t}}}^{N_t} M_t(\mathrm{d}x_t^n|x_{t-1}^{a_{t-1}^n}) \Biggr]; \end{align}\] and the second term is the law of all the random variables generated in Lines [alg:csmc95with95backward95sampling:2]–[alg:csmc95with95backward95sampling:4]: \[\begin{align} \mathrm{back}_T(\mathrm{d}[l_{1:T}, x_{1:T}']| \mathbf{z}_{1:T}) & \mathrel{\vcenter{:}}= W_T^{l_T} \biggl[ \prod_{t=1}^{T-1} B_t(\mathrm{d}l_t|\mathbf{x}_t, x_{t+1}^{l_{t+1}}; \mathbf{W}_t) \biggr] \prod_{t=1}^T \delta_{x_t^{l_t}}(\mathrm{d}x_t'). \end{align}\]

3. Finally, the law of all the random variables generated by the CSMC algorithm with ancestor sampling (Algorithm 36)—again under Assumption 1— is \[\begin{align} \mathrm{csmc}^*_T(\mathrm{d}[k_{1:T}, \mathbf{z}_{1:T}]|x_{1:T}) \mathrm{anc}_T(\mathrm{d}[l_{1:T}, x_{1:T}']| \mathbf{z}_{1:T}), \end{align}\] where the first term is the law of all the random variables generated by Lines [alg:csmc95with95ancestor95sampling:1a]–[alg:csmc95with95ancestor95sampling:1b]: \[\begin{align} \MoveEqLeft\mathrm{csmc}^*_T(\mathrm{d}[k_{1:T}, \mathbf{z}_{1:T}]|x_{1:T})\\ & \mathrel{\vcenter{:}}= \Biggl[ \lambda_1(\mathrm{d}k_1) \delta_{x_1}(\mathrm{d}x_1^{k_1}) \prod_{\substack{\mathllap{n} = \mathrlap{1}\\ \mathllap{n} \neq \mathrlap{k_1}}}^{N_1} M_1(\mathrm{d}x_1^n)\Biggr]\\[-2ex] & \quad \times \prod_{t=2}^T \Biggl[ \lambda_t(\mathrm{d}[k_t, u_{t-1}] | k_{t-1}, \mathbf{W}_{t-1}) B_{t-1}(\mathrm{d}a_{t-1}^{k_t}| \mathbf{x}_{t-1}, x_t; \mathbf{W}_{t-1})\\[-3ex] & \qquad\qquad \times \rho_{t-1}^{-k_t}(\mathrm{d}\mathbf{a}_{t-1}^{-k_t} | a_{t-1}^{k_t}, u_{t-1}; \mathbf{W}_{t-1}) \delta_{x_t}(\mathrm{d}x_t^{k_t}) \prod_{\substack{\mathllap{n} = \mathrlap{1}\\ \mathllap{n} \neq \mathrlap{k_t}}}^{N_t} M_t(\mathrm{d}x_t^n|x_{t-1}^{a_{t-1}^n}) \Biggr]; \end{align}\] and the second term is the law of all the random variables generated in Lines [alg:csmc95with95ancestor95sampling:2]–[alg:csmc95with95ancestor95sampling:4]: \[\begin{align} \mathrm{anc}_T(\mathrm{d}[l_{1:T}, x_{1:T}']|\mathbf{z}_{1:T}) \mathrel{\vcenter{:}}= W_T^{l_T} \biggl[ \prod_{t=1}^{T-1} \delta_{a_t^{l_{t+1}}}(\mathrm{d}l_t) \biggr] \prod_{t=1}^T \delta_{x_t^{l_t}}(\mathrm{d}x_t'). \end{align}\]

In the following, we use the convention that \(\mu(\mathrm{d}x) P(\mathrm{d}y| x) = [\mu \otimes P](\mathrm{d}[x, y])\) and with some abuse of notation regarding the order of the arguments. We also make the dependence of \(\smash{\widehat{\mathcal{Z}}_T = \widehat{\mathcal{Z}}_T(\mathbf{Z}_{1:T})}\) on \(\mathbf{Z}_{1:T}\) explicit in the notation.

9.3 Proof of Proposition 7↩︎

The proof slightly extends the ‘SMC as importance sampling’ argument from [2] to our setting.

By Part [lem:radon--nikodym95derivative:1] of Lemma 1, we have \(\gamma_T \otimes \mathrm{csmc}_T \ll \mathrm{smc}_T \otimes \mathrm{anc}_T\). The proof then follows by a change of measure and the tower property of conditional expectation: \[\begin{align} \gamma_T(f) & = \int f(x_{1:T}) [\gamma_T \otimes \mathrm{csmc}_T](\mathrm{d}[\mathbf{z}_{1:T}, k_{1:T}, x_{1:T}])\\ & = \mathop{\mathrm{\mathbb{E}}}\biggl[f(X_{1:T}) \dfrac{\mathrm{d}[\gamma_T \otimes \mathrm{csmc}_T]}{\mathrm{d}[\mathrm{smc}_T \otimes \mathrm{anc}_T]}(\mathbf{Z}_{1:T}, K_{1:T}, X_{1:T})\biggr]\\ & = \mathop{\mathrm{\mathbb{E}}}\biggl[\mathop{\mathrm{\mathbb{E}}}\biggl(f(X_{1:T}) \dfrac{\mathrm{d}[\gamma_T \otimes \mathrm{csmc}_T]}{\mathrm{d}[\mathrm{smc}_T \otimes \mathrm{anc}_T]}(\mathbf{Z}_{1:T}, K_{1:T}, X_{1:T})\bigg|\mathbf{Z}_{1:T}\biggr)\biggr]\\ & = \mathop{\mathrm{\mathbb{E}}}\biggl[\widehat{\mathcal{Z}}_T(\mathbf{Z}_{1:T}) \sum_{n=1}^{N_T} W_T^n f(X_{1:T}^{(n)}) \biggr]\\ & = \mathop{\mathrm{\mathbb{E}}}[ \hat{\gamma}_T(f)], \end{align}\] where the expectations in the first three lines are w.r.t.\((\mathbf{Z}_{1:T}, K_{1:T}, X_{1:T}) \sim \mathrm{smc}_T \otimes \mathrm{anc}_T\) and the expectations in the last two lines are w.r.t.\(\smash{\mathbf{Z}_{1:T} \sim \mathrm{smc}_T}\). \(\Box\)

9.4 Proof of Proposition 8↩︎

The first part of Lemma 1, gives \[\begin{align} \MoveEqLeft \pi_T(\mathrm{d}x_{1:T}) \mathrm{csmc}_T(\mathrm{d}[k_{1:T}, \mathbf{z}_{1:T}]|x_{1:T}) \mathrm{anc}_T(\mathrm{d}[l_{1:T}, x_{1:T}']| \mathbf{z}_{1:T})\\ & = \widehat{\mathcal{Z}}_T(\mathbf{z}_{1:T}) \mathrm{smc}_T(\mathrm{d}\mathbf{z}_{1:T}) \mathrm{anc}_T(\mathrm{d}[k_{1:T}, x_{1:T}]|\mathbf{z}_{1:T}) \mathrm{anc}_T(\mathrm{d}[l_{1:T}, x_{1:T}']|\mathbf{z}_{1:T})\\ & = \pi_T(\mathrm{d}x_{1:T}') \mathrm{csmc}_T(\mathrm{d}[l_{1:T}, \mathbf{z}_{1:T}]|x_{1:T}') \mathrm{anc}_T(\mathrm{d}[k_{1:T}, x_{1:T}]| \mathbf{z}_{1:T}). \end{align}\] This completes the proof. \(\Box\)

9.5 Proof of Proposition 9↩︎

The proof could be established in a number of ways, e.g., using partially-collapsed Gibbs sampling arguments [49], [50]. Here, we use the proof technique from [58] which is convenient in the sense that it simultaneously guarantees the validity of backward and ancestor sampling.

Specifically, appealing to both parts of Lemma 1, \[\begin{align} \MoveEqLeft \pi_T(\mathrm{d}x_{1:T}) \mathrm{csmc}_T(\mathrm{d}[k_{1:T}, \mathbf{z}_{1:T}]|x_{1:T}) \mathrm{back}_T(\mathrm{d}[l_{1:T}, x_{1:T}']| \mathbf{z}_{1:T})\\ & = \widehat{\mathcal{Z}}_T(\mathbf{z}_{1:T}) \mathrm{smc}_T(\mathrm{d}\mathbf{z}_{1:T}) \mathrm{anc}_T(\mathrm{d}[k_{1:T}, x_{1:T}]|\mathbf{z}_{1:T}) \mathrm{back}_T(\mathrm{d}[l_{1:T}, x_{1:T}']|\mathbf{z}_{1:T})\\ & = \pi_T(\mathrm{d}x_{1:T}') \mathrm{csmc}^*_T(\mathrm{d}[l_{1:T}, \mathbf{z}_{1:T}]|x_{1:T}') \mathrm{anc}_T(\mathrm{d}[k_{1:T}, x_{1:T}]| \mathbf{z}_{1:T}). \end{align}\] This completes the proof in the case of backward or ancestor sampling. \(\Box\)

References↩︎

---

1. School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne, UK. axel.finke@ncl.ac.uk. ORCID: 0000-0002-8379-3012↩︎

2. Department of Statistics, University of Warwick, Coventry, UK. a.m.johansen@warwick.ac.uk. ORCID: 0000-0002-3531-7628↩︎

3. School of Mathematics, University of Bristol, Bristol, UK. anthony.lee@bristol.ac.uk. ORCID: 0000-0001-7765-0616↩︎

4. Independent researcher. lawrence@indii.org. ORCID: 0000-0002-6567-6015↩︎