Nested Sequential Monte Carlo Methods

Introduction

Inference in complex and high-dimensional statistical models is a very challenging problem that is ubiquitous in applications. Examples include, but are definitely not limited to, climate informatics (Monteleoni et al., 2013), bioinformatics (Cohen, 2004) and machine learning (Wainwright and Jordan, 2008). In particular, we are interested in sequential Bayesian inference, which involves computing integrals of the form

for some sequence of probability densities

there are local dependencies among the latent variables $X_{1:k}$, both w.r.t. time $k$ and between the individual components of the (high-dimensional) vectors $X_{k}$.

One example of the type of models we consider are the so-called spatio-temporal models (Wikle, 2015; Cressie and Wikle, 2011; Rue and Held, 2005). In Figure 1 we provide a probabilistic graphical model representation of a spatio-temporal model that we will explore further in Section 6.

Sequential Monte Carlo (SMC) methods, reviewed in Section 2.1, comprise one of the most successful methodologies for sequential Bayesian inference. However, SMC struggles in high-dimensions and these methods are rarely used for dimensions, say, $d\geq 10$ (Rebeschini and van Handel, 2015). The purpose of the NSMC methodology is to push this limit well beyond $d=10$.

The basic strategy, described in Section 2.2, is to mimic the behaviour of a so-called fully adapted SMC algorithm. Full adaptation can drastically improve the efficiency of SMC in high dimensions. Unfortunately, it can rarely be implemented in practice since the fully adapted proposal distributions are typically intractable. NSMC addresses this difficulty by requiring only approximate, properly weighted, samples from the proposal distribution. The proper weighting condition ensures the validity of NSMC, thus providing a generalisation of the family of SMC methods. Furthermore, NSMC will itself produce properly weighted samples. Consequently, it is possible to use one NSMC procedure within another to construct efficient high-dimensional proposal distributions. This nesting of the algorithm can be done to an arbitrary degree. For instance, for the model depicted in Figure 1 we could use three nested samplers, one for each dimension of the “volume”.

The main methodological development is concentrated to Sections 3–4. We introduce the concept of proper weighting, approximations of the proposal distribution, and nesting of Monte Carlo algorithms. Throughout Section 3 we consider simple importance sampling and in Section 4 we extend the development to the sequential setting.

We deliberately defer the discussion of the existing body of related work until Section 5, to open up for a better understanding of the relationships to the new developments presented in Sections 3–4. We also discuss various attractive features of NSMC that are of interest in high-dimensional settings, e.g. the fact that it is easy to distribute the computation, which results in improved memory efficiency and lower communication costs. Section 6 profiles our method extensively with a state-of-the-art competing algorithm on several high-dimensional data sets. We also show the performance of inference and the modularity of the method on a $d=1\thinspace 056$ dimensional climatological spatio-temporal model (Fu et al., 2012) structured according to Figure 1. Finally, in Section 7 we conclude the paper with some final remarks.

Background and Inference Strategy

Evaluating $\bar{\pi}_{k}(f)$ as well as the normalisation constant $Z_{\pi_{k}}$ in (2) is typically intractable and we need to resort to approximations. SMC methods, or particle filters (PF), constitute a popular class of numerical approximations for sequential inference problems. Here we give a high-level introduction to the concepts underlying SMC methods, and postpone the details to Section 4. For a more extensive treatment we refer to Doucet and Johansen (2011); Cappé et al. (2005); Doucet et al. (2001). In particular, we will use the auxiliary SMC method as proposed by Pitt and Shephard (1999).

At iteration $k-1$, the SMC sampler approximates the target distribution $\bar{\pi}_{k-1}$ by a collection of weighted particles (samples) $\{(X_{1:k-1}^{i},W_{k-1}^{i})\}_{i=1}^{N}$. These samples define an empirical point-mass approximation of the target distribution

The resampling step puts emphasis on the most promising particles by discarding the unlikely ones and duplicating the likely ones. The propagation and weighting steps essentially correspond to using importance sampling when changing the target distribution from $\bar{\pi}_{k-1}$ to $\bar{\pi}_{k}$, i.e. simulating new particles from a proposal distribution and then computing corresponding importance weights.

2 Adapting the Proposal Distribution

The first working SMC algorithm was the bootstrap PF by Gordon et al. (1993), which propagates particles by sampling from the system dynamics and computes importance weights according to the observation likelihood (in the state space setting). However, it is well known that the bootstrap PF suffers from weight collapse in high-dimensional settings (Bickel et al., 2008), i.e. the estimate is dominated by a single particle with weight close to one. This is an effect of the mismatch between the importance sampling proposal and the target distribution, which typically gets more pronounced in high dimensions.

More efficient proposals, partially alleviating the degeneracy issue for some models, can be designed by adapting the proposal distribution to the target distribution (see Section 4.2). In Naesseth et al. (2014a) we make use of the fully adapted SMC method (Pitt and Shephard, 1999) for doing inference in a (fairly) high-dimensional discrete model where $x_{k}$ is a $60$-dimensional discrete vector. We can then make use of forward filtering and backward simulation, operating on the individual components of each $x_{k}$, in order to sample from the fully adapted SMC proposals. However, this method is limited to models where the latent space is either discrete or Gaussian and the optimal proposal can be identified with a tree-structured graphical model. Our development here can be seen as a non-trivial extension of this technique. Instead of coupling one SMC sampler with an exact forward filter/backward simulator (which in fact reduces to an instance of standard SMC), we derive a way of coupling multiple SMC samplers and SMC-based backward simulators. This allows us to construct procedures for mimicking the efficient fully adapted proposals for arbitrary latent spaces and structures in high-dimensional models.

Proper Weighting and Nested Importance Sampling

