Year: 2020

Abstract:
Effective variational inference crucially depends on a flexible variational family of distributions. Recent work has explored sequential Monte-Carlo (SMC) methods to construct variational distributions, which can, in principle, approximate the target posterior arbitrarily well, which is especially appealing for models with inherent sequential structure. However, SMC, which represents the posterior using a weighted set of particles, often suffers from particle weight degeneracy, leading to a large variance of the resulting estimators. To address this issue, we present a novel approach that leverages the idea of partial rejection control (PRC) for developing a robust variational inference (VI) framework. In addition to developing a superior VI bound, we propose a novel marginal likelihood estimator constructed via a dice-enterprise: a generalization of the Bernoulli factory to construct unbiased estimators for SMC-PRC. The resulting variational lower bound can be optimized efficiently with respect to the variational parameters and generalizes several existing approaches in the VI literature into a single framework. We show theoretical properties of the lower bound and report experiments on various sequential models, such as the Gaussian state-space model and variational RNN, on which our approach outperforms existing methods.

heading: 1 INTRODUCTION
text: Exact inference in latent variable models is usually intractable. Markov Chain Monte-Carlo (MCMC) (Andrieu et al., 2003) and variational inference (VI) methods (Blei et al., 2017), are commonly employed in such models to make inference tractable. While MCMC has been the traditional method of choice, often with provable guarantees, optimization-based VI methods have also enjoyed considerable recent interest due to their excellent scalability on large-scale datasets. VI is based on maximizing a lower bound constructed through a marginal likelihood estimator. For latent variable models with sequential structure, sequential Monte-Carlo (SMC) (Doucet & Johansen, 2009) returns a much lower variance estimator of the log marginal likelihood than importance sampling (Bérard et al., 2014; Cérou et al., 2011). In this work, we focus our attention on designing a low variance, unbiased, and computationally efficient estimator of the marginal likelihood.
The performance of SMC based methods is strongly dependent on the choice of the proposal distribution. Inadequate proposal distributions propose values in low probability areas under the target, leading to particle depletion (Doucet & Johansen, 2009). An effective solution is to use rejection control (Liu et al., 1998; Peters et al., 2012) which is based on an approximate rejection sampling step within SMC to reject samples with low importance weights.
In this work, we leverage the idea of partial rejection control (PRC) within the framework of SMC based VI for sequential latent variable models. To this end, we construct a novel lower bound, VSMC-PRC, and propose an efficient optimization strategy for selecting the variational parameters. Compared to other recent SMC based VI approaches (Naesseth et al., 2017; Maddison et al., 2017; Le et al., 2017), our approach consists of an inbuilt accept-reject mechanism within SMC to prevent particle depletion. The use of accept-reject within SMC makes the particle weight intractable, therefore, we use a generalization of the Bernoulli factory (Asmussen et al., 1992) to construct unbiased estimators of the marginal likelihood for SMC-PRC.
Although the idea of combining VI with an inbuilt accept-reject mechanism is not new (Salimans et al., 2015; Ruiz & Titsias, 2019; Grover et al., 2018; Gummadi, 2014), a key distinction of our approach is to incorporate an accept-reject mechanism along with a resampling framework. In contrast to standard sampling algorithms that may reject the entire stream of particles, we use a partial accept-reject on the most recent update, increasing the sampling efficiency. Further, the variational framework of SMC-PRC is interesting in itself as it combines accept-reject with particle filter methods. Therefore, our proposed bound VSMC-PRC generalizes several existing approaches for example: Variational Rejection Sampling (VRS) (Grover et al., 2018), FIVO (Maddison et al., 2014), IWAE (Burda et al., 2015), and standard variational Bayes (Blei et al., 2017).
Another key distinction is that, while existing approaches using Bernoulli factory are limited to niche one-dimensional toy examples, our proposed approach is scalable. To the best of our knowledge, there is no prior work that has used Bernoulli factories for such a general case like variational recurrent neural networks (VRNN); therefore, we believe this aspect to be a significant contribution as well. The rest of the paper is organized as follows: In Section 2, we provide a brief review on SMC, partial rejection control, and dice enterprise. In Section 3, we introduce our VSMC-PRC bound and provide new theoretical insights into the Monte-Carlo estimator and design efficient ways to optimize it. Finally, we discuss related work and present experiments on the Gaussian state-space model (SSM) and VRNN.

heading: 2 BACKGROUND
text: We denote a sequence of T real-valued observations as x1:T = (x1, x2, . . . , xT ), and assume that there is an associated sequence of latent variables z1:T = (z1, z2, . . . , zT ). We are interested in inferring the posterior distribution of the latent variables, i.e., p(z1:T |x1:T ). The task is, in general, intractable. For the rest of the paper we have used some common notations from SMC and VI literature where zit: i th particle at time t; Ait−1: ancestor variable for the i th particle at time t; θ and φ: model and variational parameters, respectively.

