Bounds all around: training energy-based models with bidirectional bounds
Year: 2022

Abstract:
Energy-based models (EBMs) provide an elegant framework for density estimation, but they are notoriously difficult to train. Recent work has established links to generative adversarial networks, where the EBM is trained through a minimax game with a variational value function. We propose a bidirectional bound on the EBM log-likelihood, such that we maximize a lower bound and minimize an upper bound when solving the minimax game. We link one bound to a gradient penalty that stabilize training, thereby provide grounding for best engineering practice. To evaluate the bounds we develop a new and efficient estimator of the Jacobi-determinant of the EBM generator. We demonstrate that these developments stabilize training and yield high-quality density estimation and sample generation. 1 Energy-based models Energy-based models (EBMs) are probabilistic models that draw inspiration from physics and have a long history in machine learning (Hopfield, 1982; Hinton and Sejnowski, 1983; Smolensky, 1986). An EBM is specified in terms of an energy function Eθ : X → R that is parameterized by θ and defines a probability distribution over X from the Gibbs density: pθ(x) = exp(−Eθ(x)) Zθ , Zθ = ∫ exp(−Eθ(x))dx, (1) where Zθ is the normalization constant or partition function. In principle, any density can be described this way for a suitable choice ofEθ. EBMs are typically learned using maximum likelihood estimation (MLE), where we wish to find a value of θ that minimized the negative data log-likelihood: L(θ) := −Ex∼pdata(x) [log pθ(x)] = Ex∼pdata [Eθ(x)] + logZθ, (2) where pdata is the data generating distribution. The fundamental challenge with EBMs for non-standard energy functions is the lack of a closed-form expression for the normalization constant, Zθ, which hinders exact pointwise density evaluation, sampling and learning. Therefore, these tasks are traditionally performed using approximate methods such as Markov chain Monte Carlo (MCMC, Metropolis et al., 1953; Neal et al., 2011). Similar to the Boltzmann learning rule (Hinton and Sejnowski, 1983; Osogami, 2017), gradient-based learning for MLE involves evaluating the gradient ∇θL(θ) = Ex∼pdata [∇θEθ(x)]− Ex∼pθ [∇θEθ(x)], where the expectation in the last term necessitates MCMC approximation. Recently, a series of papers (Kim and Bengio, 2016; Kumar et al., 2019; Abbasnejad et al., 2020; Che et al., 2020) have established links between EBMs and Wasserstein GANs (WGANs, Arjovsky et al., 2017), such that EBMs can be approximately learned through a minimax game. Xie et al. (2017, 2018c) and Wu et al. (2018) also employed the adversarial game of EBM. This allows for significant improvements in EBM learning over MCMC based methods (Grathwohl et al., 2021). ∗equal contribution 35th Conference on Neural Information Processing Systems (NeurIPS 2021). In this paper, we remark that variational bounds on the value function of the minimax game can be problematic as this is both maximized and minimized. We propose to both upper and lower bound the negative log-likelihood, and alternate between their respective minimization and maximization to alleviate this concern. Evaluation of the bounds requires evaluating the entropy of the generator, and we provide a new efficient estimator. We link the upper bound to the use of gradient penalties which are known to stabilize training. Experimentally, we demonstrate that our approach matches or surpasses state-of-the-art on diverse tasks at negligible performance increase. 1.1 Background: Variational bounds and minimax games Following Grathwohl et al. (2021), we can lower bound the log normalization constant using a proposal (variational) distribution pg(x) and Jensen’s inequality: logZθ = logEx∼pg [ exp(−Eθ(x)) pg(x) ] ≥ Ex∼pg [ log exp(−Eθ(x)) pg(x) ] = −Ex∼pg [Eθ(x)] +H[pg] (3) where H[pg] = −Epg [log pg] is the (differential) entropy of the proposal distribution. This means that we also obtain a lower bound on the negative log likelihood function given by L(θ) := Ex∼pdata(x)[Eθ(x)]− Ex∼pg(x)[Eθ(x)] +H[pg] ≤ L(θ). (4) We note that this bound is tight when pg = pθ, since the bound on the log normalization constant (3) can be equivalently expressed as logZθ ≥ logZθ − KL (pg‖pθ). To tighten the bound, we seek a proposal distribution pg that maximises L(θ), while for MLE we want to find a θ that minimizes L(θ). Accordingly, MLE can be formulated as a minimax game: min Eθ max pg {L(θ)} = min Eθ max pg { Ex∼pdata(x)[Eθ(x)]− Ex∼pg(x)[Eθ(x)] +H[pg] } . (5) This gives a tractable and MCMC-free approach to EBM learning, assuming we can evaluate the entropy (Grathwohl et al., 2021). Our key issue is that both minimizing and maximizing a lower bound is potentially unstable. By minimizing a lower bound, we run the risk of finding an optimum, where the bound is loosest, rather than where it is informative about the optima of the true objective. In particular, a minima of the lower bound may be −∞, which is rather unhelpful. Good results have been reported from using a lower bound in the minimax game (Dai et al., 2017; Kumar et al., 2019; Grathwohl et al., 2021), but the conceptual issue remains. The WGAN loss function is, as noted by Arjovsky et al. (2017), very similar to the lower boundL(θ), LWGAN = Ex∼pdata [Eθ(x)]− Ex∼pg [Eθ(x)] , (6) with the missing entropy term being the only difference. Lessons from the WGAN can, thus, be expected to apply to the EBM setting as well. For example, the success of gradient clipping (Arjovsky et al., 2017) and gradient penalties (Gulrajani et al., 2017) in WGAN training, inspired Grathwohl et al. (2021) to heuristically use gradient penalties in variational inference to stabilize training. Our work will provide justification to such heuristics. 2 Approximate minimax games through bidirectional bounds Figure 1: The bidirectional bounds ‘sandwich’ the negative data log-likelihood. Instead of solving a minimax game with a lower bounded value function (5), we propose to bound L(θ) from both above and below: bL(θ)c ≤ L(θ) ≤ dL(θ)e. (7) With this, we can now follow an optimization strategy of alternating 1. minimize dL(θ)e with respect to Eθ. 2. maximize bL(θ)c with respect to pg . This avoids the potential pitfalls of minimizing a lower bound.

