Variational Sequential Monte Carlo

Christian A. Naesseth, Scott W. Linderman, Rajesh Ranganath, David M. Blei

Introduction

Complex data like natural images, text, and medical records require sophisticated models and algorithms. Recent advances in these challenging domains have relied upon variational inference (vi) (Kingma and Welling, 2014; Hoffman et al., 2013; Ranganath et al., 2016a). Variational inference excels in quickly approximating the model posterior, yet these approximations are only useful insofar as they are accurate. The challenge is to balance faithful posterior approximation and fast optimization.

We present a new approximating family of distributions called variational sequential Monte Carlo (vsmc). vsmc blends vi and sequential Monte Carlo (smc) (Stewart and McCarty, 1992; Gordon et al., 1993; Kitagawa, 1996), providing practitioners with a flexible, accurate, and powerful approximate Bayesian inference algorithm. vsmc is an efficient algorithm that can approximate the posterior arbitrarily well.

Standard smc approximates a posterior distribution of latent variables with NN weighted particles iteratively drawn from a proposal distribution. The idea behind variational smc is to view the parameters of the proposal as indexing a family of distributions over latent variables. Each distribution in this variational family corresponds to a particular choice of proposal; to sample the distribution, we run smc to generate a set of particles and then randomly select one with probability proportional to its weight. Unlike typical variational families, the vsmc family trades off fidelity to the posterior with computational complexity: its accuracy increases with the number of particles NN, but so does its computational cost.

We develop the vsmc approximating family, derive its corresponding variational lower bound, and design a stochastic gradient ascent algorithm to optimize its parameters. We connect vsmc to the importance weighted auto-encoder (iwae) (Burda et al., 2016) and show that the iwae lower bound is a special case of the vsmc bound. As an illustration, consider approximating the following posterior with latent variables x1:Tx_{1:T} and observations y1:Ty_{1:T},

This is a toy Gaussian state space model (ssm) where the observed value at each time step depends on the square of the latent state. Figure 1c shows the approximating power of vsmc versus that of the iwae and of standard variational Bayes (vb). As the length of the sequence TT increases, naïve importance sampling effectively collapses to use only a single particle. vsmc on the other hand maintains a diverse set of particles and thereby achieves a significantly tighter lower bound of the log-marginal likelihood log⁡p(y1:T)\log p(y_{1:T}).

We focus on inference in state space and time series models, but emphasize that vsmc applies to any sequence of probabilistic models, just like standard smc (Del Moral et al., 2006; Doucet and Johansen, 2009; Naesseth et al., 2014).

In Section 5, we demonstrate the advantages of vsmc on both simulated and real data. First, we show on simulated linear Gaussian ssm data that vsmc can outperform the (locally) optimal proposal (Doucet et al., 2001; Doucet and Johansen, 2009). Then we compare vsmc with iwae for a stochastic volatility model on exchange rates from financial markets. We find that vsmc achieves better posterior inferences and learns more efficient proposals. Finally, we study recordings of macaque monkey neurons using a probabilistic model based on recurrent neural networks. vsmc reaches the same accuracy as iwae, but does so with less computation.

Much effort has been dedicated to learning good proposals for smc (Cornebise, 2009). Guarniero et al. (2017) adapt proposals through iterative refinement. Naesseth et al. (2015) uses a Monte Carlo approximation to the (locally) optimal proposal (Doucet and Johansen, 2009). Gu et al. (2015) learn proposals by minimizing the Kullback-Leibler (kl) from the posterior to proposal using SMC samples; this strategy can suffer from high variance when the initial SMC proposal is poor. Paige and Wood (2016) learn proposals by forward simulating and inverting the model. In contrast to all these methods, vsmc optimizes the proposal directly with respect to KL divergence from the smc sampling process to the posterior.

vsmc uses auxillary variables in a posterior approximation. This relates to work in vi, such as Hamiltonian VI (Salimans et al., 2015), variational Gaussian processes (Tran et al., 2016), hierarchical variational models (Ranganath et al., 2016b), and deep auxiliary variational auto-encoders (Maaløe et al., 2016). Another approach uses a sequence of invertible functions to transform a simple variational approximation to a complex one (Rezende and Mohamed, 2015; Dinh et al., 2014). All of these rich approximations can be embedded inside vsmc to build more flexible proposals.

