PaperlayerView exact revision on arXiv

Source-linked AI summary

Parallel resampling in the particle filter

Lawrence M. Murray, Anthony Lee, Pierre E. Jacob

arXiv:1301.4019v3stat.COcs.DC

TL;DR

SMC resampling is difficult to parallelise because standard schemes require collective operations across particle weights. The paper evaluates Metropolis and rejection resamplers without such operations, finding circumstances where they improve execution time and numerical stability while introducing method-specific trade-offs.

  • Problem

    Standard resampling requires cumulative operations over particle weights, creating synchronization and numerical-stability problems that hinder parallel SMC.

  • Method

    The paper introduces and evaluates Metropolis and rejection resamplers that eliminate collective operations over weights, comparing them with standard schemes.

  • Results

    The alternative resamplers are superior in certain circumstances for execution time and numerical stability, and avoid the standard schemes' single-precision instability with large particle counts.

  • Takeaways & Limitations

    The alternatives enable SMC implementation in single precision for large particle counts without the numerical instabilities of existing resampling schemes.

  • Takeaways & Limitations

    Finite-step Metropolis resampling is biased, creating a speed–bias trade-off and potentially violating assumptions needed for unbiased marginal likelihood estimates in particle MCMC.

Abstract

from arXiv · show

Modern parallel computing devices, such as the graphics processing unit (GPU), have gained significant traction in scientific and statistical computing. They are particularly well-suited to data-parallel algorithms such as the particle filter, or more generally Sequential Monte Carlo (SMC), which are increasingly used in statistical inference. SMC methods carry a set of weighted particles through repeated propagation, weighting and resampling steps. The propagation and weighting steps are straightforward to parallelise, as they require only independent operations on each particle. The resampling step is more difficult, as standard schemes require a collective operation, such as a sum, across particle weights. Focusing on this resampling step, we analyse two alternative schemes that do not involve a collective operation (Metropolis and rejection resamplers), and compare them to standard schemes (multinomial, stratified and systematic resamplers). We find that, in certain circumstances, the alternative resamplers can perform significantly faster on a GPU, and to a lesser extent on a CPU, than the standard approaches. Moreover, in single precision, the standard approaches are numerically biased for upwards of hundreds of thousands of particles, while the alternatives are not. This is particularly important given greater single- than double-precision throughput on modern devices, and the consequent temptation to use single precision with a greater number of particles. Finally, we provide auxiliary functions useful for implementation, such as for the permutation of ancestry vectors to enable in-place propagation.

1 Introduction

SMC uses weighted particles for statistical inference, but its resampling step is harder to parallelise than propagation and weighting because standard schemes require collective operations over weights. This work studies alternative resamplers that remove cumulative sums to improve parallel suitability and numerical stability.

  • SMC framework: SMC represents distributions with weighted particles, repeatedly resampling, propagating and weighting them through Markov kernels and potential functions.The framework can target a sequence of distributions, including filtering distributions and parameter posteriors.
  • 1.2 Parallelisation: Resampling is a parallelisation bottleneck because standard schemes require synchronization and collective operations across particle weights, especially on GPUs.Propagation and weighting operate independently on particles, whereas summing weights can leave threads idle and impede scaling.
  • 1.2 Parallelisation: Standard multinomial, stratified and systematic resamplers rely on cumulative sums of weights, motivating prior work on parallel prefix-sum implementations.The cumulative sum is central to these schemes and is also a generic parallel-computing problem.
  • 1.4 Contributions: Standard resampling can become numerically unstable with large particle counts or highly variable weights, while single precision is attractive because it can be substantially faster on CPUs and GPUs.The paper notes that one million particles can be realistic and that single precision offers at least a two-fold GPU speedup in the cited context.
  • 1.4 Contributions: The paper introduces Metropolis- and rejection-based resamplers that eliminate cumulative sums and compares them with multinomial, stratified and systematic schemes across CPU and GPU settings.The comparison considers numerical stability, bias, mean squared error and execution time across particle counts and weight variability.

2 Alternative resampling schemes

