Introduction

Today's paper formulates score-based generative modeling through continuous-time SDEs. Along with DDPM, the best-known diffusion approach, it is often called a foundational diffusion paper—though its difficulty hardly feels foundational. Unlike DDPM and DDIM, which rely on stepwise prediction, its sampler solves differential equations numerically. The paper's strength is that it unifies and varies many score-based approaches, enabling a discussion of more effective generative sampling.

The mathematics is much harder than DDPM, and stochastic differential equations feel more complex because the scope of a solution changes substantially with its constraints. Still, understanding this paper seems necessary to think beyond methods such as DDPM and DDIM that are confined to stepwise sampling.

Score based generative process

Suppose a dataset $x \sim \mathcal{X}$ has marginal probability distribution $p(\mathcal{\mathcal{X}})$. To start from an easy-to-sample distribution $z \sim \mathcal{Z}$ and create data $\hat{x} \sim q(\mathcal{\hat{X} \vert \mathcal{Z}})$ resembling real observations, we must find the relationship between the two distributions. Deep generative modeling aims to solve this task with deep learning.

One representative method indirectly aligns an intractable probability distribution with the maximum lower bound on likelihood, the ELBO, as in a Variational Autoencoder.

A GAN (Generative Adversarial Network) improves sampling quality through adversarial training between a decoder that creates fake images from randomly sampled noise $z$ and an encoder that distinguishes real images.

Flow-based modeling instead defines an invertible relationship between distributions as a flow and uses it to connect latent and data distributions.

Unlike these models, a score-based model learns the direction—the gradient—from a given noise sample toward the real data distribution and defines a process that gradually produces data samples from noise.

What is score?

I briefly discussed the score in the DDPM post; to restate it, the score is the gradient of a probability density function. Given a distribution $p(x)$ over the original data's continuous variables, the direction in which a noisy sample $x_t$ at step $t$ should move to become plausible is the gradient of $p(x)$. Connecting this to an energy-based approach lets us define the score as the gradient of log-likelihood.

$$\nabla_x \log p(x)$$

Two generative models built from these ideas arrived with dramatic impact: SMLD (Score Matching with Langevin Dynamics) and DDPM (Denoising Diffusion Probabilistic Model). Once generalized, both reduce to the same differential equation with different coefficients, so they can be viewed as solving the same problem. The derivation later makes this explicit.

Score estimation = Solving SDE

How can the mechanism of diffusion models be connected to stochastic differential equations? Consider the figure above. In DDPM or SMLD, data becomes noise simply by adding very small amounts of diffusion. Conventionally this noise is added over discrete states $(x_0,~x_1,~\cdots,~x_T)$. If instead the distribution's gradual transformation into noise is defined over continuous time, it can be modeled as a solution to an Itô SDE.

This becomes clearer later when it is derived in relation to the DDPM and SMLD processes, so for now we can simply accept the idea.

$$dx = f(x,~t)dt + g(t)dw$$

A stochastic differential equation is one type of differential equation, so the problem is defined through relationships among rates of change $(dx,dt,dw)$. A solution can be understood as the trajectory traced by those rates for $x(t)$. An SDE differs in that it adds a stochastic term $dw$ to the deterministic solution $x(t)$, though we can ignore that detail for the moment.

$dw$ is a diffusion term that depends only on time, independently of $x$'s actual trajectory, while $dt$ defines the direction along that trajectory. The solution $x(t)$ is therefore a function that returns the position $x(t)$ at a given time $t$, represented by the functional curve shown as the solid black line above.

Reversing the perspective, a trajectory depends on its direction. We can therefore define a new problem that predicts $x(t)$ while traveling backward from $x_T$ to $x_0$.

$$dx = \left(f(x,~t) -g^2(t)\mathbf{\nabla_x \log p_t(x)}\right) dt + g(t) d\bar{w}$$

In this expression, the diffusion term $g(t)dw$ is identical to the forward SDE. The score term, highlighted in bold, changes the trajectory associated with the drift term $dt$ from the forward direction to the reverse direction.

