Asynchronous Anytime Sequential Monte Carlo

Introduction

Particle filter based inference techniques require blocking barrier synchronization at resampling steps which limits parallel throughput and is costly in terms of memory. We introduce a new asynchronous particle filter algorithm that has statistical efficiency competitive with standard resampling algorithms, and has sufficiently higher particle throughput such that it is, on balance, more efficient per unit time. The approach uses locally-computed decision rules for each particle that do not require block synchronization of all particles, instead only requiring sharing summary statistics with particles that follow. In our algorithm each resampling point acts as a queue rather than a barrier: each particle chooses the number of its own offspring using by comparing its own weight to the weights of particles which previously reached the queue, updates its own weight, then proceeds without waiting.

An anytime algorithm is an algorithm which can be run continuously, generating progressively better solutions when afforded additional computation time. Traditional particle filter (PF) algorithms are not anytime in nature; all particles need to be propagated in lock-step to completion in order to compute expectations. Once a particle set runs to termination, inference cannot straightforwardly be continued by simply doing more computation. The naive strategy of running sequential Monte Carlo (SMC) again and merging the resulting sets of particles is suboptimal due to bias (see for explanation). More complex methods (i.e. particle Metropolis Hastings and iterated conditional sequential Monte Carlo (iCSMC) ) for correctly merging particle sets produced by additional SMC runs are closer to anytime in nature but suffer from burstiness as big sets of particles are computed then emitted at once and, fundamentally, the inner-SMC loop of such algorithms still suffers the kind of excessive synchronization performance penalty that the particle cascade directly avoids. Our asynchronous SMC algorithm, the particle cascade, is anytime in nature. The particle cascade can be run indefinitely, without resorting to merging of particle sets, and with only a fixed (tunable) memory requirement.

Our algorithm shares a superficial similarity to Bernoulli branching numbers and other search and exploration methods which have been used for particle filtering, where each particle samples some number of children to propagate to the next observation. Like the particle cascade (and in contrast to most traditional resampling algorithms), the total number of particles which exist at each generation is allowed to gradually increase and decrease. However, computing branching correction numbers is also generally a synchronous operation, requiring all particle weights to be known at each observation in order to choose an appropriate number of offspring; this also precludes use of an anytime algorithm.

Parallelizing the resampling step of sequential Monte Carlo methods has drawn increasing recent interest as the effort progresses to scale up algorithms to take advantage of high-performance computing systems and GPUs. A recent approach to removing the global collective resampling operation, quite different from the particle cascade method introduced here, can be found at .

Another recent method for running arbitrarily many particles within a fixed memory budget introduced in focuses on keeping track of random seeds used to generate proposals, allowing particular particles to be deterministically “replayed”; a first pass through all particles computes the normalizing constant of the particle weights, and a second pass re-executes those which are chosen to continue to the next generation. However, the algorithm as presented there still relies on a synchronous resampling step, and lacks the anytime property of our approach.

Background

We begin by briefly reviewing particle filtering as generally formulated on state-space models. Suppose we have a non-Markovian dynamical system with latent random variables $X_{0},\ldots,X_{N}$ and observed random variables $Y_{0},\ldots,Y_{N}$ described by the joint density

where $X_{0}$ is drawn from some initial distribution $\mu(\cdot)$, and $f$ and $g$ are conditional densities.

Given observed values $Y_{0:N}=y_{0:N}$, we approximate the posterior distribution $p(X_{0:n}|y_{0:n})$ with a weighted set of $K$ particles, with each particle $k$ denoted $x^{k}_{0:n}$ for $k=1,\dots,K$. Particles are propagated forward from proposal densities $q(x_{n}|x_{0:n-1})$ and re-weighted at each observation $n=1,\ldots,N$:

where $w^{k}_{n}$ is the weight associated with observation $y_{n}$ and $W^{k}_{n}$ is the weight of particle $k$ after observation $n$. We assume that exact evaluation of $p(x_{0:N}|y_{0:N})$ is intractable, requiring only that the conditional terms likelihoods $g(y_{n}|x^{k}_{0:n})$ can be evaluated pointwise. In many complex dynamical systems, or in black-box simulation models, evaluation of $f(x^{k}_{n}|x^{k}_{0:n-1})$ may be prohibitively costly or even impossible. As long as we are capable of simulating from the system we can set our proposal distribution $q(\cdot)\equiv f(\cdot)$, in which case the particle weights are simply $w^{k}_{n}=g(y_{n}|x^{k}_{0:n})$, eliminating the need to compute the conditional densities $f(\cdot)$ directly.