heading: 2.1 SEQUENTIAL MONTE CARLO WITH PARTIAL REJECTION CONTROL
text: An SMC sampler approximates a sequence of densities {pθ(z1:t|x1:t)}Tt=1 through a set of N weighted samples generated from a proposal distribution. Let the proposal density be
qφ(z1:T |x1:T ) = T∏ t=1 qφ(zt|x1:t, z1:t−1). (1)
Consider time t−1 at which we have uniformly weighted samples {N−1, zi1:t−1, Ait−1}Ni=1 estimating pθ(z1:t−1|x1:t−1). We want to estimate pθ(z1:t|x1:t) such that particles with a low importance weight are automatically rejected. PRC achieves this by using an approximate rejection sampling step (Liu et al., 1998; Peters et al., 2012). The overall procedure is as follows:
1. Generate zit ∼ qφ(zt|x1:t, z Ait−1 1:t−1) where i = 1, 2, . . . , N .
2. Accept zit with probability
aθ,φ(z i t|z Ait−1 1:t−1, x1:t) = 1 + M(i, t− 1)qφ(zit|x1:t, zAit−11:t−1) pθ(xt, zit|x1:t−1, z Ait−1 1:t−1) −1 , (2) where M(i, t− 1) is a hyperparameter controlling the acceptance rate (see Proposition 3 and Section 3.3 for more details). Note that PRC applies accept-reject only on zit , not on the entire trajectory.
3. If zit is rejected go to step 1. 4. The new incremental importance weight of the accepted sample is
αt(z i 1:t) = c i tZ(z Ait−1 1:t−1, x1:t), (3)
where cit is
cit = pθ(xt, z
i t|x1:t−1, z Ait−1 1:t−1)
qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|z Ait−1 1:t−1, x1:t)
, (4)
and the intractable normalization constant Z(.) (For simplicity of notation, we have ignored the dependence of Z(.) on M(i, t− 1))
Z(z Ait−1 1:t−1, x1:t) = ∫ aθ,φ(zt|z Ait−1 1:t−1, x1:t)qφ(zt|x1:t, z Ait−1 1:t−1)dzt. (5)
5. Compute Monte-Carlo estimator of unnormalized weight
w̃it = pθ(xt, z
i t|x1:t−1, z Ait−1 1:t−1) 1 K ∑K k=1 aθ,φ(δ i,k t |z Ait−1 1:t−1, x1:t)
qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|z Ait−1 1:t−1, x1:t)
, (6)
where δi,kt ∼ qφ(zt|x1:t, z Ait−1 1:t−1) and k = 1, 2, . . . ,K. Note that w̃ i t would be essential for
constructing unbiased estimator for pθ(x1:T ). 6. Generate ancestor variables Ait through dice-enterprise and set new weights w i t = N
−1 for i = 1, 2, . . . , N :
Ait ∼ Categorical
( αt(z
1 1:t)∑N
j=1 αt(z j 1:t)
, αt(z 2 1:t)∑N
j=1 αt(z j 1:t)
, . . . , αt(z N 1:t)∑N
j=1 αt(z j 1:t)
) . (7)

heading: 2.2 DICE ENTERPRISE
text: Simulation of ancestor variables in Eq. 7 is non-trivial due to intractable normalization constants in the incremental importance weight (see Eq. 3). Vanilla Monte-Carlo estimation of αt(.) yields biased samples of ancestor variables from Eq. 7. To address this issue, we leverage a generalization of Bernoulli factory, called dice-enterprise (Morina et al., 2019). Note that multinoulli extensions of Bernoulli factory (Dughmi et al., 2017) have also been used for resampling within intractable SMC before (Schmon et al., 2019), a key distinction of our approach is to design a scalable Bernoulli factory methodology especially useful for VI applications.
Suppose we can simulate Bernoulli(pit) outcomes where p i t is intractable. Bernoulli factory problem simulates an event with probability f(pit), where f(.) is some desired function. In our case, the intractable coin probability pit is the intractable normalization constant,
pit = Z(z Ait−1 1:t−1, x1:t) = ∫ aθ,φ(zt|z Ait−1 1:t−1, x1:t)qφ(zt|x1:t, z Ait−1 1:t−1)dzt. (8)
Since pit ∈ [0, 1] and we can easily simulate this coin, we obtain the dice-enterprise algorithm below.
1. Required: Constants {cit}Ni=1 see Eq. 4. 2. Sample C ∼ Categorical (
c1t∑N j=1 c j t , c2t∑N j=1 c j t , . . . , cNt∑N j=1 c j t ) 3. If C = i, generate Ui ∼ U [0, 1] and zt ∼ qφ(zt|x1:t, z Ait−1 1:t−1)
• If Ui < aθ,φ(zt|z Ait−1 1:t−1, x1:t) output i • Else go to step 2
The dice-enterprise produces unbiased ancestor variables. Note that we can easily control the efficiency of the proposed dice-enterprise through the hyper-parameter M (similar to Eq. 2) in contrast to existing Bernoulli factory algorithms (Schmon et al., 2019). For details on efficiency and correctness, please refer to Section 3.1 and Section 3.3.
Our proposed VSMC-PRC bound is constructed through a marginal likelihood estimator obtained by combining the SMC sampler with a PRC step and dice-enterprise. The variance of estimators obtained through SMC-PRC particle filter is usually low (Peters et al., 2012). Therefore, we expect VSMC-PRC to be a tighter bound compared to the standard SMC based bounds used in recent works (Maddison et al., 2017; Naesseth et al., 2017; Le et al., 2017). Algorithm 1 summarizes the generative process to simulate the VSMC-PRC bound. Please see Figure 1 to visualize the generative process for VSMC-PRC.
Algorithm 1 Estimating the VSMC-PRC lower bound
1: Required: N , K, and M 2: for t ∈ {1, 2, . . . , T} do 3: for i ∈ {1, 2, . . . , N} do 4: zit, c i t, w̃ i t ∼ PRC (q, p,M(i, t− 1)) 5: zi1:t = (z Ait−1 1:t−1, z i t) 6: end for 7: for i ∈ {1, 2, . . . , N} do 8: Ait = DICE-ENT ( {cit, zi1:t}Ni=1
) 9: end for
10: end for 11: return log ∏T t=1 ( 1 N ∑N i=1 w̃ i t ) 12: 13: PRC (q, p,M(i, t− 1)) 14: while sample not accepted do 15: Generate zit ∼ qφ(zt|x1:t, z Ait−1 1:t−1) 16: Accept zit with probability
aθ,φ(z i t|z Ait−1 1:t−1, x1:t)
17: end while 18: Sample {δi,kt }Kk=1 ∼ qφ(zt|x1:t, z Ait−1 1:t−1) 19: Calculate w̃it from Eq. 6 20: Calculate cit from Eq. 4 21: return ( zit, c i t, w̃ i t
) 22: 23: DICE-ENT ( {cit, zi1:t}Ni=1
) 24: Sample C ∼ Multinoulli ( cit∑N j=1 c j t )N i=125: if C == i then 26: Sample Ui ∼ U [0, 1] 27: zit ∼ qφ(zt|x1:t, z Ait−1 1:t−1) 28: end if 29: if Ui < aθ,φ(zit|z Ait−1 1:t−1, x1:t) then 30: return (i) 31: else 32: return DICE-ENT ( {cit, zi1:t}Ni=1
) 33: end if