The forward SDE and reverse SDE have different starting and ending points but must trace the same trajectory. If the forward SDE—the trajectory's form—is known, it can provide the solution to the reverse SDE. Predicting the score is thus equivalent to solving the reverse SDE, which generalizes generation to any sampler implemented as an SDE solver. The paper introduces two types: a PC (Predictor–Corrector) sampler and a deterministic sampler. As their names suggest, the PC sampler uses a score-based model to improve sample quality, while the deterministic sampler removes all randomness to support fast sampling and latent manipulation.

Training and sampling with SMLD/DDPM

I have not properly defined SMLD: it is an NCSN (Noise Conditional Score Network). NCSN and DDPM are both score-based generative models, but their papers develop the SDE problem through different methods. Let us examine how SMLD and DDPM train their score estimators and perform sampling.

Denoising score matching with Langevin Dynamics(SMLD)

Suppose a perturbation kernel with variance $\sigma^2$ is added to $x$. This kernel follows noise distribution $p_\sigma(\tilde{x} \vert x):=\mathcal{N}(\tilde{x};x,\sigma^2I)$. Consider increasing noise scales $\sigma_{\min}=\sigma_1<\sigma_2<\cdots<\sigma_N=\sigma_{\max}$. The smallest scale $\sigma_{\min}$ is small enough that $p_{\sigma_{\min}}(x)\approx p_\text{data}(x)$, while the largest $\sigma_{\max}$ is large enough that $p_{\sigma_{\max}}(x)\approx\mathcal{N}(x;0,\sigma_{\max}^2I)$.

Because the true distribution is intractable, denoising score matching cannot predict $\nabla_x\log p_\text{data}(x)$ directly. Instead, it predicts the conditional score of a noisy sample $\tilde{x}$ with a very small amount of added noise, $\nabla_{\tilde{x}}\log p_{\sigma_i}(\tilde{x}\vert x)$, which approximates the target. Across many noise levels, from severe corruption down to data nearly identical to the original, score estimator $s_\theta(x,\sigma)$ has the following objective.

$$\theta^\ast = \underset{\theta}{\arg\min} \sum_{i=1}^N \sigma_i^2 \mathbb{E}\_{p_\text{data}(x)} \mathbb{E}\_{p_{\sigma_i}(\tilde{x} \vert x)} \left( \parallel s_\theta(\tilde{x}, \sigma_i) - \nabla_{\tilde{x}} \log p_{\sigma_i}(\tilde{x} \vert x) \parallel_2^2 \right)$$

With enough data and model capacity, a trained score-based model predicts the scores required for sampling. Sampling is somewhat complex, but because it uses annealed Langevin dynamics, it can be understood as a nested loop. The noise level decreases from $N$ to $1$ across each $\sigma_i$, and $M$ Langevin MCMC steps run for every noise kernel—a total of $N\times M$ steps. For a distribution corrupted by noise level $\sigma_i$,

$$x_i^m = x_i^{m-1} + \epsilon_i s_{\theta^\ast}(x_i^{m-1}, \sigma_i) + \sqrt{2\epsilon_i}z_i^m$$

run the process above $M$ times with step size $\epsilon_i>0$, analogous to a gradient-ascent learning rate, and random standard-normal vector $z_i^m\sim\mathcal{N}(0,I)$. After $M$ sampling steps, lower the noise to $\sigma_{i-1}$ and run another $M$ MCMC steps, initializing $x_{i-1}^0=x_i^M$.

$$x_{i-1}^m = x_{i-1}^{m-1} + \epsilon_{i-1} s_{\theta^\ast}(x_{i-1}^{m-1}, \sigma_{i-1}) + \sqrt{2\epsilon_{i-1}}z_{i-1}^m$$

for sigma in sigma_N, ..., sigma_1:
	for step in 1, ..., M:
		noise = N(0, I) # Random noise
		x = x + epsilon * model(x, sigma) + sqrt(2*epsilon)*noise

With a sufficiently large number of samples $M$ and a very small step size $\epsilon$, the generated sample $x_1^M$ follows the original data distribution.

$$p_{\sigma_{\min}}(x) \approx p_{\text{data}}(x)$$

Denoising diffusion probabilistic models(DDPM)

DDPM defines a variance schedule $\beta_i$ whose noise scale increases step by step.

