NO MCMC FOR ME: AMORTIZED SAMPLING FOR FAST AND STABLE TRAINING OF ENERGY-BASED MODELS
Year: 2021

Abstract:
Energy-Based Models (EBMs) present a flexible and appealing way to represent uncertainty. Despite recent advances, training EBMs on high-dimensional data remains a challenging problem as the state-of-the-art approaches are costly, unstable, and require considerable tuning and domain expertise to apply successfully. In this work we present a simple method for training EBMs at scale which uses an entropy-regularized generator to amortize the MCMC sampling typically used in EBM training. We improve upon prior MCMC-based entropy regularization methods with a fast variational approximation. We demonstrate the effectiveness of our approach by using it to train tractable likelihood models. Next, we apply our estimator to the recently proposed Joint Energy Model (JEM), where we match the original performance with faster and stable training. This allows us to extend JEM models to semi-supervised classification on tabular data from a variety of continuous domains.

heading: 1 INTRODUCTION
text: Energy-Based Models (EBMs) have recently regained popularity within machine learning, partly inspired by the impressive results of Du & Mordatch (2019) and Song & Ermon (2020) on largescale image generation. Beyond image generation, EBMs have also been successfully applied to a wide variety of applications including: out-of-distribution detection (Grathwohl et al., 2019; Du & Mordatch, 2019; Song & Ou, 2018), adversarial robustness (Grathwohl et al., 2019; Hill et al., 2020; Du & Mordatch, 2019), reliable classification (Grathwohl et al., 2019; Liu & Abbeel, 2020) and semi-supervised learning (Song & Ou, 2018; Zhao et al.). Strikingly, these EBM approaches outperform alternative classes of generative models and rival hand-tailored solutions on each task.
Despite progress, training EBMs is still a challenging task. As shown in Table 1, existing training methods are all deficient in at least one important practical aspect. Markov chain Monte Carlo (MCMC) methods are slow and unstable during training (Nijkamp et al., 2019a; Grathwohl et al., 2020). Score matching mechanisms, which minimize alternative divergences are also unstable and most methods cannot work with discontinuous nonlinearities (such as ReLU) (Song & Ermon, 2019b; Hyvärinen, 2005; Song et al., 2020; Pang et al., 2020b; Grathwohl et al., 2020; Vincent, 2011). Noise contrastive approaches, which learn energy functions through density ratio estimation, typically don’t scale well to high-dimensional data (Gao et al., 2020; Rhodes et al., 2020; Gutmann & Hyvärinen, 2010; Ceylan & Gutmann, 2018).
∗Equal Contribution. Code available at github.com/wgrathwohl/VERA
In this work, we present a simple method for training EBMs which performs as well as previous methods while being faster and substantially easier to tune. Our method is based on reinterpreting maximum likelihood as a bi-level variational optimization problem, which has been explored in the past for EBM training (Dai et al., 2019). This perspective allows us to amortize away MCMC sampling into a GAN-style generator which is encouraged to have high entropy. We accomplish this with a novel approach to entropy regularization based on a fast variational approximation. This leads to the method we call Variational Entropy Regularized Approximate maximum likelihood (VERA).
Concretely, we make the following contributions:
• We improve the MCMC-based entropy regularizer of Dieng et al. (2019) with a parallelizable variational approximation.
• We show that an entropy-regularized generator can be used to produce a variational bound on the EBM likelihood which can be optimized more easily than MCMC-based estimators.
• We demonstrate that models trained in this way achieve much higher likelihoods than methods trained with alternative EBM training procedures.
• We show that our approach stabilizes and accelerates the training of recently proposed Joint Energy Models (Grathwohl et al., 2019).
• We show that the stabilization of our approach allows us to use JEM for semi-supervised learning, outperforming virtual adversarial training when little prior domain knowledge is available (e.g., for tabular data).

heading: 2 ENERGY BASED MODELS
text: An energy-based model (EBM) is any model which parameterizes a density as
pθ(x) = efθ(x)
Z(θ) (1)
where fθ : RD → R and Z(θ) = ∫ efθ(x)dx is the normalizing constant which is not explicitly modeled. Any probability distribution can be represented in this way for some fθ. The energy-based
parameterization has been used widely for its flexibility, ease of incorporating known structure, and relationship to physical systems common in chemistry, biology, and physics (Ingraham et al., 2019; Du et al., 2020; Noé et al., 2019).
The above properties make EBMs an appealing model class, but because they are unnormalized many tasks which are simple for alternative model classes become challenging for EBMs. For example, exact samples cannot be drawn and likelihoods cannot be exactly computed (or even lowerbounded). This makes training EBMs challenging as we cannot simply train them to maximize likelihood. The most popular approach to train EBMs is to approximate the gradient of the maximum likelihood objective. This gradient can be written as:
∇θ log pθ(x) = ∇θfθ(x)− Epθ(x′)[∇θfθ(x ′)]. (2)
MCMC techniques are used to approximately generate samples from pθ(x) (Tieleman, 2008). Practically, this approach suffers from poor stability and computational challenges from sequential sampling. Many tricks have been developed to overcome these issues (Du & Mordatch, 2019), but they largely still persist. Alternative estimators have been proposed to circumvent these challenges, including score matching (Hyvärinen, 2005), noise contrastive estimation (Gutmann & Hyvärinen, 2010), and variants thereof. These suffer from their own challenges in scaling to high dimensional data, and sacrifice the statistical efficiency of maximum likelihood.
In Figure 1 we visualize densities learned with our approach and Persistent Contrastive Divergence (Tieleman, 2008) (PCD) training. As we see, the sample quality of the PCD models is quite high but the learned density models do not match the true model. This is due to accrued bias in the gradient estimator from approximate MCMC sampling (Grathwohl et al., 2020). Prior work (Nijkamp et al., 2019b) has argued that this objective actually encourages the approximate MCMC samples to match the data rather than the density model. Conversely, we see that our approach (with proper entropy regularization) recovers a high quality model.

