Score matching and Langevin dynamics

The DDPM note arrived at a model by bounding the likelihood of a noising chain and simplifying the bound until it became a regression on noise. This note arrives at the same model from the other side. Instead of a density, learn its gradient, the score; instead of a decoder, sample by following the score with a little noise added, Langevin dynamics. The two routes meet: the network that predicts the noise in DDPM is, up to a known scale, the score of the noised data.

Notation. $p(x)$ is the data density on $\mathbb R^d$ and $s(x) = \nabla_x \log p(x)$ its score; $s_\theta(x)$ is the model. $q_\sigma(\tilde x) = \int p(x)\,\mathcal N(\tilde x; x, \sigma^2 I)\,dx$ is the data after adding Gaussian noise of scale $\sigma$. $\nabla \cdot s$ is the trace of the Jacobian of $s$.

Why the score#

Many models define a density only up to a constant, like the energy-based models of an earlier note: $p_\theta(x) = e^{-E_\theta(x)} / Z_\theta$ with an energy $E_\theta$ and a normalizer $Z_\theta = \int e^{-E_\theta(x)}\,dx$ that cannot be computed in high dimension. Maximum likelihood needs $\log Z_\theta$ and its gradient. The score does not:

$$\nabla_x \log p_\theta(x) = -\nabla_x E_\theta(x) - \underbrace{\nabla_x \log Z_\theta}_{=\,0} .$$

$Z_\theta$ depends on $\theta$ but not on $x$, so it disappears. The score of a Gaussian $\mathcal N(\mu, \sigma^2 I)$ is $-(x - \mu)/\sigma^2$: a vector pointing back to the mean, longer the farther away. In general the score points toward higher density, and it carries all the information about $p$ where $p$ is positive and smooth on a connected region, such as all of $\mathbb R^d$: two such densities with the same score everywhere differ by a constant factor, which normalization fixes. On two separated components the score fixes each piece only up to its own constant, and says nothing about how the mass splits between them; Gaussian smoothing removes the problem, because a noised density is positive on all of $\mathbb R^d$. When the score is parameterized through an energy, $s_\theta = -\nabla_x E_\theta$, fitting it therefore trains an energy-based model without its normalizer and without samples from the model (energy-based models). An unrestricted network $s_\theta: \mathbb R^d \to \mathbb R^d$ is a different object: a vector field whose Jacobian need not be symmetric, so it need not be the gradient of any scalar function. It is an approximate score with no energy behind it, which is how diffusion models use it.

Fitting a score: Hyvärinen’s trick#

The natural loss is the Fisher divergence, the mean squared distance between the two scores under the data:

$$J(\theta) = \tfrac12\, \mathbb E_{p}\big\lVert s_\theta(x) - \nabla_x \log p(x)\big\rVert^2 .$$

It needs the true score, which is unknown. Hyvärinen (2005) removes it by integration by parts. In one dimension, expand the square; the cross term is

$$\mathbb E_p\big[s_\theta(x)\, (\log p)'(x)\big] = \int s_\theta(x)\, p'(x)\,dx = \big[s_\theta\, p\big]_{-\infty}^{\infty} - \int s_\theta'(x)\, p(x)\,dx = -\mathbb E_p\big[s_\theta'(x)\big],$$

using $p\,(\log p)' = p'$ and assuming $s_\theta(x)\,p(x) \to 0$ at infinity. The square of the true score does not depend on $\theta$. In $d$ dimensions the same step per coordinate gives

$$J(\theta) = \mathbb E_p\Big[\nabla \cdot s_\theta(x) + \tfrac12 \lVert s_\theta(x)\rVert^2\Big] + \text{const},$$

an expectation over data that involves only the model.

Worked example. Fit a Gaussian score $s_\theta(x) = -a(x - \mu)$ with $a = 1/\sigma^2$ to the data $\{1, 2, 4, 5\}$. Then $s_\theta' = -a$ and the objective is $\mathbb E\big[\tfrac12 a^2 (x - \mu)^2 - a\big]$. Setting the derivative in $\mu$ to zero gives $\mu = \bar x = 3$; the derivative in $a$ gives $a\,\mathbb E(x - \mu)^2 = 1$, so $\sigma^2 = \frac14(4 + 1 + 1 + 4) = 2.5$. Score matching recovers the sample mean and variance, the maximum likelihood estimates for a Gaussian, without ever touching the normalizing constant.

