DOUBLY REPARAMETERIZED GRADIENT ESTIMATORS
Year: 2018

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.

heading: 1 INTRODUCTION
text: 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
Implementation of DReG estimators and code to reproduce experiments: sites.google.com/view/ dregs.
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.

heading: 2 BACKGROUND
text: 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
log pθ(x) = logEpθ(z)[pθ(x|z)] ≥ Eq(z|x) [ log pθ(x, z)
q(z|x)
] , (1)
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
Ez1:K
[ log ( 1
K K∑ i=1 pθ(x, zi) qφ(zi|x)
)] ≤ log pθ(x), (2)
with z1:K ∼ ∏ i qφ(zi|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 reparameterizable1, the standard gradient estimator of the IWAE bound is
∇θ,φEz1:K
[ log ( 1
K K∑ i=1 wi )] = ∇θ,φE 1:K [ log ( 1 K K∑ i=1 wi )] = E 1:K [ K∑ i=1 wi∑ j wj ∇θ,φ logwi ] where wi = pθ(x, zi)/qφ(zi|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
1Meaning 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 φ
E 1:K [ K∑ i=1 wi∑K j=1 wj ( − ∂ ∂φ log qφ(zi|x) + ∂ logwi ∂zi dzi dφ )] . (3)
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.

heading: 3 DOUBLY REPARAMETERIZED GRADIENT ESTIMATORS (DREGS)
text: 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 wi∑K j=1 wj ∂ ∂φ log q(zi|x) ] = K∑ i=1 E 1:K [ wi∑K j=1 wj ∂ ∂φ log q(zi|x) ] ,
so it suffices to focus on one of the K terms. Because the derivative is a partial derivative ∂∂φ , it treats zi = z( i, φ) as a constant, so we can freely change the random variable that the expectation is over to z1:K . Now,
Ez1:K [ wi∑ j wj ∂ ∂φ log qφ(zi|x) ] = Ez−iEzi [ wi∑ j wj ∂ ∂φ log qφ(zi|x) ] , (4)
where z−i = z1:i−1,i+1:K is the set of z1:K without zi. The inner expectation resembles a REINFORCE 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)
Eqφ(z|x) [ f(z) ∂
∂φ log qφ(z|x)
] = E [ ∂f(z)
∂z
∂z( , φ)
∂φ
] . (5)
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
Ezi [ wi∑ j wj ∂ ∂φ log qφ(zi|x) ] = E i [ ∂ ∂zi ( wi∑ j wj ) ∂zi ∂φ ]
= E i [( 1∑ j wj − wi ( ∑ j wj) 2 ) ∂wi ∂zi ∂zi ∂φ ] = E i [( wi∑ j wj − w 2 i ( ∑ j wj) 2 ) ∂ logwi ∂zi ∂zi ∂φ ] . (6)
This last expression can be efficiently estimated with a single Monte Carlo sample. When zi is not reparameterizable (e.g., the models in (Mnih & Rezende, 2016)), we can use a control variate (e.g.,
2logwi 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.
1 K ∂ ∂φ log qφ(zi|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, whenK > 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
∇φEz1:K
[ log ( 1
K K∑ i=1 wi )] = E 1:K  K∑ i=1 ( wi∑ j wj )2 ∂ logwi ∂zi ∂zi ∂φ  . (7) 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).

heading: 4 ALTERNATIVE TRAINING ALGORITHMS
text: Now, we review alternative training algorithms for deep generative models and derive their doubly reparameterized versions.

heading: 4.1 REWEIGHTED WAKE SLEEP (RWS)
text: 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) = Epθ(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
Epθ(z|x) [∇θ log pθ(x, z)] ≈ Ez1:K [ K∑ i=1 wi∑ j wj ∇θ log pθ(x, zi) ] ,
where z1:K ∼ ∏ i qφ(zi|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
∇φEpθ(z|x) [log pθ(z|x)− log qφ(z|x)] = −Epθ(z|x) [ ∂
∂φ log qφ(z|x)
] .
The wake update of the inference network approximates the intractable expectation by selfnormalized importance sampling
− Epθ(z|x) [ ∂
∂φ log qφ(z|x)
] ≈ −Ez1:K [ K∑ i=1 wi∑ j wj ∂ ∂φ log qφ(zi|x) ] , (8)
with zi ∼ qφ(zi|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
− Ez1:K [ K∑ i=1 wi∑ j wj ∂ ∂φ log qφ(zi|x) ] = E 1:K [ K∑ i=1 ( w2i ( ∑ j wj) 2 − wi∑ j wj ) ∂ logwi ∂zi ∂zi ∂φ ] .
(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
E 1:K [ K∑ i=1 ( α wi∑ j wj + (1− 2α) w 2 i ( ∑ j wj) 2 ) ∂ logwi ∂zi ∂zi ∂φ ] . (10)
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.

heading: 4.2 JACKKNIFE VARIATIONAL INFERENCE (JVI)
text: 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 × Ez1:K
[ log ( 1
K K∑ i=1 wi
)] − K − 1
K K∑ i=1 Ez−i log  1 K − 1 ∑ j 6=i wj  . 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.

heading: 5 RELATED WORK
text: 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).

heading: 6 EXPERIMENTS
text: 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.

heading: 6.1 TOY GAUSSIAN
text: 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, 23I), 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 ofK 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. AsK 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.

heading: 6.2 GENERATIVE MODELING
text: 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
3Defined as the mean of the estimator divided by the standard deviation. 4All dimensions of φ behaved qualitatively similarly, so for clarity, we show curves for a single randomly
chosen dimension of φ.
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).
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 RWSDReG 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).

