GENERALIZED ENERGY BASED MODELS
Year: 2021

Abstract:
We introduce the Generalized Energy Based Model (GEBM) for generative modelling. These models combine two trained components: a base distribution (generally an implicit model), which can learn the support of data with low intrinsic dimension in a high dimensional space; and an energy function, to refine the probability mass on the learned support. Both the energy function and base jointly constitute the final model, unlike GANs, which retain only the base distribution (the "generator"). GEBMs are trained by alternating between learning the energy and the base. We show that both training stages are well-defined: the energy is learned by maximising a generalized likelihood, and the resulting energy-based loss provides informative gradients for learning the base. Samples from the posterior on the latent space of the trained model can be obtained via MCMC, thus finding regions in this space that produce better quality samples. Empirically, the GEBM samples on image-generation tasks are of much better quality than those from the learned generator alone, indicating that all else being equal, the GEBM will outperform a GAN of the same complexity. When using normalizing flows as base measures, GEBMs succeed on density modelling tasks, returning comparable performance to direct maximum likelihood of the same networks.

heading: 1 INTRODUCTION
text: Energy-based models (EBMs) have a long history in physics, statistics and machine learning (LeCun et al., 2006). They belong to the class of explicit models, and can be described by a family of energies E which define probability distributions with density proportional to exp(−E). Those models are often known up to a normalizing constantZ(E), also called the partition function. The learning task consists of finding an optimal function that best describes a given system or target distribution P. This can be achieved using maximum likelihood estimation (MLE), however the intractability of the normalizing partition function makes this learning task challenging. Thus, various methods have been proposed to address this (Hinton, 2002; Hyvärinen, 2005; Gutmann and Hyvärinen, 2012; Dai et al., 2019a;b). All these methods estimate EBMs that are supported over the whole space. In many applications, however, P is believed to be supported on an unknown lower dimensional manifold. This happens in particular when there are strong dependencies between variables in the data (Thiry et al., 2021), and suggests incorporating a low-dimensionality hypothesis in the model .
Generative Adversarial Networks (GANs) (Goodfellow et al., 2014) are a particular way to enforce low dimensional structure in a model. They rely on an implicit model, the generator, to produce samples supported on a low-dimensional manifold by mapping a pre-defined latent noise to the sample space using a trained function. GANs have been very successful in generating high-quality samples on various tasks, especially for unsupervised image generation (Brock et al., 2018). The generator is trained adversarially against a discriminator network whose goal is to distinguish samples produced by the generator from the target data. This has inspired further research to extend the training procedure to more general losses (Nowozin et al., 2016; Arjovsky et al., 2017; Li et al., 2017; Bińkowski et al., 2018; Arbel et al., 2018) and to improve its stability (Miyato et al., 2018; Gulrajani et al., 2017; Nagarajan and Kolter, 2017; Kodali et al., 2017). While the generator of a GAN has effectively a low-dimensional support, it remains challenging to refine the distribution of mass on that support using pre-defined latent noise. For instance, as shown by Cornish et al. (2020) for normalizing flows, when the latent distribution is unimodal and the target distribution possesses multiple disconnected low-dimensional components, the generator, as a continuous map, compensates for this mismatch using steeper slopes. In practice, this implies the need for more complicated generators.
∗Correspondence: michael.n.arbel@gmail.com.
In the present work, we propose a new class of models, called Generalized Energy Based Models (GEBMs), which can represent distributions supported on low-dimensional manifolds, while offering more flexibility in refining the mass on those manifolds. GEBMs combine the strength of both implicit and explicit models in two separate components: a base distribution (often chosen to be an implicit model) which learns the low-dimensional support of the data, and an energy function that can refine the probability mass on that learned support. We propose to train the GEBM by alternating between learning the energy and the base, analogous to f -GAN training (Goodfellow et al., 2014; Nowozin et al., 2016). The energy is learned by maximizing a generalized notion of likelihood which we relate to the Donsker-Varadhan lower-bound (Donsker and Varadhan, 1975) and Fenchel duality, as in (Nguyen et al., 2010; Nowozin et al., 2016). Although the partition function is intractable in general, we propose a method to learn it in an amortized fashion without introducing additional surrogate models, as done in variational inference (Kingma and Welling, 2014; Rezende et al., 2014) or by Dai et al. (2019a;b). The resulting maximum likelihood estimate, the KL Approximate Lower-bound Estimate (KALE), is then used as a loss for training the base. When the class of energies is rich and smooth enough, we show that KALE leads to a meaningful criterion for measuring weak convergence of probabilities. Following recent work by Chu et al. (2020); Sanjabi et al. (2018), we show that KALE possesses well defined gradients w.r.t. the parameters of the base, ensuring well-behaved training. We also provide convergence rates for the empirical estimator of KALE when the variational family is sufficiently well behaved, which may be of independent interest.
The main advantage of GEBMs becomes clear when sampling from these models: the posterior over the latents of the base distribution incorporates the learned energy, putting greater mass on regions in this latent space that lead to better quality samples. Sampling from the GEBM can thus be achieved by first sampling from the posterior distribution of the latents via MCMC in the low-dimensional latent space, then mapping those latents to the input space using the implicit map of the base. This is in contrast to standard GANs, where the latents of the base have a fixed distribution. We focus on a class of samplers that exploit gradient information, and show that these samplers enjoy fast convergence properties by leveraging the recent work of Eberle et al. (2017). While there has been recent interest in using the discriminator to improve the quality of the generator during sampling (Azadi et al., 2019; Turner et al., 2019; Neklyudov et al., 2019; Grover et al., 2019; Tanaka, 2019; Wu et al., 2019b), our approach emerges naturally from the model we consider.
We begin in Section 2 by introducing the GEBM model. In Section 3, we describe the learning procedure using KALE, then derive a method for sampling from the learned model in Section 4. In Section 5 we discuss related work. Finally, experimental results are presented in Section 6 with code available at https://github.com/MichaelArbel/GeneralizedEBM.

heading: 2 GENERALIZED ENERGY-BASED MODELS
text: In this section, we introduce generalized energy based models (GEBM), that combine the strengths of both energy-based models and implicit generative models, and admit the first of these as a special case. An energy-based model (EBM) is defined by a set E of real valued functions called energies, where eachE∈E specifies a probability density over the data spaceX ⊂Rd up to a normalizing constant,
Q(dx)=exp(−E(x)−A)dx, A=log Å∫ exp(−E(x))dx ã . (1)
While EBMs have been shown recently to be powerful models for representing complex high dimensional data distributions, they still unavoidably lead to a blurred model whenever data are concentrated on a lower-dimensional manifold. This is the case in Figure 1(a), where the ground truth distribution is
supported on a 1-D line and embedded in a 2-D space. The EBM in Figure 1(d) learns to give higher density to a halo surrounding the data, and thus provides a blurred representation. That is a consequence of EBM having a density defined over the whole space, and can result in blurred samples for image models.
An implicit generative model (IGM) is a family of probability distributions Gθ parametrized by a learnable generator function G :Z 7→X that maps latent samples z from a fixed latent distribution η to the data space X . The latent distribution η is required to have a density over the latent space Z and is often easy to sample from. Thus, Sampling from G is simply achieved by first sampling z from η then applyingG,
x∼G ⇐⇒ x=G(z), z∼η. (2)
GANs are popular instances of these models, and are trained adversarially (Goodfellow et al., 2014). When the latent spaceZ has a smaller dimension than the input spaceX , the IGM will be supported on a lower dimensional manifold of X , and thus will not possess a Lebesgue density on X (Bottou et al., 2017). IGMs are therefore good candidates for modelling low dimensional distributions. While GANs can accurately learn the low-dimensional support of the data, they can have limited power for representing the distribution of mass on the support. This is illustrated in Figure 1(b).
A generalized energy-based model (GEBM) Q is defined by a combination of a base G and an energyE defined over a subsetX of Rd. The base component can typically be chosen to be an IGM as in (2). The generalized energy component can refine the mass on the support defined by the base. It belongs to a class E of real valued functions defined on the input spaceX , and represents the negative log-density of a sample from the GEBM with respect to the base G,
Q(dx)=exp(−E(x)−AG,E)G(dx), AG,E=log Å∫ exp(−E(x))G(dx) ã , (3)
where AG,E is the logarithm of the normalizing constant of the model w.r.t. G. Thus, a GEBM Q re-weights samples from the base according to the un-normalized importance weights exp(−E(x)). Using the latent structure of the base G, this importance weight can be pulled-back to the latent space to define a posterior latent distribution ν,
ν(z) :=η(z)exp(−E(G(z))−AG,E). (4)
Hence, the posterior latent ν can be used instead of the latent noise η for sampling from Q, as summarized by Proposition 1: Proposition 1. Sampling from Q requires sampling a latent z from ν (4) then applying the mapG,
x∼Q ⇐⇒ x=G(z), z∼ν. (5)
In order to hold, Proposition 1 does not need the generatorG to be invertible. We provide a proof in Appendix C.1 which relies on a characterization of probability distribution using generalized moments. We will see later in Section 4 how equation (5) can be used to provide practical sampling algorithms from the GEBM. Next we discuss the advantages of GEBMs.
Advantages of Generalized Energy Based Models. The GEBM defined by (3) can be related to exponential tilting (re-weighting) (Siegmund, 1976; Xie et al., 2016) of the base G. The important difference over classical EBMs is that the base G is allowed to change its support and shape in space. By learning the base G, GEBMs can accurately learn the low-dimensional support of data, just like IGMs do. They also benefit from the flexibility of EBMs for representing densities using an energyE to refine distribution of mass on the support defined by G, as seen in Figure 1(c).
Compared to EBMs, that put mass on the whole space by construction (positive density), GEBMs have the additional flexibility to concentrate the probability mass on a low-dimensional support learned by the base G, provided that the dimension of the latent spaceZ is smaller than the dimension of the ambient space X : see Figure 1(c) vs Figure 1(d). In the particular case when the dimension of Z is equal to the ambient dimension andG is invertible, the baseG becomes supported over the whole space X , and GEBM recover usual EBMs. The next proposition further shows that any EBM can be viewed as a particular cases of GEBMs, as proved in Appendix C.1. Proposition 2. Any EBM with energyE (as in (1)) can be expressed as a GEBM with baseG given as a normalizing flow with density exp(−r(x)) and a generalized energy Ẽ(x)=E(x)−r(x). In this particular case, the dimension of the latent is necessarily equal to the data dimension, i.e. dim(Z)=dim(X ).
Compared to IGMs, that rely on a fixed pre-determined latent noise distribution η, GEBMs offer the additional flexibility of learning a richer latent noise distribution. This is particularly useful when the data is multimodal. In IGMs, such a GANs, the latent noise η is usually unimodal thus requiring a more sophisticated generator to distort a unimodal noise distribution into a distribution with multiple modes, as shown by Cornish et al. (2020). Instead, GEBMs allow to sample from a posterior ν over the latent noise defined in (4). This posterior noise can be multimodal in latent space (by incorporating information from the energy) and thus can put more or less mass in specific regions of the manifold defined by the base G. This allows GEBMs to capture multimodality in data, provided the support of the base is broad enough to subsume the data support Figure 1(c). The base can be simpler, compared to GANs, as it doesn’t need to distort the input noise too much to produce multimodal samples (see Figure 8 in Appendix G.4). This additional flexibility comes at no additional training cost compared to GANs. Indeed, GANs still require another model during training, the discriminator network, but do not use it for sampling. Instead, GEBMs avoid this waist since the base and energy can be trained jointly, with no other additional model, and then both are used for sampling.

