{'fraction_non_alphanumeric': 0.05040776412413742, 'fraction_numerical': 0.011624045167718366, 'mean_word_length': 4.385532591414944, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
DVAE#: Discrete Variational Autoencoders with Relaxed Boltzmann Priors

Boltzmann machines are powerful distributions that have been shown to be an effective prior over binary latent variables in variational autoencoders (VAEs). However, previous methods for training discrete VAEs have used the evidence lower bound and not the tighter importance-weighted bound. We propose two approaches for relaxing Boltzmann machines to continuous distributions that permit training with importance-weighted bounds. These relaxations are based on generalized overlapping transformations and the Gaussian integral trick. Experiments on the MNIST and OMNIGLOT datasets show that these relaxations outperform previous discrete VAEs with Boltzmann priors. An implementation which reproduces these results is available at https://github.com/QuadrantAI/dvae .


Introduction
Advances in amortized variational inference [1,2,3,4] have enabled novel learning methods [4,5,6] and extended generative learning into complex domains such as molecule design [7,8], music [9] and program [10] generation. These advances have been made using continuous latent variable models in spite of the computational efficiency and greater interpretability offered by discrete latent variables. Further, models such as clustering, semi-supervised learning, and variational memory addressing [11] all require discrete variables, which makes the training of discrete models an important challenge.
Prior to the deep learning era, Boltzmann machines were widely used for learning with discrete latent variables. These powerful multivariate binary distributions can represent any distribution defined on a set of binary random variables [12], and have seen application in unsupervised learning [13], supervised learning [14,15], reinforcement learning [16], dimensionality reduction [17], and collaborative filtering [18]. Recently, Boltzmann machines have been used as priors for variational autoencoders (VAEs) in the discrete variational autoencoder (DVAE) [19] and its successor DVAE++ [20]. It has been demonstrated that these VAE models can capture discrete aspects of data. However, both these models assume a particular variational bound and tighter bounds such as the importance weighted (IW) bound [21] cannot be used for training.
We remove this constraint by introducing two continuous relaxations that convert a Boltzmann machine to a distribution over continuous random variables. These relaxations are based on overlapping transformations introduced in [20] and the Gaussian integral trick [22] (known as the Hubbard-Stratonovich transform [23] in physics). Our relaxations are made tunably sharp by using an inverse temperature parameter.
VAEs with relaxed Boltzmann priors can be trained using standard techniques developed for continuous latent variable models. In this work, we train discrete VAEs using the same IW bound on the log-likelihood that has been shown to improve importance weighted autoencoders (IWAEs) [21]. This paper makes two contributions: i) We introduce two continuous relaxations of Boltzmann machines and use these relaxations to train a discrete VAE with a Boltzmann prior using the IW bound. ii) We generalize the overlapping transformations of [20] to any pair of distributions with computable probability density function (PDF) and cumulative density function (CDF). Using these more general overlapping transformations, we propose new smoothing transformations using mixtures of Gaussian and power-function [24] distributions. Power-function overlapping transformations provide lower variance gradient estimates and improved test set log-likelihoods when the inverse temperature is large. We name our framework DVAE# because the best results are obtained when the power-function transformations are sharp. 2

Related Work
Previous work on training discrete latent variable models can be grouped into five main categories: i) Exhaustive approaches marginalize all discrete variables [25,26] and which are not scalable to more than a few discrete variables.
ii) Local expectation gradients [27] and reparameterization and marginalization [28] estimators compute low-variance estimates at the cost of multiple function evaluations per gradient. These approaches can be applied to problems with a moderate number of latent variables.
iii) Relaxed computation of discrete densities [29] replaces discrete variables with continuous relaxations for gradient computation. A variation of this approach, known as the straight-through technique, sets the gradient of binary variables to the gradient of their mean [30,31].
iv) Continuous relaxations of discrete distributions [32] replace discrete distributions with continuous ones and optimize a consistent objective. This method cannot be applied directly to Boltzmann distributions. The DVAE [19] solves this problem by pairing each binary variable with an auxiliary continuous variable. This approach is described in Sec. 2.
v) The REINFORCE estimator [33] (also known as the likelihood ratio [34] or score-function estimator) replaces the gradient of an expectation with the expectation of the gradient of the score function. This estimator has high variance, but many increasingly sophisticated methods provide lower variance estimators. NVIL [3] uses an input-dependent baseline, and MuProp [35] uses a first-order Taylor approximation along with an input-dependent baseline to reduce noise. VIMCO [36] trains an IWAE with binary latent variables and uses a leave-one-out scheme to define the baseline for each sample. REBAR [37] and its generalization RELAX [38] use the reparameterization of continuous distributions to define baselines.
The method proposed here is of type iv) and differs from [19,20] in the way that binary latent variables are marginalized. The resultant relaxed distribution allows for DVAE training with a tighter bound. Moreover, our proposal encompasses a wider variety of smoothing methods and one of these empirically provides lower-variance gradient estimates.