heading: 3 PARTIAL REJECTION CONTROL BASED VI FOR SEQUENTIAL LATENT VARIABLE MODELS
text: We now show how to leverage PRC to develop a robust VI framework for sequential latent variable models. Our framework is based on the VSMC-PRC bound presented below. The complete sampling distribution of Algorithm 1 is as follows.
QVSMC-PRC ( z1:N1:T , A 1:N 1:T−1, δ 1:N,1:K 1:T ) = ( K∏ k=1 N∏ i=1 qφ(δ i,k 1 |x1) T∏ t=2 N∏ i=1 K∏ k=1 qφ(δ i,k t |x1:t, z Ait−1 1:t−1) ) ×
( N∏ i=1 qφ(z i 1|x1)aθ,φ(zi1|x1) Z(x1) T−1∏ t=1 N∏ i=1 Discrete(Ait|αt) qφ(z i t+1|x1:t+1, z Ait 1:t)aθ,φ(z i t+1|z Ait 1:t , x1:t+1) Z(z Ait 1:t , x1:t+1) ) (9)
The normalization constants Z(.) in Eq. 9 are intractable and have to be estimated while calculating the weights. Therefore, we introduce an extra parameter K, denoting the number of Monte-Carlo samples used to estimate Z(.). The Monte-Carlo estimator of VSMC-PRC bound is
L̂VSMC-PRC(θ, φ;x1:T ,K) = T∑ t=1 log
( 1
N N∑ i=1 w̃it
) . (10)
We maximize the VSMC-PRC bound with respect to model parameters θ and variational parameters φ. This requires estimating the gradient the details of which are provided in Section 3.2.

heading: 3.1 THEORETICAL PROPERTIES
text: We now present properties of the Monte-Carlo estimator L̂VSMC-PRC. The key variables that affect this bound are N (number of samples), hyper-parameter M , and the number of Monte-Carlo samples used to compute the normalization constant Z(.), i.e., K. As discussed by Bérard et al. (2014); Naesseth et al. (2017), as N increases, we expect the VSMC-PRC bound to get tighter. Hence, we will focus our attention on M and K. All the proofs can be found in the appendix. Proposition 1. Dice-enterprise produces unbiased ancestor variables. Further, let Λt be the number of iterations required for generating one ancestor variable, then Λt ∼ Geom ( E[Λt]−1 ) where
E[Λt] = ∑N i=1 c
i t∑N
i=1 c i tZ(z Ait−1 1:t−1, x1:t)
.
As evident from Proposition 1, the computational efficiency of the dice-enterprise clearly relies on the normalization constantZ(.). Note that the value ofZ(.) could be interpreted as the average acceptance rate of PRC which depends on the hyper-parameter M(i, t− 1). If the average acceptance rate for PRC for all particles is γ, then we can express the expected number of iterations as E[Λit] = γ−1. Therefore, the computational efficiency of dice-enterprise is similar to the PRC step and depends crucially on the hyper-parameter M .
Proposition 2. For all K, exp(L̂VSMC-PRC) is unbiased, i.e., E [ exp(L̂VSMC-PRC) ] = pθ(x1:T ). Fur-
ther, E[L̂VSMC-PRC] is non-decreasing in K. The use of Monte-Carlo estimator in place of the true value of Z(.) creates an inefficiency, as depicted by Proposition 2. The bound monotonically increases as we increase K despite the use of resampling operation. It is important to note that Algorithm 1 produces an unbiased estimator of the marginal likelihood for all values of K. Proposition 3. Let the sampling distribution of the ith particle (generated via PRC) at time t be rθ,φ(zt|z Ait−1 1:t−1, x1:t), then
KL ( rθ,φ(zt|z Ait−1 1:t−1, x1:t) ‖ pθ(zt|z Ait−1 1:t−1, x1:t) ) ≤ KL ( qφ(zt|z Ait−1 1:t−1, x1:t) ‖ pθ(zt|z Ait−1 1:t−1, x1:t) ) .
Proposition 3 implies that the use of the accept-reject mechanism within SMC refines the sampling distribution. Instead of accepting all samples, the PRC step ensures that only high-quality samples are accepted, leading to a tighter bound for VSMC in general (not always). We show in the appendix that when M(i, t − 1)→∞, the PRC step reduces to pure rejection sampling (Robert & Casella, 2013). On the other hand, M(i, t− 1)→ 0 implies that all samples are accepted from the proposal. Recall, M(i, t− 1) is a hyperparameter that can be tuned to control the acceptance rate. For more details on tuning M , see Section 3.3.

