Parallel resampling in the particle filter

Introduction

The particle filter, and more generally Sequential Monte Carlo (SMC) methods, constitute a large class of numerical methods routinely used to perform statistical inference. Particle filters were originally developed for object tracking and time series analysis using nonlinear, non-Gaussian state-space models (Gordon et al., 1993; Doucet et al., 2001). They have been extended to accommodate general statistical models (Chopin, 2002; Del Moral, 2004; Del Moral et al., 2006), with recent applications including rare event estimation (Cérou et al., 2012), graphical models (Naesseth et al., 2014), phylogenetic inference (Bouchard-Côté et al., 2012) and variable selection (Schäfer and Chopin, 2013). SMC has shown comparative advantage over Markov chain Monte Carlo (MCMC) when the target distribution is multimodal (Chopin and Jacob, 2010; Schweizer, 2012) or when the interest lies in the normalizing constant of the target distribution (Zhou et al., 2013).

The general framework of SMC involves introducing a sequence of distributions $\pi_{0},\ldots,\pi_{T}$, where the interest might be in each distribution $\pi_{t}$ or only in the last one $\pi_{T}$. At step $t=0$, particles $\mathbf{x}_{0}^{1:N}\equiv\mathbf{x}_{0}^{1},\ldots,\mathbf{x}_{0}^{N}$ are drawn independently from $\pi_{0}$, and each weight $w_{0}^{i}$ in the vector $w_{0}^{1:N}$ is set to $1/N$. The weighted particles $(w_{0}^{1:N},\mathbf{x}_{0}^{1:N})$ constitute an empirical approximation of $\pi_{0}$. Then at any step $t\geq 1$ of the algorithm, the previous particles $(w_{t-1}^{1:N},\mathbf{x}_{t-1}^{1:N})$ approximating $\pi_{t-1}$ are propagated and weighted to obtain new particles $(w_{t}^{1:N},\mathbf{x}_{t}^{1:N})$, which approximate $\pi_{t}$. A generic way to achieve this (Del Moral, 2004) involves sequences of Markov kernels $K_{t}$ and potential functions $G_{t}$ taking values in $(0,+\infty)$, as in the algorithm in Code 1.

Before giving more details on the resampling step, let us describe two choices for $K_{t}$ and $G_{t}$. Consider, for example, a state-space model where $\mathbf{x}_{0:T}$ is an unobserved Markov chain and $\mathbf{y}_{1:T}$ are conditionally independent observations of $\mathbf{x}_{0:T}$ with additional noise, such that the joint distribution can be written

For a given data set $\mathbf{y}_{1:T}$, the interest is to draw samples from the filtering distributions $\pi_{t}(\mathbf{x}_{t})=p(\mathbf{x}_{t}\,|\,\mathbf{y}_{1:t})$ for $t=1,\ldots,T$. When the probability densities $p(\mathbf{x}_{t}\,|\,\mathbf{x}_{t-1})$ and $p(\mathbf{y}_{t}\,|\,\mathbf{x}_{t})$ are linear and Gaussian, the Kalman filter (Kalman, 1960) can be used for this purpose. When they are nonlinear and non-Gaussian, the particle filter is preferred, as it provides asymptotically consistent estimates of quantities of interest as $N\rightarrow\infty$. The bootstrap particle filter (Gordon et al., 1993) corresponds to the choice $K_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=p(\mathbf{x}_{t}|\mathbf{x}_{t-1})$ and $G_{t}(\mathbf{x}_{t})=p(\mathbf{y}_{t}|\mathbf{x}_{t})$.

Another example is that of parameter inference, where the interest is in the posterior distribution $\pi(\boldsymbol{\theta})=p(\boldsymbol{\theta}|\mathbf{y}_{1:T})$ of a parameter $\boldsymbol{\theta}$ given a data set $\mathbf{y}_{1:T}$. Let $\pi_{0}(\boldsymbol{\theta})=p(\boldsymbol{\theta})$ (the prior distribution) and introduce, for all $t=1,\ldots,T$, $\pi_{t}(\boldsymbol{\theta})\propto p(\boldsymbol{\theta})p(\mathbf{y}_{1:t}|\boldsymbol{\theta})$. In this case, a practical choice (Chopin, 2002) is to choose $K_{t}$ to be an MCMC kernel leaving $\pi_{t-1}$ invariant, such as a Metropolis–Hastings kernel, and to define $G_{t}(\boldsymbol{\theta}):=p(\mathbf{y}_{t}|\mathbf{y}_{1:t-1},\boldsymbol{\theta})$, or, in the case of independent observations, $G_{t}(\boldsymbol{\theta}):=p(\mathbf{y}_{t}|\boldsymbol{\theta})$. Since any MCMC kernel can be chosen for $K_{t}$, SMC can be seen as a framework to turn an arbitrary MCMC algorithm into a population-based, and thus parallelisable, algorithm for parameter inference.