heading: 3 LEARNING GEBMS
text: In this section we describe a general procedure for learning GEBMs. We decompose the learning procedure into two steps: an energy learning step and a base learning step. The overall learning procedure alternates between these two steps, as done in GAN training (Goodfellow et al., 2014).

heading: 3.1 ENERGY LEARNING
text: When the base G is fixed, varying the energy E leads to a family of models that all admit a density exp(−E−AG,E) w.r.t. G. When the base G admits a density exp(−r) defined over the whole space, it is possible to learn the energyE by maximizing the likelihood of the model− ∫ (E+r)dP−AG,E . However, in general G is supported on a lower-dimensional manifold so that r is ill-defined and the usual notion of likelihood cannot be used. Instead, we introduce a generalized notion of likelihood which does not require a well defined density exp(−r) for G: Definition 1 (Generalized Likelihood). The expected G-log-likelihood under a target distribution P of a GEBM model Q with base G and energyE is defined as
LP,G(E) :=− ∫ E(x)dP(x)−AG,E . (6)
To provide intuitions about the generalized likelihood in Definition 1, we start by discussing the particular case where KL(P||G) < +∞. We then present the training method in the general case where P and G might not share the same support, i.e. KL(P||G)=+∞. Special case of finite KL(P||G). When the Kullback-Leibler divergence between P and G is well defined, (6) corresponds to the Donsker-Varadhan (DV) lower bound on the KL (Donsker and Varadhan, 1975), meaning that KL(P||G)≥LP,G(E) for allE. Moreover, the following proposition holds: Proposition 3. Assume thatKL(P||G)<+∞ and 0∈E . If, in addition,E? maximizes (6), then:
KL(P||Q)≤KL(P||G). (7) In addition, we have thatKL(P||Q)=0 whenE? is the negative log-density ratio of P w.r.t. G.
We refer to Appendix C.1 for a proof. According to (7), the GEBM systematically improves over the IGM defined by G, with no further improvement possible in the limit case when G=P. Hence as long as there is an error in mass on the common support of P and G, the GEBM improves over the base G.
Estimating the likelihood in the General setting. Definition 1 can be used to learn a maximum likelihood energyE? by maximizingLP,G(E) w.r.t. E even when the KL(P||G) is infinite and when P and G don’t necessarily share the same support. Such an optimal solution is well defined whenever the set of energies is suitably constrained. This is the case if the energies are parametrized by a compact set Ψ with ψ 7→Eψ continuous over Ψ. Estimating the likelihood is then achieved using i.i.d. samples (Xn)1:N ,(Ym)1:M from P and G (Tsuboi et al., 2009; Sugiyama et al., 2012; Liu et al., 2017):
L̂P,G(E)=− 1
N N∑ n=1 E(Xn)−log
( 1
M M∑ m=1 exp(−E(Ym))
) . (8)
In the context of mini-batch stochastic gradient methods, however,M typically ranges from 10 to 1000, which can lead to a poor estimate for the log-partition functionAG,E . Moreover, (8) doesn’t exploit estimates ofAG,E from previous gradient iterations. Instead, we propose an estimator which introduces a variational parameter A ∈R meant to estimate AG,E in an amortized fashion. The key idea is to exploit the convexity of the exponential which directly implies−AG,E≥−A−exp(−A+AG,E)+1 for anyA∈R, with equality only whenA=AG,E . Therefore, (6) admits a lower-bound of the form
LP,G(E)≥− ∫ (E+A)dP− ∫ exp(−(E+A))dG+1:=FP,G(E+A),
where we introduced the functional FP,G for concision. Maximizing FP,G(E+A) over A recovers the likelihoodLP,G(E). Moreover, jointly maximizing overE andA yields the maximum likelihood energyE? and its corresponding log-partition functionA?=AG,E? . This optimization is well-suited for stochastic gradient methods using the following estimator Kanamori et al. (2011):
F̂P,G(E+A)=− 1
N N∑ n=1 (E(Xn)+A)− 1 M M∑ m=1 exp(−(E(Ym)+A))+1. (9)

heading: 3.2 BASE LEARNING
text: Unlike in Section 3.1, varying the base G does not need to preserve the same support. Thus, it is generally not possible to use maximum likelihood methods for learning G. Instead, we propose to use the generalized likelihood (6) evaluated at the optimal energyE? as a meaningful loss for learning G, and refer to it as the KL Approximate Lower-bound Estimate (KALE),
KALE(P||G)= sup (E,A)∈E×R FP,G(E+A). (10)
From Section 3.1, KALE(P||G) is always a lower bound on KL(P,G). The bound becomes tight whenever the negative log density of P w.r.t. G is well-defined and belongs to E (Appendix A). Moreover, Proposition 4 shows that KALE is a reliable criterion for measuring convergence, and is a consequence of (Zhang et al., 2017, Theorem B.1), with a proof in Appendix C.2.1: Proposition 4. Assume all energies in E are L-Lipschitz and that any continuous function can be well approximated by linear combinations of energies in E (Assumptions (A) and (B) of Appendix C.2), then KALE(P||G)≥0 with equality only if P=G and KALE(P||Gn)→0 iff Gn→P in distribution.
The universal approximation assumption holds in particular when E contains feedforward networks. In fact networks with a single neuron are enough, as shown in (Zhang et al., 2017, Theorem 2.3). The Lipschitz assumption holds when additional regularization of the energy is enforced during training by methods such as spectral normalization (Miyato et al., 2018) or additional regularization I(ψ) on the energyEψ such as the gradient penalty (Gulrajani et al., 2017) as done in Section 6.
Estimating KALE. According to Arora et al. (2017), accurate finite sample estimates of divergences that result from an optimization procedures (such as in (10)) depend on the richness of the class E ; and richer energy classes can result in slower convergence. Unlike divergences such as Jensen-Shannon, KL and the Wasserstein distance, which result from optimizing over a non-parametric and rich class of functions, KALE is restricted to a class of parametric energiesEψ . Thus, (Arora et al., 2017, Theorem 3.1) applies, and guarantees good finite sample estimates, provided optimization is solved accurately. In Appendix B, we provide an analysis for the more general case where energies are not necessarily parametric but satisfy some further smoothness properties; we emphasize that our rates do not require the strong assumption that the density ratio is bounded above and below as in (Nguyen et al., 2010).
Smoothness of KALE. Learning the base is achieved by minimizingK(θ) :=KALE(P||Gθ) over the set of parameters Θ of the generatorGθ using first order methods (Duchi et al., 2011; Kingma and Ba, 2014; Arbel et al., 2019). This requiresK(θ) to be smooth enough so that gradient methods converge to local minima and avoid instabilities during training (Chu et al., 2020). Ensuring smoothness of losses that result from an optimization procedure, as in (10), can be challenging. Results for the regularized Wasserstein are provided by Sanjabi et al. (2018), while more general losses are considered by Chu et al. (2020), albeit under stronger conditions than for our setting.Theorem 5 shows that whenE,Gθ and their gradients are all Lipschitz thenK(θ) is smooth enough.We provide a proof for Theorem 5 in Appendix C.2.1.
Theorem 5. Under Assumptions (I) to (III) of Appendix C.2, sub-gradient methods onK converge to local optima. Moreover,K is Lipschitz and differentiable for almost all θ∈Θ with:
∇K(θ)=exp(−AGθ,E?) ∫ ∇xE?(Gθ(z))∇θGθ(z)exp(−E?(Gθ(z)))η(z)dz. (11)
Estimating the gradient in (11) is achieved by first optimizing over Eψ andAusing (9), with additional regularization I(ψ). The resulting estimators Ê? and Â? are plugged in (12) to estimate ∇K(θ) using samples (Zm)1:M from η. Unlike for learning the energyE?, which benefits from using the amortized estimator of the log-partition function, we found that using the empirical log-partition for learning the base was more stable. We summarize the training procedure in Algorithm 1, which alternates between learning the energy and the base in a similar fashion to adversarial training. Algorithm 1 Training GEBM 1: Input P,N ,M , nb, ne 2: Output Trained generatorGθ and energyEψ . 3: Initialize θ , ψ andA. 4: for k=1,...,nb do 5: for j=1,...,ne do 6: Sample {Xn}1:N ∼P and {Yn}1:N ∼Gθ 7: gψ←−∇ψF̂P,Gθ (Eψ+A)+I(ψ) 8: Ã← log Ä 1 M ∑M m=1exp(−Eψ(Ym)) ä 9: gA←exp(A−Ã)−1
10: Update ψ andA using gψ and gA. 11: end for 12: Set Ê?←Eψ and Â?←A. 13: Update θ using ÷∇K(θ) from (12) 14: end for÷∇K(θ)= exp(−Â?)
M
M∑ m=1 ∇xÊ?(Gθ(Zm))∇θGθ(Zm)exp(−Ê?(Gθ(Zm))). (12)

heading: 4 SAMPLING FROM GEBMS
text: A simple estimate of the empirical distribution of observations under the GEBM is via importance sampling (IS). This consists in first sampling multiple points from the baseG, and then re-weighting the samples according to the energyE. Although straightforward, this approach can lead to highly unreliable estimates, a well known problem in the Sequential Monte Carlo (SMC) literature which employs IS extensively (Doucet et al., 2001; Del Moral et al., 2006). Other methods such as rejection sampling are known to be inefficient in high dimensions Haugh (2017). Instead, we propose to sample from the posterior ν using MCMC. Recall from (5) that a samplex fromQ is of the formx=G(z) with z sampled from the posterior latent ν of (4) instead of the priorη. While sampling fromη is often straightforward (for instance if η is a Gaussian), sampling from ν is generally harder, due to dependence of its density on complex functionsE andG. It is still possible to use MCMC methods to sample fromν, however, since we have access to its density up to a normalizing constant (4). In particular, we are interested in methods that exploit the gradient of ν, and consider two classes of samplers: Overdamped samplers and Kinetic samplers.
Overdamped samplers are obtained as a time-discretization of the Overdamped Langevin dynamics:
dzt=(∇zlogη(zt)−∇zE(G(zt)))+ √ 2dwt, (13)
wherewt is a standard Brownian motion. The simplest sampler arising from (13) is the Unadjusted Langevin Algorithm (ULA):
Zk+1 =Zk+λ(∇zlogη(Zk)−∇zE(G(Zk)))+ √ 2λWk+1, Z0∼η,
where (Wk)k≥0 are i.i.d. standard Gaussians and λ is the step-size. For large k,Zk is an approximate sample from ν (Raginsky et al., 2017, Proposition 3.3). Hence, settingX=G(Zk) for a large enough k provides an approximate sample from the GEBM Q, as summarized in Algorithm 2 of Appendix F.
Kinetic samplers arise from the Kinetic Langevin dynamics which introduce a momentum variable: dzt=vtdt, dvt=−γvtdt+u(∇logη(zt)−∇E(G(zt)))dt+ √ 2γudwt. (14)
with friction coefficient γ≥0, inverse massu≥0, momentum vector vt and standard Brownian motion wt. When the mass u−1 becomes negligible compared to the friction coefficient γ, i.e. uγ−2 ≈ 0, standard results show that (14) recovers the Overdamped dynamics (13). Discretization in time of (14)
leads to Kinetic samplers similar to Hamiltonian Monte Carlo (Cheng et al., 2017; Sachs et al., 2017). We consider a particular algorithm from Sachs et al. (2017) which we call Kinetic Langevin Algorithm (KLA) (see Algorithm 3 in Appendix F). Kinetic samplers were shown to better explore the modes of the invariant distribution ν compared to Overdamped ones (see (Neal, 2010; Betancourt et al., 2017) for empirical results and (Cheng et al., 2017) for theory), as also confirmed empirically in Appendix D for image generation tasks using GEBMs. Next, we provide the following convergence result:
Proposition 6. Assume that logη(z) is strongly concave and has a Lipschitz gradient, thatE,G and their gradients are allL-Lipschitz. Set xt=G(zt), where zt is given by (14) and call Pt the probability distribution of xt. Then Pt converges to Q in the Wasserstein sense,
W2(Pt,Q)≤LCe−cγt,
where c andC are positive constants independent of t, with c=O(exp(−dim(Z))).
Proposition 6 is proved in Appendix C.1 using (Eberle et al., 2017, Corollary 2.6), and implies that (xt)t≥0 converges at the same speed as (zt)t≥0. When the dimension q of Z is orders of magnitude smaller than the input space dimension d, the process (xt)t≥0 converges faster than typical sampling methods onX , for which the exponent controlling the convergence rate is of orderO(exp(−d)).