By normalizing the weights $W^{k}_{n}$, defining

we can approximate the posterior distribution $p(X_{0:N}|y_{0:N})$ with a weighted set of $K$ particles

In the very simple sequential importance sampling setup described here, the marginal likelihood can be estimated by

The algorithm described above suffers from a degeneracy problem wherein the normalized weights $\bar{\omega}^{1}_{n},\dots,\bar{\omega}^{K}_{n}$ become mostly very close to zero for even moderately large $n$. Traditionally this is combated by introducing a resampling step: as we progress from $n$ to $n+1$, particles with high weights are duplicated and particles with low weights are discarded. Many difference schemes for resampling particles exist; see for an overview, with discussion and theoretical results for several common approaches. One can think of a resampling scheme as a method for drawing the number of offspring particles $M^{k}_{n+1}$ that each particle $k$ will produce after stage $n$. After resampling, all outgoing particles from $n$ to $n+1$ receive a new outgoing weight $V^{k}_{n+1}$, and we have

In most traditional resampling schemes the outgoing weights of all $K$ particles are deterministically set to be equal, i.e. with $V^{k}_{n+1}=1/K$; a valid resampling scheme then must satisfy the unbiasedness condition

Introducing this resampling step prevents all the probability mass in our approximation to the posterior from accumulating on a single particle. In this version of the algorithm, where a resampling step is added at every $n$, the marginal likelihood can be estimated by

it is well-known that the estimate of the marginal likelihood is unbiased .

2 Limitations

Our goal is to scale up to very large numbers of particles, using a parallel computing architecture where each particle is simulated as a separate process or thread. In order to resample at each $n$ we must compute the normalized weights $\bar{\omega}^{k}_{n}$, requiring us to wait until all individual particles have both finished forward simulation and computed their individual weight $W^{k}_{n}$ before any can proceed. While the forward simulation itself is trivially parallelizable, the weight normalization and resampling step is a synchronous, collective operation. In practice this can lead to significant underuse of computing resources in a multiprocessor environment, hindering our ability to scale up to large numbers of particles.

Memory limitations on finite computing hardware also limit the number of simultaneous particles $K$ we are capable of running in practice. All $K$ particles must move through the system together, and all must exist simultaneously; if the total memory requirements of $K$ particles is greater than the available system RAM, then a substantial overhead will be incurred from regularly swapping memory contents to disk.

The Particle Cascade

The particle cascade algorithm we introduce addresses both these limitations: it does not require synchronization, and keeps only a bounded number of particles alive in the system at any given time. Instead of resampling, we will consider particle branching, where each particle can result in 0 or more offspring. These branching events happen asynchronously and mutually exclusively, i.e. they are processed one at a time.

At each stage $n$ of the particle filter, particles process observation $y_{n}$. Without loss of generality, we can define an ordering on the particles $1,2,\ldots$ in the order they arrive at $y_{n}$. This order need not be independent of the state of the particles $x^{k}_{0:n}$.

We keep track of the running average weight $\overline{W}^{k}_{n}$ of the first $k$ particles to arrive at observation $y_{n}$ in an online manner:

The number of children of particle $k$ should depend on the weight $W^{k}_{n}$ of particle $k$ relative to those of other particles. Particles with higher relative weight are more likely to be located in a high posterior probability part of the space, and should be allowed to spawn more child particles.

In the online asynchronous particle system as described here, we do not have access to the weights of future particles when processing $k$. Instead we will compare $W^{k}_{n}$ to the current average weight $\overline{W}^{k}_{n}$ among particles processed thus far. Specifically, the number of children, which we denote by $M^{k}_{n+1}$, will depend on the ratio

Each child of particle $k$ will be assigned a weight $V^{k}_{n+1}$ such that the total weight of all children $M^{k}_{n+1}V^{k}_{n+1}$ has expectation $W^{k}_{n}$.

A simple approach would be to sample $M^{k}_{n+1}$ independently conditioned on the weights. In such schemes we could draw each $M^{k}_{n+1}$ from some simple distribution, e.g. a Poisson distribution with mean $R^{k}_{n}$, or a discrete distribution over the integers $\{\lfloor R^{k}_{n}\rfloor,\lceil R^{k}_{n}\rceil\}$. However, one issue that arises in such approaches where the number of children for each particle is conditionally independent, is that the variance of the total number of particles at each generation can grow without bound. Suppose we start the system with $K_{0}$ particles. The number of particles at subsequent stages $n$ is given recursively as $K_{n}=\sum_{k=1}^{K_{n-1}}M^{k}_{n}$. We would like to avoid situations in which the number of particles becomes too large, or collapses to 1.