In Code 1, the resampling step is encoded by a randomised algorithm Ancestors that accepts a vector $\mathbf{w}_{t-1}\in[0,\infty)^{N}$ of weights, and returns a vector $\mathbf{a}_{t}\in\{1,\ldots,N\}^{N}$, where each $a^{i}_{t}$ is the index of the particle at time $t-1$ which is to be the ancestor of the $i$th particle at time $t$. Alternatively, the resampling step may be encoded by a randomised algorithm Offspring that also accepts a vector $\mathbf{w}_{t-1}\in[0,\infty)^{N}$ of particle weights, but instead returns a vector $\mathbf{o}_{t}\in\{0,\ldots,N\}^{N}$, where each $o^{i}_{t}$ is the number of offspring to be created from the $i$th particle at time $t-1$ for propagation to time $t$. Ancestry vectors are readily converted to offspring vectors and vice-versa (Appendix D provides functions to achieve this).

There are numerous acceptable algorithms for the resampling step. Recalling that the output is random, typically it is required only that, $\forall i\in\{1,\ldots,N\}$:

that is, the expected number of offspring of a particle should be equal to $N$ times its normalised weight. This unbiasedness condition ensures unbiased estimates of quantities such as the marginal likelihood $p(\mathbf{y}_{1:T})$, which follows from a simple extension of Proposition 7.4.1 in Del Moral (2004). Details of standard unbiased strategies, including their specific implementation in this work, are given in Appendix B. One common approach is to draw $\mathbf{o}_{t}$ according to a multinomial distribution with parameters $N$ and $\mathbf{w}_{t-1}$ (multinomial resampling). We will consider alternative strategies below, including some opportunities to relax the unbiasedness condition in exchange for a significant reduction in execution time. Figure 1 visualises both standard and alternative approaches.

2 Parallelisation

The initialisation, propagation and weighting steps of SMC are readily parallelised, being independent operations on each particle $\mathbf{x}_{t}^{i}$ and its weight $w_{t}^{i}$. Resampling, on the other hand, is a collective operation across particles and weights, so that parallelisation is more difficult. Summing the $N$ weights necessitates synchronization across threads, so that some threads may wastefully idle while waiting on others. For certain hardware that does not permit global communication between concurrently running threads, such as graphics processing units (GPUs), the elimination of collective operations can yield significant speed up. For this reason the resampling step has attracted recent attention as a potential bottleneck in the further scaling of SMC to larger systems (Murray, 2011; Whiteley et al., 2013).

At the core of most standard resampling schemes, such as the multinomial, stratified and systematic schemes, is a cumulative sum of weights, also called a prefix sum. A major theme of prior contributions has been the parallelisation of this prefix sum (Maskell et al., 2006; Hendeby et al., 2010; Chao et al., 2010; Gong et al., 2012). This is a generic problem that is also relevant to other algorithms (see e.g. Harris et al., 2007). Another theme in prior work is the partitioning of particles into disjoint subsets within which local resampling is performed (Chao et al., 2010). This is more useful in a distributed memory context, as it limits communication between processes (Whiteley et al., 2013); it has been considered in this context before (Brun et al., 2002; Bolić et al., 2005). General comments regarding the parallelisability of particle filtering algorithms are given in Lee et al. (2010) and Murray (2013).

3 Numerical precision and stability

As demonstrated later in this work, numerical instabilities can be apparent in standard resampling schemes when $N$ is large. One million particles is not an unrealistic number for some contemporary applications of SMC (see e.g. Klaas et al., 2006; Kitagawa, 2014). In double-precision, where 15 significant figures (in decimal) are expected, it is unlikely that $N$ will be sufficiently large for numerical instability to be a problem in current applications. In single-precision, however, where 8 significant figures (in decimal) are expected, numerical instabilities can manifest with this many particles. This is important because contemporary hardware has significantly faster single-precision than double-precision floating-point performance. There is reason to believe that this gap will remain; consider, for example, that single-precision is twice as fast on CPU architectures—even reasonably mature architectures—when SIMD instructions (such as those of SSE and AVX) are employed. On current GPUs, while the architecture is changing more rapidly, using single-precision also leads to at least a two-fold speedup. On other architectures, such as field-programmable gate arrays (FPGAs), custom precision is possible, allowing a trade-off between accuracy and performance (Mingas and Bouganis, 2012).

4 Contributions

As identified above, standard resampling algorithms require a cumulative sum over weights, leading to two problems:

they exhibit numerical instability for large numbers of particles or large weight variance.