The cost is the trace term. For a network $s_\theta: \mathbb R^d \to \mathbb R^d$, $\nabla \cdot s_\theta$ needs $d$ backward passes, one per diagonal entry of the Jacobian: fine for $d = 2$, impossible for an image with $d = 3 \cdot 64 \cdot 64$.

Two cheaper objectives#

Sliced score matching (Song et al., 2019) projects onto random directions $v$ with $\mathbb E[v v^\top] = I$:

$$J_{\text{SSM}}(\theta) = \mathbb E_{p}\,\mathbb E_{v}\Big[v^\top \nabla s_\theta(x)\, v + \tfrac12 \big(v^\top s_\theta(x)\big)^2\Big] .$$

In expectation over $v$ this is the objective above, because $\mathbb E[v^\top A v] = \mathrm{tr}(A)$ and $\mathbb E[(v^\top s)^2] = \lVert s\rVert^2$. One Jacobian–vector product replaces $d$ backward passes.

Denoising score matching (Vincent, 2011) avoids input derivatives of the model in the loss, with no Jacobian $\nabla_x s_\theta$ and no trace, by fitting the score of noised data; training still needs the parameter gradient $\nabla_\theta$, as any loss does. Draw $x \sim p$ and perturb it with any kernel $q(\tilde x \mid x)$ whose score is known; the noised data have the density $q(\tilde x) = \int q(\tilde x \mid x)\,p(x)\,dx$. The objective

$$J_{\text{DSM}}(\theta) = \tfrac12\, \mathbb E_{x \sim p,\ \tilde x \sim q(\tilde x \mid x)}\big\lVert s_\theta(\tilde x) - \nabla_{\tilde x} \log q(\tilde x \mid x)\big\rVert^2$$

differs from the Fisher divergence to $q(\tilde x)$ only by a constant (Song and Ermon, 2019, Eq. 2). Expand both squares. The terms $\lVert s_\theta\rVert^2$ have the same expectation, the terms without $\theta$ are constants, and the cross terms agree because $q(\tilde x)\,\nabla \log q(\tilde x) = \nabla q(\tilde x) = \int \nabla_{\tilde x} q(\tilde x \mid x)\,p(x)\,dx$, so

$$\mathbb E_{q(\tilde x)}\big[s_\theta^\top \nabla \log q(\tilde x)\big] = \mathbb E_{p(x)\,q(\tilde x \mid x)}\big[s_\theta^\top \nabla_{\tilde x} \log q(\tilde x \mid x)\big] .$$

The usual kernel is Gaussian, $\tilde x = x + \sigma\varepsilon$ with $\varepsilon \sim \mathcal N(0, I)$, which gives the density $q_\sigma$ of the notation. Its score is $\nabla_{\tilde x} \log \mathcal N(\tilde x; x, \sigma^2 I) = -(\tilde x - x)/\sigma^2 = -\varepsilon/\sigma$, and the objective becomes

$$J_{\text{DSM}}(\theta) = \tfrac12\, \mathbb E_{x,\,\varepsilon}\Big\lVert s_\theta(x + \sigma\varepsilon) + \frac{\varepsilon}{\sigma}\Big\rVert^2 .$$

Another way to see why it works: a squared-error regression learns a conditional mean: the minimizer is $s^*(\tilde x) = \mathbb E[-\varepsilon/\sigma \mid \tilde x]$, and this conditional mean is exactly $\nabla \log q_\sigma(\tilde x)$, a fact known as Tweedie’s formula.