heading: 3.2 GRADIENT ESTIMATION
text: For tuning the variational parameters, we use stochastic optimization. Algorithm 1 produces the marginal likelihood estimator by sequentially sampling the particles, ancestor variables, and particles for the normalization constant (z1:N1:T , A 1:N 1:T−1, δ 1:N,1:K 1:T ).
When the variational distribution qφ(.) is reparameterizable, we can make the sampling of δ i,k t independent of the model and variational parameters. However, the generated particles zit are not reparametrizable due to the PRC step. Finally, the ancestor variables are discrete and, therefore, cannot be reparameterized. The complete gradient can be divided into three core components (assuming qφ(.) is reparametrizable):
∇θ,φE[L̂VSMC-PRC] = EQVSMC-PRC [ ∇θ,φL̂VSMC-PRC(θ, φ;x1:T ,K) ] + gPRC + gRSAMP (11)
≈ EQVSMC-PRC [ ∇θ,φL̂VSMC-PRC(θ, φ;x1:T ,K) ] . (12)
Note that gPRC and gRSAMP denote the score gradient of PRC and resampling step, respectively. Due to high variance, we have ignored these terms for the optimization. We have derived the full gradient and explored the gradient variance issues in the appendix. Please see Figure 2 (left) comparing the convergence of biased gradient vs. unbiased gradients on a toy Gaussian SSM.

heading: 3.3 LEARNING THE M MATRIX
text: We use M as a hyperparameter for the PRC step which controls the acceptance rate of the sampler. The basic scheme of tuning M is as follows:
• Define a new random variable F (zt+1|z Ait 1:t , x1:t+1) = log
( qφ(zt+1|x1:t+1,z Ait 1:t )
pθ(xt+1,zt+1|x1:t,z Ait 1:t )
) .
• Draw zjt+1 ∼ qφ(zt+1|x1:t+1, z Ait 1:t) for j = 1, 2, . . . , J . • Evaluate γ ∈ [0, 1] quantile value of {F (zjt+1|z Ait t , x1:t+1)}Jj=1. In general for this case the
acceptance rate would be around γ for all particles.
logM(i, t) = −Q F (zt+1|z Ait 1:t ,x1:t+1) (γ). (13)
• If M matrix is very large then use a common {M(., t)}Tt=1 for every time-step. In general, for this configuration, the acceptance rate would be greater than equal to γ for all particles:
logM(., t) = min { −Q
F (zt+1|z Ait 1:t ,x1:t+1)
(γ) }N i=1 . (14)
Through γ: a user parameter, we can directly control the acceptance rate. Therefore, both diceenterprise and PRC would take around (less than) γ−1 iterations to produce a sample for M value learned from Eq. 13 (see Eq. 14). For implementation details please refer to the experiments.
Note that a similar scheme was also employed in Grover et al. (2018). We update {{M(i, t − 1)}Ni=1}Tt=1 dynamically once every F epochs to save time. To learn more on setting hyper-parameter M , see Liu et al. (1998); Peters et al. (2012).