heading: 5 RELATED WORK
text: Energy based models. Usually, energy based models are required to have a density w.r.t. to a Lebesgue measure, and do not use a learnable base measure; in other words, models are supported on the whole space. Various methods have been proposed in the literature to learn EBMs. Contrastive Divergence (Hinton, 2002) approximates the gradient of the log-likelihood by sampling from the energy model with MCMC. More recently, (Belanger and McCallum, 2016; Xie et al., 2016; 2017; 2018c; 2019; Tu and Gimpel, 2018; Du and Mordatch, 2019; Deng et al., 2020) extend the idea using more sophisticated models and MCMC sampling strategies that lead to higher quality estimators. Score Matching (Hyvärinen, 2005) calculates an alternative objective (the score) to the log-likelihood which is independent of the partition function, and was recently used in the context non-parametric energy functions to provide estimators of the energy that are provably consistent (Sriperumbudur et al., 2017; Sutherland et al., 2018; Arbel and Gretton, 2018; Wenliang et al., 2019). In Noise-Contrastive Estimation (Gutmann and Hyvärinen, 2012), a classifier is trained to distinguish between samples from a fixed proposal distribution and the target P. This provides an estimate for the density ratio between the optimal energy model and the proposal distribution. In a similar spirit, Cranmer et al. (2016) uses a classifier to learn likelihood ratios. Conversely, Grathwohl et al. (2020) interprets the logits of a classifier as an energy model obtained after marginalization over the classes. The resulting model is then trained using Contrastive Divergence. In more recent work, Dai et al. (2019a;b) exploit a dual formulation of the logarithm of the partition function as a supremum over the set of all probability distributions of some functional objective. Yu et al. (2020) explore methods for using general f-divergences, such as Jensen-Shannon, to train EBMs.
Generative Adversarial Networks. Recent work proposes using the discriminator of a trained GAN to improve the generator quality. Rejection sampling (Azadi et al., 2019) and MetropolisHastings correction (Turner et al., 2019; Neklyudov et al., 2019) perform sampling directly on the high-dimensional input space without using gradient information provided by the discriminator. Moreover, the data distribution is assumed to admit a density w.r.t. the generator. Ding et al. (2019) perform sampling on the feature space of some auxiliary pre-trained network; while Lawson et al. (2019) treat the sampling procedure as a model on its own, learned by maximizing the ELBO. In our case, no auxiliary model is needed. In the present work, sampling doesn’t interfere with training, in contrast to recently considered methods to optimize over the latent space during training Wu et al. (2019b;a). In Tanaka (2019), the discriminator is viewed as an optimal transport map between the generator and the data distribution and is used to compute optimized samples from latent space. This is in contrast to the diffusion-based sampling that we consider. In (Xie et al., 2018b;a), two independent models, a full support EBM and a generator network, are trained cooperatively using MCMC. By contrast, in the present work, the energy and base are part of the same model, and the model support is lower-dimensional than the target space X . While we do not address the mode collapse problem, Xu et al. (2018); Nguyen et al. (2017) showed that KL-based losses are resilient to it thanks to the zero-avoiding property of the KL, a good sign for KALE which is derived from KL by Fenchel duality.
The closest related approach appears in a study concurrent to the present work (Che et al., 2020), where the authors propose to use Langevin dynamics on the latent space of a GAN generator, but with a different discriminator to ours (derived from the Jensen-Shannon divergence or a Wasserstein-based divergence). Our theory results showing the existence of the loss gradient (Theorem 5), establishing weak convergence of distributions under KALE (Proposition 4), and demonstrating consistency of the KALE estimator (Appendix B) should transfer to the JS and Wasserstein criteria used in that work. Subsequent to the present work, an alternative approach has been recently proposed, based on normalising flows, to learn both the low-dimensional support of the data and the density on this support (Brehmer and Cranmer, 2020). This approach maximises the explicit likelihood of a data projection onto a learned manifold, and may be considered complementary to our approach.

heading: 6 EXPERIMENTS
text: We train the models for 150000 generator iterations using Algorithm 1. After training is completed, we rescale the energy by β=100 to get a colder version of the GEBM and sample from it using either Algorithm 2 (ULA) or Algorithm 3 (KLA) with parameters (γ=100,u=1). This colder temperature leads to an improved FID score, and needs relatively few MCMC iterations, as shown in Figure 6 of Appendix D. Sampler convergence to visually plausible modes at low tempteratures is demonstrated in Figure 2. We perform 1000 MCMC iterations with initial step-size of λ=10−4 decreased by 10 every 200 iterations. As a baseline we consider samples generated from the base of the GEBM only (without using information from the energy) and call this KALE-GAN. More details are given in Appendix G.
Results: Table 1 shows that GEBM outperforms both KALE and standard GANs when using the same networks for the base/generator and energy/critic. Moreover, KALE-GAN matches the performance of a standard GAN (with Jensen-Shannon critic), showing that the improvement of GEBM cannot be explained by the switch from Jensen-Shannon to a KALE-based critic. Rather, the improvement is largely due to incorporating the energy function into the model, and sampling using Algorithm 3.
This finding experimentally validates our claim that incorporating the energy improves the model, and that all else being equal, a GEBM outperforms a GAN with the same generator and critic architecture. Indeed, if the critic is not zero at convergence, then by definition it contains information on the remaining mismatch between the generator (base) and data mass, which the GEBM incorporates, but the GAN does not. The GEBM also outperforms an EBM even when the latter was trained using a larger network (ResNet) with supervision (S) on ImageNet, which is an easier task ( Chen et al. (2019)). More comparisons on Cifar10 and ImageNet are provided in Table 4 of Appendix D.
Table 2 shows different sampling methods using the same trained networks (generator and critic), with KALE-GAN as a baseline. All energy-exploiting methods outperform the unmodified KALE-GAN with the same architecture. That said, our method (both ULA and KLA) outperforms both (IHM) (Turner et al., 2019) and (DOT) (Tanaka, 2019), which both use the energy information.
Cifar10 LSUN CelebA ImageNet
KALE-GAN 32.03 21.67 6.91 19.37 IHM 30.47 20.63 6.39 18.15 DOT 26.35 20.41 5.93 16.21 GEBM (ULA) 23.02 16.23 5.21 14.00 GEBM (KLA) 24.29 15.25 5.38 13.94
Table 2: FID scores for different sampling methods using the same trained SNGAN (ConvNet): KALE-GAN as a baseline w/o critic information.
In Table 2, KLA was used in the high friction regime γ=100 and thus behaves like ULA. This allows to obtain sharper samples concentrated around the modes of the GEBM thus improving the FID score. If, instead, the goal is to encourage more exploration of the modes of the GEBM, then KLA with a smaller γ is a better alternative than ULA, as the former can explore multiple modes/images within the same MCMC chain, unlike (ULA): see Figures 3 to 5 of Appendix D. Moving from one mode to another results in an increased FID score while between modes, however, which can be avoided by decreasing λ.

heading: 6.2 DENSITY ESTIMATION
text: Motivation. We next consider the particular setting where the likelihood of the model is well-defined, and admits a closed form expression. This is intended principally as a sanity check that our proposed training method in Algorithm 1 succeeds in learning maximum likelihood solutions. Outside of this setting, closed form expressions of the normalizing constant are not available for generic GEBMs. While this is not an issue (since the proposed method doesn’t require a closed form expression for the normalizing constant), in this experiment only, we want to have access to closed form expressions, as they enable a direct comparison with other density estimation methods.
Experimental setting. To have a closed-form likelihood, we consider the case where the dimension of the latent space is equal to data-dimension, and choose the baseG of the GEBM to be a Real NVP (Ding et al. (2019) ) with density exp(−r(x)) and energyE(x)=h(x)−r(x). Thus, in this particular case, the GEBM has a well defined likelihood over the whole space, and we are precisely in the setting of Proposition 2, which shows that this GEBM is equal to an EBM with density proportional to exp(−h). We further require the EBM to be a second Real NVP so that its density has a closed form expression. We consider 5 UCI datasets (Dheeru and Taniskidou, 2017) for which we use the same pre-processing as in (Wenliang et al., 2019). For comparison, we train the EBM by direct maximum likelihood (ML) and contrastive divergence (CD). To train the GEBM, we use Algorithm 1, which doesn’t directly exploit the closed-form expression of the likelihood (unlike direct ML). We thus use either (8) (KALE-DV) or (9) (KALE-F) to estimate the normalizing constant. More details are given in Appendix G.2.
Results. Table 3 reports the Negative Log-Likelihood (NLL) evaluated on the test set and corresponding to the best performance on the validation set. Training the GEBM using Algorithm 1 leads to comparable performance to (CD) and (ML). As shown in Figure 7 of Appendix E, (KALE-DV) and (KALEF) maintain a small error gap between the training and test NLL and, as discussed in Section 3.1 and Appendix F, (KALE-F) leads to more accurate estimates of the log-partition function, with a relative error of order 0.1% compared to 10% for (KALE-DV).

heading: 7 ACKNOWLEDGMENTS
text: We thank Mihaela Rosca for insightful discussions and Song Liu, Bo Dai and Hanjun Dai for pointing us to important related work.