This work contributes two alternative algorithms that eliminate the cumulative sum over weights in order to remedy both problems. They are better suited to the breadth of parallelism afforded by modern hardware, and do not exhibit the numerical instability of standard schemes. The two alternative schemes are based on Metropolis and rejection samplers. In comparing these to standard schemes, we carefully consider the consequences of the choice of resampling scheme on both the CPU and GPU, considering how numerical stability, bias, mean squared error and execution time vary across the number of particles and the variability in their weights. In this respect, the work constitutes a thorough study of resampling schemes and a useful guide to the selection and implementation of the most appropriate algorithm for a given problem.

Section 2 describes the Metropolis and rejection resampling schemes. We also describe a useful permutation of ancestry vectors in Section 2.3 to prevent read and write conflicts between concurrently running threads. Empirical comparisons around bias, mean squared error and execution time are given in Section 3, with concluding remarks in Section 4.

Appendix A provides our pseudocode conventions for reference. Appendix B recalls the standard algorithms for resampling based on multinomial, stratified and systematic sampling, highlighting their use of collective operations. Appendix C gives more details on the permutation algorithm. Appendix D presents auxiliary functions for converting between offspring and ancestry vectors. Appendix E provides some implementation notes.

Alternative resampling schemes

The first approach resamples via the Metropolis algorithm (Metropolis et al., 1953) rather than direct sampling, giving a result close to that of the multinomial resampler. The approach was briefly studied by the first author in a technical report (Murray, 2011), but a more complete treatment and improved analysis is given here. Instead of the collective operation, only the ratio between pairs of weights is ever computed. Code 2 describes the approach.

The Metropolis resampler is parameterised by $B$, the number of iterations to be performed before convergence is assumed and each particle may settle on its chosen ancestor. We can view the inner for loop as iterating a Markov kernel $P$ with stationary distribution $\pi$ for $B$ steps, where

is the categorical distribution associated with the weights.

As $B$ must be finite, the algorithm produces a biased sample. Setting $B$ is a tradeoff between speed and bias, with smaller $B$ giving faster execution time but larger bias. This bias may not be much of a problem for filtering applications, but does violate the assumptions that lead to unbiased marginal likelihood estimates in a particle MCMC framework (Andrieu et al., 2010), so care should be taken.

To provide guidance as to the selection of $B$, we bound the total variation distance of the $B$-fold iterate $P^{B}(i,\cdot)$ from $\pi$, where

Such a bound can be obtained by noting that $P$ is an independent Metropolis Markov kernel with target $\pi$ and a uniform proposal on $\{1,\ldots,N\}$. By Theorem 2.1 of Mengersen and Tweedie (1996)

Because $\beta>0$ implies that the associated Markov chain is uniformly ergodic, and from Liu (1996), we know that the spectral gap of $P$ is exactly $\beta$. To ensure that $\|P^{B}(i,\cdot)-\pi(\cdot)\|_{\rm{TV}}\leq\epsilon$ for a given $\epsilon>0$ it then suffices to choose

This requires a value or lower bound on $\beta$, whose computation we would like to avoid. The bound $1/N$ in (3) is too weak, as it leads to setting $B$ roughly as a multiple of $N$ for large $N$. It is sensible instead to choose $\beta$ as some estimate of

where $\bar{w}$ and $w_{\text{max}}$ are respectively the mean of the weights and an upper bound on the weights. The serial complexity of the Metropolis resampler is $\mathcal{O}(NB)$, but $B$ may itself be a function of $N$ and the distribution of weights, as in the analysis above.

2 Rejection resampling

When an upper bound on the weights is known a priori, rejection sampling is possible. Like the Metropolis resampler, the rejection resampler avoids collective operations and associated numerical instability, but offers a couple of additional advantages:

it permits a first deterministic proposal that $a^{i}=i$, increasing the probability of this outcome, and reducing the variance in the ancestry vectors produced.

Pseudocode is given in Code 3. If line 2 is replaced with $j\sim\mathcal{U}\{1,\ldots,N\}$ (forgoing the second advantage above), rejection resampling is an alternative implementation of multinomial resampling. Its serial complexity is then $\mathcal{O}(Nw_{\text{max}}/\bar{w})$.

An issue unique to the rejection resampler is that the computational effort required to draw each ancestor varies, depending on the number of rejected proposals before acceptance. This is an example of a variable task-length problem (Murray, 2012), particularly acute in the GPU context. On GPUs, threads are grouped into warps and execute the same instructions in parallel. The threads in the same warp may trip the loop on line 3 of Code 3 a different number of times. All threads in the warp must complete before any thread can proceed beyond the loop. This is a particular case of warp divergence, which harms performance. A persistent threads strategy (Aila and Laine, 2009; Murray, 2012) might be used to mitigate the effects of this, although we have not been successful in finding such an implementation that does not lose more than it gains through additional overhead in register use and branching.