heading: 1 Energy-based models
text: Energy-based models (EBMs) are probabilistic models that draw inspiration from physics and have a long history in machine learning (Hopfield, 1982; Hinton and Sejnowski, 1983; Smolensky, 1986). An EBM is specified in terms of an energy function Eθ : X → R that is parameterized by θ and defines a probability distribution over X from the Gibbs density:
pθ(x) = exp(−Eθ(x))
Zθ , Zθ =
∫ exp(−Eθ(x))dx, (1)
where Zθ is the normalization constant or partition function. In principle, any density can be described this way for a suitable choice ofEθ. EBMs are typically learned using maximum likelihood estimation (MLE), where we wish to find a value of θ that minimized the negative data log-likelihood:
L(θ) := −Ex∼pdata(x) [log pθ(x)] = Ex∼pdata [Eθ(x)] + logZθ, (2) where pdata is the data generating distribution.
The fundamental challenge with EBMs for non-standard energy functions is the lack of a closed-form expression for the normalization constant, Zθ, which hinders exact pointwise density evaluation, sampling and learning. Therefore, these tasks are traditionally performed using approximate methods such as Markov chain Monte Carlo (MCMC, Metropolis et al., 1953; Neal et al., 2011). Similar to the Boltzmann learning rule (Hinton and Sejnowski, 1983; Osogami, 2017), gradient-based learning for MLE involves evaluating the gradient ∇θL(θ) = Ex∼pdata [∇θEθ(x)]− Ex∼pθ [∇θEθ(x)], where the expectation in the last term necessitates MCMC approximation. Recently, a series of papers (Kim and Bengio, 2016; Kumar et al., 2019; Abbasnejad et al., 2020; Che et al., 2020) have established links between EBMs and Wasserstein GANs (WGANs, Arjovsky et al., 2017), such that EBMs can be approximately learned through a minimax game. Xie et al. (2017, 2018c) and Wu et al. (2018) also employed the adversarial game of EBM. This allows for significant improvements in EBM learning over MCMC based methods (Grathwohl et al., 2021).
∗equal contribution
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In this paper, we remark that variational bounds on the value function of the minimax game can be problematic as this is both maximized and minimized. We propose to both upper and lower bound the negative log-likelihood, and alternate between their respective minimization and maximization to alleviate this concern. Evaluation of the bounds requires evaluating the entropy of the generator, and we provide a new efficient estimator. We link the upper bound to the use of gradient penalties which are known to stabilize training. Experimentally, we demonstrate that our approach matches or surpasses state-of-the-art on diverse tasks at negligible performance increase.