heading: A KL APPROXIMATE LOWER-BOUND ESTIMATE
text: We discuss the relation between KALE (10) and the Kullback-Leibler divergence via Fenchel duality. Recall that a distribution P is said to admit a density w.r.t. G if there exists a real-valued measurable function r0 that is integrable w.r.t. G and satisfies dP = r0dG. Such a density is also called the Radon-Nikodym derivative of P w.r.t. G. In this case, we have:
KL(P||G)= ∫ r0log(r0)dG. (15)
Nguyen et al. (2010); Nowozin et al. (2016) derived a variational formulation for the KL using Fenchel duality. By the duality theorem (Rockafellar, 1970), the convex and lower semi-continuous function ζ :u 7→ulog(u) that appears in (15) can be expressed as the supremum of a concave function:
ζ(u)=sup v uv−ζ?(v).
The function ζ? is called the Fenchel dual and is defined as ζ?(v) = supuuv−ζ(u). By convention, the value of the objective is set to −∞ whenever u is outside of the domain of definition of ζ?. When ζ(u) = u log(u), the Fenchel dual ζ?(v) admits a closed form expression of the form ζ?(v)=exp(v−1). Using the expression of ζ in terms of its Fenchel dual ζ?, it is possible to express KL(P||G) as the supremum of the variational objective (16) over all measurable functions h.
F(h) :=− ∫ hdP− ∫ exp(−h)dG+1. (16)
Nguyen et al. (2010) provided the variational formulation for the reverse KL using a different choice for ζ: (ζ(u) =−log(u)). We refer to (Nowozin et al., 2016) for general f -divergences. Choosing a smaller set of functionsH in the variational objective (16) will lead to a lower bound on the KL. This is the KL Approximate Lower-bound Estimate (KALE):
KALE(P||G)= sup h∈H F(h) (17)
In general, KL(P||G) ≥ KALE(P||G). The bound is tight whenever the negative log-density h0 =−logr0 belongs toH; however, we do not require r0 to be well-defined in general. Equation (17) has the advantage that it can be estimated using samples from P and G. Given i.i.d. samples (X1, ...,XN ) and (Y1, ...,YM ) from P and G, we denote by P̂ and Ĝ the corresponding empirical distributions. A simple approach to estimate KALE(P||G) is to use anM -estimator. This is achieved by optimizing the penalized objective
ĥ :=argmax h∈H “F(h)− λ 2 I2(h), (18)
where “F is an empirical version ofF and I2(h) is a penalty term that prevents overfitting due to finite samples. The penalty I2(h) acts as a regularizer favoring smoother solutions while the parameter λ determines the strength of the smoothing and is chosen to decrease as the sample size N and M increase. TheM -estimator of KALE(P||G) is obtained simply by plugging in ĥ into the empirical objective “F(h): ÷KALE(P||G) := “F(ĥ). (19) We defer the consistency analysis of (19) to Appendix B where we provide convergence rates in a setting where the set of functions H is a Reproducing Kernel Hilbert Space and under weaker assumptions that were not covered by the framework of Nguyen et al. (2010).

heading: B CONVERGENCE RATES OF KALE
text: In this section, we provide a convergence rate for the estimator in (19) whenH is an RKHS. The theory remains the same whetherH contains constants or not. With this choice, the Representer Theorem allows us to reduce the potentially infinite-dimensional optimization problem in (18) to a convex finite-dimensional one. We further restrict ourselves to the well-specified case where the density r0 of P w.r.t. G is well-defined and belongs toH, so that KALE matches the KL. While Nguyen et al.
(2010) (Theorem 3) provides a convergence rate of 1/ √ N for a related M -estimator, this requires the density r0 to be lower-bounded by 0 as well as (generally) upper-bounded. This can be quite restrictive if, for instance, r0 is the density ratio of two gaussians. In Theorem 7, we provide a similar convergence rate for the estimator defined in (19) without requiring r0 to be bounded. We start by briefly introducing some notations, the working assumptions and the statement of the convergence result in Appendix B.1 and provide the proofs in Appendix B.2.

heading: B.1 STATEMENT OF THE RESULT
text: We recall that an RKHSH of functions defined on a domain X ⊂Rd and with kernel k is a Hilbert space with dot product 〈.,.〉, such that y 7→k(x,y) belongs toH for any x∈X , and
k(x,y)=〈k(x,.),k(y,.)〉, ∀x,y∈X .
Any function h inH satisfies the reproducing property f(x)=〈f,k(x,.)〉 for any x∈X . Recall that KALE(P||G) is obtained as an optimization problem
KALE(P||G)= sup h∈H F(h) (20)
whereF is given by:
F(h) :=− ∫ hdP− ∫ exp(−h)dG+1.
Since the negative log density ratio h0 is assumed to belong to H, this directly implies that the supremum of F is achieved at h0 and F(h0) = KALE(P||G). We are interested in estimating KALE(P||G) using the empirical distributions P̂ and Ĝ,
P̂ := 1
N N∑ n=1 δXn , Ĝ := 1 N N∑ n=1 δYn ,
where (Xn)1≤n≤N and (Yn)1≤n≤N are i.i.d. samples from P and G. For this purpose we introduce the empirical objective functional,“F(h) :=−∫ hdP̂−∫ exp(−h)dĜ+1. The proposed estimator is obtained by solving a regularized empirical problem,
sup h∈H “F(h)− λ 2 ‖h‖2, (21)
with a corresponding population version,
sup h∈H F(h)− λ 2 ‖h‖2. (22)
Finally, we introduceD(h,δ) and Γ(h,δ):
D(h,δ)= ∫ δexp(−h)dG− ∫ δdP,
Γ(h,δ)=− ∫ ∫ 1
0
(1−t)δ2exp(−(h+tδ))dG.
The empirical versions of D(h,δ) and Γ(h,δ) are denoted D̂(h,δ) and Γ̂(h,δ). Later, we will show thatD(h,δ) D̂(h,δ) are in fact the gradients ofF(h) and “F(h) along the direction δ. We state now the working assumptions:
(i) The supremum ofF overH is attained at h0.
(ii) The following quantities are finite for some positive :∫ » k(x,x) dP(x),∫ » k(x,x)exp((‖h0‖+ ) » k(x,x)) dG(x),∫
k(x,x)exp((‖h0‖+ ) » k(x,x)) dG(x).
(iii) For any h∈H, ifD(h,δ)=0 for all δ then h=h0. Theorem 7. Fix any 1>η > 0. Under Assumptions (i) to (iii), and provided that λ= 1√
N , it holds
with probability at least 1−2η that |“F(ĥ)−F(h0)|≤M ′(η,h0)√ N
for a constantM ′(η,h0) that depends only on η and h0.
The assumptions in Theorem 7 essentially state that the kernel associated to the RKHSH needs to satisfy some integrability requirements. That is to guarantee that the gradient δ 7→∇F(h)(δ) and its empirical version are well-defined and continuous. In addition, the optimality condition∇F(h)=0 is assumed to characterize the global solution h0. This will be the case if the kernel is characteristic Simon-Gabriel and Scholkopf (2018). The proof of Theorem 7, in Appendix B.2, takes advantage of the Hilbert structure of the setH, the convexity of the functionalF and the optimality condition∇“F(ĥ)=λĥ of the regularized problem, all of which turn out to be sufficient for controlling the error of (19).

heading: B.2 PROOFS
text: We state now the proof of Theorem 7 with subsequent lemmas and propositions.
Proof of Theorem 7. We begin with the following inequalities:
λ 2 (‖ĥ‖2−‖h0‖2)≤ “F(ĥ)−“F(h0)≤〈∇“F(h0),ĥ−h0〉.
The first inequality is by definition of ĥ while the second is obtained by concavity of “F . For simplicity we write B=‖ĥ−h0‖ and C=‖∇“F(h0)−L(h0)‖. Using Cauchy-Schwarz and triangular inequalities, it is easy to see that
−λ 2
( B2+2B‖h0‖ ) ≤ “F(ĥ)−“F(h0)≤CB.
Moreover, by triangular inequality, it holds that
B≤‖hλ−h0‖+‖ĥ−hλ‖. Lemma 11 ensures thatA(λ)=‖hλ−h0‖ converges to 0 as λ→0. Furthermore, by Proposition 12, we have ‖ĥ−hλ‖≤ 1λD whereD(λ)=‖∇“F(hλ)−∇L(hλ)‖. Now choosing λ= 1√N and applying Chebychev inequality in Lemma 8, it follows that for any 1>η>0,we have with probability greater than 1−2η that both
D(λ)≤ C(‖h0‖η)√ N , C≤ C(‖h0‖,η)√ N ,
where C(‖h0‖,η) is defined in Lemma 8. This allows to conclude that for any η > 0, it holds with probability at least 1−2η that |“F(ĥ)−“F(h0)|≤ M ′(η,h0)√N whereM ′(η,h0) depends only on η and h0. We proceed using the following lemma, which provides an expression forD(h,δ) and D̂(h,δ) along with a probabilistic bound:
Lemma 8. Under Assumptions (i) and (ii), for any h∈H such that ‖h‖≤‖h0‖+ , there existsD(h) inH satisfying
D(h,δ)=〈δ,D(h)〉, and for any h∈H, there exists “D(h) satisfying“D(h,δ)=〈δ,“D(h)〉. Moreover, for any 0<η<1 and any h∈H such that ‖h‖≤‖h0‖+ :=M , it holds with probability greater than 1−η that
‖D(h)−“D(h)‖≤ C(M,η)√ N ,
whereC(M,η) depends only onM and η.
Proof. First, we show that δ 7→D(h,δ) is a bounded linear operator. Indeed, Assumption (ii) ensures that k(x,.) and k(x,.)exp(−h(x)) are Bochner integrable w.r.t. P and G (Retherford (1978)), hence D(h,δ) is obtained as
D(h,δ) :=〈δ,µexp(−h)G−µP〉, where µexp(−h)G = ∫ k(x,.)exp(−h(x))dG and µP = ∫ k(x,.)dP. DefiningD(h) to be =µexp(−h)G− µP leads to the desired result. “D(h) is simply obtained by taking the empirical version ofD(h). Finally, the probabilistic inequality is a simple consequence of Chebychev’s inequality.
The next lemma states thatF(h) and “F(h) are Frechet differentiable. Lemma 9. Under Assumptions (i) and (ii) , h 7→F(h) is Frechet differentiable on the open ball of radius ‖h0‖+ while h 7→ “F(h) is Frechet differentiable onH. Their gradients are given by D(h) and “D(h) as defined in Lemma 8,
∇F(h)=D(h), ∇“F(h)=“D(h) Proof. The empirical functional “F(h) is differentiable since it is a finite sum of differentiable functions, and its gradient is simply given by “D(h). For the population functional, we use second order Taylor expansion of exp with integral remainder, which gives
F(h+δ)=F(h)−D(h,δ)+Γ(h,δ).
By Assumption (ii) we know that Γ(h,δ)‖δ‖ converges to 0 as soon as ‖δ‖→ 0. This allows to directly conclude that F is Frechet differentiable, with differential given by δ 7→D(h,δ). By Lemma 8, we conclude the existence of a gradient∇F(h) which is in fact given by∇F(h)=D(h). From now on, we will only use the notation∇F(h) and∇“F(h) to refer to the gradients ofF(h) and“F(h). The following lemma states that (21) and (22) have a unique global optimum, and gives a first order optimality condition.
Lemma 10. The problems (21) and (22) admit unique global solutions ĥ and hλ inH. Moreover, the following first order optimality conditions hold: λĥ=∇“F(ĥ), λhλ=∇F(hλ). Proof. For (21), existence and uniqueness of a minimizer ĥ is a simple consequence of continuity and strong concavity of the regularized objective. We now show the existence result for (22). Let’s introduce Gλ(h)=−F(h)+ λ2 ‖h‖
2 for simplicity. Uniqueness is a consequence of the strong convexity of Gλ. For the existence, consider a sequence of elements fk ∈H such that Gλ(fk)→ infh∈HGλ(h). If h0
is not the global solution, then it must hold for k large enough that Gλ(fk)≤Gλ(h0). We also know that F(fk)≤F(h0), hence, it is easy to see that ‖fk‖≤‖h0‖ for k large enough. This implies that fk is a bounded sequence, therefore it admits a weakly convergent sub-sequence by weak compactness. Without loss of generality we assume that fk weakly converges to some element hλ ∈ H and that ‖fk‖≤‖h0‖. Hence, ‖hλ‖≤ liminfk‖fk‖≤‖h0‖. Recall now that by definition of weak convergence, we have fk(x)→k hλ(x) for all x∈X . By Assumption (ii), we can apply the dominated convergence theorem to ensure thatF(fk)→F(hλ). Taking the limit of Gλfk, the following inequality holds:
sup h∈H Gλ(h)=limsup k Gλ(fk)≤Gλ(hλ).
Finally, by Lemma 9 we know thatF is Frechet differentiable, hence we can use Ekeland and Témam (1999) (Proposition 2.1) to conclude that∇F(hλ) = λhλ. We use exactly the same arguments for (21).
Next, we show that hλ converges towards h0 inH. Lemma 11. Under Assumptions (i) to (iii) it holds that:
A(λ) :=‖hλ−h0‖→0.
Proof. We will first prove that hλ converges weakly towards h0, and then conclude that it must also converge strongly. We start with the following inequalities:
0≥F(hλ)−F(h0)≥ λ
2 (‖hλ‖2−‖h0‖2).
These are simple consequences of the definitions of hλ and h0 as optimal solutions to (20) and (21). This implies that ‖hλ‖ is always bounded by ‖h0‖. Consider now an arbitrary sequence (λm)m≥0 converging to 0. Since ‖hλm‖ is bounded by ‖h0‖, it follows by weak-compactness of balls inH that hλm admits a weakly convergent sub-sequence. Without loss of generality we can assume that hλm is itself weakly converging towards an element h∗. We will show now that h∗ must be equal to h0. Indeed, by optimality of hλm , it must hold that
λmhλm =∇F(hm).
This implies that ∇F(hm) converges weakly to 0. On the other hand, by Assumption (ii), we can conclude that ∇F(hm) must also converge weakly towards ∇F(h∗), hence ∇F(h∗) = 0. Finally by Assumption (iii) we know that h0 is the unique solution to the equation∇F(h)=0 , hence h∗=h0. We have shown so far that any subsequence of hλm that converges weakly, must converge weakly towards h0. This allows to conclude that hλm actually converges weakly towards h0. Moreover, we also have by definition of weak convergence that:
‖h0‖≤ lim inf m→∞ ‖hλm‖.
Recalling now that ‖hλm‖≤ ‖h0‖ it follows that ‖hλm‖ converges towards ‖h0‖. Hence, we have the following two properties:
• hλm converges weakly towards h0,
• ‖hλm‖ converges towards ‖h0‖.
This allows to directly conclude that ‖hλm−h0‖ converges to 0.
Proposition 12. We have that:
‖ĥ−hλ‖≤ 1
λ ‖∇F̂(hλ)−∇F(hλ)‖
Proof. By definition of ĥ and hλ the following optimality conditions hold: λĥ=∇“F(ĥ), λhλ=∇F(hλ).
We can then simply write: λ(ĥ−hλ)−(∇“F(ĥ)−∇“F(hλ))=∇“F(hλ)−∇F(hλ). Now introducing δ := ĥ−hλ andE :=∇“F(ĥ)−∇“F(hλ) for simplicity and taking the squared norm of the above equation, it follows that
λ2‖δ‖2+‖E‖2−2λ〈δ,E〉=‖∇“F(hλ)−∇F(hλ)‖2. By concavity of “F onHwe know that−〈ĥ−hλ,E〉≥0. Therefore:
λ2‖ĥ−hλ‖2≤‖∇“F(hλ)−∇F(hλ)‖2.