Background
Let x x x represent observed random variables and ζ ζ ζ continuous latent variables. We seek a generative model p(x x x, ζ ζ ζ) = p(ζ ζ ζ)p(x x x|ζ ζ ζ) where p(ζ ζ ζ) denotes the prior distribution and p(x x x|ζ ζ ζ) is a probabilistic decoder. In the VAE [1], training maximizes a variational lower bound on the marginal log-likelihood: A probabilistic encoder q(ζ ζ ζ|x x x) approximates the posterior over latent variables. For continuous ζ ζ ζ, the bound is maximized using the reparameterization trick. With reparameterization, expectations with respect to q(ζ ζ ζ|x x x) are replaced by expectations against a base distribution and a differentiable function that maps samples from the base distribution to q(ζ ζ ζ|x x x). This can always be accomplished when q(ζ ζ ζ|x x x) has an analytic inverse cumulative distribution function (CDF) by mapping uniform samples through the inverse CDF. However, reparameterization cannot be applied to binary latent variables because the CDF is not differentiable.
The DVAE [19] resolves this issue by pairing each binary latent variable with a continuous counterpart. Denoting a binary vector of length D by z z z ∈ {0, 1} D , the Boltzmann prior is p(z z z) = e −E θ θ θ (z z z) /Z θ θ θ where E θ θ θ (z z z) = −a a a T z z z − 1 2 z z z T W W Wz z z is an energy function with parameters θ θ θ ≡ {W W W , a a a} and partition function Z θ θ θ . The joint model over discrete and continuous variables is p(x x x, z z z, ζ ζ ζ) = p(z z z)r(ζ ζ ζ|z z z)p(x x x|ζ ζ ζ) where r(ζ ζ ζ|z z z) = i r(ζ i |z i ) is a smoothing transformation that maps each discrete z i to its continuous analogue ζ i . DVAE [19] and DVAE++ [20] differ in the type of smoothing transformations r(ζ|z): [19] uses spike-and-exponential transformation (Eq. (1) left), while [20] uses two overlapping exponential distributions (Eq. (1) right). Here, δ(ζ) is the (one-sided) Dirac delta distribution, ζ ∈ [0, 1], and Z β is the normalization constant: The variational bound for a factorial approximation to the posterior where q(ζ ζ ζ|x x x) = i q(ζ i |x x x) and Here q(ζ i |x x x) = zi q(z i |x x x)r(ζ i |z i ) is a mixture distribution combining r(ζ i |z i = 0) and r(ζ i |z i = 1) with weights q(z i |x x x). The probability of binary units conditioned on ζ i , q(z z z|x x x, ζ ζ ζ) = i q(z i |x x x, ζ i ), can be computed analytically. H(q(z z z|x x x)) is the entropy of q(z z z|x x x). The second and third terms in Eq. (2) have analytic solutions (up to the log normalization constant) that can be differentiated easily with an automatic differentiation (AD) library. The expectation over q(ζ ζ ζ|x x x) is approximated with reparameterized sampling.
We extend [19,20] to tighten the bound of Eq. (2) by importance weighting [21,39]. These tighter bounds are shown to improve VAEs. For continuous latent variables, the K-sample IW bound is The tightness of the IW bound improves as K increases [21].

Model
We introduce two relaxations of Boltzmann machines to define the continuous prior distribution p(ζ ζ ζ) in the IW bound of Eq. (3). These relaxations rely on either overlapping transformations (Sec. 3.1) or the Gaussian integral trick (Sec. 3.2). Sec. 3.3 then generalizes the class of overlapping transformations that can be used in the approximate posterior q(ζ ζ ζ|x x x).

