1. Foundations: Gibbs Distributions, Softmax & Energy-Based Models

1. Introduction: The Bridge Between Statistical Physics and Generative AI

At first glance, classical statistical mechanics and modern generative machine learning appear to address disparate disciplines. The former was developed in the late 19th century by Ludwig Boltzmann, James Clerk Maxwell, and Josiah Willard Gibbs to explain how macroscopic thermodynamic observables—such as temperature, pressure, and heat capacity—emerge deterministically from the stochastic interactions of $N \sim 10^{23}$ microscopic atoms. The latter was developed by computer scientists to model complex, high-dimensional probability distributions over natural language tokens, image pixels, protein structures, and neural network weights.

Yet beneath their domain-specific terminologies, both fields solve the exact same mathematical problem:

In physics, this question leads directly to the Boltzmann-Gibbs distribution $p(\sigma) \propto e^{-\beta E(\sigma)}$. In deep learning, this identical formula appears everywhere:

• The Softmax activation function with temperature $T = 1/\beta$ in autoregressive Transformer language models (GPT-4, Llama 3, Claude 3.7).
• Energy-Based Models (EBMs) and score-based generative models where probability densities are parameterized via unnormalized deep neural networks $p_\theta(x) = \frac{1}{Z_\theta} e^{-E_\theta(x)}$.
• The optimal policy in Reinforcement Learning from Human Feedback (RLHF) and Direct Preference Optimization (DPO), where policy distributions form exact Gibbs measures parameterized by reward functions.
• Modern Discrete Diffusion Models (such as SEDD and LLaDA) operating via discrete score matching over vocabulary spaces.

This chapter establishes the rigorous foundations of statistical physics required for modern deep learning. We build from microstates and macrostates to Jaynes' Maximum Entropy Principle, derive the canonical Gibbs measure from first principles, prove its equivalence to temperature-scaled softmax, analyze the training dynamics of Energy-Based Models, and explore sequence-level residual EBMs and energy-compositional controllable generation.

---

2. Microstates, Macrostates, and Phase Space

A central concept in statistical mechanics is the distinction between the complete microscopic description of a system and its low-dimensional macroscopic observation.

2.1 Discrete State Spaces

Consider a physical system composed of $N$ discrete components or a machine learning model operating on discrete sequences.

Classical Physics: The Ising Spin System

In a lattice of $N$ magnetic dipoles (spins), a microstate $\sigma$ specifies the exact orientation of every individual spin:

$$\sigma = (s_1, s_2, \dots, s_N) \in \{-1, +1\}^N$$

The microstate space $\mathcal{S} = \{-1, +1\}^N$ contains $|\mathcal{S}| = 2^N$ discrete configurations. For a macroscopic sample of matter where $N \sim 10^{23}$, the number of microstates is $2^{10^{23}}$, a quantity vastly exceeding the number of elementary particles in the observable universe ($10^{80}$).

Natural Language Processing: Token Sequences

In an autoregressive Large Language Model with vocabulary $\mathcal{V}$ (where typically $|\mathcal{V}| \approx 32{,}000$ to $128{,}000$), a sequence $\mathbf{y}$ of length $L$ represents a discrete microstate:

$$\mathbf{y} = (y_1, y_2, \dots, y_L) \in \mathcal{V}^L$$

For a context window of $L = 2048$ tokens and $|\mathcal{V}| = 32{,}000$, the discrete configuration space has cardinality:

$$|\mathcal{S}| = |\mathcal{V}|^L = (32{,}000)^{2048} \approx 10^{9{,}229}$$

No computer can enumerate, normalize, or directly evaluate all possible microstates in $\mathcal{S}$.

Macrostates and Multiplicity

A macrostate is defined by a small set of macroscopic observables $M(\sigma)$ that summarize the system without tracking microscopic coordinates:

• In physics: The total energy $E(\sigma) = -\sum_{\langle i,j \rangle} J_{ij} s_i s_j - \sum_i h_i s_i$, or the net magnetization $M(\sigma) = \sum_{i=1}^N s_i$.
• In NLP: The semantic topic, sentiment label $y \in \{+1, -1\}$, syntactic grammaticality, or total sequence cross-entropy loss.

The multiplicity (or discrete density of states) $\Omega(E)$ is the number of distinct microstates $\sigma \in \mathcal{S}$ that produce the identical macroscopic energy $E$:

$$\Omega(E) = \sum_{\sigma \in \mathcal{S}} \delta(E(\sigma) - E)$$

where $\delta(\cdot)$ is the Kronecker delta. Ludwig Boltzmann discovered that the thermodynamic entropy $S$ of an isolated macrostate is proportional to the logarithm of its multiplicity:

$$S(E) = k_B \ln \Omega(E)$$

where $k_B \approx 1.380649 \times 10^{-23} \text{ J/K}$ is Boltzmann's constant (often set to $k_B = 1$ in information theory and machine learning).

---

2.2 Continuous State Spaces and Phase Space

When physical systems or neural representations are continuous, the microstate is represented by coordinates in a continuous manifold.

Classical Phase Space $\Gamma$

For a system of $N$ classical point particles moving in $d$-dimensional Euclidean space, a microstate is a single point $\mathbf{x} = (\mathbf{q}, \mathbf{p})$ in a $2dN$-dimensional phase space $\Gamma = \mathbb{R}^{2dN}$:

$$\mathbf{q} = (\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_N) \in \mathbb{R}^{dN}, \quad \mathbf{p} = (\mathbf{p}_1, \mathbf{p}_2, \dots, \mathbf{p}_N) \in \mathbb{R}^{dN}$$

where $\mathbf{q}_i$ represents the position and $\mathbf{p}_i = m \mathbf{v}_i$ represents the linear momentum of particle $i$.

The differential volume element of phase space, nondimensionalized via Planck's constant $h$ and corrected for quantum indistinguishability of identical particles (the Gibbs factor $1/N!$), is:

$$d\Gamma = \frac{d^{dN} \mathbf{q} \, d^{dN} \mathbf{p}}{h^{dN} N!}$$

The continuous density of states $g(E)$ for a Hamiltonian $H(\mathbf{q}, \mathbf{p})$ is defined via the Dirac delta function:

$$g(E) = \int_{\Gamma} \delta(H(\mathbf{q}, \mathbf{p}) - E) \, d\Gamma$$