If line 2 of Code 3 is modified so that $j$ is sampled uniformly on $\{1,\ldots,N\}$ then the number of iterations in the while loop is a geometric random variable with success probability given by $p:=\beta\times\max_{i}w^{i}/w_{\text{max}}$, where $\beta$ is the same as that of (3), perhaps also chosen as an estimate of (5). The rejection resampler will perform poorly if this probability is small, which can occur, e.g., when $\max_{i}w^{i}\ll w_{\text{max}}$. One could use the empirical maximum, $w_{\text{max}}=\max\{w^{1},\ldots,w^{N}\}$, but this would require a collective operation over weights that would defeat the purpose of the approach and, moreover, $\beta$ could still be small.

Since $H_{N}:=\sum_{k=1}^{N}\frac{1}{k}$ is the $N$th harmonic number, we additionally have the bounds

Performance can be tuned if one is willing to concede a weighted outcome from the resampling step, rather than the usual unweighted outcome. This is the approach taken with the partial rejection control heuristic (Liu et al., 1998). To do this, choose some $v_{\text{max}}<w_{\text{max}}$, then form a categorical distribution using the weights $v^{1},\ldots,v^{N}$, where $v^{i}=\min(w^{i},v_{\text{max}})$. Clearly $v_{\text{max}}$ forms an upper bound on these new weights. One could sample from this using Code 3, with $v$ in place of $w$, and then importance weight each particle $i$ with $w^{i}\leftarrow w^{a^{i}}/v^{a^{i}}$. Note that each weight is 1 except where $w^{a^{i}}>w_{\text{max}}$. The procedure may also be used when no hard upper bound on weights exists ($w_{\text{max}}$), but where some reasonable substitute can be made ($v_{\text{max}}$).

3 Ancestor permutation for in-place propagation

A desirable feature of a resampling scheme is for it to allow in-place propagation of the particle system. This is useful for memory-intensive applications of SMC in general (Bouchard-Côté et al., 2012), but particularly for GPU implementations, since the memory available to GPU devices is typically much smaller than main memory. Instead of having an input buffer holding particles at time $t-1$, and a separate output buffer into which to write the propagated particles at time $t$, a single buffer is used with the time $t$ particles replacing the time $t-1$ particles. This in-place operation is more memory efficient by a factor of two (for a fixed number of particles), but requires guarantees that the reads and writes on the single buffer can be executed concurrently without conflicts. That is, each particle is either read from or written to, but not both.

To work in-place and prevent read and write conflicts, it is sufficient that the ancestry vector, $\mathbf{a}_{t}$, satisfies $\forall i\in\{1,\ldots,N\}$:

With this, it is possible to insert a copy step immediately before each propagation step, setting $\mathbf{x}_{t-1}^{i}\leftarrow\mathbf{x}_{t-1}^{a_{t}^{i}}$ for all $i\in\{1,\ldots,N\}$ where $a_{t}^{i}\neq i$. Each particle can then be propagated in-place by reading from and writing to the same buffer. The ancestry vector produced by resampling schemes will not typically satisfy (9), but a permutation of it will. In Appendix C we describe an algorithm to re-order the ancestry vector to achieve this, in both serial and parallel settings.

Results and discussion

The resampling algorithms are assessed empirically for bias, mean squared error and execution time. Single precision is used for all experiments, in order to highlight some of the numerical issues arising when standard resampling schemes are used with large numbers of particles. While the use of double precision eliminates the numerical artifacts in the results, the ranking of algorithms by execution time is unaffected.

Experiments are conducted on two devices. The first is an eight-core Intel Xeon E5-2650 CPU, compiling with the Intel C++ Compiler version 12.1.3, using OpenMP to parallelise over eight threads. The second device is an NVIDIA K20 GPU hosted by the same CPU, compiling with CUDA 5.0 and the same version of the Intel compiler. All compiler optimisations are applied. In particular, we use the -arch sm_35 option to the CUDA compiler to target the specific architecture of the NVIDIA K20.

Resampling algorithms are often assessed using the mean squared error (MSE, see e.g. Kitagawa (1996)), computed from the offspring vector $\mathbf{o}$ and weight vector $\mathbf{w}$. The squared error (SE) of a particular offspring vector $\mathbf{o}_{k}$ is:

Weight sets are simulated to assess the speed and accuracy of each resampling algorithm. For a number of particles $N$ and observation $y$, a weight set is generated by sampling $x^{i}\sim\mathcal{N}(0,1)$, for $i=1,\ldots,N$, and setting

The construction is analogous to having a prior distribution of $x\sim\mathcal{N}(0,1)$ and likelihood function of $y\sim\mathcal{N}(x,1)$. As $y$ increases, the relative variance in weights does too. For this set up, the maximum weight is $w_{\text{max}}=1/\sqrt{2\pi}$, and the expected weight