Archer et al. (2015); Johnson et al. (2016) develop variational inference for state space models with conjugate dynamics, while Krishnan et al. (2017) develop variational approximations for models with nonlinear dynamics and additive Gaussian noise. In contrast, vsmc is agnostic to the distributional choices in the dynamics and noise.

Importance weighted auto-encoders (Burda et al., 2016) obtain the same lower bound as variational importance sampling (vis), a special case of vsmc. However, vis provides a new interpretation that enables a more accurate variational approximation; this relates to another interpretation of iwae by Cremer et al. (2017); Bachman and Precup (2015). Variational particle approximations (Saeedi et al., 2014) also provide variational approximation that improve with the number of particles, but they are restricted to discrete latent variables.

Finally, the log-marginal likelihood lower bound (6) was developed concurrently and independently by Maddison et al. (2017) and Le et al. (2017). The difference with our work lies in how we derive the bound and the implications we explore. Maddison et al. (2017); Le et al. (2017) derive the bound using Jensen’s inequality on the smc expected log-marginal likelihood estimate, focusing on approximate marginal likelihood estimation of model parameters. Rather, we derive (6) as a tractable lower bound to the exact evidence lower bound (elbo) for the new variational family vsmc. In addition to a lower bound on the log-marginal likelihood, this view provides a new variational approximation to the posterior.

Background

We begin by introducing the foundation for variational sequential Monte Carlo (vsmc). Let p(x1:t,y1:t)p(x_{1:t},y_{1:t}) be a sequence of probabilistic models for latent (unobserved) x1:tx_{1:t} and data y1:ty_{1:t}, with t=1,…,Tt=1,\ldots,T. In Bayesian inference, we are interested in computing the posterior distribution p(x1:T ∣ y1:T)p(x_{1:T}\,|\,y_{1:T}). Two concrete examples, both from the time-series literature, are hidden Markov models and state space models (Cappé et al., 2005). In both cases, the joint density factorizes as

where ff is the prior on xx, and gg is the observation (data) distribution. For most models computing the posterior p(x1:T ∣ y1:T)p(x_{1:T}\,|\,y_{1:T}) is computationally intractable, and we need approximations such as vi and smc. Here we construct posterior approximations that combine these two ideas.

In the following sections, we review variational inference and sequential Monte Carlo, develop a variational approximation based on the samples generated by smc, and develop a tractable objective to improve the quality of the smc variational approximation. For concreteness, we focus on the state space model above. But we emphasize that vsmc applies to any sequence of probabilistic models, just like standard smc (Del Moral et al., 2006; Doucet and Johansen, 2009; Naesseth et al., 2014).

In variational inference we postulate an approximating family of distributions with variational parameters λ\lambda, q(x1:T;λ)q(x_{1:T};\lambda). Then we minimize a divergence, often the kl divergence, between the approximating family and the posterior so that q(x1:T;λ)≈p(x1:T ∣ y1:T){q(x_{1:T};\lambda)\approx p(x_{1:T}\,|\,y_{1:T})}. This minimization is equivalent to maximizing the elbo (Jordan et al., 1999),

vi turns posterior inference into an optimization problem.

Sequential Monte Carlo

smc is a sampling method designed to approximate a sequence of distributions, p(x1:t ∣ y1:t)p(x_{1:t}\,|\,y_{1:t}) for t=1…Tt=1\ldots T with special emphasis on the posterior p(x1:T ∣ y1:T)p(x_{1:T}\,|\,y_{1:T}). For a thorough introduction to smc see Doucet and Johansen (2009); Doucet et al. (2001); Schön et al. (2015).

To approximate p(x1:t ∣ y1:t)p(x_{1:t}\,|\,y_{1:t}) smc uses weighted samples,

where δX\delta_{X} is the Dirac measure at XX.

We construct the weighted set of particles sequentially for t=1,…,Tt=1,\ldots,T. At time t=1t=1 we use standard importance sampling x1i∼r(x1)x_{1}^{i}\sim r(x_{1}). For t>1t>1, we start each step by resampling auxiliary ancestor variables at−1i∈{1,…,N}{a_{t-1}^{i}\in\{1,\ldots,N\}} with probability proportional to the importance weights wt−1jw_{t-1}^{j}; next we propose new values, append them to the end of the trajectory, and reweight as follows:

