Hitchhiker’s guide on Energy-Based Models: a comprehensive review on the relation with other generative models, sampling and statistical physics
Davide Carbone
INFN, Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy
and Dipartimento di Scienze Matematiche, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy
davide.carbone@polito.it

The problem of description of data through a mathematical model is very old, being the basis of scientific method. The set of measurements use to fulfill the role that contemporary data scientists now refer to as a dataset. Presently, the model can take the form of an exceedingly complex neural network, but the underlying extrapolation remains akin to P.S. Laplace’s famous deterministic statement[laplace2012philosophical]: "An intellect which at a certain moment would know all forces [data] that set nature in motion […] would be uncertain and the future just like the past could be present before its eyes". One can easily extend this reasoning, asserting that the more data one possesses, the more robust and detailed the model that can be constructed atop them. This leads to enhanced predictions and greater stability concerning unforeseen behaviors.
This line of thought was boosted in the previous century with the advent of automatic calculators, and the velocity of development becomes astounding. For instance, consider the remarkable computational power difference between your smartphone and the computer used for the Apollo program by NASA in the 1960s [comparison]. Hence, the quest for data has become an indispensable aspect of contemporary science.
To delve deeper into this issue, let us construct a historical metaphor. One of the early modern achievements in observational astronomy is Kepler’s laws. The genesis of such results is deeply rooted in a vast collection of observational data amassed by T. Brahe [brahe]. Kepler’s formulation was, in fact, motivated by the necessity to explain these astronomical measurements. In a simplified analogy, we observe the dichotomy between the "model," embodied by Kepler, and the "dataset," represented in this narrative by Brahe. Since the 17th century, these two actors have played equally fundamental roles in the advancement of science, taking turns on the stage with the same importance. Consider, for instance, the pivotal role played by Faraday’s experiments in understanding electromagnetism [al2015birth], long before Maxwell’s laws. Or, conversely, the impact of theory of Relativity[relat] way before its experimental confirmation.
In recent years, particularly during the 2000s, we have witnessed a profound paradigm shift represented by the Big Data Era[bigdata]. Thanks to the aforementioned technological advancements in computer science, the volume of generated scientific (and not) data has dramatically increased, resulting from advancements in simulation and storage capabilities. Furthermore, there has been a growing collection of data on human activities, including images, text, sounds, and more.

Returning to the historical analogy, it is akin to Brahe suddenly providing Kepler with a thousand times the amount of data that the latter was accustomed to. This shift posed a methodological problem in what we now refer to as data science, and this is where the machine learning approach came into play[fradkov2020early]. The models required to process Big Data had already been theoretically studied since the invention of the perceptron[mcculloch1943logical]. Their application was constrained by computational power in the last century, but, as a peculiar example of convergent development, they became the primary tools in the toolbox of data scientists in the 2000s, simultaneously to the appearance of Big Data on the stage.