heading: 1.1 Background: Variational bounds and minimax games
text: Following Grathwohl et al. (2021), we can lower bound the log normalization constant using a proposal (variational) distribution pg(x) and Jensen’s inequality:
logZθ = logEx∼pg [ exp(−Eθ(x))
pg(x)
] ≥ Ex∼pg [ log
exp(−Eθ(x)) pg(x)
] = −Ex∼pg [Eθ(x)] +H[pg]
(3) where H[pg] = −Epg [log pg] is the (differential) entropy of the proposal distribution. This means that we also obtain a lower bound on the negative log likelihood function given by
L(θ) := Ex∼pdata(x)[Eθ(x)]− Ex∼pg(x)[Eθ(x)] +H[pg] ≤ L(θ). (4)
We note that this bound is tight when pg = pθ, since the bound on the log normalization constant (3) can be equivalently expressed as logZθ ≥ logZθ − KL (pg‖pθ). To tighten the bound, we seek a proposal distribution pg that maximises L(θ), while for MLE we want to find a θ that minimizes L(θ). Accordingly, MLE can be formulated as a minimax game:
min Eθ max pg {L(θ)} = min Eθ max pg
{ Ex∼pdata(x)[Eθ(x)]− Ex∼pg(x)[Eθ(x)] +H[pg] } . (5)
This gives a tractable and MCMC-free approach to EBM learning, assuming we can evaluate the entropy (Grathwohl et al., 2021).
Our key issue is that both minimizing and maximizing a lower bound is potentially unstable. By minimizing a lower bound, we run the risk of finding an optimum, where the bound is loosest, rather than where it is informative about the optima of the true objective. In particular, a minima of the lower bound may be −∞, which is rather unhelpful. Good results have been reported from using a lower bound in the minimax game (Dai et al., 2017; Kumar et al., 2019; Grathwohl et al., 2021), but the conceptual issue remains.
The WGAN loss function is, as noted by Arjovsky et al. (2017), very similar to the lower boundL(θ),
LWGAN = Ex∼pdata [Eθ(x)]− Ex∼pg [Eθ(x)] , (6)
with the missing entropy term being the only difference. Lessons from the WGAN can, thus, be expected to apply to the EBM setting as well. For example, the success of gradient clipping (Arjovsky et al., 2017) and gradient penalties (Gulrajani et al., 2017) in WGAN training, inspired Grathwohl et al. (2021) to heuristically use gradient penalties in variational inference to stabilize training. Our work will provide justification to such heuristics.

heading: 2 Approximate minimax games through bidirectional bounds
text: Instead of solving a minimax game with a lower bounded value function (5), we propose to bound L(θ) from both above and below:
bL(θ)c ≤ L(θ) ≤ dL(θ)e. (7)
With this, we can now follow an optimization strategy of alternating
1. minimize dL(θ)e with respect to Eθ. 2. maximize bL(θ)c with respect to pg .
This avoids the potential pitfalls of minimizing a lower bound.

heading: 2.1 A lower bound for maximizing pg
text: We already have a lower bound L(θ) on L(θ) in Equation (4). The first two terms of the bound can readily be evaluated by sampling from pdata and pg, respectively. The entropy of pg is less straightforward. Following standard practice in GANs, we define pg through a base variable z ∼ N (0, I), which is passed through a ‘generator’ network G : Rd → RD, i.e. x = G(z) ∼ pg . This constructs a density over a d-dimensional manifold2 in RD, which by standard change-of-variables is
log pg(G(z)) = log p0(z)− 1
2 log det(JᵀzJz), (8)
where Jz ∈ RD×d is the Jacobian ofG at z. This assume that the Jacobian exist, i.e. that the generator G has at least one derivative. The entropy of pg is then
H[pg] = −Ex∼pg [log pg(x)] = H[p0] + Ez∼p0 [ 1
2 log det(JᵀzJz)
] . (9)
The entropy of p0 trivially evaluates to H[p0] = d/2(1 + log(2π)). To avoid an expensive evaluation of the log-determinant-term, we note that it is easily bounded from below:
1 2 log det(JᵀzJz) = 1 2 d∑ i=1 log s2i ≥ d log s1, (10)
where sd ≥ . . . ≥ s1 are the singular values of the Jacobian Jz. Note that in general, this bound is not tight. With this, we have our final lower bound:
L(θ) ≥ bL(θ)c = Ex∼pdata(x)[Eθ(x)]− Ex∼pg(x)[Eθ(x)] +H[p0] + Ez∼p0 [d log s1(z)] . (11)

