marginparsep has been altered.
topmargin has been altered.
marginparwidth has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

Clarifying MCMC-based training of modern EBMs :
Contrastive Divergence versus Maximum Likelihood

Anonymous Authors1

††footnotetext:  1Anonymous Institution, Anonymous City, Anonymous Region,
Anonymous Country. Correspondence to: Anonymous Author
<EMAIL_ADDRESS>
Preliminary work. Under review by the International Conference on Machine
Learning (ICML). Do not distribute.

###### Abstract

The Energy-Based Model (EBM) framework is a very general approach to
generative modeling that tries to learn and exploit probability distributions
only defined though unnormalized scores. It has risen in popularity recently
thanks to the impressive results obtained in image generation by
parameterizing the distribution with Convolutional Neural Networks (CNN).
However, the motivation and theoretical foundations behind modern EBMs are
often absent from recent papers and this sometimes results in some confusion.
In particular, the theoretical justifications behind the popular MCMC-based
learning algorithm Contrastive Divergence (CD) are often glossed over and we
find that this leads to theoretical errors in recent influential papers Du &
Mordatch (2019); Du et al. (2020). After offering a first-principles
introduction of MCMC-based training, we argue that the learning algorithm they
use can in fact not be described as CD and reinterpret theirs methods in light
of a new interpretation. Finally, we discuss the implications of our new
interpretation and provide some illustrative experiments.

##  1 Introduction

Learning statistical models of observational data has always been one of the
primary goals in machine learning and more generally in statistics. With such
a model in hand, one can understand how the data is distributed and—of
particular importance for this work—generate new samples as if they were
coming from the data distribution.

To learn such a model, the first step is always to hypothesize a parametric
family of distributions $p_{\theta}(x)$ which we think encompasses the real
distribution—or a good approximation of it. Then we can use various
statistical tools to infer which $\theta\in\Theta$ fits the data the best.

Any choice of parametric family can be written as

$p_{\theta}(x)=\frac{1}{Z_{\theta}}\tilde{p}_{\theta}(x)$

where $\tilde{p}_{\theta}\geq 0$ is the unnormalized probability and
$Z_{\theta}=\int\tilde{p}_{\theta}(x)dx$ is the normalization constant.
Energy-based models (EBMs) consist in paremeterizing the family by defining
the unnormalized probability in the following way

$p_{\theta}(x)=\frac{1}{Z_{\theta}}e^{-E_{\theta}(x)}$

where $E_{\theta}(x)\geq 0$ is called the energy function and associates a
scalar measure of incompatibility to every point. The particular form of this
distribution comes from statistical physics and is explained in section 2. The
main advantage of EBMs is to give the most flexibility possible to how we
parameterize the parametric family : $E_{\theta}(x)$ can be any function and
does not need to be normalized. In particular, recent deep learning work have
exploited that to model image distribution by parametrizing $E_{\theta}(x)$
with a CNN : the ultimate goal being to generate new images coming from the
learned distribution.

However, we notice in these papers a growing confusion related to the
theoretical underpinnings of the MCMC-based training procedure. In particular,
we argue that two recent papers Du & Mordatch (2019); Du et al. (2020) wrongly
identify the learning algorithm they use as Contrastive Divergence, learning
to confusion and unforeseen consequences. More precisely, in this paper :

  * •

We offer a much needed first-principles introduction to EBM sampling and
training with MCMC-based objectives : Negative Log Likelihood (NLL) and
Contrastive Divergence (CD).

  * •

We perform a detailed critique of two recent paper Du & Mordatch (2019); Du et
al. (2020) and argue that objective they are using is NLL while they are
claiming it is CD. Moreover, we identify a mathematical error in the
derivation of an auxiliary loss in (Du et al., 2020) and show how, with our
new interpretation of their learning algorithm, it ends up being a good thing.

  * •

We perform some simple experiments on MNIST to illustrate our main arguments.

##  2 Physics of EBMs

###  2.1 The Boltzmann-Gibbs distribution

Suppose that we have a system $S$ which is described by a state $x\in X$ and
that this system is rapidly fluctuating due to random interactions. Next
suppose that every state is associated with an energy function $E(x)$ and that
we know that in expectation, when we observe the system, it is going to have
energy $\langle E\rangle<\infty$. In the language of physics, the second
assumption is that the system is at thermal equilibrium with a thermal bath at
a fixed temperature (i.e. $S$ interacts with another system that maintains it
a fixed average energy). Then, the Boltzmann-Gibbs distribution gives us the
probability that we observe the system $S$ in state $x$ :

$p(x)=\frac{1}{Z}e^{-E(x)/T}$

Where $T$ is the temperature and is related to $\langle E\rangle$. This tells
us that the likelihood of a state decreases exponentially when its energy goes
up : low energy states are preferred. It turns out that $p(x)$ can be derived
as the MaxENT distribution with constraint $\mathbb{E}_{x\sim p}[E(x)]=\langle
E\rangle$ Jaynes (1957). We can therefore understand $p(x)$ as the natural
tendency of the system given that no states are a priori preferred (i.e. under
the second law of thermodynamics). See Huembeli et al., 2021 for a thorough
discussion of physics and EBMs.

Importantly, notice that what we have said here is equivalent to saying that
the stationary distribution of states visited by a certain trajectory of the
system (in contact with a thermal bath at a fixed temperature) is the
Boltzmann-Gibbs distribution : if you start the system in state $x_{0}$ and
let it run to infinity, it will visit visit each state $x$ in proportion to
$p(x)$.

