Lawrence Murray

Parallel Resampling in the Particle Filter

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 parallelize, 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 analyze 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 upward 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. Supplementary materials are available online.

Bayesian learning and predictability in a stochastic nonlinear dynamical model

Bayesian inference methods are applied within a Bayesian hierarchical modeling framework to the problems of joint state and parameter estimation, and of state forecasting. We explore and demonstrate the ideas in the context of a simple nonlinear marine biogeochemical model. A novel approach is proposed to the formulation of the stochastic process model, in which ecophysiological properties of plankton communities are represented by autoregressive stochastic processes. This approach captures the effects of changes in plankton communities over time, and it allows the incorporation of literature metadata on individual species into prior distributions for process model parameters. The approach is applied to a case study at Ocean Station Papa, using particle Markov chain Monte Carlo computational techniques. The results suggest that, by drawing on objective prior information, it is possible to extract useful information about model state and a subset of parameters, and even to make useful long-term forecasts, based on sparse and noisy observations.

Particle Smoothing in Continuous Time: A Fast Approach via Density Estimation

We consider the particle smoothing problem for state-space models where the transition density is not available in closed form, in particular for continuous-time, nonlinear models expressed via stochastic differential equations (SDEs). Conventional forward-backward and two-filter smoothers for the particle filter require a closed-form transition density, with the linear-Gaussian Euler-Maruyama discretization usually applied to the SDEs to achieve this. We develop a pair of variants using kernel density approximations to relieve the dependence, and in doing so enable use of faster and more accurate discretization schemes such as Runge-Kutta. In addition, the new methods admit arbitrary proposal distributions, providing an avenue to deal with degeneracy issues. Experimental results on a functional magnetic resonance imaging (fMRI) deconvolution task demonstrate comparable accuracy and significantly improved runtime over conventional techniques.

On Disturbance State-Space Models and the Particle Marginal Metropolis-Hastings Sampler

We investigate nonlinear state-space models without a closed-form transition density and propose reformulating such models over their latent noise variables rather than their latent state variables. In doing so the tractable noise density emerges in place of the intractable transition density. For importance sampling methods such as the auxiliary particle filter, this enables importance weights to be computed where they could not be otherwise. As case studies we take two multivariate marine biogeochemical models and perform state and parameter estimation using the particle marginal Metropolis--Hastings sampler. For the particle filter within this sampler, we compare several proposal strategies over noise variables, all based on lookaheads with the unscented Kalman filter. These strategies are compared using conventional means for assessing Metropolis--Hastings efficiency, as well as with a novel metric called the conditional acceptance rate for assessing the consequences of using an estimated, and not exa...

Path storage in the particle filter

This article considers the problem of storing the paths generated by a particle filter and more generally by a sequential Monte Carlo algorithm. It provides a theoretical result bounding the expected memory cost by T+CNlogN where T is the time horizon, N is the number of particles and C is a constant, as well as an efficient algorithm to realise this. The theoretical result and the algorithm are illustrated with numerical experiments.

Sequential Monte Carlo with Highly Informative Observations

We propose sequential Monte Carlo (SMC) methods for sampling the posterior distribution of state-space models under highly informative observation regimes, a situation in which standard SMC methods can perform poorly. A special case is simulating bridges between given initial and final values. The basic idea is to introduce a schedule of intermediate weighting and resampling times between observation times, which guide particles towards the final state. This can always be done for continuous-time models, and may be done for discrete-time models under sparse observation regimes; our main focus is on continuous-time diffusion processes. The methods are broadly applicable in that they support multivariate models with partial observation, do not require simulation of the backward transition (which is often unavailable), and, where possible, avoid pointwise evaluation of the forward transition. When simulating bridges, the last characteristic cannot be avoided entirely without concessions, and we suggest an

High-performance pseudo-random number generation on graphics processing units

\...

Feynman-Kac particle integration with geometric interacting jumps

This work considers the deployment of pseudo-random number generators (PRNGs) on graphics processing units (GPUs), developing an approach based on the xorgens generator to rapidly produce pseudo-random numbers of high statistical quality. The chosen algorithm has configurable state size and period, making it ideal for tuning to the GPU architecture. We present a comparison of both speed and statistical quality with other common GPU-based PRNGs, demonstrating favourable performance of the xorgens-based approach.

Applications of heterogeneous computing in computational and simulation science

This article is concerned with the design and analysis of discrete time Feynman-Kac particle integration models with geometric interacting jump processes. We analyze two general types of model, corresponding to whether the reference process is in continuous or discrete time. For the former, we consider discrete generation particle models defined by arbitrarily fine time mesh approximations of the Feynman-Kac models with continuous time path integrals. For the latter, we assume that the discrete process is observed at integer times and we design new approximation models with geometric interacting jumps in terms of a sequence of intermediate time steps between the integers. In both situations, we provide nonasymptotic bias and variance theorems w.r.t. the time step and the size of the system, yielding what appear to be the first results of this type for this class of Feynman-Kac particle integration models. We also discuss uniform convergence estimates w.r.t. the time horizon. Our approach is based on an original semigroup analysis with first order decompositions of the fluctuation errors.

Applications of Heterogeneous Computing in Computational and Simulation Science

As the size and complexity of scientific problems and datasets grow, scientists from a broad range of discipline areas are relying more and more on computational methods and simulations to help solve their problems. This paper presents a summary of heterogeneous algorithms and applications that have been developed by a large research organization (CSIRO) for solving practical and challenging science problems faster than is possible with conventional multi-core CPUs alone. The problem domains discussed include biological image analysis, computed tomography reconstruction, marine biogeochemical models, fluid dynamics, and bioinformatics. The algorithms utilize GPUs and multi-core CPUs on a scale ranging from single workstation installations through to large GPU clusters. Results demonstrate that large GPU clusters can be used to accelerate a variety of practical science applications, and justify the significant financial investment and interest being placed into such systems.