Title: DOUBLY REPARAMETERIZED GRADIENT ESTIMATORS FOR MONTE CARLO OBJECTIVES

Abstract: Deep latent variable models have become a popular model choice due to the scalable learning algorithms introduced by (Kingma & Welling, 2013; Rezende et al., 2014) . These approaches maximize a variational lower bound on the intractable log likelihood of the observed data. Burda et al. (2015) introduced a multi-sample variational bound, IWAE, that is at least as tight as the standard variational lower bound and becomes increasingly tight as the number of samples increases. Counterintuitively, the typical inference network gradient estimator for the IWAE bound performs poorly as the number of samples increases (Rainforth et al., 2018; Le et al., 2018). Roeder et al. (2017) propose an improved gradient estimator, however, are unable to show it is unbiased. We show that it is in fact biased and that the bias can be estimated efficiently with a second application of the reparameterization trick. The doubly reparameterized gradient (DReG) estimator does not suffer as the number of samples increases, resolving the previously raised issues. The same idea can be used to improve many recently introduced training techniques for latent variable models. In particular, we show that this estimator reduces the variance of the IWAE gradient, the reweighted wake-sleep update (RWS) (Bornschein & Bengio, 2014) , and the jackknife variational inference (JVI) gradient (Nowozin, 2018) . Finally, we show that this computationally efficient, unbiased drop-in gradient estimator translates to improved performance for all three objectives on several modeling tasks.
Published as a conference paper at ICLR 2019 the bound is tighter, Rainforth et al. ( 2018) theoretically and empirically showed that the standard inference network gradient estimator for the IWAE bound performs poorly as the number of samples increases due to a diminishing signal-to-noise ratio (SNR). This motivates the search for novel gradient estimators. Roeder et al. (2017) proposed a lower-variance estimator of the gradient of the IWAE bound. They speculated that their estimator was unbiased, however, were unable to prove the claim. We show that it is in fact biased, but that it is possible to construct an unbiased estimator with a second application of the reparameterization trick which we call the IWAE doubly reparameterized gradient (DReG) estimator. Our estimator is an unbiased, computationally efficient drop-in replacement, and does not suffer as the number of samples increases, resolving the counterintuitive behavior from previous work (Rainforth et al., 2018) . Furthermore, our insight is applicable to alternative multisample training techniques for latent variable models: reweighted wake-sleep (RWS) (Bornschein & Bengio, 2014) and jackknife variational inference (JVI) (Nowozin, 2018) .
In this work, we derive DReG estimators for IWAE, RWS, and JVI and demonstrate improved scaling with the number of samples on a simple example. Then, we evaluate DReG estimators on MNIST generative modeling, Omniglot generative modeling, and MNIST structured prediction tasks. In all cases, we demonstrate substantial unbiased variance reduction, which translates to improved performance over the original estimators.

1. INTRODUCTION:

Following the influential work by (Kingma & Welling, 2013; Rezende et al., 2014) , deep generative models with latent variables have been widely used to model data such as natural images (Rezende & Mohamed, 2015; Kingma et al., 2016; Chen et al., 2016; Gulrajani et al., 2016) , speech and music time-series (Chung et al., 2015; Fraccaro et al., 2016; Krishnan et al., 2015) , and video (Babaeizadeh et al., 2017; Ha & Schmidhuber, 2018; Denton & Fergus, 2018) . The power of these models lies in combining learned nonlinear function approximators with a principled probabilistic approach, resulting in expressive models that can capture complex distributions. Unfortunately, the nonlinearities that empower these model also make marginalizing the latent variables intractable, rendering direct maximum likelihood training inapplicable. Instead of directly maximizing the marginal likelihood, a common approach is to maximize a tractable lower bound on the likelihood such as the variational evidence lower bound (ELBO) (Jordan et al., 1999; Blei et al., 2017) . The tightness of the bound is determined by the expressiveness of the variational family. For tractability, a factorized variational family is commonly used, which can cause the learned model to be overly simplistic. Burda et al. (2015) introduced a multi-sample bound, IWAE, that is at least as tight as the ELBO and becomes increasingly tight as the number of samples increases. Counterintuitively, although 2 BACKGROUND Our goal is to learn a latent variable generative model p θ (x, z) = p θ (z)p θ (x|z) where x are observed data and z are continuous latent variables. The marginal likelihood of the observed data, p θ (x) = p θ (x, z) dz, is generally intractable. Instead, we maximize a variational lower bound on log p θ (x) such as the ELBO