###  2.2 EBMs as a physical simulation

Following on this last point, consider the effect of doing MCMC sampling on a
EBM with energy function $E_{\theta}(x)$. By definition, if the MCMC algorithm
is valid, it will have as stationary distribution
$p_{\theta}(x)=\frac{1}{Z_{\theta}}e^{-E_{\theta}(x)}$. But this is exactly
the stationary distribution that a simulation an hypothetical system with
energy $E_{\theta}(x)$ would have!

This means that we can think of running a MCMC chain on an EBM (i.e. sampling)
as just simulating the hypothetical system defined by $E_{\theta}(x)$ in
contact with a heat bath at a fixed temperature (in fact we suppose that $T=1$
for the rest of the paper).

Actually, the MCMC technique that we are going to use (section 2.3) is
explicitly the simulation of a physical system.

###  2.3 Langevin sampling

Because our ultimate objective is to parametrize the energy function with a
neural network, we want a MCMC sampling algorithm that is general (unlike
Gibbs sampling which relies on knowing conditional densities of input
dimensions) and efficient (unlike random-walk Metropolis-Hastings which has
very low acceptance probability in high dimension). The generally adopted
strategy, Langevin sampling, is the discretisation of a differential equation
which describes the motion of a certain physical system.

####  2.3.1 Pollen in a glass of water

Suppose you drop a grain of pollen in a glass of still water. First, because
the particle wants to reduce its gravitational potential energy $U(x)$, it is
going to move toward the bottom of the cup. However, because it is subjected
to random molecular bombardment, a random component will be added to its
trajectory. Because of the discussion in [section 2.1], it should make sense
that the distribution of the particle’s position follows the Boltzmann-Gibbs
distribution $p(x)=\frac{1}{Z}e^{-U(x)}$. Equivalently, the stationary
distribution of the differential equation describing the motion of the pollen
grain (the Langevin equation) is the Boltzmann-Gibbs distribution.

Importantly, we can extrapolate this simple example to any potential we want :
in our case $E_{\theta}(x)$.

####  2.3.2 Discretization

Glossing over the details, the Langevin sampling algorithm is the
discretization of a certain version of the Langevin differential equation and
can be writen as :

$\displaystyle
x^{t}=x^{t-1}-\frac{\lambda}{2}\underbrace{\nabla_{x}E_{\theta}(x^{t-1})}_{\text{Move
down potential}}+\underbrace{\omega^{k}}_{\text{Random bombardment}}$ (1)

where $\omega^{k}\sim N(0,\lambda)$. It can be shown that this converges to
$p_{\theta}$ (i.e. is a valid MCMC algorithm) in the limit of
$t\rightarrow\infty$ and $\lambda\rightarrow 0$ Welling & Teh (2011). For a
given $\lambda$, Metropolis-Hastings corrections are needed to account for
discretization error : the resulting algorithm is called Metropolis-Adjusted
Langevin Algorithm (MALA) Xifara et al. (2014). In practice, these corrections
are ignored and Langevin sampling is conceptually treated like a noisy
gradient descent on the energy landscape, i.e. it is the approximation of a
valid MCMC algorithm. See Figure 1 for an example of energy landscape.

###  2.4 Generating new data with EBMs

Note that once we have an EBM that correctly models the data distribution, we
can generate new ”imagined” data points by starting a sampling chain at a
random point (i.e. random noise for images) and let it drop to a low-energy
point (i.e. a real-looking image).

##  3 Learning EBMs with maximum likelihood

Let us see what happens when we try to learn the parameters of an EBM with
maximum likelihood estimation, i.e. minimizing the expected negative log-
likelihood of the data

$\displaystyle\mathcal{L}_{\text{NLL}}$
$\displaystyle=-\mathbb{E}_{p_{\text{data}}}[\log p_{\theta}(x)]$
$\displaystyle=-\mathbb{E}_{p_{\text{data}}}[\log\frac{1}{Z_{\theta}}e^{-E_{\theta}(x)}]$
$\displaystyle=\mathbb{E}_{p_{\text{data}}}[E_{\theta}(x)]+\log Z_{\theta}$

Computing the gradient gives

$\nabla_{\theta}\mathcal{L}_{\text{NLL}}=\mathbb{E}_{p_{\text{data}}}[\nabla_{\theta}E_{\theta}(x)]+\nabla_{\theta}\log
Z_{\theta}$

where the second term can be rewritten as

$\displaystyle\nabla_{\theta}\log Z_{\theta}$
$\displaystyle=\frac{\nabla_{\theta}Z_{\theta}}{Z_{\theta}}$
$\displaystyle=\frac{\nabla_{\theta}\int e^{-E_{\theta}(x)}dx}{Z_{\theta}}$
$\displaystyle=\frac{\int\nabla_{\theta}e^{-E_{\theta}(x)}dx}{Z_{\theta}}$
$\displaystyle=\frac{\int
e^{-E_{\theta}(x)}\nabla_{\theta}(-E_{\theta}(x))dx}{Z_{\theta}}$
$\displaystyle=-\int\frac{e^{-E_{\theta}(x)}}{Z_{\theta}}\nabla_{\theta}E_{\theta}(x)dx$
$\displaystyle=-\mathbb{E}_{p_{\theta}}[\nabla_{\theta}E_{\theta}(x)]$