In this section we will lay the groundwork for the derivation of the class of NSMC algorithms. We start by considering the simpler case of importance sampling (IS), which is a fundamental component of SMC, and introduce the key concepts that we make use of. In particular, we will use a (slightly nonstandard) presentation of an algorithm as an instance of a class, in the object-oriented sense, and show that these classes can be nested to an arbitrary degree.

Interestingly, to construct a valid IS algorithm for our target $\bar{\pi}$ it is sufficient to generate samples that are properly weighted w.r.t. the proposal distribution $q$. To formalise this claim, assume that we are not able to simulate exactly from $\bar{q}$, but that it is possible to evaluate the unnormalised density $q$ point-wise. Furthermore, assume we have access to a class $\mathsf{Q}$, which works as follows. The constructor of $\mathsf{Q}$ requires the specification of an unnormalised density function, say, $q$, which will be approximated by the procedures of $\mathsf{Q}$. Furthermore, to highlight the fact that we will typically use IS (and SMC) to construct $\mathsf{Q}$, the constructor also takes as an argument a precision parameter $M$, corresponding to the number of samples used by the “internal” Monte Carlo procedure. An object is then instantiated as $\mathsf{q}=\mathsf{Q}(q,M)$. The class $\mathsf{Q}$ is assumed to have the following properties:

Let $\mathsf{q}=\mathsf{Q}(q,M)$. Assume that:

The construction of $\mathsf{q}$ results in the generation of a (possibly random) member variable, accessible as $\widehat{Z}_{q}=\mathsf{q}.\mathsf{GetZ}()$. The variable $\widehat{Z}_{q}$ is a nonnegative, unbiased estimate of the normalising constant $Z_{q}=\int q(x)dx$.

$\mathsf{Q}$ has a member function $\mathsf{Simulate}$ which returns a (possibly random) variable $X=\mathsf{q}.\mathsf{Simulate}()$, such that $(X,\widehat{Z}_{q})$ is properly weighted for $q$.

With the definition of $\mathsf{Q}$ in place, it is possible to generaliseWith $\mathsf{q}.\mathsf{GetZ}()\mapsto Z$ and $\mathsf{q}.\mathsf{Simulate}()$ returning a sample from $\bar{q}$ we obtain the standard IS method. the basic importance sampler as in Algorithm 1, which generates weighted samples $\{(X^{i},W^{i})\}_{i=1}^{N}$ targeting $\bar{\pi}$. Note that Algorithm 1 is different from a random weight IS, since it approximates the proposal distribution (and not just the importance weights).

To see the validity of Algorithm 1 we can interpret the sampler as a standard IS algorithm for an extended target distribution, defined as $\bar{\Pi}(x,u):=u\,\bar{Q}(x,u)\bar{\pi}(x)q^{-1}(x)$, where $\bar{Q}(x,u)$ is the joint PDF of the random pair $(\mathsf{q}.\mathsf{Simulate}(),\mathsf{q}.\mathsf{GetZ}())$. Note that $\bar{\Pi}$ is indeed a PDF that admits $\bar{\pi}$ as a marginal; for any measurable subset $A\subseteq\mathsf{X}$,

where the penultimate equality follows from the fact that $(X,\widehat{Z}_{q})$ is properly weighted for $q$. Furthermore, the standard unnormalised IS weight for a sampler with target $\bar{\Pi}$ and proposal $\bar{Q}$ is given by $u\,\pi/q$, in agreement with Algorithm 1.

Algorithm 1 is an example of what is referred to as an exact approximation; see e.g., Andrieu and Roberts (2009); Andrieu et al. (2010). Algorithmically, the method appears to be an approximation of an IS, but samples generated by the algorithm nevertheless target the correct distribution $\bar{\pi}$.

2 Modularity of Nested IS

To be able to implement Algorithm 1 we need to define a class $\mathsf{Q}$ with the required properties (A3.1). The modularity of the procedure (as well as its name) comes from the fact that we can use Algorithm 1 also in this respect. Indeed, let us now view $\bar{\pi}$—the target distribution of Algorithm 1—as the proposal distribution for another Nested IS procedure and consider the following definition of $\mathsf{Q}$:

Algorithm 1 is executed at the construction of the object $\mathsf{p}=\mathsf{Q}(\pi,N)$, and $\mathsf{p}.\mathsf{GetZ}()$ returns the normalising constant estimate $\widehat{Z}_{\pi}$.

where, again, we use the fact that $(X^{i},\widehat{Z}_{q}^{i})$ is properly weighted for $q$. This implies that $(X^{B},\widehat{Z}_{\pi})$ is properly weighted for $\pi$ and that our definition of $\mathsf{Q}(\pi,N)$ indeed satisfies condition (A3.1).

The Nested IS algorithm in itself is unlikely to be of direct practical interest. However, in the next section we will, essentially, repeat the preceding derivation in the context of SMC to develop the NSMC method.

Nested Sequential Monte Carlo

Let us return to the sequential inference problem. As before, let $\bar{\pi}_{k}(x_{1:k})=Z_{\pi_{k}}^{-1}\pi_{k}(x_{1:k})$ denote the target distribution at “time” $k$. The unnormalised density $\pi_{k}$ can be evaluated point-wise, but the normalising constant $Z_{\pi_{k}}$ is typically unknown. We will use SMC to simulate sequentially from the distributions $\{\bar{\pi}_{k}\}_{k=1}^{n}$. In particular, we consider the fully adapted SMC sampler (Pitt and Shephard, 1999), which corresponds to a specific choice of resampling weights and proposal distribution, chosen in such a way that the importance weights are all equal to $1/N$. Specifically, the proposal distribution (often referred to as the optimal proposal) is given by $\bar{q}_{k}(x_{k}\,|\,x_{1:k-1})=Z_{q_{k}}(x_{1:k-1})^{-1}q_{k}(x_{k}\,|\,x_{1:k-1})$, where