heading: C LATENT NOISE SAMPLING AND SMOOTHNESS OF KALE

heading: C.1 LATENT SPACE SAMPLING
text: Here we prove Proposition 6 for which we make the assumptions more precise: Assumption 1. We make the following assumption:
• logη is strongly concave and admits a Lipschitz gradient.
• There exists a non-negative constantL such that for any x,x′∈X and z,z′∈Z:
|E(x)−E(x′)|≤‖x−x′‖, ‖∇xE(x)−∇xE(x′)‖≤‖x−x′‖ |G(z)−G(z′)|≤‖z−z′‖, ‖∇zG(z)−∇zG(z′)‖≤‖z−z′‖
Throughout this section, we introduceU(z) :=−log(η(z))+E(G(z)) for simplicity.
Proof of Proposition 1 . To sample fromQG,E , we first need to identify the posterior latent distribution νG,E used to produce those samples. We rely on (23) which holds by definition of QG,E for any test function h onX : ∫
h(x)dQ(x)= ∫ h(G(z))f(G(z))η(z)dz, (23)
Hence, the posterior latent distribution is given by ν(z) = η(z)f(G(z)), and samples from GEBM are produced by first sampling from νG,E , then applying the implicit mapG,
X∼Q ⇐⇒ X=G(Z), Z∼ν.
Proof of Proposition 2. the base distribution G admits a density on the whole space denoted by exp(−r(x)) and the energy Ẽ is of the form Ẽ(x) =E(x)−r(x) for some parametric function E, it is easy to see that Q has a density proportional to exp(−E) and is therefore equivalent to a standard EBM with energyE.
The converse holds as well, meaning that for any EBM with energy E, it is possible to construct a GEBM using an importance weighting strategy. This is achieved by first choosing a base G, which is required to have an explicit density exp(−r) up to a normalizing constant, then defining the energy of the GEBM to be Ẽ(x)=E(x)−r(x) so that:
dQ(x)∝exp(−Ẽ(x))dGθ(x)∝exp(−E(x))dx (24)
Equation (24) effectively depends only on E(x) and not on G since the factor exp(r) exactly compensates for the density of G. The requirement that the base also admits a tractable implicit map G can be met by choosing G to be a normalizing flow (Rezende and Mohamed, 2015) and does not restrict the class of possible EBMs that can be expressed as GEBMs.
Proof of Proposition 6. Let πt be the probability distribution of (zt,vt) at time t of the diffusion in (14), which we recall that
dzt=vtdt, dvt =−(γvt+u∇U(zt))+ √ 2λudwt,
We call π∞ its corresponding invariant distribution given by π∞(z,v)∝exp Å −U(z)− 1 2 ‖v‖2 ã
By Lemma 13 we know thatU is dissipative, bounded from below, and has a Lipschitz gradient. This allows to directly apply (Eberle et al., 2017)(Corollary 2.6.) which implies that
W2(πt,π∞)≤Cexp(−tc),
where c is a positive constant and C only depends on π∞ and the initial distribution π0. Moreover, the constant c is given explicitly in (Eberle et al., 2017, Theorem 2.3) and is of order 0(e−q) where q is the dimension of the latent spaceZ . We now consider an optimal coupling Πt between πt and π0. Given joints samples ((zt,vt),(z,v)) from Πt, we consider the following samples in input space (xt,x) := (G(zt),G(z)). Since zt and z have marginals πt and π∞, it is easy to see that xt ∼ Pt and x∼Q. Therefore, by definition of the W2 distance, we have the following bound:
W 22 (Pt,Q)≤E [ ‖xt−x‖2 ] ≤ ∫ ‖G(zt)−G(z)‖2dΠt(zt,z)
≤L2 ∫ ‖zt−z‖2dΠt(zt,z)
≤L2W 22 (πt,π∞)≤C2L2exp(−2tc).
The second line uses the definition of (xt,x) as joint samples obtained by mapping (zt,z). The third line uses the assumption thatB isL-Lipschitz. Finally, the last line uses that Πt is an optimal coupling between πt and π∞.
Lemma 13. Under Assumption 1, there existsA>0 and λ∈(0, 14 ] such that
1 2 z>t∇U(z)≥λ
Å U(z)+ γ2 4u ‖z‖2 ã −A, ∀z∈Z, (25)
where γ and u are the coefficients appearing in (14). Moreover, U is bounded bellow and has a Lipschitz gradient.
Proof. For simplicity, let’s call u(z) = − log η(z), w(z) = E? ◦ Bθ?(z), and denote by M an upper-bound on the Lipschitz constant ofw and∇w which is guaranteed to be finite by assumption. HenceU(z)=u(z)+w(z). Equation (25) is equivalent to having
z>∇u(z)−2λu(z)− γ 2
2u ‖z‖2≥2λw(z)−z>∇w(z)−2A. (26)
Using that w is Lipschitz, we have that w(z) ≤ w(0) +M‖z‖ and −z>∇w(z) ≤M‖z‖. Hence, 2λw(z)−z>∇w(z)−2A≤2λw(0)+(2λ+1)M‖z‖−2A. Therefore, a sufficient condition for (26) to hold is
z>∇u(z)−2λu(z)− γ 2
2u ‖z‖2≥+(2λ+1)M‖z‖−2A+2λw(0). (27)
We will now rely on the strong convexity of u,which holds by assumption, and implies the existence of a positive constantm>0 such that
−u(z)≥−u(0)−z>∇u(z)+m 2 ‖z‖2,
z>∇u(z)≥−‖z‖‖∇u(0)‖+m‖z‖2.
This allows to write the following inequality,
z>∇u(z)−2λu(z)− γ 2 2u ≥(1−2λ)z>∇u(z)+λ(m+ γ 2 2u )‖z‖2−2λu(0)
≥(1−λ(m+ γ 2
2u ))‖z‖2−(1−2λ)‖z‖‖∇u(0)‖−2λu(0).
Combining the previous inequality with (27) and denotingM ′=‖∇u(0)‖ , it is sufficient to findA and λ satisfyingÅ 1−λ Å m+ γ2
2u
ãã ‖z‖2−(M+M ′+2λ(M−M ′))‖z‖−2λ(u(0)+w(0))+2A≥0.
The l.h.s. in the above equation is a quadratic function in ‖z‖ and admits a global minimum when λ< Ä m+ γ 2
2u
ä−1 . The global minimum is always positive provided thatA is large enough.
To see that U is bounded below, it suffice to note, by Lipschitzness of w, that w(z)≥w(0)−M‖z‖ and by strong convexity of u that
u(z)≥u(0)+M ′‖z‖+m 2 ‖z‖2.
Hence,U is lower-bounded by a quadratic function in ‖z‖with positive leading coefficient m2 , hence it must be lower-bounded by a constant. Finally, by assumption, u and w have Lipschitz gradients, which directly implies thatU has a Lipschitz gradient.
Proof of Proposition 3. By assumption KL(P||G)<+∞, this implies that P admits a density w.r.t. G which we call r(x). As a result P admits also a density w.r.t. Q given by:
Zexp(E?(x))r(x).
We can then compute theKL(P||Q) explicitly:
KL(P||Q)=EP[E]+log(Z)+EP[log(r)] =−LP,G(E?)+KL(P||G).
Since 0 belongs to E and by optimality ofE?, we know thatLP,G(E?)≥LP,G(0)=0. The result then follows directly.