Instead, we will allow $M^{k}_{n}$ to depend on the number of children of previous particles at $n$, in such a way that we can stabilize the total number of particles in each generation. Suppose that we wish for the number of particles to be stabilized around $K_{0}$. After $k-1$ particles have been processed, we expect the total number of children produced at that point to be approximately $k-1$, so that if the number is less than $k-1$ we should allow particle $k$ to produce more children, and vice versa. Similarly, if we already currently have more than $K$ children, we should allow particle $k$ to produce less children. We use a simple scheme which satisfies these criteria, where the number of particles is chosen at random when $R^{k}_{n}<1$, and set deterministically when $R^{k}_{n}\geq 1$:

We pause here to take note of the anytime nature of this algorithm — any given particle passing through the system needs only the previous weights $\overline{W}^{k}_{n}$ in order to make its local branching decisions, not the previous particles themselves. Thus it is possible to run this algorithm for some fixed number of initial particles $K_{0}$, inspect the output of the $K_{N}$ completed particles which have left the system, and then decide whether to continue inference by initializing additional particles.

2 Computing expectations and marginal likelihoods

Samples drawn from the particle cascade can be used to compute expectations in the same manner as samples from a standard particle filter; that is, given some function $\varphi(\cdot)$, we normalize weights $\bar{\omega}^{k}_{n}=\frac{W^{k}_{n}}{\sum_{j=1}^{K_{n}}W^{j}_{n}}$ analogously to before and approximate the posterior expectation by

We can also use the particle cascade to define an estimator of the marginal likelihood $p(y_{0:n})$,

The form of this estimate is fairly distinct from the standard SMC estimators in Section 2. In terms of predictive densities, one can think of $\hat{p}\left(y_{0:n}\right)$ as

It is interesting to note that the incrementally updated $\overline{W}^{k}_{n}$ statistics in the denominator of $R^{k}_{n}$ are very directly tied to the marginal likelihood estimate; that is, $\hat{p}(y_{0:n})=\frac{K_{n}}{K_{0}}\overline{W}^{k}_{n}$.

3 Theoretical properties, unbiasedness

Initialization at $n=0$: for $k=1,...,K_{0}$ sample $X_{0}^{k,0}\sim\mu(\cdot)$ and compute $W_{0}^{k}=g(y_{0}|X_{0}^{k,0}).$

Resampling step: resample $\left\{W_{n}^{k},X_{0:n}^{k,n}\right\}_{k=1}^{K_{n}}$ to obtain $\left\{\widetilde{W}_{n}^{k},X_{0:n}^{k,n+1}\right\}_{k=1}^{K_{n+1}}$

Forward simulation step: for $k=1,...,K_{n+1}$ sample $X_{n+1}^{k,n+1}\sim f\left(\left.\cdot\right|X_{0:n}^{k,n+1}\right)$, and set $W_{n+1}^{k}=\widetilde{W}_{n}^{k}g\left(\left.y_{n+1}\right|X_{0:n+1}^{k,n+1},y_{0:n}\right)$ and $n\leftarrow n+1.$

We denote by $B(E)$ the space of bounded real-valued functions on a space $E$. We make the following assumption on the resampling step.

Assumption R. For any $n\geq 0,$ we have $p(K_{n}>0)=1$ and for any $\varphi\in B(\mathcal{X}^{n})$

where $\mathcal{F}_{n}$ denotes the natural filtration associated to all the random variables generated by the particle algorithm before resampling at time $n$. We also denote by $\widetilde{\mathcal{F}}_{n}$ denotes the natural filtration associated to all the random variables generated by the particle algorithm just after the resampling step at time $n$.

The resampling step of the particle cascade corresponds to

i.e. each particle $X_{0:n}^{k,n}$ has $M_{n+1}^{k}$ offspring of associated weight $V_{n+1}^{k}$ so that $K_{n+1}=\sum_{k=1}^{K_{n}}M_{n+1}^{k}.$

Proof of Proposition 1. The proof follows from a backward induction. We have

Active bounding of memory usage

In an idealized computational environment, with infinite available memory, our implementation of the particle cascade could begin by launching (a very large number) $K_{0}$ particles simultaneously which then gradually propagate forward through the system. In practice, only some finite number of particles, probably much smaller than $K_{0}$, can be simultaneously simulated efficiently. Furthermore, the initial particles are not truly launched all at once, but rather in a sequence, introducing a dependency in the order in which particles arrive at each observation $n$.

