Elements of Sequential Monte Carlo

TL;DR

This tutorial provides a comprehensive overview of sequential Monte Carlo methods for approximate inference in probabilistic models, covering theory, practical design choices, and recent advances in learning proposals.

Contribution

It reviews the fundamentals of SMC, discusses recent developments in learning proposals and target distributions, and illustrates applications in machine learning models.

Findings

SMC effectively approximates intractable expectations in Bayesian inference.

Recent methods improve proposal distributions using variational inference.

SMC estimates of normalizing constants enable advanced inference techniques.

Abstract

A core problem in statistics and probabilistic machine learning is to compute probability distributions and expectations. This is the fundamental problem of Bayesian statistics and machine learning, which frames all inference as expectations with respect to the posterior distribution. The key challenge is to approximate these intractable expectations. In this tutorial, we review sequential Monte Carlo (SMC), a random-sampling-based class of methods for approximate inference. First, we explain the basics of SMC, discuss practical issues, and review theoretical results. We then examine two of the main user design choices: the proposal distributions and the so called intermediate target distributions. We review recent results on how variational inference and amortization can be used to learn efficient proposals and target distributions. Next, we discuss the SMC estimate of the normalizing constant, how this can be used for pseudo-marginal inference and inference evaluation. Throughout the tutorial we illustrate the use of SMC on various models commonly used in machine learning, such as stochastic recurrent neural networks, probabilistic graphical models, and probabilistic programs.

Figures

Tables

Table 1. Table 2.1 : Average log-probability values of the unnormalized target distribution with respect to the sampling distributions of sis and smc . The number of particles N = 10 𝑁 10 N=10 is fixed for both methods.

	$\text{mathbb}{E}_{{\hat{\gamma }}_{10}}\left[\log {\tilde{\gamma }}_{10}/10\right]$	$\text{mathbb}{E}_{{\hat{\gamma }}_{20}}\left[\log {\tilde{\gamma }}_{20}/20\right]$	$\text{mathbb}{E}_{{\hat{\gamma }}_{40}}\left[\log {\tilde{\gamma }}_{40}/40\right]$
sis	$-2.76$	$-3.35$	$-9.86$
smc	$-2.47$	$-2.51$	$-2.77$

Equations

$$p(\mathbf{x},\mathbf{y}).$$

$$p(\mathbf{x}|\mathbf{y})=\frac{p(\mathbf{x},\mathbf{y})}{p(\mathbf{y})},$$

$${\gamma }_t(x_{1:t})$$

$$x_t|x_{t-1}$$

$$y_t|x_t$$

$$p(\mathbf{x},\mathbf{y})$$

$${\tilde{\gamma }}_t(x_{1:t})$$

$$x_t|x_{1:t-1}$$

$$y_t|x_{1:t}$$

$$p(\mathbf{x},\mathbf{y})$$

$${\tilde{\gamma }}_t(x_{1:t})$$

$${\tilde{\gamma }}_1(x_1)$$

$${\tilde{\gamma }}_t(x_{1:t})$$

$${\tilde{\gamma }}_T(x_{1:T})$$

$$x_t|x_{1:t-1}$$

$$f_t(x_t|x_{t-1})$$

$$g_t(y_t|x_{1:t})$$

$$f_t(x_t|x_{1:t-1})$$

$$p(\mathbf{x},\mathbf{y})$$

$${\tilde{\gamma }}_t(x_{1:t})$$

$$p(\mathbf{x},\mathbf{y})$$

$$p(x_{t+1}=j\mid x_{1:t})$$

$${\theta }_k$$

$$p(y_t\mid x_t,\{{\theta }_k{\}}_{k\ge 1})$$

$${\gamma }_t(h_t):={\mathbb{E}}_{{\gamma }_t}\left[h_t(x_{1:t})\right].$$

$${\mathbb{E}}_{{\gamma }_t}\left[h_t(x_{1:t})\right]$$

$$w_t^i$$

$${\gamma }_t(x_{1:t})$$

$$Z_t$$

$$\mathbb{E}[{\hat{Z}}_t]=\frac{1}{N}\sum _{i=1}^N\mathbb{E}\left[\frac{{\tilde{\gamma }}_t(x_{1:t}^i)}{q_t(x_{1:t}^i)}\right]=\frac{1}{N}\sum _{i=1}^N\int \frac{{\tilde{\gamma }}_t(x_{1:t}^i)}{q_t(x_{1:t}^i)}{q}_t(x_{1:t}^i)\mathrm{d}{x}_{1:t}^i=Z_t,$$

Full text

Chapter 1 Introduction

A key strategy in machine learning is to break down a problem into smaller and more manageable parts, then process data or unknown variables recursively. Well known examples of this are message passing algorithms for graphical models and annealing for optimization or sampling. Sequential Monte Carlo (smc)is a class of methods that are tailored to solved statistical inference problems recursively. These methods have mostly received attention in the signal processing and statistics communities. With well over two decades of research in smc, they have enabled inference in increasingly complex and challenging models. Recently, there has been an emergent interest in this class of algorithms from the machine learning community. We have seen applications to probabilistic programming (Wood et al., 2014), variational inference (vi) (Maddison et al., 2017; Naesseth et al., 2018; Le et al., 2018), inference evaluation (Grosse et al., 2015; Cusumano-Towner and Mansinghka, 2017), probabilistic graphical models (pgms) (Ihler and McAllester, 2009; Naesseth et al., 2014; Paige and Wood, 2016), Bayesian nonparametrics (Fearnhead, 2004) and many other areas.

We provide a unifying view of the smc methods that have been developed since their conception in the early 1990s (Gordon et al., 1993; Stewart and McCarty, 1992; Kitagawa, 1993). In this introduction we provide relevant background material, introduce a running example, and discuss the use of code snippets throughout the tutorial.

1.1 Historical Background

smc methods are generic tools for performing approximate (statistical) inference, predominantly Bayesian inference. They use a weighted sample set to iteratively approximate the posterior distribution of a probabilistic model. Ever since the dawn of Monte Carlo methods (see e.g. Metropolis and Ulam (1949) for an early discussion), random sample-based approximations have been recognized as powerful tools for inference in complex probabilistic models. Parallel to the development of Markov chain Monte Carlo (mcmc) methods (Metropolis et al., 1953; Hastings, 1970), sequential importance sampling (sis) (Handschin and Mayne, 1969) and sampling/importance resampling (Rubin, 1987) laid the foundations for what would one day become smc.

smc methods where initially known as particle filters (Gordon et al., 1993; Stewart and McCarty, 1992; Kitagawa, 1993). Particle filters where conceived as algorithms for online inference in nonlinear state space models (ssms) (Cappé et al., 2005). Since then there has been a flurry of work applying smc and particle filters to perform approximate inference in ever more complex models. While research in smc initially focused on ssms, we will see that smc can be a powerful tool in a much broader setting.

1.2 Probabilistic Models and Target Distributions

As mentioned above, smc methods were originally developed as an approximate solution to the so called filtering problem, which amounts to online inference in dynamical models. Several overview and tutorial articles focus on particle filters, i.e. the smc algorithms specifically tailored to solve the online filtering problem (Arulampalam et al., 2002; Doucet and Johansen, 2009; Fearnhead and Künsch, 2018). In this tutorial we will take a different view and explain how smc can be used to solve more general “offline” problems. We shall see how this viewpoint opens up for many interesting applications of smc in machine learning that do not fall in the traditional filtering setup, and furthermore how it gives rise to new and interesting design choices. For complementary review and tutorial articles that treat a similar topic see also Del Moral et al. (2006b); Doucet and Lee (2018).

We consider a generic probabilistic model given by a joint probability distribution function (pdf) of latent variables $\mathbf{x}$ and observed data $\mathbf{y}$,

[TABLE]

We focus on Bayesian inference, where the key object is the posterior distribution

[TABLE]

where $p(\mathbf{y})$ is known as the marginal likelihood.

The target distributions are a sequence of probability distributions that we recursively approximate using smc. We define each target distribution $\gamma_{t}(x_{1:t})$ in the sequence as a joint pdf of latent variables $x_{1:t}=(x_{1},\ldots,x_{t})$, where $t=1,\ldots,T$. The pdf is denoted by

[TABLE]

where $\widetilde{\gamma}_{t}$ is a positive integrable function and $Z_{t}$ is the normalization constant, ensuring that $\gamma_{t}$ is indeed a pdf.

We connect the target distributions to the probabilistic model through a requirement on the final target distribution $\gamma_{T}(x_{1:T})$. We enforce the condition that $\gamma_{T}(x_{1:T})$ is either equivalent to the posterior distribution, or that it contains the posterior distribution as a marginal distribution. The intermediate target distributions, i.e. $\gamma_{t}(x_{1:t})$ for $t<T$, are useful only insofar they help us approximate the final target $\gamma_{T}(x_{1:T})$. This approach is distinct from most previous tutorials on particle filters and smc that traditionally focus on the intermediate targets, i.e. the filtering distributions. We stress that it is not necessarily the case that $x_{1:T}$ (the latent variables of the target distribution) is equal to $\mathbf{x}$ (the latent variables of the probabilistic model of interest), as the former can include additional auxiliary variables.

Below we introduce a few examples of probabilistic models and some straightforward choices of target distributions. We introduce and illustrate our running example which will be used throughout. We will return to the issue of choosing the sequence of intermediate targets in Section 3.2.

State Space Models

The state space model (or hidden Markov model) is a type of probabilistic models where the latent variables and data satisfy a Markov property. For this model we typically have $\mathbf{x}=x_{1:T}$. Often the measured data can also be split into a sequence of the same length ($T$) as the latent variables, i.e. $\mathbf{y}=y_{1:T}$. The model is defined by a transition pdf $f$ and an observation pdf $g$,

[TABLE]

The joint pdf is

[TABLE]

where $p(x_{1})$ is the prior on the initial state $x_{1}$. This class of models is especially common for data that has an inherent temporal structure such as in the field of signal processing. A common choice is to let the target distributions follow the same sequential structure as in Eq. 1.5:

[TABLE]

which means that the final normalized target distribution satisfies $\gamma_{T}(x_{1:T})=p(\mathbf{x}\,|\,\mathbf{y})$ as required. This is the model class and target distributions which are studied in the classical filtering setup.

Non-Markovian Latent Variable Models

The non-Markovian latent variable models (lvms) are characterized by either no, or higher order, Markov structure between the latent variables $\mathbf{x}$ and/or data $\mathbf{y}$. This can be seen as a non-trivial extension of the ssm, see Eq. 1.4, which has a Markov structure. Also for this class of models it is common to have $\mathbf{x}=x_{1:T}$ and $\mathbf{y}=y_{1:T}$.

Unlike the ssm, the non-Markovian lvm in its most general setting requires access to all previous latent variables $x_{1:t-1}$ to generate $x_{t},y_{t}$

[TABLE]

where we again refer to $f_{t}$ and $g_{t}$ as the transition pdf and observation pdf, respectively. The joint pdf is given by

[TABLE]

where $p(x_{1})$ is again the prior on $x_{1}$. A typical target distribution is given by

[TABLE]

with $\widetilde{\gamma}_{1}(x_{1})=p(x_{1})g_{1}(y_{1}\,|\,x_{1})$. Another option is

[TABLE]

For both these sequences of target distributions the final iteration $T$ is the posterior distribution, i.e. $\gamma_{T}(x_{1:T})=p(x_{1:T}\,|\,y_{1:T})=p(\mathbf{x}\,|\,\mathbf{y})$. However, the former one will often lead to more accurate inferences. This is because we introduce information from the data at an earlier stage in the smc algorithm.

Throughout the monograph we will exemplify the different methods using a Gaussian special case of Eq. 1.7, see Example 1.2.1. We let the prior on $x_{1:t}$, defined by the transition pdfs $f_{1},\ldots f_{t}$, be Markovian and introduce the non-Markov property instead through the observation pdfs $g_{1},\ldots,g_{t}$.

Example 1.2.1** (Non-Markovian Gaussian Sequence Model).**

As our running example for illustration purposes we use a non-Markovian Gaussian sequence model. It is

[TABLE]

with observed variables $y_{t}$ (data), and where

[TABLE]

where $\mathcal{N}\left(x\,|\,\mu,\Sigma\right)$ denotes a Gaussian distribution on $x$ with mean $\mu$ and (co)variance $\Sigma$. We let the prior at $t=1$ be $p(x_{1})=\mathcal{N}(x_{1}\,|\,0,q)$. Artificial data was generated using $(\phi,q,\beta,r)=(0.9,1,0.5,1)$. The distribution of interest is the posterior distribution $p(x_{1:T}\,|\,y_{1:T})$. We illustrate a few sample paths of $y_{1:T}$ in Fig. 1.1 for $T=100$.

We can adjust the strength of the dependence on previous latent variables in the observations, $y_{t}$, through the parameter $\beta\in[0,1]$. If we set $\beta=0$ we obtain a linear Gaussian ssm, since the data depends only on the most recent latent $x_{t}$. On the other hand if we let $\beta=1$, this signifies that $x_{k}$ for $k<t$ has equally strong effect on $y_{t}$ as does $x_{t}$.

Another example of a non-Markovian lvm often encountered in machine learning is the stochastic recurrent neural network (rnn) (Chung et al., 2015; Bayer and Osendorfer, 2014; Fraccaro et al., 2016; Maddison et al., 2017; Naesseth et al., 2018). We define the stochastic rnn below and will return to it again in Chapter 3.

Example 1.2.2** (Stochastic Recurrent Neural Network).**

A stochastic rnn is a non-Markovian lvm where the parameters of the transition and observation models are defined by rnn s. A common example is using the conditional Gaussian distribution to define the transition pdf

[TABLE]

where the functions $\mu_{t}(\cdot),\Sigma_{t}(\cdot)$ are defined by rnn s.

Conditionally independent models

A common model in probabilistic machine learning is to assume that the datapoints $y_{k}$ in the dataset $\mathbf{y}=\{y_{k}\}_{k=1}^{K}$ are conditionally independent given the latent $\mathbf{x}$. This means that the joint pdf is given by

[TABLE]

where $p(\mathbf{y}\,|\,\mathbf{x})$ is the likelihood. For this class of models it might not be immediately apparent that we can define a useful sequence of target distributions on latent variables $x_{1:t}$. There is no obvious sequence structure as in the ssm or the non-Markovian lvm. However, as we shall see, we can make use of auxiliary variables to design target distributions that can help with inference.

We will discuss two approaches to design the sequence of target distributions: using data tempering and likelihood tempering, respectively. Both of these will make use of an auxiliary variable technique, where each $x_{t}$ is a random variable on the same space as $\mathbf{x}$.

Data tempering: Using data tempering we add the data $y_{k}$ to the target distribution one by one. In this case the data index $k$ coincides with the target index $t$. We define the target distribution

[TABLE]

where the distributions $s_{k}(x_{k}\,|\,x_{k+1})$ are a design choice, known as backward kernels (Chopin, 2002; Del Moral et al., 2006a). With this choice, we have that the marginal distribution of $x_{T}$ at the final iteration is exactly the posterior distribution, i.e. $\gamma_{T}(x_{T})=p(\mathbf{x}\,|\,\mathbf{y})$. In fact, at each step we have that the marginal target distribution is a partial posterior $\gamma_{t}(x_{t})=p(\mathbf{x}\,|\,y_{1:t})$.

Likelihood tempering: With likelihood tempering, instead of adding the data one by one, we change the likelihood $p(\mathbf{y}\,|\,\mathbf{x})$ through a sequence of positive variables. We define the target distribution

[TABLE]

where $0=\tau_{1}<\ldots<\tau_{T}=1$, and again make use of user chosen backward kernels $s_{k}(x_{k}\,|\,x_{k+1})$ (Chopin, 2002; Del Moral et al., 2006a). In this setting all data is considered at each iteration. Since $\tau_{T}=1$, we have that the final marginal target distribution is again equal to the posterior $\gamma_{T}(x_{T})=p(\mathbf{x}\,|\,\mathbf{y})$.

Applying smc methods to tempered (and similar) target distributions has been studied by e.g. Chopin (2002); Del Moral et al. (2006a). If the proposal $q_{k+1}(x_{k+1}\,|\,x_{k})$ is a Markov kernel with stationary distribution $\gamma_{k+1}(x_{k+1})$, then a commonly used backward kernel is $s_{k}(x_{k}\,|\,x_{k+1})=\frac{\gamma_{k+1}(x_{k})q_{k+1}(x_{k+1}\,|\,x_{k})}{\gamma_{k+1}(x_{k+1})}$. We refer to the works of Chopin (2002); Del Moral et al. (2006a) for a thorough discussion on the choice of backward kernels $s_{k}(x_{k}\,|\,x_{k+1})$. Another well known example is annealed importance sampling by Neal (2001), which uses the backward kernel example above to define the target distributions, but relies on sis instead of smc for inference.

Models and Targets

We have seen several probabilistic models with examples of corresponding target distributions. While not limited to these, this illustrates the wide range of the applicability of smc. In fact, as long as we can design a sequence of target distributions such that $\gamma_{T}$ coincides with the distribution of interest, we can leverage smc for inference.

1.3 Applications

Sequential Monte Carlo and importance sampling methods have already seen a plethora of applications to machine learning and statistical inference problems. Before we turn to the fundamentals of the various algorithms it can be helpful to understand some of these applications. We present and discuss a few select examples of applications of smc and importance sampling (is) to probabilistic graphical models, Bayesian nonparametric models, probabilistic programming, and inference evaluation.

Probabilistic Graphical Models

Probabilistic graphical models (pgm; see e.g. Koller et al. (2009); Wainwright and Jordan (2008)) are probabilistic models where the conditional independencies in the joint pdf are described by edges in a graph. The graph structure allows for easy and strong control on the type of prior information that the user can express. The main limitation of the pgm is that exact inference is often intractable and approximate inference is challenging.

The pgm is a probabilistic model where the pdf factorizes according to an underlying graph described by a set of cliques $C\in\mathcal{C}$, i.e. fully connected subsets of the vertices $V\in\mathcal{V}$ where $\mathcal{V}$ contains all individual components of $\mathbf{x}$ and $\mathbf{y}$. The undirected graphical model can be denoted by

[TABLE]

where $x_{C}$ includes all elements of $\mathbf{x}$ and $\mathbf{y}$ in the clique $C$, and $Z$ is a normalization constant ensuring that the right hand side is a proper pdf.

Smc methods have recently been successfully applied to the task of inference in general pgms, see e.g. MacEachern et al. (1999); Chopin (2002); Del Moral et al. (2006a); Ihler and McAllester (2009); Naesseth et al. (2014); Paige and Wood (2016); Lindsten et al. (2017, 2018) for representative examples.

Probabilistic Programming

Probabilistic programming languages are programming languages designed to describe probabilistic models and to automate the process of doing inference in those models. We can think of probabilistic programming as that of automating statistical inference, particularly Bayesian inference, using tools from computer science and computational statistics. Developing a syntax and semantics to describe and denote probabilistic models and the inference problems we are interested in solving is key to the probabilistic programming language. To define what separates a probabilistic programming language from a standard programming language we quote Gordon et al. (2014): “Probabilistic programs are usual functional or imperative programs with two added constructs: (1) the ability to draw values at random from distributions, and (2) the ability to condition values of variables in a program via observations.” This aligns very well with the notion of Bayesian inference through the posterior distribution Eq. 1.2; through the syntax of the language we define $\mathbf{x}$ to be the random values we sample and $\mathbf{y}$ our observations that we condition on through the use of Bayes rule. The probabilistic program then essentially defines our joint probabilistic model $p(\mathbf{x},\mathbf{y})$.

One of the main challenges of probabilistic programming is to develop algorithms that are general enough to enable inference for any model (probabilistic program) that we could conceivably write down using the language. Recently Wood et al. (2014) have shown that smc-based approaches can be used as inference back-ends in probabilistic programs.