Worked example. Data $x \sim \mathcal N(0, 1)$ and $\sigma = 0.5$. Then $\tilde x \sim \mathcal N(0, 1.25)$, whose score is $-\tilde x / 1.25 = -0.8\,\tilde x$. Directly: $\varepsilon$ and $\tilde x$ are jointly Gaussian with $\mathrm{Cov}(\varepsilon, \tilde x) = \sigma = 0.5$, so $\mathbb E[\varepsilon \mid \tilde x] = 0.5\,\tilde x / 1.25 = 0.4\,\tilde x$ and $\mathbb E[-\varepsilon/\sigma \mid \tilde x] = -0.8\,\tilde x$: the regression target, averaged, is the score of the noised data.

The price of denoising score matching is that it estimates the score of $q_\sigma$, not of $p$. Small $\sigma$ keeps them close but makes the targets $\varepsilon/\sigma$ large and noisy; large $\sigma$ blurs the data. The way out, below, is to use many values of $\sigma$ at once.

Sampling with Langevin dynamics#

Given a score, Langevin dynamics produces samples:

$$x_{k+1} = x_k + \frac{\eta}{2}\, s(x_k) + \sqrt{\eta}\, z_k, \qquad z_k \sim \mathcal N(0, I) .$$

The drift climbs the density, the noise keeps the chain from collapsing onto the modes. For small step sizes $\eta$ and many steps the distribution of $x_k$ approaches $p$; under regularity conditions it becomes exactly $p$ as $\eta \to 0$ and the number of steps grows without bound (Song and Ermon, 2019, §2.2, citing Welling and Teh, 2011). Without the noise term the update is gradient ascent on $\log p$ and ends at a mode, not at a sample.

Worked example. For $p = \mathcal N(0, 1)$, $s(x) = -x$ and the update is $x_{k+1} = (1 - \eta/2)\, x_k + \sqrt\eta\, z_k$. Its stationary variance $v$ solves $v = (1 - \eta/2)^2 v + \eta$, so

$$v = \frac{\eta}{1 - (1 - \eta/2)^2} = \frac{1}{1 - \eta/4} .$$

With $\eta = 0.1$ the chain settles at variance 1.026 instead of 1; with $\eta = 0.5$ at 1.143. The discretization bias is of order $\eta$, and smaller steps buy accuracy with more steps; a Metropolis–Hastings correction would remove the bias, at the cost of evaluating the density, which is what the score was meant to avoid.

Shrinking the step during the run. The bias also vanishes if the step size decreases along the chain, but only if it decreases slowly. Welling and Teh (2011), Eq. 2, require

$$\sum_k \eta_k = \infty, \qquad \sum_k \eta_k^2 < \infty,$$

for example $\eta_k \propto (b + k)^{-\gamma}$ with $\gamma \in (0.5, 1]$. The first condition lets the chain travel any distance from its start; the second makes the steps, and with them the discretization error, go to zero. A geometric decay $\eta_k = \eta_0 \rho^k$ breaks the first condition: the steps add up to $\eta_0/(1 - \rho)$, and the chain freezes before it has mixed. In the Gaussian example the starting point survives in $x_K$ with the factor $\prod_k (1 - \eta_k/2)$. With $\eta_0 = 0.01$ and $\rho = 0.9$ the steps add up to at most $0.1$, and after ten steps the factor is still $0.968$: the chain has barely left its start.

Temperature. Some samplers scale the noise by a temperature $\tau$, $x_{k+1} = x_k + \frac\eta2\, s(x_k) + \tau\sqrt\eta\, z_k$. Since $\frac\eta2 \nabla \log p = \frac{\tau^2\eta}{2} \nabla \log p^{1/\tau^2}$, this is the plain sampler for the density proportional to $p^{1/\tau^2}$ with step $\tau^2\eta$. A temperature below $1$ therefore concentrates the samples near the modes. In the Gaussian example the stationary variance becomes $v = \tau^2/(1 - \eta/4)$: with $\tau = 0.5$ and $\eta = 0.1$ it is $0.256$, close to the variance $0.25$ of $p^4 \propto \mathcal N(0, 1/4)$.

Why one noise level fails#

Two problems appear as soon as the data are more than a blob (Song and Ermon, 2019).