Deep Learning: Continuous Latent Spaces and Diffusion Models

In continuous generative modeling (e.g., Score-based Diffusion Models, Continuous Normalizing Flows, and Variational Autoencoders), a microstate $\mathbf{x} \in \mathbb{R}^D$ represents an image, an audio waveform, or a continuous latent embedding $\mathbf{z} \in \mathbb{R}^d$:

• In Image Diffusion (DDPM / Score-SDE): $D = C \times H \times W$ (e.g., $3 \times 512 \times 512 \approx 7.86 \times 10^5$ dimensions).
• In Generative Flow Matching: The continuous trajectory $\mathbf{x}_t \in \mathbb{R}^D$ is governed by a time-dependent vector field $v_\theta(\mathbf{x}_t, t)$.
• The macrostate corresponds to class labels (e.g., "Image of a golden retriever"), style descriptors, or clip embeddings.

---

2.3 The Rosetta Stone: Statistical Mechanics $\longleftrightarrow$ Deep Learning

To navigate the cross-disciplinary literature, we construct the exact translation table connecting physical statistical mechanics to machine learning and computer science concepts:

---

3. The Postulate of Equal A Priori Probabilities & Jaynes' MaxEnt Principle

How do we assign probabilities to microstates when our only empirical knowledge consists of macroscopic expectations?

3.1 The Microcanonical Ensemble & Equal A Priori Probabilities

For a completely isolated physical system with fixed total energy in the interval $[E, E + \delta E]$, total volume $V$, and particle number $N$, we possess zero information distinguishing one compatible microstate from another.

In classical mechanics, this postulate is justified dynamically via Liouville's Theorem (which proves that phase space volume is preserved under Hamiltonian time evolution) and the Ergodic Hypothesis (which asserts that over long timescales, the time average of an observable along a single trajectory equals its phase-space ensemble average):

$$\lim_{\mathcal{T} \to \infty} \frac{1}{\mathcal{T}} \int_0^{\mathcal{T}} A(\mathbf{q}(t), \mathbf{p}(t)) \, dt = \int_{\Gamma} A(\mathbf{q}, \mathbf{p}) \, p_{\text{micro}}(\mathbf{q}, \mathbf{p}) \, d\Gamma$$

---

3.2 Jaynes' Principle of Maximum Entropy (MaxEnt)

In two seminal 1957 papers, Edwin T. Jaynes revolutionized statistical mechanics by reinterpreting it not as a purely physical theory of mechanical trajectories, but as a general framework for optimal statistical inference under incomplete information.

Suppose we have a discrete random variable $X$ taking values in $\mathcal{X} = \{x_1, x_2, \dots, x_K\}$. We do not know the true distribution $p = (p_1, p_2, \dots, p_K)$, but we have measured $M$ macroscopic expectation values:

$$\mathbb{E}_p [f_m(X)] = \sum_{k=1}^K p_k f_m(x_k) = \alpha_m, \quad \forall m \in \{1, 2, \dots, M\}$$

where $M < K$, meaning the system of equations is severely underdetermined with infinitely many valid probability distributions.

---

4. First-Principles Derivation of the Boltzmann/Gibbs Distribution

We now present the rigorous step-by-step mathematical derivation of the Boltzmann/Gibbs distribution using the method of Lagrange multipliers under the Jaynes MaxEnt framework.

4.1 The Constrained Optimization Problem

Let $\mathcal{S} = \{1, 2, \dots, K\}$ be the index set of all microstates. Each microstate $i \in \mathcal{S}$ has an associated microscopic energy $E_i \in \mathbb{R}$.

We seek the probability vector $\mathbf{p} = (p_1, p_2, \dots, p_K)^T$ that solves:

$$\begin{aligned} \max_{\mathbf{p}} \quad & H(\mathbf{p}) = -\sum_{i=1}^K p_i \ln p_i \\ \text{subject to} \quad & C_1: \sum_{i=1}^K p_i = 1 \quad \text{(Total probability axiom)} \\ & C_2: \sum_{i=1}^K p_i E_i = U \quad \text{(Prescribed internal energy expectation $\langle E \rangle = U$)} \\ & C_3: p_i \ge 0, \quad \forall i \in \{1, \dots, K\} \quad \text{(Non-negativity)} \end{aligned}$$

4.2 Step-by-Step Lagrange Multiplier Derivation

Step 1: Construct the Lagrangian Function

We introduce Lagrange multipliers $\lambda_0 \in \mathbb{R}$ for the normalization constraint $C_1$, and $\beta \in \mathbb{R}$ for the energy constraint $C_2$. The Lagrangian $\mathcal{L}(\mathbf{p}, \lambda_0, \beta)$ is:

$$\mathcal{L}(p_1, \dots, p_K, \lambda_0, \beta) = -\sum_{i=1}^K p_i \ln p_i - \lambda_0 \left( \sum_{i=1}^K p_i - 1 \right) - \beta \left( \sum_{i=1}^K p_i E_i - U \right)$$

Step 2: Compute Stationary Points

We take the partial derivative of $\mathcal{L}$ with respect to an arbitrary microstate probability $p_k$:

$$\frac{\partial \mathcal{L}}{\partial p_k} = \frac{\partial}{\partial p_k} \left[ -p_k \ln p_k \right] - \lambda_0 \frac{\partial p_k}{\partial p_k} - \beta E_k \frac{\partial p_k}{\partial p_k}$$

Evaluating the derivatives using the product rule:

$$\frac{\partial}{\partial p_k} \left[ -p_k \ln p_k \right] = -\ln p_k - p_k \cdot \frac{1}{p_k} = -\ln p_k - 1$$

Setting the gradient to zero for stationarity:

$$\frac{\partial \mathcal{L}}{\partial p_k} = - \ln p_k - 1 - \lambda_0 - \beta E_k = 0$$

Step 3: Solve for the Microstate Probability $p_k$

Rearranging the algebraic equation:

$$\ln p_k = -(1 + \lambda_0) - \beta E_k$$

Exponentiating both sides:

$$p_k = \exp\left( -(1 + \lambda_0) \right) \cdot \exp\left( -\beta E_k \right)$$

Notice that because $\exp(\cdot) > 0$, the non-negativity constraint $p_k > 0$ is automatically satisfied for all finite $\beta$ and $E_k$.