We refer to the final particles (samples) x1:Tix_{1:T}^{i} as trajectories. Panels (a) and (b) of Figure 1 show sets of weighted trajectories. The size of the dots represents the weights wtiw_{t}^{i} and the arrows represent the ancestors at−1ia_{t-1}^{i}. Importance sampling omits the resampling step, so each ancestor is given by the corresponding particle for the preceding time step.

The trajectories x1:Tix_{1:T}^{i} and weights wTiw_{T}^{i} define the smc approximation to the posterior. Critically, as we increase the number of particles, the posterior approximation becomes arbitrarily accurate. smc also yields an unbiased estimate of the marginal likelihood,

This estimate will play an important role in the vsmc objective.

The proposal distribution r(xt ∣ xt−1){r(x_{t}\,|\,x_{t-1})} is the key design choice. A common choice is the model prior ff—it is known as the bootstrap particle filter (bpf) (Gordon et al., 1993). However, proposing from the prior often leads to a poor approximation for a small number of particles, especially if xtx_{t} is high-dimensional. Variational smc addresses this shortcoming; it learns parameterized proposal distributions for efficient inference.

Variational Sequential Monte Carlo

We develop vsmc, a new class of variational approximations based on smc. We first define how to sample from the vsmc family and then derive its distribution. Though generating samples is straightforward, the density is intractable. To this end, we derive a tractable objective, a new lower bound to the elbo, that is amenable to stochastic optimization. Then, we present an algorithm to fit the variational parameters. Finally, we explore how to learn model parameters using variational expectation-maximization.

To sample from the vsmc family, we run smc (with the proposals parameterized by variational parameters λ\lambda) and then sample once from the empirical approximation of the posterior (2). Because the proposals r(xt∣xt−1 ;λ)r(x_{t}\mid x_{t-1}\,;\lambda) depend on λ\lambda, so does the smc empirical approximation. Algorithm 1 summarizes the generative process for the vsmc family.

The variational distribution q(x1:T ;λ)q(x_{1:T}\,;\lambda) marginalizes out all the variables produced in the sampling process, save for the output sample x1:Tx_{1:T}. This marginal comes from the joint distribution of all variables generated by vsmc,

(We have annotated this equation with the steps from the algorithm.) In this joint, the final output sample is defined by extracting the bTb_{T}-th trajectory x1:T=x1:TbT{x_{1:T}=x_{1:T}^{b_{T}}}. Note that the data y1:Ty_{1:T} enter via the weights and (optionally) the proposal distribution. This joint density is easy to calculate, but for variational inference we need the marginal distribution of x1:Tx_{1:T}. We derive this next.

Let bt≜atbt+1b_{t}\triangleq a_{t}^{b_{t+1}} for t≤T−1t\leq T-1 denote the ancestors for the trajectory x1:Tx_{1:T} returned by Algorithm 1. Furthermore, let ¬b1:T\neg b_{1:T} be all particle indices not equal to (b1,…,bT)(b_{1},\ldots,b_{T}), i.e. exactly all the particles that were not returned by Algorithm 1. Then the marginal distribution of x1:T=x1:Tb1:T=(x1b1,x2b2,…,xTbT)x_{1:T}=x_{1:T}^{b_{1:T}}=(x_{1}^{b_{1}},x_{2}^{b_{2}},\ldots,x_{T}^{b_{T}}) is given by the following proposition.

This has an intuitive form: the density of the variational posterior is equal to the exact joint times the expected inverse of the normalization constant (c.f. (3)). While we can estimate this expectation with Monte Carlo, it yields a biased estimate of log⁡q(x1:T ∣ y1:T;λ){\log q(x_{1:T}\,|\,y_{1:T};\lambda)} and the elbo (1).

To derive a tractable objective, we develop a lower bound to the elbo that is also amenable to stochastic optimization. It is

We call L~(λ)\widetilde{\mathcal{L}}(\lambda) the surrogate elbo. It is a lower bound to the true elbo for vsmc or, equivalently, an upper bound on the kl divergence. The following theorem formalizes this fact:

The surrogate elbo (6), is a lower bound to the elbo (1) when qq is defined by (5), i.e.