These are used to set the number of steps for the Metropolis resampler according to the analysis in Section 2.1. Using $\epsilon=1/100$ and $\beta=\bar{w}/w_{\text{max}}$, we set $B=B^{*}$ as defined in (4). The maximum weight $w_{\text{max}}$ is also used for the rejection resampler. This procedure is used to generate 16 different weight vectors for each combination of $N=2^{4},2^{5},\ldots,2^{22}$ and $y=0,\frac{1}{2},1,1\frac{1}{2},\ldots,4$. For each of these 16 weight vectors, each resampling algorithm is used to draw 256 offspring vectors. Results reported below are averages over the 16 weight vectors.

2 Bias results

Figure 2 plots the contribution of the empirical bias to the MSE for all algorithms. For the multinomial, stratified and systematic resamplers, this appears satisfactory until $N\geq 2^{19}$, after which the bias contribution increases rapidly. This is due to the numerical instability of the cumulative sum required by these algorithms. The instability is noticeably worse for the different multinomial algorithm used on the CPU (Code 6 in Appendix B.1) than that on the GPU (Code 5 in Appendix B.1); this is explained by the cumulative sum in the former being linear over a vector rather than recursive over a binary tree. The Metropolis and rejection algorithms do not share this instability, and otherwise empirically match the bias contribution of the other methods. This is to be expected for the rejection algorithm. For the Metropolis algorithm it suggests that the procedure for setting $B$ in Section 2.1 is appropriate. It also suggests that while the Metropolis algorithm is theoretically biased, this bias is negligible in practice compared to numerical errors.

We observe, but do not show, that the pre-sorting of weights does not fix the numerical instability of the multinomial, stratified and systematic resamplers for large numbers of particles, but that the use of double precision does. The number of particles that would be required to reproduce the same instability in double precision far surpasses, by orders of magnitude, that which would be realistic to use in SMC at present.

The selection of $B$ for the Metropolis resampler appears sufficient, but we may question whether it is too conservative. To test this empirically we compare runs of the Metropolis resampler with reduced number of steps, setting $B=B^{*}/C$ for each $C\in\{1,2,4,8\}$. The contribution of the empirical bias to the MSE is given in the leftmost plot of Figure 3. As it matches that of the multinomial resampler for $C=1$—which cannot be improved upon—and noticeably increases for $C\geq 2$, this suggests that the setting $B=B^{*}$ is indeed appropriate.

3 MSE results

Figure 4 indicates that MSE differs between methods. Little in this figure is surprising, however: it is well known that the stratified resampler reduces variance over the multinomial resampler, and that the systematic resampler can, but does not necessarily, reduce it again (Douc and Cappé, 2005). Numerical instabilities in the CPU implementation of the multinomial resampler (Code 6 in Appendix B.1) appear to increase the MSE in its outcomes for large $N$. Also of interest is that as $y$ increases, the probability of accepting the initial proposal of the rejection resampler declines, so that its MSE degrades away from that of the systematic and stratified resamplers, towards that of the multinomial and Metropolis resamplers.

4 Execution time results

Figure 5 shows the execution times for all algorithms, as well as, for context, the execution times of procedures for sorting a weight vector and computing its effective sample size (ESS) (Liu and Chen, 1995). The ESS is given by $\textnormal{Sum}(\mathbf{w})^{2}/\mathbf{w}^{T}\mathbf{w}$.

Execution times are taken until the delivery of an ancestry vector satisfying (9), and so include any of the auxiliary functions in Appendices D and C necessary to achieve this. Note that—as we would expect—the multinomial, stratified and systematic resamplers are not sensitive to $y$ (or equivalently to the variance in weights) with respect to execution time, while the Metropolis and rejection resamplers are.

In a Bayesian decision theoretic setting, we can adopt execution time as a loss function, and choose, for any combination of $N$ and $y$, the algorithm that minimises the expectation of this loss function. A more sophisticated loss function might include the bias and variance as well, but the relative weighting of the individual components is a subjective decision for the problem at hand, so we do not attempt to do this. Using execution time alone as a loss function, Figure 6 plots the resulting decision matrices across all combinations of $N$ and $y$ for which empirical results were recorded. From these matrices and Figure 5, we can conclude:

that the GPU should generally not be considered for resampling with fewer than $2^{10}$ particles,

that the systematic resampler is a good candidate overall, but

that there is a significant region of the space, especially at lower weight variances, for which the rejection or Metropolis resamplers are faster.

We observe, but do not show, that the decision boundaries in Figure 6 are not significantly affected by including the time taken to copy between GPU device and main memory. This means that the choice between the CPU or GPU device for resampling is largely independent of the choice of device for the propagation and weighting of particles. For example, use of the GPU for propagating and weighting particles does not then greatly favour the GPU for resampling: the penalty to copy the weight vector to main memory, resample using the CPU, and then copy the resulting ancestry vector back to device memory, is not significant.