heading: 2.2 An upper bound for minimizing Eθ
text: There are many ways to arrive at an upper bound for the log-likelihood. Due to the strong ties between the lower bound (4) and the WGAN objective, we take inspiration from how WGANs are optimized. The original WGAN paper (Arjovsky et al., 2017) noted optimization instabilities unless gradient clipping was applied. Gulrajani et al. (2017) noted that this clipping could be avoided by adding a gradient penalty, such that the loss become
LWGAN-GP = LWGAN + λEx̂∼g(x) [ (‖∇x̂Eθ(x̂)‖2 − 1)2 ] , (12)
where x̂ is a uniform interpolation between real data and generated samples. This is generally acknowledged as a good way to train a WGAN (Gulrajani et al., 2017). We observe that this loss is similar to our lower bound, but the gradient penalty adds a positive term, such that one could speculate that we might have an upper bound. With this intuition in hand, we prove the following statement in the supplementary material.
Theorem 1 Suppose f : X → R is L-Lipschitz continuous, g(x) is a probability density function with finite support, then there exists constants M,m ≥ 0 and p ≥ 1 such that:
logEx∼g(x) [f(x)]− Ex∼g(x) [log f(x)] ≤M(Ex∼g(x) [|∇x log f(x)|p] +m) 1/p. (13)
Note that if Ex∼g(x) [|∇x log f(x)|p] +m ≥ 1, the bound in Equation (13) can be simplified by dropping the power 1/p,
log [ Ex∼g(x) [f(x)] ] − Ex∼g(x) [log f(x)] ≤MEx∼g(x) [|∇x log f(x)|p] +Mm. (14)
For our model, we set f(x) = exp(−Eθ(x))pg(x) and g(x) = pg(x). Empirically, we observe the bound simplification in Equation (14) holds through-out most of the training, and use the following upper bound:
dL(θ)e = L(θ) +MEx∼pg(x) [|∇xEθ(x) +∇x log pg(x)| p] +Mm. (15)
For a given choice of f(x), we can view m as a constant and the term Mm can be ignored during optimization. This upper bound can directly be interpreted as the lower bound (4) plus a gradient penalty, albeit one of a different form than the traditional WGAN penalty, which is derived purely from a regularization perspective. Our upper bound can, thus, be seen as a justification of the regularization from a maximum likelihood perspective.
2Note that Eq. 8 is only the density on the spanned manifold, and that the density is zero off the manifold.
Bound tightness When pθ = pg, we have that |∇xEθ(x) +∇x log pg(x)|p = 0. In Theorem 1, m is a constant related to the Lipschitz constant of log f(x) satisfying |∇x log f(x)|p ≤ |∇y log f(y)|p+m for all x,y (see proof for details). When pθ = pg we also have |∇x log f(x)|p = 0, such that m = 0. Our upper bound is then dL(θ)e = L(θ) = L(θ), and hence it is tight.

heading: 3 Numerical evaluation of the bounds
text: Evaluating the lower bound To evaluate the lower bound in Equation (11), we need the smallest singular value of the Jacobian J = ∂zG(z). Recall that this singular value satisfy s1 = ‖Jvmin‖2 = minv 6=0
‖Jv‖2 ‖v‖2 . We can then evaluate the singular value by finding vmin with an iterative optimization
algorithm, where we opt to use the celebrated single-vector LOBPCG algorithm (Knyazev, 1998). This method performs an iterative minimization of the generalized Rayleigh quotient,
ρ(v) := vᵀJᵀJv
vᵀv , (16)
which converges to vmin. The gradient of ρ(v) is proportional to r = JᵀJv − ρ(v)v. To avoid computing the Jacobian J, we use Jacobian-vector products, which can be efficiently evaluated using automatic differentiation. To compute JᵀJv, we use the following trick (in pytorch-notation):
JᵀJv = ((Jv)ᵀJ) ᵀ = ∇z ((Jv)ᵀ.detach() ·G(z))ᵀ . (17)
The optimal learning rate for this iterative scheme can be found by maximizing the Rayleigh quotient (16). Finally, we follow the suggestions of Knyazev (2001) to improve numerical stability and accelerate convergence, which we omit here for brevity.
Evaluating the upper bound There are two challenges when evaluating Equation (15). The first is to compute ∇x log pg(x), where we empirically found that existing methods (Shi et al., 2018; Li and Turner, 2018) were too inefficient for our needs. To evaluate the term Ex∼pg(x) [|∇xEθ(x) +∇x log pg(x)|p], we further loosen the bound
|∇G(z)Eθ(G(z)) +∇G(z) log pg(G(z))|p ≤ |∇G(z)Eθ(G(z))Jz +∇G(z) log pg(G(z))Jz|p
sp1
≤ |∇G(z)Eθ(G(z))Jz +∇z log pg(G(z))|p
sp1 .
(18) where s1 is the smallest singular value of Jz. Detailed derivations are in the supplementary material. If we choose p = 2, then we can use Hutchinson’s estimator (1989):
|∇xEθ(x)Jz +∇z log pg(G(z))|2 = Ev [ (∇xEθ(x)Jzv +∇z log pg(G(z))v)2 ] , (19)
where v ∼ N (0, Id). This is easily evaluated using automatic differentiation. The second challenge is to evaluate log pg(x) which needs the Jacobian of the generator G(z) as dictated by Equation (8). Here, we opt to use our entropy estimator as described above. We could alternatively use Hutchinson’s estimator as proposed by Kumar et al. (2020). Experimentally we do not observe much difference between these two estimators.