The surrogate elbo is the expected smc log-marginal likelihood estimate. We can estimate it unbiasedly as a byproduct of sampling from the vsmc variational approximation (Algorithm 1). We run the algorithm and use the estimate to perform stochastic optimization of the surrogate elbo.

Stochastic Optimization.

While the expectations in the surrogate elbo are still not available in closed form, we can estimate it and its gradients with Monte Carlo. This admits a stochastic optimization algorithm for finding the optimal variational parameters of the vsmc family.

We assume the proposals r(xt ∣ xt−1;λ){r(x_{t}\,|\,x_{t-1};\lambda)} are reparameterizable, i.e., we can simulate from rr by setting xt=h(xt−1,εt ;λ), εt∼s(εt){x_{t}=h(x_{t-1},\varepsilon_{t}\,;\lambda),~{}\varepsilon_{t}\sim s(\varepsilon_{t})} for some distribution ss not a function of λ\lambda. With this assumption, rewrite the gradient of (6) by using the reparameterization trick (Kingma and Welling, 2014; Rezende et al., 2014),

This expansion follows from the product rule, just as in the generalized reparameterizations of Ruiz et al. (2016) and Naesseth et al. (2017). Note that all xtix_{t}^{i}, implicit in the weights wtiw_{t}^{i} and p^(y1:T)\widehat{p}(y_{1:T}) are now replaced with their reparameterizations h(⋅ ;λ)h(\cdot\,;\lambda). The ancestor variables are discrete and cannot be reparameterized—this can lead to high variance in the score function term, gscoreg_{\text{score}} from (7).

In Section 5, we empirically assess the impact of ignoring gscoreg_{\text{score}} for optimization. We empirically study optimizing with and without the score function term for a small state space model where standard variance reduction techniques, explained below, are sufficient. We lower the variance using Rao-Blackwellization (Robert and Casella, 2004; Ranganath et al., 2014), noting that the ancestor variables at−1a_{t-1} have no effect on weights prior to time tt,

Furthermore, we use the score function ∇log⁡ϕ~(a1:T−11:N ∣ ε1:T1:N ;λ)\nabla\log\widetilde{\phi}(a_{1:T-1}^{1:N}\,|\,\varepsilon_{1:T}^{1:N}\,;\lambda) with an estimate of the future log average weights as a control variate (Ranganath et al., 2014).

We found that ignoring the score function term gscoreg_{\text{score}} (8) from the ancestor variables, leads to faster convergence and very little difference in final elbo. This corresponds to approximating the gradient of L~\widetilde{\mathcal{L}} by

This is the gradient we propose to use for optimizing the variational parameters of vsmc. See the supplementary material A.3 for more details, where we also provide a general score function-like estimator and the control variates.

Algorithm.

We now describe the full algorithm to optimize the vsmc variational approximation. We form stochastic gradients ∇^L~(λ)\widehat{\nabla}\widetilde{\mathcal{L}}(\lambda) by estimating (9) using a single sample from s(⋅)ϕ~(⋅ ∣ ⋅ ;λ)s(\cdot)\widetilde{\phi}(\cdot\,|\,\cdot\,;\lambda). The sample is obtained as a byproduct of sampling vsmc (Algorithm 1). We use the step-size sequence Adam (Kingma and Ba, 2015) or ρn\rho^{n} proposed by Kucukelbir et al. (2017),

where nn is the iteration number. We set δ=10−16\delta=10^{-16} and t=0.1t=0.1, and we try different values for η\eta. Algorithm 2 summarizes this optimization algorithm.Reference implementation using Adam is available at github.com/blei-lab/variational-smc.

Variational Expectation Maximization.

Suppose the target distribution of interest p(x1:T ∣ y1:T ;θ){p(x_{1:T}\,|\,y_{1:T}\,;\theta)} has a set of unknown parameters θ\theta. We can fit the parameters using variational expectation-maximization (vem) (Beal and Ghahramani, 2003). The surrogate elbo is updated accordingly

Perspectives on Variational smc

We give some perspectives on vsmc. First, we consider the vsmc special cases of N=1{N=1} and T=1{T=1}. For N=1N=1, vsmc reduces to a structured variational approximation: there is no resampling and the variational distribution is exactly the proposal. For T=1T=1, vsmc leads to a special case we call variational importance sampling, and a reinterpretation of the iwae (Burda et al., 2016), which we explore further in the first half of this section.