heading: 4 RELATED WORK AND SPECIAL CASES
text: There is significant recent interest in developing more expressive variational posteriors for latent variable models. There are two basic schemes for constructing tighter bounds on the log marginal likelihood: sampling-based methods (MCMC, rejection sampling) (Salimans et al., 2015; Ruiz & Titsias, 2019; Hoffman, 2017; Grover et al., 2018) or multiple samples from VI distributions to increase the flexibility (IS, SMC) (Burda et al., 2015; Maddison et al., 2017; Lawson et al., 2018; Naesseth et al., 2015). In this work, we present a unified framework for combining these two approaches, utilizing the best of both worlds. Although applying sampling-based methods on VI is useful, the density ratio between the true posterior and the improved density is often intractable. Therefore, we cannot take advantage of variance-reducing schemes like resampling, which is crucial for sequential models. We solve this issue through dice-enterprise: an extension of the Bernoulli factory.
Recently, Bernoulli factory has amassed a great interest in the area of Bayesian inference (Gonçalves et al., 2017a;b; Vats et al., 2020). Although Bernoulli factory is theoretically valuable, its applicability is severely limited due to a high rejection rate. In this paper, we have presented an approach that combines SMC with dice-enterprise for efficient implementation. A closely related work from SMC literature is Schmon et al. (2019), which also utilizies Bernoulli factories to implement unbiased resampling. However, their method is not scalable and designed particularly for partially observed diffusions. Another relevant work for unbiased estimation of the marginal likelihood is that of Kudlicka et al. (2020). In contrast to our approach, this method samples one additional particle and keeps track of the number of steps required by PRC for every time-step to obtain their unbiased estimator. The weights are tractable for Kudlicka et al. (2020) as they do not take into account the effect of the normalization constant Z(.). On the other hand, we consider the effect of Z(.) on the particle’s weight, making resampling operation infeasible. To fix this intractability, we
use dice-enterprise. Comparisons of our marginal likelihood estimator versus that of Kudlicka et al. (2020) would make for interesting future work.
To provide more clarity, we will consider some special cases of VSMC-PRC bound and relate it with existing work: Note that for N = 1 our method reduces to a special case of Gummadi (2014) which uses a constraint function Ct for every time-step and restarts the particle trajectory from ∆t (if Ct is violated). Therefore, if we use the setting Ct(z1:t) = a(zt|z1:t−1, x1:t) and ∆t = t− 1, our method reduces to a specific case of Gummadi (2014). For the special case of N = 1 and T = 1, our method reduces to VRS (Grover et al., 2018). For N,T > 1, if we remove the PRC step, our bound reduces to FIVO (Maddison et al., 2017). Finally, if we remove both the PRC step and resampling, then our method effectively reduces to IWAE (Burda et al., 2015). Please refer to Figure 1 for more details.

heading: 5 EXPERIMENTS
text: In this section, we evaluate our proposed algorithm on synthetic as well as real-world datasets and compare them with relevant baselines. For the synthetic data experiment, we implement our method on a Gaussian SSM and compare our approach with VSMC (Naesseth et al., 2017). For the real data experiment, we train a VRNN (Chung et al., 2015) on the polyphonic music dataset.

heading: 5.1 GAUSSIAN STATE SPACE MODEL
text: In this experiment, we study the linear Gaussian state space model. Consider the model zt = Azt−1 + ez,
xt = Czt + ex,
where ez, ex ∼ N (0, I) and z0 = 0. We are interested in learning a good proposal for the above model. The latent variable is denoted by zt and the observed data by xt. Let the dimension of zt be dz and dimension of xt be dx. The matrix A has the elements (A)i,j = α|i−j|+1, for α = 0.42. We explore different settings of dz, dx, and matrix C. A sparse version of C matrix measures the first dx components of zt, on the other hand a dense version of C is normally distributed i.e Ci,j ∼ N (0, 1). We consider four different configurations for the experiment. For more details please refer to Figure 2.
The variational distribution is a multivariate Gaussian with unknown mean vector µ = {µd}dzd=1 and diagonal covariance matrix {log σ2d} dz d=1. We set N = 4 and T = 10 for all the cases:
q(zt|zt−1) ∼ N ( zt|Azt−1 + µ, diag(σ2) ) .
The {{M(i, t− 1)}Ni=1}Tt=1 matrix (see Eq. 13) for approximate rejection sampling is updated once every 10 epochs with acceptance rate γ ∈ {0.8, 0.4}. For estimating the intractable normalization constants, we generate K = 3 samples. Figure 2: (left) compares the convergence of biased gradient vs unbiased gradients. Note that we get a much tighter bound as compared to VSMC (Naesseth et al., 2017).