For a more thorough treatment of probabilistic programming we refer the interested reader to the recent tutorial by van de Meent et al. (2018) and the survey by Gordon et al. (2014).

Bayesian nonparametric models

Nonparametric models are characterized by having a complexity which grows with the amount of available data. In a Bayesian context this implies that the usual latent random variables (i.e., parameters) of the model are replaced by latent stochastic processes. Examples include Gaussian processes, Dirichlet processes, and Beta processes; see e.g. Hjort et al. (2010) for a general introduction.

Sampling from these latent stochastic processes, conditionally on observed data, can be done using smc. To give a concrete example, consider the Dirichlet process mixture model, which is a clustering model that can handle an unknown and conceptually infinite number of mixture components. Let $y_{t}$, $t=1,2,\ldots$ be a stream of data. Let $x_{t}$, $t=1,2,\ldots$ be a sequence of latent integer-valued random variables, such that $x_{t}$ is the index of the mixture component to which datapoint $y_{t}$ belongs. A generative representation of the mixture assignment variables is given by

[TABLE]

where $J_{t}:=\max\{x_{1:t}\}$ is the number of distinct mixture components represented by the first $t$ datapoints, and $n_{t,j}:=\sum_{k=1}^{t}\mathbb{I}\{x_{k}=j\}$ is the number of datapoints among $y_{1:t}$ that belong to the $j$th component.

The model is completed by specifying the distribution of the data, conditionally on the mixture assignment variables:

[TABLE]

where $G(\;\cdot\mid\theta)$ is an emission probability distribution parameterized by $\theta$ and $F(\theta)$ is a prior distribution over the mixture parameters $\theta$.

Note that the mixture assignment variables $x_{t}$, $t=1,2,\ldots$ evolve according to a latent stochastic process. Solving the clustering problem amounts to computing the posterior distribution of this stochastic process, conditionally on the observed data. One way to address this problem is to use smc; see Fearnhead (2004) for an efficient implementation tailored to the discrete nature of the problem.

Inference Evaluation

An important problem when performing approximate Bayesian inference is to figure out when our approximation is “good enough”? Is it possible to give practical guarantees on the approximation we obtain? We need ways to evaluate how accurate our approximate inference algorithms are when compared to the true target distribution that we are trying to approximate. We will refer to the process of evaluating and validating approximate inference methods as inference evaluation.

Inference evaluation is mainly concerned with measuring how close our approximation is to the true object we are trying to estimate, often a posterior distribution. For simulated data, Grosse et al. (2015, 2016) have shown that we can make use of smc and is to bound the symmetrized Kullback-Leibler (kl) divergence between our approximate posterior and the true posterior. In another related work Cusumano-Towner and Mansinghka (2017) have shown that smc-based methods show promise in estimating the symmetric kl divergence between the approximate posterior and a gold standard algorithm.

1.4 Example Code

We will be making use of inline Python code snippets throughout the manuscript to illustrate the algorithms and methods. Below we summarize the modules that are necessary to import to run the code snippets:

1.5 Outline

The remainder of this tutorial is organized as follows. In Chapter 2, we first introduce is, a foundational building block for smc. Then, we discuss the limitations of is and how smc resolves these. Finally, the section concludes with discussing some practical issues and theoretical results relevant to smc methods.

Chapter 3 is focused on the two key design choices of smc: the proposal and target distributions. Initially we focus on the proposal, discussing various ways of adapting and learning good proposals that will make the approximation more accurate. Then we discuss the sequence of target distributions; how we can learn intermediate distributions that help us when we try to approximate the posterior.

Chapter 4 focuses on pseudo marginal (pm) methods and other smc methods that rely on a concept known as proper weights. First, we provide a simple and straightforward proof of the unbiasedness property of the smc normalization constant estimate. Then, we describe and illustrate the combination of mcmc and smc methods through pm algorithms. We move on to detail properly weighted smc, a concept that unites and extends the random weights and nested smc algorithms. Finally, we conclude the chapter by considering a few approaches for distributed and parallel smc.

In Chapter 5 we introduce conditional sequential Monte Carlo (csmc), and related methods for simulating from and computing expectations with respect to a target distribution. First, we introduce the basic csmc algorithm and provide a straightforward proof of the unbiasedness of the inverse normalization constant estimate. Then, we show how smc and csmc can be combined to leverage multi-core and parallel computing architectures in the interacting particle Markov chain Monte Carlo (ipmcmc) algorithm. Finally, we discuss recent applications of csmc for evaluating approximate inference methods.

The tutorial concludes with a discussion and outlook in Chapter 6.

Chapter 2 Importance Sampling to Sequential Monte Carlo

Typical applications require us to be able to evaluate or sample from the target distributions $\gamma_{t}$, as well as compute their normalization constants $Z_{t}$. For most models and targets this will be intractable, and we need approximations based on e.g. Monte Carlo methods.

In this section, we first review is and some of its shortcomings. Then, we introduce the smc method, the key algorithm underpinning this monograph. Finally, we discuss some key theoretical properties of the smc algorithm.

2.1 Importance Sampling

Importance sampling (Kahn, 1950a, b; Robert and Casella, 2004) is a Monte Carlo method that constructs an approximation using samples from a proposal distribution, and corrects for the discrepancy between the target and proposal using (importance) weights.

Most applications of Bayesian inference can be formulated as computing the expectation of some generic function $h_{t}$, referred to as a test function, with respect to the target distribution $\gamma_{t}$,

[TABLE]

Examples include posterior predictive distributions, Bayesian p-values, and point estimates such as the posterior mean. Computing Eq. 2.1 is often intractable, but by a clever trick we can rewrite it as

[TABLE]

The pdf $q_{t}$ is a user chosen proposal distribution, we assume it is simple to sample from and evaluate. We can now estimate the right hand side of Eq. 2.2 using the Monte Carlo method,

[TABLE]

where $\widetilde{w}_{t}(x_{1:t}):=\nicefrac{{\widetilde{\gamma}_{t}(x_{1:t})}}{{q_{t}(x_{1:t})}}$ and $x_{1:t}^{i}$ are simulated iid from $q_{t}$. We will usually write Eq. 2.3 more compactly as

[TABLE]

where the normalized weights $w_{t}^{i}$ are defined by

[TABLE]

with $\widetilde{w}_{t}^{i}$ a shorthand for $\widetilde{w}_{t}(x_{1:t}^{i})$. The estimate in Eq. 2.4 is strongly consistent, converging (almost surely) to the true expectation as the number of samples $N$ tend to infinity. An alternate view of is is to consider it an (empirical) approximation111This should be interpreted as an approximation of the underlying probability distribution, and not of the density function itself. of $\gamma_{t}$,

[TABLE]

where $\delta_{X}$ denotes the Dirac measure at $X$. Furthermore, is provides an approximation of the normalization constant,

[TABLE]

Because the weights depend on the random samples, $x_{1:t}^{i}$, they are themselves random variables. One of the key properties is that it is unbiased, This can be easily seen by noting that $x_{1:t}^{i}$ are iid draws from $q_{t}$ and therefore

[TABLE]

since $Z_{t}=\int\widetilde{\gamma}_{t}(x_{1:t})\mathop{}\!\mathrm{d}x_{1:t}$ is nothing by the normalization constant for $\widetilde{\gamma}_{t}$. This property will be important for several of the more powerful is and smc-based methods considered in this monograph.

We summarize the importance sampling method in Algorithm 1. This algorithm is sometimes referred to as self-normalized is, because we are normalizing each individual weight using all samples.

The straightforward implementation of the is method we have described thus far is impractical for many of the example models and targets in Section 1.2. It is challenging to design good proposals for high-dimensional models. A good proposal is typically more heavy-tailed than the target; if it is not, the weights can have infinite variance. Another favorable property of a proposal is that it should cover the bulk of the target probability mass, putting high probability on regions of high probability under the target distribution. Even Markovian models, such as the ssm, can have a prohibitive computational complexity without careful design of the proposal. In the next section we will describe how we can alleviate these concerns using sis, a special case of is, with a kind of divide-and-conquer approach to tackle the high-dimensionality in $T$.

2.1.1 Sequential Importance Sampling

Sequential importance sampling (Handschin and Mayne, 1969; Robert and Casella, 2004) is a variant of is were we select a proposal distribution that has an autoregressive structure, and compute importance weights recursively. By choosing a proposal defined by

[TABLE]

we can decompose the proposal design problem into $T$ conditional distributions. This means we obtain samples $x_{1:t}^{i}$ by reusing $x_{1:t-1}^{i}$ from the previous iteration, and append a new sample, $x_{t}^{i}$, simulated from $q_{t}(x_{t}\,|\,x_{1:t-1}^{i})$. The unnormalized weights can be computed recursively by noting that

[TABLE]

We summarize the sis method in Algorithm 2, where ${q_{1}(x_{1}\,|\,x_{1:0})=q_{1}(x_{1})}$ and $\widetilde{w}_{0}=\widetilde{\gamma}_{0}=1$.

If we need to evaluate the normalization constant estimate $\widehat{Z}_{t}$, analogously to is we make use of Eq. 2.6. However, we may also obtain a (strongly) consistent estimate of the ratio of normalization constants $\nicefrac{{Z_{t}}}{{Z_{t-1}}}$

[TABLE]

While the estimate of the ratio is consistent, it is in general not unbiased. However, sis is a special case of is. This means that the sis estimate of the normalization constant for $\widetilde{\gamma}_{t}$, i.e. $\widehat{Z}_{t}$ in Eq. 2.6, is still unbiased and consistent.

In Example 2.1.1 we detail a first example proposal $q_{t}$ for the running example, and derive the corresponding weights $\widetilde{w}_{t}$. Furthermore, we include a code snippet that illustrates how to implement the sampler in Python.

Example 2.1.1** (Sequential importance sampling for Example 1.2.1).**

We revisit our running non-Markovian Gaussian example. The target distribution is

[TABLE]

with $p(x_{1})=\mathcal{N}(x_{1}\,|\,0,q)$ and

[TABLE]

A common approach is to set the proposal to be the prior (or transition) distribution $f$. A sample from the proposal $q_{t}(x_{t}\,|\,x_{t-1})=f(x_{t}\,|\,x_{t-1})$ is generated as follows

[TABLE]

We refer to this proposal simply as the prior proposal. The corresponding weight update is

[TABLE]

where $\widetilde{w}_{0}=1$. We provide LABEL:code:sis:running to illustrate how to implement sis with the prior proposal for this model in Python.

For improved numerical stability we update the log-weights $\log\widetilde{w}_{t}$ and subtract the logarithm of the sum of weights (the weight normalization) before exponentiating.

Sis can be implemented efficiently for a large class of problems, the computational complexity is usually linear in $N$ and $T$. Even so, the is methods suffer from severe drawbacks limiting their practical use for many high-dimensional problems.

2.1.2 Shortcomings of Importance Sampling

The main drawback of is is that the variance of the estimator scales unfavorably with the dimension of the problem; the variance generally increases exponentially in $T$ (Doucet and Johansen, 2009). Because sis is a special case of is it inherits this unfavorable property. To see that this is true we can illustrate using a toy example where the target factorizes completely over time and with identical distributions for each factor $\widetilde{\gamma}_{t}(x_{1:t})=\prod_{k=1}^{t}\widetilde{\gamma}(x_{k})$. We illustrate in Example 2.1.2.

Example 2.1.2** (Sequential importance sampling variance of $\widehat{Z}_{t}$).**

Assume that the target distribution of interest is given by $\widetilde{\gamma}_{t}(x_{1:t})=\prod_{k=1}^{t}\widetilde{\gamma}(x_{k})$, and we pick a proposal distribution that follows the same factorization $q_{t}(x_{1:t})=\prod_{k=1}^{t}q(x_{k})$. The normalization constant is then given by $Z_{t}=Z^{t}$ for $Z=\int\widetilde{\gamma}(x)\mathop{}\!\mathrm{d}x$. The normalization constant estimate is

[TABLE]

and its corresponding variance, normalized by $Z_{t}$, is

[TABLE]

where $\mathbb{E}_{q(x)}\left[\frac{\gamma(x)^{2}}{q(x)^{2}}\right]=\mathbb{E}_{q(x_{k})}\left[\frac{\gamma(x_{k})^{2}}{q(x_{k})^{2}}\right]$ for all $k$. This means that if the proposal is not exactly equal to the target the variance of the normalization constant estimate will in general increase exponentially with $t$.

Another way in which the unfavorable scaling in $t$ manifests itself in practice is through the normalized weights $w_{t}^{i}$. The maximum of the weights, $\max_{i}w_{t}^{i}$, will quickly approach one as $t$ increases, and consequently all other weights approach zero; a phenomena known as weight degeneracy. This means that, effectively, we approximate the target distribution using a single sample.

We illustrate weight degeneracy in Example 2.1.3 using the running example.

Example 2.1.3** (SIS weight degeneracy).**

We return to our running example, set the length $T=6$, number of particles $N=5$, and $(\phi,q,\beta,r)=(0.9,1.0,0.5,1.0)$. Fig. 2.1 shows the particles and the normalized weights $w_{t}^{i}$, where the area of the discs correspond to the size of the weights. We can see that as $t$ increases nearly all mass concentrates on the fourth particle $x_{t}^{4}$. This means that the normalized weight of the particle is almost one, $w_{t}^{4}\approx 1$. The remaining particles have normalized weights that are all close to zero, and thus have a negligible contribution to the approximation.

This concentration of mass for sis to a single particle happens very quickly. Even for very simple Markovian models the variance of e.g. our normalization constant estimator can increase exponentially fast as a function of $T$.

Sequential Monte Carlo methods tackle the weight degeneracy issue by choosing a proposal that leverages information contained in $\widehat{\gamma}_{t-1}$, the previous iteration’s target distribution approximation.

2.2 Sequential Monte Carlo

Sequential Monte Carlo methods (Gordon et al., 1993; Stewart and McCarty, 1992; Kitagawa, 1993) improve upon is by mitigating the weight degeneracy issue through a clever choice of proposal distribution. For certain sequence models the weight degeneracy issue can be resolved altogether, providing estimators to the final marginal distribution $\gamma_{T}(x_{T})$ that do not deteriorate for increasing $T$. For other sequence models, smc still tends to provide more accurate estimates in practice compared to is.

Just like in sis we need a sequence of proposal distributions $q_{t}(x_{t}\,|\,x_{1:t-1})$ for $t=1,\ldots,T$. This is a user choice that can significantly impact the accuracy of the smc approximation. For now we assume that the proposal is given and return to this choice in Chapter 3. Below, we detail the iterations (or steps) of a basic smc algorithm.

Step 1:

The first iteration of smc boils down to approximating the target distribution $\gamma_{1}$ using standard is. Simulating $N$ times independently from the first proposal

[TABLE]

and assigning corresponding weights

[TABLE]

lets us approximate $\gamma_{1}$ (cf. Eq. 2.5) by

[TABLE]

The key innovation of the smc algorithm is that it takes advantage of the information provided in $\widehat{\gamma}_{1}(x_{1})$, Eq. 2.12, when constructing a proposal for the next target distribution $\gamma_{2}$.

Step 2:

In the second iteration of smc we sample $x_{1:2}$ from the proposal $\widehat{\gamma}_{1}(x_{1})q_{2}(x_{2}\,|\,x_{1})$, rather than from $q_{1}(x_{1})q_{2}(x_{2}\,|\,x_{1})$ like sis. We sample $N$ times independently from

[TABLE]

and assign weights

[TABLE]

Simulating $x_{1:2}^{i}$, Eq. 2.13, can be broken down into parts: resampling $x_{1}^{i}\sim\widehat{\gamma}_{1}(x_{1})$, propagation $x_{2}^{i}\,|\,x_{1}^{i}\sim q_{2}(x_{2}\,|\,x_{1}^{i})$, and concatenation ${x_{1:2}^{i}=(x_{1}^{i},x_{2}^{i})}$. Note the overloaded notation for $x_{1}^{i}$. We replace the initial sample set $\{x_{1}^{i}\}_{i=1}^{N}$ from Step 1, with the resampled set $x_{1}^{i}\sim\widehat{\gamma}_{1}(x_{1})$, $i=1,\ldots,N$.

Resampling can refer to a variety of methods in statistics, for our purpose it is simple (weighted) random sampling with replacement from $x_{1}^{1:N}=\{x_{1}^{i}\}_{i=1}^{N}$ with weights $w_{1}^{1:N}=\{w_{1}^{i}\}_{i=1}^{N}$. Resampling $N$ times independently means that the number of times each particle is selected is multinomially distributed. This resampling algorithm is known as multinomial resampling, see LABEL:code:multinomial_resampling. In Sections 2.2.1 and 2.2.2 we revisit resampling and present alternative resampling algorithms, increasing efficiency by correlation and adaptation.

Propagation generates new samples independently from the proposal, i.e. $x_{2}^{i}\sim q_{2}(x_{2}\,|\,x_{1}^{i})$ for each $i=1,\ldots,N$. By concatenating $x_{1:2}^{i}=(x_{1}^{i},x_{2}^{i})$ we obtain a complete sample from the proposal $\widehat{\gamma}_{1}(x_{1})q_{2}(x_{2}\,|\,x_{1})$.

Finally, new importance weights are computed according to (2.14). This expression can be understood as a standard importance sampling weight—the target divided by the proposal. Note, however, that we use $\widetilde{\gamma}_{1}$ in the denominator: when computing the weights we “pretend” that the proposal is given by $\gamma_{1}(x_{1})q_{2}(x_{2}\,|\,x_{1})$ rather than by its approximation $\widehat{\gamma}_{1}(x_{1})q_{2}(x_{2}\,|\,x_{1})$. This is of course just an interpretation and the motivation for using this particular weight expression is given by the convergence analysis for smc; see Section 2.3.

The resulting approximation of $\gamma_{2}$ is

[TABLE]

smc can in essence be described as a synthesis of sis and resampling, which explains its alternate names sequential importance resampling (sir) or sequential importance sampling and resampling (sisr).

Step $t$:

The remaining iterations follow the recipe outlined in step $2$. First, the proposal is the product of the previous empirical distribution approximation and a conditional distribution

[TABLE]

Samples for $i=1,\ldots,N$ are generated as follows

[TABLE]

Finally, we assign the weights

[TABLE]

and approximate $\gamma_{t}$ by

[TABLE]

The normalization constant $Z_{t}$ can be estimated by

[TABLE]

We summarize the full sequential Monte Carlo sampler in Algorithm 3, where $\widehat{\gamma}_{0}=\widetilde{\gamma}_{0}=1$, and ${q_{1}(x_{1}\,|\,x_{1:0})=q_{1}(x_{1})}$.

The smc method typically achieves drastic improvements compared to sis. In Example 2.2.1 we return to our running example, using the same settings as in Example 2.1.3, to study the sample diversity and quality of the basic smc method compared to sis.

Example 2.2.1** (SMC sample diversity).**

We illustrate the weights and resampling dependencies in Fig. 2.2. The grey arrows represent what sample from iteration $t-1$ that generated the current sample at iteration $t$, referred to as its ancestor.

We can see that the weights tend to be more evenly distributed for smc. The algorithm dynamically chooses to focus computational effort on more promising samples through the resampling step. Particles with low weights tend to not be resampled, and particles with high weights are resampled more frequently.

Not only do we get more diversity in our sample set, smc also tends to find areas of higher probability. We illustrate this phenomenon in Table 2.1. We fix the number of particles $N=10$, then study the average log-probability of our target distribution $\log\widetilde{\gamma}_{T}$, under the sampling distributions $\widehat{\gamma}_{T}$, normalized by the number of iterations $T$.