The paper presents Metropolis and rejection resamplers that bypass cumulative weight sums and collective operations, targeting parallel hardware and numerical stability. It also develops guidance on their trade-offs and an ancestry permutation enabling in-place particle propagation.

  • 2.1 Metropolis resampling: Metropolis resampling replaces direct sampling with B pairwise weight-ratio iterations, producing results close to multinomial resampling without collective operations.The method applies a Markov kernel with the categorical weight distribution as its stationary target.
  • 2.1 Metropolis resampling: Finite B creates bias, trading execution speed against accuracy; this can invalidate unbiased marginal-likelihood assumptions in particle MCMC.Smaller B runs faster but produces larger bias.
  • 2.1 Metropolis resampling: Metropolis resampling has serial complexity O(NB), with B potentially depending on N and the weight distribution.Selecting B requires controlling the distance between the iterated kernel and its target, but computing the relevant convergence parameter is undesirable.
  • 2.2 Rejection resampling: Rejection resampling requires a known upper bound on weights and avoids collective operations, while a deterministic initial proposal ai = i can reduce ancestry variance.Replacing the initial proposal with uniform sampling yields an alternative implementation of multinomial resampling.
  • 2.2 Rejection resampling: Rejection resampling suffers GPU warp divergence because different ancestors require different numbers of rejected proposals before acceptance.The resulting variable task lengths can make serial complexity less attractive than parallel complexity.
  • 2.2 Rejection resampling: Rejection performance deteriorates when its acceptance probability is small, while its expected parallel complexity with N processors scales as O(log N/p).Using the empirical maximum weight would require a collective operation and defeat the approach’s purpose.
  • 2.3 Ancestor permutation for in-place propagation: A permutation of the ancestry vector permits in-place propagation by ensuring each particle is read from or written to, but not both.A preceding copy step handles non-self ancestries before propagation in the shared buffer.

3 Results and discussion

Experiments compare resampling algorithms on CPU and GPU using bias, MSE, and execution time. Standard schemes become numerically unstable at large particle counts in single precision, while alternative schemes can be faster under some conditions.

  • 3.2 Bias results: In single precision, multinomial, stratified, and systematic resamplers develop rapidly increasing bias from cumulative-sum instability at large N, whereas Metropolis and rejection resamplers do not.The instability appears after about N ≥ 2^19 and is worse for the CPU multinomial implementation; double precision eliminates the numerical artifacts.
  • 3.2 Bias results: The Metropolis resampler's empirically negligible bias supports the chosen step count B = B*, despite finite-step theoretical bias.Reducing B to B*/C increases bias for C ≥ 2, suggesting the default setting is appropriate.
  • 3.3 MSE results: Stratified resampling generally lowers MSE relative to multinomial resampling, while systematic resampling can reduce it further; rejection MSE worsens as y increases.The rejection resampler moves toward multinomial and Metropolis MSE as the probability of accepting the initial proposal declines.
  • 3.4 Execution time results: The GPU should generally not be used for resampling below 2^10 particles, while systematic is a strong overall candidate and rejection or Metropolis can be faster in substantial regions.The alternative resamplers are especially competitive at lower weight variances.
  • 3.4 Execution time results: Including CPU–GPU memory-copy time does not significantly change device-selection boundaries, so propagation on the GPU does not strongly favor GPU resampling.Copying weights to main memory, resampling on the CPU, and copying ancestry back is reported as not significant.
  • 3.3 MSE results: Metropolis execution time scales linearly with B, creating a configurable trade-off between bias and runtime for time-constrained applications.Reducing B lowers execution time but increases bias and slightly increases MSE.

4 Conclusion

The paper introduces Metropolis and rejection resamplers that avoid collective weight operations, improving parallelisation and numerical stability. Their suitability depends on particle count, weight variability, bias tolerance, available precision, and whether a usable maximum weight is available.

  • The appropriate resampling algorithm depends on particle count, weight variability, maximum-weight availability, tolerable bias, and numerical precision.
  • Metropolis and rejection resamplers can be superior in execution time and numerical stability, and support single-precision SMC with large particle counts.These advantages are demonstrated under particular conditions rather than universally.
  • Standard multinomial, stratified, and systematic resamplers require care in single precision with upwards of hundreds of thousands of particles because cumulative sums can become numerically unstable.Metropolis and rejection resamplers instead compute weight ratios, giving them better numerical properties in this setting.

Supplementary materials