heading: 4 Related work
text: In machine learning, there has been a long-standing interest in EBMs dating back to Hopfield networks (Hopfield, 1982), Boltzmann machines (Hinton and Sejnowski, 1983; Ackley et al., 1985) and restricted Boltzmann machines (Smolensky, 1986; Hinton, 2002), see e.g. reviews in the works by LeCun et al. (2006) and Scellier (2020). Learning and evaluation of these models are difficult since the normalization constant cannot be efficiently evaluated. MLE-based learning, such as the Boltzmann learning rule, relies on expensive MCMC sampling to estimate the gradient, and more advanced MCMC methods are used to reliably estimate the normalization constant (see e.g. Salakhutdinov and Murray, 2008; Grosse et al., 2013; Burda et al., 2015; Frellsen et al., 2016). For images, MCMCbased learning has been used to learn non-deep EBMs of both textures (Zhu et al., 1998; Zhu and
Mumford, 1998) and natural images (Xie et al., 2015, 2016). Learning algorithms that avoid the costly MCMC approximation have been heavily investigated. For instance, Hinton (2002) proposed k-step Contrastive Divergence (CD-k) to approximate the negative phase log-likelihood gradient, and Hyvärinen (2005) proposed an alternative method to train non-normalized graphical models using score matching. Deep versions of EBMs have subsequently been proposed, such as deep belief networks (Hinton et al., 2006) and deep Boltzmann machines (Salakhutdinov and Hinton, 2009).
In recent years, there have been renewed interest in deep generative EBMs, particularly for image generation. Grathwohl et al. (2021) gave an excellent overview of the current developments and distinct drawbacks of three classes of learning methods, which we summarize here: (1) MLE methods with MCMC sampling are slow and unstable, and while (2) score matching based methods are comparable faster, they are also unstable and do not work with discontinuous nonlinearities. (3) Noise-contrastive estimation (Gutmann and Hyvärinen, 2010) do not have these drawback, but it does not scale well with the data dimensionality.
Our proposed method belongs to a fourth class of algorithms that sidestep the costly MCMC sampling by using a simultaneously learned generator or variational distribution. Kim and Bengio (2016) proposed using a generator function and an adversarial-like training strategy similar to ours. They update the generator using the same lower bound as us, but their entropy approximation is quite different. Furthermore, their method does not have a gradient penalty like our upper bound when optimizing the energy function. Consequently, the energy function needs to be explicitly designed to prevent it from growing to infinity, limiting its potential.
The method proposed by Zhai et al. (2016) plays a min-max game to jointly optimize the energy function and generator using the same lower bound as us. However, their method relies on a specific designed bounded multi-modal energy function, which limits its potential. Furthermore, their approximation of the generator entropy comes with no theoretical guarantees, and their regularisation of the energy function does not constitute an upper bound.
Dai et al. (2017) proposed an adversarial learning framework for jointly learning the energy function and a generator. They considered two approaches to maximising the generator entropy: One, which maximizes the entropy by minimizing the conditional entropy using a variational upper bound, and another, which makes isotropic Gaussian assumptions for the data, which is not suitable for highdimensional data. Kumar et al. (2019) and Abbasnejad et al. (2019) also consider adversarial learning but different approaches to estimating the entropy of the generator. Kumar et al. (2019) estimated the entropy through its connection to mutual information, but they need an additional network to measure the entropy term. Abbasnejad et al. (2019) maximized the entropy by approximating the generator function’s Jacobian log-determinant. However, their method is impractical in high dimensions as the Jacobian is computationally expensive. Han et al. (2019) use an adversarial learning strategy in their divergence triangle loss, but their training mechanism is radically different from ours, and they rely on an extra encoder for learning the generator.
Xie et al. (2018a,b, 2021a,b) proposed cooperative learning of the energy function and a generator. However, the cooperative learning approach relies on MCMC or Langevin dynamics to draw samples from the EBM, which is expensive and difficult to tune. Like our work, VERA by Grathwohl et al. (2021) avoids the use of MCMC. VERA plays a min-max game and uses variational inference to approximate the gradient of the entropy term, which is different from ours. Furthermore, VERA uses a gradient penalty as a regularizer for the energy function, which is a heuristic unlike our upper bound, and their method is memory-consuming, and the hyperparameters are difficult to adjust.

heading: 5 Experiments
text: To demonstrate the efficiency of our energy-based model with bidirectional bounds (EBM-BB) we compare against a range of methods that represent state-of-the-art. As representatives of the GAN literature, we consider deep convolutional GANs (DCGANs) by Radford et al. (2016), spectrally normalized GANS (SNGANs) by Miyato et al. (2018) and the WGAN with zero-centered gradient penalty (WGAN-0GP) by Thanh-Tung et al. (2019). As representatives of EBMs, we consider the maximum entropy generators (MEG) by Kumar et al. (2019) and the variational entropy regularized approximate maximum likelihood (VERA) estimator (Grathwohl et al., 2021). As representatives of CoopNets, we consider two similar methods (Xie et al., 2018a, 2021a) and EBM-VAE (Xie et al., 2021b). We also consider NCSN by Song and Ermon (2019) and DDPM by Ho et al. (2020)
as representatives of a score matching model and a diffusion model respectively. To investigate the influence of our derived upper bound, we introduce another baseline in which we replace our upper bound with a zero-centered gradient penalty. We denote the resulting model EBM-0GP. Our implementation is available at https://github.com/gengcong940126/EBM-BB.