$$0 < \beta_1,~\beta_2,~\cdots,~\beta_N < 1$$

For each training data point (image) $x_0\sim p_\text{data}(x)$, it defines a perturbation kernel under the variance schedule as the following Markov chain.

$$p(x_i \vert x_{i-1}) = \mathcal{N}(x_i;~\sqrt{1-\beta_i}x_{i-1},~\beta_i I)$$

Applying this across the full sequence yields the following marginal.

$$p(x_i \vert x_0) = \mathcal{N}(x_i; \sqrt{\bar{\alpha}_i}x_0,~(1-\bar{\alpha}\_i)I),~\bar{\alpha}\_i := \prod\_{j=1}^i (1-\beta_j)$$

As in the SMLD example, this can also be expressed as a perturbed data distribution, where $p_{\bar{\alpha}_i}(\bar{x}\vert x)$ denotes the marginal obtained above.

$$p_{\bar{\alpha}\_t}(\tilde{x}) = \int p_\text{data}(x)p_{\bar{\alpha}_i}(\tilde{x} \vert x) dx$$

In DDPM, the variable's variance is fixed at $1$ throughout noise addition, so $x_N\sim\mathcal{N}(0,I)$. Although DDPM originally expresses the reverse process differently, its variational Markov chain can be written in score-matching form, like SMLD, as follows.

$$p_\theta (x_{i-1} \vert x_i) = \mathcal{N}(x_{i-1};~\frac{1}{\sqrt{1-\beta_i}}(x_i + \beta_i s_\theta(x_i,~i)), \beta_iI)$$

The evidence lower bound consequently has a different weighted sum. The expression below closely resembles the SMLD objective. Its coefficient $1-\bar{\alpha}_i$ is the variance of the $i$th perturbation kernel $p_{\bar{\alpha}_i}(\tilde{x}\vert x)$. In SMLD, too, the coefficient $\sigma_i^2$ was the perturbation kernel's variance, revealing the common structure of the two equations.

More precisely, for noisy score estimate $\nabla_x\log p(\tilde{x}\vert x)$, it is proportional to $1/\mathbb{E}\left[\parallel\nabla_x\log p(\tilde{x}\vert x)\parallel_2^2\right]$.

$$\theta^\ast = \underset{\theta}{\arg\min} \sum_{i=1}^N (1-\bar{\alpha}\_i) \mathbb{E}\_{p_\text{data}(x)} \mathbb{E}\_{p_{\bar{\alpha}\_i}(\tilde{x} \vert x)} \left( \parallel s_\theta(\tilde{x}, i) - \nabla_{\tilde{x}} \log p_{\bar{\alpha}_i}(\tilde{x} \vert x) \parallel_2^2 \right)$$

Given $\theta^\ast$ that minimizes this loss, sampling proceeds backward through the Markov chain.

$$x_{i-1} = \frac{1}{\sqrt{1-\beta_i}} (x_i + \beta_i s_{\theta^\ast}(x_i, i)) + \sqrt{\beta_i}z_i,~i=N,~N-1,~\cdots,~1$$

The authors call this ancestral sampling, naming the gradual backward sampling through reverse process $p_\theta$. DDPM's original loss has a different structure, shown below, but the reformulation above demonstrates that DDPM is in fact training a score-estimation model.

$$\mathbb{E}\_{x_0,\epsilon} \left( \frac{\beta_t^2}{2\sigma_t^2 \alpha_t(1-\bar{\alpha}\_t)} \parallel \epsilon - \epsilon_\theta(\sqrt{\bar{\alpha}\_t}x_0 + \sqrt{1-\bar{\alpha}\_t}\epsilon, t) \parallel^2\right)$$

The key point is that in both SMLD and DDPM, just as inference—the forward process—can be expressed as an SDE, the reverse process can likewise be expressed as solving an SDE through score matching.

Extending Diffusion Modeling

Existing methods can be summarized as applying perturbation kernels at several noise scales to the data. This paper seeks to generalize score-based generation from predetermined discrete noise scales to every continuous noise scale,

which amounts to solving the SDE as a continuous-variable equation, as illustrated below.