EQUATION

where q(z|x) is a variational distribution. Following the influential work by (Kingma & Welling, 2013; Rezende et al., 2014) , we consider the amortized inference setting where q φ (z|x), referred to as the inference network, is a learnable function parameterized by φ that maps from x to a distribution over z. The tightness of the bound is coupled to the expressiveness of the variational family (i.e., {q φ } φ ). As a result, limited expressivity of {q φ } φ , can negatively affect the learned model. Burda et al. (2015) introduced the importance weighted autoencoder (IWAE) bound which alleviates this coupling

EQUATION

with z 1:K ∼ i q φ (z i |x). The IWAE bound reduces to the ELBO when K = 1, is non-decreasing as K increases, and converges to log p θ (x) as K → ∞ under mild conditions (Burda et al., 2015) . When q φ is reparameterizable 1 , the standard gradient estimator of the IWAE bound is

∇ θ,φ E z 1:K log 1 K K i=1 w i = ∇ θ,φ E 1:K log 1 K K i=1 w i = E 1:K K i=1 w i j w j ∇ θ,φ log w i

where w i = p θ (x, z i )/q φ (z i |x). A single sample estimator of this expectation is typically used as the gradient estimator.

As K increases, the bound becomes increasingly tight, however, Rainforth et al. (2018) show that the signal-to-noise ratio (SNR) of the inference network gradient estimator goes to 0. This does not

1 Meaning that we can express zi as z( i, φ), where z is a deterministic, differentiable function and p( i) does not depend on θ or φ. This allows gradients to be estimated using the reparameterization trick (Kingma & Welling, 2013; Rezende et al., 2014) . Notably, the recently introduced implicit reparameterization trick has expanded the class of distributions for which we can compute reparameterization gradients to include Gaussian mixtures among many others (Graves, 2016; Jankowiak & Obermeyer, 2018; Jankowiak & Karaletsos, 2018; Figurnov et al., 2018) .

happen for the model parameters (θ). Following up on this work, Le et al. (2018) demonstrate that this deteriorates the performance of learned models on practical problems.

Because the IWAE bound converges to log p θ (x) (as K → ∞) regardless of q φ , φ's affect on the bound must diminish as K increases. It may be tempting to conclude that the SNR of the inference network gradient estimator must also decrease as K → ∞. However, low SNR is a limitation of the gradient estimator, not necessarily of the bound. Although the magnitude of the gradient converges to 0, if the variance of the gradient estimator decreases more quickly, then the SNR of the gradient estimator need not degrade. This motivates the search for lower variance inference network gradient estimators.

To derive improved gradient estimators for φ, it is informative to expand the total derivative2 of the IWAE bound with respect to φ

EQUATION

Previously, Roeder et al. (2017) found that the first term within the parentheses of Eq. 3 can contribute significant variance to the gradient estimator. When K = 1, this term analytically vanishes in expectation, so when K > 1 they suggested dropping it. Below, we abbreviate this estimator as STL. As we show in Section 6.1, the STL estimator introduces bias when K > 1.

3. DOUBLY REPARAMETERIZED GRADIENT ESTIMATORS (DREGS):

Our insight is that we can estimate the first term within the parentheses of Eq. 3 efficiently with a second application of the reparameterization trick. To see this, first note that

E 1:K K i=1 w i K j=1 w j ∂ ∂φ log q(z i |x) = K i=1 E 1:K w i K j=1 w j ∂ ∂φ log q(z i |x) ,

so it suffices to focus on one of the K terms. Because the derivative is a partial derivative ∂ ∂φ , it treats z i = z( i , φ) as a constant, so we can freely change the random variable that the expectation is over to z 1:K . Now,