Step 4: Eliminate the Multiplier $\lambda_0$ via Normalization

We enforce the total probability constraint $\sum_{k=1}^K p_k = 1$:

$$\sum_{k=1}^K p_k = \sum_{k=1}^K \exp\left( -(1 + \lambda_0) \right) \exp\left( -\beta E_k \right) = 1$$

Factoring out the constant scalar term $\exp(-(1 + \lambda_0))$:

$$\exp\left( -(1 + \lambda_0) \right) \left[ \sum_{k=1}^K \exp\left( -\beta E_k \right) \right] = 1$$

We define the sum over all unnormalized Boltzmann weights as the Canonical Partition Function $Z(\beta)$:

$$Z(\beta) \equiv \sum_{k=1}^K \exp(-\beta E_k)$$

Therefore:

$$\exp(-(1 + \lambda_0)) = \frac{1}{Z(\beta)}$$

Substituting this back into our expression for $p_k$ yields the Boltzmann/Gibbs distribution:

---

4.3 Physical and Mathematical Properties of the Partition Function

The partition function $Z(\beta)$ is not merely an arbitrary normalization constant; it is the cumulant generating function of the energy distribution and encapsulates all thermodynamic properties of the system.

1. Free Energy Connection

In thermodynamics, the Helmholtz Free Energy $F$ represents the maximum usable work obtainable from a closed system at constant temperature:

$$F \equiv U - TS = \langle E \rangle - \frac{1}{\beta} H(p)$$

Using the Gibbs distribution, we compute the entropy:

$$\begin{aligned} H(p) &= -\sum_i p_i \ln\left(\frac{e^{-\beta E_i}}{Z}\right) = -\sum_i p_i (-\beta E_i - \ln Z) \\ &= \beta \sum_i p_i E_i + \ln Z \sum_i p_i = \beta \langle E \rangle + \ln Z \end{aligned}$$

Rearranging gives the fundamental identity connecting statistical mechanics to thermodynamics:

$$-\frac{1}{\beta} \ln Z = \langle E \rangle - \frac{1}{\beta} H(p) = F$$

Hence:

$$F(\beta) = -\frac{1}{\beta} \ln Z(\beta) \iff Z(\beta) = e^{-\beta F(\beta)}$$

2. Expected Energy as the First Log-Derivative

The expected internal energy $U = \langle E \rangle$ is obtained by differentiating $\ln Z$:

$$\frac{\partial \ln Z}{\partial \beta} = \frac{1}{Z} \frac{\partial Z}{\partial \beta} = \frac{1}{Z} \frac{\partial}{\partial \beta}\left[\sum_i e^{-\beta E_i}\right] = \frac{1}{Z} \sum_i (-E_i) e^{-\beta E_i} = - \sum_i E_i p_i = -\langle E \rangle$$

Thus:

$$\langle E \rangle = -\frac{\partial \ln Z(\beta)}{\partial \beta}$$

3. Energy Fluctuations and Heat Capacity

Differentiating a second time yields the energy variance:

$$\begin{aligned} \frac{\partial^2 \ln Z}{\partial \beta^2} &= \frac{\partial}{\partial \beta} \left[ -\frac{1}{Z} \sum_i E_i e^{-\beta E_i} \right] \\ &= \frac{1}{Z^2} \left(\frac{\partial Z}{\partial \beta}\right) \sum_i E_i e^{-\beta E_i} - \frac{1}{Z} \sum_i (-E_i^2) e^{-\beta E_i} \\ &= -\langle E \rangle (-\langle E \rangle) + \sum_i E_i^2 p_i \\ &= \langle E^2 \rangle - \langle E \rangle^2 = \operatorname{Var}_p(E) \end{aligned}$$

In physics, this variance is directly proportional to the heat capacity at constant volume $C_v = \frac{\partial U}{\partial T}$:

$$\operatorname{Var}(E) = \frac{\partial^2 \ln Z}{\partial \beta^2} = -\frac{\partial \langle E \rangle}{\partial \beta} = k_B T^2 C_v$$

---

5. Exact Mathematical Equivalence: Gibbs Distributions & Softmax with Temperature

We now establish the direct, 1-to-1 isomorphism between the Boltzmann/Gibbs measure of statistical physics and the Softmax activation function used across modern neural architectures.

5.1 The Equivalence Map

Let $\mathbf{z} = (z_1, z_2, \dots, z_K)^T \in \mathbb{R}^K$ be a vector of unnormalized network outputs (logits) over $K$ discrete classes (or tokens). Let $T > 0$ denote the scalar temperature parameter.

Define the microscopic energy of state $i$ as the negative logit:

$$E_i \equiv -z_i$$

Define the thermodynamic inverse temperature as:

$$\beta \equiv \frac{1}{T} \quad (\text{setting } k_B = 1)$$

Substituting these definitions directly into the Boltzmann distribution:

$$p_i = \frac{e^{-\beta E_i}}{\sum_{j=1}^K e^{-\beta E_j}} = \frac{e^{-(1/T)(-z_i)}}{\sum_{j=1}^K e^{-(1/T)(-z_j)}} = \frac{\exp(z_i / T)}{\sum_{j=1}^K \exp(z_j / T)} \equiv \operatorname{softmax}\left(\frac{\mathbf{z}}{T}\right)_i$$

The mathematical objects are identical:

1. Softmax is the Gibbs distribution where logits are negative energies.
2. The neural network acts as an energy parameterization engine: higher logit $z_i \iff$ lower energy $E_i \iff$ higher probability $p_i$.
3. The Log-Sum-Exp operator is the negative free energy (times temperature):

$$\operatorname{LSE}(\mathbf{z}) = \ln \left( \sum_{j=1}^K e^{z_j} \right) = \ln Z = -\beta F$$

---

5.2 Asymptotic Temperature Regimes in Generative LLMs

In Large Language Models, when generating the next token $y_t \sim \operatorname{softmax}(\mathbf{z}_t / T)$, tuning $T$ controls the physical thermodynamic state of the generation process:

1. Zero Temperature Limit ($T \to 0^+ \iff \beta \to \infty$): Ground State Condensation

In physics, as $T \to 0$, all thermal fluctuations freeze out, and the system condenses with probability 1 into its ground state (the configuration with minimal energy $E_{\min} = \min_k E_k$):