The score is unknown where there are no data. Score matching weights its error by $p(x)$, so in regions of low density the learned score is barely trained, and a chain started there, as every chain started from noise is, follows arbitrary arrows. When the data lie near a low-dimensional manifold, as images do, almost all of space is such a region, and the score of the data itself is not even defined off the manifold.

Langevin dynamics cannot balance separated modes. Take $p = 0.8\,\mathcal N(-4, 1) + 0.2\,\mathcal N(4, 1)$. Near each mode the score is that of one Gaussian and does not depend on the weights; the weights show only in the gap, where the score at $x = 0$ is $4(1 - 2 \cdot 0.8) = -2.4$, but the density there is about $0.00013$ and chains rarely visit. A chain falls into the mode nearest to its start and stays. With the exact score, 4000 chains started uniformly on $[-8, 8]$ end with 0.52 of them in the left mode after 2000 steps of size 0.05, close to the share that started on the left, not the true 0.8.

Noise addresses both problems. The noised density $q_\sigma$ is positive everywhere, so its score is defined everywhere, and at large $\sigma$ the noised training points spread over a much wider region, so the network is trained where chains from noise begin; it is still trained only where noised samples actually land. For large $\sigma$, $q_\sigma$ also fills the gap between the modes, so that the weights are visible to a chain. How much this helps depends on the choices: the noise levels must overlap, each must get enough steps of a suitable size, and a learned score adds its own approximation error. The mixture example below uses the exact score.

Figure

Noise fills the gap between the modes

Noise-conditional scores and annealing#

Song and Ermon (2019) train one network for many noise levels at once. Choose a geometric sequence $\sigma_1 > \sigma_2 > \dots > \sigma_L$, condition the network on the level, and minimize a weighted sum of denoising objectives:

$$\mathcal L(\theta) = \sum_{i=1}^{L} \sigma_i^2\; \tfrac12\, \mathbb E_{x,\,\varepsilon}\Big\lVert s_\theta(x + \sigma_i\varepsilon, \sigma_i) + \frac{\varepsilon}{\sigma_i}\Big\rVert^2 .$$

The weight $\sigma_i^2$ makes the terms comparable: the target $\varepsilon/\sigma_i$ has size $1/\sigma_i$, so each weighted term is of order one. Sampling, annealed Langevin dynamics, runs a few Langevin steps at the largest noise level, where the chain moves freely, then continues at the next level from where it stopped, down to the smallest, with step sizes $\alpha_i = \epsilon\, \sigma_i^2/\sigma_L^2$ that shrink with the noise. For images the paper uses $L = 10$ levels from $\sigma_1 = 1$ to $\sigma_{10} = 0.01$, 100 steps per level and $\epsilon = 2 \cdot 10^{-5}$.

Two details make the training work. The network sees the level. In the NCSN of Song and Ermon every normalization layer has a learned scale and shift for each level, a conditional instance normalization (App. A.1). In a small fully connected network the level can simply select a learned vector that is added to each hidden layer, $h \mapsto W h + b + e_i$. Diffusion networks receive their timestep in the same ways (DDPM note). One level per example. The sum over levels has $L$ terms, and $\sum_i \lambda_i \ell_i(\theta) = L\, \mathbb E_{i \sim \mathrm{Uniform}\{1, \dots, L\}}\big[\lambda_i \ell_i(\theta)\big]$ for the weighted per-level losses $\lambda_i \ell_i$. So the average of the sampled per-level losses, one level drawn per training example, is an unbiased estimate of $\mathcal L(\theta)/L$: the objective up to the constant factor $L$, which scales the gradient but does not move the minimizer (multiply by $L$ to estimate $\mathcal L$ itself), as the random timestep does for DDPM. For example, four geometric levels from $1$ to $0.01$ are $\sigma_i = 100^{-(i-1)/3}$: $1$, $0.215$, $0.046$ and $0.01$. The figure below shows the fields such a network has to learn at these four levels.

Figure

The score field at four noise levels

In the mixture example, annealing through ten levels from $\sigma = 5$ to $\sigma = 0.1$, 100 steps each, with the exact score of each $q_\sigma$, ends with 0.79 of the chains in the left mode, against the true 0.8.

Figure

Without noise, Langevin loses the mode weights