There is indeed a discontinuity that deserves more attention: the increasing collection of data generated by humans. The term Big Data is sometimes limited to images, sounds, videos, text, and metadata resulting from human activities, not just on the internet. Unlike scientific measurements, having access to an extensive quantity of information produced by humans opened Pandora’s box, prompting the natural question: can we build artificial intelligence by leveraging Big Data? In other words, can we construct a machine capable of generating data as humans do, by training it in some smart way? Data here is to be understood in a broad sense, encompassing new theorems, art pieces, images, videos, and even novels.
Generative models represent, in this sense, the most recent breakthrough in technological advancement towards intelligent-like machines. It is complicated to provide a general definition, and there are already many available from different sources[gene, gene2]. However, if we informally focus on those already known to the general public, such as Generative Pre-Trained Transformers (GPT)[gene3], the common traits of most definitions are few. Firstly, generative models require a substantial amount of data for training, in addition to the selection of a precise architecture, which goes far beyond the original perceptron. Secondly, the training is probably not biologically inspired, i.e., we do not learn through backpropagation[grossberg1987competitive], which is the most commonly used training technique in machine learning. For completeness, it is worth noting that this thesis is still debated in neuroscience[lillicrap2020backpropagation]. Thirdly, a generative model is not necessarily informative about the data distribution; for instance, ChatGPT could achieve astounding results in text generation, but the training machine does not provide knowledge about some general features of text generated by humans.
Returning to the historical metaphor: nowadays, we are able to build "BraheGPT," which can generate and gather new plausible measurements about the orbits of planets in unobserved planetary systems after training on observed data from the solar system. However, it is not Kepler; deductive reasoning is not necessary to generate new data instances, although it remains fundamental to understanding the world. Von Neumann would certainly adapt his famous statement[dyson2004meeting] about overfitting to modern data science, cautioning against the ability to generate examples without a general picture.
Prominent data scientists, such as Yann LeCun, have recently emphasized that the use of interpretable generative models is crucial for achieving a "unified world model for AI capable of planning"[lecun]. This thesis becomes imperative in the realm of computational sciences, where qualitative generation alone is insufficient as a benchmark to evaluate model performance. In sectors like Molecular Dynamics, Biochemistry, and similar fields, the model must convey substantial information about the dataset. The generative models that excel in terms of interpretability, which form the main focus of the present work, are precisely the Energy-Based Models (EBMs). These models offer a unique advantage in their ability to provide insights into the underlying mechanisms of the data they generate. In areas such as Molecular Dynamics and Biochemistry, where understanding the intricate relationships within the dataset is crucial, the interpretability of EBMs stands out.
In adopting EBMs, researchers and practitioners gain not only the capacity to generate high-quality data but also a clearer understanding of the factors influencing the generated outputs. This interpretability is indispensable in domains where the model’s ability to convey meaningful information about the dataset is paramount. As the pursuit of a unified world model for AI continues, the emphasis on interpretable generative models, particularly EBMs, plays a pivotal role in bridging the gap between data generation and comprehensive understanding.

After providing a historical overview of generative models, this section is dedicated to exploring the origin and development of Energy-Based Models (EBMs). As we delve into this discussion, it becomes evident that the theoretical foundation of such generative models exists under different names at the intersection of various fields, including statistical physics, probability theory, computer science, and sampling, among others. In this section, we emphasize a historical perspective to shed light on the evolutionary trajectory of EBMs. While we touch upon the overarching theories, more in-depth theoretical discussions are reserved for subsequent chapters. We believe that this review serves as a valuable resource for readers across diverse fields enabling them to construct a comprehensive understanding of what constitutes an Energy-Based Model by tracing the genesis of this topic.
The first ingredient of the story is the Boltzmann-Gibbs measure, a fundamental concept in statistical mechanics, and has its origins in the works of Ludwig Boltzmann and Josiah Willard Gibbs during the late 19th century. These two influential physicists independently contributed to the development of statistical mechanics, providing a bridge between the microscopic behavior of particles and macroscopic thermodynamic properties.
Ludwig Boltzmann made significant strides in understanding the statistical nature of gases, introducing what is now known as the Boltzmann distribution[boltzmann1868studien]. Boltzmann’s statistical approach, which related the statistical weight of different microscopic configurations to their entropy, laid the groundwork for the probabilistic description of thermodynamic systems.
Josiah Willard Gibbs, in parallel with Boltzmann, extended these ideas to develop the canonical ensemble, introducing what is commonly referred to as the Gibbs measure[gibbs1902elementary]. He provides a mathematical framework for calculating thermodynamic properties based on the statistical distribution of particles in a given system. The Boltzmann-Gibbs measure, which emerged from the synthesis of these ideas, describes the probability distribution of particles in different energy states at thermal equilibrium at temperature $T$. It has become a cornerstone of statistical mechanics, applicable to diverse physical systems, including gases, liquids, and solids. We informally recall its definition: given the state of the system $x\in\Omega$, where $\Omega$ is the so-called phase space, and an energy function $U:\Omega\to\mathbb{R}^{+}$, we can express the associated probability density function

$$\rho(x)\propto e^{-\beta U(x)}$$

(1.1)