As stated earlier, both the forward and reverse diffusion processes can be expressed as SDEs.

Corrupting a Dataset Into Noise With a Differential Equation

At its most basic, diffusion modeling takes an arbitrary i.i.d. dataset $x(0)\sim p_0$ and uses a diffusion process indexed by a continuous or discrete time variable $t\in[0,T]$ to create $x(T)\sim p_T$, which approaches Gaussian noise. As introduced above, the diffusion process that transforms data distribution $p_0$ into prior distribution $p_T$ can be modeled as the solution to a stochastic differential equation.

$$dx = f(x,~t)dt + g(t) dw$$

$w$ is a standard Wiener process, also called Brownian motion. $f(x,t)$ is the drift coefficient over $x,t$, and $g(t)$ is the diffusion coefficient. Assuming $f(x,t)$ and $g(t)$ are both Lipschitz in state and time—that their derivatives are bounded—the SDE has a unique strong solution under bounded rates of change (reference).

$$\begin{aligned} &\vert \nabla_{x,~t} f(x,~t) \vert \le \epsilon_1 \newline &\vert \nabla_{t} g(t) \vert \le \epsilon_2 \end{aligned}$$

Sampling by Solving the Differential Equation in Reverse

If the prior is known or predefined, we can sample $x(T)\sim p_T$. For example, defining the prior as a multivariate Gaussian lets us obtain $x_T$ by Gaussian sampling. A derivation shows that the reverse of a diffusion SDE is itself an SDE, written as follows (reference).

$$dx = \left(f(x,~t) -g^2(t)\mathbf{\nabla_x \log p_t(x)}\right) dt + g(t) d\bar{w}$$

$\bar{w}$ is merely different notation for the same Wiener process $w$ used in the forward process. Noise generation traces the differential equation's trajectory as $t$ increases in tiny increments; the reverse equation traces it as $t$ decreases in tiny increments. If the score of marginal distribution $p_t(x)$ is known, an arbitrary noise sample $x_T$ can be transformed into $x_0$.

Predicting the Score to Solve the Reverse SDE

Sampling ultimately requires predicting the score. Since the score of the original data distribution's marginal is unavailable, we optimize time-dependent score-matching model $s_\theta(x,t)$ with a loss like those of SMLD and DDPM. SMLD and DDPM compute over a discrete perturbation sequence of length $N$; generalizing it to continuous variable $t$ gives

$$\theta^\ast = \underset{\theta}{\arg\min} \mathbb{E}\_t\left(\lambda(t) \mathbb{E}\_{x(0)} \mathbb{E}\_{x(t) \vert x(0)} \left( \parallel s_\theta(x(t), t) - \nabla_{x(t)} \log p_{0t}(x(t) \vert x(0)) \parallel_2^2 \right)\right)$$

As noted earlier, choose $\lambda(t)\propto1/\mathbb{E}\left[\parallel\nabla_{x(t)}\log p_{0t}(x(t)\vert x(0))\parallel_2^2\right]$. The score-matching method need not be denoising score matching; alternatives such as sliced score matching or finite-difference score matching also work.

Extending VE and VP SDEs

The noise-perturbation schemes of SMLD and DDPM can be viewed as discretized SDEs. With $N$ noise scales, SMLD's Markov-chain forward process is

$$x_i = x_{i-1} + \sqrt{\sigma_i^2 - \sigma_{i-1}^2} z_{i-1},~i = 1,~\cdots,~N$$

As $N$ approaches infinity, the noise scale becomes a continuous-time function $\sigma(t)$, and standard Gaussian variable $z_{i-1}$ likewise becomes $z(t)$. Replacing the Markov chain with continuous stochastic process $x(t)$ gives

$$x(t+\Delta) = x(t) + \sqrt{\sigma(t+\Delta)^2 - \sigma(t)^2}z(t)$$

Assuming $\Delta$ is very small, a first-order Taylor approximation turns this into a differential equation—the ordinary Euler method.

$$dx = \sqrt{\frac{d\sigma^2(t)}{dt}} dw$$

Likewise, DDPM follows the forward process

$$x_i = \sqrt{1-\beta_i} x_{i-1} + \sqrt{\beta_i}z_{i-1},~i=1,~\cdots,~N$$