EQUATION

where z -i = z 1:i-1,i+1:K is the set of z 1:K without z i . The inner expectation resembles a REIN-FORCE gradient term (Williams, 1992) , where we interpret wi j wj as the "reward". Now, we can use the following well-known equivalence between the REINFORCE gradient and the reparameterization trick gradient (See Appendix 8.1 for a derivation)

EQUATION

This holds even when f depends on φ. Typically, the reparameterization gradient estimator has lower variance than the REINFORCE gradient estimator because it directly takes advantage of the derivative of f . Applying the identity from Eq. 5 to the right hand side of Eq. 4 gives

EQUATION

This last expression can be efficiently estimated with a single Monte Carlo sample. When z i is not reparameterizable (e.g., the models in (Mnih & Rezende, 2016 )), we can use a control variate (e.g.,

1 K ∂ ∂φ log q φ (z i |x)).

In both cases, when K = 1, this term vanishes exactly and we recover the estimator proposed in (Roeder et al., 2017) for the ELBO. However, when K > 1, there is no reason to believe this term will analytically vanish.

Substituting Eq. 6 into Eq. 3, we obtain a simplification due to cancellation of terms

EQUATION

We call the algorithm that uses the single sample Monte Carlo estimator of this expression for the inference network gradient the IWAE doubly reparameterized gradient estimator (IWAE-DReG). This estimator has the property that when q(z|x) is optimal (i.e., q(z|x) = p(z|x)), the estimator vanishes exactly and has zero variance, whereas this does not hold for the standard IWAE gradient estimator. We provide an asymptotic analysis of the IWAE-DReG estimator in Appendix 8.2. The conclusion of that analysis is that, in contrast to the standard IWAE gradient estimator, the SNR of the IWAE-DReG estimator exhibits the same scaling behaviour of O( √ K) for both the generation and inference network gradients (i.e., improving in K).

4. ALTERNATIVE TRAINING ALGORITHMS:

Now, we review alternative training algorithms for deep generative models and derive their doubly reparameterized versions.

4.1. REWEIGHTED WAKE SLEEP (RWS):

Bornschein & Bengio (2014) introduced RWS, an alternative multi-sample update for latent variable models that uses importance sampling. Computing the gradient of the log marginal likelihood

∇ θ log p θ (x) = ∇ θ z p θ (x, z) dz p θ (x) = z p θ (x, z)∇ θ log p θ (x, z) dz p θ (x) = E p θ (z|x) [∇ θ log p θ (x, z)] ,

requires samples from p θ (z|x), which is generally intractable. We can approximate the gradient with a self-normalized importance sampling estimator

E p θ (z|x) [∇ θ log p θ (x, z)] ≈ E z 1:K K i=1 w i j w j ∇ θ log p θ (x, z i ) ,

where z 1:K ∼ i q φ (z i |x). Interestingly, this is precisely the same as the IWAE gradient of θ, so the RWS update for θ can be interpreted as maximizing the IWAE lower bound in terms of θ.

Instead of optimizing a joint objective for p and q, RWS optimizes a separate objective for the inference network. (Bornschein & Bengio, 2014) propose a "wake" update and a "sleep" update for the inference network. Le et al. ( 2018) provide empirical support for solely using the wake update for the inference network, so we focus on that update.

The wake update approximately minimizes the KL divergence from p θ (z|x) to q φ (z|x). The gradient of the KL term is

∇ φ E p θ (z|x) [log p θ (z|x) -log q φ (z|x)] = -E p θ (z|x) ∂ ∂φ log q φ (z|x) .

The wake update of the inference network approximates the intractable expectation by selfnormalized importance sampling

EQUATION

with z i ∼ q φ (z i |x). Le et al. (2018) note that this update does not suffer from diminishing SNR as K increases. However, a downside is that the updates for p and q are not gradients of a unified objective, so could potentially lead to instability or divergence.

DOUBLY REPARAMETERIZED REWEIGHTED WAKE UPDATE:

The wake update gradient for the inference network (Eq. 8) can be reparameterized