As mentioned above, at each iteration $k=1,\,\dots,\,n$, the method produces unweighted samples $\{X_{k}^{i}\}_{i=1}^{N}$ approximating $\bar{\pi}_{k}$. It also produces an unbiased estimate $\widehat{Z}_{\pi_{k}}$ of $Z_{\pi_{k}}$ (Del Moral, 2004, Proposition 7.4.1). The algorithm is expressed in a slightly non-standard form; at iteration $k$ we loop over the ancestor particles, i.e. the particles after resampling at iteration $k-1$, and let each ancestor particle $j$ generate $m_{k}^{j}$ offsprings. (The variable $L$ is just for bookkeeping.) This is done to clarify the connection with the NSMC procedure below. Furthermore, we have included a (completely superfluous) resampling step at iteration $k=1$, where the “dummy variables” $\{X_{1:0}^{i}\}_{i=1}^{N}$ are resampled according to the (all equal) weights $\{Z_{q_{1}}(X_{1:0}^{i})\}_{i=1}^{N}=\{Z_{\pi_{1}}\}_{i=1}^{N}$. The analogue of this step is, however, used in the NSMC algorithm, where the initial normalising constant $Z_{\pi_{1}}$ is estimated. We thus have to resample the corresponding initial particle systems accordingly.

2 Fully Adapted Nested SMC Samplers

In analogue with Section 3, assume now that we are not able to simulate exactly from $\bar{q}_{k}$, nor compute $Z_{q_{k}}$. Instead, we have access to a class $\mathsf{Q}$ which satisfies condition (A3.1). The proposed NSMC method is then given by Algorithm 3.

Algorithm 3 can be seen as an exact approximation of the fully adapted SMC sampler in Algorithm 2. (In Appendix A.1 we provide a formulation of NSMC with arbitrary proposals and resampling weights.) We replace the exact computation of $Z_{q_{k}}$ and exact simulation from $\bar{q}_{k}$, by the approximate procedures available through $\mathsf{Q}$. Despite this approximation, however, Algorithm 3 is a valid SMC method. This is formalised by the following theorem.

where $\{X_{1:k}^{i}\}_{i=1}^{M}$ are generated by Algorithm 3 and $\stackrel{{\scriptstyle\textrm{D}}}{{\longrightarrow}}$ denotes convergence in distribution.

The key point with Theorem 2 is that, under certain regularity conditions, the NSMC method converges at rate $\sqrt{N}$ even for a fixed (and finite) value of the precision parameter $M$. The asymptotic variance $\Sigma_{k}^{M}(f)$, however, will depend on the accuracy and properties of the approximative procedures of $\mathsf{Q}$. We leave it as future work to establish more informative results, relating the asymptotic variance of NSMC to that of the ideal, fully adapted SMC sampler.

3 Backward Simulation and Modularity of NSMC

As previously mentioned, the NSMC procedure is modular in the sense that we can make use of Algorithm 3 also to define the class $\mathsf{Q}$. Thus, we now view $\bar{\pi}_{n}$ as the proposal distribution that we wish to approximately sample from using NSMC. Algorithm 3 directly generates an estimate $\widehat{Z}_{\pi_{n}}$ of the normalising constant of $\pi_{n}$ (which indeed is unbiased, see Theorem 6). However, we also need to generate a sample $\widetilde{X}_{1:n}$ such that $(\widetilde{X}_{1:n},\widehat{Z}_{\pi_{n}})$ is properly weighted for $\pi_{n}$.

The simplest approach, akin to the Nested IS procedure described in Section 3.2, is to draw $B_{n}$ uniformly on $\{1,\,\dots,\,N\}$ and return $\widetilde{X}_{1:n}=X_{1:n}^{B_{n}}$. This will indeed result in a valid definition of the $\mathsf{Simulate}$ procedure. However, this approach will suffer from the well known path degeneracy of SMC samplers. In particular, since we call $\mathsf{q}^{j}.\mathsf{Simulate}()$ multiple times in Step 2(f)i of Algorithm 3, we risk to obtain (very) strongly correlated samples by this simple approach.

It is possible to improve the performance of the above procedure by instead making use of a backward simulator (Godsill et al., 2004; Lindsten and Schön, 2013) to simulate $\widetilde{X}_{1:n}$. The backward simulator, given in Algorithm 4, is a type of smoothing algorithm; it makes use of the particles generated by a forward pass of Algorithm 3 to simulate backward in “time” a trajectory $\widetilde{X}_{1:n}$ approximately distributed according to $\bar{\pi}_{n}$.

Algorithm 4 assumes unweighted particles and can thus be used in conjunction with the fully adapted NSMC procedure of Algorithm 2. If, however, the forward filter is not fully adapted the weights need to be accounted for in the backward simulation; see Appendix A.1.3.

The modularity of NSMC is established by the following result.

Let $\mathsf{p}=\mathsf{Q}(\pi_{n},N)$ be defined as follows:

The constructor executes Algorithm 3 with target distribution $\pi_{n}$ and with $N$ particles, and $\mathsf{p}.\mathsf{GetZ}()$ returns the estimate of the normalising constant $\widehat{Z}_{\pi_{n}}$.

$\mathsf{p}.\mathsf{Simulate}()$ executes Algorithm 4 and returns $\widetilde{X}_{1:n}$.

The class $\mathsf{Q}$ defined as in Definition 5 satisfies condition (A3.1).

Proof See Appendix A.1.3. A direct, and important, consequence of Theorem 6 is that NSMC can be used as a component of powerful learning algorithms, such as the particle Markov chain Monte Carlo (PMCMC) method (Andrieu et al., 2010) and many of the other methods discussed in Section 5. Since standard SMC is a special case of NSMC, Theorem 6 implies proper weighting also of SMC.

Practicalities and Related Work