so as $N$ approaches infinity,

$$x(t+\Delta) = \sqrt{1-\beta(t+\Delta) \cdot \Delta}x(t) + \sqrt{\beta(t+\Delta)\cdot \Delta}z(t)$$

Here $\beta$ is the noise scheduled at each original Markov step, so the infinitesimal variance change is multiplied by $\Delta$. Again applying the Euler method yields the differential equation

$$dx = -\frac{1}{2} \beta(t) x dt + \sqrt{\beta(t)} dw$$

The paper calls SMLD's differential equation a VE (Variance Exploding) SDE because variance continually increases, and DDPM's a VP (Variance Preserving) SDE because variance is preserved.

Because the coefficients have an affine form, the expression can be developed into a differential equation for variance, using the following equation.

$$\frac{d\Sigma_{\text{VP}}(t)}{dt} = \beta (t) (I - \Sigma_{\text{VP}}(t))$$

This derivation is possible because it assumes a Gaussian distribution rather than a nonlinear case. Solving the ODE gives the covariance function of $x(t)$.

$$\Sigma_\text{VP} (t) = I+\exp\left(\int_0^t -\beta(s) ds\right)(\Sigma_\text{VP}(0) - I)$$

This shows that variance $\Sigma_\text{VP}$ at every $t$ is bounded relative to $\Sigma_\text{VP}(0)$. If $\Sigma_\text{VP}(0)=I$, as in DDPM, it remains $I$ for all $t$. The authors' new sub-VP SDE is

$$dx = -\frac{1}{2} \beta(t) x dt + \sqrt{\beta(t) (1 - \exp\left(-2\int_0^t \beta(s) ds\right))}dw$$

Solving the covariance ODE in the same way gives the covariance of the new differential equation.

$$\Sigma_\text{sub-VP}(t) = I+\exp\left(-2\int_0^t \beta(s) ds \right)I + \exp\left(-\int_0^t \beta(s) ds \right) (\Sigma_\text{sub-VP}(0)-2I)$$

The VP SDE can be viewed as an upper bound on the sub-VP SDE. Since a diffusion process may take any form $dx=f(x,~t)dt+g(t)dw$, approaching it as a continuous function rather than beginning with discrete diffusion enabled the authors to propose this SDE. Sub-VP reportedly performs well in some likelihood experiments.

Stochastic Differential Equations Used in the Paper

We have introduced three different SDEs. Their perturbation kernels can be written as Gaussian marginal distributions as follows.

$$p_{0t}(x(t) \vert x(0)) = \begin{cases} \mathcal{N}(x(t); x(0), (\sigma^2(t)-\sigma^2(0))I),& \text{(VE SDE)} \newline \mathcal{N}(x(t);x(0)e^{-\frac{1}{2}\int_0^t \beta(s)ds}, I-Ie^{-\int_0^t \beta(s) ds}),&\text{(VP SDE)} \newline \mathcal{N}(x(t);x(0)e^{-\frac{1}{2}\int_0^t \beta(s)ds}, \left(1-e^{-\int_0^t \beta(s) ds}\right)^2I),&\text{(sub-VP SDE)} \end{cases}$$

The sub-VP derivation looks complicated, but it simply changes DDPM's original noise construction from

$$p(x_i \vert x_0) = \mathcal{N}(x_i; \sqrt{\bar{\alpha}_i}x_0,~(1-\bar{\alpha}_i)I),~\bar{\alpha}\_i := \prod\_{j=1}^i (1-\beta_j)$$

to a form that squares the variance term:

$$p(x_i \vert x_0) = \mathcal{N}(x_i; \sqrt{\bar{\alpha}_i}x_0,~(1-\bar{\alpha}_i)^2I),~\bar{\alpha}\_i := \prod\_{j=1}^i (1-\beta_j)$$

Solving the reverse SDE

We have repeatedly established that if a score model gives $s_\theta(\cdot)$ at every time point, the theory lets us construct the reverse SDE. Numerically solving that reverse SDE to find $x_0$ is sampling.