-E z 1:K K i=1 w i j w j ∂ ∂φ log q φ (z i |x) = E 1:K K i=1 w 2 i ( j w j ) 2 - w i j w j ∂ log w i ∂z i ∂z i ∂φ .

(9) We call the algorithm that uses the single sample Monte Carlo estimator of this expression as the wake update for the inference network RWS-DReG.

Interestingly, the inference network gradient estimator from (Roeder et al., 2017) can be seen as the sum of the IWAE gradient estimator and the wake update of the inference network (as the wake update minimizes, we add the negative of Eq. 9). Their positive results motivate further exploration of convex combinations of IWAE-DReG and RWS-DReG

EQUATION

)

We refer to the algorithm that uses the single sample Monte Carlo estimator of this expression as DReG(α). When α = 1, this reduces to RWS-DReG, when α = 0, this reduces to IWAE-DReG and when α = 0.5, this reduces STL.

4.2. JACKKNIFE VARIATIONAL INFERENCE (JVI):

Alternatively, Nowozin (2018) reinterprets the IWAE lower bound as a biased estimator for the log marginal likelihood. He analyzes the bias and introduces a novel family of estimators, Jackknife Variational Inference (JVI), which trade off reduction in bias for increased variance. This additional flexibility comes at the cost of no longer being a stochastic lower bound on the log marginal likelihood. The first-order JVI has significantly reduced bias compared to IWAE, which empirically results in a better estimate of the log marginal likelihood with fewer samples (Nowozin, 2018) . For simplicity, we focus on the first-order JVI estimator

K × E z1:K log 1 K K i=1 w i - K -1 K K i=1 E z-i   log   1 K -1 j =i w j     .

It is straightforward to apply our approach to higher order JVI estimators.

DOUBLY REPARAMETERIZED JACKKNIFE VARIATIONAL INFERENCE (JVI):

The JVI estimator is a linear combination of K and K -1 sample IWAE estimators, so we can use the doubly reparameterized gradient estimator (Eq. 7) for each term.

5. RELATED WORK:

Mnih & Rezende (2016) introduced a generalized framework of Monte Carlo objectives (MCO). The log of an unbiased marginal likelihood estimator is a lower bound on the log marginal likelihood by Jensen's inequality. In this view, the ELBO can be seen as the MCO corresponding to a single importance sample estimator of the marginal likelihood with q θ as the proposal distribution. Similarly, IWAE corresponds to the K-sample estimator. Maddison et al. (2017) show that the tightness of an MCO is directly related to the variance of the underlying estimator of the marginal likelihood.

However, Rainforth et al. (2018) point out issues with gradient estimators of multi-sample lower bounds. In particular, they show that although the IWAE bound is tighter, the standard IWAE gradient estimator's SNR scales poorly with large numbers of samples, leading to degraded performance. Le et al. (2018) experimentally investigate this phenomenon and provide empirical evidence of this degradation across multiple tasks. They find that RWS (Bornschein & Bengio, 2014) does not suffer from this issue and find that it can outperform models trained with the IWAE bound. We conclude that it is not sufficient to just tighten the bound; it is important to understand the gradient estimators of the tighter bound as well.

Wake-sleep is an alternative approach to fitting deep generative models, first introduced in (Hinton et al., 1995) as a method for training Hemholtz machines. It was extended to the multi-sample setting by (Bornschein & Bengio, 2014) and the sequential setting in (Gu et al., 2015) . It has been applied to generative modeling of images (Ba et al., 2015) .

6. EXPERIMENTS:

To evaluate DReG estimators, we first measure variance and signal-to-noise ratio (SNR)3 of gradient estimators on a toy example which we can carefully control. Then, we evaluate gradient variance and model learning on MNIST generative modeling, Omniglot generative modeling, and MNIST structured prediction tasks.

6.1 TOY GAUSSIAN We reimplemented the Gaussian example from (Rainforth et al., 2018) . Consider the generative model with z ∼ N (θ, I) and x|z ∼ N (z, I) and inference network q φ (z|x) ∼ N (Ax + b, 2 3 I), where φ = {A, b}. As in (Rainforth et al., 2018) , we sample a set of parameters for the model and inference network close to the optimal parameters (perturbed by zero-mean Gaussian noise with standard deviation 0.01), then estimate the gradient of the inference network parameters for increasing number of samples (K).