Overlapping Relaxations
We obtain a continuous relaxation of p(z z z) through the marginal p(ζ ζ ζ) = z p(z z z)r(ζ ζ ζ|z z z) where r(ζ ζ ζ|z z z) is an overlapping smoothing transformation [20] that operates on each component of z z z and ζ ζ ζ independently; i.e., r(ζ ζ ζ|z z z) = i r(ζ i |z i ). Overlapping transformations such as mixture of exponential in Eq. (1) may be used for r(ζ ζ ζ|z z z). These transformations are equipped with an inverse temperature hyperparameter β to control the sharpness of the smoothing transformation. As β → ∞, r(ζ ζ ζ|z z z) approaches δ(ζ ζ ζ − z z z) and p(ζ ζ ζ) = z p(z z z)δ(ζ ζ ζ − z z z) becomes a mixture of 2 D delta function distributions centered on the vertices of the hypercube in R D . At finite β, p(ζ ζ ζ) provides a continuous relaxation of the Boltzmann machine.
To train an IWAE using Eq. (3) with p(ζ ζ ζ) as a prior, we must compute log p(ζ ζ ζ) and its gradient with respect to the parameters of the Boltzmann distribution and the approximate posterior. This computation involves marginalization over z z z, which is generally intractable. However, we show that this marginalization can be approximated accurately using a mean-field model.

Computing log p(ζ ζ ζ) and its Gradient for Overlapping Relaxations
Since overlapping transformations are factorial, the log marginal distribution of ζ ζ ζ is The first term in Eq. (4) is the log partition function of the Boltzmann machinep(z z z) with augmented energy functionÊ β θ θ θ,ζ ζ ζ (z z z) Estimating the log partition function accurately can be expensive, particularly because it has to be done for each ζ ζ ζ. However, we note that each ζ i comes from a bimodal distribution centered at zero and one, and that the bias b b b β (ζ ζ ζ) is usually large for most components i (particularly for large β). In this case, mean field is likely to provide a good approximation ofp(z z z), a fact we demonstrate empirically in Sec. 4.
To compute log p(ζ ζ ζ) and its gradient, we first fit a mean-field distribution m(z z z) = i m i (z i ) by minimizing KL(m(z z z)||p(z z z)) [40]. The gradient of log p(ζ ζ ζ) with respect to β, θ θ θ or ζ ζ ζ is: where m m m T = [m 1 (z 1 = 1) · · · m D (z D = 1)] is the mean-field solution and where the gradient does not act on m m m. The first term in Eq. (5) is the result of computing the average energy under a factorial distribution. 3 The second expectation corresponds to the negative phase in training Boltzmann machines and is approximated by Monte Carlo sampling from p(z z z).
To compute the importance weights for the IW bound of Eq. (3) we must compute the value of log p(ζ ζ ζ) up to the normalization; i.e. the first term in Eq. (4). Assuming that KL m(z z z)||p(z z z) ≈ 0 and using the first term of Eq. (4) is approximated as H m(z z z) −Ê β θ θ θ,ζ ζ ζ (m m m).

The Gaussian Integral Trick
The computational complexity of log p(ζ ζ ζ) arises from the pairwise interactions z z z T W W Wz z z present in E θ θ θ (z z z). Instead of applying mean field, we remove these interactions using the Gaussian integral trick [41]. This is achieved by defining Gaussian smoothing: for an invertible matrix A A A and a diagonal matrix βI I I with β > 0. Here, β must be large enough so that [41]. However, neither of these choices places the modes of p(ζ ζ ζ) on the vertices of the hypercube in R D . Instead, we take A A A = (W W W + βI I I) −1 giving the smoothing transformation r(ζ ζ ζ|z z z) = N (ζ ζ ζ|z z z, (W W W + βI I I) −1 ). The joint density is then where 1 1 1 is the D-vector of all ones. Since p(z z z, ζ ζ ζ) no longer contains pairwise interactions z z z can be marginalized out giving where c i is the i th element of (W W W + βI I I)ζ ζ ζ.
The marginal p(ζ ζ ζ) in Eq. (7) is a mixture of 2 D Gaussian distributions centered on the vertices of the hypercube in R D with mixing weights given by p(z z z). Each mixture component has covariance Σ Σ Σ = (W W W + βI I I) −1 and, as β gets large, the precision matrix becomes diagonally dominant. As β → ∞, each mixture component becomes a delta function and p(ζ ζ ζ) approaches z p(z z z)δ(ζ ζ ζ − z z z). This Gaussian smoothing allows for simple evaluation of log p(ζ ζ ζ) (up to Z θ θ θ ), but we note that each mixture component has a nondiagonal covariance matrix, which should be accommodated when designing the approximate posterior q(ζ ζ ζ|x x x).
The hyperparameter β must be larger than the absolute value of the most negative eigenvalue of W W W to ensure that W W W + βI I I is positive definite. Setting β to even larger values has the benefit of making the Gaussian mixture components more isotropic, but this comes at the cost of requiring a sharper approximate posterior with potentially noisier gradient estimates.