Solving an SDE numerically amounts to predicting the form—the trajectory—of its solution function. Methods include Euler–Maruyama and stochastic Runge–Kutta, among others, though I do not know their details. The important point is that given a score predictor, any SDE solver can generate samples.

Ancestral sampling in DDPM is a special case and cannot be applied when the SDE changes even slightly, making it difficult to generalize. To relax this limitation, the authors propose the reverse-diffusion sampler below, which reportedly performs slightly better than ordinary ancestral sampling in SMLD and DDPM.

$$x_i = x_{i+1} - f_{i+1}(x_{i+1}) + G_{i+1}G_{i+1}^\top s_{\theta^\ast}(x_{i+1}, i+1) + G_{i+1}z_{i+1}$$

Defining an infinitesimal time change $\Delta$ ever closer to zero should reduce the gap between ancestral sampling and the true SDE solution. The original ancestral sampler seems weaker because it depends on $p_\theta(x_i\vert x_{i+1})$ without accounting for the SDE's original form.

Predictor-corrector samplers

Ordinarily, solving an SDE performs only score prediction. Predictor–Corrector sampling adds an MCMC approach, namely Langevin dynamics. After the numerical SDE predicts the solution at the next time step, correction sampling adjusts that point using the score estimate. The numerical SDE solver acts as the predictor, and the score model acts as the corrector.

For both SMLD and DDPM, combining Predictor and Corrector substantially improves sampling efficiency over using either one alone.

Probability flow and Neural ODE

A score-based model lets us solve the SDE numerically. In fact, every diffusion process has an ODE with the same marginal likelihood $p(x_{0:T})$. If the original diffusion SDE is

$$dx = f(x, t)dt + g(t)dw$$

then the ODE with the same marginal likelihood is

$$dx = \left( f(x, t) - \frac{1}{2}g(t)^2 \nabla_x \log p_\theta(x) \right) dt$$

as shown above. The proof is too complicated for me to understand yet. In any case, define $\tilde{f}(x,t)=f(x,t)-\frac{1}{2}g(t)^2\nabla_x\log p_\theta(x)$. The ODE with Wiener-process term $dw$ removed,

$$dx = \tilde{f}(x, t)dt$$

has the same marginal likelihood as the original SDE, meaning that sampling yields the same result. Rather than forcing an approximation of discrete likelihood as DDPM does below,

$$p_\theta(x_0 \vert x_1) = \prod_{i=1}^D \int_{\delta_{-}(x_0^i)}^{\delta_+(x_0^i)} \mathcal{N}(x; \mu_\theta^i (x_1, 1), \sigma_1^2) dx$$ $$\delta_+(x) = \begin{cases} \infty,&\text{if }x=1 \newline x+\frac{1}{255},&\text{if }x<1 \end{cases},~~\delta_-(x) = \begin{cases} -\infty,&\text{if }x=-1 \newline x-\frac{1}{255},&\text{if }x>-1 \end{cases}$$

we can calculate the actual likelihood.

$$\log p_0(x(t)) = \log p_T(x(T)) + \int_0^T \nabla \cdot \tilde{f}_\theta(x(t),t)dt$$

Because an ODE can ignore the random term $dw$ entirely, it has the advantage of a one-to-one mapping between a specific data point $x(0)$ and its endpoint $x(T)$ along the path. This resembles an invertible flow-based model or neural ODE, making image editing and interpolation easy through latent representations. Each image maps to a unique latent, so the forward SDE can also be viewed as an encoder; the forward process is in fact independent of parameter learning.

Finally, removing stochasticity makes sampling with a neural ODE fast.

Conditional sampling

Another benefit of solving the continuous SDE is that it can generate not only unconditional samples but also samples from conditional probability $p_0(x(0)\vert y)$. This requires class probability $p_t(y\vert x(t))$ at each time step.

$$dx = \left( f(x, t) - g(t)^2(\nabla_x \log p_t (x) + \nabla_x \log p_t(y \vert x)) \right)dt + g(t)d\bar{w}$$

The class-conditional probability is learned by training an encoder under supervised learning on noisy samples, just as in the class-guidance paper Diffusion Beats GANs. The remaining methods appear in the appendix; I did not think they needed a detailed discussion here.