While smc methods do not suffer from weight degeneracy as $t$ increases, the fundamental difficulty of simulating from a high-dimensional target still remains. For smc this instead manifests itself in another form of weight degeneracy and something known as path degeneracy. We discuss weight degeneracy below and illustrate path degeneracy in Example 2.2.2.

Weight degeneracy in SMC:

While introducing resampling can alleviate the weight degeneracy stemming from sis’s unfavorable scaling in $t$, smc can still suffer from weight degeneracy due to a high-dimensional $x_{t}$. If at each iteration of Algorithm 3 the maximum of the normalized importance weight $\max_{i}w_{i}$ is close to one, and the remaining close to zero, we say that the smc sampler suffers from weight degeneracy. This occurs when the variance of the unnormalized importance weights $\widetilde{w}_{t}^{i}$ is high. Snyder et al. (2008); Bickel et al. (2008); Bengtsson et al. (2008); Snyder et al. (2015) showed that to avoid weight degeneracy the number of particles $N$ must scale at least exponentially with the variance of the log-weights, which in simple examples is proportional to the dimension $n_{x}$ of the latent variable $x_{t}$. These results rely on a Gaussian target, are only valid as both $N$ and $n_{x}$ tend to infinity, and focus on the variability of the log-weights. Naesseth et al. (2019, Proposition 3) showed that asymptotic variance of the resulting estimator of simple test functions can also scale exponentially in the dimension $n_{x}$, relying on asymptotics only in $N$ and not $n_{x}$. We return to this example as we discuss the central limit theorem (clt) in Section 2.3.

The weight degeneracy in smc caused by high-dimensional latent variables is still an open research problem, see e.g. Djurić and Bugallo (2013); Briggs et al. (2013); Beskos et al. (2014); Rebeschini and Van Handel (2015); Naesseth et al. (2019) and the references therein. In this tutorial we focus on how choosing good proposals and target distributions can ameliorate, but not completely solve, this important issue.

Example 2.2.2** (SMC path degeneracy).**

In Fig. 2.3 we reiterate the result from our previous example, Example 2.2.1. However, this time we only include the arrows corresponding to samples that have been consistently resampled and form our final approximation $\widehat{\gamma}_{T}$.

We can see that our approximation for early iterations collapses back to a single ancestor, e.g. we have $N=5$ identical copies of $x_{1}^{1}$ in Fig. 2.3. This phenomena is known as path degeneracy and occurs because of the resampling mechanism. In Jacob et al. (2015) the authors study the impact this has on the memory requirements. They show that for state space models, under suitable conditions on the observation pdf $g$, the expected distance from the current iteration to a coalescence of the paths is bounded from above by $\mathcal{O}(N\log N)$.

In Fig. 2.4 we study the impact that increasing dependence on earlier iterations has on our smc estimate of the log-normalization constant. We let $N=20$, $T=100$ be fixed and vary the value of $\beta\in(0,1)$, where increasing values of $\beta$ correspond to more long-range dependencies in our non-Markovian lvm. We can see that for modest values of $\beta$ the smc method achieves accurate estimates, whereas for higher values the drop-off in efficiency is significant.

This manifests itself as an increase in the Monte Carlo (mc) variance in the estimate (width of bars), as well as a negative bias. As we will discuss in Section 2.3, this negative bias in the estimate of $\log Z_{T}$ is typical for smc and is.

In Sections 2.2.1 and 2.2.2 we will discuss two standard practical approaches that help alleviate (but not solve) the issue of path degeneracy and improve overall estimation accuracy. The first is low-variance resampling — we lower the variance in the resampling step by correlating random variables. The second is adaptive resampling, which essentially means that we do not resample at each iteration but only on an as-needed basis. In Chapter 3 we go a step further and show how to choose, or learn, proposal and target distributions that lead to even further improvements.

2.2.1 Low-variance Resampling

The resampling step introduces extra variance in the short-term, using $\widehat{\gamma}_{t-1}$ to estimate $\gamma_{t-1}$ is better than using the resampled particle set. However, discarding improbable particles, and putting more emphasis on promising ones is crucial to the long-term performance of smc. To keep long-term performance but minimize the added variance, we can employ standard variance reduction techniques based on correlating samples drawn from $\widehat{\gamma}_{t-1}$. First, we explain a common technique for implementing multinomial resampling based on inverse transform sampling. Then, we explain two low-variance resampling alternatives, stratified and systematic resampling.

To sample from $\widehat{\gamma}_{t-1}$ we use inverse transform sampling based on the weights $w_{t-1}^{i}$. We have that

[TABLE]

where $a$ is an integer random variable on $\{1,\ldots,N\}$, such that the following is true

[TABLE]

for $u\sim\mathcal{U}(0,1)$. Repeating the above process independently $N$ times gives multinomial resampling. That is, we draw $u^{i}\sim\mathcal{U}(0,1)$ and find the corresponding $a^{i}$ for each $i=1,\ldots,N$. In the context of smc the index variables $\{a^{i}\}_{i=1}^{N}$ determine the ancestry of the particles and they are therefore often referred to as ancestor indices.

Stratified Resampling

One way to improve resampling is by stratification on the uniform random numbers $u^{i}$ (Kitagawa, 1996; Fearnhead, 1998). This means we divide the unit interval into $N$ strata,

${(0,\nicefrac{{1}}{{N}}),(\nicefrac{{1}}{{N}},\nicefrac{{2}}{{N}}),\ldots,(1-\nicefrac{{1}}{{N}},1)}$. Then, we generate $u^{i}\sim\mathcal{U}(\nicefrac{{i-1}}{{N}},\nicefrac{{i}}{{N}})$ for each strata $i=1,\ldots,N$. Finally, the corresponding ancestor indices $a^{i}$ are given by studying Eq. 2.19.

LABEL:code:stratified_resampling shows how this can be implemented in Python. The main change compared to multinomial resampling, LABEL:code:multinomial_resampling, is the way the $u^{i}$’s are generated.

Systematic Resampling

We can take correlation to the extreme by only generating a single uniform number $u\sim\mathcal{U}(0,1)$ to set all the $u^{i}$’s. Systematic resampling (Carpenter et al., 1999; Whitley, 1994) means that we let

[TABLE]

where the random number $u\sim\mathcal{U}(0,1)$ is identical for all $i=1,\ldots,N$. Then, we find corresponding indices $a^{i}$ by again studying Eq. 2.19. Note that just like in stratified resampling, systematic resampling also generates one $u^{i}$ in each strata $(\nicefrac{{i-1}}{{N}},\nicefrac{{i}}{{N}})$. However, in this case the $u^{i}$’s are based on a single random value $u$. This means that systematic resampling will be more computationally efficient than stratifies and multinomial resampling.

The code change to implement systematic resampling is simple. We only need to change a single line in LABEL:code:stratified_resampling. We replace line $3$ by: u = (np.arange(N) + npr.rand())/N.

Both systematic and stratified resampling are heavily used in practice. In many cases systematic resampling achieves slightly better results (Hol et al., 2006). However, systematic resampling can sometimes lead to non-convergence (Gerber et al., 2017), depending on the ordering of the samples.

For a more thorough discussion on these, and other, resampling methods see e.g. Douc and Cappé (2005); Hol et al. (2006).

2.2.2 Effective Sample Size and Adaptive Resampling

The resampling step introduces extra variance by eliminating particles with low weights, and replicating particles with high weights. If the variance of the normalized weights is low, this step might be unnecessary. By tracking the variability of the normalized weights, and trigger a resampling step only when the variability crosses a pre-specified threshold, we can alleviate this issue. This is known as adaptive resampling in the smc literature. Often we study the effective sample size (ess) (Liu, 2004) to gauge the variability of the weights, which for iteration $t$ is

[TABLE]

The ess is a positive variable taking values in the continuous range between $1$ and $N$. For is the ess is an approximation of the number of exact samples from $\gamma_{t}$ that we would need to achieve a comparable estimation accuracy (Doucet and Johansen, 2009). A common approach is to resample only at iterations when the ess falls below $\nicefrac{{N}}{{2}}$.

Adaptive resampling can be implemented by slight alterations to the proposal and weight updates in Eqs. 2.15 and 2.16. If $\textrm{ESS}_{t-1}$ is above the prespecified threshold we simply omit the resampling step in Eq. 2.15 and obtain

[TABLE]

The fact that resampling is omitted at some iterations needs to be accounted for when computing the weights. For iterations where we do not resample the weight update is,

[TABLE]

This can be thought of as adding an extra importance correction, where the previous weights $\{w_{t-1}^{i}\}_{i=1}^{N}$ define a “target distribution” for the ancestor indices and the factor $\nicefrac{{1}}{{N}}$ is the “proposal” corresponding to not resampling. When the ess falls below our pre-set threshold, we use the standard update Eqs. 2.15 and 2.16 with resampling instead.

Adaptive resampling is usually combined with the low-variance resampling techniques explained above for further variance reduction.

Remark 2.2.1**.**

The constant factors in the weight expression will cancel when normalizing the weights. However, these constants still need to be included for the expression for the normalizing constant estimate Eq. 2.18 to be valid. That is, the extra importance correction added at those iterations when we do not resample should be the ratio of the normalized weights from the previous iteration, divided by the constant $\nicefrac{{1}}{{N}}$. An alternative approach appearing in the literature is to neglect the constants when computing the weights, but then modify the expression for the normalizing constant estimate Eq. 2.18 instead.

2.3 Analysis and Convergence

Since its conception in the 1990s, significant effort has been spent on studying the theoretical properties of smc methods. We will review and discuss a few select results in this section. For an early review of the area, see e.g. Del Moral (2004). The theorems we discuss below all hold for a number of conditions on the proposal, probabilistic model, and test functions. For brevity we have omitted these exact conditions and refer to the cited proofs for details.

Unbiasedness

One of the key properties of smc approximations is that they provide unbiased approximations of integrals of functions $h_{t}$ with respect to the unnormalized target distribution $\widetilde{\gamma}_{t}$. We formalize this in Theorem 2.3.1.

Theorem 2.3.1** (Unbiasedness).**

[TABLE]

Proof.

See Del Moral (2004, Theorem 7.4.2). For the special case $h_{t}\equiv 1$ see also Section 4.1 and Appendix 4.A. ∎

A particularly important special case is when $h_{t}\equiv 1$ and we approximate the normalization constant of $\widetilde{\gamma}_{t}$. We have that $\mathbb{E}[\widehat{Z}_{t}]=Z_{t}$, where the expectation is taken with respect to all the random variables generated by the smc algorithm. If we instead consider the more numerically stable $\log\widehat{Z}_{t}$, we have by Jensen’s inequality that $\mathbb{E}\left[\log\widehat{Z}_{t}\right]\leq\log Z_{t}$. This means that the estimator of the log-normalization constant is negatively biased. This is illustrated by the violin plot discussed in Example 2.2.2.

We will delve deeper into the applications of the unbiasedness property and its consequences in Chapter 4.

Laws of Large Numbers

While integration with respect to unnormalized distributions can be estimated unbiasedly, this is unfortunately not true when estimating expectations with respect to the normalized target distribution $\gamma_{t}$. However, smc methods are still strongly consistent, leading to exact solutions when the number of particles $N$ tend to infinity. We formalize this law of large numbers in Theorem 2.3.2.

Theorem 2.3.2** (Law of Large Numbers).**

[TABLE]

Proof.

See Del Moral (2004, Theorem 7.4.3). ∎

Central Limit Theorem

While the law of large numbers from the previous section shows that smc approximations are exact in the limit of infinite computation, it tells us nothing about the quality of our estimate. The clt in Theorem 2.3.3 tells us about the limiting distribution of our smc estimate and its asymptotic variance. This gives us a first approach to get an understanding for the precision of our smc approximation to $\gamma_{t}(h_{t})$.

Theorem 2.3.3** (Central Limit Theorem).**

[TABLE]

where $\textrm{V}_{t}(\cdot)$ is defined recursively for a measurable function $h$,

[TABLE]

initialized with $\widetilde{\textrm{V}}_{1}(h)=\textrm{Var}_{q_{1}}(h)$ and where

[TABLE]

Proof.

See Chopin (2004). ∎

Example 2.3.1** (State Space Model).**

We study the following example of a state space model

[TABLE]

where $n_{x}$ denotes the dimension of $x_{l}$ and $y_{l}$. For simplicity assume that $y_{l,k}=y_{l,m}$, $\forall k,m$ and that $\mathbb{E}_{p(\mathbf{x}\,|\,\mathbf{y})}\left[\mathbf{x}\right]=0$. Let the function of interest be $h_{t}(x_{1:t})=\sum_{k=1}^{n_{x}}x_{t,k}$, i.e. the sum of the components of the state at the most recent iteration. Studying the asymptotic variance in Theorem 2.3.3 for the so called fully adapted smc, explained in Section 3.2.2, we get

[TABLE]

where the constants $A_{l}$, $B_{l}$, and $C_{l}$ are defined by

[TABLE]

For the proof see Naesseth et al. (2019, Proposition 3).

This result tells us that the asymptotic variance is exponential in $n_{x}$, the dimension of $x_{t}$. Furthermore, it tells us that for the asymptotic variance to be stable in $t$, i.e. to not be increasing with increasing $t$, we need that our model forgets its initial conditions. Concretely, we need the inner expectations in $A_{l}$ and $C_{l}$ to tend to [math] (the posterior mean) quickly enough to be able to ensure that we can bound the sum by a (finite) constant that does not depend on $t$. This holds more generally (under various assumptions on the model and proposal) for approximations of the filtering distribution, $\gamma_{t}(x_{t})=p(x_{t}\,|\,y_{1:t})$, of ssms. We formalize the result in Proposition 2.3.1.

Proposition 2.3.1** (Asymptotic Variance for Filtering Approximation in SSM).**

Under appropriate forgetting conditions on the model

[TABLE]

for some finite constant $\mathrm{c}$.

Proof.

See Whiteley (2013). ∎

Proposition 2.3.1 shows that the asymptotic variance can be stable as a function of $t$, this is in general only true for test functions that depend only on the most recent latent variable $x_{t}$ and under various assumptions on the model and proposal (Del Moral and Guionnet, 2001; Del Moral, 2004; Whiteley, 2013). When considering the variance of $\nicefrac{{\widehat{Z}_{t}}}{{Z_{t}}}$ of ssms, it is instead possible to show that the variance increases only linearly with $t$ (Whiteley, 2013; Del Moral, 2004). This implies that to obtain a good estimate of the normalization constant, we must choose the number of particles $N$ to scale at least linearly in $t$.

Sample Bounds

Another way of looking at the smc method, disregarding test functions, is as a direct approximation to the target distribution itself. The weights and samples generated by the smc algorithm provide an approximation $\widehat{\gamma}_{t}\left(x_{1:t}\right)$, Eq. 2.17, to the true target distribution $\gamma_{t}(x_{1:t})$.

With this point of view we can study bounds on the difference between the distribution of a sample drawn from the smc approximation compared to that of a sample drawn from the target distribution. Specifically, assume that we generate a sample $x_{1:t}^{\prime}$ by first running an smc sampler to generate an approximation $\widehat{\gamma}_{t}$ of $\gamma_{t}$, and then simulate $x_{1:t}^{\prime}\sim\widehat{\gamma}_{t}$. Then, the marginal distribution of $x_{1:t}^{\prime}$ is $\mathbb{E}[\widehat{\gamma}_{t}]$, where the expectation is taken with respect to all random variables generated by the smc algorithm. It is worth noting that this distribution, $\mathbb{E}[\widehat{\gamma}_{t}]$, may be continuous despite the fact that $\widehat{\gamma}_{t}$ is a point-mass distribution by construction. In Theorem 2.3.4 we restate a generic sample bound on the kl divergence from the expected smc approximation $\mathbb{E}\left[\widehat{\gamma}_{t}\right]$ to the target distribution $\gamma_{t}$. The kl divergence from a distribution $q$ to another distribution $p$ is defined by $\textrm{KL}\,\left(p(x)\|q(x)\right):=\int p(x)\left(\log p(x)-\log q(x)\right)\mathop{}\!\mathrm{d}x$.

Theorem 2.3.4** (Sample Bound).**

[TABLE]

for a finite constant $\mathcal{C}$.

Proof.

See Del Moral (2004, Theorem 8.3.2). ∎

From this result we can conclude that in fact the smc approximation tends to the true target in distribution as the number of particles increase.

Recently there has been an increased interest for using smc as an approximation to the target distribution, rather than just as a method for estimating expectations with respect to test functions. This point of view has found applications in e.g. probabilistic programming and variational inference (Wood et al., 2014; Naesseth et al., 2018; Huggins and Roy, 2015).

Chapter 3 Learning Proposals and Twisting Targets

The two main design choices of sequential Monte Carlo methods are the proposal distributions $q_{t}$, $t\in\{1,\dots,T\}$ and the intermediate (unnormalized) target distributions $\widetilde{\gamma}_{t}$, $t\in\{1,\dots,T-1\}$. Carefully choosing these can drastically improve the efficiency of the algorithm and the accuracy of our estimation. Adapting both the proposal and the target distribution is especially important if the latent variables are high-dimensional.

In the first section below we discuss how to choose, or learn, the proposal distribution for a fixed target distribution. Then, in the final section we discuss how we can design a good sequence of intermediate target distributions for a given probabilistic model.

3.1 Designing the Proposal Distribution

The choice of the proposal distribution is perhaps the most important design choice for an efficient sequential Monte Carlo algorithm. A common choice is to propose samples from the model prior; this is simply known as the prior proposal. When using the prior as proposal the algorithm is commonly known as the bootstrap particle filter or bootstrap smc (Gordon et al., 1993). However, using the prior can lead to poor approximations for a small number of particles, especially if the latent space is high-dimensional.

We will in this section derive the locally (one-step) optimal proposal distribution, or optimal proposal for short. Because it is typically intractable, we further discuss various ways to either emulate it or approximate it directly. Finally, we discuss alternatives to learn efficient proposal distributions using an end-to-end variational perspective.

3.1.1 Locally Optimal Proposal Distribution

The locally optimal proposal distribution (Doucet et al., 2000) is the distribution we obtain if we assume that we have a perfect approximation at iteration $t-1$. Then, choose the proposal $q_{t}$ that minimizes the kl divergence from the joint distribution $\gamma_{t-1}(x_{1:t-1})q_{t}(x_{t}\,|\,x_{1:t-1})$ to $\gamma_{t}(x_{1:t})$. We formalize this result in Proposition 3.1.1. Equivalently, we can view this as the proposal minimizing the variance of the incremental weights, i.e. $\nicefrac{{\widetilde{w}_{t}^{i}}}{{\widetilde{w}_{t-1}^{i}}}$, with respect to the newly generated samples $x_{t}^{i}$.

Proposition 3.1.1** (Locally optimal proposal distribution).**

The optimal proposal $q_{t}^{\star}(x_{t}\,|\,x_{1:t-1})$ minimizing $\textrm{KL}\left(\gamma_{t-1}(x_{1:t-1})q_{t}(x_{t}\,|\,x_{1:t-1})\|\gamma_{t}(x_{1:t})\right)$ is given by

[TABLE]

where $\widetilde{\gamma}_{t}(x_{1:t-1})=\int\widetilde{\gamma}_{t}(x_{1:t})\mathop{}\!\mathrm{d}x_{t}$.

Proof.

If we let “const” denote terms constant with respect to the proposal distribution $q_{t}(x_{t}\,|\,x_{1:t-1})$, we get

[TABLE]

where the inner (conditional) Kullback-Leibler divergence is zero if and only if $q_{t}(x_{t}\,|\,x_{1:t-1})\equiv\gamma_{t}(x_{t}\,|\,x_{1:t-1})$. ∎

In Example 3.1.1 we show that the optimal proposal is analytically tractable for our running non-Markovian Gaussian example.

Example 3.1.1** (Optimal proposal for Example 1.2.1).**