There has been much recent interest in using SMC within SMC in various ways. The SMC2 by Chopin et al. (2013) and the recent method by Crisan and Míguez (2013) are sequential learning algorithms for state space models, where one SMC sampler for the parameters is coupled with another SMC sampler for the latent states. Johansen et al. (2012) and Chen et al. (2011) address the state inference problem by splitting the state variable into different components and run coupled SMC samplers for these components. These methods differ substantially from NSMC; they solve different problems and the “internal” SMC sampler(s) is constructed in a different way (for approximate marginalisation instead of for approximate simulation). Another related method is the random weights PF of Fearnhead et al. (2010a), requiring exact samples from $\bar{q}$ and where the importance weights are estimated using a nested Monte Carlo algorithm.

The method most closely related to NSMC is the space-time particle filter (ST-PF) (Beskos et al., 2014a), which has been developed independently and in parallel with our work. The ST-PF is also designed for solving inference problems in high-dimensional models. It can be seen as a island PF (Vergé et al., 2015) implementation of the method presented by Naesseth et al. (2014b). Specifically, for a spatio-temporal models they run an island PF over both spatial and temporal dimensions. However, the ST-PF does not generate an approximation of the fully adapted SMC sampler.

Another key distinction between NSMC and ST-PF is that in the latter each particle in the “outer” SMC sampler comprises a complete particle system from the “inner” SMC sampler. For NSMC, on the other hand, the particles will simply correspond to different hypotheses about the latent variables (as in standard SMC), regardless of how many samplers that are nested. This is a key feature of NSMC, since it implies that it is easily distributed over the particles. The main computational effort of Algorithm 3 is the construction of $\{\mathsf{q}^{j}\}_{j=1}^{N}$ and the calls to the $\mathsf{Simulate}$ procedure, which can be done independently for each particle. This leads to improved memory efficiency and lower communication costs. Furthermore, we have found (see Section 6) that NSMC can outperform ST-PF even when run on a single machine with matched computational costs.

Another strength of NSMC methods are their relative ease of implementation, which we show in Section 6.3. We use the framework to sample from what is essentially a cubic grid Markov random field (MRF) model just by implementing three nested samplers, each with a target distribution defined on a simple chain.

There are also other SMC-based methods designed for high-dimensional problems, e.g., the block PF studied by Rebeschini and van Handel (2015), the location particle smoother by Briggs et al. (2013) and the PF-based methods reviewed in Djuric and Bugallo (2013). However, these methods are all inconsistent, as they are based on various approximations that result in systematic errors.

The previously mentioned PMCMC (Andrieu et al., 2010) is a related method, where SMC is used as a component of an MCMC algorithm. We make use of a very similar extended space approach to motivate the validity of our algorithm. Note that our proposed algorithm can be used as a component in PMCMC and most of the other algorithms mentioned above, which further increases the scope of models it can handle.

Experimental Results

We illustrate NSMC on three high-dimensional examples, both with real and synthetic data. We compare NSMC with standard (bootstrap) PF and the ST-PF of Beskos et al. (2014a) with equal computational budgets on a single machine (i.e., neglecting the fact that NSMC is more easily distributed). These methods are, to the best of our knowledge, the only other available consistent online methods for full Bayesian inference in general sequential models. For more detailed explanations of the models and additional results, see Appendix A.3Code available at https://github.com/can-cs/nestedsmc.

We start by considering a high-dimensional Gaussian state space model, where we have access to the true solution through belief propagation. The latent variables and measurements $\{X_{1:k},Y_{1:k}\}$, with $\{X_{k},Y_{k}\}=\left\{X_{k,l},Y_{k,l}\right\}_{l=1}^{d}$, are modeled by a $d\times k$ lattice Gaussian MRF. The true data is simulated from a nearly identical state space model (see Appendix A.3.1). We run a 2-level NSMC sampler. The outer level is fully adapted, i.e. the proposal distribution is $q_{k}=p(x_{k}\,|\,x_{k-1},y_{k})$, which thus constitute the target distribution for the inner level. To generate properly weighted samples from $q_{k}$, we use a bootstrap PF operating on the $d$ components of the vector $x_{k}$. Note that we only use bootstrap proposals where the actual sampling takes place, and that the conditional distribution $p(x_{k}\,|\,x_{k-1},y_{k})$ is not explicitly used.

We simulate data from this model for $k=1,\ldots,100$ for different values of $d=\text{dim}(x_{k})\in\{50,100,200\}$. The exact filtering marginals are computed using belief propagation.We compare with both the ST-PF and standard (bootstrap) PF.

The results are evaluated based on the effective sample size (ESS, see e.g. Fearnhead et al. (2010b)) defined as,

where $\widehat{x}_{k,l}$ denote the mean estimates and $\mu_{k,l}$ and $\sigma_{k,l}^{2}$ denote the true mean and variance of $x_{k,l}\,|\,y_{1:k}$ obtained from belief propagation. The expectation in (5) is approximated by averaging over $100$ independent runs of the involved algorithms. The ESS reflects the estimator accuracy, obvious by the definition which is tightly related to the mean-squared-error. Intuitively the ESS corresponds to the equivalent number of i.i.d. samples needed for the same accuracy.

We also consider the effective resample size (ERS, Kong et al. (1994)), which is based on the resampling weights at the top levels in the respective SMC algorithms,

The ERS is an estimate of the effective number of unique particles (or particle systems in the case of ST-PF) available at each resampling step.

We use $N=500$ and $M=2\cdot d$ for NSMC and match the computational time for ST-PF and bootstrap PF. We report the results in Figure 2. The bootstrap PF is omitted from $d=100$, $200$ due to its poor performance already for $d=50$ (which is to be expected). Each dimension $l=1,\ldots,d$ provides us with a value of the ESS, so we present the median (lines) and $15$–$85$% percentiles (shaded regions) in the first row of Figure 2. The ERS is displayed in the second row of Figure 2. Note that ESS gives a better reflection of estimation accuracy than ERS.