Figure 1: Illustration of a complex energy landscape taken from the internet
to help with intuition. Sampling (noisy gradient descent) will spend most of
its time in low energy : the stationary distribution is
$p(x)=\frac{1}{Z}e^{-E(x)}$

Therefore, we can write the learning gradient as

$\displaystyle\nabla_{\theta}\mathcal{L}_{\text{NLL}}$
$\displaystyle=\underbrace{\mathbb{E}_{p_{\text{data}}}[\nabla_{\theta}E_{\theta}(x)]}_{\text{Positive
phase}}\underbrace{-\mathbb{E}_{p_{\theta}}[\nabla_{\theta}E_{\theta}(x)]}_{\text{Negative
phase}}$ (2)
$\displaystyle=\nabla_{\theta}\Big{(}\underbrace{\mathbb{E}_{p_{\text{data}}}[E_{\theta}(x)]}_{\mathcal{L}_{+}}\underbrace{-\mathbb{E}_{p_{\theta}}[E_{\theta}(x)]}_{\mathcal{L}_{-}}\Big{)}$
(3)

We can separate the loss in two phases of learning called positive and
negative. The positive phase consists in minimizing the energy of points drawn
from $p_{\text{data}}$ and is straightforward to compute assuming that
$E_{\theta}(x)$ is differentiable (like a CNN). This has the effect of making
the data points more likely (i.e. push down on their energy). However, this
alone would lead to a flat energy landscape (i.e. 0 everywhere). To counteract
this, the negative phase pushes up on the energy of points drawn from the
model distribution $p_{\theta}$ : the energy of the points that are unlikely
in $p_{\text{data}}$ will only get pushed up while those that are likely are
going to get balanced out by the positive phase and will stay relatively low.
It is the balancing of these two forces that shapes the energy landscape.

The naïve way to go about computing the negative phase is to collect samples
from $p_{\theta}$ using MC chains and then make a Monte Carlo estimate of the
expectation. The problem is that in practice the MCMC chains mix very poorly
on complex data. For example, in the space of pixels, plausible images under
the data distribution are separated by big distances of highly unlikely
images.

In other words, low energy regions in the energy landscape are often separated
by immense mountains of high energy which Langevin chains have no chance of
traversing in finite time. Hence the samples produced by a chain are not
representative of $p_{\theta}$ and the expectation is strongly biased. The
next section proposes an alternative method for training EBMs that sidesteps
this problem.

##  4 Contrastive Divergence training

First, consider the NLL objective written in terms of the KL-divergence (the
constant term corresponding to the entropy of the data distribution is omitted
because it does not change the objective):

$\mathcal{L}_{\text{NLL}}=D_{\text{KL}}[p_{\text{data}}(x)\|p_{\theta}(x)]$

The Contrastive Divergence (CD) objective, introduced by Hinton, 2002, can be
expressed as follows

$\displaystyle\mathcal{L}_{\text{CD}}$
$\displaystyle=\underbrace{D_{\text{KL}}[p_{\text{data}}(x)\|p_{\theta}(x)]}_{\mathcal{L}_{\text{NLL}}}-D_{\text{KL}}[q_{\theta}^{t}(x)\|p_{\theta}(x)]$

where $q_{\theta}^{t}$ is the distribution resulting from running $t$ steps of
the MC chain starting from the data. Notice that the added term is the
negative of another NLL loss : it tries to minimize
$\mathbb{E}_{q_{\theta}^{t}}[\log p_{\theta}]$, i.e. push $q_{\theta}^{t}$
away from $p_{\theta}$.

###  4.1 Equivalence to NLL

Because $q^{0}_{\theta}(x)=p_{\text{data}}(x)$ and
$q^{\infty}_{\theta}(x)=p_{\theta}(x)$, the CD loss is often expressed like

$\mathcal{L}_{\text{CD}}=D_{\text{KL}}[q_{\theta}^{0}(x)\|q_{\theta}^{\infty}(x)]-D_{\text{KL}}[q_{\theta}^{t}(x)\|q_{\theta}^{\infty}(x)]$

Assuming the MCMC sampling algorithm is valid (the stationary distribution of
the chains is $p_{\text{data}}$), $q_{\theta}^{t}$ is $t$ step closer to the
equilibrium than $q_{\theta}^{0}=p_{\text{data}}$ and we are guaranteed that
$\mathcal{L}_{\text{CD}}$ non-negative. Also the two terms are equal if and
only if $q_{\theta}^{t}=q_{\theta}^{0}=p_{\text{data}}$ which in turn imply
that $q_{\theta}^{\infty}=p_{\text{data}}$, i.e. $\mathcal{L}_{\text{NLL}}=0$.
Therefore,

$\mathcal{L}_{\text{CD}}=0\iff\mathcal{L}_{\text{NLL}}=0$

so the two objectives have the same solution, but they get there in different
ways. The next sections elaborate on this fact. First we motivate this new
loss mathematically and then discuss how to think about it versus the NLL
loss.

###  4.2 Derivation of the gradient

We can write the gradient of the added term like :