Then, we think of sampling from vsmc as sampling a highly optimized smc approximation. This means many of the theoretical smc results developed over the past 25 years can be adapted for vsmc. We explore some examples in the second half of this section.

The case where T=1{T=1} is smc without any resampling, i.e., importance sampling. The corresponding special case of vsmc is vis. The surrogate elbo for vis is exactly equal to the iwae lower bound (Burda et al., 2016).

This equivalence provides new intuition behind the iwae’s variational approximation on the latent variables. If we want to make use of the approximation q(x1:T ;λ⋆)q(x_{1:T}\,;\lambda^{\star}) learned with the iwae lower bound, samples from the latent variables should be generated with Algorithm 1, i.e. vis.

For vis it is possible to show that the surrogate elbo is always tighter than the one obtained by standard vb (equivalent to vis with N=1N=1) (Burda et al., 2016). This result does not carry over to vsmc, i.e. we can find cases when the resampling creates a looser bound compared to standard vb or vis. However, in practice the vsmc lower bound outperforms the vis lower bound.

Figure 2 provides a simple example of vis applied to a multimodal p(x ∣ y)∝N(x ;0,1) N(y ;x2/2,ex/2){p(x\,|\,y)\propto\mathcal{N}(x\,;0,1)\,\mathcal{N}(y\,;x^{2}/2,e^{x/2})} with a normal proposal r(x ;λ)=N(x ;μ,σ2){r(x\,;\lambda)=\mathcal{N}(x\,;\mu,\sigma^{2})} and a kernel density estimate of the corresponding variational approximation q(x ;λ)q(x\,;\lambda). The number of particles is N=10N=10. Standard vb with a Gaussian approximation only captures one of the two modes; which one depends on the initialization. We see that even a simple proposal can lead to a very flexible posterior approximation. This property is also inherited by the more general T>1T>1 case, vsmc.

Theoretical Properties.

This fact means that the gap in Theorem 1 disappears and the distribution of the trajectory returned by vsmc will tend to the true target distribution p(x1:T ∣ y1:T){p(x_{1:T}\,|\,y_{1:T})}. A bound on the kl divergence gives us the rate

for some constant c(λ)<∞c(\lambda)<\infty. This is a special case of a “propagation of chaos” result from Del Moral (2004, Theorem 8.3.2).

We can arrive at this result informally by studying (5): as the number of particles increases, the marginal likelihood estimate will converge to the true marginal likelihood and the variational posterior will converge to the true posterior. Huggins and Roy (2017) provide further bounds on various divergences and metrics between smc and the target distribution.

vsmc and T𝑇T.

Like smc, variational sequential Monte Carlo scales well with TT. Bérard et al. (2014) show a central limit theorem for the smc approximation log⁡p^(y1:T)−log⁡p(y1:T){\log\widehat{p}(y_{1:T})-\log p(y_{1:T})} with N=bT{N=bT}, where b>0{b>0}, as T→∞{T\to\infty}. Under the same conditions as in that work, and assuming that log⁡p^(y1:T){\log\widehat{p}(y_{1:T})} is uniformly integrable, we can show that

The implication for vsmc is significant. We can make the variational approximation arbitrarily accurate by setting N∝T{N\propto T}, even as TT goes to infinity. The supplement shows that this holds in practice; see A.4 for the toy example from Figure 1. We emphasize that neither standard vb nor iwae (vis) have this property.

Empirical Study

The linear Gaussian ssm is a ubiquitous model of time series data that enjoys efficient algorithms for computing the exact posterior. We use this model to study the convergence properties and impact of biased gradients for vsmc. We further use it to confirm that we learn good proposals. We compare to the bootstrap particle filter (bpf), which uses the prior as a proposal, and the (locally) optimal proposal that tilts the prior with the likelihood.

where vt∼N(0,Q)v_{t}\sim\mathcal{N}(0,Q), et∼N(0,R)e_{t}\sim\mathcal{N}(0,R), and x1∼N(0,I)x_{1}\sim\mathcal{N}(0,I). The log-marginal likelihood log⁡p(y1:T)\log p(y_{1:T}) can be computed using the Kalman filter.