heading: 6.3 STRUCTURED PREDICTION ON MNIST
text: 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.
5We 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).
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.

heading: 7 DISCUSSION
text: 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.

heading: ACKNOWLEDGMENTS
text: We thank Ben Poole and Diederik P. Kingma for helpful discussion and comments on drafts of this paper. We thank Sergey Levine and Jascha Sohl-Dickstein for insightful discussion.

heading: 8 APPENDIX

heading: 8.1 EQUIVALENCE BETWEEN REINFORCE GRADIENT AND REPARAMETERIZATION TRICK GRADIENT
text: Given a function f(z, φ), we have
Eqφ(z) [ f(z, φ) ∂ log qφ(z)
∂φ
] = E [ ∂f(z, φ)
∂z
∂z( , φ)
∂φ
] ,
for a reparameterizable distribution qφ(z). To see this, note that
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 + Eqφ(z) [ ∂f(z, φ) ∂φ ] = Eqφ(z) [ f(z, φ) ∂ log qφ(z)
∂φ
] + Eqφ(z) [ ∂f(z, φ)
∂φ
] ,
via the REINFORCE gradient. On the other hand,
d
dφ ∫ z qφ(z)f(z, φ) dz = d dφ Eqφ(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( , φ)
∂φ
] + Eqφ(z) [ ∂f(z, φ)
∂φ
] ,
via the reparameterization trick. Thus, we conclude that Eqφ(z) [ f(z, φ) ∂ log qφ(z)
∂φ
] +Eqφ(z) [ ∂f(z, φ)
∂φ
] = E [ ∂f(z, φ)
∂z
∂z( , φ)
∂φ
] +Eqφ(z) [ ∂f(z, φ)
∂φ
] ,
from which the identity follows.

heading: 8.2 ASYMPTOTIC ANALYSIS
text: 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 rateO((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 : R2 → R. Define the following notation for the partials of g evaluated at the mean of random variables X,Y ,
gx(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),
Var(g(X,Y )) ≈ gx(X,Y )2Var(X) + 2gx(X,Y )gy(X,Y )Cov(X,Y ) + gy(X,Y )2Var(Y ) Now, assume without loss of generality that φ is a single real-valued parameter. Let ui = w2i ∂ logwi ∂zi ∂zi ∂φ , X = ∑K i=1 ui, and Y = ∑K i=1 wi. 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(wi) and U = E(ui) we get in this case after cancellations,
Var(g(X,Y )) ≈ 1 Z4 Var(X) K4 − 4U Z5 Cov(X,Y ) K4 +
4U2
Z6 Var(Y ) K4
Because wi are all mutually independent, we get Var(Y ) = KVar(wi). Similarly for Var(X) and ui. Because the wi and ui are identically distributed and independent for i 6= j, we have Cov(X,Y ) = KCov(wi, ui). All together we can see that Var(g(X,Y )) scales like K−3. Thus, the standard deviation scales like K−3/2.

heading: 8.3 UNIFIED SURROGATE OBJECTIVES FOR ESTIMATORS
text: 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 respect 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, zi), qθ(zi|x), and wi which are functions of θ and zi = z(θ, i). We use X̃ to meanX with stopped gradients with respect to zi, 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:
IWAE:
LIWAE(θ) = E 1:K [ K∑ i=1 w̄i∑ j w̄j logwi ] (11)
DReG IWAE:
LDReG−IWAE(θ) = E 1:K  K∑ i=1 w̄i∑ j w̄j log p̃θ(x, zi) + ( w̄i∑ j w̄j )2 log ŵi  (12) RWS:
LRWS(θ) = E 1:K [ K∑ i=1 w̄i∑ j w̄j (log p̃θ(x, zi) + log q̃θ(zi|x)) ] (13)
DReG RWS:
LDReG−RWS(θ) = E 1:K  K∑ i=1 w̄i∑ j w̄j log p̃θ(x, zi) +  w̄i∑ j w̄j − ( w̄i∑ j w̄j )2 log ŵi  (14)
STL:
LSTL(θ) = E 1:K [ K∑ i=1 w̄i∑ j w̄j (log p̃θ(x, zi) + log ŵi) ] (15)
DReG(α):
LDReG(α)(θ) = E 1:K  K∑ i=1 w̄i∑ j w̄j log p̃θ(x, zi) + α w̄i∑ j w̄j + (1− 2α) ( w̄i∑ j w̄j )2 log ŵi 
(16)
The only subtle difference is that DReG(α = 0.5) does not correspond exactly to STL due to the scaling between terms:
LDReG(α=0.5)(θ) = E 1:K [ K∑ i=1 w̄i∑ j w̄j (log p̃θ(x, zi) + 0.5 log ŵi) ] (17)

Review summaries:
This paper applies a reparameterization trick to estimate the gradients objectives encountered in variational autoencoder based frameworks with continuous latent variables.  Especially the authors use this double reparameterization trick on Importance Weighted Auto-Encoder (IWAE) and Reweighted Wake-Sleep (RWS)  methods. Compared to IWAE, the developed method's SNR does not go to zero with increasing the number of particles.