$\displaystyle\nabla_{\theta}D_{\text{KL}}[q_{\theta}^{t}\|p_{\theta}(x)]=\frac{\partial
D_{\text{KL}}[q_{\theta}^{t}\|p_{\theta}]}{\partial
q_{\theta}^{t}}\frac{\partial q_{\theta}^{t}}{\partial\theta}+\frac{\partial
D_{\text{KL}}[q_{\theta}^{t}\|p_{\theta}]}{\partial p_{\theta}}\frac{\partial
p_{\theta}}{\theta}$
$\displaystyle=\frac{\partial}{\partial\theta}D_{\text{KL}}[q_{\theta}^{t}\|p_{\Omega(\theta)}]+\frac{\partial}{\partial\theta}D_{\text{KL}}[q_{\Omega(\theta)}^{t}\|p_{\theta}]$
$\displaystyle=\frac{\partial}{\partial\theta}D_{\text{KL}}[q_{\theta}^{t}\|p_{\Omega(\theta)}]+\Big{(}\underbrace{\mathbb{E}_{q_{\Omega(\theta)}^{t}}[\nabla_{\theta}E_{\theta}(x)]}_{\text{Divergence
phase}}\underbrace{-\mathbb{E}_{p_{\theta}}[\nabla_{\theta}E_{\theta}(x)]}_{\text{Negative
phase}}\Big{)}$

where $\Omega(\theta)$ means we consider $\theta$ as constant. The last line
is obtained essentially the same way as (2), replacing $p_{\text{data}}$ with
$q_{\Omega(\theta)}^{t}$. This means that the expectation over $p_{\theta}$
(the problematic negative phase) cancel in the full loss :

$\displaystyle\nabla_{\theta}\mathcal{L}_{\text{CD}}=\nabla_{\theta}\mathcal{L}_{\text{NLL}}-\nabla_{\theta}D_{\text{KL}}[q_{\theta}^{t}\|p_{\theta}(x)]$
$\displaystyle=\underbrace{\mathbb{E}_{p_{\text{data}}}[\nabla_{\theta}E_{\theta}(x)]}_{\text{Positive
phase}}\underbrace{-\mathbb{E}_{q_{\Omega(\theta)}^{t}}[\nabla_{\theta}E_{\theta}(x)]}_{\text{Divergence
phase}}-\underbrace{\nabla_{\theta}D_{\text{KL}}[q_{\theta}^{t}\|p_{\Omega(\theta)}]}_{\text{KL
term}}$
$\displaystyle=\nabla_{\theta}\Big{(}\underbrace{\mathbb{E}_{p_{\text{data}}}[E_{\theta}(x)]}_{\mathcal{L}_{+}}\underbrace{-\mathbb{E}_{q_{\Omega(\theta)}^{t}}[E_{\theta}(x)]}_{\mathcal{L}_{d}}\underbrace{-D_{\text{KL}}[q_{\theta}^{t}\|p_{\Omega(\theta)}]}_{\mathcal{L}_{\text{KL}}}\Big{)}$

Observe that the new gradient doesn’t require us to get any samples from the
stationary distribution : what we couldn’t do anyway. So one way to view CD is
like a neat cancellation trick.

###  4.3 Divergence phase vs Negative phase

Notice that the divergence phase is very similar to—and in some sense an
approximation of—the negative phase of the NLL. However it has quite a
different interpretation, mainly because in the divergence phase the chains
always start from the data distribution. First of all, notice that the
negative phase comes from the minimization of
$\mathcal{L}_{\text{NLL}}=D_{\text{KL}}[p_{\text{data}}(x)\|p_{\theta}(x)]$
while the divergence phase comes from the maximization of
$D_{\text{KL}}[q_{\Omega(\theta)}^{t}\|p_{\theta}]$ in
$\mathcal{L}_{\text{CD}}$.

