Diffusion probabilistic models (Sohl-Dickstein 2015): Deep Unsupervised Learning using Nonequilibrium Thermodynamics
Add noise to data in many small steps until nothing is left, then train a network to undo the steps one at a time.
Each small noising step has a Gaussian reverse, so the network only has to predict a mean and a variance per step. Run the learned reverse steps from pure noise and you get samples; average the log of the ratio of reverse to forward step probabilities along a noising path and you get a lower bound on the log-likelihood.
Explaining the paperDeep Unsupervised Learning using Nonequilibrium ThermodynamicsA drop of ink spreads through a glass of water until the glass is evenly gray. The spreading is easy to simulate and very hard to run backward. This 2015 paper shows that if the spreading happens in small enough steps, each backward step is simple enough for a neural network to learn.
This is the paper that introduced diffusion models. DDPM, score-based SDEs, Stable Diffusion and the diffusion policies used in robotics all keep its three parts: a fixed process that turns data into noise, a learned process that turns noise back into data, and an objective that ties the two together. Several details that are now standard (predicting the noise, the squared-error loss, the formula for jumping to any noise level in one step) came five years later, and this page marks which is which.
The argument has five pieces, taken in order below: why flexible probability models are usually intractable; the forward process that dissolves data into a Gaussian; why the reverse of a small noising step is itself Gaussian, so a network can learn it by regression; the bound that trains the network and also reports a likelihood; and how the same chain can be conditioned on an observation to inpaint or denoise an image.
Flexible or tractable
A generative model is a probability distribution fitted to data, say 32×32 color images. You want to do three things with it: draw samples, evaluate for a given image, and fit it to a dataset. The paper calls a family tractable when those three operations are cheap and exact, and flexible when it can match the structure of real data.
A Gaussian is tractable: its density is a formula, sampling takes one call to a random number generator, and fitting it means computing a mean and a covariance. It cannot represent a dataset of handwritten digits, because digits are not an ellipsoidal cloud. The usual route to flexibility is an energy model, any non-negative function divided by its total mass:
For an image, is an integral over a 3072-dimensional space, and nobody can compute it for an expressive . Without you cannot evaluate , and fitting or sampling falls back on Markov chain Monte Carlo, which is slow and only approximate.
The paper gets both properties by defining the model as the end of a long Markov chain. Each link in the chain is a Gaussian, so each link can be sampled and evaluated in closed form. The flexibility comes from the number of links: a chain of 1000 slightly different Gaussian steps can carry a standard Gaussian onto a very non-Gaussian distribution, and the abstract describes the result as a deep generative model "with thousands of layers or time steps".
The forward process: how diffusion destroys data
The chain is built from two processes. The forward process is fixed and needs no training: it takes a data point and applies noising steps. Each step scales the current point down slightly and adds a little Gaussian noise:
The number is the variance of the noise added at step , and it is small. The mean is the old point multiplied by , which is slightly less than 1. Shrinking the signal while adding noise keeps the variance fixed: if has variance 1 per coordinate, then has variance . The paper scales every dataset except the binary one to variance 1, and the chain ends at the standard Gaussian instead of growing without bound. Score-SDE (2021) later named this kind of chain "variance preserving"; the 2015 paper gives the kernel without a name.
With a constant , one step multiplies the signal by . After 100 steps the original signal has been multiplied by , after 500 steps by , and after 1000 steps by , at which point the sample is Gaussian noise with a trace of the image. The paper uses steps for its image models, except the bark textures (), and 40 steps for its two-dimensional swiss roll.
The schedule does not have to be constant. For Gaussian diffusion the paper learns by gradient ascent on the training objective, with fixed to a small constant to prevent overfitting; a footnote adds that a fixed schedule seems to work as well. For binary data it uses the fixed schedule , which removes the same 1/T fraction of the original signal at every step, so after of steps a fraction of it is left: half at step 500.
The figure runs the forward process on the paper's toy dataset. Drag the step from 0 to 1000 and watch the readout: the signal coefficient falls and the noise coefficient rises while the total variance stays at 1.00.
The figure jumps straight to step . Composing Gaussian steps gives one Gaussian whose signal coefficient is the product , written since DDPM. The paper's text never writes this marginal or the notation; it gives only the single step (1). The reference implementation released with it does use the cumulative product to sample in one jump during training, as the code later on this page shows.
Why small steps make the reverse step Gaussian
To generate, the model has to run the chain backward: start from a draw of and undo one step at a time. Undoing step means sampling from , the distribution over where a point was one step earlier given where it is now. That distribution depends on the data, so it is unknown and has to be learned.
The paper rests on a classical result it attributes to Feller (1949): for a diffusion with small enough steps, the reverse step has the same functional form as the forward step. The forward step is Gaussian, so for small the reverse step is also Gaussian, and a Gaussian is described completely by a mean and a covariance. Learning the reverse process therefore reduces to learning, for every step, a function that outputs a mean and a function that outputs a covariance. In the paper's words, "the task of estimating a probability distribution has been reduced to the task of performing regression on the functions which set the mean and covariance of a sequence of Gaussians."
The step has to be small because of Bayes' rule. The reverse step is proportional to the forward likelihood times the distribution of the earlier point:
The first factor is a Gaussian bump of width about around . When is small, the bump is narrow, and over such a short range the log of the second factor is close to a linear function. A Gaussian times the exponential of a linear function is another Gaussian, shifted, so the product is still a single Gaussian bump. When is large, the bump is wide enough to cover several clusters of the data, and the product inherits their peaks.
The figure shows this for one-dimensional data with two clusters. The teal curve is the exact reverse distribution over the clean value given one noisy observation; the dashed curve is the closest single Gaussian; the readout is the KL divergence between them. At the reverse is one narrow bump and the KL is 0.01 nats. At the default it has split into two peaks, the best Gaussian sits in the empty valley between them, and the KL is 0.48.
This sets the trade-off behind the method. Small steps make every reverse step learnable by regression but require a long chain; large steps shorten the chain but make each reverse step multimodal, and a Gaussian reverse kernel then fits it badly. Diffusion takes the long chain. The model's reverse step is
where and are outputs of one network that reads the noisy point and the step index. The network outputs the reverse Gaussian's mean and its covariance (diagonal, for images) directly. It does not predict the added noise ; that parameterization, with its simpler loss, is DDPM's.
Sampling: run the learned chain backward
Once and are trained, sampling is a loop. Draw from . At each step, evaluate the network to get a mean and a covariance, and draw from that Gaussian. After steps, is a sample from the model. A fresh Gaussian draw is made at every step, which makes this an ancestral sampler: each step samples from the conditional distribution given the step before, the way you sample a Bayesian network from its roots down.
The figure runs this loop on the swiss roll with 40 steps, the length the paper used for this dataset. The reverse steps here are the exact ones for this toy distribution, standing in for a trained network. The teal arrows are the drift , the change in mean each step applies to a point, which the paper plots in the bottom row of its Figure 1. Press Play and watch the arrows and the readout. At every arrow points at the center of the data and the longest one is 4.1 units, almost the full distance from the corner of the grid to the center. By the arrows over the roll have shrunk, the outer ones still point inward, and the particles have spread around the spiral; by most of them sit on its arms.
The drift changes because of how much signal is left. At the noisy point says almost nothing about which part of the roll it came from, so the best estimate of the earlier point is close to the mean of the whole dataset and the step mostly replaces the input with fresh noise. As the noise level falls, a point carries enough information to place it near one arm, and the mean of the step moves it toward that arm instead of toward the center. The coarse layout of a sample is settled at high noise and the fine detail at low noise.
Switch the fresh noise off and every particle follows the mean of each step. The first reverse step from forgets almost everything about its input, so all 110 particles start the second step at nearly the same point, follow the same means, and end together on one spot. The spread of the samples comes from the noise drawn with variance at every step.
How diffusion models are trained: the bound K
Training needs an objective. The natural one is the log-likelihood of the data, but the model's probability of a clean point is an integral over every path the reverse chain could have taken to reach it:
For images that integral runs over 1000 × 3072 dimensions. The paper borrows a trick from annealed importance sampling: multiply and divide inside the integral by the forward path probability , which is known exactly. The likelihood becomes an average, over forward paths started at the data point, of the ratio between reverse and forward step probabilities. Forward paths are cheap to sample. The log of an average is at least the average of the log (Jensen's inequality), which gives a lower bound:
is a lower bound on the log-likelihood, so training maximizes ; the loss handed to the optimizer is . The paper notes that the derivation parallels the bound used in variational Bayesian methods, the same shape of argument that trains a variational autoencoder.
Appendix B of the paper applies Bayes' rule to the forward process and rewrites (3) as a sum of KL divergences plus entropy terms, all of which can be computed in closed form:
The sum starts at and carries a minus sign, so maximizing means shrinking every KL term. The three entropies at the end depend only on the forward schedule and the prior, not on the network. Inside each KL, the network's reverse step is compared with , the reverse of the forward process when the clean starting point is also known. That conditioning makes the target computable: given , the pair , is jointly Gaussian, so the target is a Gaussian with a closed-form mean and variance. The KL between two Gaussians is also closed form, so every term in (4) is analytic, and one training step can score a single randomly chosen .
The figure computes exactly for a one-dimensional dataset of two clusters, scaled to variance 1, where the true reverse steps can be written down. Each bar is one step's KL between the true reverse step and the model's Gaussian, split in two. The violet part is the error of an imperfect kernel, which training removes. The amber part is what remains when the kernel is the best possible Gaussian, the cost of approximating a non-Gaussian reverse step by a Gaussian. The bottom bar measures against the baseline the paper reports in its Table 1, , the log-likelihood of the plain Gaussian . For this dataset no model can beat that baseline by more than 0.740 bits.
Read three settings off the figure. With , falls short of its ceiling, plus 0.740 bits, by the full 0.740 bits at every training level: a single Gaussian step from noise to data can at best produce a Gaussian, and the best Gaussian for variance-1 data is itself. With full training the shortfall is 0.060 bits at , 0.006 at , and under 0.001 by . With training at 0 the shortfall is 0.740 bits for every , because that untrained chain reproduces exactly.
One step of training on a real example looks like this. Take a CIFAR-10 image as a vector of numbers, add uniform noise to undo the 8-bit pixel quantization, and scale to variance 1. Pick of . Under the reference code's schedule about half the original variance is left at step 500, so
The network reads and the step index and returns a 3072-dimensional mean and a 3072-dimensional diagonal variance for . The closed-form target gives another mean and variance. The per-pixel Gaussian KL between the two, averaged over pixels and images and multiplied by to stand for the sum over steps, is the estimate of that gets backpropagated:
# One training step, CIFAR-10, after model.py in the reference code.
# x0: (n, 3072) images with uniform dequantization noise, scaled to variance 1.
t = randint(1, T) # ONE step per minibatch, T = 1000
abar = prod(1 - beta[:t + 1]) # signal variance left at step t
eps = randn_like(x0) # frozen noise: held fixed for the gradient
xt = sqrt(abar) * x0 + sqrt(1 - abar) * eps # jump straight to step t
mu_q, var_q = forward_posterior(x0, xt, t) # q(x_{t-1} | x_t, x0), closed form
z_mu, z_var = net(xt, t) # multi-scale conv net, (n, 3072) each
mu_p = xt * sqrt(1 - beta[t]) + z_mu * sqrt(beta[t])
var_p = sigmoid(z_var + logit(beta[t])) # reverse variance near beta[t]
kl = gauss_kl(mu_q, var_q, mu_p, var_p) # (n, 3072), analytic per pixel
neg_K = T * kl.mean() + entropy_terms # one-sample estimate of -K
neg_K.backward(); opt.step() # descend -K = ascend the bound KThe noise is drawn once and held fixed while the gradient is taken, which the paper calls frozen noise and credits to Kingma and Welling. Holding fixed makes a differentiable function of the schedule, and that is how the paper's Gaussian runs get a gradient for . In the reference code the schedule is fixed, so the gradient reaches only the network. Sampling needs no data and no target; it runs the learned chain from noise:
# Sampling: ancestral rollout of the learned reverse chain (sampler.py).
x = randn(n, 3072) # start from the prior pi = N(0, I)
for t in range(T - 1, 0, -1): # walk the noise level down
mu, sigma = net_mu_sigma(x, t) # learned reverse mean and std
x = mu + sigma * randn(n, 3072) # fresh noise at every step
return x # one sample per rowThe network and the reference code
The paper used three networks. The swiss roll model is a radial basis function network with one hidden layer of 16 units, shared across all 40 steps, with a separate readout layer per step; its variance output passes through a sigmoid to stay between 0 and 1. The binary sequence model is a multilayer perceptron with three hidden layers of 50 units that outputs the Bernoulli rates for 2000 steps of binomial diffusion. The image models use a multi-scale convolutional network: each layer mean-pools the image by powers of two, convolves at every scale, upsamples, and sums, which gives a full-resolution output that still sees long-range structure.
The image network does not output the mean and variance directly. For each pixel it outputs two vectors of coefficients that weight a set of bump functions of , producing per-pixel values and that vary smoothly across steps. Appendix D then turns them into the reverse step as a perturbation of the forward kernel:
The function is the logistic sigmoid. With the reverse variance equals the forward , so the network starts from a sensible guess and learns a correction. The mean moves the pixel a fraction of the way toward : a pixel at with and gets mean . The paper describes as the mean of the distribution you would reach by applying the same reverse step many times.
The public reference implementation, which the paper says produced the CIFAR-10 results, differs from this in a few places. Its active mean is , with a form close to (5) left commented out. It fixes at a schedule that removes a constant fraction of the original variance per step, with learning it marked as a to-do. It samples one step per minibatch and multiplies that step's KL by , as in the code above (file and line numbers are in the Provenance panel at the end).
Where the thermodynamics comes in
The title refers to nonequilibrium statistical physics. The forward process is a diffusion in the physical sense: it moves a system from a structured state to the featureless equilibrium . The paper borrows two results about such processes.
The first is the likelihood estimator. The average of a reverse-to-forward path ratio in (3) is the weight used by annealed importance sampling (Neal, 2001), a method for estimating an intractable normalizing constant by moving gradually from a simple distribution to a hard one. The Jarzynski equality (1997) is the physics version of the same identity: averaging the exponentiated work over many fast, irreversible runs of a process gives the equilibrium free-energy difference exactly. The paper cites both as the source of its approach. This page adds an analogy the paper does not make: Jensen's inequality in (3) plays the part of the second law of thermodynamics, which says the average work done on a system is at least its free-energy change.
The second is when the bound becomes exact. A process run slowly enough to stay in equilibrium at every moment is called quasi-static, and it is reversible: the forward and reverse path distributions are identical. If they are identical, the ratio inside (3) is the same for every path, so a single sampled path gives the exact likelihood and the inequality becomes an equality. For the model, that limit is infinitely many infinitely small steps with a perfect reverse kernel, which is what Figure 4 approaches as grows at full training. Small steps alone are not enough; the reverse kernels also have to match. The paper also uses the physics to frame the schedule: in thermodynamics the schedule taken between equilibrium states determines how much free energy is lost, and in the model the schedule affects how tight can get; the paper learns it by gradient ascent on .
Inpainting and denoising: multiplying distributions
Many uses of a generative model are conditional: fill in a missing region of an image, or remove noise from a corrupted one. Both mean multiplying the model by a second function that encodes the observation, and renormalizing:
For inpainting, is a delta function that pins the known pixels to their observed values and is constant over the missing ones. For denoising, is the Gaussian likelihood of the noisy image given a candidate original. The renormalization needs a new normalizing constant, the same kind of integral that makes energy models intractable, and the paper notes that multiplying distributions is costly for variational autoencoders, GSNs, NADEs and most graphical models.
The figure shows what (6) produces for the swiss roll when observes only the horizontal coordinate. A vertical line crosses the spiral several times, so the conditioned distribution has several separate arcs.
Under a diffusion model the product can be built into the sampler. The paper multiplies every intermediate distribution of the reverse chain by , using the same at every step in its experiments, which changes each reverse step to
The chain also starts from instead of ; in the inpainting figure of the paper, that start is the image with its missing region replaced by Gaussian noise. Each reverse step is a narrow Gaussian. If is smooth over the width of that Gaussian, a first-order expansion of around the mean (Appendix C) shifts the mean and leaves the covariance alone:
The shift has a plus sign, it is the covariance times the gradient of , and the gradient is evaluated at the unshifted mean , not at . Steering a diffusion sampler with the gradient of a log-likelihood at every step is the mechanism Dhariwal and Nichol named classifier guidance in 2021, with a classifier's probability of the target class and a scale on the gradient; classifier-free guidance (2022) removed the separate classifier.
The approximation (7) holds only when is smooth relative to . For a Gaussian observation of width , the shift moves the mean a fraction of the way to the observed value, so the remaining offset is multiplied by each step. When that factor has magnitude above 1: the step lands further from the observed value than it started, on the other side. If enough steps do this, the offsets compound and the chain diverges. The exact product of the two Gaussians moves the mean a fraction , which is always below 1. When can be multiplied with a Gaussian in closed form, the paper multiplies exactly instead, and the reference sampler does so for both demonstrations: for denoising it combines the two Gaussians into one with precision equal to the sum of their precisions, and for inpainting it sets each known pixel to its value with zero variance at every step.
The figure runs the conditioned 40-step chain both ways, with the same noise draws. At the default width , exceeds 2 on 22 of the 40 steps, the product of over the chain is about , and the first-order samples are thrown off the plot. The exact product puts about 90% of the same samples within of the band's center, on the arcs that Figure 5 lights up. At the first-order version also works: the first five, noisiest steps have up to 5.9 and overshoot, but the remaining steps have and pull the samples back. Push the slider to with the exact product and some samples stay in the band but sit between the two arcs, off the roll. Multiplying each step by and renormalizing step by step is itself an approximation: the paper replaces the exact modified step, whose normalizer depends on , with a normalized one, and notes the two agree when is a small perturbation of each step.
The paper's inpainting result uses the exact case. A diffusion model trained on bark textures () fills in the central 100×100-pixel region of a bark image by sampling the missing pixels from the posterior given the rest, and a crack entering from the left continues into the filled region. The denoising demonstration corrupts CIFAR-10 test images with Gaussian noise of variance 1 (a signal-to-noise ratio of 1) and samples clean images from the posterior, with the Gaussian noise likelihood.
Results: what the 2015 numbers mean
The paper reports the bound on held-out data for every dataset, and also , the improvement over an isotropic Gaussian (or independent binomial) model of the same data. Higher is better for both.
- Swiss roll: bits, 6.45 bits above the Gaussian baseline.
- Binary heartbeat (length-20 sequences with a 1 in every fifth bin): bits per sequence. The true distribution has , so the model is within 0.092 bits of perfect.
- Dead leaves (synthetic images of overlapping discs that mimic natural-image statistics): 1.489 bits per pixel, against 1.244 for MCGSM, the previous best, so the paper claims state of the art. Since 1.489 is a lower bound, the model's true log-likelihood is at least that.
- CIFAR-10: bits per pixel, 11.5 above the Gaussian baseline, on data dequantized with uniform noise. An earlier arXiv version reported higher numbers that came from the model learning the 8-bit quantization of the pixels; this is the corrected figure.
- Bark: bits per pixel, 1.5 above the baseline.
MNIST is reported differently. Most earlier MNIST likelihoods were Parzen-window estimates, which fit a kernel density to model samples and score test digits under it, so the paper does the same with the code from Goodfellow et al. (2014). The diffusion model scores 317 ± 2.7 bits, against 325 ± 2.9 for the adversarial net, 309 ± 1.6 for a deep GSN, 199 ± 2.9 for a DBN, 174 ± 2.3 for a stacked contractive autoencoder, and 349 ± 3.3 for the training data itself. The paper calls its result comparable to these methods, and it trails the GAN.
No single row is the main result. The paper's claim is that one method, trained the same way, fits a 2-D toy distribution, binary sequences, digits, natural images and textures, and for all of them can draw samples, report a likelihood bound, and compute posteriors such as an inpainting.
What later papers changed
DDPM (Ho, Jain and Abbeel, 2020) kept the forward process and the bound but changed the parameterization: the network predicts the noise that was added, the reverse variance is fixed instead of learned, and the weighted sum of KL terms is replaced by an unweighted squared error on . That version produced sharp images. A year earlier, NCSN (Song and Ermon, 2019) learned the score at several noise levels and sampled with Langevin dynamics, and Score-SDE (2021) showed that both are discretizations of one stochastic differential equation. The continuous-time version of Feller's reverse-process result that Score-SDE builds on is Anderson (1982), which this paper does not cite.
Later work moved further from the 2015 chain: DDIM sampled in tens of steps instead of a thousand, latent diffusion ran the chain in an autoencoder's latent space, and flow matching replaced the chain with a learned velocity field. The guidance used in text-to-image models descends from the multiplying-distributions step in (7).
Questions you might still have
Is this the same thing as DDPM?
Same model family, different training target. DDPM (2020) predicts the added noise instead of the reverse mean, fixes the reverse variance instead of learning it, and replaces the full bound with a plain squared-error loss. This 2015 paper learns the reverse Gaussian's mean and variance and maximizes the bound K. The DDPM explainer covers the later version.
Why thousands of small steps instead of a few big ones?
The reverse of a noising step is close to Gaussian only when the step is small. A big step from a two-cluster distribution has a two-peaked reverse that no single Gaussian fits, and the bound pays for the mismatch. The paper used 1000 steps for most image models, 500 for bark, and 40 for the 2-D swiss roll.
Where is the score function?
Not in this paper. It samples by ancestral rollout of the learned reverse chain and never computes the gradient of log p. Score matching with Langevin sampling arrived with NCSN (2019), and Score-SDE (2021) connected the two lines.
Did the paper learn the noise schedule?
For Gaussian diffusion, yes: beta_2 to beta_T are trained by gradient ascent on K, with beta_1 fixed to a small constant. A footnote adds that a fixed schedule seems to work as well, and the released reference code keeps the schedule fixed. Binomial diffusion uses the fixed schedule beta_t = 1/(T - t + 1).
Why predict a variance at all?
The reverse step is a Gaussian, and the bound scores both its mean and its variance. The reference code parameterizes the variance as a learned perturbation of the forward beta_t. DDPM later found that fixing the variance and dropping the bound gave sharper samples, and later work went back to learning it.
How does it inpaint?
It multiplies the model by a second distribution r at every reverse step. For inpainting, r pins the known pixels, so the sampler sets those pixels to their observed values with zero variance at each step and samples only the missing region. For a smooth r the step mean shifts by the variance times the gradient of log r, which is the mechanism later called classifier guidance.
Is the reported likelihood exact?
No. Every number except the MNIST row is the lower bound K, so the true log-likelihood is at least that high. The bound is exact only when the forward and reverse trajectory distributions are identical, the quasi-static limit of infinitely many small steps. The MNIST row is a Parzen-window estimate from samples.
Footnotes & further reading
- The paper: Sohl-Dickstein, Weiss, Maheswaranathan, Ganguli, Deep Unsupervised Learning using Nonequilibrium Thermodynamics (ICML 2015, PMLR 37). Reference implementation: Diffusion-Probabilistic-Models (
model.py,sampler.py). - The reverse-of-a-diffusion result: W. Feller, On the theory of stochastic processes, with particular reference to applications (Berkeley Symposium, 1949). The continuous-time reverse process, a reverse-time SDE with a score-modified drift: B. D. O. Anderson, Reverse-time diffusion equation models (1982), not cited by the paper.
- Ho, Jain, Abbeel, Denoising Diffusion Probabilistic Models (2020): noise prediction, the notation, and the simplified loss. See the DDPM explainer.
- Song & Ermon, Generative Modeling by Estimating Gradients of the Data Distribution (NCSN, 2019), and Song et al., Score-Based Generative Modeling through Stochastic Differential Equations (2021), the source of the names "variance preserving" and "variance exploding".
- Kingma & Welling, Auto-Encoding Variational Bayes (arXiv 2013, ICLR 2014), cited for frozen noise. See the VAE explainer.
- Neal, Annealed Importance Sampling (2001), and C. Jarzynski, Nonequilibrium equality for free energy differences (Physical Review Letters, 1997).
- Dhariwal & Nichol, Diffusion Models Beat GANs on Image Synthesis (classifier guidance, 2021), and Ho & Salimans, classifier-free guidance (2022). See the classifier guidance and classifier-free guidance explainers.
How could this explainer be improved? Found an error, or something unclear? I read every message.