In addition to signal-to-noise ratio (SNR), we plot the squared bias and variance of the gradient estimators4 in Fig. 1 . The bias is computed relative to the expected value of the IWAE gradient estimator. As a result, although the average of K ELBO gradient estimators is an unbiased estimator of the ELBO gradient, it is a biased gradient estimator of the IWAE objective. Importantly, SNR does not penalize estimators that are biased, so trivial constant estimators can have infinite SNR. Thus, it is important to consider additional evaluation measures as well. As K increases, the SNR of the IWAE-DReG estimator increases, whereas the SNR of the standard gradient estimator of IWAE goes to 0, as previously reported. Furthermore, we can see the bias present in the STL estimator. As a check of our implementation, we verified that the observed "bias" for IWAE-DReG was statistically indistinguishable from 0 with a paired t-test. For the biased estimators (e.g., STL), we could easily reject the null hypothesis with few samples. Bias Squared 10 0 10 1 10 2 10 3 K 10 -10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0

Variance K ELBOs IWAE STL (Roeder et al. 2017 ) IWAE-DReG

Figure 1 : Signal-to-noise ratios (SNR), bias squared, and variance of gradient estimators with increasing K over 10 random trials with 1000 measurement samples per trial (mean in bold). The observed "bias" for IWAE-DReG is not statistically significant under a paired t-test (as expected because IWAE-DReG is unbiased). IWAE-DReG is unbiased, its SNR increases with K, and it has the lowest variance of the estimators considered here.

6.2. GENERATIVE MODELING:

Training generative models of the binarized MNIST digits dataset is a standard benchmark task for latent variable models. For this evaluation, we used the single latent layer architecture from (Burda et al., 2015) . The generative model used 50 Gaussian latent variables with an isotropic prior and passed z through two deterministic layers of 200 tanh units to parameterize factorized Bernoulli outputs. The inference network passed x through two deterministic layers of 200 tanh units to parameterize a factorized Gaussian distribution over z. Because our interest was in improved gradient estimators and optimization performance, we used the dynamically binarized MNIST dataset, which minimally suffers from overfitting. We used the standard split of MNIST into train, validation, and test sets.

We trained models with the IWAE gradient, the RWS wake update, and with the JVI estimator. In all three cases, the doubly reparameterized gradient estimator reduced variance5 and as a result substantially improved performance (Fig. 2 ).

0 1000 2000 3000 4000 5000

Steps (thousands) Steps (thousands) We found similar behavior with different numbers of samples (Fig. 3 and Appendix Fig. 8 ). Interestingly, the biased gradient estimators STL and RWS-DReG perform best on this task with RWS-DReG slightly outperforming STL. As observed in (Le et al., 2018) , RWS increasingly outperforms IWAE as K increases. Finally, we experimented with convex combinations of IWAE-DReG and RWS-DReG (right Fig. 3 ). On this dataset, convex combinations that heavily weighted RWS-DReG had the best performance. However, as we show below, this is task dependent.

Next, we performed the analogous experiment with the dynamically binarized Omniglot dataset using the same model architecture. Again, we found that the doubly reparameterized gradient estimator reduced variance and as a result improved test performance (Figs. 5 and 6 in the Appendix).

6.3. STRUCTURED PREDICTION ON MNIST:

Structured prediction is another common benchmark task for latent variable models. In this task, our goal is to model a complex observation x given a context c (i.e., model the conditional distribution p(x|c)). We can use a conditional latent variable model p θ (x, z|c) = p θ (x|z, c)p θ (z|c), however, as before, computing the marginal likelihood is generally intractable. It is straightforward to adapt the bounds and techniques from the previous section to this problem. To evaluate our method in this context, we use the standard task of modeling the bottom half of a binarized MNIST digit from the top half. We use a similar architecture, but now learn a conditional prior distribution p θ (z|c) where c is the top half of the MNIST digit. The conditional prior feeds c to two deterministic layers of 200 tanh units to parameterize a factorized Gaussian distribution over z. To model the conditional distribution p θ (x|c, z), we concatenate z with c and feed it to two deterministic layers of 200 tanh units to parameterize factorized Bernoulli outputs.