where $\beta=1/k_{B}T$, $k_{B}$ being the Boltzmann constant. A detailed mathematical description will be provided in the next chapters.
The analysis of the impact of Boltzmann-Gibbs ensemble on physics would require a full monography per se; for the sake of the present work, we directly advance to 1924, when E. Ising presented his PhD thesis[isi]. The so called Ising model is a fundamental mathematical model in statistical mechanics. It serves as a simplified yet powerful representation of magnetic systems, particularly in understanding the behavior of spins in a lattice $\Lambda$ — for simplicity, we can imagine a graph with $N$ nodes. In the Ising model, each lattice site is associated with a magnetic spin, which can take two possible values, usually denoted as "up" or "down", that is $\Omega=\{-1,1\}^{N}$. The interactions between spins are typically modeled using a simple energy function, namely

$$\displaystyle U_{Is}(x)=-\sum_{\langle ij\rangle}J_{ij}x_{i}x_{j}-\mu\sum_{j}h_{j}x_{j}$$

(1.2)

Let us briefly clarify the notation: $i,j\in\Lambda$ are indexes of sites in the lattice; $\langle ij\rangle$ indicates that the sum is restricted to first neighbours and $J_{ij}$ is the strength of the interaction. The field $h_{i}$ instead individually acts on each site and $\mu$ is just a constant that traditionally corresponds to magnetic moment. In laymen terms, each magnetic spin interacts with its first neighbours and with an external field. The alignment of spins is encouraged.
In considering (1.1) as associated to $U_{Is}$, the primary focus is often on the behavior of the system as a function of temperature. In a nutshell, at high temperatures, thermal fluctuations dominate, and the system exhibits no long-range order. As the temperature decreases, there is a critical point at which the system undergoes a phase transition, leading to spontaneous magnetization and the emergence of long-range order.
For some decades the interest for Ising model and its extensions was confined to physics. The motivation for invoking such a model in the present work is the following: in the 80s a fundamental connection between Ising model and data science manifested through Hopfield networks[hopfield1982neural] and Boltzmann machines[ackley1985learning]. Both can be viewed as an Ising lattice where interactions are not confined to first neighbors. Apart from the initial summation, which, for the former, extends to $\forall i,j\in\Lambda$ rather than just $\langle ij\rangle$, the energy function bears resemblance to (1.2). From a statistical physics standpoint, the distinction between a Hopfield network and Boltzmann machines lies solely in the temperature value.
The purpose of the former is pattern recognition and associative memory tasks. A distinctive feature of Hopfield networks is their proficiency in storing and retrieving patterns through symmetric connections between neurons, that is Ising sites, in the network. In practice, when provided with a set of network configurations $y^{\lambda}\in\Omega$ representing patterns, denoted by $\lambda=1,\dots,n$, one constructs the coupling $J$ as follows:

$$\displaystyle J_{ij}={\frac{1}{n}}\sum_{\lambda=1}^{n}y_{i}^{\lambda}y_{j}^{\lambda}$$

(1.3)

This involves employing the Hebbian rule[hebb1949organization] "neurons wire together if they fire together"[lowel1992selection], but further specifications[storkey1997increasing] are available. This phase is commonly referred to as the training of the network. Subsequently, one can define a retrieval iterative dynamics starting from any configuration $x^{k=0}\in\Omega$, as exemplified by the equation:

$$x^{k+1}=\text{sgn}(Jx^{k}+h)\quad\quad k\in\mathbb{N}$$

(1.4)