$$\lim_{T \to 0^+} p_i = \lim_{\beta \to \infty} \frac{e^{-\beta E_i}}{\sum_j e^{-\beta E_j}} = \lim_{\beta \to \infty} \frac{e^{-\beta (E_i - E_{\min})}}{\sum_j e^{-\beta (E_j - E_{\min})}}$$

Let $\mathcal{G} = \{k \mid E_k = E_{\min}\}$ be the set of degenerate ground states. For any state $i \notin \mathcal{G}$, $E_i - E_{\min} > 0$, so $e^{-\beta(E_i - E_{\min})} \to 0$ as $\beta \to \infty$. Thus:

$$\lim_{T \to 0^+} p_i = \begin{cases} \frac{1}{|\mathcal{G}|} & \text{if } i \in \mathcal{G} \iff z_i = \max_j z_j \\ 0 & \text{otherwise} \end{cases}$$

In LLM generation, this corresponds to Greedy Decoding ($\operatorname{argmax}_i z_i$). The entropy collapses to $S \to 0$, eliminating output diversity.

2. Infinite Temperature Limit ($T \to \infty \iff \beta \to 0^+$): Pure Maximum Entropy

In physics, as $T \to \infty$, thermal fluctuations completely overwhelm the energetic potential barrier. Expanding the exponential to first order:

$$e^{-\beta E_i} = 1 - \beta E_i + \mathcal{O}(\beta^2) \xrightarrow{\beta \to 0} 1$$

Therefore:

$$\lim_{T \to \infty} p_i = \frac{1}{\sum_{j=1}^K 1} = \frac{1}{K}$$

In LLMs, $T \to \infty$ produces uniform random sampling across the entire vocabulary $\mathcal{V}$, yielding complete semantic incoherence and maximum entropy $H = \ln |\mathcal{V}|$.

3. Intermediate Sampling Techniques: Truncated Canonical Support

To prevent sampling from the high-energy (low-probability) tail of the Gibbs distribution while preserving local diversity, practitioners use:

• Top-$k$ Sampling: Truncates the state space to the $k$ lowest-energy microstates: $\mathcal{S}_k = \{i \mid \operatorname{rank}(E_i) \le k\}$.
• Top-$p$ (Nucleus) Sampling (Holtzman et al., 2020): Restricts sampling to the minimal set of microstates $\mathcal{S}_p \subset \mathcal{S}$ whose cumulative Gibbs mass satisfies:

$$\sum_{i \in \mathcal{S}_p} p_i \ge p$$

---

6. Energy-Based Models (EBMs) in Machine Learning

Having established that parameterized distributions over configurations $\mathbf{x}$ naturally take the form of a Gibbs measure, we now examine Energy-Based Models (EBMs) in machine learning.

6.1 Formulation of Unnormalized Generative Density

An Energy-Based Model parameterizes a continuous or discrete probability distribution over data $\mathbf{x} \in \mathbb{R}^D$ via an unnormalized scalar energy function $E_\theta(\mathbf{x}): \mathbb{R}^D \to \mathbb{R}$ parameterized by neural network weights $\theta$:

Why EBMs? Architectural Freedom and Compositionality

Unlike other generative frameworks, EBMs impose zero architectural restrictions:

1. Autoregressive Models require causal masking and chain-rule factorization $p(\mathbf{x}) = \prod_{t} p(x_t \mid x_{<t})$, imposing fixed generation orderings.
2. Normalizing Flows require invertible transformations with tractable Jacobian determinants $\det(\frac{\partial f}{\partial \mathbf{x}})$, restricting expressivity.
3. VAEs require tractable approximate posteriors $q_\phi(\mathbf{z} \mid \mathbf{x})$ and suffer from blurry reconstructions due to loose variational bounds.
4. GANs require unstable adversarial min-max optimization ($\min_G \max_D V(D,G)$).
5. EBMs allow $E_\theta(\mathbf{x})$ to be any arbitrary deep neural network mapping inputs to a scalar energy value. Crucially, EBMs allow exact compositionality: combining independent constraints amounts to simply summing their energies:

$$p(\mathbf{x}) \propto e^{-(E_1(\mathbf{x}) + E_2(\mathbf{x}))} = e^{-E_1(\mathbf{x})} e^{-E_2(\mathbf{x})}$$

---

6.2 The Fundamental Challenge: The Intractable Partition Function

The primary challenge of EBMs is that the partition function $Z(\theta) = \int_{\mathbb{R}^D} e^{-E_\theta(\mathbf{x})} d\mathbf{x}$ is a high-dimensional integral. In dimensions $D \ge 100$, computing $Z(\theta)$ numerically via grid quadrature requires $\mathcal{O}(K^D)$ evaluations, which is fundamentally computationally intractable (#P-hard).

---

6.3 Maximum Likelihood Learning and the Positive / Negative Phase

Suppose we are given an empirical dataset $\mathcal{D} = \{\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N\}$ sampled i.i.d. from true data distribution $p_{\text{data}}(\mathbf{x})$. We wish to train the parameters $\theta$ of $E_\theta(\mathbf{x})$ by maximizing the average log-likelihood:

$$\mathcal{L}(\theta) = \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} [\ln p_\theta(\mathbf{x})] = \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} \left[ -E_\theta(\mathbf{x}) - \ln Z(\theta) \right]$$

We now compute the exact gradient $\nabla_\theta \mathcal{L}(\theta)$.

Derivation of the Log-Likelihood Gradient

The gradient decomposes into two terms:

$$\nabla_\theta \mathcal{L}(\theta) = -\mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} [\nabla_\theta E_\theta(\mathbf{x})] - \nabla_\theta \ln Z(\theta)$$

We evaluate the gradient of the log partition function:

$$\begin{aligned} \nabla_\theta \ln Z(\theta) &= \frac{1}{Z(\theta)} \nabla_\theta Z(\theta) \\ &= \frac{1}{Z(\theta)} \nabla_\theta \left[ \int_{\mathcal{X}} e^{-E_\theta(\mathbf{x}')} \, d\mathbf{x}' \right] \\ &= \frac{1}{Z(\theta)} \int_{\mathcal{X}} \nabla_\theta \left[ e^{-E_\theta(\mathbf{x}')} \right] \, d\mathbf{x}' \\ &= \frac{1}{Z(\theta)} \int_{\mathcal{X}} \left( -\nabla_\theta E_\theta(\mathbf{x}') \right) e^{-E_\theta(\mathbf{x}')} \, d\mathbf{x}' \\ &= -\int_{\mathcal{X}} \left[ \frac{e^{-E_\theta(\mathbf{x}')}}{Z(\theta)} \right] \nabla_\theta E_\theta(\mathbf{x}') \, d\mathbf{x}' \\ &= -\int_{\mathcal{X}} p_\theta(\mathbf{x}') \nabla_\theta E_\theta(\mathbf{x}') \, d\mathbf{x}' \\ &= -\mathbb{E}_{\mathbf{x}' \sim p_\theta} [\nabla_\theta E_\theta(\mathbf{x}')] \end{aligned}$$

Substituting this back yields the central theorem of Energy-Based Learning:

Physical and Intuitive Interpretation:

1. The Positive Phase: Takes true observed data points $\mathbf{x} \sim p_{\text{data}}$ and updates $\theta$ to push their energy down ($E_\theta(\mathbf{x}) \downarrow \implies p_\theta(\mathbf{x}) \uparrow$).
2. The Negative Phase: Generates ungrounded "fantasies" or "hallucinations" $\mathbf{x}' \sim p_\theta$ from the current model distribution and updates $\theta$ to pull their energy up ($E_\theta(\mathbf{x}') \uparrow \implies p_\theta(\mathbf{x}') \downarrow$).
3. Stationary Equilibrium: Learning converges when the expectations match:

$$\mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} [\nabla_\theta E_\theta(\mathbf{x})] = \mathbb{E}_{\mathbf{x}' \sim p_\theta} [\nabla_\theta E_\theta(\mathbf{x}')] \iff p_\theta(\mathbf{x}) = p_{\text{data}}(\mathbf{x})$$

---

6.4 Sampling Algorithms: SGLD & Contrastive Divergence

To compute the negative phase gradient $\mathbb{E}_{\mathbf{x}' \sim p_\theta}[\nabla_\theta E_\theta(\mathbf{x}')]$, we must sample configurations $\mathbf{x}' \sim p_\theta(\mathbf{x}) \propto e^{-E_\theta(\mathbf{x})}$.

Stochastic Gradient Langevin Dynamics (SGLD)

In continuous spaces $\mathbf{x} \in \mathbb{R}^D$, we sample from the Gibbs measure using Langevin dynamics (Welling & Teh, 2011; Du & Mordatch, 2019):

$$\mathbf{x}_{k+1} = \mathbf{x}_k - \frac{\epsilon}{2} \nabla_\mathbf{x} E_\theta(\mathbf{x}_k) + \sqrt{\epsilon} \mathbf{\xi}_k, \quad \mathbf{\xi}_k \sim \mathcal{N}(0, I_D)$$

where:

• $-\frac{\epsilon}{2} \nabla_\mathbf{x} E_\theta(\mathbf{x}_k)$ is the deterministic drift pulling the sample toward low-energy basins.
• $\sqrt{\epsilon} \mathbf{\xi}_k$ is the Brownian thermal noise injecting entropy to explore the entire phase space.

As step size $\epsilon \to 0$ and steps $K \to \infty$, the stationary distribution of the Markov chain converges exactly to the Gibbs distribution $p_\theta(\mathbf{x}) \propto e^{-E_\theta(\mathbf{x})}$.

Contrastive Divergence (CD-$k$) and Replay Buffers

Running full SGLD to equilibrium at every gradient step is prohibitively slow. Contrastive Divergence ($CD_k$, Hinton 2002) initializes the Langevin Markov chain at data points $\mathbf{x}_0 \sim p_{\text{data}}$ and runs only $k$ steps ($k \approx 5 - 50$). Modern deep EBMs (Du & Mordatch, 2019) maintain a persistent Replay Buffer of past negative samples to allow chains to mix across training iterations.

---

6.5 Circumventing the Partition Function: Score Matching

Can we train an Energy-Based Model without evaluating or approximating $Z(\theta)$ at all?

Aapo Hyvärinen (2005) introduced Score Matching, which optimizes the vector field of log-density derivatives (the score function):

$$\mathbf{s}_\theta(\mathbf{x}) \equiv \nabla_\mathbf{x} \ln p_\theta(\mathbf{x})$$

Applying the definition of an EBM:

$$\mathbf{s}_\theta(\mathbf{x}) = \nabla_\mathbf{x} \left[ -E_\theta(\mathbf{x}) - \ln Z(\theta) \right] = -\nabla_\mathbf{x} E_\theta(\mathbf{x}) - \underbrace{\nabla_\mathbf{x} \ln Z(\theta)}_{= 0 \text{ (indep. of } \mathbf{x})} = -\nabla_\mathbf{x} E_\theta(\mathbf{x})$$

Explicit Score Matching (ESM) Objective

Score matching minimizes the expected squared Fisher divergence between the model score and true data score:

$$J_{\text{ESM}}(\theta) = \frac{1}{2} \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} \left[ \left\| \nabla_\mathbf{x} \ln p_\theta(\mathbf{x}) - \nabla_\mathbf{x} \ln p_{\text{data}}(\mathbf{x}) \right\|_2^2 \right]$$

Using integration by parts (assuming $p_{\text{data}}(\mathbf{x}) \nabla_\mathbf{x} \ln p_\theta(\mathbf{x}) \to 0$ at the domain boundary), Hyvärinen proved this is equivalent to:

$$J_{\text{SM}}(\theta) = \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} \left[ \frac{1}{2} \left\| \nabla_\mathbf{x} E_\theta(\mathbf{x}) \right\|_2^2 - \operatorname{Tr}\left( \nabla_\mathbf{x}^2 E_\theta(\mathbf{x}) \right) \right]$$

where $\operatorname{Tr}(\nabla_\mathbf{x}^2 E_\theta(\mathbf{x})) = \sum_{d=1}^D \frac{\partial^2 E_\theta(\mathbf{x})}{\partial x_d^2}$ is the trace of the Hessian (the Laplacian $\Delta E_\theta(\mathbf{x})$).

This remarkable identity gave birth to modern Score-based Generative Models and Diffusion SDEs (Song & Ermon, 2019; Song et al., 2021; Ho et al., 2020), which we will explore in detail in Chapter 6.

