SUMO: UNBIASED ESTIMATION OF LOG MARGINAL PROBABILITY FOR LATENT VARIABLE MODELS
Year: 2020

Abstract:
Standard variational lower bounds used to train latent variable models produce biased estimates of most quantities of interest. We introduce an unbiased estimator of the log marginal likelihood and its gradients for latent variable models based on randomized truncation of infinite series. If parameterized by an encoder-decoder architecture, the parameters of the encoder can be optimized to minimize its variance of this estimator. We show that models trained using our estimator give better test-set likelihoods than a standard importance-sampling based approach for the same average computational cost. This estimator also allows use of latent variable models for tasks where unbiased estimators, rather than marginal likelihood lower bounds, are preferred, such as minimizing reverse KL divergences and estimating score functions.

heading: 1 INTRODUCTION
text: Latent variable models are powerful tools for constructing highly expressive data distributions and for understanding how high-dimensional observations might possess a simpler representation. Latent variable models are often framed as probabilistic graphical models, allowing these relationships to be expressed in terms of conditional independence. Mixture models, probabilistic principal component analysis (Tipping & Bishop, 1999), hidden Markov models, and latent Dirichlet allocation (Blei et al., 2003) are all examples of powerful latent variable models. More recently there has been a surge of interest in probabilistic latent variable models that incorporate flexible nonlinear likelihoods via deep neural networks (Kingma & Welling, 2014). These models can blend the advantages of highly structured probabilistic priors with the empirical successes of deep learning (Johnson et al., 2016; Luo et al., 2018). Moreover, these explicit latent variable models can often yield relatively interpretable representations, in which simple interpolation in the latent space can lead to semantically-meaningful changes in high-dimensional observations (e.g., Higgins et al. (2017)).
It can be challenging, however, to fit the parameters of a flexible latent variable model, since computing the marginal likelihood of the data requires integrating out the latent variables. Typical approaches to this problem include the celebrated expectation maximization algorithm (Dempster et al., 1977), Markov chain Monte Carlo, and the Laplace approximation. Variational inference generalizes expectation maximization by forming a lower bound on the aforementioned (log) marginal likelihood, using a tractable approximation to the unmanageable posterior over latent variables. The maximization of this lower bound—rather than the true log marginal likelihood—is often relatively straightforward when using automatic differentiation and Monte Carlo sampling. However, a
∗Equal contribution.
lower bound may be ill-suited for tasks such as posterior inference and other situations where there exists an entropy maximization objective; for example in entropy-regularized reinforcement learning (Williams & Peng, 1991; Mnih et al., 2016; Norouzi et al., 2016) which requires minimizing the log probability of the samples under the model.
While there is a long history in Bayesian statistics of estimating the marginal likelihood (e.g., Newton & Raftery (1994); Neal (2001)), we often want high-quality estimates of the logarithm of the marginal likelihood, which is better behaved when the data is high dimensional; it is not as susceptible to underflow and it has gradients that are numerically sensible. However, the log transformation introduces some challenges: Monte Carlo estimation techniques such as importance sampling do not straightforwardly give unbiased estimates of this quantity. Nevertheless, there has been significant work to construct estimators of the log marginal likelihood in which it is possible to explicitly trade off between bias against computational cost (Burda et al., 2016; Bamler et al., 2017; Nowozin, 2018). Unfortunately, while there are asymptotic regimes where the bias of these estimators approaches zero, it is always possible to optimize the parameters to increase this bias to infinity.
In this work, we construct an unbiased estimator of the log marginal likelihood. Although there is no theoretical guarantee that this estimator has finite variance, we find that it can work well in practice. We show that this unbiased estimator can train latent variable models to achieve higher test log-likelihood than lower bound estimators at the same expected compute cost. More importantly, this unbiased estimator allows us to apply latent variable models in situations where these models were previously problematic to optimize with lower bound estimators. Such applications include latent variable modeling for posterior inference and for reinforcement learning in high-dimensional action spaces, where an ideal model is one that is highly expressive yet efficient to sample from.

heading: 2 PRELIMINARIES

heading: 2.1 LATENT VARIABLE MODELS
text: Latent variable models (LVMs) describe a distribution over data in terms of a mixture over unobserved quantities. Let pθ(x) be a family of probability density (mass) functions on a data space X , indexed by parameters θ. We will generally refer to this as a “density” for consistency, even when the data should be understood to be discrete; similarly we will use integrals even when the marginalization is over a discrete set. In a latent variable model, pθ(x) is defined via a space of latent variables Z , a family of mixing measures on this latent space with density denoted pθ(z), and a conditional distribution pθ(x | z). This conditional distribution is sometimes called an “observation model” or a conditional likelihood. We will take θ to parameterize both pθ(x | z) and pθ(z) in the service of determining the marginal pθ(x) via the mixture integral:
pθ(x) := ∫ Z pθ(x | z)pθ(z) dz = Ez∼pθ(z) [pθ(x | z)] . (1)
This simple formalism allows for a large range of modeling approaches, in which complexity can be baked into the latent variables (as in traditional graphical models), into the conditional likelihood (as in variational autoencoders), or into both (as in structured VAEs). The downside of this mixing approach is that the integral may be intractable to compute, making it difficult to evaluate pθ(x)— a quantity often referred to in Bayesian statistics and machine learning as the marginal likelihood or evidence. Various Monte Carlo techniques have been developed to provide consistent and often unbiased estimators of pθ(x), but it is usually preferable to work with log pθ(x) and unbiased estimation of this quantity has, to our knowledge, not been previously studied.