heading: 5.2 VARIATIONAL RNN
text: VRNN (Chung et al., 2015) comprises of three core components: the observation xt, stochastic latent state zt, and a deterministic hidden state ht(zt−1, xt−1, ht−1), which is modeled through a RNN. For the experiments, we use a single-layer LSTM for modeling the hidden state. The conditional distributions pt(zt|.) and qt(zt|.) are assumed to be factorized Gaussians, parametrized by a single layer neural net. The output distribution gt(xt|.) depends on the dataset. For a fair comparison, we use the same model setting as employed in FIVO (Maddison et al., 2017). We evaluate our model on four polyphonic music datasets: Nottingham, JSB chorales, Musedata, and Piano-midi.de.
Each observation xt is represented as a binary vector of 88 dimensions. Therefore, we model the observation distribution gt(xt|.) by a set of 88 factorized Bernoulli variables. We split all four datasets into the standard train, validation, and test sets. For tuning the learning rate, we use the validation test set. For a fair comparison, we use the same learning rate and iterations for all the models. Let the dimension of hidden state (learned by single layer LSTM) be dh and dimension of latent variable be dz . We choose the setting dz = dh = 64 for all the data-sets except JSB. For modeling JSB, we use dz = dh = 32. For VSMC-PRC we have considered N ∈ {4, 6} Further, for each N , we consider four settings (K, γ) ∈ {(1, 0.9), (1, 0.8), (3, 0.9), (3, 0.8)}. The M hyper-parameter for PRC step is learned from Eq. 14 due to large size. We have updated M value once every 50 epoches. Note that in this scenario, the acceptance rate for all particles would be greater than equal to γ. For more details on experiments, please refer to the appendix.
As discussed in Section 3.1, the PRC step and dice-enterprise have time complexity O(N/γ) for producing N samples (assuming average acceptance rate γ). Therefore, we consider dNγ−1e particles for IWAE and FIVO to ensure effectively the same number of particles, where N ∈ {4, 6} and γ = 0.8. Note, however, that the acceptance rate is ≥ γ, so this adjustment actually favors the other approaches more. For FIVO, we perform resampling when ESS falls below N/2. Table 1 summarizes the results which show whether rejecting samples provide us with any benefit or not, and as the results show, our approach, even with the aforementioned adjustment, outperforms the other approaches in terms of test log-likelihoods, while still having a similar computational cost.
In Sec. 3.1, we discussed the effect of K and PRC rejection rate on VSMC-PRC bound. We expect a performance improvement when K and the rejection rate is increased. Although the results for VSMC-PRC’s different configurations are almost the same, we still get the best average ranking for (K = 3, γ = 0.8). Overall, for most cases, VSMC-PRC bound performs better than FIVO (Maddison et al., 2017) and IWAE (Burda et al., 2015) for a variety of configurations.
In VSMC-PRC, improvement in the bound value comes at the cost of estimating the normalization constant Z(.), i.e., K. On further inspection, we can clearly see that increasing K doesn’t provide us
with any substantial benefits despite the increase in computational cost. Therefore, to maintain the computational trade-off (K = 1, γ > 0.8) seems to be a reasonable choice for VI practitioners.
Table 1 signifies that rejecting samples with low importance weight is better instead of keeping a large number of particles (at least for a reasonably high acceptance rate γ). The proposed bound uses more particles (PRC step and dice-enterprise) than existing approaches like FIVO and IWAE due to intractability. Future work aims at designing a scalable implementation for VSMC-PRC bound that consumes fewer particles.

heading: 6 CONCLUSION
text: We introduced VSMC-PRC, a novel bound that combines SMC and partial rejection sampling with VI in a synergistic manner. This results in a robust VI procedure for sequential latent variable models. Instead of using standard sampling algorithms, we have employed a partial sampling scheme suitable for high dimensional sequences. Our experimental results clearly demonstrate that VSMC-PRC outperforms existing bounds like IWAE (Burda et al., 2015) and standard particle filter bounds (Maddison et al., 2017; Naesseth et al., 2017; Le et al., 2017). The future work aims to explore partial versions of powerful sampling algorithms like Hamiltonian Monte Carlo (Neal et al., 2011) instead of rejection sampling.