The same model as DDPM#

DDPM trains a network $\varepsilon_\theta(x_t, t)$ to predict the noise in $x_t = \sqrt{\bar\alpha_t}\, x_0 + \sqrt{1 - \bar\alpha_t}\,\varepsilon$ (forward process). The score of this perturbation kernel is

$$\nabla_{x_t} \log q(x_t \mid x_0) = -\frac{x_t - \sqrt{\bar\alpha_t}\, x_0}{1 - \bar\alpha_t} = -\frac{\varepsilon}{\sqrt{1 - \bar\alpha_t}} ,$$

so denoising score matching at level $t$ regresses $s_\theta(x_t, t)$ onto $-\varepsilon/\sqrt{1 - \bar\alpha_t}$. Writing $s_\theta(x_t, t) = -\varepsilon_\theta(x_t, t)/\sqrt{1 - \bar\alpha_t}$ turns the DSM loss into $\lVert\varepsilon_\theta - \varepsilon\rVert^2 / (1 - \bar\alpha_t)$, and the weight $1 - \bar\alpha_t$, the variance of the noise, is the $\sigma^2$ weighting above. The DDPM objective $L_{\text{simple}}$ is noise-conditional denoising score matching with the noise variance as weight, and the DDPM sampling step is a Langevin-like step with the learned score. Song et al. (2021) make the correspondence exact in continuous time: both models discretize a stochastic differential equation that adds noise, and generation runs its time reversal, whose drift contains the score $\nabla_x \log p_t(x)$.

What changed#

	Likelihood (VAE, flows)	Score matching
What is learned	a density or a bound on it	the gradient of the log-density
Normalizing constant	needed, or avoided by the architecture	never needed
Training signal	log-likelihood or ELBO	Fisher divergence, via a trace, random projections, or denoising
Sampling	one decoder pass	many Langevin steps
Weak spot	the architecture restricts the model	regions without data; separated modes, mitigated by noise levels

Limitations#

• Many sampling steps. Annealed Langevin runs hundreds to thousands of network evaluations per sample; faster samplers are the subject of later work.
• Discretization bias. Unadjusted Langevin with a finite step samples a slightly wrong distribution, as the variance 1.026 instead of 1 shows. A decreasing step removes the bias only if the steps still add up to infinity.
• The noise schedule is a design choice. The largest level must cover the distances between modes, the smallest must be below the detail that matters, and the steps between levels must overlap.
• No likelihood for free. A score model does not give $\log p(x)$ directly; it can be computed through the probability flow ODE of Song et al. (2021), at the cost of solving it.

Summary#

• The score $\nabla_x \log p(x)$ does not depend on the normalizing constant and, where $p$ is positive on a connected region, determines $p$ up to it.
• Hyvärinen’s integration by parts turns the Fisher divergence into $\mathbb E_p[\nabla\cdot s_\theta + \frac12 \lVert s_\theta\rVert^2]$; sliced score matching estimates the trace with random projections; denoising score matching regresses onto the score of a perturbation kernel, $-\varepsilon/\sigma$ for Gaussian noise, and learns the score of the noised data.
• Langevin dynamics samples with the score; its step size trades bias against the number of steps, a decaying step must sum to infinity, and a temperature $\tau$ samples $p^{1/\tau^2}$ instead of $p$.
• One noise level fails in empty regions and between modes; noise-conditional scores with annealing mitigate both when the levels and step sizes suit the data. The network receives the level in every layer and is trained with one random level per example.
• DDPM’s noise prediction is noise-conditional denoising score matching: $s_\theta = -\varepsilon_\theta / \sqrt{1 - \bar\alpha_t}$.

Notebook#

The score matching notebook of the HSE-LAMBDA Deep Generative Models repository works on a two-dimensional Swiss roll and is based on the ACIDS-IRCAM diffusion tutorial (cell 1; cell numbers count every cell, markdown included, from 1). It has no saved outputs, so it documents the procedure rather than results.