If we let $\widetilde{\gamma}_{t}(x_{1:t})=\widetilde{\gamma}_{t-1}(x_{1:t-1})f(x_{t}\,|\,x_{t-1})g(y_{t}\,|\,x_{1:t})$, then the (locally) optimal proposal for our running example is analytically tractable. It is

[TABLE]

In most practical cases the optimal proposal distribution is not a feasible alternative. The resulting importance weights are intractable, or simulating random variables from the optimal proposal is too computationally costly. Below we discuss various common approaches for approximating it.

3.1.2 Approximations to the Optimal Proposal Distribution

There have been a number of suggestions over the years on how to approximate the optimal distribution. We will review three analytic approximations based on a Gaussian assumption, as well as briefly describe an exact approximation.

Laplace Approximation

We obtain the Laplace approximation to the optimal proposal by a second-order Taylor approximation of the log-pdf around a point $\bar{x}$ (Doucet et al., 2000). If we let $l_{t}(x_{t}):=\log\gamma_{t}(x_{t}\,|\,x_{1:t-1})$, supressing the dependence on $x_{1:t-1}$, then

[TABLE]

where $\nabla l_{t}$ and $\nabla^{2}l_{t}$ are the gradient and the Hessian of the log-pdf with respect to $x_{t}$, respectively. A natural approximation to the optimal proposal is then

[TABLE]

The mode of the distribution can be a good choice for the linearization point $\bar{x}$ if the distribution is unimodal. With this choice the mean simplifies to just $\bar{x}$. However, the mode is usually unknown and will depend on the value of $x_{1:t-1}$. This means that we are required to run a separate optimization for each particle $x_{1:t-1}^{i}$ and iteration $t$ to find the mode and Hessian. This can outweigh the benefits of the improved proposal distribution.

Extended and Unscented Kalman Filter Approximations

The extended Kalman filter (ekf) and the unscented Kalman filter (ukf) (Anderson and Moore, 1979; Julier and Uhlmann, 1997) are by now standard solutions to non-linear and non-Gaussian filtering problems. We can leverage these ideas to derive another class of Gaussian approximations to the optimal proposal distribution (Doucet et al., 2000; Van Der Merwe et al., 2001).

The two methods depend on a structural equation representation of our probabilistic model,

[TABLE]

where $v_{t}$ and $e_{t}$ are random variables. The representation implies a joint distribution on $x_{t}$ and $y_{t}$ conditional on $x_{1:t-1}$, i.e. $p(x_{t},y_{t}\,|\,x_{1:t-1})$. The locally optimal proposal distribution corresponding to this representation is given by $q_{t}^{\star}(x_{t}\,|\,x_{1:t-1})=p(x_{t}\,|\,x_{1:t-1},y_{t})$.

The ekf and ukf approximations rely on a Gaussian approximation to the joint conditional distribution of $x_{t}$ and $y_{t}$,

[TABLE]

where we have suppressed the dependence on $x_{1:t-1}$ for clarity. We block the mean $\widehat{\mu}$ and covariance $\widehat{\Sigma}$ as follows

[TABLE]

Under the assumptions of the approximation in Eq. 3.6, the distribution of $x_{t}\,|\,x_{1:t-1},y_{t}$ is tractable:

[TABLE]

with

[TABLE]

The key difference between the ekf- and ukf-based approximations is how to compute the estimates $\widehat{\mu}$ and $\widehat{\Sigma}$. We leave the details of these procedures to Appendix 3.A, and focus here on the intuition behind them.

The ekf uses a first-order Taylor approximation to the non-linear functions $a,c$ to derive an approximation to the full posterior distribution based on exact (analytical) updates of the approximate model. Following this line of thinking we can similarly linearize $a,c$ locally. Then under a Gaussian assumption on the noise $v_{t},e_{t}$, we compute the distribution $\widehat{p}(x_{t},y_{t}\,|\,x_{1:t-1})$ exactly for the linearized model. The ekf-based approximations $\widehat{\mu}$ and $\widehat{\Sigma}$ can be found in Eq. 3.53 in the appendix.

The ukf on the other hand uses so-called sigma points to reach the Gaussian approximation $\widehat{p}(x_{t},y_{t}\,|\,x_{1:t-1})$. The key idea is to choose a set of points, pass them through the nonlinear functions $a$ and $c$, and then estimate the mean and variance of the transformed point-set. Unlike the ekf-based approximation, the ukf-based approximation does not require that the noises $v_{t}$ and $e_{t}$ are Gaussian distributed. However, we do require that the mean and variance for these random variables are available. The ukf-based approximations $\widehat{\mu}$ and $\widehat{\Sigma}$ can be found in Eq. 3.59 in the appendix.

Analytic Gaussian Approximations

We refer to the Laplace-, ekf- and ukf-based approximations as analytic Gaussian approximations to the locally optimal proposal distribution. We summarize these three approaches in Algorithm 4. Which one works best depends heavily on the model being studied.

In Fig. 3.1 we study a simple example when $x$ and $y$ are one-dimensional random variables.

Example 3.1.2** (Gaussian Proposal Approximations).**

We illustrate the analytic Gaussian approximations from Algorithm 4 based on a simple scalar example

[TABLE]

The prior on the latent variable $x$ is a standard normal, but the measurement model $c$ is a polynomial in $x$ with additive Gaussian noise. In Fig. 3.1 we show the true (bimodal) posterior $p(x\,|\,y)$, and the corresponding approximations based on the Laplace, ekf, and ukf methods. Neither of these methods can capture the bimodality of the true normalized distribution, since they are all based on the basic simplifying assumption that the posterior is Gaussian. However, as previously discussed a good proposal will cover the bulk of probability mass of the target. This means it is likely that the ekf and ukf proposals will outperform the Laplace proposal in this case.

Exact Approximations

In exact approximations we use another level of mc. Instead of drawing exact samples from $q_{t}^{\star}$, we approximate draws from it using a nested mc algorithm (Naesseth et al., 2015, 2019). This nested algorithm is chosen such that it does not alter the asymptotic exactness of the outer smc method. We have discussed how smc and is can be used as distribution approximators. Nested mc methods takes this one step further, using a separate smc (or is) approximation $\widehat{q}_{t}^{\star}$ for each sample from $q_{t}^{\star}$ we would like. This can lead to substantial benefits for e.g. high-dimensional latent variables $x_{t}$. One of the key differences between the analytic Gaussian approximations and the exact approximations is the fixed form distribution assumption on the analytic approximation. The exact approximation is not limited to a standard parametric distribution, and we can obtain arbitrarily accurate approximations with enough compute. The trade-off is the additional effort needed to generate each sample $x_{t}^{i}$, which requires us to run independent smc algorithms for each particle $i$ at each iteration $t$. While this might seem wasteful, it can in fact improve accuracy compared to smc methods with standard proposals (Naesseth et al., 2015), even for equal compute.

We are interested in approximating the locally optimal proposal distribution, i.e. $q_{t}^{\star}(x_{t}\,|\,x_{1:t-1})\propto\widetilde{\gamma}_{t}\left(x_{1:t}\right)$, separately for each particle $i$ using an mc method. A first approach is to use a nested is sampler with a proposal $r_{t}(x_{t}\,|\,x_{1:t-1})$ targeting $q_{t}^{\star}(x_{t}\,|\,x_{1:t-1})$. We construct an approximation of the optimal proposal independently for each particle $x_{1:t-1}^{i}$ as follows

[TABLE]

where the weights $\widetilde{\nu}_{t}$ are given by

[TABLE]

We replace $q_{t}(x_{t}\,|\,x_{1:t-1})$ in Eq. 2.15 with samples from the distribution defined in Eq. 3.15. The corresponding weight update, under the assumption that we have resampled $x_{1:t-1}^{i}$, is given by

[TABLE]

We summarize the nested is approach to approximate the optimal proposal in Algorithm 5. The algorithm generates a single sample, and must be repeated for each $i=1,\ldots,N$.

By what we know from the theory of smc, see Section 2.3, we would expect that $\widehat{q}_{t}^{\star}\stackrel{{\scriptstyle\textrm{d}}}{{\longrightarrow}}q_{t}^{\star}$ if we let $M\to\infty$. This is true, and as shown by Naesseth et al. (2019) the asymptotic variance in the limit of $M\to\infty$ is equal to that of using the locally optimal proposal distribution. The main implication is that we can obtain arbitrarily accurate optimal proposal approximations at the cost of increasing the number of samples $M$ for the nested mc method. We can also combine the analytic Gaussian approximations above with exact approximations by choosing $r_{t}$ as either the Laplace, ekf, ukf, or any other analytic approximation.

The computational complexity of the nested is method is increased by a factor of $M$ compared to standard smc. A natural question is therefore whether or not we are better of running standard smc with $r_{t}$ as a proposal and using $NM$ particles. However, the memory requirement is lower for nested smc and there are opportunities for speeding up through parallelization. There are even variants of nested mc that can outperform on the same computational budget (Naesseth et al., 2015), i.e. when compared to standard smc with $NM$ particles without taking parallelization into account. We return to this class of methods in Section 4.3.2.

3.1.3 Learning a Proposal Distribution

Rather than relying on approximating the locally optimal proposal distribution, we can instead directly learn a good proposal distribution that can take into account global information of the target distribution. By parameterizing a suitable class of distributions and choosing a cost function to optimize, we can use standard optimization tools (such as stochastic gradient descent) to adapt our proposal to the problem we are trying to solve. When combined with the smc approximation, we refer to these types of approaches as adaptive or variational smc methods.

We will in this section focus on proposal distributions, $q_{t}$, parameterized by $\lambda$, we denote this as $q_{t}(x_{t}\,|\,x_{1:t-1}\,;\,\lambda)$. A common choice is to use a conditional Gaussian distribution

[TABLE]

where $\mu_{\lambda}(\cdot),\sigma_{\lambda}^{2}(\cdot)$ are neural networks parameterized by $\lambda$. The proposal is also a function of the data, either explicitly as a part of the input to the functions or implicitly through the cost function for the parameters. These types of conditional distributions have recently been used with success in everything from image and speech generation (Gulrajani et al., 2017; Maddison et al., 2017) to causal inference (Louizos et al., 2017; Krishnan et al., 2017).

Below we discuss two classes of methods to learn the parameters $\lambda$, adaptive and variational methods, respectively.

Adaptive Sequential Monte Carlo

Adaptive methods are characterized by choosing a cost function that tries to fit directly the proposal distribution $q_{t}(x_{1:T})$ to the target distribution $\gamma_{t}(x_{1:T})$. This can either be done locally at each iteration for $q_{t}(x_{t}\,|\,x_{1:t-1}\,;\,\lambda)$ as in e.g. Cornebise et al. (2008), or globally for $q_{T}(x_{1:T}\,;\,\lambda)$ as in e.g. Gu et al. (2015); Paige and Wood (2016). Because we are mainly interested in the approximation of the final target distribution $\gamma_{T}$, we will here focus our attention to methods that take a global approach to learning proposals. Below we focus on recent methods developed in the machine learning literature. However, there is also an important body of work on adaptive smc samplers (Del Moral et al., 2006a) stemming from the statistics community. For a comprehensive description of those methods we refer to Jasra et al. (2011); Fearnhead and Taylor (2013); Schäfer and Chopin (2013); Buchholz et al. (2018).

Since we generally need the proposal distribution to cover areas of high probability under the target distribution, we need the proposal to be more diffuse than the target. This will ensure that the weights we assign will have finite variance and that our approximation gets more accurate. This means that using the cost function $\textrm{KL}(q_{T}(x_{1:T}\,;\,\lambda)\|\gamma_{T}(x_{1:T}))$, like in standard variational inference (e.g. Blei et al., 2017), would not result in a good proposal distribution. This cost function will encourage a proposal that is more concentrated than the target, and not less. Instead we focus our attention on the inclusive kl divergence $\textrm{KL}(\gamma_{T}(x_{1:T})\|q_{T}(x_{1:T}\,;\,\lambda))$ like in e.g. Gu et al. (2015) or Paige and Wood (2016). The inclusive kl divergence encourages a $q_{T}$ that covers the high-probability regions of $\gamma_{T}$.

We focus our exposition on adaptive methods to the case when $\gamma_{T}(x_{1:T})$ is the posterior distribution $p(x_{1:T}\,|\,y_{1:T})$. The cost function to minimize is the inclusive kl divergence

[TABLE]

where “const” includes all terms constant with respect to the proposal distribution. If we could compute the gradients of Eq. 3.17 with respect to $\lambda$ we could employ a standard gradient descent method to optimize it. Unfortunately the gradients, given by

[TABLE]

requires us to compute expectations with respect to the posterior distribution. Computing these types of expectations is the problem we try to solve with smc in the first place! Below we will detail two approaches to solve this problem. First, we consider a stochastic gradient method that uses smc to estimate the gradients. Then, we consider an alternative solution that uses inference amortization combined with stochastic gradient methods.

In the first approach, we use our current best guess for the parameters $\lambda$ and then the smc procedure itself to approximate the gradient in Eq. 3.18. Given this estimate of the gradient we can do an update step based on standard stochastic gradient descent methods. Suppose we have the iterate $\lambda^{n-1}$, we estimate the gradient in Eq. 3.18 using Algorithm 3 with proposals $q_{t}(x_{t}\,|\,x_{1:t-1}\,;\,\lambda^{n-1})$, and get

[TABLE]

With stochastic gradient descent we update our iterate by

[TABLE]

where $\alpha_{n}$ are a set of positive step-sizes, typically chosen such that $\sum_{n}\alpha_{n}=\infty$ and $\sum_{n}\alpha_{n}^{2}<\infty$. Note that unlike standard stochastic optimization methods, the gradient estimate is biased and convergence to a local minima of the cost function is not guaranteed. However, this approach has still shown to deliver useful proposal adaptation in practice (Gu et al., 2015).

The other approach relies on amortizing inference (Gershman and Goodman, 2014; Paige and Wood, 2016). One view of inference amortization is as a procedure that does not just learn a single optimal $\lambda^{\star}$ for the observed dataset $y_{1:T}$, but rather learns a mapping $\lambda^{\star}(\cdot)$ from the data space to the parameter space. This mapping is optimized to make sure $q_{T}(x_{1:T}\,;\,\lambda^{\star}(y_{1:T}))$ is a good approximation to $p(x_{1:T}\,|\,y_{1:T})$ for any $y_{1:T}$ that is likely under the probabilistic model $p(y_{1:T})$. Note that the way we parameterize $q_{T}$ may differ from the above approach even though we use the same notation for the parameters $\lambda$.

In this setting instead of the kl divergence in Eq. 3.17, we consider minimizing

[TABLE]

If we let $\lambda_{\eta}(\cdot)$ be a parametric function with parameters $\eta$, we can compute stochastic gradients of Eq. 3.21 w.r.t. $\eta$ with no need to resort to smc. A common choice is to let $\lambda_{\eta}$ be defined by a neural network where $\eta$ are the weights and biases. We define the gradient

[TABLE]

and estimate it, for the current iterate $\eta^{n}$, using mc

[TABLE]

Previous methods we have considered focus on proposals that try to emulate the posterior or the locally optimal proposal, both conditionally on the observed data $y_{1:T}$. However, the amortized inference approach, in this setting, learns proposals based on simulated data from the model $p(x_{1:T},y_{1:T})$. This is performed offline before using the learned proposal and the real dataset for inference.

The amortized inference method follows the same procedure as above in Eq. 3.20 when updating $\eta^{n}$, replacing the gradient with the expression from Eq. 3.23

[TABLE]

Unlike the above approach that uses smc to estimate the gradients, the amortized inference approach results in an unbiased approximation of the gradient. This means that using the update Eq. 3.24 will ensure convergence to a local minima of its cost function Eq. 3.21 by standard stochastic approximation results (Robbins and Monro, 1951). On the other hand, this approach requires us to learn a proposal that works well for any dataset that could be generated by our model. This puts extra stress on the model to be accurate for the actual observed dataset to be able to learn a good proposal. For more thorough discussion of these topics, see Paige and Wood (2016).

We summarize these two approaches to optimize the proposal in Algorithm 6.

Variational Sequential Monte Carlo

The key idea in variational sequential Monte Carlo (vsmc) (Naesseth et al., 2018; Le et al., 2018; Maddison et al., 2017) is to use the parametric distribution for our proposals, $q_{t}(x_{t}\,|\,x_{1:t-1}\,;\,\lambda)$, and optimize the fit in kl divergence from the expected smc approximation $\mathbb{E}\left[\widehat{\gamma}_{T}\right]$ to the target distribution $\gamma_{T}$. The dependence on $\lambda$ enters implicitly in the expected smc approximation $\mathbb{E}\left[\widehat{\gamma}_{T}\right]$ through the proposed samples, weights, and resampling step. In contrast to mimicking the locally optimal proposal distribution, this cost function takes into account the complete smc approximation and defines a coherent global objective function for it. Compared to the previously described adaptive smc methods, vsmc optimizes the fit of the final smc distribution approximation to the true target, rather than the fit between the proposal to the target distribution. This means that we explicitly take into account the resampling steps that are a key part of the smc algorithm.

Studying the marginal distribution of a single sample from $\widehat{\gamma}_{T}$, i.e. $\mathbb{E}\left[\widehat{\gamma}_{T}\right]$, it is possible to show that (Naesseth et al., 2018) the kl divergence from this distribution to the target distribution $\gamma_{T}$ can be bounded from above

[TABLE]

where the log-normalization constant estimate, cf. Eq. 2.18, is

[TABLE]

The expectation is taken with respect to all random variables generated by the smc algorithm. Because $\log Z_{T}$ does not depend on the parameters $\lambda$, minimizing Eq. 3.25 is equivalent to maximizing

[TABLE]

The kl divergence $\textrm{KL}(\mathbb{E}[\widehat{\gamma}_{T}]\|\gamma_{T})$ is non-negative, which means that Eq. 3.27 is a lower bound to the log-normalization constant $\log Z_{T}$. This is why Eq. 3.27 is typically known as an evidence lower bound in the variational inference literature. By maximizing the cost function in Eq. 3.27, with respect to $\lambda$, we can find a proposal that fits the complete smc distribution to the target distribution. The main issue is that evaluating and computing the gradient of this cost function is intractable, so we need to resort to approximations.

We assume that the proposal distributions $q_{t}$ are reparameterizable (Kingma and Welling, 2014; Rezende et al., 2014; Ruiz et al., 2016; Naesseth et al., 2017). This means we can replace simulating $x_{t}\,|\,x_{1:t-1}\sim q_{t}(\cdot\,|\,x_{1:t-1})$ by

[TABLE]

where the distribution of the random variable $\varepsilon_{t}$ is independent of the parameters $\lambda$. If we further assume that $h_{t}$ is differentiable we can use

[TABLE]

where $\nabla_{\lambda}\log\widetilde{w}_{t}(\cdot)$ can be computed using e.g. automatic differentiation, replacing $x_{1:t}^{i}$ with its defintion through Eq. 3.28. The approximation follows from ignoring the gradient part that results from the resampling step, which has been shown to work well in practice (Naesseth et al., 2018; Le et al., 2018; Maddison et al., 2017).

Just like for the adaptive smc methods, Eq. 3.29 suggests a stochastic gradient method to optimize $\lambda$. Given an iterate $\lambda^{n-1}$, we compute the gradient estimate by running smc with proposals $q_{t}(x_{t}\,|\,x_{1:t-1}\,;\,\lambda^{n-1})$ and evaluate

[TABLE]

With stochastic gradient ascent we update our iterate by

[TABLE]

where $\alpha_{n}$ are positive step-sizes, chosen such that $\sum_{n}\alpha_{n}=\infty$ and $\sum_{n}\alpha_{n}^{2}<\infty$.