As in the previous tasks, the doubly reparameterized gradient estimator improves across all three updates (IWAE, RWS, and JVI; Appendix Fig. 7 ). However, on this task, the biased estimators (STL and RWS) underperform unbiased IWAE gradient estimators (Fig. 4 ). In particular, RWS becomes unstable later in training. We suspect that this is because RWS does not directly optimize a consistent objective. 

7. DISCUSSION:

In this work, we introduce doubly reparameterized estimators for the updates in IWAE, RWS, and JVI. We demonstrate that across tasks they provide unbiased variance reduction, which leads to improved performance. Furthermore, DReG estimators have the same computational cost as the original estimators. As a result, we recommend that DReG estimators be used instead of the typical gradient estimators.

Variational Sequential Monte Carlo (Maddison et al., 2017; Naesseth et al., 2018; Le et al., 2018) and Neural Adapative Sequential Monte Carlo (Gu et al., 2015) extend IWAE and RWS to sequential latent variable models, respectively. It would be interesting to develop DReG estimators for these approaches as well.

We found that a convex combination of IWAE-DReG and RWS-DReG performed best, however, the weighting was task dependent. In future work, we intend to apply ideas from (Baydin et al., 2017) to automatically adapt the weighting based on the data.

Finally, the form of the IWAE-DReG estimator (Eq. 7) is surprisingly simple and suggests that there may be a more direct derivation that is applicable to general MCOs.

8. APPENDIX:

0 1000 2000 3000 4000 5000

Steps (thousands) for a reparameterizable distribution q φ (z). To see this, note that