We have conducted additional experiments with different model parameters and different choices for $N$ and $M$ (some additional results are given in Appendix A.3.1). Overall the results seem to be in agreement with the ones presented here, however ST-PF seems to be more robust to the trade-off between $N$ and $M$. A rule-of-thumb for NSMC is to generally try to keep $N$ as high as possible, while still maintaining a reasonably large ERS.

2 Non-Gaussian State Space Model

Next, we consider an example with a non-Gaussian SSM, borrowed from Beskos et al. (2014a) where the full details of the model are given. The transition probability $p(x_{k}\,|\,x_{k-1})$ is a localised Gaussian mixture and the measurement probability $p(y_{k}\,|\,x_{k})$ is t-distributed. The model dimension is $d=1\thinspace 024$. Beskos et al. (2014a) report improvements for ST-PF over both the bootstrap PF and the block PF by Rebeschini and van Handel (2015). We use $N=M=100$ for both ST-PF and NSMC (the special structure of this model implies that there is no significant computational overhead from

running backward sampling) and the bootstrap PF is given $N=10\thinspace 000$. In Figure 3 we report the ESS (5), estimated according to Carpenter et al. (1999). The ESS for the bootstrap PF is close to , for ST-PF around 1–2, and for NSMC slightly higher at 7–8. However, we note that all methods perform quite poorly on this model, and to obtain satisfactory results it would be necessary to use more particles.

3 Spatio-Temporal Model – Drought Detection

In this final example we study the problem of detecting droughts based on measured precipitation data (Jones and Harris, 2013) for different locations on earth. We look at the situation in North America during the years $1901$–$1950$ and the Sahel region in Africa during the years $1950$–$2000$. These spatial regions and time frames were chosen since they include two of the most devastating droughts during the last century, the so-called Dust Bowl in the US during the 1930s (Schubert et al., 2004) and the decades long drought in the Sahel region in Africa starting in the 1960s (Foley et al., 2003; Hoerling et al., 2006).

We consider the spatio-temporal model defined by Fu et al. (2012) and compare with the results therein. Each location in a region is modelled to be in either a normal state or in an abnormal state $1$ (drought). Measurements are given by precipitation (in millimeters) for each location and year. At every time instance $k$ our latent structure is described by a rectangular $2$D grid $X_{k}=\{X_{k,i,j}\}_{i=1,j=1}^{I,J}$; in essence this is the model showcased in Figure 1. Fu et al. (2012) considers the problem of finding the maximum aposteriori configuration, using a linear programming relaxation. We will instead compute an approximation of the full posterior filtering distribution $\bar{\pi}_{k}(x_{k})=p(x_{k}\,|\,y_{1:k})$.

The rectangular structure is used to instantiate an NSMC method that on the first level targets the full posterior filtering distribution. To sample from $X_{k}$ we run, on the second level, an NSMC procedure that operates on the “columns” $X_{k,1:I,j}$, $j=1,\,\dots,\,J$. Finally, to sample each column $X_{k,1:I,j}$ we run a third level of SMC, that operates on the individual components $X_{k,i,j}$, $i=1,\,\dots,\,I$, using a bootstrap proposal. The structure of our NSMC method applied to this particular problem is illustrated in Figure 4.

Figure 5 gives the results on the parts of North America that we consider. The first row shows the number of locations where the estimate of $p(x_{k,i,j}=1)$ exceeds $\{0.5,0.7,0.9\}$, for both regions. These results seems to be in agreement with Fu et al. (2012, Figures 3, 6). However, we also receive an approximation of the full posterior and can visualise uncertainty in our estimates, as illustrated by the three different levels of posterior probability for drought. In general, we obtain a rich sample diversity from the posterior distribution. However, for some problematic years the sampler degenerates, with the result that the three credibility levels all coincide. This is also visible in the second row of Figure 5, where we show the posterior estimates $p(x_{k,i,j}\,|\,y_{1:k})$ for the years 1939–1941, overlayed on the regions of interest. For year 1940 the sampler degenerates and only reports 0-1 probabilities for all sites. Naturally, one way to improve the estimates is to run the sampler with a larger number of particles, which has been kept very low in this proof-of-concept.

Conclusions

We have shown that a straightforward NSMC implementation with fairly few particles can attain reasonable approximations to the filtering problem for dimensions in the order of hundreds, or even thousands. This means that NSMC methods takes the SMC framework an important step closer to being viable for high-dimensional statistical inference problems. However, NSMC is not a silver bullet for solving high-dimensional inference problems, and the approximation accuracy will be highly model dependent. Hence, much work remains to be done, for instance on combining NSMC with other techniques for high-dimensional inference such as localisation (Rebeschini and van Handel, 2015) and annealing (Beskos et al., 2014b), in order to solve even more challenging problems.

Acknowledgments

This work was supported by the projects: Learning of complex dynamical systems (Contract number: 637-2014-466) and Probabilistic modeling of dynamical systems (Contract number: 621-2013-5524), both funded by the Swedish Research Council.

A Appendix

In this appendix we start out in Section A.1 by providing a more general formulation of the NSMC method and proofs of the central limit and proper weighting theorems of the main manuscript. We also detail (Section A.2) a straightforward extension of nested IS to a sequential version. We show that a special case of this nested sequential IS turns out to be more or less equivalent to the importance sampling squared algorithm by Tran et al. (2013). This relationship serves as evidence that illustrates that the NSMC framework being more widely applicable than the scope of problems considered in this article. Finally, in Section A.3 we give more details and results on the experiments considered in the main manuscript.

We start by presenting a general formulation of a nested auxiliary SMC sampler in Algorithm 5. In this formulation, $q_{k}(x_{k}\,|\,x_{1:k-1})$ is an arbitrary (unnormalised) proposal, normalised by