---

6.6 Noise Contrastive Estimation (NCE)

Another approach to bypass $Z(\theta)$ is Noise Contrastive Estimation (NCE) (Gutmann & Hyvärinen, 2010). NCE treats the unnormalized model as a binary classifier trained to distinguish true data $\mathbf{x} \sim p_{\text{data}}$ from synthetic noise generated by a known noise distribution $\mathbf{y} \sim q(\mathbf{x})$:

$$p(D=1 \mid \mathbf{x}) = \frac{p_\theta(\mathbf{x})}{p_\theta(\mathbf{x}) + \nu q(\mathbf{x})} = \frac{e^{-E_\theta(\mathbf{x}) - c}}{e^{-E_\theta(\mathbf{x}) - c} + \nu q(\mathbf{x})}$$

where the log partition function $c = \ln Z(\theta)$ is simply treated as a single scalar learnable parameter alongside $\theta$. As the number of noise samples $\nu \to \infty$, the parameter estimates $\hat{\theta}$ and $\hat{c}$ converge to their true maximum likelihood values.

---

7. Case Studies in NLP: Residual EBMs & Controllable Energy Generation

How do these statistical mechanics principles manifest in state-of-the-art Natural Language Processing?

7.1 Sequence-Level Residual EBMs for Text (Deng et al., ICLR 2020)

Standard autoregressive language models $p_{\text{AR}}(\mathbf{y}) = \prod_{t=1}^L p(y_t \mid y_{<t})$ decompose text generation into a sequence of local, myopic token predictions. While highly efficient for sampling, autoregressive models struggle with global, sequence-level constraints (e.g., long-range discourse coherence, holistic semantic consistency).

Deng et al. (ICLR 2020) proposed Residual Energy-Based Models for text. Instead of training an EBM from scratch, they define the joint sequence distribution as a product of a pre-trained autoregressive base model $p_{\text{AR}}(\mathbf{y} \mid \mathbf{x})$ and a non-autoregressive residual energy function $E_\theta(\mathbf{x}, \mathbf{y})$:

Why Residual Parameterization Succeeds

1. Division of Labor: $p_{\text{AR}}(\mathbf{y} \mid \mathbf{x})$ captures local grammar, syntax, and short-range transitions. $E_\theta(\mathbf{x}, \mathbf{y})$ (parameterized via bidirectional Transformers like BERT or RoBERTa) inspects the entire sequence simultaneously to evaluate holistic coherence.
2. Tractable Importance Sampling: Because $p_{\text{AR}}$ is already close to the data distribution, we can draw candidate sequences $\mathbf{y}^{(1)}, \dots, \mathbf{y}^{(K)} \sim p_{\text{AR}}(\mathbf{y} \mid \mathbf{x})$ and compute self-normalized importance sampling estimates of the conditional partition function:

$$Z_\theta(\mathbf{x}) = \mathbb{E}_{\mathbf{y}' \sim p_{\text{AR}}} \left[ \exp(-E_\theta(\mathbf{x}, \mathbf{y}')) \right] \approx \frac{1}{K} \sum_{k=1}^K \exp(-E_\theta(\mathbf{x}, \mathbf{y}^{(k)}))$$

Training is conducted stably via Noise Contrastive Estimation (NCE), yielding state-of-the-art results on text summarization and machine translation.

---

7.2 Controllable Generation via Energy Composition: "Mix and Match" (ACL 2022)

In controllable text generation, we wish to sample text $\mathbf{y}$ that satisfies multiple simultaneous attribute constraints:

1. High fluency / grammaticality (base language model $E_{\text{LM}}$).
2. Positive sentiment ($E_{\text{sentiment}}$).
3. Specific topic / domain, e.g. "Quantum Computing" ($E_{\text{topic}}$).
4. Strict safety / anti-toxicity constraint ($E_{\text{safety}}$).

Fine-tuning a separate 70B parameter model for every possible permutation of attributes is computationally prohibitive. Statistical mechanics provides an elegant, zero-shot solution through the composition of energy functions (Product of Experts, Hinton 2002; Mireshghallah et al., ACL 2022):

Discrete Gibbs Sampling / MCMC in Token Space

Because text is discrete, continuous SGLD cannot be applied directly. Instead, "Mix and Match" applies Discrete Gibbs / Metropolis-Hastings Sampling:

1. Initialize with an arbitrary sequence $\mathbf{y}^{(0)} = (y_1, \dots, y_L)$.
2. At step $t$, select a random token position $j \in \{1, \dots, L\}$.
3. Mask token $y_j$ and propose a replacement token $v \in \mathcal{V}$ using a masked language model (e.g. BERT) proposal distribution $q(v \mid \mathbf{y}_{-j})$.
4. Accept the candidate sequence $\mathbf{y}' = (y_1, \dots, v, \dots, y_L)$ with Metropolis-Hastings acceptance probability:

$$\alpha(\mathbf{y} \to \mathbf{y}') = \min\left(1, \, \frac{\exp(-E_{\text{total}}(\mathbf{y}'))}{\exp(-E_{\text{total}}(\mathbf{y}))} \cdot \frac{q(\mathbf{y}_j \mid \mathbf{y}'_{-j})}{q(v \mid \mathbf{y}_{-j})} \right)$$

This allows non-autoregressive, plug-and-play controllable steering of pre-trained models without modifying a single model weight.

---

7.3 Direct Preference Optimization (DPO): An Exact Gibbs Partition Cancellation

A preview of Lecture 2: In Reinforcement Learning from Human Feedback (RLHF), the goal is to find an optimal policy $\pi_\theta(y \mid x)$ maximizing expected human reward $r(x,y)$ regularized by KL-divergence to a reference model $\pi_{\text{ref}}(y \mid x)$:

$$\max_{\pi} \mathbb{E}_{x \sim \mathcal{D}, y \sim \pi(\cdot \mid x)} [r(x,y)] - \beta D_{\text{KL}}(\pi(y \mid x) \,\|\, \pi_{\text{ref}}(y \mid x))$$

Solving this constrained optimization problem analytically (as we will prove in Chapter 2) yields an exact Gibbs distribution:

$$\pi^*(y \mid x) = \frac{1}{Z(x)} \pi_{\text{ref}}(y \mid x) \exp\left( \frac{1}{\beta} r(x,y) \right)$$

where $Z(x) = \sum_{y} \pi_{\text{ref}}(y \mid x) \exp(\frac{1}{\beta} r(x,y))$ is the intractable partition function.

In Direct Preference Optimization (DPO) (Rafailov et al., NeurIPS 2023), taking the logarithm of this exact Gibbs relationship allows expressing the latent reward function directly in terms of policy probabilities:

$$r(x,y) = \beta \ln \frac{\pi^*(y \mid x)}{\pi_{\text{ref}}(y \mid x)} + \beta \ln Z(x)$$

When evaluating the pairwise Bradley-Terry preference probability $p(y_w \succ y_l \mid x) = \sigma(r(x, y_w) - r(x, y_l))$, the partition functions cancel each other out identically:

$$r(x, y_w) - r(x, y_l) = \beta \ln \frac{\pi_\theta(y_w \mid x)}{\pi_{\text{ref}}(y_w \mid x)} - \beta \ln \frac{\pi_\theta(y_l \mid x)}{\pi_{\text{ref}}(y_l \mid x)}$$

This allows optimizing human preferences without training a reward model and without sampling from the negative phase of the Gibbs distribution during training.

---

8. Summary & Key Mathematical Identities

To synthesize the theoretical foundations established in Lecture 1, we summarize the core mathematical formulas and their computational roles:

---

9. Analytical Exercise Problems

To master the mathematical tools derived in this chapter, work through the following complete analytical problems.

Problem 1: Maximum Entropy with Multiple Feature Constraints

Prompt: Let $X$ be a discrete random variable taking values in $\mathcal{X}$. Suppose we measure $M$ independent macroscopic feature functions $\mathbf{f}(X) = (f_1(X), f_2(X), \dots, f_M(X))^T$ such that $\mathbb{E}[f_m(X)] = \alpha_m$ for $m=1, \dots, M$. Prove that the Maximum Entropy distribution is a member of the general exponential family:

$$p_{\boldsymbol{\theta}}(x) = \frac{1}{Z(\boldsymbol{\theta})} \exp\left( \sum_{m=1}^M \theta_m f_m(x) \right)$$

where $\boldsymbol{\theta} = (\theta_1, \dots, \theta_M)^T$ are the Lagrange multipliers.

Solution: We formulate the Lagrangian with multipliers $\lambda_0$ for normalization and $\boldsymbol{\theta} = (\theta_1, \dots, \theta_M)^T$ for the $M$ expectation constraints:

$$\mathcal{L}(p, \lambda_0, \boldsymbol{\theta}) = -\sum_{x} p(x) \ln p(x) - \lambda_0 \left( \sum_x p(x) - 1 \right) + \sum_{m=1}^M \theta_m \left( \sum_x p(x) f_m(x) - \alpha_m \right)$$

Differentiating with respect to $p(x)$ and setting to zero:

$$\frac{\partial \mathcal{L}}{\partial p(x)} = -\ln p(x) - 1 - \lambda_0 + \sum_{m=1}^M \theta_m f_m(x) = 0$$

Solving for $p(x)$:

$$p(x) = \exp(-(1 + \lambda_0)) \exp\left( \sum_{m=1}^M \theta_m f_m(x) \right)$$

Enforcing normalization $\sum_x p(x) = 1$ gives $e^{1+\lambda_0} = \sum_x \exp\left( \sum_{m=1}^M \theta_m f_m(x) \right) \equiv Z(\boldsymbol{\theta})$. Hence:

$$p_{\boldsymbol{\theta}}(x) = \frac{1}{Z(\boldsymbol{\theta})} \exp\left( \boldsymbol{\theta}^T \mathbf{f}(x) \right) \quad \blacksquare$$

---

Problem 2: Score Matching Hessian Trace Identity

Prompt: Prove that for any smooth probability density $p_{\text{data}}(\mathbf{x})$ and model $p_\theta(\mathbf{x}) = \frac{1}{Z_\theta} e^{-E_\theta(\mathbf{x})}$, the expected squared score error:

$$J(\theta) = \frac{1}{2} \int_{\mathbb{R}^D} p_{\text{data}}(\mathbf{x}) \left\| \nabla_\mathbf{x} \ln p_\theta(\mathbf{x}) - \nabla_\mathbf{x} \ln p_{\text{data}}(\mathbf{x}) \right\|_2^2 d\mathbf{x}$$

is equal (up to a constant independent of $\theta$) to:

$$J_{\text{SM}}(\theta) = \int_{\mathbb{R}^D} p_{\text{data}}(\mathbf{x}) \left[ \frac{1}{2} \|\nabla_\mathbf{x} E_\theta(\mathbf{x})\|_2^2 - \operatorname{Tr}(\nabla_\mathbf{x}^2 E_\theta(\mathbf{x})) \right] d\mathbf{x}$$

Solution: Expanding the squared $L_2$ norm inside the integral:

$$\begin{aligned} J(\theta) &= \frac{1}{2} \int p_{\text{data}}(\mathbf{x}) \left( \|\nabla_\mathbf{x} \ln p_\theta(\mathbf{x})\|_2^2 - 2 \langle \nabla_\mathbf{x} \ln p_\theta(\mathbf{x}), \nabla_\mathbf{x} \ln p_{\text{data}}(\mathbf{x}) \rangle + \|\nabla_\mathbf{x} \ln p_{\text{data}}(\mathbf{x})\|_2^2 \right) d\mathbf{x} \end{aligned}$$

The third term $\|\nabla_\mathbf{x} \ln p_{\text{data}}(\mathbf{x})\|_2^2$ does not depend on $\theta$ and can be dropped as a constant $C$.

Now consider the cross term. Notice that $p_{\text{data}}(\mathbf{x}) \nabla_\mathbf{x} \ln p_{\text{data}}(\mathbf{x}) = \nabla_\mathbf{x} p_{\text{data}}(\mathbf{x})$. Thus:

$$I_{\text{cross}} = -\int_{\mathbb{R}^D} \langle \nabla_\mathbf{x} \ln p_\theta(\mathbf{x}), \nabla_\mathbf{x} p_{\text{data}}(\mathbf{x}) \rangle d\mathbf{x} = -\sum_{d=1}^D \int_{\mathbb{R}^D} \frac{\partial \ln p_\theta(\mathbf{x})}{\partial x_d} \frac{\partial p_{\text{data}}(\mathbf{x})}{\partial x_d} d\mathbf{x}$$

Applying multidimensional integration by parts to each coordinate $x_d$:

$$\int_{-\infty}^\infty \frac{\partial \ln p_\theta(\mathbf{x})}{\partial x_d} \frac{\partial p_{\text{data}}(\mathbf{x})}{\partial x_d} dx_d = \left[ p_{\text{data}}(\mathbf{x}) \frac{\partial \ln p_\theta(\mathbf{x})}{\partial x_d} \right]_{-\infty}^\infty - \int_{-\infty}^\infty p_{\text{data}}(\mathbf{x}) \frac{\partial^2 \ln p_\theta(\mathbf{x})}{\partial x_d^2} dx_d$$

Assuming boundary conditions $p_{\text{data}}(\mathbf{x}) \to 0$ as $\|\mathbf{x}\| \to \infty$, the boundary term vanishes. Summing over all dimensions $d \in \{1, \dots, D\}$:

$$I_{\text{cross}} = \int_{\mathbb{R}^D} p_{\text{data}}(\mathbf{x}) \sum_{d=1}^D \frac{\partial^2 \ln p_\theta(\mathbf{x})}{\partial x_d^2} d\mathbf{x} = \int_{\mathbb{R}^D} p_{\text{data}}(\mathbf{x}) \operatorname{Tr}\left( \nabla_\mathbf{x}^2 \ln p_\theta(\mathbf{x}) \right) d\mathbf{x}$$

Substituting $\ln p_\theta(\mathbf{x}) = -E_\theta(\mathbf{x}) - \ln Z_\theta$:

• $\nabla_\mathbf{x} \ln p_\theta(\mathbf{x}) = -\nabla_\mathbf{x} E_\theta(\mathbf{x}) \implies \|\nabla_\mathbf{x} \ln p_\theta(\mathbf{x})\|_2^2 = \|\nabla_\mathbf{x} E_\theta(\mathbf{x})\|_2^2$
• $\nabla_\mathbf{x}^2 \ln p_\theta(\mathbf{x}) = -\nabla_\mathbf{x}^2 E_\theta(\mathbf{x}) \implies \operatorname{Tr}(\nabla_\mathbf{x}^2 \ln p_\theta(\mathbf{x})) = -\operatorname{Tr}(\nabla_\mathbf{x}^2 E_\theta(\mathbf{x}))$

Combining terms gives:

$$J_{\text{SM}}(\theta) = \mathbb{E}_{\mathbf{x} \sim p_{\text{data}}} \left[ \frac{1}{2} \|\nabla_\mathbf{x} E_\theta(\mathbf{x})\|_2^2 - \operatorname{Tr}(\nabla_\mathbf{x}^2 E_\theta(\mathbf{x})) \right] \quad \blacksquare$$

---

Problem 3: Exact Bradley-Terry Partition Function Cancellation in DPO

Prompt: Consider two sequences $y_w$ (preferred) and $y_l$ (dispreferred) conditioned on prompt $x$. Under the Bradley-Terry preference model, the probability that $y_w \succ y_l$ is given by:

$$p(y_w \succ y_l \mid x) = \sigma(r^*(x, y_w) - r^*(x, y_l)) = \frac{1}{1 + e^{-(r^*(x, y_w) - r^*(x, y_l))}}$$

Assume the optimal Gibbs policy is $\pi^*(y \mid x) = \frac{1}{Z(x)} \pi_{\text{ref}}(y \mid x) \exp\left( \frac{1}{\beta} r^*(x, y) \right)$. Show that the Bradley-Terry preference likelihood can be expressed entirely in terms of the ratio of policies $\pi^*$ and $\pi_{\text{ref}}$, without any dependence on the partition function $Z(x)$.

Solution: From the optimal Gibbs policy formula:

$$\pi^*(y \mid x) = \frac{\pi_{\text{ref}}(y \mid x) \exp\left(\frac{1}{\beta} r^*(x,y)\right)}{Z(x)}$$

Taking the natural logarithm of both sides:

$$\ln \pi^*(y \mid x) = \ln \pi_{\text{ref}}(y \mid x) + \frac{1}{\beta} r^*(x,y) - \ln Z(x)$$

Rearranging to isolate the ground-truth reward $r^*(x,y)$:

$$r^*(x,y) = \beta \ln\left(\frac{\pi^*(y \mid x)}{\pi_{\text{ref}}(y \mid x)}\right) + \beta \ln Z(x)$$

Now compute the reward difference between the winning completion $y_w$ and losing completion $y_l$:

$$\begin{aligned} r^*(x, y_w) - r^*(x, y_l) &= \left[ \beta \ln\left(\frac{\pi^*(y_w \mid x)}{\pi_{\text{ref}}(y_w \mid x)}\right) + \beta \ln Z(x) \right] - \left[ \beta \ln\left(\frac{\pi^*(y_l \mid x)}{\pi_{\text{ref}}(y_l \mid x)}\right) + \beta \ln Z(x) \right] \\ &= \beta \ln\left(\frac{\pi^*(y_w \mid x)}{\pi_{\text{ref}}(y_w \mid x)}\right) - \beta \ln\left(\frac{\pi^*(y_l \mid x)}{\pi_{\text{ref}}(y_l \mid x)}\right) + \underbrace{\left( \beta \ln Z(x) - \beta \ln Z(x) \right)}_{= 0} \\ &= \beta \ln \left[ \frac{\pi^*(y_w \mid x)}{\pi_{\text{ref}}(y_w \mid x)} \cdot \frac{\pi_{\text{ref}}(y_l \mid x)}{\pi^*(y_l \mid x)} \right] \end{aligned}$$

Substituting this reward difference directly into the sigmoid function:

$$p(y_w \succ y_l \mid x) = \sigma\left( \beta \ln \frac{\pi^*(y_w \mid x)}{\pi_{\text{ref}}(y_w \mid x)} - \beta \ln \frac{\pi^*(y_l \mid x)}{\pi_{\text{ref}}(y_l \mid x)} \right)$$

The partition function $Z(x)$ has cancelled identically, enabling direct preference optimization without reward modeling or MCMC sampling. $\blacksquare$

---

10. References & Further Reading