heading: A PROOF OF THEORETICAL RESULTS
text: Proof of proposition 1 : Dice-enterprise produces unbiased ancestor variables. Let’s evaluate the probability that Dice-enterprise would output i as the ancestor index. Assume that the algorithm terminates after r steps where r ∈ {1, 2, . . . ,∞}, then the probability of output i is given as follows:
Pr(output = i) = ∞∑ r=1 Pr(output = i|after r steps)
= citZ(z Ait−1 1:t−1, x1:t)∑N m=1 c m t ∞∑ r=1  N∑ j=1 cjt (1− Z(z Ajt−1 1:t−1, x1:t)∑N m=1 c m t r−1
= citZ(z Ait−1 1:t−1, x1:t)∑N
j=1 c j tZ(z Ajt−1 1:t−1, x1:t)
.
It is easy to see that Λt is geometrically distributed with success probability given by Pr(getting output in a loop) = ∑N i=1 c i tZ(z Ait−1 1:t−1, x1:t)∑N
m=1 c m t
Proof of proposition 2 : Before explaining the proof, we will first introduce Ẑ: Monte-Carlo estimator of the unknown normalization constant Z(.). Since we are using K samples
Ẑ(x1:t, z Ait−1 1:t−1;K) = 1
K K∑ k=1 aθ,φ(δ i,k t |x1:t, z Ait−1 1:t−1). (15)
Algorithm 1 is producing an unbiased estimator of the marginal likelihood. We will first integrate out δ1:N,1:K1:T from the marginal likelihood estimator.
EQVSMC-PRC [ exp(L̂VSMC-PRC) ] =
∑ A1:N1:T−1 ∫ T∏ t=1 1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1)Ẑ(x1:t, z Ait−1 1:t−1;K) qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1)
QVSMC-PRC ( z1:N1:T , A 1:N 1:T−1, δ 1:N,1:K 1:T ) dz1:N1:T dδ 1:N,1:K 1:T
= ∑
A1:N1:T−1
∫ T∏ t=1 1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1) ∫ Ẑ(x1:t, z Ait−1 1:t−1;K) ∏K k=1 qφ(δ i,k t |x1:t, z Ait−1 1:t−1)dδ i,k t qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1)
QVSMC-PRC ( z1:N1:T , A 1:N 1:T−1 ) dz1:N1:T
= ∑
A1:N1:T−1
∫ T∏ t=1 1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1)Z(x1:t, z Ait−1 1:t−1) qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1) QVSMC-PRC ( z1:N1:T , A 1:N 1:T−1 ) dz1:N1:T
= ∑
A1:N1:T−1
∫ T∏ t=1 1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1) rθ,φ(zit|x1:t, z Ait−1 1:t−1) QVSMC-PRC ( z1:N1:T , A 1:N 1:T−1 ) dz1:N1:T
= ∑
A1:N1:T−1
∫ T∏ t=1 1 N N∑ i=1 wit ( N∏ i=1 rθ,φ(z i 1|x1) T−1∏ t=1 N∏ i=1 Discrete(Ait|αt)rθ,φ(zit+1|x1:t+1, z Ait 1:t) ) dz1:N1:T
= E [ T∏ t=1 1 N N∑ i=1 wit ] . (16)
Note that rθ,φ(zit|.) is the sampling density of PRC. It is easy to see that Eq. 16 is a standard SMC estimator for the marginal likelihood pθ(x1:T ). The proof is quite standard in SMC literature and can be found in (Naesseth et al., 2017; 2019). The key factor that makes our bound unbiased is the ability to produce unbiased ancestor samples despite the presence of intractable normalization constant. Using Jensen inequality we can easily show that VSMC-PRC bound is smaller then log marginal likelihood.
E [ L̂VSMC-PRC;K ] = ∫ T∑ t=1 log  1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1)Ẑ(x1:t, z Ait−1 1:t−1;K) qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1)  QVSMC-PRC ( z1:N1:T , A 1:N 1:T−1, δ 1:N,1:K 1:T ) dz1:N1:T dδ 1:N,1:K 1:T dA 1:N 1:T−1 ≤ log pθ(x1:T ) .
(17)
We will show that E[L̂VSMC-PRC;K] is non-decreasing with K. Let’s define a collection of subsets {{Ii,t}Ni=1}Tt=1 ⊂ {1, 2, . . . ,K} having elements {i1, i2, . . . , im} randomly drawn from the set {1, 2, . . . ,K} having length m ≤ K. We can easily show the following expectation:
Ẑ(x1:t, z Ait−1 1:t−1;K) = EIi,t={i1,i2,...,im} [ Ẑ(x1:t, z Ait−1 1:t−1;m) ] = 1
K K∑ k=1 aθ,φ(δ i,k t |x1:t, z Ait−1 1:t−1)
(18)
Substitute the values of above expectation in equation 17. Use Jensen’s inequality i.e E[logX] ≤ logE[X] to complete the proof .
E [ L̂VSMC-PRC;K ]
= E  T∑ t=1 logE{{Ii,t}Ni=1}Tt=1  1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1)Ẑ(x1:t, z Ait−1 1:t−1;m) qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1)  ≥ E
E{{Ii,t}Ni=1}Tt=1 T∑ t=1 log  1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1)Ẑ(x1:t, z Ait−1 1:t−1;m) qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1)  ≥ E [ L̂VSMC-PRC;m ]
Now we will see what happens when the limit K →∞. Using dominated convergence theorem we can write down the estimator as follows:
limK→+∞ E[L̂VSMC-PRC]
= lim K→+∞
E  T∑ t=1 log  1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1)Ẑ(x1:t, z Ait−1 1:t−1;K) qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1)  = E
 T∑ t=1 log  1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1) limK→+∞ Ẑ(x1:t, z Ait−1 1:t−1;K) qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1)  = E
 T∑ t=1 log  1 N N∑ i=1 pθ(xt, z i t|x1:t−1, z Ait−1 1:t−1)Z(x1:t, z Ait−1 1:t−1) qφ(zit|x1:t, z Ait−1 1:t−1)aθ,φ(z i t|x1:t, z Ait−1 1:t−1)  ≤ log pθ(x1:T )
Proof of proposition 3 : We will show that PRC step refines the learned distribution.
KL ( r(zt|z Ait−1 1:t−1, x1:t)||p(zt|z Ait−1 1:t−1, x1:t) )
= ∫ log r(zt|zAit−11:t−1, x1:t) p(zt|z Ait−1 1:t−1, x1:t)  r(zt|zAit−11:t−1, x1:t)dzt = ∫ log
q(zt|zAit−11:t−1, x1:t)a(zt|zAit−11:t−1, x1:t) p(zt|z Ait−1 1:t−1, x1:t)Z(x1:t, z Ait−1 1:t−1)  r(zt|zAit−11:t−1, x1:t)dzt ≤ ∫ log
q(zt|zAit−11:t−1, x1:t)a(zt|zAit−11:t−1, x1:t) p(zt|z Ait−1 1:t−1, x1:t)Z(x1:t, z Ait−1 1:t−1)  q(zt|zAit−11:t−1, x1:t)dzt ≤ KL ( q(zt|z Ait−1 1:t−1, x1:t)||p(zt|z Ait−1 1:t−1, x1:t) ) +
∫ log a(zt|zAit−11:t−1, x1:t) Z(x1:t, z Ait−1 1:t−1)  q(zt|zAit−11:t−1, x1:t)dzt ≤ KL ( q(zt|z Ait−1 1:t−1, x1:t)||p(zt|z Ait−1 1:t−1, x1:t)
) First, we will use the property of negatively correlated random variables. Note that the following two random variables
X = log
( r(zt|z Ait−1 1:t−1 ,x1:t)
p(zt|z Ai t−1
1:t−1 ,x1:t)
) and Y = a(zt|z Ait−1 1:t−1, x1:t),
are negatively correlated. We know that for negatively correlated variables following identity holds E[XY ] ≤ E[X]E[Y ]. Further we have used Jensen’s inequality to show that
E[log a(zt|z Ait−1 1:t−1, x1:t)] ≤ logZ(x1:t, z Ait−1 1:t−1)
Case 1: M(i, t− 1)→ 0 implies r(zt|z Ait−1 1:t−1, x1:t)→ q(zt|z Ait−1 1:t−1, x1:t)
In this situation all samples would be accepted. Hence, we can express the sampling distribution as:
r(zt|z Ait−1 1:t−1, x1:t) = q(z i t|x1:t, z Ait−1 1:t−1), (19)
Case 2: M(i, t− 1)→∞ implies r(zt|z Ait−1 1:t−1, x1:t)→ p(zt|z Ait−1 1:t−1, x1:t)
In this case, the acceptance probability would reduce to standard rejection sampling. Therefore, the sampling distribution would become equal to the true posterior.
r(zt|z Ait−1 1:t−1, x1:t) =
q(zit|x1:t, z Ait−1 1:t−1)
p(zit,xt|x1:t−1,z Ait−1 1:t−1 )
M(i,t−1)q(zit|x1:t,z Ai t−1
1:t−1 )∫ q(zit|x1:t, z Ait−1 1:t−1) p(zit,xt|x1:t−1,z Ai t−1 1:t−1 )
M(i,t−1)q(zit|x1:t,z Ai t−1 1:t−1 )
dzit
= p(zit|x1:t, z Ait−1 1:t−1)