The choice of $B$ for the Metropolis algorithm permits a trade off between bias and execution time. This may be particularly useful in applications with hard execution time constraints, such as real-time object tracking. Recall that execution time is linear in $B$. Execution time results for various choices of $B$ are given in the rightmost plot of Figure 3, with associated biases in the leftmost plot. Recall that the rejection resampler is also somewhat configurable by using an approximate maximum weight. To do this, one must be willing to accept a still-weighted output from the resampling step, and the cumulative implications of this within SMC are problem-specific and not overly clear. We leave this for future work.

A further consideration is that the execution time of both the Metropolis and rejection resamplers depends on the PRNG used. This dependence is by a constant factor, but can be substantial. Here, robust PRNGs for Monte Carlo work have been used (see Appendix E), but conceivably cheaper, if less robust, PRNGs might be considered. This represents another trade-off between execution time and bias.

Conclusion

This work has presented two alternative resampling schemes for SMC that eliminate collective operations over weights. Consequently, they are more readily parallelised and are more numerically stable than standard resampling schemes.

The appropriate choice of resampling algorithm depends on a number of problem-specific factors, including:

the number of particles required and the typical variability in their associated weights,

whether a maximum weight exists, or can be approximated sufficiently accurately, to configure and use the Metropolis and rejection algorithms,

the tolerable level of bias in resampling outcomes, and

are demonstrably, in certain circumstances, superior in terms of execution time and numerical stability (see e.g. Figures 6 and 4, respectively), and

allow a practitioner to implement SMC in single precision for large numbers of particles without the numerical instabilities of existing resampling schemes.

With respect to numerical stability, especially in single precision, great care should be taken when using the standard multinomial, stratified or systematic resamplers with upwards of hundreds of thousands of particles, even though these schemes are unbiased in theory. This is due to numerical instability in the cumulative sum operation that these algorithms require. The alternative Metropolis and rejection resamplers have better numerical properties, as they compute only ratios of weights. This is important in light of the temptation to use single-precision, or even custom-precision floating point to improve execution times with modern computer architectures.

Supplementary materials

LibBi package Resampling providing scripts to reproduce the numerical results of this article using all algorithms as implemented in the LibBi software, available at www.libbi.org.

Acknowledgements

The third author acknowledges EPSRC for funding this research through grant EP/K009362/1.

References

Appendix A Pseudocode conventions

The algorithms presented in this work are described using pseudocode with a number of conventions. We distinguish between the for each and for constructs. The former is used where the body of the loop is to be executed for each element of a set, with the order unimportant. The latter is used where the body of the loop is to be executed for each element of a sequence, where the order must be preserved. The intended implication is that for each loops may be parallelised, while for loops cannot be. The atomic keyword is used to indicate that a line must be executed as if it constitutes one instruction (i.e. an atomic operation) in order to avoid read and write conflicts between concurrently running threads.

A number of primitive operations such as searches, transformations, reductions, sorts and prefix sums are used throughout pseudocode. These are specified in Code 4. Such operations will be familiar to users of, for example, the C++ standard template library (STL) or Thrust library (Hoberock and Bell, 2010), and their implementation on GPUs has been well-studied (see e.g. Harris et al., 2007; Satish et al., 2009). The advantage of describing algorithms in this way is that we can specify intent without prescribing implementation; the efficient implementation of these primitives in both serial and parallel contexts is well understood, and a single pseudocode description that uses primitives will often suffice for both serial and parallel contexts.

Appendix B Standard resampling schemes

Multinomial resampling proceeds by drawing each $a^{i}$ independently from the categorical distribution over $\mathcal{C}=\{1,\ldots,N\}$, where $P(a^{i}=j)=w^{j}/\textnormal{Sum}(\mathbf{w})$. Pseudocode is given in Code 5. The algorithm is dominated by the $N$ calls of Lower-Bound, which if implemented with a binary search, will give a serial complexity of $\mathcal{O}(N\log_{2}N)$ overall.

The Inclusive-Prefix-Sum operation on line 5 of Code 5 is not numerically stable, as large values may be added to relatively insignificant ones during the procedure (an issue intrinsic to any large summation). With large $N$, assigning the weights to the leaves of a binary tree and summing with a depth-first recursion over this will help. With large variance in weights, pre-sorting may also help. While log-weights are often used in the implementation of SMC, these need to be exponentiated (perhaps after rescaling) for the Inclusive-Prefix-Sum operation, so this does not alleviate the issue.

Serially, the same approach may be used, although a single-pass approach of complexity $\mathcal{O}(N)$ is enabled by generating sorted uniform random variates (Bentley and Saxe, 1979). Code 6 details this approach. A drawback is the use of relatively expensive logarithm functions. There is scope for a small degree of parallelism in this new algorithm by dividing $N$ among a handful of threads. Each thread must still step through all $N$ weights, however, so that the complexity is not improved with parallelism. We find it faster than Code 5 when on CPU, but slower when on GPU.