heading: 3 VARIATIONAL MAXIMUM LIKELIHOOD
text: We seek the energy function which maximizes likelihood given in Equation 1. We examine the intractable component of the log-likelihood, the log partition-function logZ(θ) = log ∫ efθ(x)dx. We can re-write this quantity as the optimum of
logZ(θ) = max q Eq(x)[fθ(x)] +H(q) (3)
where q is a distribution and H(q) = −Eq(x)[log q(x)] denotes its entropy 1 (see the Appendix A.1 for the derivation). Plugging this into our original maximum likelihood statement we obtain:
θ̂ = argmax θ [ Epdata(x)[fθ(x)]−maxq [ Eq(x)[fθ(x)] +H(q) ]] (4)
which gives us an alternative method for training EBMs. We introduce an auxiliary sampler qφ which we train online to optimize the inner-loop of Equation 4. This objective was used for EBM training in Kumar et al. (2019); Abbasnejad et al. (2019); Dai et al. (2017), Dai et al. (2019) (motivated by Fenchel Duality (Wainwright & Jordan, 2008)). Abbasnejad et al. (2019) use an implicit generative model and Dai et al. (2019) propose to use a sampler which is inspired by MCMC sampling from pθ(x) and whose entropy can be computed exactly.
Below we describe our approach which utilizes the same objective with a simpler sampler and a new approach to encourage high entropy. We note that when training pθ(x) and q(x) online together, the inner maximization will not be fully optimized. This leads our training objective for pθ(x) to be an upper bound on log pθ(x). In Section 5.1 we explore the impact of this fact and find that the bound is tight enough to train models that achieve high likelihood on high-dimensional data.
1For continuous spaces, this would be the differential entropy, but we simply use entropy here for brevity.

heading: 4 METHOD
text: We now present a method for training an EBM pθ(x) = efθ(x)/Z(θ) to optimize Equation 4. We introduce a generator distribution of the form qφ(x) = ∫ z qφ(x|z)q(z)dz such that:
qφ(x|z) = N (gψ(z), σ2I), q(z) = N (0, I) (5) where gψ is a neural network with parameters ψ and thus, φ = {ψ, σ2}. This is similar to the decoder of a variational autoencoder (Kingma & Welling, 2013). With this architecture it is easy to optimize the first and second terms of Equation 4 with reparameterization, but the entropy term requires more care.

heading: 4.1 ENTROPY REGULARIZATION
text: Estimating entropy or its gradients is a challenging task. Multiple, distinct approaches have been proposed in recent years based on Mutual Information estimation (Kumar et al., 2019), variational upper bounds (Ranganath et al., 2016), Denoising Autoencoders (Lim et al., 2020), and nearest neighbors (Singh & Póczos, 2016).
The above methods require the training of additional auxiliary models or do not scale well to high dimensions. Most relevant to this work are Dieng et al. (2019); Titsias & Ruiz (2019) which present a method for encouraging generators such as ours to have high entropy by estimating∇φH(qφ).The estimator takes the following form:
∇φH(qφ) = ∇φEqφ(x)[log qφ(x)] = ∇φEp(z)p( )[log qφ(x(z, ))] (Reparameterize sampling) = Ep(z)p( )[∇φ log qφ(x(z, ))] = Ep(z)p( )[∇x log qφ(x(z, ))T∇φx(z, )] (Chain rule) (6)
where we have written x(z, ) = gψ(z) + σ . All quantities in Equation 6 can be easily computed except for the score-function ∇x log qφ(x). The following estimator for this quantity can be easily derived (see Appendix A.2):
∇x log qφ(x) = Eqφ(z|x)[∇x log qφ(x|z)] (7) which requires samples from the posterior qφ(z|x) to estimate. Dieng et al. (2019); Titsias & Ruiz (2019) generate these samples using Hamiltonian Monte Carlo (HMC) (Neal et al., 2011), a gradientbased MCMC algorithm. As used in Dieng et al. (2019), 28 sequential gradient computations must be made per training iteration. Since a key motivation of this work is to circumvent the costly sequential computation of MCMC sampling, this is not a favourable solution. In our work we propose a more efficient solution that we find works just as well empirically.