Here, $J$ represents the coupling matrix defined element-wise in (1.3), and $h$ is a bias vector that influences the preferences for ’up’ or ’down’. It is noteworthy that in a Hopfield network, there is no use of the Boltzmann-Gibbs ensemble; the objective is to construct a dynamical system with prescribed attractors, which are the minima of $U(x)$ by design.
Boltzmann Machines share the same structure and energy function but the goal extends beyond the mere retrieval of patterns; it is to model their overall distribution. To illustrate this concept, consider a finite set of $n$ natural images of cats and dogs. A meticulously designed Hopfield Network could perfectly retrieve any of these examples. On the contrary, a trained Boltzmann Machine aspires to generate new instances of cats and dogs, capturing, in a sense, the distribution of such images. The objective appears to be on a different level of difficulty: although possibly big, the cardinality of the set of patterns is finite; the number of possible variations of cats and dogs is not. Thus, one can immediately guess why the training and generation phases (n.b. it is no more just a retrieval) are completely different w.r.t. Hopfield Networks. The take home message is the hypothesis that the distribution of the given patterns can be described by a Boltzmann-Gibbs ensemble associated to the energy of the Hopfield Network at temperature $T$.
It is convenient to consider Boltzmann Machines as a specific instance of Energy-Based Models, a term introduced by Hinton et al. [teh2003energy], to describe both training and generation phases. EBMs differ from Boltzmann Machines in the use of a generic parametric energy $U_{\theta}(x)$ instead of the usual choice made for the latter. Here, $\theta\in\Theta$ needs to be selected and trained so that the Boltzmann-Gibbs ensemble $\rho_{\theta}$ associated with $U_{\theta}(x)$ "fits well" the distribution of the given patterns, which we refer to as $\rho_{*}$. After training the EBM, the generative phase involves sampling equilibrium configurations from $\rho_{\theta}$. Specifically, a Boltzmann Machine corresponds to an EBM with the choice of $U(x)$ as the energy of a Hopfield Network and $\theta=J$.
Despite their conceptual simplicity, both training and generation represent fundamental open problems that intersect multiple research fields. In essence, sampling from a Boltzmann-Gibbs ensemble is a challenging task in general, and unfortunately, it is necessary even during the training phase. For this reason, the use of Boltzmann Machines was limited to toy models until the proposal of the Contrastive Divergence algorithm by Hinton [hinton2002training].
This procedure, along with its generalizations, made it possible to apply EBMs to practical problems. Moreover, thanks to the adoption of a deep neural network [xie2016theory] as $U_{\theta}$, the interest towards this class of generative models critically increased and in 2010s the use of EBMs for state-of-the-art tasks became standard. However, all that glitters is not gold. Despite its success in generating high-quality individual samples, the use of Contrastive Divergence is known to be biased. For instance, it could happen that individual images are correctly generated, but ensemble properties as the relative proportion of the two species is incorrect. Although Hinton et al. originally claimed that this bias is generally small [carreira2005contrastive], numerous counterexamples have been shown in the more than 20 years since their original paper. The absence of novel paradigm shifts, coupled with the rise of alternative generative models (e.g., diffusion-based ones [song2020score]), has reduced attention on EBMs and consequently on Boltzmann Machines.
The main purpose of the present work is to review the main features of EBMs, in particular in relation with generative modelling, sampling and Statistical Physics. The idea is to provide a useful orientation guide for people coming from very different communities that would like to investigate such topics; we believe that a resource of this kind could be beneficial to favor further developments in the field of energy-based modelling. We will present instances of fruitful interplay between Physics and generative modelling throughout the next sections.
The structure of the review will be the following:

• •

  in Section 2 we define Energy-Based Models and we highlight the main difficulty related to its training through cross-entropy minimization;

• •

  in Section 3 we present a review of the principal generative models as opposed to EBMs. Then, we conclude the section with a comparative scheme between all the presented models, i.e. EBMs, Generative Adversarial Networks, Variational Auto-Encoders, Normalizing Flow and Diffusion-Based Generative Models;

• •

  in Section 4 we review the main MCMC methods used to sample from a Boltzmann-Gibbs ensemble, hence used in the context of EBM training to generate the necessary sample used to perform gradient descent on cross-entropy;

• •

  in Section 5 we summarize the derivation of Boltzmann-Gibbs equilibrium ensemble, motivating the importance of concepts as free energy in the context of statistical learning;

• •

  in Section 6 we present the state of the art about EBM training, that is Constrastive Divergence algorithm, and we highlight his limitation; in conclusion, we provide some references about recent works trying to overcome the main issues of such procedure.

Typical training methods, e.g. based on the so-called Constrastive Divergence[hinton2002training] and its specifications (see Subsection 6), resort to various approximations to calculate the expectation $\mathbb{E}_{\theta}[\partial_{\theta}U_{\theta}]$. While these approaches have proven successful in many situations, they are prone to training instabilities that limit their applicability. The cross-entropy is more stringent, and therefore better, than objectives like the Fisher divergence used to train other generative models: for example, unlike the latter, it is sensitive to the relative probability weights of modes on $\rho_{*}$ separated by low-density regions [song2021train].