Instead of pushing up on the model distribution, the divergence phase attempts
to stop MC chains from diverging (hence the name Contrastive Divergence) from
the data distribution after running for $t$ steps (by setting up big mountains
on the borders of the data manifold). This is in some sense a different way of
doing the same thing : balancing out the positive phase (remember discussion
in section 3.

This alternative way of approaching the negative phase is the main way in
which CD differs from the NLL objective. Note that both objectives have the
same solution. Pragmatically, the CD objective is a slight modification of the
NLL objective which reframes the problem so that samples from the stationary
distribution are no longer needed.

###  4.4 KL term and approximate CD

Notice that the CD formulation introduces another term called the KL term. It
tries to stop MC chains from diverging (i.e. maximize
$D_{\text{KL}}[q_{\theta}^{t}\|p_{\theta}]$) but only with respect to
$q_{\theta}^{t}$ (i.e. the MC chain). Notice that the divergence term does the
same thing but with respect to $p_{\theta}$.

It is ignored in most implementations of CD because is difficult to compute
and because in some cases it has been shown to have negligible impact. In
practice the loss that is most often optimized is an approximation of the CD
loss where this term is dropped :

$\mathcal{L}_{\text{CD*}}=\mathcal{L}_{+}+\mathcal{L}_{d}\approx\mathcal{L}_{\text{CD}}$

It has been shown that this loss doesn’t follow the gradient of any fixed
function and that it is biased Hinton (2012); Fischer & Igel (2011).

###  4.5 Small note on $t$

The value of $t$ varies a lot between different implementations. The general
idea is that large values of $t$ work better but are more expensive. Original
implementation by Hinton, 2002 uses $t=1$. In implementation we are going to
use $t=10$.

###  4.6 Spurious minima and persistent initializations

Because classical CD starts the MC chains from the data and runs them only for
a few steps $t$, only regions near the data manifold are ever explored (and
optimized). Regions far away from the data manifold—that have vanishingly low
likelihood under the data—can therefore be assigned very high probability
under the model. These points are called spurious minimas.

Persistent CD Tieleman (2008) is a variation of the classical algorithm that
tries to alleviate this problem. At the beginning of training, the MC chains
are initialized on the data points like as usual. However, for the subsequent
updates, the initial points are taken to be the result of the previous MC
chains. The reasoning is that if the learning rate is not too big, the model
doesn’t change that much between updates so that it is reasonable to say that
(at least for a few updates) the multiple chains count as one big chain, i.e.

$p_{\theta_{k}}^{2t}(x)\approx p_{\theta_{k}}^{t}(p_{\theta_{k-1}}^{t}(x))$

In words, running many subsequent chains with persistent initialisation and
updates of the parameter in between is approximately equal to running one big
chain on the last parameter. Of course the more updates you make the less
sensible this approximation becomes. For this reason, it can be a good idea to
reset the chains to the data after a certain number of updates.

In practice, PCD has been shown to work well even if the chains are very
rarely reset Tieleman (2008).

##  5 ConvNet EBMs in practice

In this section, we review two recent results : Implicit Generation and
Modeling with Energy Based Models Du & Mordatch (2019) and Improved
Contrastive Divergence Training of Energy-Based Model Du et al. (2020). They
are deep learning contributions that make the training and sampling of ConvNet
EBMs work in practice. The first paper is one of the first method to
successfully train an EBM on complex images and the second paper improves the
method by introducing additional tricks. See Figure 2 for a overview of their
framework. We start by explaining why the theoretical foundation of their
framework is confused and how it can be clarified. Then we describe and
justify the different tricks they use in a more principled way.

###  5.1 NLL or CD?

As the name of the second paper indicates, the authors interpret their
framework as being an improved version of CD. However, as we will argue,
because of an important design decision they made, their learning procedure is
actually optimizing the NLL objective. This distinction is subtle but
important when reasoning about the algorithm.

####  5.1.1 Initialization of sampling chains

Taking inspiration from persistent CD, Du & Mordatch, 2019 propose to have a
reservoir of images from which to sample the initialization of the chains for
the divergence phase. When a chain finishes, the final image is put back into
the reservoir. The reservoir has a limited size (e.g. 10000 in Du et al.,
2020) and adding a new image to it pushes out the oldest element.

However, at the beginning of training, the buffer is initialized with random
noise images. Moreover, every time the buffer is sampled, the image is
replaced with random noise with a small probability ($1\%$ in Du et al.,
2020).

Because of that and what we discussed in section 4.3, it doesn’t make sense to
call this CD. Indeed, starting the chains on the data and optimizing them to
not diverge is the central point of CD.

The training scheme proposed is much more in line with the NLL objective
stance where the samples are actually supposed to come from the $p_{\theta}$
(whereas they are only expected to come from $q_{\theta}^{t}$ in CD). Indeed,
having a persistent reservoir can be seen as a way to improve mixing :
starting the chains on random noise encourage the samples to cover more of the
energy landscape and starting subsequent chains from previous ending points is
approximately like running very long chains (which is required to mix).

However, the more the images are complex, the more the energy landscape in
pixel space is complex and irregular and the harder it is for chains to mix.

We now discuss each of the tricks used by the authors and justify them as ways
of helping the samples to approach the stationary distribution.

###  5.2 Energy function

The energy function is implemented by a convolutional neural network connected
to a fully connected layer that outputs a scalar value. The specific
architecture of the CNN depends on the complexity of the images in the
dataset. For complex images like faces, state of the art CNN including
residual block and self-attention modules are used.

####  5.2.1 Compositional Multi-scale Generation

Figure 2: Illustration of the overall proposed framework for training EBMs in
Du et al., 2020.

An important architectural component introduced by Du et al., 2020 is the
usage of many CNN at different resolutions. The different CNNs act as
independent energy functions and each take as input the original image
downsampled to certain size (see Figure 2). Then the final energy of an image
is the sum of the CNNs output. Importantly, the sum of many energy function
defines the product of their underlying distribution :

$\displaystyle p_{1}(x)\cdot p_{2}(x)$
$\displaystyle=\frac{e^{-E^{(1)}_{\theta}(x)}}{Z^{(1)}_{\theta}}\cdot\frac{e^{-E^{(2)}_{\phi}(x)}}{Z^{(2)}_{\phi}}$
$\displaystyle=\frac{e^{-(E^{(1)}_{\theta}(x)+E^{(2)}_{\phi}(x))}}{Z^{(1)}_{\theta}Z^{(2)}_{\phi}}$

Therefore, for the energy of an image to be low, only one of the CNN must
assign it low energy. In particular, the low resolution CNN can rapidly learn
to assign low energy to real images (since they are not trying to model all
the details), which means that the sampling will spend more time sensible
regions and help the higher resolution CNN to learn.

####  5.2.2 Spectral normalization

In (Du & Mordatch, 2019), authors sometimes found it useful to use spectral
normalization on the weights of the CNN. This has the effect of controlling
the Lipschitz constant of the energy function, i.e. the maximum absolute value
of the slope of $E_{\theta}(x)$. Doing this makes for a smoother energy
landscape : one where sampling is more stable and doesn’t fall as often in
steep holes.

###  5.3 Sampling

As discussed in section 2.3 the sampling algorithm used is Langevin sampling.
However, only using classic Langevin sampling for images leads to especially
poor mixing (see end of section 3). Essentially, the sparsity of real images
in the space of pixels coupled with the very small steps normally taken by
Langevin sampling means MC chains always get trapped in a local modes.

####  5.3.1 Step size

Even if in theory the step size is supposed to be related to with the variance
of the noise (see equation 2.3), in practice both are fixed independently and
are considered hyperparameters. Notably the step size is sometimes fixed to
high values, between 10 and 1000, while the variance of the noise is usually
set to small values like 0.005. High step size are meant to help the sampling
chains move faster on the energy landscape, however, this has the effect of
making them a poor approximation of the valid MCMC algorithm.

####  5.3.2 Data Augmentation Transitions

In the second paper, the authors propose to periodically perform semantically
meaningful transformations to images so that sampling chains can sometimes
teleport to very far away regions in pixel space while remaining sensible (see
Figure 2). Example of transformations are : shift in color, rescalling,
Gaussian blur, deformation, etc..

###  5.4 Estimation of the KL term

One of the main contributions of Du et al., 2020 is the addition of auxiliary
loss which is meant to estimate the normally ignored KL term. Because of our
discussion in section 5.1.1, this should not be a good idea. However, the
authors made a error in the derivation of this auxiliary loss (flipped a sign)
and this resulted in them minimizing the KL term, when it technically should
be maximized (see section 4.2). As we will see, it turns out that this ”doubly
wrong” auxiliary loss is a helpful training signal in the context of NLL.

####  5.4.1 Optimizing the KL loss

Notice that we can write the gradient of the KL loss as follows :

$\displaystyle\nabla_{\theta}\mathcal{L}_{\text{KL}}$
$\displaystyle=\nabla_{\theta}\Big{(}-D_{\text{KL}}[q_{\theta}^{t}\|p_{\Omega(\theta)}]\Big{)}$
$\displaystyle=\nabla_{\theta}\Big{(}-\mathbb{E}_{q_{\theta}^{t}}[\log
q_{\theta}^{t}]+\mathbb{E}_{q_{\theta}^{t}}[\log p_{\Omega(\theta)}]\Big{)}$
$\displaystyle=\nabla_{\theta}\Big{(}\mathcal{H}(q_{\theta}^{t})+\mathbb{E}_{q_{\theta}^{t}}\left[\log\frac{e^{-E_{\Omega(\theta)}(x)}}{Z_{\Omega(\theta)}}\right]\Big{)}$
$\displaystyle=\nabla_{\theta}\Big{(}\underbrace{\mathcal{H}(q_{\theta}^{t})}_{\mathcal{L}_{\text{ent}}}\underbrace{-\mathbb{E}_{q_{\theta}^{t}}\left[E_{\Omega(\theta)}(x)\right]}_{\mathcal{L}_{\text{opt}}}\Big{)}$

so that the KL loss can be further decomposed in two terms

$\mathcal{L}_{\text{KL}}=\mathcal{L}_{\text{ent}}+\mathcal{L}_{\text{opt}}$

where $\mathcal{L}_{\text{ent}}$ corresponds to minimizing the entropy of
$q_{\theta}^{t}$ and $\mathcal{L}_{\text{opt}}$ corresponds to pushing up on
the points resulting from sampling but with respect to the sampling process
$\bm{q_{\theta}^{t}}$.

However, because of some sign error made in the derivation, the authors
flipped the sign of the objective and therefore optimize the exact opposite.
That is, they maximize the entropy of $q_{\theta}^{t}$ and minimize the energy
of the sampling process. We give more details and explain why and how it makes
sense (even if it is wrongly justified).

####  5.4.2 Minimizing the energy of sampling

To minimize the energy of the MCMC samples with respect to the sampling
process, one must propagate the gradient through the steps of Langevin
sampling and compute second order derivatives. This is a very straightforward
task using modern machine learning programming tools since you can just let
the gradient flow through the sampling process and let automatic
differentiation do the job. However, backpropagating through all the steps is
very costly so the authors propose only backpropagating through a few steps.
They show that backpropagating through the whole chain leads to very similar
results.

####  5.4.3 Maximizing the entropy of sampling

To compute the entropy, the authors propose to use a non-parametric nearest
neighbor entropy estimator Beirlant et al. (1997). The entropy $\mathcal{H}$
of a distribution $p(x)$ can be estimated through a set
$X=x_{1},x_{2},\ldots,x_{n}$ of $n$ different points sampled from $p(x)$ as

$\mathcal{H}(p(x))=\frac{1}{n}\sum_{i=1}^{n}\log(n\cdot\text{NN}(x_{i},X))+O(1)$

where $\text{NN}(x_{i},X)$ is the minimum distance between $x_{i}$ and $X$.
Therefore, we can maximize the entropy $\mathcal{H}(q_{\theta}^{t})$ by
minimizing the loss

$\mathcal{L}_{\text{ent}}^{\prime}=-\mathbb{E}_{q_{\theta}^{t}}[\log\text{NN}(x,B)]$

where $B$ is a set containing some number of past samples from MC chains. By
repelling the MCMC samples from each other, we increase entropy, i.e.
diversity.

####  5.4.4 Why does it work?

In the CD loss, the $\mathcal{L}_{\text{KL}}$ loss was meant to make the
sampling distribution $q_{\theta}^{t}$ diverge from the model distribution
$p_{\theta}$, which makes sense in the context of CD but not so much in the
context of NLL. Indeed, in NLL, for the negative phase to be accurate, we need
the sampling distribution to converge to the model distribution!

From this new point of view, minimizing the energy of the sampling process
makes sense as an auxiliary objective because it encourages the samples to end
up in low points (i.e. which are more likely under $p_{\theta}$). Also,
maximizing the entropy has the effect of diversifying the negative samples,
which could all end up in the same hole and give a very biased samples of
$p_{\theta}$.

Overall, the “doubly wrong” auxiliary loss they introduce has the effect of
encouraging sampling to closely approximate the stationary distribution
$p_{\theta}$.

##  6 Experiments

We perform some experiments to answer the following questions : 1) does actual
CD learning lead to good results? 2) what is the effect of better negative
sampling on the energy landscape? We perform the experiments on MNIST.