heading: 2.2 TRAINING LATENT VARIABLE MODELS
text: Fitting a parametric distribution to observed data is often framed as the minimization of a difference between the model distribution and the empirical distribution. The most common difference measure is the forward Kullback-Leibler (KL) divergence; if pdata(x) is the empirical distribution and pθ(x) is a parametric family, then minimizing the KL divergence (DKL) with respect to θ is equivalent to maximizing the likelihood:
DKL(pdata || pθ) = ∫ X pdata(x) log pdata(x) pθ(x) dx = −Edata [log pθ(x)] + const . (2)
Equivalently, the optimization problem of finding the MLE parameters θ comes down to maximizing the expected log probability of the data:
θMLE = arg min θ DKL(pdata || pθ) = arg max θ Edata [log pθ(x)] . (3)
Since expectations can be estimated in an unbiased manner using Monte Carlo procedures, simple subsampling of the data enables powerful stochastic optimization techniques, with stochastic gradient descent in particular forming the basis for learning the parameters of many nonlinear models. However, this requires unbiased estimates of ∇θ log pθ(x), which are not available for latent variable models. Instead, a stochastic lower bound of log pθ(x) is often used and then differentiated for optimization.
Though many lower bound estimators (Burda et al., 2016; Bamler et al., 2017; Nowozin, 2018) are applicable, we focus on an importance-weighted evidence lower bound (Burda et al., 2016). This lower bound is constructed by introducing a proposal distribution q(z;x) and using it to form an importance sampling estimate of the marginal likelihood:
pθ(x) = ∫ Z pθ(x | z) pθ(z) dz = ∫ Z q(z;x) pθ(x | z) pθ(z) q(z;x) dz = Ez∼q [ pθ(x | z) pθ(z) q(z;x) ] . (4)
IfK samples are drawn from q(z;x) then this provides an unbiased estimate of pθ(x) and the biased “importance-weighted autoencoder” estimator IWAEK(x) of log pθ(x) is given by
IWAEK(x) := log 1
K K∑ k=1 pθ(x | zk) pθ(zk) q(zk;x) , zk iid∼ q(z;x) . (5)
The special case of K = 1 generates an unbiased estimate of the evidence lower bound (ELBO), which is often used for performing variational inference by stochastic gradient descent. While the IWAE lower bound acts as a useful replacement of log pθ(x) in maximum likelihood training, it may not be suitable for other objectives such as those that involve entropy maximization. We discuss tasks for which a lower bound estimator would be ill-suited in Section 3.4.
There are two properties of IWAE that will allow us to modify it to produce an unbiased estimator: First, it is consistent in the sense that as the number of samples K increases, the expectation of IWAEK(x) converges to log pθ(x). Second, it is also monotonically non-decreasing in expectation:
log pθ(x) = lim K→∞
E[IWAEK(x)] and E[IWAEK+1(x)] ≥ E[IWAEK(x)] . (6)
These properties are sufficient to create an unbiased estimator using the Russian roulette estimator.

heading: 2.3 RUSSIAN ROULETTE ESTIMATOR
text: In order to create an unbiased estimator of the log probability function, we employ the Russian roulette estimator (Kahn, 1955). This estimator is used to estimate the sum of infinite series, where each sample of the estimator almost surely requires only a finite amount of computation. Intuitively, the Russian roulette estimator relies on a randomized truncation and upweighting of each term to account for the possibility of not computing the remaining terms.
To illustrate the idea, let ∆̃k denote the k-th term of an infinite series. Assume the partial sum of the series ∑∞ k=1 ∆̃k converges to some quantity we wish to obtain. We can construct a simple estimator by always computing the first term then flipping a coin b ∼ Bernoulli(q) to determine whether we stop or continue evaluating the remaining terms. With probability 1− q, we compute the rest of the series. By reweighting the remaining future terms by 1/(1−q), we obtain an unbiased estimator:
Ỹ = ∆̃1 + (∑∞ k=2 ∆̃k 1− q ) 1b=0 + (0)1b=1 E[Ỹ ] = ∆̃1 + ∑∞ k=2 ∆̃k 1− q (1− q) = ∞∑ k=1 ∆̃k.
To obtain the “Russian roulette” (RR) estimator (Forsythe & Leibler, 1950), we repeatedly apply this trick to the remaining terms. In effect, we make the number of terms a random variable K, taking values in 1, 2, . . . to use in the summation (i.e., the number of coin flips) from some distribution
with probability mass function p(K) = P(K = K) with support over the positive integers. With K drawn from p(K), the estimator takes the form:
Ŷ (K) = K∑ k=1 ∆̃k P(K ≥ k)
EK∼p(K)[Ŷ (K)] = ∞∑ k=1 ∆̃k . (7)
The equality on the right hand of equation 7 holds so long as (i) P(K ≥ k) > 0, ∀k > 0, and (ii) the series converges absolutely, i.e., ∑∞ k=1 |∆̃k| <∞ (Chen et al. (2019); Lemma 3). This condition ensures that the average of multiple samples will converge to the value of the infinite series by the law of large numbers. However, the variance of this estimator depends on the choice of p(K) and can potentially be very large or even infinite (McLeish, 2011; Rhee & Glynn, 2015; Beatson & Adams, 2019).

heading: 3 SUMO: UNBIASED ESTIMATION OF LOG PROBABILITY FOR LVMS