E q φ (z) f (z, φ) ∂ log q φ (z) ∂φ = E ∂f (z,

d dφ z q φ (z)f (z, φ) dz = z ∂ ∂φ q φ (z)f (z, φ) dz = z f (z, φ) ∂ ∂φ q φ (z) + q φ (z) ∂ ∂φ f (z, φ) dz = z f (z, φ)q φ (z) ∂ log q φ (z) ∂φ dz E q φ (z) ∂f (z, φ) ∂φ = E q φ (z) f (z, φ) ∂ log q φ (z) ∂φ + E q φ (z) ∂f (z, φ) ∂φ ,

via the REINFORCE gradient. On the other hand,

d dφ z q φ (z)f (z, φ) dz = d dφ E q φ (z) [f (z, φ)] = d dφ E [f (z( , φ), φ)] = E d dφ f (z( , φ), φ) = E ∂f (z, φ) ∂z ∂z( , φ) ∂φ + E ∂f (z, φ) ∂φ z=z( ,φ) = E ∂f (z, φ) ∂z ∂z( , φ) ∂φ + E q φ (z) ∂f (z, φ) ∂φ ,

via the reparameterization trick. Thus, we conclude that

E q φ (z) f (z, φ) ∂ log q φ (z) ∂φ +E q φ (z) ∂f (z, φ) ∂φ = E ∂f (z, φ) ∂z ∂z( , φ) ∂φ +E q φ (z) ∂f (z, φ) ∂φ ,

from which the identity follows.

8.2. ASYMPTOTIC ANALYSIS:

At a high level, Rainforth et al. (2018) show that the expected value of the IWAE gradient of the inference network collapses to zero with rate 1/K, while its standard deviation is only shrinking at a rate of 1/ √ K. This is the essence of the problem that results in the SNR (expectation divided by standard deviation) of the inference network gradients going to zero at a rate O(

(1/K)/(1/ √ K)) = O(1/ √ K)

, worsening with K. In contrast, Rainforth et al. (2018) show that the generation network gradients scales like O( √ K), improving with K.

Because the IWAE-DReG estimator is unbiased, we cannot hope to change the scaling of the expected value in K, but we can hope to change the scaling of the variance. In particular, in this subsection, we provide an informal argument, via the delta method, that the standard deviation of IWAE-DReG scales like K -3/2 , which results in an overall scaling of O( √ K) for the inference network gradient's SNR (i.e., increasing with K). Thus, the SNR of the IWAE-DReG estimator improves similarly in K for both inference and generation networks.

We will appeal to the delta method on a two-variable function g : R 2 → R. Define the following notation for the partials of g evaluated at the mean of random variables X, Y ,

g x (X, Y ) = ∂g(x, y) ∂x (x,y)=(E(X),E(Y ))

The delta method approximation of Var(g(X, Y )) is given by (Section 5.5 of Casella & Berger), Cov(X, Y )

EQUATION

K 4 + 4U 2 Z 6

Var(Y ) K 4 Because w i are all mutually independent, we get Var(Y ) = KVar(w i ). Similarly for Var(X) and u i . Because the w i and u i are identically distributed and independent for i = j, we have Cov(X, Y ) = KCov(w i , u i ). All together we can see that Var(g(X, Y )) scales like K -3 . Thus, the standard deviation scales like K -3/2 .

8.3. UNIFIED SURROGATE OBJECTIVES FOR ESTIMATORS:

In the main text, we assumed that θ and φ were disjoint, however, it can be helpful to share parameters between p and q (e.g., (Fraccaro et al., 2016) ). With the IWAE bound, we differentiate a single objective with to both the p and q parameters. Thus it is straightforward to adapt IWAE and IWAE-DReG to the shared parameter setting. In this section, we discuss how to deal with shared parameters in RWS. Suppose that both p and q are parameterized by θ. If we denote the unshared parameters of q by φ, then we can restrict the RWS wake update to only φ. Alternatively, with a modified RWS wake update, we can derive a single surrogate objective for each scenario such that taking the gradient with respect to θ results in the proper update. For clarity, we introduce the following modifier notation for p θ (x, z i ), q θ (z i |x), and w i which are functions of θ and z i = z(θ, i ). We use X to mean X with stopped gradients with respect to z i , X to mean X with stopped gradients with respect to θ (but not θ is not stopped in z(θ, i )), and X to mean X with stopped gradients for all variables. Then, we can use the following surrogate objectives: The only subtle difference is that DReG(α = 0.5) does not correspond exactly to STL due to the scaling between terms:

L DReG(α=0.5) (θ) = E 1:K K i=1 wi j wj (log pθ (x, z i ) + 0.5 log ŵi ) (17)

:

Implementation of DReG estimators and code to reproduce experiments: sites.google.com/view/ dregs.

log wi depends on φ in two ways: φ and zi = z( i, φ). The total derivative accounts for both sources of dependence and the partial derivative ∂ ∂φ considers zi as a constant.

Defined as the mean of the estimator divided by the standard deviation.

All dimensions of φ behaved qualitatively similarly, so for clarity, we show curves for a single randomly chosen dimension of φ.

We summarized the Covariance matrix of the gradient estimators by reporting the trace of the Covariance normalized by the number of parameters. To estimate this quantity, we estimated moments from exponential moving averages with decay=0.99 (we found that the results were robust to the exact value of the decay).

Figures and Tables:

Figure 2: MNIST generative modeling trained according to IWAE (left), RWS (middle), and JVI (right). The top row compares the variance of the original gradient estimator (dashed) with the variance of the doubly reparameterized gradient estimator (solid). The bottom row compares test performance. The left and middle plots show the IWAE (stochastic) lower bound on the test set. The right plot shows the JVI estimator (which is not a bound) on the test set. The bold lines are the average over three trials, and individual trials are displayed as semi-transparent). All methods used K = 64.

Figure3: Log-likelihood lower bounds for generative modeling on MNIST. The left and middle plots compare performance with different number of samples K = 32, 256. For clarity the legend is shared between the plots. The bold lines are the average over three trials, and individual trials are displayed as semi-transparent). The right plot compares performance as the convex combination between IWAE-DReG and RWS-DReG is varied (Eq. 10). To highlight differences, we plot the difference between the test IWAE bound and the test IWAE bound IWAE-DReG achieved at that step.