###  6.1 How does true CD perform?

We look at what happens if we modify their framework so that the learning
algorithm is true CD. To do so, we make a small modification to the
initialization buffer : we initialize it with real images and every time a
point is sampled from it, it has some chance (e.g. $10\%$) of instead being
picked from the dataset (thus keeping the approximation discussed in 4.6
reasonable). No tricks employed in Du et al., 2020 were used. Moreover, we
tested if the ”right” KL loss could help in that scenario.

We were not able to make a model trained this way generate satisfying images
at test time : while a chain initialized on a real image successfully stays on
it (i.e. it doesn’t diverge, as the objective encouraged), a chain initialized
on random noise always fails to converge to a sensible image (Figure 3). This
indicates that the model didn’t really learn to fit the data distribution but
rather to erect high mountains on the border of the data manifold (as
discussed in section 4.6) and leaving the energy landscape malformed
elsewhere. Moreover adding the right KL loss $\mathcal{L}_{\text{KL}}$ did not
help. However, maybe careful tuning of hyperparameters could lead to better
results.

Figure 3: Top two are sampling chains on a model trained with true CD. The
first starts on random noise and the second starts on a data point. Bottom
image is a chain starting from the same data point but on the EBM trained with
the framework in Du & Mordatch, 2019. All are run for the same number of steps
(10000) and the same noise (0.005).