B.2 Stratified resampling

The variance in outcomes produced by the multinomial resampler may be reduced (Douc and Cappé, 2005) by stratifying the cumulative probability function of the same categorical distribution, and randomly drawing one particle from each stratum. This stratified resampler (Kitagawa, 1996) most naturally delivers not the ancestry vector $\mathbf{a}$ or offspring vector $\mathbf{o}$, but the cumulative offspring vector, which we denote $\mathbf{O}$, and define as $\mathbf{O}=\textnormal{Inclusive-Prefix-Sum}(\mathbf{o})$. Pseudocode is given in Code 7. The algorithm is of serial complexity $\mathcal{O}(N)$.

As for multinomial resampling, the Inclusive-Prefix-Sum operation on line 7 of Code 7 is not numerically stable. The same strategies to ameliorate the problem apply. Line 7 of Code 7 is more problematic. Consider that there may be a $j$ such that, for $i\geq j$, $u^{k^{i}}$ is not significant against $r^{i}$ under the floating-point model, so that the result of $r^{i}+u^{k^{i}}$ is just $r^{i}$. For such $i$, no random sample is being made within the strata. Furthermore, rounding up on the same line might easily deliver $O^{N}=N+1$, not $O^{N}=N$ as required, if not for the quick-fix use of $\min$. Given that single precision has about seven significant figures in decimal, consider that, with $N$ around one million, almost certainly no $u^{k^{i}}\in[0,1)$ is significant against $r^{i}$ at high $i$. Note that while pre-sorting weights and summing over a binary tree can help with the numerical stability of the Inclusive-Prefix-Sum operation, it does not help with this latter issue.

B.3 Systematic resampling

The variance in outcomes of the stratified resampler may often, but not always (Douc and Cappé, 2005), be further reduced by using the same random offset within each stratum. This is the systematic resampler (equivalent to the deterministic method described in the appendix of Kitagawa (1996)). Pseudocode is given in Code 8, which is a simple modification to Code 7. The same complexity and numerical caveats apply to the systematic resampler as for the stratified resampler.

The resampling algorithms presented here do not constitute an exhaustive list of those in use, for instance residual resampling has been omitted (Liu and Chen, 1998). However they are reasonably representative, and can form the building blocks of more elaborate schemes.

Appendix C Ancestor permutation for in-place propagation

An ancestry vector may be permuted to satisfy (9) in the main article. A serial algorithm to achieve this is straightforward and given in Code 9. This $\mathcal{O}(N)$ algorithm makes a single pass through the ancestry vector with pair-wise swaps to satisfy the condition.

The simple algorithm is complicated in a parallel context as the pair-wise swaps are not readily serialised without heavy-weight mutual exclusion. In parallel we propose Code 10. This algorithm does not perform the permutation in-place, but instead produces a new vector $\mathbf{c}\in\{1,\ldots,N\}^{N}$ that is the permutation of the input vector $\mathbf{a}$. It introduces a new vector $\mathbf{d}\in\{1,\ldots,N+1\}^{N}$, through which, ultimately, $c^{i}=a^{d^{i}}$. In the first stage of the algorithm, Prepermute, the thread for element $i$ attempts to claim position $a^{i}$ in the output vector by setting $d^{a^{i}}\leftarrow i$. By virtue of the $\min$ function on line 10, the element of lowest index always succeeds in this claim while all others contesting the same place fail, and the outcome of the whole permutation procedure is deterministic. This is desirable so that the results of a particle filter are reproducible for the same pseudorandom number seed. For each element $i$ that is not successful in its claim, the thread for $i$ instead attempts to claim $d^{i}$, if unsuccessful again then $d^{d^{i}}$, then recursively $d^{d^{d^{i}}},\ldots$ etc, until an unclaimed place is found.

We offer a proof of the termination of Code 10. First note that Prepermute leaves $\mathbf{d}$ in a state where, excluding all values of $N+1$, the remaining values are unique. Furthermore, in Permute the conditional on line 10 means that the loop on line 10 is only entered for values of $i$ that are not represented in $\mathbf{d}$.

For each such $i$, the while loop traverses the sequence $x_{0}=i$, $x_{n}=d^{x_{n-1}}$, until $d^{x_{n}}=N+1$. For the procedure to terminate this sequence must be finite. Because each $x_{n}$ is an element of the finite set $\{1,\ldots,N\}$, to show that the sequence is finite it is sufficient to show that it never revisits the same value twice. The proof is by induction.

As no value of $\mathbf{d}$ is $i$, the sequence cannot revisit its initial value $\mathbf{x}_{0}=i$. The element $x_{0}$ is therefore unique.

For $k\geq 1$, assume that the elements of $x_{0:k-1}$ are unique.