heading: 5.1 Training details
text: On toy data and MNIST, all models are realized with multi-layer perceptrons (MLPs), while for natural images, we use the convolutional architecture from StudioGAN (Kang and Park, 2020). Details of the network architecture are given in the supplementary material. All experiments are conducted on a single 12GB NVIDIA Titan GPU using a pytorch (Paszke et al., 2017) implementation. To improve training, we use a positive margin for our energy function to balance the bounds. Specifically, we let
dL(θ)e = bL(θ)c+ [ MEx∼pg(x) [|∇xEθ(x) +∇x log pg(x)| p]− ζ ] + , (20)
where [·]+ = max(0, ·) is the usual hinge. In most experiments, ζ = 1 as this allows us to apply the simplified bound in Equation (15). But in general we recommend starting with ζ = 0 and trying successively larger values to stabilize training.
5.2 Toy data
Figure 2 shows estimated densities on the 25-Gaussians and swiss-roll datasets using both our methods, WGAN-0GP and two baselines. We observe that WGAN’s discriminator is not suitable as a density estimator. This is unsurprising as WGAN is not supposed to provide a density estimate. MEG and EBM-0GP can have inaccurate density information in some edge and peripheral areas. This may be due to insufficient or excessive maximization of entropy since the zero-centered gradient penalty is not a principled objective for maximum
likelihood. VERA and our EBM-BB can learn a sharp distribution, but our method is more stable in some inner regions of the 25-Gaussians and edge regions of the swiss-roll.

heading: 5.3 MNIST
text: Mode counting Mode collapse is a frequently occurring phenomenon in which the generator function maps all latent input to a small number of points in the observation space. In particular, GANs are plagued by this problem. Since the generator is trained to maximize the entropy of the generated distribution, several EBM-based models have been shown to capture all the modes of the data distribution. To empirically verify that our model also captures a variety of modes in the data distribution, we follow the test procedure of Kumar et al. (2019). We train our generative
model on the StackedMNIST dataset, a synthetic dataset created by stacking MNIST on different channels. The true total number of modes is 1,000, and they are counted using a pretrained MNIST classifier. The KL divergence is calculated empirically between the generated mode distribution and the data distribution. The results appear in Table 1. As expected, GAN-based methods suffer from mode collapse except for WGAN-0GP. All the EBM-based methods capture all or nearly all the modes of the data distribution. Our model EBM-0GP captures all modes and report the smallest KL measures. It worth mentioned that WGAN-0GP obtains comparable results, which was also observed by Grathwohl et al. (2021).
Entropy Estimation In Figure 3, we explore the quality of our entropy estimator measured on MNIST (LeCun, 1998). We use networks with fully connected layers as the Jacobian is then easily derived in closed-form, giving us a ground truth for the entropy. We compare our estimator (red) with the ground truth (black) and Hutchinson’s estimator (blue), as was proposed by Kumar et al. (2020). Finally, as we only run our iterative estimator for a few steps, we also compare with an estimator that uses the smallest singular value computed with high precision (green). We observe that Hutchinson’s estimator is reasonable close to the ground truth but provides an upper bound, making it inapplicable for our lower
bound. We further observe a noticeable gap between our estimator and the high precision singular value estimator. Finally, we see that all estimators follow roughly the same trend, which suggests that our estimator provides a suitable target for optimization.