Generalizing Overlapping Transformations
The previous sections developed two r(ζ ζ ζ|z z z) relaxations for Boltzmann priors. Depending on this choice, compatible q(ζ ζ ζ|x x x) parameterizations must be used. For example, if Gaussian smoothing is used, then a mixture of Gaussian smoothers should be used in the approximate posterior. Unfortunately, the overlapping transformations introduced in DVAE++ [20] are limited to mixtures of exponential or logistic distributions where the inverse CDF can be computed analytically. Here, we provide a general approach for reparameterizing overlapping transformations that does not require analytic inverse CDFs. Our approach is a special case of the reparameterization method for multivariate mixture distributions proposed in [42].
Assume q(ζ|x x x) = (1 − q)r(ζ|z = 0) + qr(ζ|z = 1) is the mixture distribution resulting from an overlapping transformation defined for one-dimensional z and ζ where q ≡ q(z = 1|x x x). Ancestral sampling from q(ζ|x x x) is accomplished by first sampling from the binary distribution q(z|x x x) and then sampling ζ from r(ζ|z). This process generates samples but is not differentiable with respect to q.
To compute the gradient (with respect to q) of samples from q(ζ|x x x), we apply the implicit function theorem. The inverse CDF of q(ζ|x x x) at ρ is obtained by solving: where ρ ∈ [0, 1] and R(ζ|z) is the CDF for r(ζ|z). Assuming that ζ is a function of q but ρ is not, we take the gradient from both sides of Eq. (8) with respect to q giving which can be easily computed for a sampled ζ if the PDF and CDF of r(ζ|z) are known. This generalization allows us to compute gradients of samples generated from a wide range of overlapping transformations. Further, the gradient of ζ with respect to the parameters of r(ζ|z) (e.g. β) is computed similarly as With this method, we can apply overlapping transformations beyond the mixture of exponentials considered in [20]. The inverse CDF of exponential mixtures is shown in Fig. 1(a) for several β. As β increases, the relaxation approaches the original binary variables, but this added fidelity comes at the cost of noisy gradients. Other overlapping transformations offer alternative tradeoffs: Uniform+Exp Transformation: We ensure that the gradient remains finite as β → ∞ by mixing the exponential with a uniform distribution. This is achieved by defining r (ζ|z) = (1 − )r(ζ|z) + where r(ζ|z) is the exponential smoothing and ζ ∈ [0, 1]. The inverse CDF resulting from this smoothing is shown in Fig. 1(b).
Power-Function Transformation: Instead of adding a uniform distribution we substitute the exponential distribution for one with heavier tails. One choice is the power-function distribution [24]: otherwise for ζ ∈ [0, 1] and β > 1.
The conditionals in Eq. (10) correspond to the Beta distributions B(1/β, 1) and B(1, 1/β) respectively. The inverse CDF resulted from this smoothing is visualized in Fig. 1(c). In the second row, the gradient of the inverse CDF with respect to q is visualized. Each column corresponds to a different smoothing transformation. As the transition region sharpens with increasing β, a sampling based estimate of the gradient becomes noisier; i.e., the variance of ∂ζ/∂q increases. The uniform+exp exponential has a very similar inverse CDF (first row) to the exponential but has potentially lower variance (bottom row). In comparison, the power-function smoothing with β = 40 provides a good relaxation of the discrete variables while its gradient noise is still moderate. See the supplementary material for a comparison of the gradient noise.