Furthermore, the resampling weights are obtain by multiplying the importance weights with the arbitrary adjustment multipliers $\nu_{k-1}(x_{1:k-1},Z_{q_{k}})$, which may depend on both the state sequence $x_{1:k-1}$ and the normalising constant (estimate). The fully adapted NSMC sampler (Algorithm 3 in the main document) is obtained as a special case if we choose

and $\nu_{k-1}(x_{1:k-1},Z_{q_{k}})=Z_{q_{k}}$, in which case the importance weights are indeed given by $W_{k}^{i}\equiv 1$.

The validity of Algorithm 5 can be established by interpreting the algorithm as a standard SMC procedure for a sequence of extended target distributions. If $\widehat{Z}_{q_{k}}$ is computed deterministically, proper weighting (i.e., unbiasedness) ensures that $\widehat{Z}_{q_{k}}=Z_{q_{k}}$ and it is evident that the algorithm reduces to a standard SMC sampler. Hence, we consider the case when the normalising constant estimates $\widehat{Z}_{q_{k}}$ are random.

For $k=1,\,\dots,\,n+1$, let us introduce the random variable $U_{k-1}$ which encodes the complete internal state of the object $\mathsf{q}$ generated by $\mathsf{q}=\mathsf{Q}(q_{k}(\cdot\,|\,x_{1:k-1}),M)$. Let the distribution of $U_{k-1}$ be denoted as $\bar{\psi}_{k-1}^{M}(u_{k-1}\,|\,x_{1:k-1})$. To put Algorithm 5 into a standard (auxiliary) SMC framework, we shall interpret steps 2a–2b of Algorithm 5 as being the last two steps carried out during iteration $k-1$, rather than the first two steps carried out during iteration $k$. This does not alter the algorithm per se, but it results in that the resampling step is conducted first at each iteration, which is typically the case for standard auxiliary SMC formulations.

The estimator of the normalising constant is computable from the internal state of $\mathsf{q}$, so that we can introduce a function $\tau_{k}$ such that $\widehat{Z}_{q_{k}}=\tau_{k}(U_{k-1})$. Furthermore, note that the simulation of $X_{k}$ via $X_{k}=\mathsf{q}.\mathsf{Simulate}()$ is based solely on the internal state $U_{k-1}$, and denote by $\bar{\gamma}_{k}^{M}(x_{k}\,|\,U_{k-1})$ the distribution of $X_{k}$.

Assume that $\mathsf{Q}$ satisfies condition (A3.1) in the main manuscript. Then,

Proof The pair $(X_{k},\tau_{k}(U_{k-1}))$ are properly weighted for $q_{k}$. Hence, for a measurable function $f$,

Since $f$ is arbitrary, the result follows.

We can now define the sequence of (unnormalised) extended target distributions for the Nested SMC sampler as,

and $\Pi_{0}(u_{0})=\bar{\psi}_{0}^{M}(u_{0})$. We write $\Theta_{k}=\mathsf{X}_{k}\times\mathsf{U}_{k}$ for the domain of $\Pi_{k}$.

Assume that $\mathsf{Q}$ satisfies condition (A3.1) in the main manuscript. Then,

where the penultimate equality follows by applying Lemma 7 and the induction hypothesis to the two integrals, respectively.

As a corollary to Lemma 8, it follows that

Consequently, $\Pi_{k}$ is normalised by the same constant $Z_{\pi_{k}}$ as $\pi_{k}$, and by defining $\bar{\Pi}_{k}(x_{1:k},u_{0:k}):=Z_{\pi_{k}}^{-1}\Pi_{k}(x_{1:k},u_{0:k})$ we obtain a probability distribution which admits $\bar{\pi}_{k}$ as a marginal (note that $\bar{\Pi}_{0}=\Pi_{0}$, which is normalised by construction). This implies that we can use $\bar{\Pi}_{k}$ as a proxy for $\bar{\pi}_{k}$ in a Monte Carlo algorithm, i.e., samples drawn from $\bar{\Pi}_{k}$ can be used to compute expectations w.r.t. $\bar{\pi}_{k}$. This is precisely what Algorithm 5 does; it is a standard auxiliary SMC sampler for the (unnormalised) target sequence $\Pi_{k}$, $k=0,\,\dots,\,n$, with adjustment multiplier weights $\nu_{k-1}(x_{1:k-1},\tau_{k}(u_{k-1}))$ and proposal distribution $\bar{\gamma}_{k}^{M}(x_{k}\,|\,u_{k-1})\bar{\psi}_{k}^{M}(u_{k}\,|\,x_{1:k})$. The (standard) weight function for this sampler is thus given by

which is the same as the expression on line 2(f)ii of Algorithm 5.

A.1.2 Central Limit Theorem – Proof of Theorem 2 in the Main Manuscript

Now that we have established that Nested SMC is in fact a standard auxiliary SMC sampler, albeit on an extended state space, we can reuse existing convergence results from the SMC literature; see e.g., Johansen and Doucet (2008); Douc and Moulines (2008); Douc et al. (2009); Chopin (2004) or the extensive textbook by Del Moral (2004).

Here, in order to prove Theorem 2 of the main manuscript, we make use of the result for the auxiliary SMC sampler by Johansen and Doucet (2008), which in turn is based on the central limit theorem by Chopin (2004). The technique used by Johansen and Doucet (2008) is to reinterpret (as detailed below) the auxiliary SMC sampler as a sequential importance sampling and resampling (SISR) particle filter, by introducing the modified (unnormalised) target distribution

The auxiliary SMC sampler described in the previous section can then be viewed as a SISR algorithm for (9). Indeed, if we write $\bar{Q}_{k}^{M}(x_{k},u_{k}\,|\,x_{1:k-1},u_{k-1}):=\bar{\psi}_{k}^{M}(u_{k}\,|\,x_{1:k})\bar{\gamma}_{k}^{M}(x_{k}\,|\,u_{k-1})$ for the joint proposal distribution of $(x_{k},u_{k})$, then the weight function for this SISR sampler is given by