We summarize the vsmc algorithm to adapt the proposals in Algorithm 7.

3.2 Adapting the Target Distribution

A less well-known design aspect of smc algorithms is the target distribution itself. If we are mainly interested in the final distribution $\gamma_{T}$, we are free to choose the intermediate target distributions $\gamma_{t}$ to maximize the accuracy of our estimate $\widehat{\gamma}_{T}$. By making use of information such as future observations, we can change the target distribution to take this into account. This leads to so-called auxiliary or twisted smc methods (Guarniero et al., 2017; Heng et al., 2017).

3.2.1 Twisting the Target

Even if we could simulate exactly from the locally optimal proposal distribution Eq. 3.1, we still would not be getting exact samples from our target distribution. The reason for this is that when sampling and weighting our particles at iteration $t$ we do not take into account potential future iterations. These future iterations can add new information regarding earlier latent variables. For instance, in an lvm (cf. Eq. 1.9) a natural choice of target distributions is

[TABLE]

However, with this choice we do not take future observations $y_{t+1},\ldots,y_{T}$ into account at iteration $t$. The smc approximation for earlier iterations has finite support (represented by the particles $x_{1:t-1}^{i})$ and so the only way we can incorporate this is by reweighting and resampling. These two operations will typically impoverish our approximation for the full states $x_{1:t}$, a phenomena related to the path degeneracy discussed in Chapter 2. The key idea in twisting (or tilting) the target distribution is to change our intermediate target distributions to take into account information from future iterations already at the current step. These types of approaches also have a connection to so-called lookahead strategies (Lin et al., 2013) and block sampling (Doucet et al., 2006) for smc.

The optimal target distribution, under the assumption that we are using the locally optimal proposal, is to let the target at each iteration be the marginal distribution of the final iteration’s target, i.e.

[TABLE]

With this choice all samples from $x_{1:T}\sim\widehat{\gamma}_{T}$ are perfect samples from $\gamma_{T}$. We can easily see this if we write out the optimal proposal for this sequence of targets

[TABLE]

where the final equality corresponds to the case when $\gamma_{T}$ is the posterior distribution $p(x_{1:T}\,|\,y_{1:T})$ for an lvm. This means that the locally optimal proposal distribution is no longer only locally optimal, it is in fact optimal in a global sense.

In the following discussion we will assume that the final unnormalized target distribution $\widetilde{\gamma}_{T}(x_{1:T})$ can be split into a product of factors:

[TABLE]

With this form of the joint distribution of data and latent variables, we choose the following structure for our unnormalized target distributions:

[TABLE]

where $\psi_{t}\geq 0$ are our so-called twisting potentials, with $\psi_{0}\equiv\psi_{T}\equiv 1$. It is easy to check that this recursive definition is consistent with our assumed factorization of the final target $\widetilde{\gamma}_{T}(x_{1:T})$ in Eq. 3.33. The target distribution structure postulated in Eq. 3.34 can be directly applied in our basic smc method in Algorithm 3 without any changes.

To deduce the optimal twisting potentials $\psi_{t}^{\star}$ we can make use of the property that the optimal targets must fulfill

[TABLE]

This follows from Eq. 3.32, each target $\gamma_{t}^{\star}$ should be the marginal distribution of $x_{1:t}$ under the final target distribution $\gamma_{T}$. By replacing $\widetilde{\gamma}_{t}^{\star}$ in Eq. 3.35 with the definition from Eq. 3.34 we get

[TABLE]

While Eq. 3.36 is typically not available analytically for practical applications, we will see in Section 3.2.2 that it can serve as a guideline for designing tractable twisting potentials.

Below we give a few concrete examples on what the optimal potentials might look like for both non-Markovian lvm s and conditionally independent models with tempering.

Non-Markovian Latent Variable Model

The non-Markovian lvm can be described by a transition pdf $f_{t}$ and an observation pdf $g_{t}$ and follows the structure in Eq. 3.33 by defining

[TABLE]

For this model we can rewrite Eq. 3.36 as follows

[TABLE]

where the final expression is the predictive likelihood of $y_{t+1:T}$ given the latent variables $x_{1:t}$. In Example 3.2.1 we derive the optimal twisting potentials for our Gaussian non-Markovian lvm.

Example 3.2.1** (Analytical Twisting Potential).**

Our running example is a non-Markovian lvm with an analytical optimal twisting potential $\psi_{t}^{\star}$. This is because the joint distribution is Gaussian, and thus the integrals we need to compute to construct it (cf. Eq. 3.37) are tractable. We can rewrite the equation for the observations $y_{t}$ for each $t$ only as a function of $x_{1:l}$, $e_{t}$ and $v_{l+1:t}$, denoted by $y_{t\,|\,l}$, as follows

[TABLE]

where $\gamma_{k}$ is given by

[TABLE]

Because we know that $\psi_{t}^{\star}(x_{1:t})=p(y_{t+1:T}\,|\,x_{1:t})$ we get the optimal value by considering the distribution of $(y_{t+1\,|\,t},\ldots,y_{T\,|\,t})^{\top}$. Since $e_{t}$ and $v_{t}$ are independent Gaussian, by Eq. 3.39 the predictive likelihood is a correlated multivariate Gaussian.

Conditionally Independent Models

Another class of probabilistic models we discussed in Section 1.2 was conditionally independent models. In the notation of Eq. 3.34 we can express this class of models as

[TABLE]

where $p(x)$ is the prior distribution on the (static) latent variable, $s_{t}$ is a backward kernel, and $g_{t}$ is a potential depending on the observed data $\mathbf{y}$. We initialize by $\widetilde{f}_{1}(x_{1})=p(x_{1})g_{1}(x_{1},\mathbf{y})$. The potential $g_{t}$ usually takes either of the two following forms:

[TABLE]

where $0=\tau_{1}<\tau_{2}<\ldots<\tau_{T}=1$. The backward kernel $s_{t}$ is an artificially introduced distribution to enable the use of smc methods, requiring approximating a space of increasing dimension. See Section 1.2 for more examples of useful sequence models in annealing methods, or Del Moral et al. (2006a) for a more thorough treatment of the subject.

When we replace $\widetilde{f}_{t}$ in Eq. 3.36 with its definition from Eq. 3.40 we get

[TABLE]

where we can see that because of the (reverse) Markov structure in the annealing model from Eq. 3.40, the optimal twisting potential only depends on the current value of $x_{t}$, i.e. $\psi_{t}^{\star}(x_{1:t})=\psi_{t}^{\star}(x_{t})$. This will hold true also for the ssm, the Markovian special case of Eq. 3.37 where $f_{t}(x_{t}\,|\,x_{1:t-1})=f_{t}(x_{t}\,|\,x_{t-1})$ and $g_{t}(y_{t}\,|\,x_{1:t},y_{1:t-1})=g_{t}(y_{t}\,|\,x_{t})$.

3.2.2 Designing the Twisting Potentials

The optimal twisting potentials are typically not tractable, requiring us to solve integration problems that might be as difficult to compute as finding the posterior itself. However, they can give insight into how to design and parameterize approximate twisting functions. First, we study a certain class of models that admit a locally optimal twisting potential, which gives rise to the so-called fully adapted smc (Pitt and Shephard, 1999; Johansen and Doucet, 2008). After discussing the locally optimal choice, we then move on to discuss a general approach to learning twisting potentials for smc algorithms.

Locally Optimal Twisting Potential

If we assume that the twisting potentials only satisfy the optimal twisting equation, Eq. 3.36, for one iteration we get a locally optimal twisting potential

[TABLE]

The target distribution based on these potentials will allow us to adapt for information from one step into the future.

Furthermore, we combine the locally optimal twisting potentials with the proposals based only on factors up until iteration $t$, i.e.

[TABLE]

With these choices of twisting functions and proposals, based on the twisted target structure Eq. 3.34, we obtain the fully adapted smc (Pitt and Shephard, 1999; Johansen and Doucet, 2008) with corresponding weights

[TABLE]

Note that this perspective on fully adaptive smc approximates the “predictive target”

[TABLE]

including information from iteration $t+1$ already at iteration $t$. In some situations we might be interested in the “filtering target” $\prod_{k=1}^{t}\widetilde{f}_{k}(x_{1:k})$. If so, we can make use of the approach described above to generate the particles $x_{1:t}^{i}$. Then to estimate the filtering target we simply average the particles $x_{1:t}^{i}$, ignoring the weights in Eq. 3.44, accounting for the discrepancy between the filtering and predictive targets.

Example 3.2.2** (Locally Optimal Twisting for Non-Markovian LVM).**

Recall the non-Markovian lvm with factors given by (see Eq. 3.37) $\widetilde{f}_{t}(x_{1:t})=f_{t}(x_{t}\,|\,x_{1:t-1})g_{t}(y_{t}\,|\,x_{1:t},y_{1:t-1})=p(y_{t},x_{t}\,|\,x_{1:t-1},y_{1:t-1}).$ For this model the locally optimal twisting potential Eq. 3.42 is given by

[TABLE]

The proposal in Eq. 3.43 becomes

[TABLE]

with weights according to Eq. 3.44 given by $\widetilde{w}_{t}(x_{1:t})=p(y_{t+1}\,|\,x_{1:t}).$ This can be recognized as a one-step look-ahead strategy, incorporating the observation $y_{t+1}$ already at iteration $t$. This is indeed how the fully adaptive smc is often presented for ssms and non-Markovian lvms.

Remark 3.2.1**.**

The proposal used above is not the locally optimal proposal distribution for our (predictive) target distribution. The locally optimal proposal (see Eq. 3.1) for the target Eq. 3.34 is

[TABLE]

which would result in weights

[TABLE]

Just like in the non-twisted case the weights for the locally optimal proposal are independent of the actual samples we generate at iteration $t$ of the smc algorithm. However, even if the potentials in Eq. 3.42 are available, the integration in Eq. 3.45 is not necessarily tractable.

Except for a few special cases the locally optimal twisting potential is not available, therefore we resort to approximating it. We can extend the ideas from Section 3.1 and either make an approximation to the locally optimal potential, or we can learn a full twisting potential.

Learning Twisting Potentials

By parameterizing a class of positive functions and a cost function, we can effectively learn useful twisting potentials for the application at hand. From Eq. 3.36 we know that the optimal twisting potentials must satisfy a recursive relationship. We can make use of this property as we define the cost function and the twisting potentials.

We explain the approach by Guarniero et al. (2017), but generalize slightly to allow for non-Markovian lvm s. This means that our model is defined by $\widetilde{f}_{t}=f_{t}(x_{t}\,|\,x_{1:t-1})g_{t}(y_{t}\,|\,x_{1:t})$, cf. Eq. 3.37. Furthermore, we define the twisting potentials $\psi_{t}(x_{1:t})$ implicitly as expectations with respect to $f_{t+1}$ of positive functions $\bar{\psi}_{t+1}(x_{1:t+1}\,;\,\rho_{t+1})$ with parameters $\rho_{t+1}$, i.e.

[TABLE]

We further assume that this expectation can be computed analytically. With this setup we focus on learning the $\bar{\psi}_{t}$’s instead. Using the optimality condition from Eq. 3.36 with twisting potentials defined by Eq. 3.46, we get the recursive approximation criteria for the $\bar{\psi}_{t}$’s

[TABLE]

Because we assumed that $\mathbb{E}_{f_{t}(x_{t}\,|\,x_{1:t-1})}\left[\bar{\psi}_{t}(x_{1:t}\,;\,\rho_{t})\right]$ is available analytically, and Eq. 3.47 is telling us to view it as an approximation of $p(y_{t:T}\,|\,x_{1:t-1})$, it makes sense to use as the proposal

[TABLE]

With this choice of proposal we tie together the parameters for the proposal and the potentials. The corresponding weights in the smc algorithm are

[TABLE]

What remains is how to use Eq. 3.47 to actually learn a useful set of potentials. To do this we follow the approach by Guarniero et al. (2017) which uses an iterative refinement process: initialize parameters, run smc with the current parameters, update parameters based on the current smc approximation $\widehat{\gamma}_{t}$, repeat. The update step solves for $\rho_{t}^{n}$ by recursively minimizing

[TABLE]

where the last step is for $t=T-1,\ldots,1$. With the updated set of parameters, we re-run smc and repeat the updates in Eq. 3.50 if a stopping criteria has not been met. Guarniero et al. (2017) proposes a stopping criteria based on the normalization constant estimate $\widehat{Z}_{T}$.

We summarize the iterated auxiliary smc method explained above in Algorithm 8, and give an example setup for a stochastic rnn in Example 3.2.3.

Example 3.2.3** (Stochastic Recurrent Neural Network).**

A stochastic rnn is a non-Markovian lvm where the parameters of the transition and observation models are defined by rnn s. A common example is using the conditional Gaussian distribution to define the transition pdf

[TABLE]

where the functions $\mu_{t}(\cdot),\Sigma_{t}(\cdot)$ are defined by rnn s. This model together with a Gaussian-like definition for the potentials satisfies the criteria for using Algorithm 8. We define $\bar{\psi}_{t}$ as follows

[TABLE]

where the functions $\Lambda_{t}(\cdot),\iota_{t}(\cdot),c_{t}(\cdot)$ depend on the parameters $\rho_{t}$—for notational brevity we have not made this dependence explicit. The proposal (cf. Eq. 3.48) is given by

[TABLE]

and the twisting potential

[TABLE]

Twisting the target distribution is an area of active research (Guarniero et al., 2017; Heng et al., 2017; Lawson et al., 2018; Lindsten et al., 2018). The iterated auxiliary smc (Guarniero et al., 2017) makes strong assumptions on the model for easier optimization. End-to-end optimization of a priori independent twisting potentials $\psi_{t}$ and proposals $q_{t}$ might lead to even more accurate algorithms.

Appendix 3.A Taylor and Unscented Transforms

In this appendix we detail the ekf- and ukf-based approximations to the distribution $p(x_{t},y_{t}\,|\,x_{1:t-1})$. We refer to the two approaches as Taylor and Unscented transforms, respectively.

Taylor Transform

If we assume that $v_{t}$ and $e_{t}$ are zero-mean and linearize $a$ around $v_{t}=0$ and $c$ around $x_{t}=a(x_{1:t-1},0),e_{t}=0$ we get

[TABLE]

where $\bar{a}=a(x_{1:t-1},0)$, $\bar{c}=c\left((x_{1:t-1},\bar{a}),0\right)$, and $J_{v_{t}}^{a}(\bar{a})$ denotes the Jacobian of the function $a$ with respect to $v_{t}$ evaluated at the point $\bar{a}$. We rewrite Eq. 3.51:

[TABLE]

For the approximation in Eq. 3.52 $x_{t},y_{t}\,|\,x_{1:t-1}$ is indeed Gaussian, and with $v_{t}\sim\mathcal{N}(0,Q),e_{t}~{}\sim\mathcal{N}(0,R)$, we can identify the blocks of $\widehat{\mu}$ and $\widehat{\Sigma}$:

[TABLE]

Using e.g. automatic differentiation we can easily compute the Jacobians needed for the ekf-based proposal. However, we still have to be able to evaluate the densities defined by $a(x_{1:t-1},v_{t})$ and $c(x_{1:t},e_{t})$ for the weight updates $\widetilde{w}_{t}$. For notational convenience we have omitted the dependence on previous latent states in the expressions above. However, in general because of this dependence we will have to compute the proposals separately for each particle, which includes inverse matrix computations that can be costly if $x_{t}$ is high-dimensional.

Unscented Transform

We rewrite Eq. 3.3 as a single joint random variable, and such that the functions $a,c$ only depend on $x_{1:t-1},v_{t}$, and $e_{t}$,

[TABLE]

where $z_{t}=(v_{t},e_{t})^{\top}$. By approximating the conditional distribution of $(x_{t},y_{t})^{\top}$ given $x_{1:t-1}$ as Gaussian, we can again compute the distribution of $x_{t}\,|\,x_{1:t-1},y_{t}$ for the approximation. We choose sigma points based on the two first moments of $z_{t}$. The mean and variance of $z_{t}$ is

[TABLE]

where the scalars $\sigma_{l}$ and vectors $\mathrm{u}_{l}$ correspond to the singular value decomposition of $\Sigma_{z}$. We then choose the $2n_{x}+2n_{y}+1$ sigma points for $z$ as follows,

[TABLE]

where $\lambda=\alpha^{2}(n_{x}+n_{y}+\kappa)-n_{x}-n_{y}$ (Julier and Uhlmann, 1997). We will discuss the choice of design parameters $\alpha$ and $\kappa$ below. The mean and variance of $(x_{t},y_{t})^{\top}$ is estimated by

[TABLE]

where the coefficients $\mathrm{\omega}^{l}$ are

[TABLE]

and $\beta$ is another design parameter.

Like in the ekf-based approximation we still need access to the densities corresponding to Eq. 3.3. Unlike in the ekf case the, ukf-based approximation does not need access to derivatives. However, we need to choose the three design parameters: $\alpha,\beta,\kappa$. Both $\alpha$ and $\kappa$ control the spread of the sigma points, and $\beta$ is related to the distribution of $z_{t}$ (Julier, 2002). Typical values of these parameters are $(\alpha,\kappa,\beta)=(10^{-3},0,2)$.

Chapter 4 Nested Monte Carlo: Algorithms and Applications

In this chapter we will discuss nesting Monte Carlo (mc)algorithms, using several layers of mc to construct composite algorithms. Usually, this is done to approximate some intractable idealized algorithm. The intractability can stem from e.g. a likelihood or weight that we cannot evaluate in closed form, or a proposal distribution that we cannot sample exactly from. Either way, the key idea of these methods is to introduce additional auxiliary variables to tackle the intractability in the idealized method. The approximate method is constructed in such a way that it retains some of the favorable properties of the idealized mc method that it tries to emulate. We will in this section study particle-based pseudo marginal (pm)methods for several idealized mc methods, ranging from marginal Metropolis-Hastings (mh) (Andrieu and Roberts, 2009; Andrieu et al., 2010) to distributed smc (Vergé et al., 2015) and various nested smc methods (Fearnhead et al., 2010; Chopin et al., 2013; Naesseth et al., 2015).

In the first section below we discuss the unbiasedness property of the smc normalization constant estimate, a key property underpinning several nested mc methods. Then, in the next section we discuss particle Metropolis-Hastings (pmh). Next, we introduce the concept of proper weights and discuss various smc methods relying on random (but proper) weights. Finally, we introduce a few perspectives on distributed smc.

4.1 The Unbiasedness of the Normalization Constant Estimate

We briefly touched upon one of the key theoretical properties of smc in Section 2.3, the unbiasedness of the normalization constant estimate $\widehat{Z}_{t}$. This will be important for several of the pm methods that we discuss in this chapter. We repeat the definition of $\widehat{Z}_{t}$ again for clarity:

[TABLE]

In Proposition 4.1.1 we provide a result that formalizes the unbiasedness property. Furthermore, we provide a straightforward self-contained proof of this important property of the smc algorithm.

Proposition 4.1.1**.**

The smc normalization constant estimate $\widehat{Z}_{t}$ is non-negative and unbiased

[TABLE]

Proof.

See Appendix 4.A. ∎

Often it will prove useful to think of the $\widehat{Z}_{t}$ simply as a function of the random seed used to generate all the random variables in the smc algorithm. Below we formalize this point-of-view, which will be useful for e.g. pmh and nested sequential Monte Carlo (nsmc) in Section 4.2 and Section 4.3, respectively.

Auxiliary Variables and Random Seeds

