Energy-based models, score matching and Langevin dynamics

Before this lesson: lesson 01 (maximum likelihood as a KL divergence) and lesson 05 (critics that learn unnormalized density ratios). Gradient descent (optimization lesson 01). From the probability track: Gaussians, Markov chains and the Metropolis–Hastings algorithm. Stochastic differential equations (SDEs) and the Fokker–Planck equation are explained here as far as we need them; the ML for science track treats them in full. Time: about 100 minutes, plus the notebook.

---

Why this matters #

An energy-based model (EBM) is the most general way to write a density: no invertible layers (as in flows), no ordering of the variables (as in autoregressive models). The ideas of this lesson appear in many places:

• Diffusion models learn the score of noisy data with denoising score matching, and sample with steps that look like Langevin dynamics. NCSN (Song and Ermon 2019) is annealed Langevin dynamics; predictor–corrector samplers add Langevin "corrector" steps.
• Bayesian deep learning uses stochastic gradient Langevin dynamics (SGLD): SGD plus the right amount of noise samples a posterior over weights.
• Physics and chemistry. The Boltzmann distribution $p \propto e^{-U/(k_B T)}$ is an EBM with a known energy, and Langevin dynamics is Brownian dynamics. Slow crossing of energy barriers is the rare-event problem of molecular simulation.
• Inverse problems. A posterior $p(x \mid y) \propto p(y \mid x)\,p(x)$ is an EBM whose energy is the data misfit plus the prior energy. Its score is the sum of two scores, and the unknown $p(y)$ drops out.

1. Predict first #

2. The picture #

2.1 An energy is a landscape #

An energy $E(x)$ is a landscape over the data space. The model density is

Valleys of $E$ are regions of high probability (modes). Hills are regions of low probability.

• $E(x) = x^2/2$ gives $p(x) \propto e^{-x^2/2}$: the standard normal $\mathcal N(0, 1)$, with $Z = \sqrt{2\pi}$.
• With a parameter, $E_\theta(x) = \theta x^2/2$ with $\theta > 0$ gives $Z_\theta = \sqrt{2\pi/\theta}$ and $p_\theta = \mathcal N(0, 1/\theta)$.
• The double well $E(x) = x^4/4 - x^2$ has two valleys at $x = \pm\sqrt 2$ with $E = -1$, and a hill at $x = 0$ with $E = 0$. The density has two modes, and the barrier between them has height $E(0) - E(\pm\sqrt2) = 1$.

Only differences of energy matter. Adding a constant $c$ to $E$ multiplies $e^{-E}$ and $Z$ by the same factor $e^{-c}$, so $p$ does not change. The ratio $p(x)/p(x') = e^{E(x') - E(x)}$ needs no $Z$.

Langevin dynamics is noisy gradient descent. The update

is a gradient-descent step on the energy with step size $h$, plus a random kick. Without the kick, $x$ falls into the nearest valley and stops: you get a mode, not a sample. With the kick of size $\sqrt{2h}$, falling and kicking balance, and the points visit each region in proportion to $e^{-E}$. Picture a ball in the landscape that is shaken by heat.

Interactive figure: loading…

The figure shows two facts that section 3 proves. First, the particles settle to the right density for small $h$. Second, a finite step $h$ gives a slightly wrong density (a bias); a Metropolis–Hastings correction (MALA) removes it.

2.2 Learning pushes the energy down at the data and up at model samples #

How do we fit $E_\theta$ to data? Maximum likelihood (as in lesson 01) makes $\log p_\theta(x)$ large on the data. Section 3.2 shows that its gradient has two parts:

• the positive phase lowers the energy at the data points;
• the negative phase raises the energy at samples from the model.

If the model puts probability where there are no data, its samples land there, and the negative phase raises the energy there. Learning stops when the two pushes balance: the model samples look like the data, as far as the features of $E_\theta$ can tell.

Interactive figure: loading…

2.3 The score field and many noise levels #

The score of a density is the vector field $s(x) = \nabla_x \log p(x)$. It points in the direction where the density grows fastest. For an EBM, $\log p = -E - \log Z$, and $\log Z$ does not depend on $x$, so

So Langevin dynamics needs only the score: $x \leftarrow x + h\,s(x) + \sqrt{2h}\,\xi$. If we can learn the score directly (section 3.7), we never need $Z$ or even an energy.

Noise smooths the data. Let $x = x_0 + \sigma\varepsilon$ with $x_0$ from the data and $\varepsilon \sim \mathcal N(0, I)$. The density $p_\sigma$ of $x$ is the data density blurred with a Gaussian of width $\sigma$. For small $\sigma$ it is close to the data. For large $\sigma$ the modes merge into one wide bump.

The figure uses the 8 Gaussians of this track: means $\mu_k = 2(\cos\frac{2\pi k}{8}, \sin\frac{2\pi k}{8})$, standard deviation $s = 0.15$. After noising, each component has variance $s^2 + \sigma^2$, so the score of $p_\sigma$ is known in closed form. The arrows are the exact score, not a learned one.

Interactive figure: loading…

The figure shows the two problems of a single noise level, and the fix:

1. The score does not see the mixture weights. Near a mode, the score points to that mode, whatever its weight. Particles go to the closest mode and stay there.
2. Particles do not cross between modes. At small σ the region between modes has almost no probability, and Langevin needs an exponentially long time to cross it (section 3.9).
3. The fix: start with large noise. At large σ the modes overlap, particles move freely, and they split according to the weights. As σ shrinks, each particle stays in the mode it chose.

3. Definitions and main results #

3.1 Energy-based models #

An energy-based model is a density of the form

where $E_\theta: \R^d \to \R$ is any function with parameters $\theta$ (often a neural network). The partition function (normalizing constant) $Z_\theta$ must be finite: $E_\theta$ must grow fast enough as $\|x\| \to \infty$, or $x$ must live in a bounded set. For $E_\theta(x) = \theta x^2/2$ this needs $\theta > 0$.