heading: C.2 TOPOLOGICAL AND SMOOTHNESS PROPERTIES OF KALE
text: Topological properties of KALE. Denseness and smoothness of the energy class E are the key to guarantee that KALE is a reliable criterion for measuring convergence. We thus make the following assumptions on E :
(A) For all E ∈ E , −E ∈ E and there is CE > 0 such that cE ∈ E for 0 ≤ c ≤ CE . For any continuous function g, any compact support K in X and any precision > 0, there exists a finite linear combination of energiesG= ∑r i=1aiEi such that supx∈K |f(x)−G(x)|≤ .
(B) All energiesE in E are Lipschitz in their input with the same Lipschitz constantL>0.
Assumption (A) holds in particular when E contains feedforward networks with a given number of parameters. In fact networks with a single neuron are enough, as shown in (Zhang et al., 2017, Theorem 2.3). Assumption (B) holds when additional regularization of the energy is enforced during training by methods such as spectral normalization Miyato et al. (2018) or gradient penalty Gulrajani et al. (2017) as done in Section 6. Proposition 4 states the topological properties of KALE ensuring that it can be used as a criterion for weak convergence. A proof is given in Appendix C.2.1 and is a consequence of (Zhang et al., 2017, Theorem B.1). Proposition 14. Under Assumptions (A) and (B) it holds that:
1. KALE(P||G)≥0 with KALE(P||G)=0 if and only if P=G.
2. KALE(P||Gn)→0 if and only if Gn→P under the weak topology.

heading: C.2.1 TOPOLOGICAL PROPERTIES OF KALE
text: In this section we prove Proposition 4. We first start by recalling the required assumptions and make them more precise:
Assumption 2. Assume the following holds:
• The setX is compact.
• For all E ∈ E , −E ∈ E and there is CE > 0 such that cE ∈ E for 0 ≤ c ≤ CE . For any continuous function g, any compact support K in X and any precision > 0, there exists a finite linear combination of energiesG= ∑r i=1aiEi such that |f(x)−G(x)|≤ onK.
• All energiesE in E are Lipschitz in their input with the same Lipschitz constantL>0.
For simplicity we consider the setH=E+R, i.e.: H is the set of functions h of the form h=E+c where E ∈E and c∈R. In all what follows P1 is the set of probability distributions with finite first order moments. We consider the notion of weak convergence on P1 as defined in (Villani, 2009, Definition 6.8) which is equivalent to convergence in the Wasserstein-1 distanceW1.
Proof of Proposition 4 . We proceed by proving the separation properties (1st statement), then the metrization of the weak topology (2nd statement).
Separation. We have by Assumption 2 that 0∈E , hence by definition KALE(PP ||G)≥FP,G(0)=0. On the other hand, whenever P=G, it holds that:
FP,G(h)=− ∫ (exp(−h)+h−1)dP, ∀h∈H.
Moreover, by convexity of the exponential, we know that exp(−x)+x−1≥0 for all x∈R. Hence, FP,G(h)≤FP,G(0) = 0 for all h∈H. This directly implies that KALE(P|G) = 0. For the converse, we will use the same argument as in the proof of (Zhang et al., 2017, Theorem B.1). Assume that KALE(P|G)=0 and let h be inH. By Assumption 2, there existsCh>0 such that ch∈H and we have:
F(ch)≤KALE(P||G)=0. Now dividing by c and taking the limit to 0, it is easy to see that− ∫ hdP+ ∫ hdG≤0.Again, by Assump-
tion 2, we also know that−h∈H, hence, ∫ hdP− ∫ hdG≤0. This necessarily implies that ∫ hdP−∫
hdG=0 for allh∈H. By the density ofH in the set continuous functions on compact sets, we can conclude that the equality holds for any continuous and bounded function, which in turn implies that P=G.
Metrization of the weak topology. We first show that for any P and G with finite first moment, it holds that KALE(P|G)≤LW1(P,G), whereW1(P,G) is the Wasserstein-1 distance between P and G. For any h∈H the following holds:
F(h)=− ∫ hdP− ∫ exp(−h)dG+1
= ∫ h(x)dG(x)−h(x′)dP(x′)
− ∫
(exp(−h)+h−1)︸ ︷︷ ︸ ≥0 dG
≤ ∫ h(x)dG(x)−h(x′)dP(x′)≤LW1(P,G)
The first inequality results from the convexity of the exponential while the last one is a consequence of h being L-Lipschitz. This allows to conclude that KALE(P||G) ≤ LW1(P,G) after taking the supremum over all h∈H. Moreover, sinceW1 metrizes the weak convergence onP1 (Villani, 2009, Theorem 6.9), it holds that whenever a sequence Gn converges weakly towards P inP1 we also have W1(P,Gn)→0 and thus KALE(P||Gn)→0. The converse is a direct consequence of (Liu et al., 2017, Theorem 10) since by assumptionX is compact.
Well-defined learning. Assume that for any > 0 and any h and h′ in E there exists f in 2E such that ‖h+h′−f‖∞≤ then there exists a constantC such that:
KALE(P,Q)≤CKALE(P,G)
This means that the proposed learning procedure which first finds the optimal energy E? given a base G by maximum likelihood then minimizes KALE(P,G) ensures ends up minimizing the distance between the data end the generalized energy-based model Q.
KALE(P,Q)=sup h∈E LP,QG(h)
=−KALE(P,G)+sup h∈E LP,G(h+E?)
Let’s choose =KALE(P,G) and let h∈ 2E such that ‖h+E?−f‖∞≤ . We have by concavity of the function (α,β) 7→LP,G(α(h+E?−f)+βf) we have that:
LP,G(h+E?)≤2LP,G( 1
2 f)−LP,G(h+E?−f)
By assumption, we have that ‖h+E?−f‖∞≤ , thus |LP,G(h+E?−f)|≤2 . Moreover, we have thatLP,G( 12f)≤KALE(P,G) since 1 2f ∈E . This ensures that:
LP,G(h+E?)≤3KALE(P,G).
Finally, we have shown that:
KALE(P,Q)≤2KALE(P,G).
Hence, minimizing KALE(P,G) directly minimizes KALE(P,Q).