Experiments
In this section we compare the various relaxations for training DVAEs with Boltzmann priors on statically binarized MNIST [43] and OMNIGLOT [44] datasets. For all experiments we use a generative model of the form p(x x x, ζ ζ ζ) = p(ζ ζ ζ)p(x x x|ζ ζ ζ) where p(ζ ζ ζ) is a continuous relaxation obtained from either the overlapping relaxation of Eq. (4) or the Gaussian integral trick of Eq. (7). The underlying Boltzmann distribution is a restricted Boltzmann machine (RBM) with bipartite connectivity which allows for parallel Gibbs updates. We use a hierarchical autoregressively-structured q(ζ ζ ζ|x x x) = G g=1 q(ζ ζ ζ g |x x x, ζ ζ ζ <g ) to approximate the posterior distribution over ζ ζ ζ. This structure divides the components of ζ ζ ζ into G equally-sized groups and defines each conditional using a factorial distribution conditioned on x x x and all ζ ζ ζ from previous groups.
The smoothing transformation used in q(ζ ζ ζ|x x x) depends on the type of relaxation used in p(ζ ζ ζ). For overlapping relaxations, we compare exponential, uniform+exp, Gaussian, and power-function. With the Gaussian integral trick, we use shifted Gaussian smoothing as described below. The decoder p(x x x|ζ ζ ζ) and conditionals q(ζ ζ ζ g |x x x, ζ ζ ζ <g ) are modeled with neural networks. Following [20], we consider both linear (-) and nonlinear (∼) versions of these networks. The linear models use a single linear layer to predict the parameters of the distributions p(x x x|ζ ζ ζ) and q(ζ ζ ζ g |x x x, ζ ζ ζ <g ) given their input. The nonlinear models use two deterministic hidden layers with 200 units, tanh activation and batchnormalization. We use the same initialization scheme, batch-size, optimizer, number of training iterations, schedule of learning rate, weight decay and KL warm-up for training that was used in [20] (See Sec. 7.2 in [20]). For the mean-field optimization, we use 5 iterations. To evaluate the trained models, we estimate the log-likelihood on the discrete graphical model using the importance-weighted bound with 4000 samples [21]. At evaluation p(ζ ζ ζ) is replaced with the Boltzmann distribution p(z z z), and q(ζ ζ ζ|x x x) with q(z z z|x x x) (corresponding to β = ∞).
For DVAE, we use the original spike-and-exp smoothing. For DVAE++, in addition to exponential smoothing, we use a mixture of power-functions. The DVAE# models are trained using the IW bound in Eq. (3) with K = 1, 5, 25 samples. To fairly compare DVAE# with DVAE and DVAE++ (which can only be trained with the variational bound), we use the same number of samples K ≥ 1 when estimating the variational bound during DVAE and DVAE++ training.
With the Gaussian integral trick, each mixture component in the prior contains off-diagonal correlations and the approximation of the posterior over ζ ζ ζ should capture this. We recall that a multivariate Gaussian N (ζ ζ ζ|µ µ µ, Σ Σ Σ) can always be represented as a product of Gaussian conditionals i N ζ i |µ i + ∆µ i (ζ ζ ζ <i ), σ i where ∆µ i (ζ ζ ζ <i ) is linear in ζ ζ ζ <i . Motivated by this observation, we provide flexibility in the approximate posterior q(ζ ζ ζ|x x x) by using shifted Gaussian smoothing where r(ζ i |z i ) = N (ζ i |z i + ∆µ i (ζ ζ ζ <i ), 1/β i ), and ∆µ i (ζ ζ ζ <i ) is an additional parameter that shifts the distribution. As the approximate posterior in our model is hierarchical, we generate ∆µ i (ζ ζ ζ <g ) for the i th element in g th group as the output of the same neural network that generates the parameters of q(ζ ζ ζ g |x x x, ζ ζ ζ <g ). The parameter β i for each component of ζ ζ ζ g is a trainable parameter shared for all x x x.
Training also requires sampling from the discrete RBM to compute the θ θ θ-gradient of log Z θ θ θ . We have used both population annealing [45] with 40 sweeps across variables per parameter update and persistent contrastive divergence [46] for sampling. Population annealing usually results in a better generative model (see the supplementary material for a comparison). We use QuPA 4 , a GPU implementation of population annealing. To obtain test set log-likelihoods we require log Z θ θ θ , which we estimate with annealed importance sampling [47,48]. We use 10,000 temperatures and 1,000 samples to ensure that the standard deviation of the log Z θ θ θ estimate is small (∼ 0.01).
We compare the performance of DVAE# against DVAE and DVAE++ in Table 1. We consider four neural net structures when examining the various smoothing models. Each structure is denoted "G -/∼" where G represent the number of groups in the approximate posterior and -/∼ indicates linear/nonlinear conditionals. The RBM prior for the structures "1 -/∼" is 100 × 100 (i.e. D = 200) and for structures "2/4 ∼" the RBM is 200 × 200 (i.e. D = 400).
We make several observations based on Table 1: i) Most baselines improve as K increases. The improvements are generally larger for DVAE# as they optimize the IW bound. ii) Power-function smoothing improves the performance of DVAE++ over the original exponential smoothing. iii) DVAE# and DVAE++ both with power-function smoothing for K = 1 optimizes a similar variational bound with same smoothing transformation. The main difference here is that DVAE# uses the marginal p(ζ ζ ζ) in the prior whereas DVAE++ has the joint p(z z z, ζ ζ ζ) = p(z z z)r(z z z|ζ ζ ζ). For this case, it can be seen that DVAE# usually outperforms DVAE++ . iv) Among the DVAE# variants, the Gaussian integral trick and Gaussian overlapping relaxation result in similar performance, and both are usually inferior to the other DVAE# relaxations. v) In DVAE#, the uniform+exp smoothing performs better than exponential smoothing alone. vi) DVAE# with the power-function smoothing results in the best generative models, and in most cases outperforms both DVAE and DVAE++.
Given the superior performance of the models obtained using the mean-field approximation of Sec. 3.1.1 top(ζ ζ ζ), we investigate the accuracy of this approximation. In Fig. 2(a), we show that the mean-field model converges quickly by plotting the KL divergence of Eq. (6) with the number of mean-field iterations for a single ζ ζ ζ. To assess the quality of the mean-field approximation, in Fig. 2(b) we compute the KL divergence for randomly selected ζ ζ ζs during training at different iterations for exponential and power-function smoothings with different βs. As it can be seen, throughout the  training the KL value is typically < 0.2. For larger βs, the KL value is smaller due to the stronger bias that b b b β (ζ ζ ζ) imposes on z z z.
Lastly, we demonstrate that the lower variance of power-function smoothing may contribute to its success. As noted in Fig. 1, power-function smoothing potentially has moderate gradient noise while still providing a good approximation of binary variables at large β. We validate this hypothesis in Fig. 2(c) by measuring the variance of the derivative of the variational bound (with K = 1) with respect to the logit of q during training of a 2-layer nonlinear model on MNIST. When comparing the exponential (β = 10) to power-function smoothing (β = 30) at the β that performs best for each smoothing method, we find that power-function smoothing has significantly lower variance.