We study the impact of the biased gradient (9) for optimizing the surrogate elbo (6). First, consider a simple scalar model with A=0.5{A=0.5}, Q=1{Q=1}, C=1{C=1}, R=1{R=1}, and T=2{T=2}. For the proposal we use r(xt ∣ xt−1 ;λ)=N(xt ;λ+0.5xt−1,1){r(x_{t}\,|\,x_{t-1}\,;\lambda)=\mathcal{N}(x_{t}\,;\lambda+0.5x_{t-1},1)}, with x0≡0{x_{0}\equiv 0}. Figure 3 (left) shows the mean and spread of estimates of gscoreg_{\text{score}} (8), with control variates, and grepg_{\text{rep}} (9), as a function of λ\lambda for four randomly generated datasets. The optimal setting of λ\lambda is where the sum of the means is equal to zero. Ignoring the score function term gscoreg_{\text{score}} (8) will lead to a perturbation of the optimal λ\lambda. However, even for this simple model, the variance of the score function term (red) is several orders of magnitude higher than that of the reparameterization term (blue), despite the variance reduction techniques of Section 3. This variance has a significant impact on the convergence speed of the stochastic optimization.

Next, we study the magnitude of the perturbation, and its effect on the surrogate elbo. We generate data with T=10T=10, (A)ij=α∣i−j∣+1(A)_{ij}=\alpha^{|i-j|+1} for α=0.42\alpha=0.42, Q=IQ=I, and R=IR=I. We explored several settings of dx=dim⁡(xt)d_{x}=\dim(x_{t}), dy=dim⁡(yt)d_{y}=\dim(y_{t}), and CC. Sparse CC measures the first dyd_{y} components of xtx_{t}, and dense CC has randomly generated elements Cij∼N(0,1)C_{ij}\sim\mathcal{N}(0,1). Figure 3 (right) shows the true log-marginal likelihood and elbo as a function of iteration. It shows vsmc with biased gradients (blue) and unbiased gradients (red). We choose the proposal

with λ={μt,βt,σt2}t=1T\lambda=\left\{\mu_{t},\beta_{t},\sigma_{t}^{2}\right\}_{t=1}^{T}, and set the number of particles to N=4N=4. Note that while the gradients are biased, the resulting elbo is not. We can see that the final vsmc elbo values are very similar, regardless of whether we train with biased or unbiased gradients. However, biased gradients converge faster. Thus, we use biased gradients in the remainder of our experiments.

Next, we study the effect of learning the proposal using vsmc compared with standard proposals in the smc literature. The most commonly used is the bpf, sampling from the prior ff. We also consider the so-called optimal proposal, r∝f⋅gr\propto f\cdot g, which minimizes the variance of the incremental importance weights (Doucet and Johansen, 2009). Table 1 shows results for a linear Gaussian ssm when T=25{T=25}, Q=0.12I{Q=0.1^{2}I}, R=1{R=1}, dx=10{d_{x}=10}, and dy=1{d_{y}=1}. Because of the relatively high-dimensional state, bpf exhibits significant bias whereas the optimal proposal smc performs much better. vsmc outperforms them both, learning an accurate proposal that results in an elbo only 0.90.9 nats lower than the true log-marginal likelihood. We further emphasize that the optimal proposal is unavailable for most models.

Stochastic Volatility

A common model in financial econometrics is the (multivariate) stochastic volatility model (Chib et al., 2009). The model is

where vt∼N(0,Q){v_{t}\sim\mathcal{N}(0,Q)}, et∼N(0,I){e_{t}\sim\mathcal{N}(0,I)}, x1∼N(μ,Q){x_{1}\sim\mathcal{N}(\mu,Q)}, and θ=(μ,ϕ,Q,β){\theta=(\mu,\phi,Q,\beta)}. (In the multivariate case, multiplication is element-wise.) Computing log⁡p(y1:T ;θ){\log p(y_{1:T}\,;\theta)} and its gradients for this model is intractable, we study the vem approximation to find the unknown parameters θ\theta. We compare vsmc with iwae and structured vi. For the proposal in vsmc and iwae we choose