heading: 5.4 Natural Images
text: Image generation The data studied thus far are simple and perhaps do not challenge our upper bound and entropy estimator. We, therefore, train our model on the standard benchmark 32 × 32 CIFAR-10 (Krizhevsky et al., 2009) dataset and the 64× 64 cropped ANIMEFACE3 dataset, which both represent a significant increase in complexity. Similar to recent work on GANs (Miyato et al., 2018), we report Inception Score (IS), Fréchet Inception Distance (FID) scores and two Fβ scores (Sajjadi et al., 2018). We compare with competitive GANs, EBM baselines, CoopNets, NCSN and DDPM (Table 3). All GANs and EBMs are reproduced using the same network architecture (see supplements for details) except DCGAN and IGEBM, where we respectively borrow the results from StudioGAN (Kang and Park, 2020) and the original implementation (Du and Mordatch, 2019). We use the DCGAN (Radford et al., 2016) network architecture for CIFAR-10 and a Resnet architecture (Kang and Park, 2020) for ANEMIFACE. Parameters are chosen as in the original papers. For VERA (Grathwohl et al., 2021), we were unable to reproduce the reported performance, so we choose hyper-parameters according to an extensive grid search. Like the original paper, we choose the entropy weight to be 0.0001. For CoopNets and NCSN, we report the results from the original papers. For DDPM, we used a public available implementation.4
From Table 3, we see that our EBM-BB model is the best performing EBM on CIFAR-10, though it is surpassed by NCSN and DDPM. This is not surprising as NCSN and DDPM focus on sample quality rather than optimising a data likelihood. Further note that sampling NCSN and DDPM is significantly more expensive than our method. On ANIMEFACE our method is highly competitive. Figures 4 and 5 show samples from different models. For ANIMEFACE, our method generates more diverse samples in terms of face parts than the baselines. All models predominantly generate female faces. This suggests that while our approach captures more density modes than the baselines, all approaches still misses several modes. We further draw attention to several corrupt samples generated by VERA, despite extensive parameter search. For CIFAR-10, we observe no immediate differences in the generative capabilities between models.
Capacity usage For the proposal distribution pg to be adaptive, the associated generator network G(z) should be able to use as much of its available capacity as possible. To compare different methods, we consider an implicit measure of capacity usage Cz that locally measures the intrinsic dimensionality of G(z). In particular, we use the anisotropy index (Wang and Ponce, 2021) that, for
3https://www.kaggle.com/splcher/animefacedataset 4https://github.com/rosinality/denoising-diffusion-pytorch
a given z, measure the standard deviation of the norm of the directional derivative of G(z) along the input dimensions, i.e.
Cz = std ({‖Jzei‖}i=1...d) , (21)
where std(·) computes the standard deviation of the input, and ei denotes the ith standard basis vector. A small Cz indicates that the different input dimensions contribute equally to the output of G, which imply good capacity usage. We measure the mean value of Cz for z ∼ N (0, Id) and report mean and standard deviation of the
result running for several times in Table 2. We observe that most models, with the exception of VERA, perform well on CIFAR-10, while on ANIMEFACE there is more diversity. In both cases, EBM-BB has the best capacity usage. The large difference between EBM-0GP and EBM-BB on
ANIMEFACE
WGAN-0GP 2.22 ± 0.030 9.76 ± 0.674 0.95 ± 0.005 0.98 ± 0.005 MEG 2.20 ± 0.020 9.31 ± 0.007 0.95 ± 0.005 0.98 ± 0.001
VERA 2.15 ± 0.001 41.00 ± 1.072 0.515 ± 0.078 0.78 ± 0.013 DDPM 2.18 8.81 0.94 0.98
EBM-0GP (ours) 2.26 ± 0.017 20.53 ± 0.524 0.889 ± 0.008 0.909 ± 0.019 EBM-BB (ours) 2.26 ± 0.005 12.75 ± 0.045 0.94 ± 0.001 0.96 ± 0.005
ANIMEFACE suggest that our proposed upper bound helps increase the entropy of the proposal distribution pg .
Out-of-distribution detection As EBMs are density estimators, they should, in principle, assign low likelihood to out-of-distribution (OOD) observations. OOD detection performance can then be seen as a proxy for the quality of the estimated density. To test how well our model fares in this regard, we follow previous work (Hendrycks and Gimpel, 2017; Hendrycks et al., 2019; Alemi et al., 2018; Choi et al., 2018; Ren et al., 2019; Havtorn et al., 2021) and report the threshold independent evaluation metrics of Area Under the Receiver Operator Characteristic (AUROC↑), Area Under the Precision Recall Curve (AUPRC↑) and False Positive Rate at 80% (FPR80↓), where the arrow indicates the direction of improvement of the metrics. The results are reported in Table 4. Each column of the table takes the form ‘In-distribution / Out-of-distribution’ in reference to the training and test set, respectively. We observe that on CIFAR-10 / SVHN, EBM-BB is the top performer. On CIFAR-10 / CIFAR-100 the overall performance degraded significantly, EBM-0GP fared noticeable better than other models. The overall degradation may be caused by the strong similarity between CIFAR-10 and CIFAR-100. For ANIMEFACE / Bedroom a similar situation occurs, even if these datasets are highly dissimilar, but we observe EBM-0GP is much better than other models.
No clear winner can be found from this study. We note that our method consistently performs well, but surprisingly, so does WGAN-0GP even if it is not a density estimator. As OOD detection is a task that comes with many subtle pitfalls (Havtorn et al., 2021), we suggest that these results should be taken with a grain of salt even if our model performs well.
WGAN-0GP 100000 13h MEG 100000 20h VERA 200000 44h
EBM-0GP 100000 38h EBM-BB 100000 40h
Running time Table 5 gives an overview of the time required to train the proposed models and the baselines. All models were trained on a single 12GB Titan GPU. We observe that our model should be trained with a lower learning rate, and therefore may increase the total number of epochs.