heading: 3.1 RUSSIAN ROULETTE TO DEBIAS LOWER BOUNDS
text: We can turn any absolutely convergent series into a telescoping series and apply the Russian roulette randomization to form an unbiased stochastic estimator. We focus here on the IWAE bound described in Section 2.2. Let ∆k(x) = IWAEk+1(x)− IWAEk(x), then since Eq[∆k(x)] converges absolutely, we apply equation 7 to construct our estimator, which we call SUMO (Stochastically Unbiased Marginalization Objective). The detailed derivation of SUMO is in Appendix A.1.
SUMO(x) = IWAE1(x) + K∑ k=1 ∆k(x) P(K ≥ k) where K ∼ p(K) . (8)
The randomized truncation of the series using the Russian roulette estimator means that this is an unbiased estimator of the log marginal likelihood, regardless of the distribution p(K):
E [SUMO(x)] = log pθ(x) , (9)
where the expectation is taken over p(K) and q(z;x) (see Algorithm 1 for our exact sampling procedure). Furthermore, under some conditions, we have E [∇θSUMO(x)] = ∇θE [SUMO(x)] = ∇θ log pθ(x) (see Appendix A.4).
3.2 OPTIMIZING VARIANCE-COMPUTE PRODUCT BY CHOICE OF p(K)
To efficiently optimize a limit, one should choose an estimator to minimize the product of the second moment of the gradient estimates and the expected compute cost per evaluation. The choice of p(K) effects both the variance and computation cost of our estimator. Denoting Ĝ := ∇θŶ and ∆gk := ∇θ[IWAEk+1(x)− IWAEk(x)], the Russian roulette estimator is optimal across a broad family of unbiased randomized truncation estimators if the ∆gk are statistically independent, in which case it has second moment E||Ĝ||22 = ∑∞ k=1 E||∆gk|| 2 2/P(K≥k) (Beatson & Adams, 2019). While the ∆gk are not in fact strictly independent with our sampling procedure (Algorithm 1), and other estimators within the family may perform better, we justify our choice by showing that E∆i∆j for i 6= j converges to zero much faster than E∆2k (Appendices A.2 & A.3). In the following, we assume independence of ∆gk and choose p(K) to minimize the product of compute and variance.
Algorithm 1 Computing SUMO, an unbiased estimator of log p(x).
Input: x, m ≥ 1, encoder q(z;x), decoder p(x, z), p(K), reverse cdf(·) = P(K ≥ ·) 1: Sample K ∼ p(K) 2: Sample {zk}K+mk=1
iid∼ q(z;x) 3: logwk ← log p(x, zk)− log q(zk;x) 4: ks← [1, . . . ,K +m] 5: cum iwae← log_cumsum_exp(logwk)− log(ks[:k+1]) 6: inv weights = 1/reverse cdf(ks)
return cum iwae[m-1] + sum(inv weights * (cum iwae[m:] - cum iwae[m-1:-1]))
We first show that E||∆gk||22 is O(1/k2) (Appendix A.5). This implies the optimal compute-variance product (Rhee & Glynn, 2015; Beatson & Adams, 2019) is given by P(K ≥ k) ∝ √ E||∆gk||22). In our case, this gives P(K ≥ k) = 1/k, which results in an estimator with infinite expected computation and no finite bound on variance. In fact, any p(K) which gives rise to provably finite variance requires a heavier tail than P(K ≥ k) = 1/k and so will have infinite expected computation. Though we could not theoretically show that our estimator and gradients have finite variance, we empirically find that gradient descent converges —even in the setting of minimizing log probability. We plot ||∆k||22 for the toy variational inference task used to assess signal to noise ratio in Tucker et al. (2018) and Rainforth et al. (2018b), and find that they converge faster than 1k2 in practice (Appendix A.6). While this indicates the variance is better than the theoretical bound, an estimator having infinite expected computation cost will always be an issue as it indicates significant probability of sampling arbitrarily large K. We therefore modify the tail of the sampling distribution such that the estimator has finite expected computation:
P(K ≥ k) = {
1/k if k < α 1/α · (1− 0.1)k−α if k ≥ α (10)
We typically choose α = 80, which gives an expected computation cost of approximately 5 terms.

heading: 3.2.1 TRADING VARIANCE AND COMPUTE
text: One way to improve the RR estimator is to construct it so that some minimum number of terms (denoted here as m) are always computed. This puts a lower bound on the computational cost, but can potentially lower variance, providing a design space for trading off estimator quality against computational cost. This corresponds to a choice of RR estimator in which P(K = K) = 0 for K ≤ m. This computes the sum out to m terms (effectively computing IWAEm) and then estimates the remaining difference with Russian roulette:
SUMO(x) = IWAEm(x) + K∑ k=m ∆k(x) P(K ≥ k) , K ∼ p(K) (11)
In practice, instead of tuning parameters of p(K), we setm to achieve a given expected computation cost per estimator evaluation for fair comparison with IWAE and related estimators.
3.3 TRAINING q(z;x) TO REDUCE VARIANCE
The SUMO estimator does not require amortized variational inference, but the use of an “encoder” to produce an approximate posterior q(z;x) has been shown to be a highly effective way to perform rapid feedforward inference in neural latent variable models. We use φ to denote the parameters of the encoder qφ(z;x). However, the gradients of SUMO with respect to φ are in expectation zero precisely because SUMO is an unbiased estimator of log pθ(x), regardless of our choice of qφ(z;x). Nevertheless, we would expect the choice of qφ(z;x) significantly impacts the variance of our estimator. As such, we optimize qφ(z;x) to reduce the variance of the SUMO estimator. We can obtain unbiased gradients in the following way (Ruiz et al., 2016; Tucker et al., 2017):
∇φVar[SUMO] = ∇φE[SUMO2]−((((( ((∇φ(E[SUMO])2 = E[∇φSUMO2] . (12)
Notably, the expectation of this estimator depends on the variance of SUMO, which we have not been able to bound. In practice, we observe gradients which are sometimes very large. We apply gradient clipping to the encoder to clip gradients which are excessively large in magnitude. This helps stabilize the training progress but introduces bias into the encoder gradients. Fortunately, the encoder itself is merely a tool for variance reduction, and biased gradients with respect to the encoder can still significantly help optimization.

heading: 3.4 APPLICATIONS OF UNBIASED LOG PROBABILITY
text: Here we list some applications where an unbiased log probability is useful. Using SUMO to replace existing lower bound estimates allows latent variable models to be used for new applications where a lower bound is inappropriate. As latent variable models can be both expressive and efficient to sample from, they are frequently useful in applications where the data is high-dimensional and samples from the model are needed.
Minimizing log pθ(x). Some machine learning objectives include terms that seek to increase the entropy of the learned model. The “reverse KL” objective—often used for training models to perform approximate posterior inference—minimizes Ex∼pθ(x)[log pθ(x)− log π(x)] where π(x) is a target density that may only be known up a normalization constant. Local updates of this form are the basis of the expectation propagation procedure (Minka, 2001). This objective has also been used for distilling autoregressive models that are inefficient at sampling (Oord et al., 2018). Moreover, reverse KL is connected to the use of entropy-regularized objectives (Williams & Peng, 1991; Ziebart, 2010; Mnih et al., 2016; Norouzi et al., 2016) in decision-making problems, where the goal is to encourage the decision maker toward exploration and prevent it from settling into a local minimum.
Unbiased score function ∇θ log pθ(x). The score function is the gradient of the log-likelihood with respect to the parameters and has uses in estimating the Fisher information matrix and performing stochastic gradient Langevin dynamics (Welling & Teh, 2011), among other applications. Of particular note, the REINFORCE gradient estimator (Williams, 1992)—generally applicable for optimizing objectives of the form maxθ Ex∼pθ(x)[R(x)]—is estimated using the score function. This can be replaced with the gradient of SUMO which itself is an estimator of the score function∇θ log pθ(x).
∇θEx∼pθ(x)[R(x)] = Ex∼pθ(x)[R(x)∇θ log pθ(x)] = Ex∼pθ(x)[R(x)∇θE[SUMO(x)]] = Ex∼pθ(x)[E[R(x)∇θSUMO(x)]]
(13)
where the inner expectation is over the stochasticity of the SUMO estimator. Such estimators are often used for reward maximization in reinforcement learning where pθ(x) is a stochastic policy.

heading: 4 RELATED WORK
text: There is a long history in Bayesian statistics of marginal likelihood estimation in the service of model selection. The harmonic mean estimator (Newton & Raftery, 1994), for example, has a long (and notorious) history as a consistent estimator of the marginal likelihood that may have infinite variance (Murray & Salakhutdinov, 2009) and exhibits simulation psuedo-bias (Lenk, 2009). The Chib estimator (Chib, 1995), the Laplace approximation, and nested sampling (Skilling, 2006) are alternative proposals that can often have better properties (Murray & Salakhutdinov, 2009). Annealed importance sampling (Neal, 2001) probably represents the gold standard for marginal likelihood estimation. These, however, turn into consistent estimators at best when estimating the log marginal probability (Rainforth et al., 2018a). Bias removal schemes such as jackknife variational inference (Nowozin, 2018) have been proposed to debias log-evidence estimation, IWAE in particular. Hierarchical IWAE (Huang et al., 2019) uses a joint proposal to induce negative correlation among samples and connects the convergence of variance of the estimator and the convergence of the lower bound.
Russian roulette also has a long history. It dates back to unpublished work from von Neumann and Ulam, who used it to debias Monte Carlo methods for matrix inversion (Forsythe & Leibler, 1950) and particle transport problems (Kahn, 1955). It has gained popularity in statistical physics (Spanier & Gelbard, 1969; Kuti, 1982; Wagner, 1987), for unbiased ray tracing in graphics and rendering (Arvo & Kirk, 1990), and for a number of estimation problems in the statistics community (Wei & Murray, 2017; Lyne et al., 2015; Rychlik, 1990; 1995; Jacob & Thiery, 2015; Jacob et al., 2017). It has also been independently rediscovered many times (Fearnhead et al., 2008; McLeish, 2011; Rhee & Glynn, 2012; Tallec & Ollivier, 2017).
The use of Russian roulette estimation in deep learning and generative modeling applications has been gaining traction in recent years. It has been used to solve short-term bias in optimization problems (Tallec & Ollivier, 2017; Beatson & Adams, 2019). Wei & Murray (2017) estimates the reciprocal normalization constant of an unnormalized density. Han et al. (2018) uses a similar random truncation approach to estimate the distribution of eigenvalues in a symmetric matrix. Along similar motivations with our work, Chen et al. (2019) uses this estimator to construct an unbiased estimator of the change of variables equation in the context of normalizing flows (Rezende & Mohamed, 2015), and Xu et al. (2019) uses it to construct unbiased log probability for a nonparameteric distribution in the context of variational autoencoders (Kingma & Welling, 2014).
Though we extend latent variable models to applications that require unbiased estimates of log probability and benefit from efficient sampling, an interesting family of models already fulfill these requirements. Normalizing flows (Rezende & Mohamed, 2015; Dinh et al., 2017) offer exact log probability and certain models have been proven to be universal density estimators (e.g. Huang et al. (2018)). However, these models often require restrictive architectural choices with no dimensionality-reduction capabilities, and make use of many more parameters to scale up (Kingma & Dhariwal, 2018) than alternative generative models. Discrete variable versions of these models are still in their infancy and make use of biased gradients (Tran et al., 2019; Hoogeboom et al., 2019), whereas latent variable models naturally extend to discrete observations.

heading: 5 DENSITY MODELING EXPERIMENTS
text: We first compare the performance of SUMO when used as a replacement to IWAE with the same expected cost on density modeling tasks. We make use of two benchmark datasets: dynamically binarized MNIST (LeCun et al., 1998) and binarized OMNIGLOT (Lake et al., 2015).
We use the same neural network architecture as IWAE (Burda et al., 2016). The prior p(z) is a 50- dimensional standard Gaussian distribution. The conditional distributions p(xi|z) are independent Bernoulli, with the decoder parameterized by two hidden layers, each with 200 tanh units. The approximate posterior q(z;x) is also a 50-dimensional Gaussian distribution with diagonal covariance, whose mean and variance are both parameterized by two hidden layers with 200 tanh units. We reimplemented and tuned IWAE, obtaining strong baseline results which are better than those previously reported. We then used the same hyperparameters to train with the SUMO estimator. We find clipping very large gradients can help performance, as large gradients may be infrequently sampled. This introduces a small amount of bias into the gradients while reducing variance, but can nevertheless help achieve faster convergence and should still result in a less-biased estimator. A posthoc study of the effect on final test performance as a function of this bias-variance tradeoff mechanism is discussed in Appendix A.7. We note that gradient clipping is only done for the density modeling experiments.
The averaged test log-likelihoods and standard deviations over 3 runs are summarized in Table 1. To be consistent with existing literature, we evaluate our model using IWAE with 5000 samples. In all the cases, SUMO achieves slightly better performance than IWAE with the same expected cost. We also bold the results that are statistically insignificant from the best performing model according to an unpaired t-test with significance level 0.05. However, we do see diminishing returns as we increase k, suggesting that as we increase compute, the variance of our estimator may impact performance more than the bias of IWAE.

heading: 6 LATENT VARIABLES MODELS FOR ENTROPY MAXIMIZATION
text: We move on to our first task for which a lower bound estimate of log probability would not suffice. The reverse KL objective is useful when we have access to a (possibly unnormalized) target
distribution but no efficient sampling algorithm.
min θ DKL(pθ(x)||p∗(x)) = min θ Ex∼pθ(x)[log pθ(x)− log p ∗(x)] (14)
A major problem with fitting latent variables models to this objective is the presence of an entropy maximization term, effectively a minimization of log pθ(x). Estimating this log marginal probability with a lower bound estimator could result in optimizing θ to maximize the bias of the estimator instead of the true objective. Our experiments demonstrate that this causes IWAE to often fail to optimize the objective unless we use a large amount of computation.
Modifying IWAE. The bias of the IWAE estimator can be interpreted as the KL between an importance-weighted approximate posterior qIW (z;x) implicitly defined by the encoder and the true posterior p(z|x) (Domke & Sheldon, 2018). Both the encoder and decoder parameters can therefore affect this bias. In practice, we find that the encoder optimization proceeds at a faster timescale than the decoder optimization: i.e., the encoder can match qIW (z;x) to the decoder’s posterior p(z|x) more quickly than the latter can match an objective. For this reason, we train the encoder to reduce bias and use a minimax training objective
min p(x,z) max q(z;x)
Ex∼p(x)[IWAEK(x)− log p∗(x)] (15)
Though this is still a lower bound with unbounded bias, it makes for a stronger baseline than optimizing q(z;x) in the same direction as p(x, z). We find that this approach can indeed work in practice, but requires setting k extremely high.
We choose a “funnel” target distribution (Figure 1) similar to the distribution used as a benchmark for inference in Neal et al. (2003), where p∗ has support in R2 and is defined p∗(x1, x2) = N (x1; 0, 1.352)N (x2; 0, e2x1) We use neural networks with one hidden layer of 200 hidden units and tanh activations for both the encoder and decoder networks. We use 20 latent variables, with p(z), pθ(x|z), and qφ(z;x) all being Gaussian distributed.
Figure 2 shows the learning curves when using IWAE and SUMO. Unless k is set very large, IWAE will at some point start optimizing the bias instead of the actual objective. The reverse KL is a non-negative quantity, so any estimate significantly below zero can be attributed to the unbounded bias. On the other hand, SUMO
correctly optimizes for the objective even with a small expected cost. Increasing the expected cost k for SUMO simply reduces variance. We also found that if k is set sufficiently large, then IWAE can work when we train using the minimax objective in equation 15, suggesting that a sufficiently debiased estimator can also work in practice. However, this requires much more compute and likely does not scale compared to SUMO. We also visualize the contours of the resulting models in Figure 1. For IWAE, we visualize the model a few iterations before it reaches numerical instability.

heading: 7 LATENT VARIABLE POLICIES FOR COMBINATORIAL OPTIMIZATION
text: Let us now consider the problem of finding the maximum of a non-differentiable function, a special case of reinforcement learning without an interacting environment. Variational optimization (Staines & Barber, 2012) can be used to reformulate this as the optimization of a parametric distribution,
max x R(x) ≥ max θ Ex∼pθ(x)[R(x)], (16)
which is now a differentiable function with respect to the parameters θ, whose gradients can be estimated using a combination of the REINFORCE gradient estimator and the SUMO estimator (equation 13). Furthermore, entropy regularized reinforcement learning—where we maximize R(x) + λH(pθ) with H(pθ) being the entropy of pθ(x)—encourages exploration and is inherently related to minimizing a reverse KL objective with the target being an exponentiated reward (Norouzi et al., 2016).
For concreteness, we focus on the problem of quadratic pseudo-Boolean optimization (QPBO) where the objective is to maximize
R(x) = ∑ i=1 wi(xi) + ∑ i<j wij(xi, xj) (17)
where {xi}di=1 ∈ {0, 1} are binary variables. Without further assumptions, QPBO is NPhard (Boros & Hammer, 2002). As there exist complex dependencies between the binary variables and optimization of equation 16 requires sampling from the policy distribution pθ(x), a model that is both expressive and allows efficient sampling would be ideal. For this reason, we motivate the use of latent variable models with independent conditional distributions, trained using the SUMO estimator. Our baselines are an autoregressive policy, which captures dependencies but for which sampling must be performed sequentially, and an independent policy, which is easy to sample from but captures no dependencies.
pLVM(x) := ∫ d∏ i=1 pθ(xi|z)p(z)dz pAutoreg(x) := d∏ i=1 p(xi|x<i) pIndep(x) := d∏ i=1 p(xi)
We note that Haarnoja et al. (2018) also argued for latent variable policies in favor of learning diverse strategies but ultimately made use of normalizing flows which do not require marginalization.
We constructed one problem instance for each d ∈ {100, 500}, which we note are already intractable for exact optimization. For each instance, we randomly sampled the weights wi and wij uniformly from the interval [−1, 1]. Figure 3 shows the performance of each policy model. In general, the independent policy is quick to converge to a local minima and is unable to explore different regions, whereas more complex models have a better grasp of the “frontier” of reward distributions during optimization. The autoregressive model works well overall but is much slower to train due to its sequential sampling procedure; with d = 500, it is 19.2× slower per iteration than SUMO. Surprisingly, we find that estimating the REINFORCE gradient with IWAE results in decent performance when no entropy regularization is present. With entropy regularization, all policies improve
significantly; however, training with IWAE in this setting results in performance similar to the independent model, suggesting that it may not be using the latent state. On the other hand, SUMO works with both REINFORCE gradient estimation and entropy regularization, albeit at the cost of slower convergence due to variance.

heading: 8 CONCLUSION
text: We introduced SUMO, a new unbiased estimator of the log probability for latent variable models, and demonstrated tasks for which this estimator performs better than standard lower bounds. Specifically, we investigated applications involving entropy maximization where a lower bound performs poorly, but our unbiased estimator can train properly with relatively smaller amount of compute.
In the future, we plan to investigate new families of gradient-based optimizers which can handle heavy-tailed stochastic gradients. It may also be fruitful to investigate the use of convex combination of consistent estimators within the SUMO approach, as any convex combination is unbiased, or to apply variance reduction methods to increase stability of training with SUMO.

heading: ACKNOWLEDGEMENTS
text: This work was partially funded by NSF IIS-1421780. Y.L and J.Z were supported by the NSF China Project (No. 61620106010), Beijing NSF Project (No. L172037), the JP Morgan Faculty Research Program and the NVIDIA NVAIL Program with GPU/DGX Acceleration.

heading: A APPENDIX
text: A.1 DERIVATION OF SUMO
Let
Eq[IWAEk(x)] = Ez1,...,zk∼q(z;x)
[ log 1
K K∑ k=1 pθ(x | zk) pθ(zk) q(zk;x)
] ,
where z1, .., zk are sampled independently from q(z;x). And we define the k-th term of the infinite series ∆̃k(x) := Eq[IWAEk+1(x)] − Eq[IWAEk(x)]. Using the properties of IWAE in equation 6, we have ∆̃k(x) ≥ 0, and
∞∑ k=1 |∆̃k(x)| = ∞∑ k=1 ∆̃k(x) = lim k→∞ Eq[IWAEk(x)]− Eq[IWAE1(x)]
= log pθ(x)− Eq[IWAE1(x)] <∞, (18)
which means the series converges absolutely. This is a sufficient condition for finite expectation of the Russian roulette estimator (Chen et al. (2019); Lemma 3). Applying equation 7 to the series:
log pθ(x) = Eq[IWAE1(x)] + ∞∑ k=1 ∆̃k(x) (19)
= Eq[IWAE1(x)] + EK∼p(K) [ K∑ k=1 ∆̃k(x) P(K ≥ k) ] (20)
= EK∼p(K) [ Eq[IWAE1(x)] +
K∑ k=1 Eq [IWAEk+1(x)− IWAEk(x)] P(K ≥ k)
] (21)
= EK∼p(K) [ Eq [ IWAE1(x) +
K∑ k=1 IWAEk+1(x)− IWAEk(x) P(K ≥ k)
]] . (22)
Let ∆k(x) := IWAEk+1(x)− IWAEk(x), Hence our estimator is constructed:
SUMO(x) = IWAE1(x) + K∑ k=1 ∆k(x) P(K ≥ k) , K ∼ p(K), zk iid∼ q(z;x) . (23)
And it can be easily seen from equation 22 and equation 23 that SUMO is an unbiased estimator of the log marginal likelihood:
EK∼p(K),z1,...,zK∼q(z;x) [SUMO(x)] = log pθ(x). (24)
A.2 CONVERGENCE OF ∆k
We follow the analysis of JVI (Nowozin, 2018), which applied the delta method for moments to show the asymptotic results on the bias and variance of IWAEk both at a rate of O( 1k ). We build on this analysis to analyze the convergence of ∆k.
Let wi = p(x|zi)p(zi) q(zi;x) and we define Yk := 1k ∑k i=1 wi as the sample mean and we have E[Yk] = E[w] = µ.
IWAEk = log Yk = log [µ+ (Yk − µ)]
= logµ− ∞∑ t=1 (−1)t tµt (Yk − µ)t
(25)
We note that we rely on ||Yk−µ|| < 1 for this power series to converge. This condition was implicitly assumed, but not explicitly noted, in (Nowozin, 2018). This condition will hold for sufficiently
large k so long as the moments ofwi exist: one could bound the probability ||Yk−µ|| ≥ 1 by Chebyshev’s inequality or by the Central Limit Theorem. We use the central moments γt := E[(Yk − µ)t] and µt := E[(w − µ)t] for t ≥ 2.
E∆2k = E(IWAEk+1 − IWAEk)2 (26)
= E [ logµ−
∞∑ t=1 (−1)t tµt (Yk+1 − µ)t − logµ+ ∞∑ t=1 (−1)t tµt (Yk − µ)t
]2 (27)
= E [ ∞∑ t=1 (−1)t tµt [ (Yk − µ)t − (Yk+1 − µ)t ]]2 (28)
Expanding Eq. 28 to order two gives E∆2k = E [ − 1 µ (Yk − µ− Yk+1 + µ) + 1 2µ2 [ (Yk − µ)2 − (Yk+1 − u)2 ]]2 + o(k−2) (29)
= 1 µ2 E [ Yk+1 − Yk + 1 2µ (Yk + Yk+1 − 2µ)(Yk − Yk+1) ]2 + o(k−2) (30)
= 1 µ2 E [ 2(Yk+1 − Yk) + 1 2µ (Yk + Yk+1)(Yk − Yk+1) ]2 + o(k−2) (31)
Since we use cumulative sum to compute Yk and Yk+1, we obtain Yk+1 = kYk+wk+1
k+1 .
=⇒ E∆2k = 1
µ2 E [ 2 wk+1 − 1 k + 1 + (wk+1 + 2k+1k+1 ∑ki=1 wk 2kµ )(wk+1 − 1 k + 1 )]2 + o(k−2) (32)
We note that wk+1−1k+1 = O( 1 k ) and
wk+1+ 2k+1 k+1 ∑k i=1 wk
2kµ = O(1). Therefore ∆k isO( 1 k ), and E∆ 2 k =
O( 1k2 ).
A.3 CONVERGENCE OF ∆k∆j
Without loss of generality, suppose j ≥ k + 1,
E∆k∆j=E [( ∞∑ t=1 (−1)t tµt [ (Yk − µ)t − (Yk+1 − µ)t ])( ∞∑ t=1 (−1)t tµt [ (Yj − µ)t − (Yj+1 − µ)t ])] (33)
For clarity, let Ck = Yk − µ be the zero-mean random variable. Nowozin (2018) gives the relations
E[C2k ] = γ2 = µ2 k
(34)
E[C3k ] = γ3 = µ3 k2
(35)
E[C4k ] = γ4 = 3µ22 k2 + µ4 − 3µ22 k3 (36)
E∆k∆j = E [( ∞∑ t=1 (−1)t tµt (Ctk − Ctk+1) )( ∞∑ t=1 (−1)t tµt (Ctj − Ctj+1) )] (37)
Expanding both the sums inside the brackets to order two:
E∆k∆j ≈ E 1
µ2 (Ck+1 − Ck)(Cj+1 − Cj) (1)
− E 1 2µ3 (C2k+1 − C2k)(Cj+1 − Cj) (2)
− E 1 2µ3 (Ck+1 − Ck)(C2j+1 − C2j ) (3)
+ E 1
4µ4 (C2k+1 − C2k)(C2j+1 − C2j ) (4)
We will proceed by bounding each of the terms (1), (2), (3), (4). First, we decompose Cj . Let Bk,j := 1 j ∑j i=k+1(wi − µ).
Cj = 1
j
( kCk +
j∑ i=k+1 (wi − µ)
) = k
j Ck +
1
j j∑ i=k+1 (wi − µ) = k j Ck +Bk,j (38)
We know that Bk,j is independent of Ck and E[Bk,j ] = 0, implying E[CkBk,j ] = 0. Note C2j = k2 j2C 2 k + 2 k jCkBk,j +B 2 k,j .
Now we show that (1) is zero:
E[ 1
µ2 (Ck+1 − Ck)(Cj+1 − Cj)] =
1 µ2 E [ Ck+1 k + 1 j + 1 + Ck+1Bj+1,k+1
− k + 1 j C2k+1 −Bj,k+1Ck+1 − Ck k j + 1
− CkBj+1,k + k
j C2k + CkBj,k ] = 1
µ2 E[ − k + 1 j(j + 1) C2k+1 + C 2 k
k
j(j + 1) ]
= 1 µ2 [− k + 1 j(j + 1) µ2 k + 1 + µ2 k
k
j(j + 1) ] = 0
We now investigate (2):
E[− 1 2µ3 (C2k+1 − C2k)(Cj+1 − Cj)] = 1 2µ3 E [ C3k k j + 1 + C2kBj+1,k − C3k k j − C2kBj,k
+ C3k+1 k + 1
j + 1 + C2k+1Bj,k +
k j C2k + CkBj+1,k ] = 1
µ2 E[− k + 1 j(j + 1) C2k+1 + C 2 k
k
j(j + 1) ]
= 1 2µ3 [− µ3 kj(j + 1) + µ3 (k + 1)j(j + 1) ] = − µ3 2µ3 [ 1 k(k + 1)j(j + 1) ]
We now show that (3) is zero:
E[ 1
2µ3 (Ck+1 − Ck)(C2j − C2j+1)] =
1
2µ3 E[Ck+1C2j − Ck+1C2j+1 − CkC2j + CkC2j+1)]
= 1 2µ3 [ µ3 j2 − µ3 (j + 1)2 − µ3 j2 − µ3 (j + 1)2 ]
= 0
Finally, we investigate (4):
E[ 1
4µ4 (C2k+1 − C2k)(C2j+1 − C2j )]
Using the relation in equation 36, we have
E[C2kC2j ] = E[C2k( k2
j2 C2k +
2k
j CkBj,k +B
2 j,k)] (39)
= k2
j2 γ4 + γ2 (j − k)µ2 j2
(40)
= (2k + j − 3)µ22 + µ4
j2k (41)
E[ 1
4µ4 (C2k+1 − C2k)(C2j+1 − C2j )] = (2k + j − 3)µ22 + µ4 j2k − (2k + j − 2)µ 2 2 + µ4 (j + 1)2k
− (2k + j − 1)µ 2 2 + µ4
j2(k + 1) + (2k + j)µ22 + µ4 (j + 1)2(k + 1)
= (j2 − 5j − 3)µ22 j2(j + 1)2k(k + 1) − µ4 (j + 1)2k(k + 1)
= j2(µ22 − µ4)− (5j + 3)µ22
j2(j + 1)2k(k + 1)
=O(j−2k−2)

heading: In summary, E∆k∆j is O(k−2j−2).
text: A.4 GRADIENT OF SUMO
Assume that ∇θSUMO is bounded: it is sufficient that ∇θIWAE1 is bounded and that the sampling probabilities are chosen such that the partial sums of ∇θ∆kP(K≥k) converge, i.e. P(K ≥ k) > ck||∇θ∆k|| for some constant c. Then we have E [∇θSUMO(x)] = ∇θE [SUMO(x)] = ∇θ log pθ(x) directly by the dominated convergence theorem, as long as SUMO is everywhere differentiable, which is satisfied by all of our experiments. If ReLU neural networks are to be used, one may be able to show the same property using Theorem 5 of Bikowski et al. (2018), assuming finite higher moments and Lipschitz continuity. A.5 CONVERGENCE OF ∇ ( IWAEk+1 − IWAEk )
The IWAE log likelihood estimate is:
Lk = log (1 k k∑ i=1 pθ(x, zi) qψ(zi|x) ) The gradient of this with respect to λ, where λ is either θ or ψ, is
dLk dλ = 1∑k
i=1 pθ(x,zi) qψ(zi|x) k∑ i=1 d dλ pθ(x, zi) qψ(zi|x)
We abbreviate wi := pθ(x,zi) qψ(zi|x) , and νi = dwi dλ . In both λ = ψ and λ = θ cases, it suffices to treat the wi and νi as i.i.d. random variables with finite variance and expectation. Being a likelihood ratio, wi could be ill behaved when the importance sampling distribution qψ(zi|x) is is particularly mismatched from the true posterior p(zi|x) = pθ(x,ziEz∼p(z)pθ(x,z) . However, the analysis from IWAE (Burda et al., 2016) requires assuming that the likelihood ratios wi = pθ(x,zi) qψ(zi|x) are bounded, and we adopt this assumption. Reasoning about when this assumption holds, and the behavior of IWAE-like estimators when it does not, is an interesting area for future work.
Consider the differences between two gradients: we label ∆g as follows:
∆gk := dLk+1 dλ − dLk dλ
We have:
∆gk = 1∑k+1
i=1 wi νk+1 + ( 1∑k+1 i=1 wi − 1∑k i=1 wi ) k∑ i=1 νi
= 1∑k+1
i=1 wi νk+1 + wk+1(∑k+1 i=1 wi )(∑k i=1 wi ) k∑ i=1 νi
We again let Yk denote the kth sample mean 1k ∑ i wi. Then:
∆gk = 1
kYk νk+1 + wk+1 (k + 1)YkYk+1 ν̄k
The sample means Yk and µ̄k have finite expectation and variance. The variance vanishes as k →∞ (but the expectation does not change).
E||∆gk|| 2 2 =
1 k2 E||νk+1 Yk + k k + 1 wk+1ν̄k YkYk+1 ||22
Let νk+1 Yk + k k + 1 wk+1ν̄k YkYk+1 := φk
=⇒ E||∆gk|| 2 2 =
1
k2 ||Eφk||22 +
1
k2 Var(φk)
The second term vanishes at a rate strictly faster than 1k2 : the variance of φk goes to zero as k →∞. But the first term does not: φk is a biased estimator of φ∞ so Eφk does change with k, but it does not necessarily go to zero:
Eφ∞ = E [ ν Ew + k k + 1 wEν (Ew)2 ] = Eν Ew
Thus, E||∆gk||22 is at most O( 1 k2 ).
A.6 EMPIRICAL CONFIRMATION ON THE CONVERGENCE OF E∆2k AND E||∆ g k||22
We measure the ∆2k and ||∆ g k||22 on a toy example to verify the convergence rates empirically. We re-implement the toy Gaussian example from Rainforth et al. (2018b); Tucker et al. (2018). The generative model is pθ(x, z) = N (z|θ, I)N (x|z, I), where both x and z are in RD. The encoder is qφ(z|x) = N (z|Ax+ b, 23I), where φ = (A, b). The synthetic dataset was generated with D = 20 andN = 1000 data points using the true model parameter θtrue from a standard Gaussian. Alongside ||∆k||22, we plot several reference convergence rates such as O(1/kc), c > 1, and O(ck), c < 1, as a visual guide. The results are shown in Figure 5. Following the setup in Rainforth et al. (2018b), we sample a group of model parameters close to the optimal values which are perturbed by Gaussian noise from N (0, 0.012). The gradient ∆gk is taken w.r.t. the model parameter θ.
A.7 BIAS-VARIANCE TRADEOFF VIA GRADIENT CLIPPING
While SUMO is unbiased, its variance is extremely high or potentially infinite. This property leads to poor performance compared to lower bound estimates such as IWAE when maximizing loglikelihood. In order to obtain models with competitive log-likelihood values, we can make use of gradient clipping. This allows us to ignore rare gradient samples with extremely large values due to the heavy-tailed nature of its distribution.
Gradient clipping introduces bias in favor of reduced variance. Figure 6 shows how the performance changes as a function of the clipping value, and more importantly, the percentage of clipped gradients. As shown, neither full clipping nor no clipping are desirable. We performed this experiment after reporting the results in Table 1, so this grid search was not used to tune hyperparameters for our experiments. As bias is introduced, we do not use gradient clipping for entropy maximization or policy gradient (REINFORCE).
A.8 EXPERIMENTAL SETUP
In density modeling experiments, all the models are trained using a batch size of 100 and the AMSGrad optimizer (Reddi et al., 2018) with parameters lr = 0.001, β1 = 0.9, β2 = 0.999 and = 10−4. The learning rate is reduced by factor 0.8 if the validation likelihood does not improve for 50 epochs. We use gradient norm scaling in both the inference and generative networks. We train SUMO using the same architecture and hyperparameters as IWAE except the gradient clipping norm. We set the gradient norm to 5000 for encoder and {20, 40, 60} for decoder in SUMO. For IWAE, the gradient norm is fixed to 10 in all the experiments. We report the performance of models with early stopping if no improvements have been observed for 300 epochs on the validation set.
We add additional plots of the test NLL against the norm and percentage of gradients clipped for the decoder in Figure 6. The plot is based on MNIST with expected number of compute k = 15. Gradient clipping was not used in the other experiments except the density modeling ones, where it was used as a simple tool to obtain a better bias-variance trade-off.
A.8.1 REVERSE KL AND COMBINATORIAL OPTIMIZATION
These two tasks use the same encoder and decoder architecture: one hidden layer with tanh nonlinearities and 200 hidden units. We set the latent state to be of size 20. The prior is a standard Gaussian with diagonal covariance, while the encoder distribution is a Gaussian with parameterized diagonal covariance. For reverse KL, we used independent Gaussian conditional likelihoods for p(x|z), while for combinatorial optimization we used independent Bernoulli conditional distributions. We found it helps stablize training for both IWAE and SUMO to remove momentum and used RMSprop with learning rate 0.00005 and epsilon 1e-3 for fitting reverse KL. We used Adam with learning rate 0.001 and epsilon 1e-3, plus standard hyperparameters for the combinatorial optimization problems. SUMO used an expected compute of 15 terms.

Review summaries:
This paper proposes an unbiased estimator of $\log p_\theta(x)$. Many unbiased estimators of $p_\theta$ exist, but $\log p_\theta$ is needed in many other settings, some of which are not well-served by standard estimators of $p_\theta$. The SUMO estimator is essentially a Russian roulette-based extension of IWAE; it is exactly unbiased, but takes a random and unbounded number of samples.