heading: 4.2 VARIATIONAL APPROXIMATION WITH IMPORTANCE SAMPLING
text: We propose to replace HMC sampling of qφ(z|x) with a variational approximation ξ(z | z0) ≈ qφ(z | x) where z0 is a conditioning variable we will define shortly. We can use this approximation with self-normalized importance sampling to estimate
∇x log qφ(x) = Eqφ(z|x)[∇x log qφ(x | z)] = Ep(z0)qφ(z|x)[∇x log qφ(x | z)]
= Ep(z0)ξ(z|z0) [ qφ(z | x) ξ(z | z0) ∇x log qφ(x | z) ]
= Ep(z0)ξ(z|z0) [ qφ(z, x)
qφ(x)ξ(z | z0) ∇x log qφ(x | z) ] ≈
k∑ i=1 wi∑k j=1 wj ∇x log qφ(x | zi) ≡ ∇̃x log qφ(x; {zi}ki=1, z0) (8)
where {zi}ki=1 ∼ ξ(z | z0) and wi ≡ qφ(zi,x) ξ(zi|z0) . We use k = 20 importance samples for all experiments in this work. This approximation holds for any conditioning information we would like to use.
To choose this, let us consider how samples x ∼ qφ(x) are drawn. We first sample z0 ∼ N (0, I) and then x ∼ qφ(x | z0). In our estimator we want a variational approximation to qφ(z | x) and by construction, z0 is a sample from this distribution. For this reason we let our variational approximation be
ξη(z | z0) = N (z | z0, η2I), (9)
or simply a diagonal Gaussian centered at the z0 which generated x. For this approximation to be useful we must tune the variance η2. We do this by optimizing the standard Evidence Lower-Bound at every training iteration
LELBO(η; z0, x) = Eξη(z|z0) [log(qφ(x | z)) + log(qφ(z))] +H(ξη(z | z0)). (10)
We then use ξη(z | z0) to approximate ∇x log qφ(x) which we plug into Equation 6 to estimate ∇φH(qφ) for training our generator. A full derivation and discussion can be found in Appendix A.3. Combining the tools presented above we arrive at our proposed method which we call Variational Entropy Regularized Approximate maximum likelihood (VERA), outlined in Algorithm 1. We found it helpful to further add an `2-regularizer with weight 0.1 to the gradient of our model’s likelihood as in Kumar et al. (2019). In some of our larger-scale experiments we reduced the weight of the entropy regularizer as in Dieng et al. (2019). We refer to the entropy regularizer weight as λ.
Algorithm 1: VERA Training Input : EBM pθ(x) ∝ efθ(x), generator qφ(x, z), approximate posterior ξη(z|z0), entropy weight λ, gradient penalty γ = .1 Output: Parameters θ such that pθ ≈ p while True do
Sample mini-batch x, and generate mini-batch xg, z0 ∼ qφ(x, z) Compute LELBO(η; z0, xg) and update η // Update posterior Compute log fθ(x)− log fθ(xg) + γ||∇x log fθ(x)||2 and update θ // Update EBM Sample {zi}ki=1 ∼ ξη(z|z0) Compute s = ∇̃x log qφ(x; {zi}ki=1, z0) // Estimate score fn (Eq.8) Compute g = sT∇φxg // Estimate entropy gradient (Eq.6) Update φ with∇φ log fθ(xg) + λg // Update generator
end

heading: 5 EBM TRAINING EXPERIMENTS
text: We present results training various models with VERA and related approaches. In Figure 1 we visualize the impact of our generator’s entropy on the learned density model and compare this with MCMC sampling used in PCD learning. In Section 5.1, we explore this quantitatively by training tractable models and evaluating with test-likelihood. In Section 5.2 we explore the bias of our entropy gradient estimator and the estimator’s effect on capturing modes.

heading: 5.1 FITTING TRACTABLE MODELS
text: Optimizing the generator in VERA training minimizes a variational upper bound on the likelihood of data under our model. If this bound is not sufficiently tight, then training the model to maximize this bound will not actually improve likelihood. To demonstrate the VERA bound is tight enough to train large-scale models we train NICE models (Dinh et al., 2014) on the MNIST dataset. NICE is a normalizing flow (Rezende & Mohamed, 2015) model – a flexible density estimator which enables both exact likelihood computation and exact sampling. We can train this model with VERA (which does not require either of these abilities), evaluate the learned model using likelihood, and generate exact samples from the trained models. Full experimental details2 can be found in Appendix B.3 and additional results can be found in Appendix C.1.
2This experiment follows the NICE experiment in Song et al. (2020) and was based on their implementation.
We compare the performance of VERA with maximum likelihood training as well as a number of approaches for training unnormalized models; Maximum Entropy Generators (MEG) (Kumar et al., 2019), Persistent Contrastive Divergence (PCD), Sliced Score Matching (SSM) (Song et al., 2020), Denoising Score Matching (DSM) (Vincent, 2011), Curvature Propagation (CP-SM) (Martens et al., 2012), and CoopNets (Xie et al., 2018). As an ablation we also train VERA with the HMCbased entropy regularizer of Dieng et al. (2019), denoted VERA-(HMC). Table 2 shows that VERA outperforms all approaches that do not require a normalized model.
Figure 2 shows exact samples from our NICE models. For PCD we can see (as observed in Figure 1) that while the approximate MCMC samples resemble the data distribution, the true samples from the model do not. This is further reflected in the reported likelihood value which falls behind all methods besides DSM. CoopNets perform better than PCD, but exhibit the same behavior of generator samples resembling the data, but not matching true samples. We attribute this behavior to the method’s reliance on MCMC sampling.
Conversely, models trained with VERA generate coherent and diverse samples which reasonably capture the data distribution. We also see that the samples from the learned generator
much more closely match true samples from the NICE model than PCD and MEG. When we remove the entropy regularizer from the generator (λ = 0.0) we observe a considerable decrease in likelihood and we find that the generator samples are far less diverse and do not match exact samples at all. Intriguingly, entropy-free VERA outperforms most other methods. We believe this is because even without the entropy regularizer we are still optimizing a (weak) bound on likelihood. Conversely, the score-matching methods minimize an alternative divergence which will not necessarily correlate well with likelihood. Further, Figure 2 shows that MEG performs on par with entropy-free VERA indicating that the Mutual Information-based entropy estimator may not be accurate enough in high dimensions to encourage high entropy generators.

heading: 5.2 UNDERSTANDING OUR ENTROPY REGULARIZER
text: In Figure 3, we explore the quality of our score function estimator on a PCA (Tipping & Bishop, 1999) model fit to MNIST, a setting where we can compute the score function exactly (see Appendix B.4 for details). The importance sampling estimator (with 20 importance samples) has somewhat larger variance than the HMC-based estimator but has a notably lower bias of .12. The HMC estimator using 2 burn-in steps (recommended in Dieng et al. (2019)) has a bias of .48. Increasing the burn-in steps to 500 reduces the bias to .20 while increasing variance. We find the additional variance of our estimator is remedied by mini-batch averaging and the reduced bias helps explain the improved performance in Table 2.
Further, we compute the effective sample size (Kong, 1992) (ESS) of our importance sampling proposal on our CIFAR10 and MNIST models and achieve an ESS of 1.32 and 1.29, respectively using 20 importance samples. When an uninformed proposal (N (0, I)) is used, the ESS is 1.0 for both models. This indicates our gradient estimates are informative for training. More details can be found in Appendix B.6.
Next, we count the number of modes captured on a dataset with 1,000 modes consisting of 3 MNIST digits stacked on top of one another (Dieng et al., 2019; Kumar et al., 2019). We find that both VERA and VERA (HMC) recover 999 modes, but training with no entropy regularization recovers 998 modes. We conclude that entropy regularization is unnecessary for preventing mode collapse in this setting.

heading: 6 APPLICATIONS TO JOINT ENERGY MODELS
text: Joint Energy Models (JEM) (Grathwohl et al., 2019) are an exciting application of EBMs. They reinterpret standard classifiers as EBMs and train them as such to create powerful hybrid generative/discriminative models which improve upon purely-discriminative models at outof-distribution detection, calibration, and adversarial robustness.
Traditionally, classification tasks are solved with a function fθ : RD → RK which maps from the data to K unconstrained real-valued outputs (where k is the number of classes). This function parameterizes a distribution over labels y given data x: pθ(y|x) = ef(x)[y]/ ∑ y′ e fθ(x)[y ′]. The same function fθ can be used to define an EBM for the joint distribution over x and y as: pθ(x, y) = efθ(x)[y]/Z(θ). The label y can be marginalized out to give an unconditional model pθ(x) = ∑ y e
fθ(x)[y]/Z(θ). JEM models are trained to maximize the factorized likelihood:
log pθ(x, y) = α log pθ(y|x) + log pθ(x) (11) where α is a scalar which weights the two terms. The first term is optimized with cross-entropy and the second term is optimized using EBM training methods. In Grathwohl
et al. (2019) PCD was used to train the second term. We train JEM models on CIFAR10, CIFAR100, and SVHN using VERA instead of PCD. We examine how this change impacts accuracy, generation, training speed, and stability. Full experimental details can be found in Appendix B.7.
Speed and Stability While the results presented in Grathwohl et al. (2019) are promising, training models as presented in this work is challenging. MCMC sampling can be slow and training can easily diverge. Our CIFAR10 models train 2.8x faster than the official JEM implementation3 with the default hyper-parameters. With these default parameters JEM models would regularly diverge. To train for the reported 200 epochs training needed to be restarted multiple times and the number of MCMC steps needed to be quadrupled, greatly increasing run-time.
Conversely, we find that VERA was much more stable and our models never diverged. This allowed us to remove the additive Gaussian noise added to the data which is very important to stabilize MCMC training (Grathwohl et al., 2019; Nijkamp et al., 2019a; Du & Mordatch, 2019).
Hybrid Modeling In Tables 3 and 4 we compare the discriminative and generative performance of JEM models trained with VERA, PCD (JEM), and HDGE (Liu & Abbeel, 2020). With α = 1 we find that VERA leads to models with poor classification performance but strong generation performance. With α = 100 VERA obtains stronger classification performance than the original JEM model while still having improved image generation over JEM and HDGE (evaluated with FID (Heusel et al., 2017)).
3https://github.com/wgrathwohl/JEM
Model CIFAR10 CIFAR100 SVHN
Classifier 95.8 74.2 97.7
JEM 92.9 72.2 96.6 HDGE 96.7 74.5 N/A
VERA α = 100 93.2 72.2 96.8 VERA α = 1 76.1 48.7 94.2
Table 3: Classification on image datasets.
Model FID ↓ JEM 38.4 HDGE 37.6 SNGAN (Miyato et al., 2018) 25.50 NCSN (Song & Ermon, 2019b) 23.52
VERA α = 100 30.5 VERA α = 1 27.5
Table 4: FID on CIFAR10.
Unconditional samples can be seen from our CIFAR10 and CIFAR100 models in Figure 4. Samples are refined through a simple iterative procedure using the latent space of our generator, explained in Appendix B.7.1. Additional conditional samples can be found in Appendix C.5
Out-of-Distribution Detection JEM is a powerful approach for out-of-distribution detection (OOD), greatly outperforming tractable likelihood models like VAEs and flows (Nalisnick et al., 2018). In Table 5, reporting AUROC (Hendrycks & Gimpel, 2016), we see that for all but 1 dataset, VERA outperforms JEM with PCD training but under-performs contrastive training (HDGE). Intriguingly, VERA performs poorly on CelebA. This result, along with the unreliable performance of DSM models at this task (Li et al., 2019) leads to questions regarding special benefits of MCMC training that are lost in our method as well as DSM. We leave this to future work.

heading: 6.1 TABULAR DATA
text: Training with VERA is much more stable and easy to apply to domains beyond images where EBM training has been extensively tuned. To demonstrate this we show that JEM models trained with VERA can provide a benefit to semi-supervised classification on datasets from a variety of domains. Considerable progress has been made in semi-supervised learning but the most impressive results require considerable domain knowledge (Chen et al., 2020). In domains like images, text, and audio such knowledge exists but for data from particle accelerators, gyroscopes, and satellites, such intuition may not be available and these techniques cannot be applied. In these settings there are far fewer options for semi-supervised learning.
We present VERA as one such option. We train semi-supervised JEM models on data from a variety of continuous domains. We perform no data augmentation beyond removing redundant features and standardizing the remaining features. To further demonstrate the versatility of VERA we use an identical network for each dataset and method. Full experimental details can be found in Appendix B.8.
In Table 6, we find on each dataset tested, we find that VERA outperforms the supervised baseline and outperforms VAT which is the strongest domain agnostic semi-supervised learning method we are aware of.

heading: 7 RELATED WORK
text: Kumar et al. (2019) train EBMs using entropy-regularized generators, attempting to optimize the same objective as our own. The key difference is how the generator’s entropy is regularized. Kumar
4We treat MNIST as a tabular dataset since we do not use convolutional architectures.
et al. (2019) utilize a Mutual Information estimator to approximate the generator’s entropy whereas we approximate the gradients of the entropy directly. The method of Kumar et al. (2019) requires the training of an additional MI-estimation network, but our approach only requires the optimization of the posterior variance which has considerably fewer parameters. As demonstrated in Section 5.1, their approach does not perform as well as VERA for training NICE models and their generator collapses to a single point. This is likely due to the notorious difficulty of estimating MI in high dimensions and the unreliability of current approaches for this task (McAllester & Stratos, 2020; Song & Ermon, 2019a).
Dai et al. (2019) train EBMs using the same objective as VERA. The key difference here is in the architecture of the generator. Dai et al. (2019) use generators with a restricted form, inspired by various MCMC sampling methods. In this setting, a convenient estimator for the generator’s entropy can be derived. In contrast, our generators have an unconstrained architecture and we focus on entropy regularization in the unconstrained setting.
Abbasnejad et al. (2019) train EBMs and generators as well to minimize the reverse KL-divergence. Their method differs from ours in the architecture of the generator and the method for encouraging high entropy. Their generator defines an implicit density (unlike ours which defines a latent-variable model). The entropy is maximized using a series approximation to the generator function’s Jacobian log-determinant which approximates the change-of-variables for injective functions.
Gao et al. (2020) train EBMs using Noise Contrastive Estimation where the noise distribution is a normalizing flow. Their training objective differs from ours and their generator is restricted to having a normalizing flow architecture. These architectures do not scale as well as the GAN-style architectures we use to large image datasets.
As well there exist CoopNets (Xie et al., 2018) which cooperatively train an EBM and a generator network. Architecturally, they are similar to VERA but are trained quit differently. In CoopNets, the generator is trained via maximum likelihood on its own samples refined using MCMC on the EBM. This maximum likleihood step requires MCMC as well to generate posterior samples as in Pang et al. (2020a). In contrast, the generator in VERA is trained to minimize the reverse KL-divergence. Our method requires no MCMC and was specifically developed to alleviate the difficulties of MCMC sampling.
The estimator of Dieng et al. (2019) was very influential to our work. Their work focused on applications to GANs. Our estimator could easily be applied in this setting and to implicit variational inference (Titsias & Ruiz, 2019) as well but we leave this for future work.

heading: 8 CONCLUSION
text: In this work we have presented VERA, a simple and easy-to-tune approach for training unnormalized density models. We have demonstrated our approach learns high quality energy functions and models with high likelihood (when available for evaluation). We have further demonstrated the superior stability and speed of VERA compared to PCD training, enabling much faster training of JEM (Grathwohl et al., 2019) while retaining the performance of the original work. We have shown that VERA can train models from multiple data domains with no additional tuning. This enables the applications of JEM to semi-supervised classification on tabular data – outperforming a strong baseline method for this task and greatly outperforming JEM with PCD training.

heading: A KEY DERIVATIONS

heading: A.1 DERIVATION OF VARIATIONAL LOG-PARTITION FUNCTION
text: Here we show that the variational optimization given in equation 3 recovers logZ(θ).
max q Eq(x)[fθ(x)] +H(q)
= max q ∫ x q(x)fθ(x)dx− ∫ x q(x) log(q(x))dx
= max q ∫ x q(x) log ( exp(fθ(x)) q(x) ) dx
= max q ∫ x q(x) log ( exp(fθ(x)) q(x) ) dx− logZ(θ) + logZ(θ)
= max q ∫ x q(x) log ( exp(fθ(x))/Z(θ) q(x) ) dx+ logZ(θ) = max q −KL(q(x)‖pθ(x)) + logZ(θ)
= logZ(θ).

heading: A.2 SCORE FUNCTION ESTIMATOR
text: Here we derive the equivalent expression for∇x log qφ(x) given in Equation 7.
∇x log qφ(x) = ∇xqφ(x) qφ(x)
= ∇x ∫ z qφ(x, z)dz
qφ(x)
= ∫ z ∇xqφ(x, z)dz qφ(x)
= ∫ z ∇xqφ(x, z) qφ(x) dz
= ∫ z ∇xqφ(x | z)qφ(z) qφ(x) dz
= ∫ z (∇x log qφ(x | z))qφ(x | z)qφ(z) qφ(x) dz
= Eqφ(z|x)[∇x log qφ(x|z)].

heading: A.3 ENTROPY GRADIENT ESTIMATOR
text: From Equation 6 we have
∇φH(qφ) = Ep(z0)p( )[∇x log qφ(x) T∇φx(z0, )].
Plugging in our score function estimator gives
∇φH(qφ) = Ep(z0)p( )[∇x log qφ(x) T∇φx(z0, )] (12) = Ep(z0)p( ) [ Eqφ(z|x) [∇x log qφ(x | z)] T ∇φx(z0, ) ]
= Ep(z0)p( ) [ Eξη(z|z0) [ qφ(z | x) ξη(z | z0) ∇x log qφ(x | z) ]T ∇φx(z0, ) ]
= Ep(z0)p( ) [ Eξη(z|z0) [ qφ(x, z)
qφ(x)ξη(z | z0) ∇x log qφ(x | z)
]T ∇φx(z0, ) ]
≈ Ep(z0)p( ) [ k∑ i=1 wi∑k j=1 wj ∇x log qφ(x | zi) ]T ∇φx(z0, )  where {zi}ki=1 ∼ ξ(z | z0) and wi ≡ qφ(zi,x) ξ(zi|z0) .

heading: A.3.1 DISCUSSION
text: We discuss when the approximations in Equation 13 will hold. The importance sampling estimator will be biased when qψ(z | x) differs greatly from ξ(z | z0). Since the generator function gψ(z) is a smooth, Lipschitz function (as are most neural networks) and the output Gaussian noise is small, the space of z values which could have generated x should be concentrated near z0. In these settings, z0 should be useful for predicting q(z|x). The accuracy of this approximation is based on the dimension of z compared to x and the Lipschitz constant of gψ . In all settings we tested, dim(z) dim(x) where this approximation should hold. If dim(z) dim(x) the curse of dimensionality would take effect and z0 would be less and less informative about q(z|x). In settings such as these, we do not believe our approach would be as effective. Thankfully, almost all generator architectures we are aware of have dim(z) dim(x). The approximation could also break down if the Lipschitz constant blew up. We find this does not happen in practice, but this can be addressed with many forms of regularization and normalization.

heading: B EXPERIMENTAL DETAILS

heading: B.1 HYPERPARAMETER RECOMMENDATIONS
text: *When α > 1, learning rates were divided by α. †We found λ = 10−4 to work best on large image datasets, but in general we recommend starting with λ = 1 and trying successively smaller values of λ until training is stable.
We give some general tips on how to set hyperparameters when training VERA in Table 7. In all VERA experiments, we use the gradient norm penalty with weight 0.1. This was not tuned during our experiments. When using VERA and MEG we train with Adam (Kingma & Ba, 2014) and set β1 = 0, β2 = .9 as is standard in the GAN literature (Miyato et al., 2018). In general, we recommend setting the learning rate for the generator to twice the learning rate of the energy function and equal to the learning rate of the approximate posterior sampler.
B.1.1 IMPACT OF λ
Let us rewrite the generator’s training objective L(q;λ) = Eq(x)[fθ(x)] + λH(q). (13)
We can easily see that this objective is equivalent (up to a multiplicative constant) to
Eq(x)
[ fθ(x)
λ
] +H(q). (14)
From this, it is clear that maximizing Equation 14 is the same as minimizing the KL-divergence between q and a tempered version of pθ(x) defined as
efθ(x)/λ
Z . (15)
Tempering like this is standard practice in EBM training and is done in many recent works. Tempering has the effect of increasing the weight of the gradient signal in SGLD sampling relative to the added Gaussian noise. In all of Du & Mordatch (2019); Grathwohl et al. (2019); Nijkamp et al. (2019a;b) the SGLD samplers used a temperature of 1/λ = 20,000. This value is near the value of 10, 000 that we use in this work.
Thus we can see that our most important hyper-parameter is actually the temperature of the sampler’s target distribution. Ideally this temperature would be set to 1, but this can lead to unstable training (as it does with SGLD). To train high quality models, we recommend setting λ as close to 1 as possible, decreasing if training becomes unstable.

heading: B.2 TOY DATA VISUALIZATIONS
text: We train simple energy functions on a 2D toy data. The EBM is a 2-layer MLP with 100 hidden units per layer using the Leaky-ReLU nonlinearity with negative slope .2. The generator is a 2-layer
MLP with 100 hidden units per layer and uses batch normalization (Ioffe & Szegedy, 2015) with ReLU nonlinearities. All models were trained for 100,000 iterations with all learning rates set to .001 and used the Adam optimizer (Kingma & Ba, 2014).
The PCD models were trained using an SGLD sampler and a replay buffer with 10,000 examples, reinitialized every iteration with 5% probability. We used 20 steps of SGLD per training iteration to make runtime consistent with VERA. We tested σ values outside of the presented range but smaller values did not produce decent samples or energy functions and for larger values, training diverged.

heading: B.3 TRAINING NICE MODELS
text: The NICE models were exactly as in Song et al. (2020). They have 4 coupling layers and each coupling layer had 5 hidden layers. Each hidden layer has 1000 units and uses the Softplus nonlinearity. We preprocessed the data as in Song et al. (2020) by scaling the data to the range [0, 1], adding uniform noise in the range [−1/512, 1/512], clipping to the range [.001, .999] and applying the logit transform log(x) − log(1 − x). All models were trained for 400 epochs with the Adam optimizer (Kingma & Ba, 2014) with β1 = 0 and β2 = .9. We use a batch size of 128 for all models. We re-ran the score matching model of Song et al. (2020) to train for 400 epochs as well and found it did not improve performance as its best test performance happens very early in training.
For all generator-based training methods we use the same fixed generator architecture. The generator has a latent-dimension of 100 and 2 hidden layers with 500 units each. We use the Softplus nonlinearity and batch-normalization (Ioffe & Szegedy, 2015) as is common with generator networks.
For VERA the hyper-parameters we searched over were the learning rates for the NICE model and for the generator. Compared to Song et al. (2020) we needed to use much lower learning rates. We searched over learning rates in {.0003, .00003, .000003} for both the generator and energy function. We found .000003 to work best for the energy function and .0003 to work best for the generator. This makes intuitive sense since the generator needs to be fully optimized for the bound on likelihood to be tight. When equal learning rates were used (.0003, .0003) we observed high sample quality from the generator but exact samples and likelihoods from the NICE model were very poor.
For PCD we search over learning rates in {.0003, .00003, .000003}, the number of MCMC steps in {20, 40} and the SGLD noise standard-deviation in {1.0, 0.1}. All models with 20 steps and SGLD standard-deviation 1.0 quickly diverged. Our best model used learning rate .000003, stepsize 0.1, and 40 steps. We tested the gradient-norm regularizer from Kumar et al. (2019) and found it decreased performance for PCD trained models. Most models with 20 MCMC steps diverged early in training.
For review the MCMC sampler we use is stochastic gradient Langevin dynamics (Welling & Teh, 2011). This sampler updates its samples by
xt = xt−1 + σ2
2 ∇xfθ(x) + σ, ∼ N (0, I) (16)
where σ is the noise standard-deviation and is a parameter of the sampler.
For Maximum Entropy Generators (MEG) (Kumar et al., 2019) we must choose a mutualinformation estimation network. We follow their work and use an MLP with LeakyReLU nonlinearities with negative slope .2. Our network mirrors the generator and has 3 hidden layers with 500 units each. We searched over the same hyper-parameters as VERA. We found MEG to perform almost identically to training with no entropy regularization at all. We believe this has to do with the challenges of estimating MI in high dimensions Song & Ermon (2019a).
For CoopNets (Xie et al., 2018) we use the same flow and generator architectures as VERA. Following the MNIST experiments in Xie et al. (2018) we train using 10 SGLD steps per iteration. We tried the recommended learning rates of .007 and .0001 for the flow and generator, respectively but found this to lead quick divergence. For this reason, we search over learning rates for the flow and generator from {.0003, .00003, .000003} as we did for VERA and found the best combination to be .000003 for the flow and .0003 for the generator. Other combinations resulted in higher quality generator samples but much worse likelihood values. We tested the recommended SGLD step-size of .002 and found this to lead to divergence as well in this setup. Thus, we searched over larger values of {.002, .02, .1} found .1 to perform the best, as with PCD.

heading: B.4 ESTIMATION OF BIAS OF ENTROPY REGULARIZER
text: If we restrict the form of our generator to a linear function
x = Wz + µ+ σ , z, ∼ N (0, I) then we have
q(x|z) = N (Wz + µ, σ2I), q(x) = N (µ,WTW + σ2I) meaning we can exactly compute log q(x), and∇x log q(x) which is the quantity that VERA (HMC) approximates with the HMC estimator from Dieng et al. (2019) and we approximate with VERA. To explore this, we fit a PCA model on MNIST and recover the parameters W,µ of the linear generator and the noise parameter σ. Samples from this model can be seen in Figure 5.
Both VERA and VERA (HMC) have some parameters which are tuned automatically to improve the estimator throughout training. For VERA this is the posterior variance which is optimized according to Equation 10 with Adam with default hyperparameters. For VERA (HMC) this is the stepsize of HMC which is tuned automatically as outlined in (Dieng et al., 2019). Both estimators were trained with a learning rate of .01 for 500 iterations with a batch size of 5000. Samples from the estimators during training were taken with the default parameters, in particular the number of burn-in steps or the number of importance samples was not varied as they were during evaluation.
The bias of the estimators were evaluated on a batch of 10 samples from the generator. For each example in the batch, 5000 estimates were taken from the estimator and averaged to be taken as an estimate of the score function for this batch example. This estimate of the score function was subtracted from the true score function and then averaged over all dimensions and examples in the batch and taken as an estimate of the bias per dimension.
For VERA (HMC) we varied the number of burn-in steps used for samples to evaluate the bias of the estimator. We also tried to increase the number of posterior samples taken off the chain, but we found that this did not clearly reduce the bias of this estimator as the number of samples increased. For VERA we computed the bias on the default of 20 importance samples.

heading: B.5 MODE COUNTING
text: We train VERA without labels on the STACKEDMNIST dataset. This dataset consists of 60,000 examples of 3 stacked MNIST images sampled with replacement from the original 60,000 MNIST train set. As in (Dieng et al., 2019) we resize the images to 64 × 64 and use the DCGAN (Radford et al., 2015) architecture for both the energy-function and generator. We use a latent code of 100 dimensions. We train for 17 epochs with a learning rate of .001 and batch size 100.
We estimate the number of modes captured by taking S = 10, 000 samples as in (Dieng et al., 2019) and classifying each digit of the 3 stacked images separately with a pre-trained classifier on MNIST.

heading: B.6 EFFECTIVE SAMPLE SIZE
text: When performing importance sampling, the quality of the proposal distribution has a large impact. If the proposal is chosen poorly, then typically 1 sample will dominate in the expectation. This can
be quantified using the effective sample size (Kong, 1992) (ESS) which is defined as
wi = p(x)
q(x)
w̃i = wi∑N j=1 wj
ESS = 1∑N
i=1 w̃i 2
(17)
where p(x) is the target distribution, q(x) is the proposal distribution and N is the number of samples. If the self-normalized importance weights are dominated by one weight close to 1 then the ESS will be 1. If the proposal distribution is identical to the target, so that the self-normalized importance weights are then uniform, then the ESS will be N . When ESS = N , importance sampling is as efficient as using the target distribution.
To understand the effect of the proposal distribution on ESS we plot the ESS when doing importance sampling with 20 importance samples from a 128-dimensional Gaussian target distribution with µ = 0 and Σ = I . We use a proposal which is a 128-dimensional Gaussian with µ increasing from 0 to 5. Results can be seen in Figure 6. We see when the means differ by greater than 2, the ESS is approximately 1.0 and importance sampling has effectively failed.

heading: B.7 JEM MODELS
text: Our energy function used the same architecture as in Grathwohl et al. (2019); Liu & Abbeel (2020), a Wide ResNet Zagoruyko & Komodakis (2016) 28-10 with batch-norm removed and with dropout. Our generator architecture is identical to Miyato et al. (2018). As in Grathwohl et al. (2019); Liu & Abbeel (2020) we set the learning rate for the energy function equal to 0.0001. We set the learning rate for the generator equal to 0.0002. We train for 200 epochs using the Adam optimizer with β1 = 0 and β2 = .9. We set the batch size to 64. Results presented are from the models after 200 epochs with no early stopping. We believe better results could be obtained from further training.
We trained models with α ∈ {1, 30, 100} and found classification to be best with α = 100 and generation to be best with α = 1.
Prior work on PCD EBM training (Grathwohl et al., 2019; Du & Mordatch, 2019; Nijkamp et al., 2019a;b) recommends adding Gaussian noise to the data to stabilize training. Without this, PCD training of JEM models very quickly diverges. Early in our experiments we found training with VERA was stable without the addition of Gaussian noise so we do not use it.
As mentioned in Dieng et al. (2019), when the strength of the entropy regularizer λ is too high, the generator may fall into a degenerate optimum where it just outputs high-entropy Gaussian noise. To combat this, as suggested in Dieng et al. (2019), we decrease the strength of λ to .0001 for all JEM experiments. This value was chosen by decreasing λ from 1.0 by factors of 10 until learning took place (quantified by classification accuracy).

heading: B.7.1 MCMC SAMPLE REFINEMENT
text: Our generator qφ(x) is trained to approximate our EBM pθ(x). After training, the samples from the generator are of high quality (see Figures 9 and 10, left) but they are not exactly samples from pθ(x). We can use MCMC sampling to improve the quality of these samples. We use a simple MCMC refinement procedure based on the Metropolis Adjusted Langevin Algorithm (Besag, 1994) applied to an expanded state-space defined by our generator and perform the Accept/Reject step in the data space.
We can reparameterize a generator sample x ∼ qφ(x) as a function x(z, ) = gψ(z) + σ, and we can define an unnormalized density over {z, }, log h(z, ) ≡ fθ(gψ(z) + σ)− logZ(θ, φ) which is the density (under pθ(x)) of the generator sample.
Starting with an initial sample z0, 0 ∼ N (0, I) we define the proposal distribution p(zt|zt−1, t−1) = N ( zt−1 + δ
2 ∇zt−1 log h(zt−1, t−1), δ2I ) p( t|zt−1, t−1) = N ( zt−1 + δ
2 ∇ t−1 log h(zt−1, t−1), δ2I ) p(zt, t|zt−1, t−1) = p(zt|zt−1, t−1)p( t|zt−1, t−1)
and accept a new sample with probability
min [ h(zt, t)p(zt, t|zt−1, t−1) h(zt−1, t)p(zt−1, t−1|zt, t) , 1 ] .
Here, δ is the step-size and is a parameter of the sampler. We tune δ with a burn-in period to target an acceptance rate of 0.57.
A similar sampling procedure was proposed in Kumar et al. (2019) and Che et al. (2020) and in both works was found to improve sample quality. In all experiments,
We clarify that this procedure is not a valid MCMC sampler for pθ(x) due to the augmented variables and the change in density of gψ which are not corrected for. The density of the samples will be a combination of pθ(x) and qφ(x). As the focus of this work was training and not sampling/generation, we leave the development of more correct generator MCMC-sampling to future work. Regardless, we find this procedure to improve visual sample quality. In Figure 7 we visualize a sampling chain using the above method applied to our JEM model trained on SVHN.
B.7.2 IMAGE QUALITY METRICS
We present results on Inception Score (Salimans et al., 2016) and Frechet Inception Distance (Heusel et al., 2017). These metrics are notoriously fickle and different repositories are known to give very different results (Grathwohl et al., 2019). For these evaluations we generate 12,800 samples from the model and (unless otherwise stated) refine the samples with 100 steps of our latent-space MALA procedure (Appendix B.7.1). The code to generate our reported FID comes from this publicly available repository. The code to generate our reported Inception Score can be found here.

heading: B.8 SEMI-SUPERVISED LEARNING ON TABULAR DATA
text: B.8.1 DATA
We provide details about each of the datasets used for the experiments in Section 6.1. HEPMASS5 is a dataset obtained from a particle accelerator where we must classify signal from background noise. CROP6 is a dataseset for classifying crop types from optical and radar sensor data. HUMAN7 is a dataset for human activity classification from gyroscope data. MNIST is an image dataset of handwritten images, treated here as tabular data.
For HEPMASS and HUMAN data we remove features which repeat the same exact same value more than 5 times. For CROP data we remove features which have covariation greater than 1.01. For MNIST we linearly standardize features to the interval [−1, 1]. We take a random 10% subset of the data to use as a validation set.

heading: B.8.2 TRAINING
text: We use the same architecture for all experiments and baselines. It has 6 layers of hidden units with dimensions [1000, 500, 500, 250, 250, 250] and a Leaky-ReLU nonlinearity with negative-slope .2 between each layer of hidden units. The only layers which change between datasets are the input layer and the output layer which change according to the number of features and number of classses respectively.
The training process for semi-supervised learning is similar to JEM with an additional objective commonly used in semi-supervised learning:
log pθ(x, y) = α log pθ(y | x) + log pθ(x) + βH(pθ(y | x)) (18)
where H(p(y|x)) is the entropy of the predictive distribution over the labels. For all models we report the accuracy the model converged to on the held-out validation set. We report the average taken over three training runs with different seeds.
We use equal learning rates for the energy model, generator, and the entropy estimator. We tune the learning rate and decay schedule for supervised models on the full-set of labels and 10 labels per class.
5http://archive.ics.uci.edu/ml/datasets/HEPMASS 6https://archive.ics.uci.edu/ml/datasets/Crop+mapping+using+fused+
optical-radar+data+set 7https://archive.ics.uci.edu/ml/datasets/Human+Activity+Recognition+ Using+Smartphones
On VERA we tune the learning rate and decay schedule, the weighting of the entropy regularization λ and the weighting of the entropy of classification outputs β.
For VAT we tune the perturbation size ∈ {.01, .1, 1, 3, 10}. All other hyperparameters were fixed according to tuning on VERA.
For MEG we used the hyperparameters tuned according to VERA.
For JEM we tune the number of MCMC steps in the range κ ∈ {20, 40, 80}. We generate samples using SGLD with stepsize 1 and noise standard deviation 0.01 as in Grathwohl et al. (2019).

heading: C ADDITIONAL RESULTS

heading: C.1 TRAINING MIXTURES OF GAUSSIANS
text: We present additional results training mixure of Gaussian models using VERA and PCD. Each model consits of 100 diagonal Gaussians and mixing weights. Our experimental setup and hyperparameter search was identical to that presented in Appendiex B.2. We see in Table 10 that VERA outperforms PCD.

heading: C.2 SAMPLES FROM SSL MODELS
text: We present some samples from our semi-supervised MNIST models in Figure 8.

heading: C.3 HYBRID MODELING
text: We present an extended Table 4 with Inception Score Salimans et al. (2016) and which includes more comparisons.
Model FID IS
JEM 38.4 8.76 HDGE 37.6 9.19 SNGAN (Miyato et al., 2018) 25.50 8.59 NCSN (Song & Ermon, 2019b) 23.52 8.91 ADE (Dai et al., 2019) N/A 7.55 IGEBM (Du & Mordatch, 2019) 37.9 8.30 Glow (Kingma & Dhariwal, 2018) 48.9 3.92 FCE (Gao et al., 2020) 37.3 N/A
VERA α = 100 30.5 8.11 VERA α = 1 27.5 8.00 VERA α = 1 (generator samples) 32.4 7.34

heading: C.5 UNCONDITIONAL IMAGE GENERATION: CIFAR100

heading: C.4 UNCONDITIONAL IMAGE GENERATION: CIFAR10

heading: C.6 CONDITIONAL IMAGE GENERATION: CIFAR10 AND CIFAR100

Review summaries:
This paper proposes a new method on training energy-based models with maximum likelihood. Instead of using MCMC approaches to sample from the EBM, authors follow previous work on training neural generators for faster sample generation. In particular, authors consider a special generator where the output is convolved with Gaussian noise. The score function of this generator can be estimated with self-normalized importance sampling, which is then used to estimate the entropy term through the reparameterization trick. Authors demonstrate that their new method is able to train EBMs efficiently, and improves the stability and performance of JEMs compared to MCMC-based training approaches.

This paper presents a method for improving training of energy-based models. Rather than drawing samples using persistent contrastive divergence / MCMC, this approach parameterizes a separate model, which is trained to directly output samples. This effectively adds an additional KL divergence to the objective. The authors use a particular form of sampling model (a latent Gaussian model), borrowing a few tricks for getting entropy estimates out of the model. Results are demonstrated on a few qualitative setups, but most of the results are centered on improved JEM on sample quality, out-of-distribution detection, and semi-supervised learning. The main benefit of the approach seems to be speed and stability, however, the authors also claim that minimal tuning is needed.