Thanks to such results, an estimate of the log-likelihood can be obtained using the Expectation-Maximization (EM) algorithm[dempster1977maximum]: (E) step corresponds to solve the unconstrained variational problem at fixed $\theta$

$$q_{*}(z)=\operatorname*{arg\,max}_{q\in\mathcal{F}}\mathcal{L}(\theta,q(z);x)$$

(3.9)

while (M) step to maximization of ELBO w.r.t. $\theta$ at fixed $q(z)$. To be precise, the output of the (E) steps is conditioned on $x$, which is $q(z)=q(z|x)$. It can be theoretically proven that under suitable condition such an algorithm converges to the optimum and satisfies the equality in (3.7).
For now there is no evident advantage: solving an explicit variational optimization problem can be unfeasible as the computation of (3.3). But further simplifications are possible: in so-called fixed-form variational inference[honkela2010approximate], the variational family $\mathcal{F}$ is constrained to be any parametric family of PDFs $q_{\lambda}(z|x)$ dependent on $\lambda\in\Lambda$; e.g. for the gaussian family $q_{\lambda}(z|x)=\mathcal{N}(z;\boldsymbol{\mu},\boldsymbol{\Sigma})$ we have $\lambda=\{\boldsymbol{\mu},\boldsymbol{\Sigma}\}$. The advantage is that one can perform the (E) step as optimizing $\lambda$ and not in a function class, and possibly find

$$\lambda^{*}=\operatorname*{arg\,max}_{\lambda}\mathcal{L}(\theta,\lambda;x)$$

(3.10)

Since we have to deal with a dataset of $N$ data point, we rewrite

$$\mathcal{L}(\theta,\lambda;\mathcal{X})=\sum_{i=1}^{N}\mathcal{L}(\theta,\lambda_{i};x_{i})$$

(3.11)

and ideally perform gradient-based optimization routines both in (E) and (M) step. But we immediately notice that optimizing the "local" $\lambda_{i}$ for each sample if $N$ is big is very impractical: for instance, for the gaussian class in dimension $d$ we should update $N$ means and covariance matrices, that is $Nd^{2}(d+1)/2$ scalars.
Thus, a last assumption is necessary to practically train the generative model, leading to the so-called amortized variational inference. It corresponds to assume that there exists a parametric map $f_{\phi}$ such that $\lambda_{i}=f_{\phi}(x_{i})$. In this way, the definitive learning objective for EM algorithm is

$$\mathcal{L}(\theta,\lambda;\mathcal{X})=\sum_{i=1}^{N}\mathcal{L}(\theta,\phi;x_{i})=\sum_{i=1}^{N}\mathbb{E}_{q_{\phi}(z_{i}|x_{i})}[\log\rho_{\theta}(x_{i},z_{i})-\log q_{\phi}(z_{i}|x_{i})]$$

(3.12)

Summarizing this first part, we started from the problem of training the decoder network $\rho_{\theta}(x|z)$ and we had to face the issue of computing the marginal log-likelihood. Thanks to a reformulation of the problem, we could explicit an equivalent objective (3.12). Given $q_{\phi}(z|x)$, that is the approximation of the intractable posterior $\rho_{\theta}(z|x)$, $\mathcal{L}(\theta,\lambda;x)$ can be attacked via EM algorithm, i.e. alternatively optimizing $\theta$ and $\phi$.
VAE approach can be seen as a particular case of amortized variational inference in which $q_{\phi}(z|x)$ is approximated via a neural network, which by analogy with AE is denoted as encoder network. Similarly to the decoder network, a widely used choice is a gaussian, i.e.

$$\rho_{\phi}(z|x)=\mathcal{N}(z,\boldsymbol{\mu}_{\phi}(x),\operatorname{diag}\{\boldsymbol{\sigma}^{2}_{\phi}(x)\})$$

(3.13)