The complete set of random variables generated in the smc algorithm consists of each of the particles $x_{t}^{i}$, for $i\in\{1,\ldots,N\}$ and $t\in\{1,\ldots,T\}$, and ancestor variables $a_{t}^{i}$ for $i\in\{1,\ldots,N\}$ and $t\in\{1,\ldots,T-1\}$. Both $x_{t}^{i}$ and $a_{t}^{i}$ are often generated by sampling base auxiliary random numbers such as Gaussian or uniforms, then the particles and ancestor variables are simply a function of these auxiliary variables. We will denote the collection of all base auxiliary random variables, also referred to as the random seed, by $\mathbf{u}$. Each individual particle $x_{t}^{i}$ and ancestor $a_{t}^{i}$ are then a transformation from the random seed. For example we might have $\mathbf{u}=u_{x,1:T}^{1:N}\cup u_{a,1:T-1}^{1:N}$, where each $u_{x,t}^{i},u_{a,t}^{i}\sim\mathcal{U}(0,1)$. The function that connects $u_{x,t}^{i}$ to $x_{t}^{i}$ and $u_{a,t}^{i}$ to $a_{t}^{i}$ are then the inverse cumulative distribution functions (cdfs) of the proposal $q_{t}$ and resampling mechanism from Section 2.2.1, respectively. A common tool in pm is to think of the normalization constant estimate in Eq. 4.1 as a function of all the random variables $\mathbf{u}$, effectively the random seed of the smc algorithm. Then we can write Eq. 4.1 as $\widehat{Z}_{t}(\mathbf{u})$, i.e.

[TABLE]

where the weights $\widetilde{w}_{k}^{i}$ are functions of the random seed $\mathbf{u}$. We denote the distribution of the random seed $\mathbf{u}$ by $\tau(\mathbf{u})$. The main result in Proposition 4.1.1 can then be written in terms of the random seed $\mathbf{u}$: $\mathbb{E}_{\tau(\mathbf{u})}\left[\widehat{Z}_{t}(\mathbf{u})\right]=Z_{t}$.

Example 4.1.1** ($\widehat{Z}_{t}(\mathbf{u})$ for Example 1.2.1).**

We revisit our running non-Markovian Gaussian example and consider a standard smc estimator with the prior proposal and no twisting of the target. This means that for $\mathbf{u}=\varepsilon_{1:T}^{1:N}\cup u_{1:T-1}^{1:N}$, where $\varepsilon_{t}^{i}\sim\mathcal{N}(0,1)$ and $u_{t}^{i}\sim\mathcal{U}(0,1)$, we have

[TABLE]

where $x_{1:k}^{i}=(x_{1:k-1}^{a_{k-1}^{i}},x_{k}^{i})$, the ancestor variable $a_{k-1}^{i}$ is defined by the resampling mechanism Section 2.2.1 using $u_{k-1}^{i}$, and

[TABLE]

In LABEL:code:rs:running we provide a code snippet, definitions left out for clarity, that illustrates how $\widehat{Z}_{t}(\mathbf{u})$ can be computed. For numerical stability purposes we compute the log-normalization constant estimate, $\log\widehat{Z}_{t}(\mathbf{u})$, instead.

4.2 Particle Metropolis-Hastings

Pmh melds particle-based methods, such as is and smc, with Metropolis-Hastings, providing the practitioners with powerful approximate Bayesian inference for a large class of intractable likelihood problems. The particle mcmc methods discussed in this section were mainly developed in Andrieu et al. (2010).

Mh is a mcmc method (Robert and Casella, 2004) for approximate simulation from a pdf. The typical example is the posterior distribution of unknown model parameters $\theta$ given the observed data $\mathbf{y}$,

[TABLE]

where here $\mathbf{x}$ are latent (unknown) variables that we wish to integrate out. The likelihood $p(\mathbf{y}\,|\,\theta)=\int p(\mathbf{x},\mathbf{y}\,|\,\theta)\mathop{}\!\mathrm{d}\mathbf{x}$ is often intractable. A common solution is to use data augmentation, sampling instead (approximately) from the posterior distribution of both $\theta$ and $\mathbf{x}$ given the data $\mathbf{y}$, i.e. $p(\theta,\mathbf{x}\,|\,\mathbf{y})$. Marginal mh, see Algorithm 9, on the other hand generates approximate samples directly from the marginal distribution $p(\theta\,|\,\mathbf{y})$, provided that the likelihood can be evaluated exactly. The key idea is to generate a Markov chain $\theta_{1},\theta_{2},\ldots$ with a stationary distribution equal to the target distribution of interest Eq. 4.4. The accept-reject step in Algorithm 9, accepting the proposed value $\theta^{\prime}$ as the new iterate with probability $A$, ensures that the stationary distribution is equal to Eq. 4.4. The acceptance probability

[TABLE]

is the minimum of one and the ratio of the posterior and the proposal evaluated at the proposed value, divided by the same ratio evaluated at the previous iterate. For a more thorough treatment of the mh algorithm see e.g. Robert and Casella (2004).

4.2.1 Particle Marginal Metropolis-Hastings

The pseudo marginal method replaces the exact likelihood with an estimate of it, usually constructed by is or smc like in Eq. 4.1. For a fixed value of the parameter $\theta$ we can use the tools from the previous chapters to construct a non-negative and unbiased approximation to $p(\mathbf{y}\,|\,\theta)$, we will use the notation from (4.3) and label this estimate $\widehat{Z}(\mathbf{u},\theta)$. From Proposition 4.1.1 we have that

[TABLE]

Simply replacing $p(\mathbf{y}\,|\,\theta)$ with this unbiased and non-negative approximation, $\widehat{Z}(\mathbf{u},\theta)$, results in the particle marginal Metropolis-Hastings (pmmh) algorithm in Algorithm 10 (Andrieu et al., 2010). In a practical implementation we do not need to save the entire seed $\mathbf{u}$ for each iteration, it suffices to save the scalar value for our normalization constant estimate $\widehat{Z}(\theta,\mathbf{u})$. We illustrate the mh and pmmh algorithms applied to our running example in Fig. 4.1.

Example 4.2.1** (Particle marginal Metropolis-Hastings for Example 1.2.1).**

We revisit our running non-Markovian Gaussian example. To keep it simple we focus on the is special case, i.e. with $T=1$, which means we are interested in inference for $\theta=(q,r)$. We assign a log-normal prior to the variance parameters $q$ and $r$, i.e. $\log q\sim\mathcal{N}(0,1)$ and $\log r\sim\mathcal{N}(0,1)$. The corresponding marginal likelihood is then given by $p(\mathbf{y}\,|\,\theta)=\mathcal{N}(\mathbf{y}\,|\,0,q+r)$. We use a normal distribution random walk proposal with variance $0.5$ for the logarithm of the parameters. In Fig. 4.1 we can see the kernel density estimate based on the samples of a $20\,000$ iterations long run of the two respective algorithms. The two algorithms obtain the same marginal distribution as expected.

The motivation for why this works is fairly straightforward. We can think of the pmmh algorithm as a standard mh on the extended space of $(\theta,\mathbf{u})$, i.e. adding the random seed as an auxiliary variable. The proposal on this extended space is

[TABLE]

and the distribution we are trying to draw samples from, i.e. the target distribution, is

[TABLE]

which by the unbiasedness of $\widehat{Z}(\mathbf{u},\theta)$ has the posterior distribution as a marginal

[TABLE]

This means that if we can generate samples approximately distributed as the distribution in Eq. 4.7, we get samples approximately distributed as $p(\theta\,|\,\mathbf{y})$ as a by-product by considering only the samples for $\theta$.

To show that pmmh is a standard mh on an extended space we first note that the proposal $q(\theta^{\prime}\,|\,\theta)\tau(\mathbf{u}^{\prime})$ is the distribution of the proposed value $(\theta^{\prime},\mathbf{u}^{\prime})$. Then studying the corresponding acceptance probability $A$ for the above choice of proposal and target distributions,

[TABLE]

we get exactly the expression from Algorithm 10. These two properties ensure that pmmh is an mh algorithm on the extended space with Eq. 4.7 as its stationary distribution.

4.2.2 Sampling the Latent Variables

If we are interested not just in the parameters but in the latent variables also, we can adapt the particle mh framework with a simple extension of the pmmh algorithm (Andrieu et al., 2010). We can use particle methods to generate approximate samples from the complete target with both parameters and latent variables, $p(\theta,\mathbf{x}\,|\,\mathbf{y})$. We simply use the final smc approximation for $p(\mathbf{x}\,|\,\theta,\mathbf{y})$ that we obtain as a by-product of computing $\widehat{Z}(\mathbf{u},\theta)$ as a proposal in the mh algorithm. The smc approximation to the conditional distribution $\mathbf{x}\,|\,\theta,\mathbf{y}$ is

[TABLE]

where we suppress the dependence on the parameters $\theta$ and data $\mathbf{y}$ for clarity. Just like $\widehat{Z}(\mathbf{u},\theta)$ the conditional distribution estimator is also completely defined by the random seed $\mathbf{u}$. Using the above distribution as our proposal for $\mathbf{x}$ lets us derive the particle Metropolis-Hastings algorithm, see Algorithm 11.

Clearly, this algorithm can be applied for sampling the latent variables also when the model parameters $\theta$ are all known. In this case we simply drop all $\theta$-dependent terms and variables in Algorithm 11. The proposal on the extended space is then

[TABLE]