The supplementary materials provide reproducibility scripts and decision matrices for choosing resampling algorithms by particle count and weight variability.

  • The LibBi package provides scripts to reproduce the article’s numerical results using all implemented algorithms.The package is available through libbi.org.
  • Figure 6 maps the fastest algorithm across particle counts and weight variability for CPU-only, GPU-only, and combined CPU/GPU choices.The figure also considers whether copying weights and ancestry vectors changes the combined-device decision boundaries.

A Pseudocode conventions

The pseudocode conventions distinguish loops that may be parallelised from ordered loops and define reusable primitive operations without prescribing a particular implementation.

  • “For each” loops indicate bodies whose execution order is unimportant and may therefore be parallelised, whereas “for” loops preserve sequence order.
  • The pseudocode uses familiar primitives such as searches, transformations, reductions, sorts, and prefix sums for both serial and parallel contexts.
  • Lower-Bound returns the lowest insertion position for a value in an ascending-sorted weight vector.

B.1 Multinomial resampling

Multinomial resampling uses cumulative weight sums and lower-bound searches, while a single-pass variant offers linear serial complexity but has different CPU and GPU performance.

  • Multinomial resampling draws each ancestor independently from a categorical distribution defined by normalised weights, using a prefix sum and lower-bound search.With binary search, the serial complexity is O(N log2 N).
  • The inclusive prefix sum is numerically unstable because large values can absorb relatively insignificant additions, especially as particle counts increase.Tree-based summation and pre-sorting may help in some settings, but exponentiating log-weights does not remove the issue.
  • A sorted-uniform, single-pass multinomial variant has O(N) serial complexity but is faster than the parallel approach on CPU and slower on GPU.Its implementation uses an exclusive prefix sum and a sequential search procedure.

B.2 Stratified resampling

Stratified resampling reduces multinomial variance by sampling one particle from each stratum, but relies on a serial O(N) prefix sum and has numerical-instability risks.

  • Stratified resampling reduces multinomial variance by randomly drawing one particle from each stratum of the cumulative probability function.
  • The algorithm naturally returns a cumulative offspring vector O, defined as the inclusive prefix sum of the offspring vector o.
  • The serial complexity is O(N), and the procedure is specified in Code 7.
  • Floating-point rounding can prevent sampling within later strata and may produce O_N = N + 1.

B.3 Systematic resampling

Systematic resampling uses one shared random offset across strata to often reduce variance further, while retaining stratified resampling’s complexity and numerical caveats.

  • Systematic resampling often further reduces stratified resampling’s variance by using the same random offset within every stratum.
  • The method is a simple modification of Code 7 and is specified in Code 8.
  • Systematic resampling retains the stratified method’s complexity and numerical caveats.
  • The resampling schemes discussed here are representative rather than exhaustive because residual resampling is omitted.

C Ancestor permutation for in-place propagation

The paper replaces difficult in-place parallel swaps with a deterministic parallel ancestry permutation that resolves competing claims and is proved to terminate.

  • The serial O(N) swap algorithm is difficult to parallelize because pair-wise swaps require heavy-weight mutual exclusion.
  • The proposed parallel permutation creates a new vector c rather than performing ancestry swaps in place.
  • Prepermute assigns competing output positions to the lowest-index claimant, making the resulting permutation deterministic and reproducible for a fixed pseudorandom seed.
  • The termination proof shows that the traversal sequence never revisits a value and therefore remains finite.
  • Ancestry, offspring, and cumulative offspring vectors can be converted using prefix-sum, adjacent-difference, and the stated conversion procedures.

E Implementation

The implementation combines CUDA-oriented ancestry permutation and vector-conversion routines with numerically careful weight handling, prefix-sum standard resamplers, and PRNG-dependent alternative resamplers.

  • Code 10’s parallel permutation uses Prepermute to claim slots atomically and Permute to construct a new vector c satisfying the ancestry condition.
  • Log-weights are renormalized before sums or prefix sums, whereas Metropolis and rejection resampling require only weight ratios or log-weight differences.
  • Standard resampler performance depends largely on prefix sums, while Metropolis and rejection performance depends mostly on PRNG selection.
  • The implementation supports conversions between cumulative offspring, offspring, and ancestry vectors through prefix sums, adjacent differences, and atomic counting.
  • GPU implementation uses XORWOW, accepts inherent random memory access in alternative resamplers, and appends Prepermute to avoid a separate kernel launch.
Loading 1301.4019v3…