Figure4: Log-likelihood lower bounds for structured prediction on MNIST. The left plot uses K = 64 samples and the right plot uses K = 256 samples. For clarity the legend is shared between the plots. The bold lines are the average over three trials, and individual trials are displayed as semi-transparent). The right plot compares performance as the convex combination between IWAE-DReG and RWS-DReG is varied (Eq. 10). To highlight differences, we plot the difference between the test IWAE bound and the test IWAE bound IWAE-DReG achieved at that step.

Figure 5: Omniglot generative modeling trained according to IWAE (left), RWS (middle), and JVI (right). The top row compares the variance of the original gradient estimator (dashed) with the variance of the doubly reparameterized gradient estimator (solid). The bottom row compares test performance. The left and middle plots show the IWAE (stochastic) lower bound on the test set. The right plot shows the JVI estimator (which is not a bound) on the test set. The bold lines are the average over three trials, and individual trials are displayed as semi-transparent). All methods used K = 64.

Figure 7: Structured prediction on MNIST according to IWAE (left), RWS (middle), and JVI (right).The top row compares the variance of the original gradient estimator (dashed) with the variance of the doubly reparameterized gradient estimator (solid). The bottom row compares test performance. The left and middle plots show the IWAE (stochastic) lower bound on the test set. The right plot shows the JVI estimator (which is not a bound) on the test set. The bold lines are the average over three trials, and individual trials are displayed as semi-transparent). All methods used K = 64.

Y )g y (X, Y )Cov(X, Y ) + g y (X, Y ) 2 Var(Y )Now, assume without loss of generality that φ is a single real-valued parameter. Let u i = . Let g(X, Y ) = X/Y 2 , then g(X, Y ) is the IWAE-DReG estimator whose variance we seek to understand. Letting Z = E(w i ) and U = E(u i ) we get in this case after cancellations, Var(g(X, Y ))

(x, z i ) + log qθ (z i |x))

Peer review (iclr2019):
 Overall: This paper works on improving the gradient estimator of the ELBO. Author experimentally found that the estimator of the existing work(STL) is biased and proposed to reduce the bias by using the technique like REINFORCE. The problem author focused on is unique and the solution is simple, experiments show that proposed method seems promising. Clarity: The paper is clearly written in the sense that the motivation of research is clear, the derivation of the proposed method is easy to understand. Significance: I think this kind of research makes the variational inference more useful, so this work is significant. But I cannot tell the proposed method is really useful, so I gave this score. The reason I doubt the reason is that as I written in the below, the original STL can handle the mixture of Gaussians as the latent variable but the proposed method cannot. So I do not know which is better and whether I should use this method or use the original STL with flexible posterior distribution to tighten the evidence lower bound. I think additional experiments are needed. I know that motivation is a bit different for STL and proposed method but some comparisons are needed. Question and minor comments: In the original paper of STL, the author pointed out that by freezing the gradient of variational parameters to drop the score function term, we can utilize the flexible variational families like the mixture of Gaussians. In this work, since we do not freeze the variational parameters, we cannot utilize the mixture of Gaussians as in the STL. IWAE improves the lower bound by increasing the samples, but we can also improve the bound by specifying the flexible posteriors like the mixture of Gaussians in STL. Faced on this, I wonder which strategy is better to tighten the lower bound, should we use the STL with the mixture of Gaussians or use the proposed method? To clarify the usefulness of this method, I think the additional experimental comparisons are needed. About the motivation of the paper, I think it might be better to move the Fig.1 about the Bias to the introduction and clearly state that the author found that the STL is biased "experimentally". Followings are minor comments. In experiment 6.1, I*m not sure why the author present the result of K ELBO estimator in the plot of Bias and Variance. I think author want to point that when K=1, STL is unbiased with respect to the 1 ELBO, but when k>1, it is biased with respect to IWAE estimator. However, the objective of K ELBO and IWAE are different, it may be misleading. So this should be noted in the paper. In Figure 3, the left figure, what each color means? Is the color assignment is the same with the middle figure? (Same for Figure 4)