heading: 6 Discussion
text: The first observation in our work is that current methods for training energy-based models (EBMs) interchangeably minimize and maximize a lower bound. As this may be a potential source of training instability, we propose to bound the negative log-likelihood from above and below and switch between bounds when minimizing and maximizing. The lower bound, bL(θ)c, is similar to existing ones, but we provide a new algorithm for its realization. Unlike past work, this algorithm does not need additional networks and does not rely on hard-to-tune parameters; our only parameter controls the number of iterations for the Jacobian approximation and represents a trade-off between the accuracy of the bound and computational budget. The upper bound, dL(θ)e, is new to the literature but similar to the common regularization practice of introducing gradient penalties. To the best of our knowledge, this is the first time that gradient penalties have been derived from the perspective of bounding the log-likelihood. It is rewarding that current better engineering practice can be justified from a probabilistic perspective. Empirically, we find our proposed model generally performs as well or better than some of the current state-of-the-art on a variety of tasks. The evidence suggests that our bidirectional bounds allow the generator to increase its entropy and capacity usage. We see this both directly and through improved sample quality on diverse datasets.
Limitations The current drawbacks of the method are mainly three concerns. First, we find that a smaller-than-usual learning rate helps our model. While we expect this to be a matter of implementation rather than a more profound concern, the current implication is that training is approximately twice that of a EBM baseline. Second, our bounds rely on computing the smallest singular value of the Jacobian of the generator. We provide an efficient implementation of a method for this task, but it needs to run for several iterations to guarantee convergence. In practice, we stop after a fixed, low number of iterations, which we found to work well, but technically this violates the bound. This seems to be an unavoidable aspect of our approach, but we find that the benefits significantly compensate for this issue. The availability of a fast (approximate) entropy estimator that does not require additional networks and parameters is highly valuable when training EBMs. Finally, for the upper bound, we need an estimate of the volume of support M of the proposal distribution pg . We do not have a viable method for this, and, in practice, we treat M as a hyper-parameter. We have not found it difficult to tune this parameter, but it nonetheless constitutes a limitation.
Negative societal impact All high-capacity generative models carry the risk of being used for misinformation, and our model is no exception. The value of EBMs over GANs is that they come with a likelihood function, which is more valuable in data analysis, than e.g. for creating deepfakes.
Acknowledgments and Disclosure of Funding
This work was funded in part by the Novo Nordisk Foundation through the Center for Basic Machine Learning Research in Life Science (NNF20OC0062606). It also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research, innovation programme (757360), National science foundation of China under grant 61771305 and Shanghai Municipal Science and Technology Major Project (2021SHZDZX0102). JF was supported in part by the Novo Nordisk Foundation (NNF20OC0065611) and the Independent Research Fund Denmark (9131- 00082B). SH was supported in part by a research grant (15334) from VILLUM FONDEN. The authors are grateful to anonymous reviewers and the handling area chair for valuable discussions and feedback on an early version of this manuscript.

Review summaries:
This paper proposes a bidirectional bound on the EBM log-likelihood, such that it can maximize a lower bound and minimize an upper bound when solving the minimax game of EBM. The paper demonstrates the effectiveness of the proposed method on density estimation and sample generation. The major contribution of the current paper is to use the bidirectional bounds to train the EBMs.

This work proposes a new method for density estimation using Energy-based model. Classical variational bounds replace the negative log-likelihood by a computable lower bound, and phrase the maximum likelihood problem as a minimax game on this lower bound. A possible algorithm to solve the problem is then to alternate between minimization and maximization of this lower bound (with respect to different variables). Although this approach has shown good performance, the authors claim that the concept of minimizing such a lower bound can lead to instabilities. In particular, they claim that the WGAN objective tends to be unstable to optimize, and that while clipping the gradient in this case solves the instability issue, the origin of the instability lies in minimizing an upper bound. To circumvent this issue, the current work constructs both an upper and a lower bound to the negative log-likelihood, and alternate between minimizing the upper bound and maximizing the lower bound. The upper bound is inspired by what is done in WGAN, i.e., by adding a gradient penalty. The algorithm is tested on several benchmarks, and is shown to be rather competitive compared to existing algorithms.

This paper proposes a method to train EBMs with lower and upper variational bounds. Previous approaches used only one side of the bound. This causes significant training instabilities as the training objective involves alternately maximizing and minimizing the same quantity.  The lower bound is estimated by lower-bounding the log determinant of the jacobian multiplied with itself with the smallest singular value of the jacobian.   The upper bound is inspired by tricks used to make WGAN and related EBM training approaches stable by regularizing the gradient.   The minimum singular value, needed for computing the lower bound, is computed using an iterative algorithm run for a small number of steps. Automatic differentiation is used cleverly to ensure full Jacobians don’t need to be instantiated.  The upper bound is loosened, and Hutchinson’s trace estimator is employed to efficiently compute.  

The authors propose to sandwich the EBM log-likelihood with newly introduced upper and lower bounds, with which a minimax objective is formed to train EBMs.  Specifically, based on the variational lower bound from Grathwohl et al. (2021), the authors approximate the entropy term therein to form the presented lower bound, and then adding a gradient-penalty term mimicking WGAN to form the new upper bound.  The lower and upper bounds are optimized wrt different groups of parameters, as, intuitively, both minimizing and maximizing a lower bound wrt the same group of parameters is potentially unstable (this is the main motivation for this submission).