While the resampling scheme in Eq. 14 is designed to stabilize the number of particles over time, we can still see an explosion in the number of particles $K_{n}$. The degree to which the particle count becomes unstable depends on the extent to which the ordering of the particles is permuted as we progress to each $n$. In Fig. 1 we compare a best-case situation where the ordering of particles at $n$ is completely independent of the ordering of particles at $n+1$, to a worst-case situation where the ordering of particles is completely preserved from $n$ to $n+1$. In practice, a naïve implementation of the incremental resampling scheme will have a very strong dependence in ordering across $n$ — a particle which is one of the first to reach stage $n$ is quite likely one of the first to reach stage $n+1$ as well.

Our implementation of the particle cascade addresses these issues by explicitly injecting randomness into the execution order of particles, and by imposing a machine-dependent hard cap on the number of simultaneous extant processes. This permits us to run our particle filter system indefinitely, for arbitrarily large initial particle counts $K_{0}$, while consuming only a fixed computational budget.

Each particle in our implementation runs as an independent operating system process. In order to efficiently run a large number of particles, we impose a hard limit limit $\rho$ on the total number of particles which can simultaneously exist in the particle system; most of these will generally be sleeping processes. The ideal choice for this number will vary based on hardware capabilities, but in general should be made as large as possible.

Scheduling across particles is managed via a global first-in random-out process queue of length $\rho$; this can equivalently be conceptualized as a random-weight priority queue. Each particle corresponds to a single live process, augmented by a single additional control process which is responsible only for spawning additional initial particles (i.e. incrementing the initial particle count $K_{0}$). When any particle $k$ arrives at any likelihood evaluation $n$, it computes its target number of child particles $M_{n+1}^{k}$ and outgoing particle weight $V_{n+1}^{k}$. If $M_{n+1}^{k}=0$ it immediately terminates; otherwise it enters the queue. Once this particle either enters the queue or terminates, some other process continues execution — this process is chosen uniformly at random, and as such may be a sleeping particle at any stage $n<N$, or it may instead be the control process which then launches a brand new particle. At any given time, there are some number of particles $K_{\rho}<\rho$ currently in the queue, and so the probability of resuming any particular individual particle, or of launching a new particle, is $\frac{1}{K_{\rho}+1}$. If the particle released from the queue has exactly one child to spawn, it advances to the next observation and repeats the resampling process. If, however, a particle has more than one child particle to spawn, rather than launching all child particles at once it launches a single particle to simulate forward, decrements the total number of particles left to launch by one, and itself re-enters the queue.

In the event that the process count is fully saturated (i.e. the process queue is full), then we forcibly prevent particles from duplicating themselves and creating new children. If we release a particle from the queue which seeks to launch $m>1$ additional particles when the queue is full, we instead collapse all the remaining particles into a single particle; this single particle represents a “virtual” set of particles, but does not actually create a new process and requires no additional CPU or memory resources. We keep track of a particle count multiplier $C^{k}_{n}$ that we propagate forward along with the particle. All particles are initialized with $C^{k}_{0}=1$, and then when a particle collapse takes place, update their multiplier at $n+1$ to $mC^{k}_{n}$.

This affects the way in which running weight averages are computed; suppose a new particle $k$ arrives with multiplier $C^{k}_{n}$ and weight $W^{k}_{n}$. We incorporate all these values into the average weight immediately, and update $\overline{W}^{k}_{n}$ taking into account the multiplicity, with

This does not affect the computation of the ratio $R^{k}_{n}$. We preserve the particle multiplier, until we reach the final $n=N$; then, after all forward simulation is complete, we re-incorporate the particle multiplicity when reporting the final particle weight $W^{k}_{N}=C^{k}_{N}V^{k}_{N}w^{k}_{N}$. The system is initialized by seeding the system with a number of initial particles $\rho_{0}<\rho$ at $n=0$, creating $\rho_{0}$ active initial processes.

The ideal choice for the process count constraint $\rho$ may vary across operating systems and hardware configurations; online optimization of this parameter remains an avenue for future work.

Experiments

We run a preliminary set of experiments on two simple state space models, each with $N=50$ observations, with the goal of demonstrating the overall validity and utility of the particle cascade algorithm. Results are presented here on two simple models. The first is a hidden Markov model (HMM) with 10 latent discrete states, each with an associated Gaussian emission distribution; the second is a one-dimensional linear Gaussian model. In both models we can use an exact algorithm to compute posterior marginals at each $n$ and compute the marginal likelihood $Z=p(y_{1:N})$.