which corresponds to first running the smc algorithm (sample $\mathbf{u}^{\prime}\sim\tau(\mathbf{u})$) and then sample one particle trajectory from the resulting smc approximation (sample $\mathbf{x}^{\prime}\sim\widehat{p}(\mathbf{x}\,|\,\mathbf{u}).$ Since this is done independently from the current state of the Markov chain, it is referred to as an independent proposal and the resulting algorithm is know as particle independent Metropolis-Hastings (pimh) (Andrieu et al., 2010). The acceptance probability for pimh is simply the minimum of one and the ratio of normalization constant estimators, i.e.

[TABLE]

Example 4.2.2** (Particle independent Metropolis-Hastings for Example 1.2.1).**

We revisit our running non-Markovian Gaussian example. We set $T=20$ and run pimh for $N=5$, $10$ and $20$. The kernel density estimates of the marginal posterior distributions of $x_{1}$ and $x_{T}$ based on the $10,000$ samples generated can be seen in Fig. 4.2. We can see that, while in theory exact in the limit of infinite number of samples, the approximations for the lower number of particles $N=5$, $10$, are distinctly non-Gaussian. It is only when we use a large enough number of particles, e.g. $N=20$, that we obtain more satisfactory results for a finite number of mcmc iterations.

Analogously to pmmh, we can motivate pmh as a standard mh algorithm on an extended space. We focus on the the setting where the model parameters are known and the extended space is given by $(\mathbf{u},\mathbf{x})$. The proposal on this extended space is

[TABLE]

and the target distribution is

[TABLE]

By the unbiasedness property of smc, i.e. Theorem 2.3.1, we have that the posterior distribution $p(\mathbf{x}\,|\,\mathbf{y})$ is a marginal of Eq. 4.11

[TABLE]

where $\stackrel{{\scriptstyle\textrm{d}}}{{=}}$ denotes equal to in distribution. Studying the corresponding acceptance probability $A$ for the above proposal and target distribution we obtain

[TABLE]

which is identical to the expression from Algorithm 11 (when the model parameters are known).

4.2.3 Correlated Pseudo-Marginal Methods

One of the main limitations of the pm methods that we have covered so far is that they require that the normalization constant estimate $\widehat{Z}$ is sufficiently accurate to achieve good performance. This requirement means we have to focus much of the computational effort on the smc algorithm by using a high number of particles $N$, rather than generating more samples from the posterior $p(\theta\,|\,\mathbf{y})$. By studying the pmmh algorithm we can see that the proposed parameter $\theta^{\prime}$ at each iteration is correlated with the previous parameter $\theta^{j-1}$, encouraging local exploration and higher acceptance rates. However, this is not true for the random seed $\mathbf{u}$ that is generated independently at each iteration. Introducing correlation also for this random seed, i.e. replacing the independent proposal $\tau(\mathbf{u})$ with a correlated proposal $\tau(\mathbf{u}^{\prime}\,|\,\mathbf{u})$, can allow us to effectively lower the computational cost of the smc estimate through choosing a lower number of particles $N$. Since the dimension of $\mathbf{u}$ typically increases linearly with $N$, we need a proposal distribution which scales favorably to high dimensions. Examples include preconditioned Crank-Nicolson (Deligiannidis et al., 2018; Dahlin et al., 2015), elliptical slice sampling (Murray and Graham, 2016) and Hamiltonian Monte Carlo (Lindsten and Doucet, 2016).

4.3 Proper Weights and Sequential Monte Carlo

So far in this manuscript we have focused on smc algorithms where the particles are sampled exactly from the proposal $q_{t}$ and where the weights $\widetilde{w}_{t}$ are deterministic functions of the particles. We can in fact relax both of these restrictions significantly by considering the weights and particles as joint random variables that satisfy a proper weighting condition. This will lead us to powerful algorithms such as the random-weights particle filter (Fearnhead et al., 2010), smc 2 (Chopin et al., 2013), and nsmc methods (Naesseth et al., 2015, 2019).

To formally justify these algorithms we can make use of the concept of proper weighting. This is a property on the expected value of the weights and samples when applied to functions.

Definition 4.3.1** (Proper Weighting).**

We say the (random) pair $(x,\widetilde{w})$ are properly weighted for the (unnormalized) distribution $\widetilde{\gamma}$ if $\widetilde{w}\geq 0$ a.s. and for all measurable functions $h$

[TABLE]

for some positive constant $\mathcal{C}>0$ that is independent of $x$ and $\widetilde{w}$.

This is highly related to the unbiasedness property of the smc method that we saw already in Theorem 2.3.1. However, while the basic smc has deterministic weights when conditioned on the samples, here we allow also for random weights within the confines of Definition 4.3.1. See also Liu (2004) for an early discussion on this topic in the context of random weights and is.

We will denote the distribution of the random pair $(\widetilde{w},x)$ by $Q(\mathop{}\!\mathrm{d}\widetilde{w},\mathop{}\!\mathrm{d}x)$. By simulating from $Q$ independently

[TABLE]

we can construct an estimate of $\mathbb{E}_{\gamma}\left[h(x)\right]$ that satisfies the law of large numbers,

[TABLE]

Proper weighting ensures that the numerator converges to $Z\mathcal{C}\mathbb{E}_{\gamma}[h(x)]$ and the denominator to $Z\mathcal{C}$, and thus the right hand side converges to the true expectation $\mathbb{E}_{\gamma}[h(x)]$. Recall that the constant $\mathcal{C}$ is as defined in Definition 4.3.1.

We have already seen one example of a properly weighted is method in Section 3.1, namely nested is (Algorithm 5). Because standard is is a special case of nested is it follows that also standard is leads to weight-sample pairs that are properly weighted. Below we will see how this can be used to motivate more powerful algorithms for approximate inference.

4.3.1 Random Weights Sequential Monte Carlo

Random weights smc and is focuses on solving the problem with intractable weights. Because we pick the proposal ourselves it is often easy to sample from, whereas the target distribution depends largely on the model we want to study and can lead to intractable importance weights. A concrete example of this is partially observed stochastic differential equations (sdes) as studied by e.g. Beskos et al. (2006); Fearnhead et al. (2010). In this section we focus exclusively on random weights smc algorithms for sdes. The corresponding target distribution, and thus the weights, contains an exponentiated integral that is intractable. If we replace the intractable weights in the standard smc framework, $\widetilde{w}(x)$, with (non-negative) unbiased approximations of them, $\widehat{\widetilde{w}}(x)$, we obtain random weights smc. While we focus the exposition on the basic is case, the result and method follows for smc when we have unbiased estimates of the weights conditional on previous samples. It is straightforward to prove that the random pair $(\widehat{\widetilde{w}}(x),x)$ are properly weighted for $\widetilde{\gamma}$,

[TABLE]

A simplified version of the target distribution in an sde has the following form of its intractable part

[TABLE]

where $\Phi(x)$ is an intractable functional, often an integral. We assume that we can generate unbiased approximations $\widehat{\Phi}(x)$ of the functional for a given value of $x$. Using the Poisson estimator (Wagner, 1987) we can obtain an unbiased estimate of the exponentiated functional in Eq. 4.15.

[TABLE]

where the expectation is with respect to $\kappa\sim\textrm{Poisson}(\lambda)$. A mc estimate of Eq. 4.16 leads to the Poisson estimator

[TABLE]

where $\widehat{\Phi}^{j}(x)$ are iid unbiased estimators of $\Phi(x)$. While Eq. 4.17 is unbiased, it is not necessarily non-negative and might be unsuitable for weights estimation in an smc algorithm. With proper modification to the estimate we can sometimes ensure that the estimate is non-negative (Fearnhead et al., 2010; Jacob and Thiery, 2015). Replacing the intractable part in the weights with its unbiased estimate, like Eq. 4.17 in the sde case, results in the random weights smc algorithm Algorithm 12. The samples $x_{t}^{i}$ are generated by a standard user chosen proposal $q_{t}(x_{t}\,|\,x_{1:t-1})$ just like in standard smc.

Rather than interpreting random weights smc as a properly weighted smc algorithm, we can instead motivate it as a standard smc on an extended space.

4.3.2 Nested Sequential Monte Carlo

We already briefly touched on nsmc in Section 3.1, illustrating how a nested is sampler can be used to approximate the locally optimal proposal distribution. Here we will instead consider the case when the nested sampler is itself smc rather than is. This usually allows the user to derive more efficient proposal approximations, pushing the boundaries of the type of problems that smc methods can be effectively applied to.

Nsmc (Naesseth et al., 2015, 2019) is specifically targeting models where the latent state $x_{t}$ is of dimension $n_{x}>1$ (typically in the order of 10’s to 100’s). The goal of nsmc is to approximate the locally optimal proposal distribution by running a separate internal (or nested) smc algorithm over the components of $\mathbf{x}_{t}=(x_{t,1},\ldots,x_{t,n_{x}})$ for each sample we need. At first this might seem wasteful of our computational resources, running a nested smc sampler with $M$ internal particles for each particle in the outer algorithm results in $\mathcal{O}(NM)$ operations. However, for many models it results in much more accurate estimate than, for instance, a corresponding standard smc algorithm with $NM$ particles and the prior proposal.

Unlike in the random weights smc, nsmc changes not only the way we compute weights, but also the way we generate the samples $x_{t}^{i}$. The key idea is to construct a smc approximation to the locally optimal proposal $q_{t}^{\star}(x_{t}\,|\,x_{1:t-1})$. For each particle $x_{t}^{i}$ we want to simulate we construct an smc approximation to $q_{t}^{\star}$

[TABLE]

where $(\widetilde{\nu}_{n_{x}}^{j,i},\bar{x}_{1:n_{x}}^{j,i})$ are properly weighted for $q_{t}^{\star}(x_{t}\,|\,x_{1:t-1}^{i})$. We assume that Eq. 4.18 is constructed using a nested smc sampler on the components of $x_{t}$; first targeting $x_{t,1}$, then $x_{t,1:2}$, etc, where the final target is $q_{t}^{\star}(x_{t,1:n_{x}}\,|\,x_{1:t-1}^{i})$. With this approximation nsmc replaces the sampling step Eq. 2.15 by simulating from Eq. 4.18, and the weighting step Eq. 2.16 by

[TABLE]

where $\widetilde{\nu}_{k}^{j,i}$ are the corresponding weights in the nested smc sampler. We summarize nsmc in Algorithm 13 and give examples below of ways to construct $\widehat{q}_{t}^{\star}(x_{t}\,|\,x_{1:t-1}^{i})$.

Example 4.3.1** (Nested sequential Monte Carlo for Example 1.2.1).**

We revisit our running non-Markovian Gaussian example. Because $x_{t}$ is scalar for this example, we make a straightforward extension to the model to illustrate nsmc. We replicate the model from Example 1.2.1 $n_{x}$ times, i.e.

[TABLE]

To construct $\widehat{q}_{t}^{\star}(x_{t}\,|\,x_{1:t-1}^{i})$ we run an internal smc sampler with $M$ samples targeting

[TABLE]

for $l=1,\ldots,n_{x}$. We denote the particles and weights for the nested smc by $\bar{x}_{1:l}^{j,i}$ and $\widetilde{\nu}_{l}^{j,i}$, respectively. A straightforward proposal for the above target distribution is $\bar{x}_{1:l}^{j,i}\sim\widehat{\eta}_{l-1}(x_{t,1:l-1}\,|\,x_{1:t-1}^{i})f(x_{t,l}\,|\,x_{t-1,l}^{i})$, i.e. using the prior proposal. We have that $\widehat{\eta}_{l}$ is

[TABLE]

where the weights are given by

[TABLE]

We obtain the approximation to the locally optimal proposal $\widehat{q}_{t}^{\star}(x_{t}\,|\,x_{1:t-1}^{i})\equiv\widehat{\eta}_{n_{x}}(x_{t,1:n_{x}}\,|\,x_{1:t-1}^{i})$. This defines the nsmc method in Algorithm 13 for the multivariate extension to our running example. We illustrate the nested target distributions in Fig. 4.3.

Just like in the standard smc sampler, there is room for improvement in the nested smc by e.g. choosing a better internal proposal distributions and target distributions under the proper weighting constraint. It is possible to further improve performance by correlating samples through a combination of fully adapted smc and a forward-filtering backward-sampling procedure as discussed by Naesseth et al. (2015, 2019).

4.4 Distributed Sequential Monte Carlo

Developing computational methods that can take advantage of parallel and multi-core computing architectures is important for efficient approximate inference algorithms. In this section we discuss how the idea of nesting mc algorithms can open up for distributed implementation of smc. First we discuss independent importance weighted smc samplers, a way to combine the outputs from multiple smc samplers using an outer level importance sampler. Then we discuss the island particle filter (Vergé et al., 2015), where the outer level sampling is done using smc as well. It should be noted that the purpose of this section is to illustrate how the idea of nested mc opens up for parallelization of smc, and not to provide a comprehensive overview of distributed smc algorithms. Indeed, many useful methods have been proposed for parallelizing the smc algorithm itself. Often this is done at the granularity of individual particles. The key challenge is then to parallelize the resampling step, since this step requires interaction between particles. One way to address this challenge is to make use of mh or rejection sampling to implement the resampling (Murray et al., 2016). Another approach is to modify the resampling step in order to limit the interaction between particles, thereby making a distributed implementation easier (Whiteley et al., 2016).

4.4.1 Importance Weighted Independent SMC Samplers

Assume that we have access to $M$ computing nodes. These could correspond to different machines, different cores, or even $M$ consecutive runs on the same machine. A (very) simple way of distributing computation across these nodes is to run separate independent smc samplers on each node. Using the notation from above, let $\mathbf{u}^{m}\sim\tau(\mathbf{u})$ denote all the random variables generated during the run of an smc sampler with $N$ particles in node $m\in\{1,\dots,M\}$. The variable $\mathbf{u}^{m}$ can also be thought of as the random seed used to initialize the $m$:th sampler. We then get $m$ independent approximations of the target $\gamma_{T}$ given by

[TABLE]

where we have assumed that each node uses the same number of particles $N$ and the same smc algorithm. Assume that we want to compute the expected value $\mathbb{E}_{\gamma_{T}}[h(x_{1:T})]$ of some test function $h$. Each node provides an estimate,

[TABLE]

and since these estimates are independent and identically distributed a natural approach is to take the average $\frac{1}{M}\sum_{m=1}^{M}\widehat{h}(\mathbf{u}^{m})$ as the final estimate. However, since each smc estimate is biased for finite $N$, the resulting average will also be biased regardless of the number of nodes $M$ used in the computation. A better approach can therefore be to weight the independent smc samplers according to their respective normalizing constant estimates:

[TABLE]

This aggregated estimator is obtained by viewing the smc sampler in each node as a proposal for an outer-level importance sampler. Specifically, the target distribution for this importance sampler is given by

[TABLE]

and the proposal is given by $\tau(\mathbf{u})$, resulting in weights $\Pi(\mathbf{u})/\tau(\mathbf{u})\propto\widehat{Z}_{T}(\mathbf{u})$. By computations analogous to those in Section 4.2 it holds that $\mathbb{E}_{\Pi}[\widehat{h}(\mathbf{u})]=\mathbb{E}_{\gamma}[h(x_{1:T})].$ Consequently, by a standard importance sampling argument it holds that Eq. 4.21 is a consistent estimate of $\mathbb{E}_{\gamma}[h(x_{1:T})]$ as $M$ becomes large, regardless of the number $N$ of particles used in each smc sampler.

It is not expected that Eq. 4.21 will be a better estimate than what we would obtain from a single smc sampler with a total of $NM$ particles. However, as pointed out above, the computation of the estimators Eq. 4.21 can be distributed, or even performed sequentially on the same machine in situations when memory constraints limit the number of particles than can be used.

4.4.2 The Island Particle Filter

In the preceding section we considered importance sampling of smc estimators. It is possible to take this one step further and develop smc sampling of smc estimators. This is the idea behind the island particle filter by Vergé et al. (2015).

Similarly to above, the island particle filter consists of $M$ nodes, each running an smc sampler with $N$ particles. However, these $M$ samplers are allowed to interact at each iteration according to the standard resampling, propagation, and weighting steps of smc. We provide pseudo-code in Algorithm 14, where we assume that resampling at the outer level only occurs when the effective sample size (ess)drops below some threshold (see Section 2.2.2). This type of adaptive resampling is recommended for a distributed implementation of the island particle filter. Indeed, resampling at the outer level at each iteration would result in very high communication costs, since it involves copying complete smc states (particles and weights) from node to node. Furthermore, we can control the outer level ess by increasing $N$, thereby ensuring that outer level resampling happens rarely when done adaptively.

The island particle filter fits into the standard smc framework. To make this explicit we start by writing $\mathbf{u}=\mathbf{u}_{1:T}$ where $\mathbf{u}_{t}$ corresponds to all the basic random variables generated at iteration $t$ of an (internal) smc sampler with $N$ particles. For example, $\mathbf{u}_{t}$ could correspond to a collection of independent uniform random variables needed to implement the resampling and propagation step at iteration $t$. Equivalently, we can view $\mathbf{u}_{t}$ as the random seed used for generating these variables. By making the (non restrictive) assumption that these auxiliary variables are drawn independently at each iteration we can write

[TABLE]

The island particle filter is then a standard smc sampler with $M$ particles, targeting the sequence of distributions

[TABLE]

and using $\tau_{t}(\mathbf{u}_{t})$ as a proposal distribution at iteration $t$. Similarly to the preceding section, targeting these “extended distributions” is motivated by the fact that $\mathbb{E}_{\Pi_{t}}[\widehat{h}_{t}({\mathbf{u}_{1:t}})]=\mathbb{E}_{\gamma}[h_{t}(x_{1:t})],$ where $\widehat{h}_{t}({\mathbf{u}_{1:t}})$ is defined analogously to (4.20).

If we denote the outer level (island level) unnormalized smc weights by $\{\widetilde{\Omega}_{t}^{m}\}_{m=1}^{M}$ we get

[TABLE]

at iterations when the islands are resampled, and

[TABLE]

at iterations when the outer level resampling is omitted; see Eq. 2.21. For implementation purposes we can use the latter expression at all iterations if we overwrite $\Omega_{t-1}^{m}\leftarrow\nicefrac{{1}}{{M}}$ whenever the islands are resampled. In the above expressions, $\{\widetilde{w}_{t}^{m,i}\}_{i=1}^{M}$ denote the unnormalized (inner level) smc weights at node $M$.

The derivation above is based on viewing the complete random seed as an auxiliary variable. In practice it is usually more convenient to simply keep track of the resulting particle trajectories and weights, which can be computed deterministically given $\mathbf{u}_{1:t}$. Thus, in Algorithm 14 we have identified $\mathbf{u}_{1:t}$ with $\mathbf{v}_{t}:=\left\{\left(x_{1:t}^{i},\widetilde{w}_{t}^{i}\right)\right\}_{i=1}^{N}$ and express the steps of the algorithm using the latter variable.

Appendix 4.A Proof of Unbiasedness

For clarity we prove the unbiasedness of the normalizing constant estimate for the version of smc presented in Algorithm 3, which involves multinomial resampling at each iteration of the algorithm. We note, however, that the property holds more generally, for instance when using low-variance and adaptive resampling. The proof presented here can be straightforwardly extended to such cases.

The distribution of all random variables generated by the smc method in Algorithm 3 is given by

[TABLE]

where the particles are $x_{1:T}^{1:N}=\bigcup_{t=1}^{T}\left\{x_{t}^{i}\right\}_{i=1}^{N}$, and the ancestor variables from the resampling step are $a_{1:T-1}^{1:N}=\bigcup_{t=1}^{T-1}\left\{a_{t}^{i}\right\}_{i=1}^{N}$. Assuming that $\widetilde{w}_{t}^{i}\geq 0$ for all $i$ and $t$, we have that $\widehat{Z}_{T}\geq 0$. We are now going to prove that

[TABLE]

from which the result follows. Note that the results holds for any value of $T$. To prove the result we introduce another set of variables $b_{T}^{i}=i$ and $b_{t}^{i}=a_{t}^{b_{t+1}^{i}}$, for $t=1,\ldots,T-1$. The notation describes the ancestor index at iteration $t$ for the final particle $x_{1:T}^{i}$. This means that we can write $x_{1:T}^{i}=x_{1:T}^{b_{1:T}^{i}}=(x_{1}^{b_{1}^{i}},\ldots,x_{T}^{b_{T}^{i}})$. Using this notation we rewrite the integrand in Eq. 4.23,

[TABLE]

where the first equality is the definition, the second equality separates the dependence of the particle $x_{1:T}^{i}=x_{1:T}^{b_{1:T}^{i}}$ from the rest, and the third equality simplifies the product between weights and proposals for particle $x_{1:T}^{i}=x_{1:T}^{b_{1:T}^{i}}$.

Now, using Eq. 4.23 and the above rewrite of the integrand we obtain

[TABLE]

where in the first equality we note that the sum over $i$ generates $N$ equal values, and thus we can arbitrarily choose one of these (in this case $i=1$) to evaluate and multiply the result by $N$. In the second equality we marginalize the variables not involved in the particle $x_{1:T}^{b_{1:T}^{1}}$. The two final equalities follow because we are averaging $N^{T-1}$ values that are all equal to $1$. This concludes the proof.

Chapter 5 Conditional SMC: Algorithms and Applications

Conditional sequential Monte Carlo (csmc)is a recent algorithm combining smc and mcmc (Andrieu et al., 2010), originally developed as a means to approximate and mimic idealized Gibbs samplers for local latent variable models. The intractability of the exact Gibbs sampler stems from the challenge in simulating from the exact conditional of latent variables given the data and the parameters. The key idea in csmc is to simulate from this conditional distribution approximately using a slightly modified smc approximation, retaining a single particle at each iteration. By iterating this procedure we obtain a mcmc method that has the desired target distribution as its stationary distribution. We will see how this can be used not only for simulating from a target distribution, but also how we can use it to generate unbiased inverse normalization constant estimates, speed up inference for e.g. probabilistic programming, and evaluating other approximate inference methods.

In the first section below we introduce the csmc algorithm and illustrate it with an example. Then, we discuss the unbiasedness property of the csmc inverse normalization constant estimate, a key property underpinning the rest of the section. Using csmc as a tool we move on by deriving the distribution of the expected smc empirical distribution approximation. Next, we introduce the ipmcmc algorithm, leveraging csmc and smc to derive a method that can take advantage of distributed and multi-core computational architectures. Finally, we discuss ways in which smc and csmc can be used for evaluating approximate inference.

5.1 Conditional Sequential Monte Carlo

One way to look at csmc is to view it as Markov kernel with the target distribution $\gamma_{T}$ as its stationary distribution. Each iteration takes the previous sample as an input, a reference trajectory, and outputs an updated sample with no need for an accept–reject step. If we keep iterating, applying csmc over and over, we obtain the iterated csmc algorithm. The iterated csmc sampler is a type of Gibbs sampler on an extended space of all random variables generated at each step (Andrieu et al., 2010). We will sometimes refer to the previous sample $x_{1:T}^{\prime}$, i.e. the reference trajectory, as the retained particle.

The csmc algorithm differs from standard smc described in Section 2.2 only in that a single particle is set deterministically to the retained particle from the previous iteration. This means that we condition on the event that a specific particle survives and is equal to the retained particle $x_{1:T}^{\prime}$, hence the name conditional smc. We obtain the csmc algorithm by a simple modification of the smc proposal and weight update equations. At each iteration, for the first $N-1$ particles, $i=1,\ldots,N-1$, we follow the standard updates according to Eqs. 2.15 and 2.16. The proposal is

[TABLE]

and the weights

[TABLE]

The last particle $i=N$ we set deterministically equal to the retained particle $x_{1:t}^{N}=x_{1:t}^{\prime}$. Its weight is computed according to Eq. 5.2, that is we use the same weight expression for all $N$ particles. The empirical distribution approximation of $\gamma_{T}$ is, analogously to smc, equal to

[TABLE]

The updated retained particle is obtained by drawing a single sample from the empirical distribution approximation in Eq. 5.3, i.e. $x_{1:T}^{\prime}\sim\widehat{\gamma}_{T}(x_{1:T})$.

We summarize the csmc algorithm in Algorithm 15. Iterated csmc is obtained by letting the retained particle at the next iteration be the output of the previous step. This algorithm constructs a Markov chain, consisting of the retained particles at each iteration, with the target distribution $\gamma_{T}$ as its stationary distribution. For a thorough treatment of the theory of iterated csmc, particle Gibbs (pg) and particle Markov chain Monte Carlo (pmcmc), i.e. the extension to sampling also model parameters, see Andrieu et al. (2010).

The path degeneracy issue discussed in Section 2.2 is particularly prominent in iterated csmc. Because we enforce that the retained particle survives to the very end, we simultaneously force that all the particles collapse for early steps $t$ to $x_{1:t}^{\prime}$. This means that the output from Algorithm 15, with high probability, is equal to the input for early values of $x_{1:T}^{\prime}$. The samples we get are highly correlated and any estimator will require a significant number of samples to reasonably approximate expectation with respect to the target distribution $\gamma_{T}$. By modifying the way we update particle $N$, the retained particle, we can alleviate this problem. One way to do this is using ancestor sampling (Lindsten et al., 2014). With ancestor sampling, instead of setting $x_{1:t}^{N}=x_{1:t}^{\prime}$ deterministically as in Algorithm 15, we sample at each time step

[TABLE]

where the auxiliary weights $\widetilde{\nu}_{t}^{i}$ are given by

[TABLE]

Using ancestor sampling does not avoid path degeneracy but it can allow for the particle trajectory we degenerate to to be different than the retained particle from a previous iteration. For an in-depth treatment of the path degeneracy issue and its other potential solutions we refer to Lindsten and Schön (2013).

Example 5.1.1** (Conditional sequential Monte Carlo for Example 1.2.1).**

We revisit our running non-Markovian Gaussian example. We let $T=20$ and study the performance of iterated csmc without (blue) and with ancestor sampling (red) in Fig. 5.1. We set $N=5$, $10$ and study the kernel density estimates based on $10,000$ mcmc samples of the marginal distributions $\gamma_{T}(x_{1})$ and $\gamma_{T}(x_{T})$, respectively.

We can see the striking effect that path degeneracy has the iterated csmc approximation in Fig. 5.1(a); effectively we approximate $\gamma_{T}(x_{1})$ using a single sample when $N=5$, whereas the approximation for $\gamma_{T}(x_{T})$ is more diverse. When using ancestor sampling on the other hand, the approximation is already accurate for $N=5$. This accuracy comes at a cost however, the ancestor sampling step requires that we evaluate Eq. 5.5 at each iteration. This extra step has a computational cost that in general may be proportional to $T^{2}$, whereas the cost of standard csmc is often only on the order of $T$. For some target distributions, such as the ssm or our running example, there is only a small constant computational overhead associated with evaluating Eq. 5.5. However, whether the benefit of ancestor sampling will outweigh its computational cost will have to be evaluated on a case by case basis.

5.2 The Unbiasedness of the Inverse Normalization Constant Estimate

Unlike in the standard smc algorithm, the normalization constant estimate for csmc based on the product of the average weights is biased. We restate the normalization constant estimate below:

[TABLE]

Because we condition on a reference trajectory, this normalization constant estimate has different properties than does the smc estimate. However, it is still consistent, converging almost surely to $Z_{T}$ as the number of particles tends to infinity.

While the above estimator is biased for a finite number of particles, interestingly the csmc inverse normalization constant estimate $\nicefrac{{1}}{{\widehat{Z}_{T}}}$ is unbiased. We will see that this is particularly relevant for the inference evaluation methods that we will study in Section 5.5. The unbiasedness holds under the assumption that the reference trajectory is distributed according to the target distribution $x_{1:T}^{\prime}\sim\gamma_{T}(x_{1:T})$. This assumption relies on the Markov chain having already converged to its stationary distribution. However, in some applications it is possible to make use of csmc with the problem set up in such a way that this is always true.

We formalize and prove the unbiasedness result in Proposition 5.2.1.

Proposition 5.2.1**.**

Assuming that the retained particle $x_{1:T}^{\prime}\sim\gamma_{T}(x_{1:T})$, then the csmc inverse normalization constant estimate $\frac{1}{\widehat{Z}_{T}}$ is non-negative and unbiased

[TABLE]

Proof.

See Appendix 5.A. ∎

5.3 The Expected SMC Empirical Distribution

Another method for representing the expected smc empirical distribution $\displaystyle\mathbb{E}\left[\widehat{\gamma}_{T}(x_{1:T})\right]$ is using csmc. This can be useful for some applications of smc, such as in variational smc or in inference evaluation which we will see later in this chapter. The marginal distribution of a single sample from the smc empirical distribution $\widehat{\gamma}_{T}$ is

[TABLE]

where $\textrm{Q}_{\text{SMC}}$ (see Eq. 4.22 for its definition) is the distribution of all particles and ancestor variables generated by a single run of the smc algorithm.

We formalize the alternate representation of the expected smc empirical distribution in Proposition 5.3.1.

Proposition 5.3.1**.**

The expected smc empirical distribution satisfies

[TABLE]

where $\textrm{Q}_{\text{CSMC}}$ is the distribution of all random variables generated by a single run of the csmc algorithm with $x_{1:T}$ as its retained particle.

Proof.

See Appendix 5.B. ∎

5.4 Parallelization

In Section 4.4 we briefly discussed a few approaches for leveraging multi-core and distributed computational architectures. Here we will see how we can combine smc and csmc to derive ipmcmc (Rainforth et al., 2016), an algorithm that takes advantage of parallelization to speed up inference for amongst other things probabilistic programming.

The ipmcmc method, like the iterated csmc, can be thought of as a type of Gibbs sampler on an extended space that draws approximate samples from the target distribution $\gamma_{T}$. However unlike iterated csmc, the ipmcmc algorithm generates not only one sample at each iteration but several. The method is based on a pool of several interacting standard and conditional sequential Monte Carlo samplers. These pools of samplers are run as parallel processes and will be referred to as nodes below. Each node either runs an iteration of csmc like in Algorithm 15, or an iteration of smc like in Algorithm 3. After each run of this pool, we update the node indices running csmc by applying successive Gibbs updates. This means that the nodes running csmc can change from iteration to iteration. This approach can let the user trade off exploration (smc) and exploitation (csmc) to achieve a faster and more accurate approximate inference method.

To simplify notation we will let $\mathbf{x}$ denote the full trajectory $x_{1:T}$. We assume that we have a total of $M$ nodes, indexed by $j=1,\ldots,M$, of which $P$ run csmc and the remaining run standard smc. The conditional nodes, i.e. the nodes running csmc at any given time, are denoted by a set of indices $c_{1:P}\in\{1,\ldots,M\}^{P}$ where each $c_{j}$ is distinct. The conditional node $c_{j}$’s retained particle is referred to by $\mathbf{x}_{j}^{\prime}$. After each run all nodes produce a normalization constant estimate $\widehat{Z}_{T}^{j}$, i.e. the product of the average unnormalized weights in that node. Ipmcmc proceeds by successively updating the conditional node indices and retained particles by Gibbs steps. For $j=1,\ldots,P$ we sample

[TABLE]

where $\mathbbm{1}\left(\cdot\right)$ is the identity function evaluating to one if its argument is true and zero otherwise. Furthermore, we have used the node index shorthand $c_{\neg j}:=\{c_{1},\ldots,c_{j-1},c_{j+1},\ldots,c_{P}\}$. After selecting new nodes to become conditional nodes through Eq. 5.10 we select their new corresponding retained particle by simulating

[TABLE]

where $\widehat{\gamma}_{T}^{c_{j}}$ is the standard target distribution approximation from node $c_{j}$.

We summarize one iteration of ipmcmc in Algorithm 16. Iterating the update steps, letting the retained particles at the previous step be the input for the subsequent step, generates a Markov chain on an extended space where each retained particle is approximately distributed according to the target distribution.

Ipmcmc is suited to distributed and multi-core architectures. The main computational cost is running the smc and csmc iterations, which can be done in parallel by the nodes. At each iteration, only a single normalization constant estimate needs be communicated between nodes and the calculation of the updates of $c_{1:P}$ is negligible. Retained particles do not need to be communicated and may be stored locally in each node until needed. Furthermore, ipmcmc may be amenable to asynchronous adaptation under certain assumptions (Rainforth et al., 2016). The algorithm is also suited for use in probabilistic programming where we are often only able to simulate from the prior and evaluate the likelihood. Interacting particle Markov chain Monte Carlo can allow for tackling more challenging inference problems, i.e. more complex probabilistic programs, compared to standard pmcmc methods due to leveraging distributed and multi-core computational architectures (Rainforth et al., 2016).

5.5 Evaluating Inference

Approximate Bayesian inference is central to a large part of probabilistic machine learning and its application. It is a coherent approach to treat uncertainty in both the data and the model with applications to everything from robotics to health care. A notoriously difficult problem is to measure the accuracy of our posterior approximation compared to the true posterior distribution. Inference checking, or inference evaluation, refers to the process of comparing our approximate inference result to what we would have obtained using the true posterior. In practice whether an approximation is accurate or not may differ between different applications, we might be interested in different aspects of the posterior. Some examples are tail probabilities in rare event estimation or posterior mean for minimum mean square error estimation. In this section we will focus on comparing directly the posterior approximation to the exact posterior distribution in terms of symmetrized kl divergence, also known as the Jeffreys divergence, following Grosse et al. (2016); Cusumano-Towner and Mansinghka (2017).

The posterior distribution that we want to approximate is denoted by $p(\mathbf{x}\,|\,\mathbf{y})$, and the approximate posterior will be denoted by $q(\mathbf{x})$. The approximation $q$ could for example be the expected smc empirical distribution approximation, similar to vsmc, or it could be an optimized variational approximation, or any of the other examples detailed in Grosse et al. (2016); Cusumano-Towner and Mansinghka (2017). In either case, the symmetrized kl divergence from $q$ to $p$ is given by

[TABLE]

where the standard Kullback-Leibler divergence from $q(\mathbf{x})$ to $p(\mathbf{x})$ is defined by $\textrm{KL}\left(q(\mathbf{x})\,\|\,p(\mathbf{x})\right):=\mathbb{E}_{q(\mathbf{x})}\left[\log q(\mathbf{x})-\log p(\mathbf{x})\right]$. We can rewrite Eq. 5.12 to eliminate the need for computing the marginal likelihood $p(\mathbf{y})$ as follows

[TABLE]

To be able to estimate Eq. 5.13 using a standard mc estimator we require that we can sample $q(\mathbf{x})$ and $p(\mathbf{x}\,|\,\mathbf{y})$, respectively. Furthermore, we require that we are able to evaluate the pdfs $q(\mathbf{x})$ and $p(\mathbf{x},\mathbf{y})$ pointwise. These requirements can be restrictive for the types of problems we encounter in practice, and Grosse et al. (2016); Cusumano-Towner and Mansinghka (2017) make various relaxations and assumptions to make the problem tractable. Below we will detail the two approaches to inference checking known as bounding divergences with reverse annealing (bread) (Grosse et al., 2016) and auxiliary inference divergence estimator (aide) (Cusumano-Towner and Mansinghka, 2017).

Bounding Divergences with Reverse Annealing

The bread approach focuses on bounding the symmetrized kl for a posterior approximation based on annealed importance sampling (ais) (Neal, 2001) for simulated data. By using simulated data we ensure that we can generate an exact sample from the posterior distribution $\mathbf{x}\sim p(\mathbf{x}\,|\,\mathbf{y})$. We can think of ais as a single sample, $N=1$, smc method that uses as its proposal a Markov kernel with stationary distribution $\gamma_{t}(x_{t})$. The target distribution is typically chosen according to a data or likelihood tempering scheme as discussed in Section 1.2. Per definition the posterior distribution is obtained as the marginal at the final iteration, i.e. $p(\mathbf{x}\,|\,\mathbf{y})=\gamma_{T}(x_{T})$ where $\mathbf{x}\equiv x_{T}$. The posterior distribution approximation is then defined by $q(\mathbf{x}):=\mathbb{E}[\widehat{\gamma}_{T}(x_{T})]$ and the expectation is with respect to all the proposal distributions.

Because of this set up the intractability in estimating, or bounding, the Jeffreys divergence in Eq. 5.13 comes from having to evaluate the pdf $q(\mathbf{x})$ pointwise. We can bound Eq. 5.13 from above by noting that marginalizing out variables cannot increase the kl divergence

[TABLE]

where the equality follows from Proposition 5.3.1. We can further bound Eq. 5.14 using Jensen’s inequality and simplifying, leading to

[TABLE]

The inequality in Eq. 5.15 holds generally even if we would use a standard smc algorithm rather than ais as our posterior approximation. The details are provided in Appendix 5.C. However, when ais is used we can further simplify to obtain the same bound as Grosse et al. (2016)

[TABLE]

where $q(x_{1:T})=q_{1}(x_{1})\prod_{t=1}^{T}q_{t}(x_{t}\,|\,x_{t-1})$ are the proposals (forward kernels), $s(x_{1:T-1}\,|\,x_{T})=\prod_{t=1}^{T-1}s_{t}(x_{t}\,|\,x_{t+1})$ are the corresponding backward kernels (cf. Section 1.2, Conditionally independent models), and the weights are in the likelihood tempering variant given by

[TABLE]

with $\widetilde{w}_{1}=1$. This means it is possible to obtain a stochastic upper bound, i.e. a mc estimate of the upper bound in Eq. 5.16, by running a standard (forward) ais evaluating the negative log-weights, and running a reverse ais evaluating the log-weights; the left and right expressions in Eq. 5.16, respectively.

Auxiliary Inference Divergence Estimator

The aide approach, on the other hand, treats a more general class of posterior approximations such as variational inference, smc, or mcmc. However, rather than assuming simulated data, Cusumano-Towner and Mansinghka (2017) propose to study the Jeffreys divergence from the approximation to that of a gold standard inference approximation. Typically this gold standard approximation would consist of a long run of an mcmc method or smc with a large number of particles that we know approximates the posterior well.

In this setting not only is the pdf for the approximation $q(\mathbf{x})$ difficult to evaluate, so is the pdf for our gold standard algorithm. Rather than introducing new notation for the gold standard algorithm we will keep the posterior distribution notation $p(\mathbf{x}\,|\,\mathbf{y})$. The approach that the aide takes is to use meta-inference methods, based on smc or csmc, to obtain approximations of the pdfs and their reciprocals. A meta-inference algorithm takes as input a sample from the distribution of interest and outputs an unbiased estimate of either its pdf or the reciprocal. In fact it suffices that the meta-inference algorithms provide estimates that are unbiased up to a constant of proportionality because any constant is canceled by the symmetry in the Jeffrey’s divergence.

We note that the approach suggested by Cusumano-Towner and Mansinghka (2017) is essentially equivalent to the bound derived in Eqs. 5.14, 5.15 and 5.16. However, instead of the exact target distribution $\gamma_{T}(x_{1:T})$ we have another smc approximation to represent the gold standard algorithm. To differentiate between the two smc approximations we will use superscript $q$ for the approximation and $p$ for the gold standard. With this the corresponding Jeffrey’s divergence is

[TABLE]

where $\textrm{Q}_{\text{CSMC}}^{q}$, $\widehat{Z}_{T}^{q}$ and $\textrm{Q}_{\text{CSMC}}^{p}$, $\widehat{Z}_{T}^{p}$ are the csmc distributions and normalization constant estimates for the approximation and gold standard algorithm, respectively. The first equality follows from the definition of $q(\mathbf{x})$ and $p(\mathbf{x}\,|\,\mathbf{y})$. The second equality follows from Proposition 5.3.1 for the two distributions, respectively. The divergence in Eq. 5.17 can be bounded, analogously to the bread case, as follows

[TABLE]

where $q(\mathbf{x})$ and $p(\mathbf{x}\,|\,\mathbf{y})$ are the two corresponding expected smc empirical distributions.

The bound above differs slightly from that derived by Cusumano-Towner and Mansinghka (2017). The aide procedure uses a form of importance weighted version of Eq. 5.18, in their notation it is equivalent to the bound based on a single run of the meta-inference for the approximation ($M_{t}=1$) and gold standard ($M_{g}=1$). For more details we refer the interested reader to that paper.

Appendix 5.A Proof of Unbiasedness

When introducing the csmc method we consistently set the $N$:th particle to the reference trajectory. In practice it does not matter which particle is deterministically set to the reference trajectory. It will prove simplifying to let the particle that is conditioned on, i.e. the reference trajectory, have its own index $b_{1:T}=(b_{1},\ldots,b_{T})\in\{1,\ldots,N\}^{T}$. For clarity we focus on using multinomial resampling at each iteration of the algorithm. The distribution of all random variables generated by the csmc method in Algorithm 15 is then given by

[TABLE]

where the particles are $x_{1:T}^{\neg b_{1:T}}=\bigcup_{t=1}^{T}\left\{x_{t}^{i}\right\}_{i=1}^{N}\backslash\{x_{t}^{b_{t}}\}_{t=1}^{T}$, and the ancestor variables from the resampling step are $a_{1:T-1}^{\neg b_{1:T-1}}=\bigcup_{t=1}^{T-1}\left\{a_{t}^{i}\right\}_{i=1}^{N}\backslash\{b_{t}\}_{t=1}^{T-1}$. Assuming that $\widetilde{w}_{t}^{i}\geq 0$ for all $i$ and $t$, we have that $\widehat{Z}_{T}\geq 0$ and thus $\nicefrac{{1}}{{\widehat{Z}_{T}}}\geq 0$. We are now going to prove that when $\left(x_{1:T}^{b_{1:T}},b_{1:T}\right)\sim\gamma_{T}(x_{1:T})\cdot N^{-T}$

[TABLE]

from which the result follows. We rewrite the integrand in Eq. 5.20,

[TABLE]

where the first line is the integrand. In the first equality we have divided and multiplied by $q_{1}\prod_{t=2}^{T}w_{t-1}^{b_{t-1}}q_{t}$ to be able to identify the difference between $\textrm{Q}_{\text{SMC}}$ and $\textrm{Q}_{\text{CSMC}}$. Note that we have replaced $b_{t}=a_{t}^{b_{t+1}}$ in the notation for $\textrm{Q}_{\text{SMC}}\left(x_{1:T}^{1:N},a_{1:T-1}^{1:N}\right)$ to identify the smc distribution. For the second equality we have rewritten the normalized weights in the denominator, together with $N^{-T}$ and dividing and multiplying by $\sum_{i}\widetilde{w}_{T}^{i}$, to identify the normalization constant estimator. For the final equality we see that the normalization constant estimators cancel, and the product in the denominator simplifies together with $\widetilde{\gamma}_{T}$ to become $\widetilde{w}_{T}^{b_{T}}$.

Using this simplified expression for the integrand in Eq. 5.20 we obtain

[TABLE]

where again we first replace the original sum over $b_{1:T-1}$ with the ancestor variable $a$ notation. The first equality marginalizes $b_{T}$ and the final equality follows because $\textrm{Q}_{\text{SMC}}$ is a distribution on $x_{1:T}^{1:N},a_{1:T-1}^{1:N}$. This concludes the proof.

Appendix 5.B Proof of SMC Distribution Representation

We will prove the following alternate representation of the expected smc empirical distribution

[TABLE]

were we for clarity again focus on using multinomial resampling at each iteration of the algorithm. Note that the left hand side is a probability measure and the right hand side a probability distribution function. To be precise the equality means that the corresponding probability measure of the right hand side (usually adding either the Lebesgue our counting measure) is equal almost everywhere to the left hand side.

To prove the result we need to show that the expectation of any bounded and continuous function $h(x_{1:T})$ coincides for the left and right hand sides of Eq. 5.21. We have that the expectation of $h(x_{1:T})$ with respect to the right hand side of Eq. 5.21 is

[TABLE]

where the first equality is the definition. In the second equality we have expanded the terms and can identify the expression derived in Appendix 5.A above. In the third equality we have replaced it with the corresponding $\textrm{Q}_{\text{SMC}}$ distribution, again using $b_{t}=a_{t}^{b_{t+1}}$ in the notation for $\textrm{Q}_{\text{SMC}}$. In the penultimate equality we re-arrange a few terms and replace the sum over $b_{T}$ with an integration with respect to a point mass distribution. In the final equality we identify the expectation with respect to $\mathbb{E}_{\textrm{Q}_{\text{SMC}}}\left[\widehat{\gamma}_{T}(x_{1:T})\right]$. This concludes the proof.

Appendix 5.C Jeffreys Divergence Bound

We detail the derivation of the bound in Eq. 5.15 below under the assumption of multinomial resampling at each iteration of the algorithm,

[TABLE]

where the first equality is the definition. The second equality is a simplification of the previous expression where the terms cancelling, $\widetilde{\gamma}_{T},\gamma_{T},Z_{T}$, have been removed. Next, the inequality follows from the (conditional) Jensen’s inequality for the convex functions $x\log x$ and $-\log x$, respectively. The final equality follows from Proposition 5.3.1, noting that the first term is equivalent to an expectation with respect to the full smc distribution $\textrm{Q}_{\text{SMC}}$.

The Jeffrey’s divergence bound between two generic expected smc empirical distributions in Eq. 5.17 is a straightforward extension to the above derivation.

Chapter 6 Discussion

We have described sequential Monte Carlo, methods that use random samples to approximate the posterior. The goal is to approximate the distribution of the latent (unobserved) variables given the data. The idea is to draw samples from a (simpler) proposal distribution, then correct for the discrepancy between the proposal and the posterior using weights. smc is generally applicable and can be considered to be an alternative, or complement, to mcmc methods.

First, we discussed the roots of smc in (sequential) importance sampling and sampling/importance resampling. Next, we introduced the smc algorithm in its full generality. The algorithm was illustrated with synthetic data from a non-Markovian latent variable model. Then, we discussed practical issues and theoretical results. The smc estimate satisfies favorable theoretical properties; under appropriate assumptions it is both consistent and has a central limit theorem. Then, we discussed and explained how to learn a good proposal and target distribution. As a motivating example, we described how to go about learning a proposal- and target distribution for a stochastic recurrent neural network. Next, we discussed proper weighting and how to use smc to construct nested estimators based on several layers of mc approximations to approximate idealized mc algorithms. We provided a straightforward proof to the unbiasedness property of the smc normalization constant estimate. We studied pmh algorithms, methods that construct high-dimensional proposals for mh based on smc samplers. We studied random and nested smc methods that deal with intractability in the weights and proposal distributions with applications to stochastic differential equations and optimal proposal approximation. We illustrated how to leverage proper weighting and smc to design algorithms that makes use of distributed and multi-core computing architectures. Finally, we introduced conditional sequential Monte Carlo algorithms and applications of these. We provided proof that the csmc sampler obtains unbiased estimates of the inverse normalization constant. We showed how the ipmcmc method combines smc and csmc to design a powerful mcmc method that can take advantage of parallel computing. As a motivating example we show how smc and csmc can be used to evaluate inference in the context of approximate Bayesian inference.

There are several avenues open for further work and research in smc methods and their application. We provide a few examples below. The proposal distribution is crucial for performance, and finding efficient ways to learn useful proposals is an area of continued research. The twisting potential can drastically improve accuracy. However, jointly learning twisting potentials and proposal distributions for complex models is a challenging and largely unsolved problem.

A rather obvious—but still very interesting and important—direction consists in adapting the smc algorithm to parallel and distributed computing environments. This is indeed an active area of research, see e.g. Lee and Whiteley (2016); Lee et al. (2010). However, we have reason to believe that there is still a lot that can be done in this direction. Just to give a few concrete examples we mention the fact that we can assemble smc algorithm in new ways that we directly designed for parallel and distributed computing. As a concrete example we mention the smoother developed by Jacob et al. (2019) which is based on the clever debiasing scheme of Glynn and Rhee (2014). Perhaps the divide and conquer SMC algorithm (Lindsten et al., 2017) can be useful in this effort. Its successful implementation within a probabilistic programming environment with support for parallel computation might give rise to an interesting scalable use of smc.

An important observation is that smc is increasingly being used as a component in various larger (often iterative) algorithms. Just to give a few examples of this we have Bayesian learning, where the particle filter is used within a Markov chain Monte Carlo algorithm (Andrieu et al., 2010; Lindsten et al., 2014), iterated filtering (Ionides et al., 2015) and maximum likelihood using expectation maximization (Schön et al., 2011). The idea of using smc within itself has occurred in several forms, including SMC2 by Chopin et al. (2013) and nested smc by Naesseth et al. (2015). There are still most likely many interesting algorithms to uncover by following this idea in new directions.

There has recently been a very interesting line of research when it comes to merging deep learning architectures—such as recurrent neural networks, temporal convolutional networks and gated recurrent units—with stochastic latent models—such as the state space model—for sequence modelling, see e.g. (Bai et al., 2018; Aksan and Hilliges, 2019; Osendorfer and Bayer, 2014; Fraccaro et al., 2016). An important idea here is to combine the expressiveness of the deep neural networks with the stochastic latent variable model’s ability to represent uncertainty. These models are often trained by optimizing the ELBO objective. It has recently been shown that we can exploit the fact that the smc estimator of the likelihood is unbiased to construct tighter lower bounds of this objective (Le et al., 2018; Maddison et al., 2017; Naesseth et al., 2018). This can most likely open up for better training of combinations of deep learning architectures and stochastic latent models.

Given that we have now understood that the smc methods are much more generally applicable than most of us first thought we foresee an interesting and creative future for these algorithms. We hope that this paper can serve as a catalyst for further research on these methods and in particular on their application in machine learning.

Acknowledgments

This research was financially supported by the Swedish Foundation for Strategic Research (SSF) via the projects ASSEMBLE (contract number: RIT15-0012) and Probabilistic Modeling and Inference for Machine Learning (contract number: ICA16-0015), and by the Swedish Research Council via the projects NewLEADS - New Directions in Learning Dynamical Systems (contract number: 621-2016-06079) and Learning of Large-Scale Probabilistic Dynamical Models (contract number: 2016-04278).