where mean and covariance are modelled by a NN. The proposal to train a VAE[DBLP:journals/corr/KingmaW13] is to perform gradient-based optimization on the joint set of parameters $\{\theta,\phi\}$ with (3.12) as objective. Since the encoder and decoder are in cascade, the joint training can be suboptimal[he2018lagging] with respect to the alternating routine in EM algorithm.
Despite this drawback, using the definition of KL divergence and conditional probability, we rewrite (3.12) as

$$\mathcal{L}(\theta,\lambda;\mathcal{X})=\sum_{i=1}^{N}\mathbb{E}_{q_{\phi}(z_{i}|x_{i})}[\log\rho_{\theta}(x_{i}|z_{i})]-\sum_{i=1}^{N}D_{\text{KL}}[q_{\phi}(z_{i}|x_{i})\parallel\rho(z)]$$

(3.14)

The two summations can be easily interpreted: the first one is related to reconstruction accuracy and measures the fidelity of encoding and decoding chain; the second one is a regularization term that forces the posterior (encoder) to be close to the prior, which is a set of independent gaussians — ideally, each $z$ entry should encode an independent characteristic of the data.
Regarding the actual implementation of a gradient routine, the whole point of ELBO reformulation was the intractability of the marginal likelihood. Thus, we have to ensure to not have the same issue for $\mathcal{L}$. The regularization term has an analytical expression for the usual mentioned choices for $q_{\phi}(z|x)$ and $\rho(z)$ (e.g. if it is the KL divergence between gaussian densities). Thus, the computation of the gradient of that summation w.r.t. to $\theta$ or $\phi$ is not a problem for backpropagation algorithm (n.b. we are dealing with NN). On the other hand, the first summation in analytically intractable: the only possibility is the use of a Monte Carlo estimate using samples $\{z^{(r)_{i}}\}_{r=1}^{R}$ drawn from $q_{\phi}(z_{i}|x_{i})$. Sampling from a gaussian encoder is an easy task, but unfortunately it is not a differentiable operation and it poses an obstacle to perform backpropagation w.r.t. $\phi$. The solution to this last issue is the following reparametrization trick:

$$z^{(r)}_{i}=\boldsymbol{\mu}_{\phi}(x_{i})+\operatorname{diag}\{\boldsymbol{\sigma}^{2}_{\phi}(x)\}^{\frac{1}{2}}\epsilon^{(r)}\quad\quad\epsilon^{(r)}\sim\mathcal{N}(\boldsymbol{0}_{L},\boldsymbol{I}_{L})$$

(3.15)

which allows to effectively compute the gradient w.r.t. $\phi$. The resulting empirical estimate of $\mathcal{L}(\theta,\lambda;\mathcal{X})$ is

$$\hat{\mathcal{L}}(\theta,\lambda;\mathcal{X})=\sum_{i=1}^{N}\frac{1}{R}\sum_{i=r}^{R}\log\rho_{\theta}(x_{i}|z_{i}^{(r)})-\sum_{i=1}^{N}D_{\text{KL}}[q_{\phi}(z_{i}|x_{i})\parallel\rho(z)],$$

(3.16)

which is the objective for the joint optimization of $\theta$ and $\phi$.
After some manipulation, we conclude that VAEs can be trained on log-likelihood objective. The main strength appears to be the ease of generation, since for common choices of encoder and decoder such task reduces to sample from a gaussian distribution. In fact, the main drawbacks[wei2020variations, oussidi2018deep, sengupta2020review] of VAEs lays in the training phase. First of all, VAEs have several hyperparameters (e.g., the choice of prior, a possible imbalanced weighting of the reconstruction and regularization terms) that can significantly impact their performance. Finding the optimal set of hyperparameters can be a challenging task. The assumed simple structure of the latent space in VAEs might not capture the complex dependencies present in the data, limiting the expressiveness of the learned representations. Plus, achieving perfect disentanglement remains a challenge. The latent variables might still be entangled, making it challenging to control specific factors independently. Empirically, it is observed that VAEs sometimes generate blurry samples or suffer from mode collapse, where the model focuses on capturing only a few modes of the data distribution, neglecting others. In general it seems to be an issue related to their limited capacity: they might struggle with capturing complex and high-dimensional data distributions effectively, especially when compared to other generative models.