Overall, we observe that training the EBM with NLL leads to a better energy
landscape. For example the last image in Figure 3 shows that the chain
actually mixes. Also, more pragmatically, if our end goal is to generate new
images starting from random points, then it makes sense to make chains start
from random points at training (i.e. with NLL).

###  6.2 How do the tricks impact the energy landscape?

####  6.2.1 Degeneracy of the energy landscape without the tricks

Notice that if at test-time we increase the magnitude of the noise, the
sampling chains end up in deeper minima. Indeed, noisy gradient descent can
more easily escape shallow minima when the noise is high. We observe that when
we add enough noise, all chains converge very fast to the same digit (or few
digits). This indicate that the energy landscape has a funnel shape and
doesn’t fit $p_{\theta}$ well (all digits should be equally likely).

Figure 4: Sampling chains with varying noise level. Top three have the same
starting point. First two are run for the same number of steps (200) but with
(first) low and (second) high noise. Third and fourth are run with high noise
for many more steps (10000) but the fourth has a different starting point.

For example, in the model illustrated in Figure 4, we observed that with high
noise every chain ended up on a 9 no matter the initialization : suggesting 9
is a global energy minimum. We hypothesize that this is because the 9 shares
the most features with all other digits : e.g. when we lower the energy of a
1,2,4,7 or 8 the CNN implicitly lowers the energy of all 9.

###  6.3 Effect of better negative sampling

We implemented the procedure of data augmentation transition (section 5.3.2)
but with a new type of transformation that is better suited to the MNIST
dataset. We used no other tricks (we found that data augmentation was the most
impactful). We use elastic deformations, transformations that warp the image
in a continuous fashion. In training, every time an image is sampled from the
buffer it is slightly warped. Then when generating images at test-time, we
apply the transformation every 100 step of Langevin sampling. We observe that
this leads to sampling chains that mix much better at test time (i.e. don’t
fall in a deep hole). See Figure 5.

Figure 5: Long sampling chains (20000 steps). First is trained and sampled
with deformation. Second is trained and sampled without deformation. Third has
been trained without deformations but sampled with deformations.

We see from the figure that the deformations don’t just help when generating
images at test time but more importantly help during training and lead to a
better energy landscape. Indeed, at test time, the EBM trained with regular
sampling falls in the inescapable deep hole (in this case a ”1”) even when
sampled with deformations.

This is because during training regular chains don’t have time to reach the
”1” quickly enough and so the negative phase doesn’t balance it properly. More
formally, increasing mixing during training means that the negative phase is
less biased and the learning gradient of better quality.

##  7 Conclusion and future work

In this work, we offer an in-depth presentation of EBMs sampling and training.
Using this knowledge, we showed that (Du & Mordatch, 2019; Du et al., 2020)
have mistaken their training objective for CD while it is in fact NLL. We then
discussed different aspect of their training algorithm in light of this
insight : in particular, we showed that a mathematical error they made ends up
not mattering because of our new interpretation.

We showed –theoretically and experimentally– that in order to improve current
MCMC-based EBM systems, the focus should be on making negative sampling chains
mix better, e.g. the tricks in (Du & Mordatch, 2019; Du et al., 2020). As we
discussed throughout the paper, doing Langevin sampling on a complicated high-
dimensional energy landscape (e.g. pixels) necessarily leads to slow mixing
because of the small steps.

Because of that we think that future research should focus on ways to make
large steps during the sampling process. This can be done by leveraging
knowledge about the structure of the data : for example the Data Augmentation
Transition trick makes big steps on the energy landscape based on known
symmetries of the data (e.g. a crop of a face is still a face). More
generally, we think that operating at higher levels of abstraction –where the
statistical relationships are lower-dimensional and simpler– are the key to
resolve the inescapable problems of high-dimensional sampling. An exciting
recent line of work exploring this idea is the series of papers on GFlowNets
Zhang et al. (2022).