These experiments are not designed to stress-test the particle cascade; rather, they are designed to show that performance of the particle branching scheme closely approximates that of the fully synchronous particle filter, even in a small-data small-complexity regime where we expect particle filter performance to be very good. In addition to comparing to a particle filter which resamples synchronously, we also compare to a worst-case particle filter in which we never resample, instead propagating particles forward deterministically with a single child particle at every $n$. While the statistical (per-sample) efficiency of this approach is quite poor, it is fully parallelizable with no blocking operations in the algorithm at all, and thus provides a ceiling estimate of the raw sampling speed attainable in our overall implementation.

We also benchmark against what we believed to be the most practically competitive similar approach, iterated conditional SMC . Iterated conditional SMC corresponds to the particle Gibbs algorithm in the case where parameter values are known; by using a particle filter sweep as a step within a larger MCMC algorithm, iCSMC provides a statistically valid approach to sampling from a posterior distribution by repeatedly running sequential Monte Carlo sweeps each with a fixed number of particles. One downside to iCSMC is that it does not provide an estimate of the marginal likelihood.

On both these models we see the statistical efficiency of the particle cascade is approximately in line with the true particle filter, slightly outperforming the iCSMC algorithm and significantly outperforming the fully parallelized non-resampling approach. This suggests that the approximations made by computing weights at each $n$ based on only the previously observed particles, and the total particle count limit imposed by $\rho$, do not have an adverse effect on overall performance. In Fig. 2 we plot convergence per particle to the true posterior distribution, as well as convergence in our estimate of the normalizing constant.

Although values will be implementation-dependent, we are ultimately interested not in per-sample efficiency but rather in our rate of convergence over time. We record wall clock time for each algorithm for both of these models; the results for convergence of our estimates of values and marginal likelihood are shown in Fig. 3. These particular experiments were all run on Amazon EC2, in an 8-core environment with Intel Xeon E5-2680 v2 processors. The particle cascade provides a much faster and more accurate estimate of the marginal likelihood than the competing methods, in both models. Convergence in estimates of values is quick as well, faster than the iCSMC approach. We note that for very small numbers of particles, running a simple particle filter is faster than the particle cascade, despite the blocking nature of the resampling step. This is due to the overhead incurred by the particle cascade in sending an initial flurry of $\rho_{0}$ particles into the system before we see any particles progress to the end; this initial speed advantage diminishes as the number of samples increases. Furthermore, in stark contrast to the simple SMC method, there are no barriers to drawing more samples from the particle cascade indefinitely. On this fixed hardware environment, our implementation of SMC, which aggressively parallelizes all forward particle simulations, exhibits a dramatic loss of performance as the number of particles increases from $10^{4}$ to $10^{5}$, to the point where simultaneously running $10^{5}$ particles is simply not possible in a feasible amount of time.

We are also interested in how the particle cascade scales up to larger hardware, or down to smaller hardware. A comparison across 5 different hardware configurations is shown in Fig. 4.

Discussion

The particle cascade has broad applicability, appropriate for all SMC and particle filtering inference applications. For example, constructing an appropriate sequence of densities for SMC is possible in arbitrary probabilistic graphical models, including undirected graphical models; see e.g. the sequential decomposition approach of . We are particularly motivated by the SMC-based probabilistic programming systems that have recently appeared in the literature . In both references it was suggested that primary performance bottleneck in their inference algorithms was barrier synchronization, something we have done away with entirely. What is more, while particle MCMC methods are particularly appropriate when there is a clear boundary that can be exploited between between parameters of interest and nuisance state variables, in a growing number of applications in the probabilistic programming and SMC communities, parameters values are generated as part of the state trajectory itself, leaving no explicitly denominated latent parameter variables per se. The particle cascade is particularly relevant to such approaches.

Finally, an attractive property of this algorithm is that it yields an unbiased estimate of the marginal likelihood, and thus can be plugged directly into PIMH, SMC2 , and other so-called pseudomarginal methods.

Yee Whye Teh’s research leading to these results has received funding from EPSRC (grant EP/K009362/1) and the ERC under the EU’s FP7 Programme (grant agreement no. 617411). Frank Wood is supported under DARPA PPAML. This material is based on research sponsored by DARPA through the U.S. Air Force Research Laboratory under Cooperative Agreement number FA8750-14-2-0004. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation heron. The views and conclusions contained herein are those of the authors and should be not interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of DARPA, the U.S. Air Force Research Laboratory of the U.S. Government.

References

References