Now, the elements of $x_{0:k}$ are not unique if there exists some $j\in\{1,\ldots,k-1\}$ such that $x_{k}=d^{x_{k-1}}=x_{j}=d^{x_{j-1}}$, with $x_{j-1}\neq x_{k-1}$ by the uniqueness of $x_{0:k-1}$. But this contradicts the uniqueness of the (non $N+1$) values of $\mathbf{d}$. Thus the elements of $x_{0:k}$ are unique, the sequence is finite, and the program must terminate.

Appendix D Auxiliary functions

The multinomial, Metropolis and rejection resamplers most naturally return the ancestry vector $\mathbf{a}$, while the stratified and systematic resamplers return the cumulative offspring vector $\mathbf{O}$. Conversion between these is reasonably straightforward. An offspring vector $\mathbf{o}$ may be converted to a cumulative offspring vector $\mathbf{O}$ via the Inclusive-Prefix-Sum primitive, and back again via Adjacent-Difference. A cumulative offspring vector may be converted to an ancestry vector via Code 11, and an ancestry vector to an offspring vector via Code 12. These functions perform well on both CPU and GPU. An alternative approach to Cumulative-Offspring-To-Ancestors, using a binary search for each ancestor, was found to be slower.

Appendix E Implementation

All the algorithms described in the article and the appendices have been implemented as part of the LibBi software (www.libbi.org, Murray (2013)) for performing methods such as the particle filter on high-performance computing devices. We enumerate the most important considerations of the implementation here, and avoid painstaking detail of the remainder so as not to oversell their importance relative to these. It is worth emphasising that some important decisions, such as the choice of pseudorandom number generator (PRNG), depend on the particular problem at hand.

The weight vector may contain many very small values. Because of this, a typical implementation will store log-weights rather than weights for numerical accuracy. The log-weights may be large and negative, and one should avoid taking a floating point exponential of these large negative numbers, which is often zero. All of the algorithms presented are robust to the scaling of weights by a constant factor, however. When computing sums or prefix sums, a vector of log-weights can therefore be renormalised using, say, the maximum value, denoted $\log w_{\text{max}}$. For example, the logarithm of the sum of weights, stored as log-weights, is accurately computed using the identity:

Renormalisation is not required for the Metropolis and rejection algorithms, as they feature only pairwise ratios between weights, or pairwise differences between log-weights.

The performance of the multinomial, stratified and systematic resamplers depends largely on the implementation of the prefix sum operation. We defer to existing work for these operations, in particular to that invested in the Thrust library (Bell and Hoberock, 2012), which the implementation uses. Conceptually, the implementation in the Thrust library is based on up- and down-sweeps of a balanced binary tree (Harris et al., 2007).

The performance of the Metropolis and rejection resamplers is dependent mostly on the selection of PRNG. Performance is not the only consideration in this selection, however. PRNGs are assessed both on execution speed and the statistical quality of the pseudorandom number sequence that they produce, typically using test suites such as DIEHARD (Marsaglia, 1996) or TestU01 (L’Ecuyer and Simard, 2007). Among clients of PRNGs, Monte Carlo algorithms, such as the particle filter, have high demands for statistical quality. To this end, our CPU code uses the Mersenne Twister PRNG (Matsumoto and Nishimura, 1998) as implemented in the Boost.Random library (www.boost.org). This is standard for Monte Carlo applications. Our GPU code uses the XORWOW PRNG (Marsaglia, 2003) from the CURAND library (NVIDIA Corporation, 2012). This particular PRNG belongs to a family that is readily shaped to the GPU architecture (Nandapalan et al., 2012). Faster but lower quality PRNG may be used. This would constitute a relaxing of the unbiasedness condition (2). As any such decision is problem-specific, it is not investigated in this work.

The Metropolis and rejection algorithms use random access patterns to memory. Spatiotemporally local access patterns are preferred for good cache performance on CPU, and streaming, or at least coalesced access, is preferred on GPU. The random access pattern is, unfortunately, inherent to the algorithms, and we can only rely on the presence of a large cache to mitigate associated latencies. On GPU, juditious use of shared memory may help, but there is no reason to believe that this can achieve better results than the hardware-controlled cache found on more recent architectures; we rely on the latter. As such, the GPU is configured to use 48 KB of L1 cache and 16 KB of shared memory. This maximises the size of the cache for random access patterns, but still provides sufficient shared memory for all kernels.

Our Metropolis and rejection resampler kernels compile to 32 registers per thread, as reported by the CUDA compiler. This is satisfactory with respect to occupancy of the device, and we do not seek further reductions.

Finally, the auxiliary algorithms presented in §D pose little challenge. Implemented using CUDA, they compile to kernels using no shared memory and fewer than 16 registers per thread, which is of no hindrance to occupancy of the device. On GPU, we append the Prepermute procedure of Code 10 to the end of any procedure that produces an ancestry vector. This saves the launch of a separate kernel and the associated overhead of doing so.