• Cell 5: 10,000 points of a Swiss roll with noise 0.5, the first and third coordinates divided by 10.
• Cells 7–12: standard score matching with a two-layer Softplus network, the exact Jacobian built by one backward pass per output coordinate, and 2000 full-batch Adam steps.
• Cells 13–14: plot_gradients draws the score as arrows on a $50 \times 50$ grid over $[-1.5, 2]^2$, each rescaled to $v \log(1 + \lVert v\rVert)/\lVert v\rVert$, the length rule of the score-field figure, so that the huge scores far from the data stay readable. Cells 19, 23, 27 and 32–33 reuse the same plot.
• Cells 15–19: gradient ascent with the score (20 steps of $10^{-3}$) against Langevin dynamics (10 steps from $\epsilon = 10^{-2}$, multiplied by 0.9 after each step), both from the fixed point $(1.5, -1.5)$; the text says the start is drawn from $\mathcal N(0, I)$. Cell 17 credits the convergence of Langevin dynamics to Welling and Teh (2011). sample_langevin (cell 18) has a temperature argument, default 1, that scales the noise. Its geometric decay of the step breaks the condition $\sum_k \eta_k = \infty$: the ten steps add up to $0.065$, so the chain barely moves from its start.
• Cells 20–23: sliced score matching on the first 1000 points. The random directions are normalized to unit length, so $\mathbb E[v v^\top] = I/d$, and both terms are averaged over the $d$ coordinates: the Jacobian term becomes $v^\top \nabla s_\theta\, v / d$, with expectation $\mathrm{tr}(\nabla s_\theta)/d^2$, and the second term is the mean of $(v_i\, s_i)^2/2$ rather than $(v^\top s)^2/2$, with expectation $\lVert s_\theta\rVert^2/(2d^2)$. Both carry the same factor $1/d^2$, so the objective is proportional to Hyvärinen’s and has the same minimizer.
• Cells 24–27: denoising score matching. The markdown of cell 24 gives the general objective for any perturbation kernel before the Gaussian one; the code uses $\sigma = 0.01$ and 5000 steps. At this noise level the targets are of size $1/\sigma = 100$ and points off the roll get almost no training signal, which is the failure that motivates the next cells.
• Cells 28–33: a noise-conditional network with the four levels $1$, $0.215$, $0.046$ and $0.01$ and the $\sigma^2$-weighted loss (anneal_power = 2). ConditionalLinear adds a learned vector per level, an nn.Embedding(4, 128) initialized uniformly on $[-0.1, 0.1]$, to the output of each of the two hidden linear layers, the additive variant described above. Each of the 5000 full-batch steps draws one level per point with torch.randint. Cell 32 plots the learned score field at each level. Annealed Langevin sampling is not implemented, and the last plot assigns a random noise level to every grid point, so its arrows mix four different fields.

Next#

Diffusion and score models work in pixel space: every step of the reverse chain processes the full image, and most of its bits describe details no one perceives. The next note runs the same process in the latent space of an autoencoder and adds the conditioning that turns it into text-to-image generation.

References#

• A. Hyvärinen. Estimation of Non-Normalized Statistical Models by Score Matching. Journal of Machine Learning Research, 6, 2005.
• P. Vincent. A Connection Between Score Matching and Denoising Autoencoders. Neural Computation, 23(7), 2011.
• Y. Song, S. Garg, J. Shi, S. Ermon. Sliced Score Matching: A Scalable Approach to Density and Score Estimation. UAI, 2019.
• Y. Song, S. Ermon. Generative Modeling by Estimating Gradients of the Data Distribution. NeurIPS, 2019. Noise-conditional score networks and annealed Langevin dynamics.
• Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, B. Poole. Score-Based Generative Modeling through Stochastic Differential Equations. ICLR, 2021.
• B. Efron. Tweedie’s Formula and Selection Bias. Journal of the American Statistical Association, 106(496), 2011.
• M. Welling, Y. W. Teh. Bayesian Learning via Stochastic Gradient Langevin Dynamics. ICML, 2011. Eq. 2 (step-size conditions).
• G. O. Roberts, R. L. Tweedie. Exponential Convergence of Langevin Distributions and Their Discrete Approximations. Bernoulli, 2(4), 1996.