Conclusions
We have introduced two approaches for relaxing Boltzmann machines to continuous distributions, and shown that the resulting distributions can be trained as priors in DVAEs using an importance-weighted bound. We have proposed a generalization of overlapping transformations that removes the need for computing the inverse CDF analytically. Using this generalization, the mixture of power-function smoothing provides a good approximation of binary variables while the gradient noise remains moderate. In the case of sharp power smoothing, our model outperforms previous discrete VAEs.

A Population Annealing vs. Persistence Contrastive Divergence
In this section, we compare population annealing (PA) to persistence contrastive divergence (PCD) for sampling in the negative phase. In Table 2, we train DVAE# with the power-function smoothing on the binarized MNIST dataset using PA and PCD. As shown, PA results in a comparable generative model when there is one group of latent variables and better models in other cases.

B On the Gradient Variance of the Power-function Smoothing
Our experiments show that power-function smoothing performs best because it provides a better approximation of the binary random variables. We demonstrate this qualitatively in Fig. 1 and quantitatively in Fig. 2(c) of the paper. This is also visualized in Fig. 3. Here, we generate 10 6 samples from q(ζ) = (1 − q)r(ζ|z = 0) + qr(ζ|z = 1) for q = 0.5 using both the exponential and power smoothings with different values of β (β ∈ {8, 9, 10, . . . , 15} for exponential, and β ∈ {10, 20, 30, . . . , 80} for power smoothing). The value of β is increasing from left to right on each curve. The mean of |ζi − zi| (for zi = 1 [ζ i >0.5] ) vs. the variance of ∂ζi/∂q is visualized in this figure. For a given gradient variance, power function smoothing provides a closer approximation to the binary variables.  Figure 3: Average distance between ζ and its binarized z vs. variance of ∂ζ/∂q measured on 10 6 samples from q(ζ). For a given gradient variance, power function smoothing provides a closer approximation to the binary variables.