## Acknowledgements

This work was done as a final project in the class IFT 6269 (Probabilistic
Graphical Models) given by Simon Lacoste-Julien in Fall 2021 at Mila.

## References

  * Beirlant et al. (1997) Beirlant, J., Dudewicz, E. J., Györfi, L., and Denes, I.  Nonparametric entropy estimation. an overview.  1997\.
  * Du & Mordatch (2019) Du, Y. and Mordatch, I.  Implicit generation and generalization in energy-based models.  _CoRR_ , abs/1903.08689, 2019.  URL http://arxiv.org/abs/1903.08689.
  * Du et al. (2020) Du, Y., Li, S., Tenenbaum, J. B., and Mordatch, I.  Improved contrastive divergence training of energy based models.  _CoRR_ , abs/2012.01316, 2020.  URL https://arxiv.org/abs/2012.01316.
  * Fischer & Igel (2011) Fischer, A. and Igel, C.  Bounding the bias of contrastive divergence learning.  _Neural Computation_ , 23(3):664–673, 2011.  ISSN 0899-7667.  doi: 10.1162/NECO˙a˙00085.
  * Hinton (2002) Hinton, G. E.  Training products of experts by minimizing contrastive divergence.  _Neural Computation_ , 14(8):1771–1800, 2002\.
  * Hinton (2012) Hinton, G. E.  A practical guide to training restricted boltzmann machines.  In Montavon, G., Orr, G. B., and Müller, K.-R. (eds.), _Neural Networks: Tricks of the Trade (2nd ed.)_ , volume 7700 of _Lecture Notes in Computer Science_ , pp. 599–619. Springer, 2012.  ISBN 978-3-642-35288-1.  URL http://dblp.uni-trier.de/db/series/lncs/lncs7700.html#Hinton12.
  * Huembeli et al. (2021) Huembeli, P., Arrazola, J. M., Killoran, N., Mohseni, M., and Wittek, P.  The physics of energy-based models.  (28), 2021.  URL https://physicsofebm.github.io/.
  * Jaynes (1957) Jaynes, E. T.  Information theory and statistical mechanics.  _Phys. Rev._ , 106(4):620–630, May 1957.  doi: 10.1103/PhysRev.106.620.  URL http://prola.aps.org/abstract/PR/v106/i4/p620_1.
  * Kingma & Ba (2017) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization, 2017.
  * Tieleman (2008) Tieleman, T.  Training restricted boltzmann machines using approximations to the likelihood gradient.  In Cohen, W. W., McCallum, A., and Roweis, S. T. (eds.), _Machine Learning, Proceedings of the Twenty-Fifth International Conference (ICML 2008), Helsinki, Finland, June 5-9, 2008_ , volume 307 of _ACM International Conference Proceeding Series_ , pp. 1064–1071. ACM, 2008.  doi: 10.1145/1390156.1390290.  URL https://doi.org/10.1145/1390156.1390290.
  * van Tulder (2021) van Tulder, G.  Elastic deformations for n-dimensional images (python, scipy, numpy, tensorflow, pytorch), 2021.  URL https://github.com/gvtulder/elasticdeform.
  * Welling & Teh (2011) Welling, M. and Teh, Y. W.  Bayesian learning via stochastic gradient langevin dynamics.  In _Proceedings of the 28th International Conference on International Conference on Machine Learning_ , ICML’11, pp. 681–688, Madison, WI, USA, 2011. Omnipress.  ISBN 9781450306195.
  * Xifara et al. (2014) Xifara, T., Sherlock, C., Livingstone, S., Byrne, S., and Girolami, M.  Langevin diffusions and the metropolis-adjusted langevin algorithm.  _Statistics & Probability Letters_, 91:14–19, Aug 2014\.  ISSN 0167-7152.  doi: 10.1016/j.spl.2014.04.002.  URL http://dx.doi.org/10.1016/j.spl.2014.04.002.
  * Zhang et al. (2022) Zhang, D., Malkin, N., Liu, Z., Volokhova, A., Courville, A., and Bengio, Y.  Generative flow networks for discrete probabilistic modeling, 2022.

##  8 Implementation details

###  8.1 Dataset

We used the classic MNIST dataset because it is very simple and doesn’t
require expensive state-of-the-art CNNs, yet it has enough structure to be
interesting to study. It consists in a total of 70000 $28\times 28$ images of
hand-written digits, preprocessed to they have a nice standard form.

###  8.2 Code

Our code is based on van Tulder (2021), implementing a basic PyTorch version
of the method proposed in the first paper covered Du & Mordatch (2019). We
made substantial modification to the code in order to become familiar with it
and make it more modular. Additionally, we implemented most of the tricks from
Du et al., 2020 in order to test their effectiveness.

###  8.3 CNN architecture

We used the CNN architecture from the code base : a simple 4-layers CNN
producing a latent vector of length 128 which is then passed through two
linear layers to output a single scalar. We experimented with more powerful
architecture but improvement in performance didn’t justify the huge different
in running time.

###  8.4 Optimizer

To optimize our model, we used the Adam optimizer Kingma & Ba (2017) with
learning rate $0.0001$ and standard hyperparameters except for the momentum
which is set to 0. We also used gradient clipping of $0.1$ to help with
training stability.