with variational parameters λ=(μ1,…,μT,Σ1,…,ΣT){\lambda=(\mu_{1},\ldots,\mu_{T},\Sigma_{1},\ldots,\Sigma_{T})}. We define the variational approximation for structured vi to be q(x1:T ;λ,θ)=∏t=1Tr(xt ∣ xt−1 ;λ,θ)q(x_{1:T}\,;\lambda,\theta)=\prod_{t=1}^{T}r(x_{t}\,|\,x_{t-1}\,;\lambda,\theta).

We study 1010 years of monthly returns (9/20079/2007 to 8/20178/2017) for the exchange rate of 2222 international currencies with respect to US dollars. The data is from the Federal Reserve System. Table 2 reports the optimized elbo (higher is better) for different settings of the number of particles/samples N={4,8,16}N=\{4,8,16\}. vsmc outperforms the competing methods with almost 0.20.2 nats per time-step.

In theory we can improve the bound of both iwae and vsmc by increasing the number of samples NN. This means we can first learn proposals using only a few particles NN, for computational efficiency. Then, at test time, we can increase NN as needed for improved accuracy. We study the impact of increasing the number of samples for vsmc and iwae using fix θ⋆\theta^{\star} and λ⋆\lambda^{\star} optimized with N=16N=16. Figure 4 shows that the gain for iwae is limited, whereas for vsmc it can be significant.

Deep Markov Model

An important problem in neuroscience is understanding dynamics of neural circuits. We study a population of 105 motor cortex neurons simultaneously recorded in a macaque monkey as it performed reaching movements (c.f. Gao et al., 2016). In each trial, the monkey reached toward one of fourteen targets; each trial is T=21{T=21} time steps long. We train on 700{700} trials and test on 8484.

We use recurrent neural networks to model both the dynamics and observations. The model is

where vt∼N(0,I)v_{t}\sim\mathcal{N}(0,I), x0≡0x_{0}\equiv 0, and μ,σ,η\mu,\sigma,\eta are neural networks parameterized by θ\theta. The multiplication in the transition dynamics is element-wise. This is a deep Markov model (Krishnan et al., 2017).

For inference we use the following proposal for both vsmc and iwae,

where μx,σx,μy,σy\mu^{x},\sigma^{x},\mu^{y},\sigma^{y} are neural networks parameterized by λ\lambda, and the proposal factorizes over the components of xtx_{t}.

Figure 5 illustrates the result for dx={3,5,10}{d_{x}=\{3,5,10\}} with N=8{N=8}. vsmc gets to the same elbo faster.

Conclusions

We introduced the variational sequential Monte Carlo (vsmc) family, a new variational approximating family that provides practitioners with a flexible, accurate, and powerful approximate Bayesian inference algorithm. vsmc melds variational inference (vi) and sequential Monte Carlo (smc). This results in a variational approximation that lets us trade-off fidelity to the posterior with computational complexity.

Acknowledgements

Christian A. Naesseth is supported by CADICS, a Linnaeus Center, funded by the Swedish Research Council (VR). Scott W. Linderman is supported by the Simons Foundation SCGB-418011. This work is supported by ONR N00014-11-1-0651, DARPA PPAML FA8750-14-2-0009, the Alfred P. Sloan Foundation, and the John Simon Guggenheim Foundation.

References

Appendix A Variational Sequential Monte Carlo – Supplementary Material

We start by noting that the distribution of all random variables generated by the vsmc algorithm is given by

We insert the above expression in (13) and we get

A.2 Proof of Theorem 1

The evidence lower bound (elbo), using the above result about the distribution of q(x1:T ;λ)q(x_{1:T}\,;\lambda), is given by

where the last step follows because q(x1:T ;λ)q(x_{1:T}\,;\lambda) is the marginal of ϕ~(x1:T1:N,a1:T−11:N ;λ)\widetilde{\phi}(x_{1:T}^{1:N},a_{1:T-1}^{1:N}\,;\lambda). □\square

A.3 Stochastic Optimization

In practice we use a stochastic estimate of ctc_{t}.

For T=2T=2 we can use a leave-one-out estimator of the ancestor variable score function gradient

Below we provide the derivation of a score function-like estimator that is applicable in very general settings. However, we have found that in practice the variance tends to be quite high.

A.4 Scaling With Dimension

In this section we study how the methods compare on a simple toy model defined by