heading: C.2.2 SMOOTHNESS PROPERTIES OF KALE
text: We will now prove Theorem 5. We begin by stating the assumptions that will be used in this section:
(I) E is parametrized by a compact set of parameters Ψ. (II) Functions in E are jointly continuous w.r.t. (ψ,x) and are L-lipschitz and L-smooth w.r.t.
the input x:
‖Eψ(x)−Eψ(x′)‖≤Le‖x−x′‖, ‖∇xEψ(x)−∇xEψ(x′)‖≤Le‖x−x′‖.
(III) (θ,z) 7→Gθ(z) is jointly continuous in θ and z, with z 7→Gθ(z) uniformly Lipschitz w.r.t. z:
‖Gθ(z)−Gθ(z′)‖≤Lb‖z−z′‖, ∀z,z′∈Z,θ∈Θ.
There exists non-negative functions a and b defined from Z to R such that θ 7→Gθ(z) are a-Lipschitz and b-smooth in the following sense:
‖Gθ(z)−Gθ′(z)‖≤a(z)‖θ−θ′‖, ‖∇θGθ(z)−∇θGθ′(z)‖≤b(z)‖θ−θ′‖.
Moreover, a and b are integrable in the following sense:∫ a(z)2exp(2LeLb‖z‖)dη(z)<∞, ∫ exp(LeLb‖z‖)dη(z)<∞,
∫ b(z)exp(LeLb‖z‖)dη(z)<∞.
To simplify notation, we will denote byLθ(f) the expected Gθ log-likelihood under P. In other words, Lθ(E) :=LP,Gθ (E)=− ∫ EdP−log ∫ exp(−E)dGθ.
We also denote by pE,θ the density of the model w.r.t. Gθ,
pE,θ= exp(−E) ZGθ,E , ZGθ,E=
∫ exp(−E)dGθ.
We writeK(θ) :=KALE(P||Gθ) to emphasize the dependence on θ.
Proof of Theorem 5. To show that sub-gradient methods converge to local optima, we only need to show thatK is Lipschitz continuous and weakly convex. This directly implies convergence to local optima for sub-gradient methods, according to Davis and Drusvyatskiy (2018); Thekumparampil et al. (2019). Lipschitz continuity ensures thatK is differentiable for almost all θ∈Θ, and weak convexity simply means that there exits some positive constant C ≥ 0 such that θ 7→K(θ)+C‖θ‖2 is convex. We now proceed to show these two properties.
We will first prove that θ 7→K(θ) is weakly convex in θ. By Lemma 15, we know that for anyE∈E , the function θ 7→ Lθ(E) is M -smooth for the same positive constant M . This directly implies that it is also weakly convex and the following inequality holds:
Lθt(E)≤ tLθ(E)+(1−t)Lθ′(E)+ M
2 t(1−t)‖θ−θ′‖2.
Taking the supremum w.r.t. E, it follows that
K(θt)≤ tK(θ)+(1−t)K(θ′)+ M
2 t(1−t)‖θ−θ′‖2.
This means precisely thatK is weakly convex in θ. To prove that K is Lipschitz, we will also use Lemma 15, which states that Lθ(E) is Lipschitz in θ uniformly on E . Hence, the following holds:
Lθ(E)≤Lθ(E)+LC‖θ−θ′‖.
Again, taking the supremum overE, it follows directly that
K(θ)≤K(θ′)+LC‖θ−θ′‖.
We conclude that K is Lipschitz by exchanging the roles of θ and θ′ to get the other side of the inequality. Hence, by the Rademacher theorem,K is differentiable for almost all θ. We will now provide an expression for the gradient ofK. By Lemma 16 we know that ψ 7→Lθ(Eψ) is continuous and by Assumption (I) Ψ is compact. Therefore, the supremum supE∈ELθ(E) is achieved for some function E?θ . Moreover, we know by Lemma 15 that Lθ(E) is smooth uniformly on E , therefore the family (∂θLθ(E))E∈E is equi-differentiable. We are in position to apply Milgrom and Segal (2002)(Theorem 3) which ensures thatK(θ) admits left and right partial derivatives given by
∂+e K(θ)= lim t>0 t→0 ∂θLθ(E?θ+te)>e, ∂−e K(θ)= lim t<0 t→0 ∂θLθ(E?θ+te)>e, (28)
where e is a given direction in Rr. Moreover, the theorem also states that K(θ) is differentiable iff t 7→E?θ+te is continuous at t=0. Now, recalling thatK(θ) is actually differentiable for almost all θ, it must hold thatE?θ+te→t→0E?θ and ∂+e K(θ)=∂−e K(θ) for almost all θ. This implies that the two limits in (28) are actually equal to ∂θLθ(E?θ )>e. The gradient ofK, whenever defined, in therefore given by
∇θK(θ)=Z−1Gθ,E?θ
∫ ∇xE?θ (Gθ(z))∇θGθ(z)exp(−E?θ (Gθ(z)))η(z)dz.
Lemma 15. Under Assumptions (I) to (III), the functional Lθ(E) is Lipschitz and smooth in θ uniformly on E:
|Lθ(E)−Lθ′(E)|≤LC‖θ−θ′‖, ‖∂θLθ(E)−∂θLθ′(E))‖≤2CL(1+L)‖θ−θ′‖.
Proof. By Lemma 16, we have thatLθ(E) is differentiable, and that ∂θLθ(E) := ∫ (∇xE◦Gθ)∇θGθ(pE,θ◦Gθ)dη.
Lemma 16 ensures that ‖∂θLθ(E)‖ is bounded by some positive constantC that is independent from E and θ. This implies in particular thatLθ(E) is Lipschitz with a constantC. We will now show that it is also smooth. For this, we need to control the difference
D :=‖∂θLθ(E)−∂θLθ′(E)‖.
We have by triangular inequality: D≤ ∫ ‖∇xE◦Gθ−∇xE◦Gθ′‖‖∇θGθ‖(pE,θ◦Gθ)dη︸ ︷︷ ︸
I
+ ∫ ‖∇xE◦Gθ‖‖∇θGθ−∇θGθ′‖(pE,θ◦Gθ)dη︸ ︷︷ ︸
II
+ ∫ ‖∇xE◦Gθ∇θGθ‖|pE,θ◦Gθ−pE,θ′ ◦Gθ′ |dη︸ ︷︷ ︸
III
.
The first term can be upper-bounded usingLe-smoothness ofE and the fact thatGθ is Lipschitz in θ: I≤Le‖θ−θ′‖ ∫ |a|2(pE,θ◦Gθ)dη
≤LeC‖θ−θ′‖.
The last inequality was obtained by Lemma 17. Similarly, using that∇θGθ is Lipschitz, it follows by Lemma 17 that
II≤Le‖θ−θ′‖ ∫ |b|(pE,θ◦Gθ)dη
≤LeC‖θ−θ′‖.
Finally, for the last term III , we first consider a path θt= tθ+(1−t)θ′ for t∈ [0,1], and introduce the function s(t) :=pE,θt ◦Gθt . We will now control the difference pE,θ◦Gθ−pE,θ′ ◦Gθ′ , also equal to s(1)−s(0). Using the fact that st is absolutely continuous we have that s(1)−s(0)= ∫ 1 0 s′(t)dt. The derivative s′(t) is simply given by s′(t) = (θ−θ′)>(Mt−M̄t)s(t) whereMt= (∇xE◦Bθt)∇θGθt and M̄t= ∫ MtpE,θt ◦Gθtdη. Hence,
s(1)−s(0)=(θ−θ′)> ∫ 1
0
(Mt−M̄t)s(t)dt.
We also know thatMt is upper-bounded byLa(z),which implies III≤L2e‖θ−θ′‖ ∫ 1
0
Ç∫ |a(z)|2s(t)(z)dη(z)+ Å∫ a(z)s(t)(z)dη(z) ã2å ≤L2e(C+C2)‖θ−θ′‖,
where the last inequality is obtained using Lemma 17. This allows us to conclude thatLθ(E) is smooth for anyE∈E and θ∈Θ.
Lemma 16. Under Assumptions (II) and (III), it holds that ψ 7→ Lθ(Eψ) is continuous, and that θ 7→Lθ(Eψ) is differentiable in θ with gradient given by
∂θLθ(E) := ∫ (∇xE◦Gθ)∇θGθ(pE,θ◦Gθ)dη.
Moreover, the gradient is bounded uniformly in θ andE: ‖∇θLθ(E)‖≤Le Å∫ exp(−LeLb‖z‖)dη(z) ã−1∫ a(z)exp(LeLb‖z‖)dη(z).
Proof. To show that ψ 7→Lθ(Eψ) is continuous, we will use the dominated convergence theorem. We fixψ0 in the interior of Ψ and consider a compact neighborhoodW ofψ0. By assumption, we have that (ψ,x) 7→Eψ(x) and (ψ,z) 7→Eψ(Gθ(z)) are jointly continuous. Hence, |Eψ(0)| and |Eψ(Gθ(0))| are bounded onW by some constantC. Moreover, by Lipschitz continuity of x 7→Eψ , we have
|Eψ(x)|≤|Eψ(0)|+Le‖x‖≤C+Le‖x‖, exp(−E(Gθ(z)))≤exp(−E(Gθ(0)))exp(LeLb‖z‖)≤exp(C)exp(LeLb‖z‖).
Recalling that P admits a first order moment and that by Assumption (III), exp(LeLb‖z‖) is integrable w.r.t. η, it follows by the dominated convergence theorem and by composition of continuous functions that ψ 7→Lθ(Eψ) is continuous in ψ0. To show that θ 7→ Lθ(Eψ) is differentiable in θ, we will use the differentiation lemma in (Klenke, 2008, Theorem 6.28). We first fix θ0 in the interior of Θ, and consider a compact neighborhood V of θ0. Since θ 7→ |E(Gθ(0))| is continuous on the compact neighborhood V it admits a maximum valueC; hence we have using Assumptions (II) and (III) that
exp(−E(Gθ(z)))≤exp(−E(Gθ(0)))exp(LeLb‖z‖)≤exp(C)exp(LeLb‖z‖). Along with the integrability assumption in Assumption (III), this ensures that z 7→exp(−E(Gθ(z))) is integrable w.r.t η for all θ in V . We also have that exp(−E(Gθ(z))) is differentiable, with gradient given by
∇θexp(−E(Gθ(z)))=∇xE(Gθ(z))∇θGθ(z)exp(−E(Gθ(z))). Using thatE is Lipschitz in its inputs andGθ(z) is Lipschitz in θ, and combining with the previous inequality, it follows that
‖∇θexp(−E(Gθ(z)))‖≤exp(C)Lea(z)exp(LeLb‖z‖), where a(z) is the location dependent Lipschitz constant introduced in Assumption (III). The r.h.s. of the above inequality is integrable by Assumption (III) and is independent of θ on the neighborhood V . Thus (Klenke, 2008, Theorem 6.28) applies, and it follows that
∇θ ∫ exp(−E(Gθ0(z)))dη(z)= ∫ ∇xE(Gθ0(z))∇θGθ0(z)exp(−E(Gθ0(z)))dη(z).
We can now directly compute the gradient ofLθ(E), ∇θLθ(E)= Å∫ exp(−E(Gθ0))dη ã−1∫ ∇xE(Gθ0)∇θGθ0exp(−E(Gθ0))dη.
Since E and Gθ are Lipschitz in x and θ respectively, it follows that ‖∇xE(Gθ0(z))‖ ≤ Le and ‖∇θGθ0(z)‖≤a(z). Hence, we have
‖∇θLθ(E)‖≤Le ∫ a(z)(pE,θ◦Gθ(z))dη(z).
Finally, Lemma 17 allows us to conclude that ‖∇θLθ(E)‖ is bounded by a positive constant C independently from θ andE.
Lemma 17. Under Assumptions (II) and (III), there exists a constantC independent from θ andE such that ∫
ai(z)(pE,θ◦Gθ(z))dη(z)<C, (29)∫ b(z)(pE,θ◦Gθ(z))dη(z)<C,
for i∈1,2.
Proof. By Lipschitzness of E and Gθ, we have exp(−LeLb‖z‖)≤ exp(E(Gθ(0))−E(Gθ(z))≤ exp(LeLb‖z‖), thus introducing the factor exp(E(Bθ0(0)) in (29) we get∫ ai(z)(pE,θ◦Gθ(z))dη(e)≤Le Å∫ exp(−LeLb‖z‖)dη(z) ã−1∫
a(z)iexp(LeLb‖z‖)dη(z),∫ b(z)(pE,θ◦Gθ(z))dη(z)≤Le Å∫ exp(−LeLb‖z‖)dη(z) ã−1∫ b(z)exp(LeLb‖z‖)dη(z).
The r.h.s. of both inequalities is independent of θ andE, and finite by the integrability assumptions in Assumption (III).
D IMAGE GENERATION
Figures 3 and 4 show sample trajectories using Algorithm 3 with no friction γ=0 for the 4 datasets. It is clear that along the same MCMC chain, several image modes are explored. We also notice the transition from a mode to another happens almost at the same time for all chains and corresponds to the gray images. This is unlike Langevin or when the friction coefficient γ is large as in Figure 5. In that case each chain remains within the same mode.
Table 4 shows further comparisons with other methods on Cifar10 and ImageNet 32x32.

heading: E DENSITY ESTIMATION
text: Figure Figure 7 (left) shows the error in the estimation of the log-partition function using both methods (KALE-DV and KALE-F). KALE-DV estimates the negative log-likelihood on each batch of size 100 and therefore has much more variance than KALE-F which maintains the amortized estimator of the log-partition function.
Figure Figure 7 (right) shows the evolution of the negative log-likelihood (NLL) on both training and test sets per epochs for RedWine and Whitewine datasets. The error decreases steadily in the case of KALE-DV and KALE-F while the error gap between the training and test set remains controlled.
Larger gaps are observed for both direct maximum likelihood estimation and Contrastive divergence although the training NLL tends to decrease faster than for KALE.

heading: F ALGORITHMS
text: Estimating the variational parameter. Optimizing (9) exactly over A yields (8), with the optimal A equal to Ã= log( 1M ∑M m=1exp(−E(Ym))). However, to maintain an amortized estimator of the log-partition we propose to optimize (9) iteratively using second order updates:
Ak+1 =Ak−λ(exp(Ak−Ãk+1)−1), A0 =Ã0 (30)
where λ is a learning rate and Ãk+1 is the empirical log-partition function estimated from a batch of new samples. By leveraging updates from previous iterations, A can yield much more accurate estimates of the log-partition function as confirmed empirically in Figure 7 of Appendix E.
Tempered GEBM. It can be preferable to sample from a tempered version of the model by rescaling the energyE by an inverse temperature parameter β, thus effectively sampling from Q. High temperature regimes (β→0) recover the base model G while low temperature regimes (β→∞) essentially sample from minima of the energyE. As shown in Section 6, low temperatures tend to produce better sample quality for natural image generation tasks.
Training In Algorithm 1, we describe the general algorithm for training a GEBM which alternates between gradient steps on the energy and the generator. An additional regularization, denoted by I(ψ) is used to ensure conditions of Proposition 4 and Theorem 5 hold. I(ψ) can includeL2 regularization over the parameters ψ, a gradient penalty as in Gulrajani et al. (2017) or Spectral normalization Miyato et al. (2018). The energy can be trained either using the estimator in (8) (KALE-DV) or the one in (9) (KALE-F) depending on the variable C.
Sampling In Algorithm 3, we describe the MCMC sampler proposed in Sachs et al. (2017) which is a time discretization of (14).
Algorithm 2 Overdamped Langevin Algorithm 1: Input λ, γ, u,η,E,G 2: OuputXT 3: Z0∼η // Sample Initial latent from η. 4: for t=0,...,T do 5: Yt+1←∇zlogη(Zt)−∇zE◦B(Zt) // Evaluating∇zlog(ν(Zt+1)) using (4). 6: Wt+1∼N (0,I) // Sample standard Gaussian noise 7: Zt+1←Zt+λYt+1+ √ 2λWt+1
8: end for 9: XT←G(ZT )
Algorithm 3 Kinetic Langevin Algorithm 1: Input λ, γ, u,η,E,G 2: OuputXT 3: Z0∼η // Sample Initial latent from η. 4: for t=0,...,T do 5: Zt+1←Zt+ λ2Vt 6: Yt+1←∇zlogη(Zt+1)−∇zE◦B(Zt+1) // Evaluating∇zlog(ν(Zt+1)) using (4). 7: Vt+1←Vt+ uλ2 Yt+1. 8: Wt+1∼N (0,I) // Sample standard Gaussian noise 9: Ṽt+1←exp(−γλ)Vt+ 12 + √ u(1−exp(−2γλ))Wt+1 10: Vt+1← Ṽt+1+ uλ2 Yt+1 11: Zt+1←Zt+1+ λ2Vt+1 12: end for 13: XT←G(ZT )

heading: G EXPERIMENTAL DETAILS
text: In all experiments, we use regularization which is a combination of L2 norm and a variant of the gradient penalty Gulrajani et al. (2017). For the image generation tasks, we also employ spectral normalization Miyato et al. (2018). This is to ensure that the conditions in Proposition 4 and Theorem 5 hold. We pre-condition the gradient as proposed in Simsekli et al. (2020) to stabilize training, and to avoid taking large noisy gradient steps due to the exponential terms in (8) and (9). We also use the second-order updates in (30) for the variational constant cwhenever it is learned.
G.1 IMAGE GENERATION
z∈R100∼N (0,I) dense→Mg×Mg×512
4×4, stride=2 deconv. BN 256 ReLU 4×4, stride=2 deconv. BN 128 ReLU 4×4, stride=2 deconv. BN 64 ReLU
3×3, stride=1 conv. 3 Tanh
RGB image x∈RM×M×3 3×3, stride=1 conv 64 lReLU 4×4, stride=2 conv 64 lReLU 3×3, stride=1 conv 128 lReLU 4×4, stride=2 conv 128 lReLU 3×3, stride=1 conv 256 lReLU 4×4, stride=2 conv 256 lReLU 3×3, stride=1 conv 512 lReLU
dense→1.
Table 6: Energy/Discriminator of SNGAN ConvNet: M=32.
z∈R100∼N (0,I) dense, 4×4×256 ResBlock up 256 ResBlock up 256 ResBlock up 256 BN, ReLu, 3×3 conv, Tanh
Table 8: Base/Generator of SNGAN ResNet.
Network Architecture Table 5 and Table 6 show the network architectures used for the GEBM in the case of SNGAN ConvNet. Table 5 and Table 6 show the network architectures used for the GEBM in the case of SNGAN ResNet. The residual connections of each residual block consists of two convolutional layers proceeded by a BatchNormalization and ReLU activation: BN+ReLU+Conv+BN+ReLU+Conv as in (Miyato et al., 2018, Figure 8).
Training: We train both base and energy by alternating 5 gradient steps to learn the energy vs 1 gradient step to learn the base. For the first two gradient iterations and after every 500 gradient iterations on base, we train the energy for 100 gradient steps instead of 5. We then train the model up to 150000 gradient iterations on the base using a batch-size of 128 and Adam optimizer Kingma and Ba (2014) with initial learning rate of 10−4 and parameters (0.5,.999) for both energy and base.
Scheduler: We decrease the learning rate using a scheduler that monitors the FID score in a similar way as in Bińkowski et al. (2018); Arbel et al. (2018). More precisely, every 2000 gradient iterations on the base, we evaluate the FID score on the training set using 50000 generated samples from the base and check if the current score is larger than the score 20000 iterations before. The learning rate is decreased by a factor of 0.8 if the FID score fails to decrease for 3 consecutive times.
Sampling: For (DOT) Tanaka (2019), we use the following objective:
z 7→‖z−zy+ ‖+ 1
keff E◦G(z) (31)
where zy is sampled from a standard Gaussian, is a perturbation meant to stabilize sampling and keff is the estimated Lipschitz constant of E ◦B. Note that (31) uses a flipped sign for the E ◦B compared to Tanaka (2019). This is becauseE plays the role of−D whereD is the discriminator in Tanaka (2019). Introducing the minus sign in (31) leads to a degradation in performance. We perform 1000 gradient iterations with a step-size of 0.0001 which is also decreased by a factor of 10 every 200 iterations as done for the proposed method. As suggested by the authors of Tanaka (2019) we perform the following projection for the gradient before applying it:
g←g− (g >z) √ q z.
We set the perturbation to 0.001 and keff to 1 which was also shown in Tanaka (2019) to perform well. In fact, we found that estimating the Lipschitz constant by taking the maximum value of ‖∇E◦G(z)‖ over 1000 latent samples according to η lead to higher values for keff : ( Cifar10: 9.4, CelebA : 7.2, ImageNet: 4.9, Lsun: 3.8). However, those higher values did not perform as well as setting keff =1.
For (IHM) Turner et al. (2019) we simply run the MCMC chain for 1000 iterations.

heading: G.2 DENSITY ESTIMATION
text: Pre-processing We use code and pre-processing steps from Wenliang et al. (2019) which we describe here for completeness. For RedWine and WhiteWine, we added uniform noise with support equal to the median distances between two adjacent values. That is to avoid instabilities due to the quantization of the datasets. For Hepmass and MiniBoone, we removed ill-conditioned dimensions as also done in Papamakarios et al. (2017). We split all datasets, except HepMass into three splits. The test split consists of 10% of the total data. For the validation set, we use 10% of the remaining data with an upper limit of 1000 to reduce the cost of validation at each iteration. For HepMass, we used the sample splitting as done in Papamakarios et al. (2017). Finally, the data is whitened before fitting and the whitening matrix was computed on at most 10000 data points.
Regularization: We set the regularization parameter to 0.1 and use a combination ofL2 norm and a variant of the gradient penalty Gulrajani et al. (2017):
I(ψ)2 = 1
dψ ‖ψ‖2+E
î ‖∇xfψ(‹X)‖2ó
Network Architecture. For both base and energy, we used an NVP Dinh et al. (2016) with 5 NVP layers each consisting of a shifting and scaling layer with two hidden layers of 100 neurons. We do not use Batch-normalization.
Training: In all cases we use Adam optimizer with learning rate of 0.001 and momentum parameters (0.5,0.9). For both KALE-DV and KALE-F, we used a batch-size of 100 data samples vs 2000 generated samples from the base in order to reduce the variance of the estimation of the energy. We alternate 50 gradient steps on the energy vs 1 step on the base and further perform 50 additional steps on the energy for the first two gradient iterations and after every 500 gradient iterations on base. For Contrastive divergence, each training step is performed by first producing 100 samples from the model using 100 Langevin iterations with a step-size of 10−2 and starting from a batch of 100 data-samples. The resulting samples are then used to estimate the gradient of the of the loss.
For (CD), we used 100 Langevin iterations for each learning step to sample from the EBM. This translates into an improved performance at the expense of increased computational cost compared to the other methods. All methods are trained for 2000 epochs with batch-size of 100 (1000 on Hepmass and Miniboone datasets) and fixed learning rate 0.001, which was sufficient for convergence.
G.3 ILLUSTRATIVE EXAMPLE IN FIGURE 1
We consider parametric functionsG(1)θ andG (2) θ from R to R of the form:
G (1) θ (x)=sin(8πWx)/(1+4πBx), G (2) θ (x)=4πW ′x+b
with θ= (W,B,W ′,b). we also call θ?= (1,1,1,0). In addition, we consider a sigmoid like function h from [0,1] to [0,1] of the form:
z̃= tan(π(z− 1 2 ), h(z)= 1 2
Å z+
1
1+exp(−9z̃)
ã .
Data generation : To generate a data point X = (X1,X2), we consider the following simple generative model:
• Sample a uniform r.v. Z from [0,1]. • Apply the distortion function h to get a latent sample Y =h(Z).
• Generate pointX usingX1 =G(1)θ? (Y ) andX2 =G (2) θ? (Y ).
Hence, the data are supported on the 1-d line defined by the equationX2 =G (2) θ? (X1).
GAN For the generator we sample Z uniformly from [0, 1] then generate a sample (X1,X2)=(G (1) θ (Z),G (2) θ (Z)). The goal is to learn θ.
For the discriminator, we used an MLP with 6 layers and 10 hidden units.
GEBM For the base we use the same generator as in the GAN model. For the energy we use the same MLP as discriminator of the GAN model.
EBM To ensure tractability of the likelihood, we use the following model:
X2|X1∼N (G(2)θ (X1),σ0) X1∼MoG((µ1,σ1),(µ2,σ2))
MoG((µ1,σ1),(µ2,σ2)) refers to a Mixture of two gaussians with mean and variances µi and σi. We learn each of the parameters (θ,σ0,µ1,σ1,µ2,σ2) by maximizing the likelihood.
Both GAN and GEBM have the capacity to recover the the exact support by finding the optimal parameter θ?. For the EBM, when θ = θ?, the mean Gθ?(X1) of the conditional gaussian X2|X1 draws a line which matches the data support exactly, i.e.: X2 =G (2) θ? (X1).

heading: G.4 BASE/GENERATOR COMPLEXITY
text: To investigate the effect of model complexity of the performance gap between GANs and GEBMs, we performed additional experiments using the setting of Figure 1. Now we allow the generator/base network to better model the hidden transformation h that produces the first coordinate X1 given the latent noise. We choose G(1)θ to be either a one hidden layer network or an MLP with 3 hidden layers both with leaky ReLU activation, instead of a simple linear transform as previously done in Appendix G.3. The network has universal approximation capability that depends on the number of units. This provides a direct control over the complexity of the generator/base. We then varied the number of hidden units from 1 to 5∗104 units for the one hidden layer network and from 10 to 5∗103 units per layer for the MLP. Note that the MLP with 5∗103 units per layer stores a matrix of size 2.5∗107 and thus contains 2 orders of magnitudes more parameters than the widest shallow network with 5∗104 units. We then compared the performance of the GAN and GEBM using the Sinkhorn divergence Feydy et al. (2019) between each model and the data distribution. In all experiments, we used the same discriminator/energy network described in Appendix G.3. Results are provided in Figure 8.
Estimating the Sinkhorn divergence. The Sinkhorn is computed using 6000 samples from the data and the model, with squared euclidean distance as a ground cost and using a regularization =1e−3. We then repeat the procedure 5 times and average the result to get the final estimate of the Sinkhorn distance for a given run.
Training. Each run optimizes the parameters of the model using Adam optimizer (β1 = .5,β2 = .99), learning rate lr=1e−4 for the energy/discriminator and lr=1e−5 for the base/generator and weight decay of 1e−2 for the base/generator. Training is performed using KALE for 2000 epochs using a batch size of 5000 and 10 gradient iterations for the energy/discriminator per base/generator iteration. We use the gradient penalty for the energy/discriminator with a penalty parameter of 0.01. We then perform early stopping and retain the best performing model on a validation set.
Observations We make the following observations from Figure 8: the GAN generator indeed improves when we increase the number of hidden units. The performance of the GEBM remains stable as the number of hidden units increases. The performance of the GEBM is always better than the GAN, although we can see the GAN converging towards the GEBM. GEBM with a simpler base already outperforms the GAN with more powerful generators. The gap between the GEBM and the GAN reduces as the GAN becomes more expressive. Using a deeper network further reduces the gap compared to a shallow network.
These observations support the prior discussion that the energy witnesses a remaining difference between the generator and training samples, as long as it is not flat. This information allows the GEBM to perform better than a GAN that ignores it. The performance gap between the GEBM and the GAN reduces as the generator becomes more powerful and forces the energy to be more flat. This is consistent with the result in Proposition 3.

Review summaries:
This paper proposes a framework called GEBM that combine an implicit generator and an EBM to define a probabilistic model on low dimensional manifold. Specifically, the implicit generator defines the base distribution, and the EBM refines the base. Finally, this method is equivalent to define an EBM on the latent space of the implicit generator together with a mapping from the latent space to the data space. The authors propose to use the KALE to train the generator, and provide a theoretical guarantee about the validness of the KALE.