where $W_{k}$ is defined in (8). This weight expression thus accounts for both the importance weights and the adjustment multipliers of the auxiliary SMC sampler formulation.

Since this SISR algorithm does not target $\bar{\Pi}_{k}$ (and thus not $\bar{\pi}_{k}$) directly, we use an additional IS step to compute estimators of expectations w.r.t. to $\bar{\pi}$. The proposal distribution for this IS procedure is given by

Note that we obtain an approximation of (11) after the propagation Step 2(f)i of Algorithm 5, but before the weighting step. The resulting IS weights, for target distribution $\bar{\Pi}_{k}(x_{1:k},u_{0:k})$ and with proposal distribution (11), are given by

which, again, is in agreement with Algorithm 5.

We have now reinterpreted the NSMC algorithm; first as a standard auxiliary SMC sampler, and then further as a standard SISR method. Consequently, we are now in the position of directly applying, e.g., the central limit theorem by Chopin (2004, Theorem 1). The conditions and the statement of the theorem are reproduced here for clarity.

A.1.3 Nested SMC Generates Properly Weighted Samples – Proof of Theorem 6 in the Main Manuscript

In the previous two sections we showed that the NSMC procedure is a valid inference algorithm for $\bar{\pi}_{n}$. Next, we turn our attention to the modularity of the method and the validity of using the algorithm as a component in another NSMC sampler. Let us start by stating a more general version of the backward simulator in Algorithm 6. Clearly, if the forward NSMC procedure is fully adapted $W_{k}^{i}\equiv 1$, Algorithm 6 reduces to the backward simulator stated in the main manuscript.

We will now show that the pair $(\widehat{Z}_{\pi_{n}},\widetilde{X}_{1:n})$ generated by Algorithms 5 and 6 is properly weighted for $\pi_{n}(x_{1:n})$, and thereby prove Theorem 6 in the main manuscript.

The proof is based on the particle Markov chain Monte Carlo (PMCMC) construction (Andrieu et al., 2010). The idea used by Andrieu et al. (2010) was to construct an extended target distribution, incorporating all the random variables generated by an SMC sampler as auxiliary variables. This opened up for using SMC approximations within MCMC in a provably correct way; these seemingly approximate methods simply correspond to standard MCMC samplers for the (nonstandard) extended target distribution. Here we will use the same technique to prove the proper weighing property of the NSMC procedure.

We start by introducing some additional notation for the auxiliary variables of the extended target construction. While Algorithm 5 is expressed using multinomial random variables $m_{k}^{1:N}$ in the resampling step, it is more convenient for the sake of the proof to explicitly introduce the ancestor indices $\{A_{k}^{i}\}_{i=1}^{N}$; see e.g., Andrieu et al. (2010). That is, $A_{k}^{i}$ is a categorical random variable on $\{1,\,\dots,\,N\}$, such that $X_{1:k-1}^{A_{k}^{i}}$ is ancestor particle at iteration $k-1$ of particle $X_{k}^{i}$. The resampling Step 2d of Algorithm 5 can then equivalently be expressed as: simulate independently $\{A_{k}^{i}\}_{i=1}^{N}$ from the categorical distribution with probabilities

Let $\mathbf{X}_{k}=\{X_{k}^{1},\,\dots,\,X_{k}^{N}\}$, $\mathbf{U}_{k}=\{U_{k}^{1},\,\dots,\,U_{k}^{N}\}$, and $\mathbf{A}_{k}=\{A_{k}^{1},\,\dots,\,A_{k}^{N}\}$, denote all the particles, internal states of the proposals, and ancestor indices, respectively, generated at iteration $k$ of the NSMC algorithm. We can then write down the joint distribution of all the random variables generated in executing Algorithm 5 (up to an irrelevant permutation of the particle indices) as,

where we interpret $\widehat{\nu}_{k}^{i}$ and $W_{k}^{i}$ as deterministic functions of $(x_{1:k}^{i},u_{0:k}^{i})$.

Let $B_{n}$ denote a random variable defined on $\{1,\,\dots,\,N\}$. The extended target distribution for PMCMC samplers corresponding to (14) is then given by

where $\widehat{Z}_{\pi_{n}}$ is a deterministic function of $(\mathbf{x}_{1:n},\mathbf{u}_{0:n},\mathbf{a}_{1:n})$. We know from Andrieu et al. (2010) that $\bar{\Phi}$ is a probability distribution which admits $\bar{\Pi}_{n}$ as its marginal distribution for $(X_{1:n}^{b_{n}},U_{0:n}^{b_{n}})$. Consequently, by (7) it follows that the marginal distribution of $X_{1:n}^{b_{n}}$ is $\bar{\pi}_{n}$. For later reference we define recursively $b_{k-1}:=a_{k}^{b_{k}}$ for $k=1,\,\dots,\,n$, the particle indices for the trajectory obtained by tracing backward the genealogy of the $b_{n}$’th particle at iteration $n$.

We now turn our attention to the backward simulator in Algorithm 6. Backward simulation has indeed been used in the context of PMCMC, see e.g. Whiteley (2010); Lindsten and Schön (2013); Lindsten et al. (2014). The strategy used for combining PMCMC with backward simulation is to show that each step of the backward sampler corresponds to a partially collapsed Gibbs sampling step for the extended target distribution $\bar{\Phi}$. This implies that the backward sampler leaves $\bar{\Phi}$ invariant.

We use the same approach here, but we need to be careful in how we apply the existing results, since the PMCMC distribution $\bar{\Phi}$ is defined w.r.t. to $\bar{\Pi}_{n}$, whereas the backward simulator of Algorithm 6 works with the original target distribution $\bar{\pi}_{n}$. Nevertheless, from the proof of Lemma 1 by Lindsten et al. (2014) it follows that we can write the following collapsed conditional distribution of $\bar{\Phi}$ as:

By Lemma 7 we know that each factor of the product (in brackets) on the second line integrates to $1$ over $u_{s-1}$. Hence, plugging (17) into (A.1.3) and integrating over $u_{k:n}^{b_{k:n}}$ yields

which coincides with the expression used to simulate the index $B_{k}$ in Algorithm 6. Hence, simulation of $B_{k}$ indeed corresponds to a partially collapsed Gibbs sampling step for $\bar{\Phi}$ and it will thus leave $\bar{\Phi}$ invariant. (Note that, in comparison with the PMCMC sampler derived by Lindsten et al. (2014) we further marginalise over the variables $u_{k:n}^{b_{k:n}}$ which, however, still results in a valid partially collapsed Gibbs step.)

We now have all the components needed to prove proper weighting of the combined NSMC/backward simulation procedure. For notational simplicity, we write

for the distribution of $B_{k}$ in Algorithm 6. Let $(\widehat{Z}_{\pi_{n}},\widetilde{X}_{1:n})$ be generated by Algorithms 5 and 6. Let $f$ be a measurable function and consider

A.2 Nested Sequential Importance Sampling

Here we give the definition of the nested sequential importance sampler and we show that a special case of this is the importance sampling squared (IS2) method by Tran et al. (2013).

We present a straightforward extension of the Nested IS class to a sequential IS version. Consider the following definition of the Nested SIS $\mathsf{Q}$:

Algorithm 7 is executed at the construction of the object $\mathsf{p}=\mathsf{Q}(\pi_{n},N)$, and $\mathsf{p}.\mathsf{GetZ}()$ returns the normalising constant estimate $\widehat{Z}_{\pi_{n}}$.

Note that we do not require that the procedure $\mathsf{Q}$ is identical for each individual proposal $q_{k}$, thus we have a flexibility in designing our algorithm as can be seen in the example in Section A.2.2. We can motivate the algorithm in the same way as for Nested IS and similar theoretical results hold, i.e. Nested SIS is properly weighted for $\pi_{n}$ and it admits $\bar{\pi}_{n}$ as a marginal.

A.2.2 Relation to IS2

Here we will show how IS2, proposed by Tran et al. (2013), can be viewed as a special case of Nested SIS. We are interested in approximating the posterior distribution of parameters $\theta$ given some observed values $y$

We assume that the data likelihood $p(y\,|\,\theta)$ can, by introducing a latent variable $x$, be computed as an integral

Now, let our target distribution in Nested SIS be $\bar{\pi}_{2}(\theta,x)=\bar{\pi}_{2}(x\,|\,\theta)\bar{\pi}_{1}(\theta)=\frac{p(y\,|\,x,\theta)p(x\,|\,\theta)}{p(y\,|\,\theta)}p(\theta)$. We set our proposal distributions to be

First, $\mathsf{Q}(q_{1}(\cdot),1)$ runs an exact sampler from the proposal $g_{\text{IS}}$. Then at iteration $k=2$ we let the nested procedure $\mathsf{Q}(q_{2}(\cdot\,|\,\theta^{i}),M)$ be a standard IS algorithm with proposal $h(x\,|\,y,\theta)$, giving us properly weighted samples for $q_{2}$. Putting all this together gives us samples $\theta^{i}$ distributed according to $g_{\text{IS}}(\theta)$ and weighted by

A.3 Further Details on the Experiments

We provide some further details and results for the experiments presented in the main manuscript.

We generate data from a synthetic $d$-dimensional (dim$(x_{k})=d$) dynamical/spatio-temporalNote that in a previous version this was erraneously stated as equivivalent to the Gaussian MRF we use for sequential inference. Thus this example actually illustrates a problem where we have a misspecified model. However, this misspecification does not lead to any discernible difference in the MSE results. This because the exact filtering marginals for the two different models (LGSS, GMRF) with the parameters chosen differs with orders of magnitudes much lower than the Monte Carlo errors. model defined by

where $\Sigma$ and $\mu_{k}$ are given as follows

Alternatively, in a more standard state space model notation, we have

where $A=a\tau_{\rho}\Sigma$, $Q=\Sigma$ and $R=\tau_{\phi}^{-1}I$. We assume that the parameters $\theta=(\tau_{\psi},a,\tau_{\rho},\tau_{\phi})=(1,0.5,1,10)$ are known.

To do inference with this generated data-set $\{y_{k}\}$ we propose to target the following slightly different model

where the observation potential $\boldsymbol{\phi}$ and interaction potentials $\boldsymbol{\rho}$ and $\boldsymbol{\psi}$ are given by

This can be visualised as a Gaussian rectangular ($d\times k$) lattice MRF, i.e. it grows with “time” $k$. The goal is to estimate the filtering distribution $p(x_{k}\,|\,y_{1:k})$. Note that this model has almost identical filtering marginals as the data generating distribution and leads to a simpler implementation of NSMC and ST-PF.

Results (mean-squared-error, MSE) comparing NSMC and ST-PF for different settings of $N$ and $M$ can be found in the first row of Figure 6 and the second row displays the results when comparing ST-PF to the SMC method by Naesseth et al. (2014b) for equal computational budgets. We show median (over dimensions $d$) MSE for posterior marginal mean and variance estimates of the respective algorithms. True values are obtained using belief propagation. Note that setting $N=1$ in ST-PF can be viewed as a special case of the SMC method by Naesseth et al. (2014b).

A.3.2 Spatio-Temporal Model – Drought Detection

We present the full model for drought detection in our notation, this is essentially the model by Fu et al. (2012) adapted for estimating the filtering distribution. The latent variables for each location on a finite world grid, $x_{k,i,j}$, are binary, i.e. being normal state and $1$ being the abnormal (drought) state. Measurements, $y_{k,i,j}$, are available as real valued precipitation values in millimeters. The probabilistic model for filtering is given as,

References