heading: B GRADIENT ESTIMATION
text: In this section, we will derive the unbiased gradients for the Monte-Carlo estimate E[L̂VSMC-PRC]. Note that we can express the complete gradient into three core components (assuming q() is reparametriz-
able):
∇θ,φE[L̂VSMC-PRC] = grep + gPRC + gRSAMP grep = EQVSMC-PRC [ ∇θ,φL̂VSMC-PRC(θ, φ;x1:T ,K) ] ,
gPRC = EQVSMC-PRC L̂VSMC-PRC(θ, φ;x1:T ,K)∇θ,φ  N∑ i=1 T∑ t=1 log a(zit|x1:t, z Ait−1 1:t−1) Z(x1:t, z Ait−1 1:t−1)  , gRSAMP = EQVSMC-PRC [ L̂VSMC-PRC(θ, φ;x1:T ,K)∇θ,φ
N∑ i=1 T−1∑ t=1 Ait logαt ] Note that for both gPRC and gRSAMP unbiased score gradient estimates are not available due to intractability. Therefore, we have used Monte-Carlo samples to estimate the log density for figure 3.

heading: C EXPERIMENTAL SETUP
text: For the real data experiment, we train a VRNN (Chung et al., 2015) on the polyphonic music dataset. Polyphonic music comprises of four datasets: : Nottingham, JSB chorales, Musedata, and Piano-midi.de. Each dataset was divided into standard train, validation, and test datset.
The validation data was used to tune the learning rate: we picked the learning rate from the following set: {3× 10−4, 1× 10−4, 3× 10−5, 1× 10−5}, instead of optimizing for each method, we picked the learning rate at which FIVO (Maddison et al., 2017) validation performance is the best. Once the learning rate is decided, we ran every method for the same number of iterations to ensure uniformity. We use a single-layer LSTM for modeling the hidden state having dimension dh. For a length T sequence, the variational distribution and joint data likelihood are defined as follows
r(z1:T |x1:T ) = ∏T t=1 qt(zt|ht(zt−1, xt−1, ht−1), xt)at(zt|ht(zt−1, xt−1, ht−1), xt)∏T
t=1 Zt(ht(zt−1, xt−1, ht−1), xt)
Note that at(zt|.) is the acceptance probability for the PRC step and Zt(.) is the intractable normalization constant. Similarly, we can write down the joint data likelihood as
p(z1:T , x1:T ) = T∏ t=1 pt(zt|ht(zt−1, xt−1, ht−1), xt)gt(xt|ht(zt−1, xt−1, ht−1), zt)
The conditional distributions pt(zt|.) and qt(zt|.) are assumed to be factorized Gaussians, where dimension of the latent variable zt is dz . Note that the conditional densities are parametrized by fully connected neural networks with a single layer having size dh. The output distribution gt(xt|.) is modeled by a set of 88 iid Bernoulli variables. Please see Table 2 for more details regarding implementation details. The unknown weights and biases are initialized using a Xavier initialization. For setting up the optimization we used a batch size of 4 with adam optimizer. The unknown hyperparameter M is updated once every 50 epoches through Eq. 14 to save time.

Review summaries:
This paper describes an SMC algorithm to sample the posterior distribution of latent states $p_\theta(z_{1:T}|x_{1:T})$ in a latent variable models $p_\theta(x_{1:T}, z_{1:T})$. The authors consider a completely general setting (the authors assume Eq.(1) but clearly there is nothing to assume here, this the standard Bayes rule). It is well known that the vanilla SMC sampler is a good candidate for ELBO because it provides an unbiased estimator of the likelihood. But the authors prefer here to use a more sophisticated version of the  SMC algorithm, which features a partial rejection algorithm, which amounts to eliminate proposed particles which are "large enough" likelihood.