In high dimension, $Z_\theta$ is an integral that we cannot compute. So $\log p_\theta(x)$ is known only up to the unknown constant $\log Z_\theta$. We can compare two points ($p_\theta(x)/p_\theta(x')$), but we cannot evaluate the likelihood itself.

3.2 The gradient of the log-likelihood #

Theorem (EBM log-likelihood gradient). For one data point,

For the average log-likelihood $L(\theta) = \mathbb E_{x \sim p_{\text{data}}}[\log p_\theta(x)]$:

Proof.

1. Take the log: $\log p_\theta(x) = -E_\theta(x) - \log Z_\theta$.
2. Differentiate $\log Z_\theta$. Use $\nabla_\theta \log Z = \nabla_\theta Z / Z$ and move the gradient inside the integral (this needs a mild regularity condition, as for dominated convergence):

3. Combine steps 1 and 2 for one $x$, then average over the data. $\blacksquare$

Reading the result. Gradient ascent on $L$ moves $\theta$ in the direction of $-\nabla_\theta E_\theta$ at the data, which lowers the energy there, and of $+\nabla_\theta E_\theta$ at model samples, which raises the energy there. The gradient is zero when

Exponential families. If the energy is linear in $\theta$, $E_\theta(x) = \theta^\top \phi(x)$ with fixed features $\phi$, then $\nabla_\theta E_\theta = \phi(x)$ and the condition is moment matching: $\mathbb E_{\text{data}}[\phi] = \mathbb E_{p_\theta}[\phi]$. The Hessian of $\log Z_\theta$ is $\operatorname{Cov}_{p_\theta}[\phi] \succeq 0$, so $L$ is concave and the optimum is unique when the covariance is positive definite. The widget of section 2.2 is this case with $\phi(x) = (x^4/4, -x^2/2)$: at the end, $\overline{x^2}$ and $\overline{x^4}$ of the chains equal those of the data.

In code. Draw model samples $\tilde x_1, \dots, \tilde x_M$ with an MCMC sampler and detach them from the computation graph. Then minimize the surrogate loss

Its gradient is a Monte Carlo estimate of $-\nabla_\theta L$. The value of $\ell$ is not a log-likelihood; only its gradient matters.

3.3 Langevin dynamics #

To get model samples, we need a sampler that uses only $\nabla E$. A stochastic differential equation (SDE) $dX_t = b(X_t)\,dt + \sqrt 2\,dW_t$ describes a point that moves with velocity $b$ and is kicked by Brownian motion $W_t$: over a time $h$, the kick $W_{t+h} - W_t$ is $\mathcal N(0, h\,I)$ and independent of the past. The simplest way to simulate it is the Euler–Maruyama step:

The noise scales with $\sqrt h$, not $h$: the variance of a sum of independent kicks grows linearly with time.

The (overdamped) Langevin SDE for a target $p$ uses $b = \nabla \log p$:

Theorem (stationarity). If $X_0 \sim p$, then $X_t \sim p$ for all $t$: $p$ is a stationary density of the Langevin SDE. Moreover, if $\rho_t$ is the density of $X_t$, then $\mathrm{KL}(\rho_t \,\|\, p)$ never increases.

Why this is true (intuition). The density of $X_t$ changes for two reasons: the drift moves probability uphill toward the modes, and the noise spreads it out (diffusion). At $\rho = p$ the two flows of probability cancel exactly at every point.

The Fokker–Planck argument

1. For $dX = b\,dt + \sqrt 2\,dW$, the density $\rho_t$ of $X_t$ solves the Fokker–Planck equation $\partial_t \rho = -\nabla\cdot(\rho\, b) + \Delta\rho$. The first term moves probability with velocity $b$; the second term (the Laplacian $\Delta = \nabla\cdot\nabla$) spreads it like heat.
2. Write $\Delta\rho = \nabla\cdot(\nabla\rho) = \nabla\cdot(\rho\,\nabla\log\rho)$, because $\nabla \rho = \rho\,\nabla\log\rho$. With $b = \nabla\log p$:

$$\partial_t\rho = \nabla\cdot\big(\rho\,\nabla\log\rho - \rho\,\nabla\log p\big) = \nabla\cdot\Big(\rho\,\nabla\log\frac{\rho}{p}\Big).$$

3. At $\rho = p$ the term $\log(\rho/p)$ is 0, so $\partial_t\rho = 0$: $p$ is stationary.
4. $\frac{d}{dt}\mathrm{KL}(\rho_t\,\|\,p) = \int \partial_t\rho\,\big(\log\frac{\rho}{p} + 1\big)\,dx$. The "$+1$" part integrates to 0, because the total mass stays 1. Insert step 2 and integrate by parts (the boundary terms vanish when $\rho$ decays fast):

$$\frac{d}{dt}\mathrm{KL}(\rho_t\,\|\,p) = -\int \rho\,\Big\|\nabla\log\frac{\rho}{p}\Big\|^2 dx \le 0,$$

with equality only when $\rho \propto p$, that is, $\rho = p$. $\blacksquare$

Only $\nabla \log p$ appears, so the dynamics never need $Z$. The theorem says nothing about speed. For a target with separated modes, the KL can decrease very slowly (section 3.9).

3.4 ULA is biased; MALA removes the bias #

The unadjusted Langevin algorithm (ULA) is the Euler–Maruyama step of the Langevin SDE with a fixed step $h$:

The continuous SDE has stationary density $p$. The discrete chain does not: it has an error of order $h$. A Gaussian target shows this exactly.

Theorem (ULA bias for a Gaussian). For the target $\mathcal N(\mu, \sigma^2)$ in 1D, ULA is stable for $0 < h < 2\sigma^2$, its stationary distribution is Gaussian with mean $\mu$, and its stationary variance is

Proof.

1. The score is $\nabla\log p(x) = -(x - \mu)/\sigma^2$. So $x_{k+1} - \mu = a\,(x_k - \mu) + \sqrt{2h}\,\xi_k$ with $a = 1 - h/\sigma^2$.
2. This is a linear (AR(1)) chain. It forgets its start when $|a| < 1$, that is, $0 < h < 2\sigma^2$.
3. Taking means: $\mathbb E[x_{k+1}] - \mu = a\,(\mathbb E[x_k] - \mu) \to 0$. No bias in the mean.
4. Taking variances, with $\xi_k$ independent of $x_k$: $v_{k+1} = a^2 v_k + 2h$. At the fixed point, $v = a^2 v + 2h$, so

For $\mathcal N(0, 1)$ and $h = 0.1$: $v = 1/0.95 = 1.0526$, about 5% too large. The bias is $\approx h/(2\sigma^2)$ for small $h$. A smaller $h$ reduces the bias but makes the chain slower: $a = 1 - h/\sigma^2$ is close to 1, so successive samples are strongly correlated.

Several dimensions and conditioning

For a target $\mathcal N(\mu, \Sigma)$, the same computation holds in each eigendirection of $\Sigma$, with $\sigma^2$ replaced by the eigenvalue $\lambda_i$. Stability needs $h < 2\lambda_{\min}$, the bias is largest in the narrowest direction, and the widest direction needs about $\lambda_{\max}/h$ steps to mix. So the cost grows with the condition number $\lambda_{\max}/\lambda_{\min}$, as for gradient descent (optimization lesson 02). Example: $\Sigma = \operatorname{diag}(1.75, 0.25)$ and $h = 0.1$ give ULA variances $1.8015$ (+2.9%) and $0.3125$ (+25%); stability needs $h < 0.5$; the wide direction has $a = 0.943$.

MALA. The Metropolis-adjusted Langevin algorithm uses the ULA step as a proposal, $q(x' \mid x) = \mathcal N\big(x';\ x + h\nabla\log p(x),\ 2h\,I\big)$, and accepts $x'$ with probability

If it rejects, the chain stays at $x$. This is the Metropolis–Hastings rule, so detailed balance makes $p$ exactly stationary for every $h$: the bias is gone. Three facts:

• $p(x')/p(x) = e^{E(x) - E(x')}$: $Z$ cancels.
• MALA needs the energy $E$ itself, not only its gradient. A network that outputs only a score gives no acceptance ratio.
• Large steps are often rejected, so the chain moves less. In high dimension, $h$ must shrink like $d^{-1/3}$ to keep a good acceptance rate (Roberts and Rosenthal 1998).

3.5 Contrastive divergence and persistent chains #

The negative phase needs samples from $p_\theta$, and $\theta$ changes at every update. Running a chain to equilibrium after each update is too expensive. Two approximations are common.

• Contrastive divergence (CD-$k$) (Hinton 2002): start the chains at the data batch, run $k$ MCMC steps (often $k = 1$), and use the end points as negative samples. It is cheap but biased: after $k$ steps the chains are still close to the data, not distributed as $p_\theta$. Low-energy regions far from the data (spurious modes) are never visited, so they are never corrected. In general, CD-$k$ is not the gradient of any objective.
• Persistent CD (PCD; Tieleman 2008): keep the chains from one update to the next instead of restarting them at the data. $\theta$ changes slowly, so the chains stay close to the current $p_\theta$ and the bias is smaller. The widget of section 2.2 uses persistent chains.
• Deep EBMs (Du and Mordatch 2019) run tens of Langevin steps per update, started from a replay buffer of earlier samples (about 95%) or from fresh noise (about 5%), and add a small penalty on the size of the energies. In practice the step size and the noise scale are often tuned separately, with much less noise than $\sqrt{2h}$. The chains are then close to noisy gradient descent, and the samples are not exact samples of $p_\theta$.

Two methods avoid negative samples altogether: score matching (section 3.7), and noise-contrastive estimation, which trains a logistic classifier to separate data from a known noise distribution and learns $\log Z$ as one more parameter.

3.6 The score needs no partition function #

Definition. The score of a density $p$ on $\R^d$ is $s(x) = \nabla_x \log p(x)$, a vector field. (Do not confuse it with the score-function estimator of lesson 02, which uses $\nabla_\phi \log q_\phi$, the gradient with respect to parameters.)

For an EBM, $\log p_\theta(x) = -E_\theta(x) - \log Z_\theta$, and $\log Z_\theta$ does not depend on $x$:

More generally, if $p = \tilde p / Z$ for any unnormalized $\tilde p$, then $\nabla_x \log p = \nabla_x \log \tilde p$. Examples: the score of $\mathcal N(\mu, \sigma^2 I)$ is $-(x - \mu)/\sigma^2$; the score of a mixture $\sum_k \pi_k p_k$ is $\sum_k w_k(x)\,\nabla\log p_k(x)$, where $w_k(x) = \pi_k p_k(x)/p(x)$ are the responsibilities (compute them with log-sum-exp).

So a score model $s_\theta(x) \approx \nabla_x \log p(x)$ can be trained and used for Langevin sampling without ever touching $Z$, or even an energy.

3.7 Score matching #

We want $s_\theta \approx \nabla \log p$ from samples of $p$. In this lesson every score-matching loss has the factor $\tfrac12$ (as in Hyvärinen 2005 and Vincent 2011); the minimizers do not depend on this choice.

Explicit score matching (ESM) minimizes the Fisher divergence

We cannot compute it: the target $\nabla\log p$ is unknown. Two tricks remove the target.

Implicit score matching (Hyvärinen 2005). Integration by parts gives

where $C = \tfrac12\mathbb E\|\nabla \log p\|^2$ does not depend on $\theta$. Proof in 1D: expand the square; the cross term is $\mathbb E[s_\theta\,(\log p)'] = \int s_\theta\,p'\,dx = -\int s_\theta'\,p\,dx = -\mathbb E[s_\theta']$, when $s_\theta p \to 0$ at $\pm\infty$. In $d$ dimensions, do this for each coordinate. The loss needs only samples of $p$. But the trace of the Jacobian $\nabla_x s_\theta$ costs $d$ backward passes, which is too slow for images.

Sliced score matching

Sliced score matching (Song, Garg, Shi and Ermon 2019) replaces the trace by a random projection: $\operatorname{tr}(J) = \mathbb E_v[v^\top J v]$ for a random $v$ with $\mathbb E[vv^\top] = I$, the same idea as Hutchinson's estimator in lesson 04. The loss $\mathbb E_v \mathbb E_{x}\big[v^\top \nabla_x s_\theta(x)\, v + \tfrac12 \|s_\theta(x)\|^2\big]$ needs one Jacobian-vector product per $v$, not $d$ backward passes.

Denoising score matching (Vincent 2011). Add Gaussian noise: $x = x_0 + \sigma\varepsilon$, with $x_0 \sim p_{\text{data}}$ and $\varepsilon \sim \mathcal N(0, I)$. The noising kernel is $q(x \mid x_0) = \mathcal N(x;\ x_0,\ \sigma^2 I)$, and its score is known:

The DSM loss regresses the model on this known target:

The noisy point $x$ has the density $p_\sigma(x) = \int p_{\text{data}}(x_0)\,q(x \mid x_0)\,dx_0$, the data density blurred with $\mathcal N(0, \sigma^2 I)$.

Theorem (DSM = ESM + constant). Assume $\mathbb E\|s_\theta(x)\|^2 < \infty$ for $x \sim p_\sigma$. Then

where $C$ does not depend on $\theta$. So DSM trains $s_\theta$ toward the score of the noisy density $p_\sigma$, not of the clean data.

Why this is true (intuition). One noisy point $x$ can come from many clean points $x_0$. Each of them gives a different target $-\varepsilon/\sigma = (x_0 - x)/\sigma^2$. A least-squares regression learns the average target at each $x$. That average is exactly the score of $p_\sigma$ (step 3 below).

Proof.

1. Expand both squared norms with $\|s - y\|^2 = \|s\|^2 - 2\langle s, y\rangle + \|y\|^2$.
2. The terms $\tfrac12\mathbb E\|s_\theta(x)\|^2$ are the same in both losses, because in both $x$ has the density $p_\sigma$.
3. Key identity. Differentiate $p_\sigma(x) = \int p_{\text{data}}(x_0)\,q(x \mid x_0)\,dx_0$ under the integral and use $\nabla q = q\,\nabla\log q$:

Divide by $p_\sigma(x)$. The factor $p_{\text{data}}(x_0)\,q(x \mid x_0)/p_\sigma(x)$ is the posterior $p(x_0 \mid x)$, so

The marginal score is the posterior mean of the conditional score. 4. The cross terms are equal. By step 3 and the tower property, $\mathbb E\langle s_\theta(x), \nabla\log p_\sigma(x)\rangle = \mathbb E\big\langle s_\theta(x), \mathbb E[-\varepsilon/\sigma \mid x]\big\rangle = \mathbb E\langle s_\theta(x), -\varepsilon/\sigma\rangle$. 5. The remaining terms do not contain $\theta$. So $C = \tfrac12\big(\mathbb E\|\varepsilon/\sigma\|^2 - \mathbb E\|\nabla\log p_\sigma(x)\|^2\big) = \tfrac12\big(d/\sigma^2 - \mathbb E\|\nabla\log p_\sigma(x)\|^2\big)$. $\blacksquare$

Consequences.

• Over all functions, the minimizer is the conditional mean of the target: $s^*(x) = \mathbb E[-\varepsilon/\sigma \mid x] = \nabla\log p_\sigma(x)$. This is the general rule "the best least-squares prediction is the conditional mean".
• The loss does not go to 0. Its minimum is $C = \tfrac12\mathbb E\,\|\varepsilon/\sigma - \mathbb E[\varepsilon/\sigma \mid x]\|^2 \ge 0$, the conditional variance of the target. No function of $x$ can predict the part of the noise that $x$ does not reveal. Do not judge a score model by how close its loss is to 0.
• Noise prediction. Multiply the loss by $\sigma^2$: $\sigma^2\|s_\theta + \varepsilon/\sigma\|^2 = \|\sigma s_\theta + \varepsilon\|^2$. So a network $\hat\varepsilon_\theta(x)$ that predicts the added noise gives the score $s_\theta = -\hat\varepsilon_\theta/\sigma$. The DDPM loss of lesson 07 is this loss over many noise levels.

Interactive figure: loading…

3.8 Tweedie's formula #

The DSM optimum also gives a denoiser. From $x_0 = x - \sigma\varepsilon$ and step 3 of the proof:

This is Tweedie's formula (Robbins 1956, Efron 2011). The score at a noisy point, times $\sigma^2$, is the step back to the average clean point. So "learn the score of $p_\sigma$" and "learn the best denoiser at noise level $\sigma$" are the same task. The same formula gives the posterior covariance: $\operatorname{Cov}[x_0 \mid x] = \sigma^2\big(I + \sigma^2\nabla^2\log p_\sigma(x)\big)$.

Check with a Gaussian. Data $\mathcal N(m, s^2)$ in 1D. Then $p_\sigma = \mathcal N(m, s^2 + \sigma^2)$ and Tweedie gives $\mathbb E[x_0 \mid x] = x - \sigma^2 (x - m)/(s^2 + \sigma^2) = (s^2 x + \sigma^2 m)/(s^2 + \sigma^2)$. This is the posterior mean from Gaussian conditioning: a weighted average of the observation and the prior mean.

3.9 Two problems of a single score, and many noise levels #

A score learned from clean data (or from $p_\sigma$ with small $\sigma$), used in Langevin dynamics, fails in three ways (Song and Ermon 2019).

1. Inaccurate scores in low-density regions. The score is trained only where there are data. Between modes and far from the data, the learned score can be wrong, and the chains start exactly there.

2. The score does not see the mixture weights. Let $p = \pi p_1 + (1 - \pi)\,p_2$, where $p_1$ and $p_2$ have disjoint supports. On the support of $p_1$, $\log p = \log\pi + \log p_1$, so

Langevin chains that cannot cross between the supports keep whatever split they started with. Nearly disjoint supports (well-separated modes) behave almost the same.

3. Slow mixing. To leave a mode, a chain must climb a barrier of height $\Delta E = E(\text{saddle}) - E(\text{mode})$. For overdamped Langevin at unit temperature, the mean time to cross grows like $e^{\Delta E}$ (the Arrhenius law). Kramers' formula makes this precise for a 1D double well: the mean time to go from one well to the other is about

Interactive figure: loading…

The fix: many noise levels. Let $p_\sigma = p * \mathcal N(0, \sigma^2 I)$. For large $\sigma$:

• the noised data fill space, so a learned score is accurate where the chains start;
• the modes overlap, the barriers are low, and the chains move freely between modes;
• the score far from the data now depends on the weights, so the chains split correctly.

Example: for $\tfrac12\mathcal N(-4, 0.5^2) + \tfrac12\mathcal N(4, 0.5^2)$, the barrier $\log p_\sigma(4) - \log p_\sigma(0)$ is $31.31$ at $\sigma = 0$, $5.71$ at $\sigma = 1$, $1.19$ at $\sigma = 2$, $0.21$ at $\sigma = 3$, and $0$ at $\sigma = 4$ (one mode).

Annealed Langevin dynamics (NCSN, Song and Ermon 2019) uses noise levels $\sigma_1 > \sigma_2 > \dots > \sigma_L$ (geometric), and one network $s_\theta(x, \sigma)$ trained with DSM at all of them.

1. Start from a broad distribution (uniform or a wide Gaussian).
2. For $i = 1, \dots, L$: set the step $h_i = \epsilon\,\sigma_i^2$ (proportional to the variance at that level) and run $T$ Langevin steps $x \leftarrow x + h_i\,s_\theta(x, \sigma_i) + \sqrt{2h_i}\,\xi$.
3. Each level starts from the samples of the previous level. Return the samples of the last level.

At large $\sigma$ the chains find the right mode weights; small $\sigma$ then sharpens the samples. Diffusion models do the same with a continuous noise level and a better-designed sampler (lessons 07–09); predictor–corrector samplers still add Langevin corrector steps at each level.

Answers to section 1. (1) Larger: $1/(1 - 0.1/2) = 1.0526$, about 5% too large (section 3.4). (2) No: at $\sigma = 0$ the barrier is $31.3$, so the crossing time is of order $e^{31.3} \approx 4 \times 10^{13}$ time units, times a prefactor. Near $+4$ the score is $(4 - x)/0.25$ for any weights (section 3.9). (3) It learns the score of the noisy density $p_\sigma$, not of the clean data; the loss does not go to 0, its minimum is the constant $C > 0$ (section 3.7).

4. Worked example #

(a) The likelihood gradient of a Gaussian EBM. $E_\theta(x) = \theta x^2/2$, so $\nabla_\theta E_\theta(x) = x^2/2$ and $p_\theta = \mathcal N(0, 1/\theta)$. Data: $\{-1, 0.5, 2\}$.

• Second moment of the data: $\overline{x^2} = (1 + 0.25 + 4)/3 = 1.75$.
• Positive phase: $-\overline{x^2}/2 = -0.875$.
• Negative phase at $\theta = 1$: $\mathbb E_{p_\theta}[x^2/2] = 1/(2\theta) = 0.5$.
• $\nabla_\theta L = -0.875 + 0.5 = -0.375 < 0$: decrease $\theta$, that is, make the model wider.
• The gradient $-\overline{x^2}/2 + 1/(2\theta)$ is zero at $\hat\theta = 1/1.75 = 0.5714$. The model variance $1/\hat\theta = 1.75$ equals the second moment of the data: moment matching.

(b) ULA on the fitted model $\mathcal N(0, 1.75)$.

• With $h = 0.7$: $a = 1 - 0.7/1.75 = 0.6$, and the update is $x' = 0.6\,x + \sqrt{1.4}\,\xi = 0.6\,x + 1.1832\,\xi$.
• Stationary variance: $1.4/(1 - 0.36) = 2.1875 = 1.75/0.8$, so 25% too large. Stable for $h < 2 \cdot 1.75 = 3.5$.
• With $h = 0.07$: $a = 0.96$, the variance is $0.14/(1 - 0.9216) = 1.7857$, only 2% too large.
• The price: the integrated autocorrelation time of an AR(1) chain is $(1 + a)/(1 - a)$. It grows from $1.6/0.4 = 4$ steps to $1.96/0.04 = 49$ steps. About twelve times less bias costs about twelve times more steps per independent sample.

(c) One MALA step. Target $\mathcal N(0, 1)$, $h = 0.5$, current point $x = 1$, noise $\xi = 1.5$.

• Proposal: $x' = x - h\,x + \sqrt{2h}\,\xi = 1 - 0.5 + 1 \cdot 1.5 = 2$.
• Target ratio: $\log p(x') - \log p(x) = -\tfrac{4}{2} + \tfrac{1}{2} = -1.5$.
• Proposal densities (both have variance $2h = 1$, so the constants cancel): $\log q(x' \mid x) = -\tfrac12(2 - 0.5)^2 = -1.125$ and $\log q(x \mid x') = -\tfrac12\big(1 - (2 - 1)\big)^2 = 0$.
• $\log\alpha = -1.5 + 0 - (-1.125) = -0.375$, so $\alpha = e^{-0.375} = 0.687$. The move is accepted with probability 0.687.

(d) Score, Tweedie and the DSM floor for Gaussian data. Data $\mathcal N(2, 1)$, noise $\sigma = 1$. Then $p_\sigma = \mathcal N(2, 2)$.

• Score at $x = 4$: $-(4 - 2)/2 = -1$.
• Tweedie: $\mathbb E[x_0 \mid x = 4] = 4 + 1 \cdot (-1) = 3$. Gaussian conditioning gives the same: $(1 \cdot 4 + 1 \cdot 2)/(1 + 1) = 3$.
• Posterior variance: $\sigma^2(1 + \sigma^2 \cdot (-\tfrac12)) = \tfrac12$, which is $s^2\sigma^2/(s^2 + \sigma^2)$ with $s^2 = 1$.
• DSM floor: $\mathbb E[(\nabla\log p_\sigma)^2] = \mathbb E[(x - 2)^2]/4 = 2/4 = \tfrac12$, so $C = \tfrac12(1/\sigma^2 - \tfrac12) = 0.25$. Even the exact score has $J_{\text{DSM}} = 0.25$.

(e) A mixture. Data $\tfrac12\mathcal N(-2, 0.5^2) + \tfrac12\mathcal N(2, 0.5^2)$ (the "narrow" data of the DSM figure), $\sigma = 0.5$. Noising adds $\sigma^2$ to each variance: $p_\sigma = \tfrac12\mathcal N(-2, 0.5) + \tfrac12\mathcal N(2, 0.5)$. At $x = 1$:

• The two component scores are $-(1 - 2)/0.5 = 2$ and $-(1 + 2)/0.5 = -6$.
• The responsibility of the right component: the ratio of the two densities is $e^{-(1-2)^2/1} / e^{-(1+2)^2/1} = e^{8}$, so $w_+ = 1/(1 + e^{-8}) = 0.999665$.
• The score is the responsibility-weighted average: $0.999665 \cdot 2 + 0.000335 \cdot (-6) = 1.99732$.
• Tweedie: $\mathbb E[x_0 \mid x = 1] = 1 + 0.25 \cdot 1.99732 = 1.4993$. This is almost the posterior mean inside the right component alone, $(0.25 \cdot 1 + 0.25 \cdot 2)/0.5 = 1.5$.

5. In code #

NumPy and SciPy check the numbers of sections 3 and 4. The notebook has these cells and the exercises.

import numpy as np
from scipy.special import logsumexp
from scipy.stats import norm

rng = np.random.default_rng(0)

# (a) EBM gradient for E_theta(x) = theta x^2 / 2: positive + negative phase
data = np.array([-1.0, 0.5, 2.0])
theta = 1.0
grad = -np.mean(data**2 / 2) + 1 / (2 * theta)       # E_p[x^2/2] = 1/(2 theta)
print(round(grad, 3), round(1 / np.mean(data**2), 4))   # -0.375 0.5714

# (b) ULA on N(0, 1.75) with h = 0.7: run 20000 chains for 200 steps
s2, h = 1.75, 0.7
x = np.zeros(20_000)
for _ in range(200):
    x = x - h * x / s2 + np.sqrt(2 * h) * rng.standard_normal(x.size)
print(round(x.var(), 3), s2 / (1 - h / (2 * s2)))  # 2.172 2.1875 (sampling error about 0.02)

# (c) one MALA step: target N(0, 1), h = 0.5, from x = 1 to x' = 2
x0, x1, h = 1.0, 2.0, 0.5
log_alpha = (norm.logpdf(x1) - norm.logpdf(x0)
             + norm.logpdf(x0, x1 - h * x1, np.sqrt(2 * h))    # log q(x | x')
             - norm.logpdf(x1, x0 - h * x0, np.sqrt(2 * h)))   # log q(x' | x)
print(round(log_alpha, 3), round(np.exp(log_alpha), 3))         # -0.375 0.687

# (d) DSM on Gaussian data N(2, 1) with sigma = 1: p_sigma = N(2, 2)
sigma = 1.0
x0 = 2 + rng.standard_normal(1_000_000)
eps = rng.standard_normal(x0.size)
xs = x0 + sigma * eps
exact_score = -(xs - 2) / 2
print(round(np.mean(0.5 * (exact_score + eps / sigma) ** 2), 3))   # 0.25: the DSM floor
print(4 + sigma**2 * (-(4 - 2) / 2))                               # Tweedie at x = 4: 3.0

# (e) mixture 0.5 N(-2, 0.25) + 0.5 N(2, 0.25), sigma = 0.5: score at x = 1
def log_p(x, sigma, w=0.5, m=2.0, s=0.5):
    sd = np.sqrt(s**2 + sigma**2)
    return logsumexp([np.log(w) + norm.logpdf(x, -m, sd),
                      np.log(1 - w) + norm.logpdf(x, m, sd)], axis=0)
d = 1e-6
score = (log_p(1 + d, 0.5) - log_p(1 - d, 0.5)) / (2 * d)
print(round(score, 5), round(1 + 0.25 * score, 4))   # 1.99732 1.4993

# (f) barrier of 0.5 N(-4, 0.25) + 0.5 N(4, 0.25) after noising with sigma
grid = np.linspace(0, 8, 80_001)
for sg in [0, 1, 2, 3, 4]:
    lp = log_p(grid, sg, m=4.0)
    print(sg, round(lp.max() - log_p(0.0, sg, m=4.0), 2))
# 0 31.31 / 1 5.71 / 2 1.19 / 3 0.21 / 4 0.0

The mixture log-density uses logsumexp: far from both components, each Gaussian density underflows to 0 in floating point, and the responsibilities would become $0/0$.

6. Where it shows up #

• Score-based diffusion models. NCSN (Song and Ermon 2019) trains $s_\theta(x, \sigma)$ with DSM and samples with annealed Langevin dynamics. Score SDEs and DDPM use the same DSM (noise-prediction) loss with a continuous noise level (lessons 07–09). Predictor–corrector samplers add Langevin corrector steps.
• Denoisers as priors. By Tweedie's formula, any good denoiser gives a score. "Plug-and-play" and "regularization by denoising" methods for image restoration use this, and diffusion posterior sampling (lesson 15) adds a likelihood score: $\nabla_x\log p(x \mid y) = \nabla_x\log p(x) + \nabla_x\log p(y \mid x)$, where $p(y)$ drops out.
• Bayesian deep learning. SGLD (Welling and Teh 2011) adds Langevin noise to mini-batch SGD with a decreasing step size to sample a posterior over network weights.
• Classifiers as EBMs. The logits $f(x)$ of a classifier define an energy $E(x) = -\log\sum_y e^{f(x)_y}$ (JEM, Grathwohl et al. 2020). Such energies are also used for out-of-distribution detection.
• Physics and chemistry. Overdamped Langevin dynamics is Brownian dynamics, and the score of a Boltzmann density is the force divided by $k_B T$. Slow crossing of barriers is the rare-event problem (protein folding, phase transitions). Boltzmann generators use normalizing flows to sample all modes without crossing barriers.
• Data assimilation. With a Gaussian likelihood, the posterior energy is $E(x) = \|y - Ax\|^2/(2s^2) + E_{\text{prior}}(x)$. 4D-Var minimizes this energy (the MAP estimate); Langevin with a learned prior score samples it. Score-based data assimilation (Rozet and Louppe 2023) uses a few Langevin corrector steps per noise level.

7. Common mistakes #

• Noise $\sqrt h$ instead of $\sqrt{2h}$. The update $x + h\nabla\log p + \sqrt h\,\xi$ samples $p^2$ (half the temperature): a Gaussian target comes out with half its variance as $h \to 0$.
• No noise at all. Then it is gradient descent on $E$: it returns a mode, not samples.
• Ignoring the ULA bias. Test a sampler on a Gaussian with known variance. Steps above $2\lambda_{\min}$ diverge.
• MH with a score network. The MALA ratio needs $E(x) - E(x')$; a score field alone does not give it.
• Sign errors in the EBM loss. Minimize $\overline{E(x_{\text{data}})} - \overline{E(\tilde x)}$, with the negative samples $\tilde x$ detached from the graph.
• Trusting short-run chains. CD-1 or 20 Langevin steps from noise do not sample $p_\theta$. The learned energy can work well with its own sampler and still be a poor density.
• Thinking DSM learns the clean score. It learns the score of $p_\sigma$. The clean score is only the limit $\sigma \to 0$, where the targets have variance $d/\sigma^2$ and the regression is hopeless.
• Expecting the DSM loss to reach 0. Its minimum is $C(\sigma) > 0$, which depends on $\sigma$ and on the data. Compare with a closed-form score when you can.
• One chain on a multimodal target. A chain stuck in one mode looks converged. Run many chains from spread-out starts, and anneal the noise level.
• Confusing the two "scores". $\nabla_x\log p(x)$ (this lesson) is a gradient with respect to the data. The score function of statistics and of REINFORCE, $\nabla_\theta\log p_\theta(x)$, is a gradient with respect to the parameters.

8. Check yourself #

1. $E_\theta(x) = \theta|x|$ with $\theta > 0$. Compute $Z_\theta$ and $\nabla_\theta \log Z_\theta$, check the theorem of section 3.2, and give the maximum-likelihood $\theta$.
   Answer

   $Z_\theta = 2\int_0^\infty e^{-\theta x}\,dx = 2/\theta$, so $\nabla_\theta\log Z_\theta = -1/\theta$. The theorem says $\nabla_\theta \log Z_\theta = -\mathbb E_{p_\theta}[|x|]$, and for this Laplace density $\mathbb E|x| = 1/\theta$. ✓ The gradient of $L$ is $-\overline{|x|} + 1/\theta$, which is zero at $\hat\theta = 1/\overline{|x|}$.
2. ULA on $\mathcal N(0, 4)$ with $h = 1$: write the update, the stationary variance and the largest stable step.
   Answer

   $x' = 0.75\,x + \sqrt2\,\xi$. $v = 2/(1 - 0.5625) = 4.571 = 4/(1 - 1/8)$, about 14% too large. Stable for $h < 2\sigma^2 = 8$.
3. Show that $dX = \nabla\log p\,dt + \sqrt{2T}\,dW$ has a stationary density proportional to $p^{1/T}$. Which $T$ does the "noise $\sqrt h$" mistake correspond to?
   Answer

   The Fokker–Planck equation is $\partial_t\rho = \nabla\cdot\big(\rho\,(T\nabla\log\rho - \nabla\log p)\big)$. It is 0 when $T\nabla\log\rho = \nabla\log p$, that is, $\rho \propto p^{1/T}$. Noise $\sqrt h = \sqrt{2Th}$ means $T = \tfrac12$: the chain samples $p^2$.
4. For $p = \mathcal N(0, 1)$ and the model $s_a(x) = -ax$, compute the implicit score-matching loss $\mathbb E[s_a'(x) + \tfrac12 s_a(x)^2]$ and its minimizer.
   Answer

   $s_a' = -a$ and $\tfrac12 s_a^2 = \tfrac12 a^2x^2$, so the loss is $-a + \tfrac12 a^2$ (using $\mathbb E[x^2] = 1$). It is minimal at $a = 1$, the true score $-x$. Check: $J_{\text{ESM}} = \tfrac12\mathbb E[(1 - a)^2x^2] = \tfrac12(1 - a)^2 = (-a + \tfrac12 a^2) + \tfrac12$, so $C = \tfrac12$.
5. Data $\mathcal N(0, 4)$ in 1D and noise $\sigma = 3$. What is the minimum value of the DSM loss (with the factor $\tfrac12$)? Which function reaches it?
   Answer

   $p_\sigma = \mathcal N(0, 13)$ and $\mathbb E[(\nabla\log p_\sigma)^2] = 1/13$. So $C = \tfrac12(1/9 - 1/13) = 2/117 \approx 0.0171$. It is reached by the exact score $s(x) = -x/13$.
6. Why can you not use MALA with the score network of a diffusion model? What could you do instead?
   Answer

   MALA needs $\log p(x') - \log p(x)$, and the network gives only $\nabla\log p$. You could parameterize the model as an energy and take its gradient (then the energy gives the ratio), integrate the score along the segment from $x$ to $x'$ (costly, and a learned field need not be a true gradient), or accept the ULA bias with a small step.
7. For the mixture of section 1, why does annealed Langevin recover the 50/50 weights when plain Langevin started near $+4$ does not?
   Answer

   At $\sigma = 0$ the barrier is 31.3, so chains never cross, and the score near each mode does not contain the weights. At $\sigma \ge 4$ the noised mixture has one mode (barrier 0), so the chains mix and split according to the weights. Lowering $\sigma$ slowly keeps each chain in the mode it chose, so the split survives as the modes separate.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Write the gradient of $\log p_\theta(x)$ for an energy-based model.

   Show answer

2. Write one step of the unadjusted Langevin algorithm.

   Show answer

3. Stationary variance of ULA for the target $\mathcal N(\mu, \sigma^2)$ with step $h$?

   Show answer

4. How does MALA remove the bias of ULA, and what does it need that ULA does not?

   Show answer

5. Why does the score of an EBM not depend on $Z_\theta$?

   Show answer

6. Write the denoising score matching loss and its target.

   Show answer

7. What does DSM learn at its optimum, and why?

   Show answer

8. State Tweedie's formula for $x = x_0 + \sigma\varepsilon$.

   Show answer

9. Why does contrastive divergence (CD-$k$) give a biased gradient?

   Show answer

10. Why does Langevin with many noise levels (annealing) work better than one small noise level?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why are EBMs hard to train?" — $Z_\theta$ is intractable, so the likelihood gradient needs model samples from MCMC; short chains are biased, long chains are expensive, and mixing between modes is slow.
• "What is the score, and why do diffusion models learn it instead of the density?" — $\nabla_x\log p(x)$; it does not contain $Z$, and it can be learned by denoising (DSM).
• "Langevin versus gradient descent?" — The same step plus noise $\sqrt{2h}\,\xi$; the noise turns mode finding into sampling. Without noise you get a mode.
• "What is the bias of ULA and how do you remove it?" — $O(h)$; for $\mathcal N(0, \sigma^2)$ the variance is $\sigma^2/(1 - h/(2\sigma^2))$. Remove it with MALA or a decreasing step size.
• "Prove that denoising score matching works." — The marginal score is the posterior mean of the conditional score (differentiate under the integral); so the cross terms of DSM and ESM agree, and the losses differ by a constant.
• "Why many noise levels?" — Scores are accurate only where noised data have mass; large noise merges modes, so Langevin can cross and find the weights; small noise sharpens.
• "Contrastive divergence in one sentence?" — Approximate the negative phase with $k$ MCMC steps started at the data: cheap, biased, and not the gradient of a fixed objective in general.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution. All of them use NumPy and exact formulas; no network is trained.

1. ULA on $\mathcal N(0, 1)$: measure the variance for several $h$, compare with $1/(1 - h/2)$, then write MALA and check that the bias disappears.
2. Crossing time against barrier height on the double well $a(x^2 - 1)^2$: simulate and compare with the exact formula and Kramers' approximation.
3. Denoising score matching with a model that is linear in fixed features, fitted by least squares on 1D mixture data; compare with the exact noised score.
4. Annealed versus fixed-σ Langevin on the 8 Gaussians with unequal weights, with exact scores: compare the mode proportions.
5. Persistent contrastive divergence for a 1D Gaussian EBM; compare with the closed-form maximum likelihood estimate.

11. Watch and read #

Tutorial 8: Deep Energy-Based Generative Models (Part 1)Phillip Lippe · UvA Deep Learning course · 15 min

A deep EBM trained with contrastive divergence, Langevin sampling and a replay buffer.

UofT GenAI Course -- Lecture 22: Langevin Algorithm and Contrastive Divergence for EBM TrainingAli Bereyhi · University of Toronto · 31 min

The likelihood gradient of an EBM, Langevin sampling and contrastive divergence, step by step.

MIT 6.S184: Flow Matching and Diffusion Models - Lecture 03A - Score Functions (2026)Peter Holderrieth · MIT 6.S184 · 46 min

Scores of noisy distributions, denoising score matching and Langevin dynamics, in the notation of modern diffusion models.

Diffusion and Score-Based Generative ModelsYang Song · MIT CBMM · 1 h 32 min

The author of NCSN explains score matching, the problems of a single noise level and annealed Langevin dynamics.

• Song and Kingma, How to Train Your Energy-Based Models (2021): §2 (EBMs), §3 (maximum likelihood with MCMC, Langevin, CD and persistent CD), §4 (score matching, including denoising and sliced score matching).
• Song and Ermon, Generative Modeling by Estimating Gradients of the Data Distribution (NeurIPS 2019): §2 (score matching, Langevin), §3 (the problems of low density and slow mixing), §4 (noise levels and annealed Langevin, Algorithm 1).
• Yang Song, Generative Modeling by Estimating Gradients of the Data Distribution (blog post): a readable tour from scores to score SDEs.
• Vincent, A Connection Between Score Matching and Denoising Autoencoders (2011): the original DSM = ESM + constant result.
• Hyvärinen, Estimation of Non-Normalized Statistical Models by Score Matching (JMLR 2005): implicit score matching.
• Song, Garg, Shi and Ermon, Sliced Score Matching (UAI 2019).
• UvA Deep Learning tutorial 8, Deep Energy-Based Generative Models: a practical EBM with contrastive divergence and a replay buffer.
• Murphy, Probabilistic Machine Learning: Advanced Topics, free PDF: chapter 24 (energy-based models: MCMC training, score matching, noise-contrastive estimation).