Diffusion Models: The Mathematics of Image Generation

Picture a Central Pollution Control Board monitoring station on the Yamuna, watching for illegal effluent discharge. A sensor at Wazirabad reports a spike in dissolved chromium at 6 a.m.; by 9 a.m. the plume has spread and diluted across three kilometres of river. The inverse problem CPCB actually needs solved is: given the spread-out, noisy concentration reading at 9 a.m., work backward to the concentrated point source at 6 a.m. and its exact location. This is a genuinely hard problem because diffusion destroys information — many different source configurations could have produced the same diluted reading, and the physically correct answer requires knowing not just the current concentration field but its gradient, the direction in which concentration is increasing fastest, at every point in space. That gradient of a density field is the single mathematical object this entire chapter is about. Diffusion generative models — the machinery underneath Stable Diffusion, Imagen, and DALL·E's later versions — solve the exact same class of problem on pixels instead of pollutants: given a field of pure noise, they walk backward along the gradient of a probability density to recover something that looks like it was sampled from the distribution of natural images. This chapter builds that machinery from first principles: the forward noising process as a Markov chain, why the reverse of that chain cannot be computed in closed form, score matching as the technique that sidesteps this, the training objective Ho, Jain, and Abbeel actually optimize, and the continuous-time view that unifies discrete DDPM with faster deterministic samplers like DDIM.

The forward process as a Markov chain

Sohl-Dickstein, Weiss, Maheswaranathan, and Ganguli's 2015 paper "Deep Unsupervised Learning using Nonequilibrium Thermodynamics" is where the name "diffusion model" comes from: they borrowed the mathematics of a substance diffusing through a medium — the same physics governing that chromium plume — and pointed it at data distributions instead of physical space. Start with a clean image (or, for the worked arithmetic below, a single scalar) x₀ drawn from the true data distribution q(x₀). Define a fixed, non-learned Markov chain that corrupts it over T steps by adding a small amount of Gaussian noise at each step:

q(x_t | x_{t-1}) = N( x_t ; sqrt(1 - β_t) · x_{t-1} ,  β_t · I )

Here β_t ∈ (0,1) is a variance schedule, typically increasing with t so that later steps inject more noise. Because each step only depends on the immediately preceding state, this is a first-order Markov chain — the defining property that makes the whole framework tractable is that q(x_t | x_{t-1}, x_{t-2}, …, x_0) = q(x_t | x_{t-1}). Crucially, this forward chain is fixed by design, not learned: there is no network involved in going from x₀ to x_T, only a schedule of numbers.

The chain has a closed form that lets you jump straight from x₀ to any x_t without simulating the intermediate steps, which is what makes DDPM training practical. Define α_t = 1 − β_t and the cumulative product ᾱ_t = ∏_{s=1}^{t} α_s. Because the sum of independent Gaussians is itself Gaussian, repeatedly applying the reparameterization trick collapses the whole chain into one shot:

q(x_t | x_0) = N( x_t ; sqrt(ᾱ_t) · x_0 ,  (1 − ᾱ_t) · I )
x_t = sqrt(ᾱ_t) · x_0 + sqrt(1 − ᾱ_t) · ε,   ε ~ N(0, I)

As t → T, ᾱ_t → 0, so x_T converges to pure isotropic noise regardless of what x_0 was — the chain provably forgets everything about the original image. That forgetting is precisely the "nonequilibrium thermodynamics" framing: entropy strictly increases along the forward direction, exactly as it does in the chromium plume spreading through the river. The generative task is to run time backward and locally decrease entropy — which physics forbids for an isolated system, but a learned model can approximate given enough training data, because it has access to the full population of images, not just one sample's physical trajectory.

Why you cannot simply run the chain backward

Common misconception: a student who has just learned that q(x_t | x_{t-1}) is a simple Gaussian often assumes the reverse step q(x_{t-1} | x_t) must also be some fixed, computable Gaussian — as if you could just algebraically invert the forward formula and run the same machinery backward. This is false, and seeing exactly why is the whole reason score matching exists. Bayes' rule gives the true reverse conditional as

q(x_{t-1} | x_t) = q(x_t | x_{t-1}) · q(x_{t-1}) / q(x_t)

The numerator's first factor is the known forward Gaussian, but q(x_{t-1}) and q(x_t) are marginal densities over the entire population of natural images at that noise level — integrals over the whole data distribution, which is exactly the object we do not have an analytic formula for. (If we did, we would not need a generative model in the first place.) The reverse conditional is only tractable in closed form when it is additionally conditioned on x_0, i.e. q(x_{t-1} | x_t, x_0), and even then only because you know which particular x_0 produced this trajectory — information a sampler starting from noise does not have. Anderson's 1982 result on reverse-time diffusion equations shows that in the continuous-time limit the reverse process is itself a diffusion process of the same mathematical family, but its drift term requires exactly one extra piece of information not present in the forward equations: the gradient of the log-density, ∇x log q_t(x). That gradient is called the score function, and learning it — not "inverting a formula" — is what makes generation possible.

Score matching: learning a gradient you cannot integrate

Score matching originates with Hyvärinen's 2005 paper on estimating non-normalized statistical models. The core difficulty it solves: for almost any interesting distribution, the density p(x) = exp(−E(x)) / Z has an intractable normalizing constant Z (the integral of exp(−E(x)) over all of image space). You cannot maximize log-likelihood directly because you cannot compute Z. But the score, ∇x log p(x) = −∇x E(x), does not depend on Z at all — the constant vanishes under differentiation. So instead of learning the density, you learn its gradient field directly, by minimizing the expected squared error between a network's output and the true score:

J(θ) = E_{x~p(x)} [ ‖ s_θ(x) − ∇x log p(x) ‖² ]

This objective still looks unusable, since it references the unknown ∇x log p(x) on the right. Vincent's 2011 paper made the connection that turns this into something trainable: if you add Gaussian noise to a clean sample to get x̃ = x₀ + σε, the score of the resulting noised conditional distribution has a closed form, because it is just a Gaussian centered at x₀:

q(x̃ | x_0) = N(x̃ ; x_0, σ²I)
∇x̃ log q(x̃ | x_0) = −(x̃ − x_0) / σ² = −ε / σ

Vincent's denoising score-matching theorem proves that minimizing the network's error against this tractable conditional score, averaged over noise draws, minimizes the same objective as matching the true marginal score of the noised data distribution, up to a constant that does not affect the optimum. In other words: you never need to know the true score of natural images. You only need to know the score of a Gaussian you added yourself, and a network trained to predict it will, in expectation, learn the real thing. This is the identical trick CPCB's inverse problem exploits in reverse: because the forward diffusion physics is known exactly (Fick's law, known diffusion coefficient), you can compute what any hypothetical source would look like after diffusing, and search for the source whose forward simulation best matches the observed plume — you never need to "invert" the physics symbolically, only run it forward and compare.

DDPM's training objective is score matching in disguise

Ho, Jain, and Abbeel's 2020 "Denoising Diffusion Probabilistic Models" paper parameterizes the reverse process as a chain of learned Gaussians, p_θ(x_{t-1} | x_t) = N(x_{t-1}; μ_θ(x_t, t), σ_t²I), and derives a variational lower bound on log p_θ(x_0) by writing the joint over the whole reverse trajectory and matching it, step by step, to the tractable forward posterior q(x_{t-1} | x_t, x_0). That posterior — the one term Bayes' rule made tractable above — turns out to be Gaussian with a known mean:

μ̃_t(x_t, x_0) = [ sqrt(ᾱ_{t-1})·β_t / (1 − ᾱ_t) ] · x_0 + [ sqrt(α_t)·(1 − ᾱ_{t-1}) / (1 − ᾱ_t) ] · x_t

The variational bound reduces to a sum of KL divergences between this known Gaussian and the model's predicted Gaussian at each step — and since matching two Gaussians of equal variance is just matching their means, the loss becomes a weighted squared error between μ̃_t and μ_θ. The key design choice is how to parameterize μ_θ. Rather than predicting x_0 or μ directly, Ho et al. substitute the forward relation x_0 = (x_t − sqrt(1−ᾱ_t)·ε) / sqrt(ᾱ_t) into μ̃_t and have the network predict the noise ε instead:

μ_θ(x_t, t) = (1/sqrt(α_t)) · [ x_t − (β_t / sqrt(1 − ᾱ_t)) · ε_θ(x_t, t) ]

Substituting this back, every per-step KL term collapses to a constant times ‖ε − ε_θ(x_t, t)‖². Ho et al. then found empirically that dropping the theoretical weighting constant entirely and training on the unweighted sum over all t — the now-famous simplified objective — produces better samples than the exact variational weighting:

L_simple(θ) = E_{t, x_0, ε} [ ‖ ε − ε_θ(x_t, t) ‖² ],   x_t = sqrt(ᾱ_t)x_0 + sqrt(1−ᾱ_t)ε

Compare this to the denoising score-matching target derived above: ∇x̃ log q(x̃|x_0) = −ε/σ. With σ = sqrt(1 − ᾱ_t), this says the true score at noise level t is −ε / sqrt(1 − ᾱ_t). A network trained to predict ε is therefore, up to the known scale factor −1/sqrt(1 − ᾱ_t), a network trained to predict the score. ε_θ and a score network s_θ are the same object under a linear reparameterization:

ε_θ(x_t, t) = −sqrt(1 − ᾱ_t) · s_θ(x_t, t)

This is not a loose analogy — it is an exact algebraic identity, and it is why the U-Net inside Stable Diffusion, trained with the plain "predict the noise" objective, can be plugged directly into score-based samplers (annealed Langevin dynamics, the probability-flow ODE) with no retraining. The two research lineages — Song and Ermon's 2019 score-based models trained with denoising score matching, and Ho et al.'s noise-prediction DDPM — converge on the identical network.

Worked example: tracing the chain and verifying the identity by hand

Take a toy one-dimensional "image" x₀ = 2.0 and a four-step schedule β = [0.1, 0.2, 0.3, 0.4], so α = [0.9, 0.8, 0.7, 0.6]. Fix a single noise draw ε = 0.5 for the whole trajectory, since the closed-form jump formula only needs one. First compute the cumulative products:

t   β_t   α_t   ᾱ_t     sqrt(ᾱ_t)   1−ᾱ_t   sqrt(1−ᾱ_t)
1   0.1   0.9   0.9000   0.9487      0.1000   0.3162
2   0.2   0.8   0.7200   0.8485      0.2800   0.5292
3   0.3   0.7   0.5040   0.7099      0.4960   0.7043
4   0.4   0.6   0.3024   0.5499      0.6976   0.8352

Now generate each noised sample directly from x_0 using x_t = sqrt(ᾱ_t)·x_0 + sqrt(1−ᾱ_t)·ε:

x_1 = 0.9487(2.0) + 0.3162(0.5) = 2.05548
x_2 = 0.8485(2.0) + 0.5292(0.5) = 1.96163
x_3 = 0.7099(2.0) + 0.7043(0.5) = 1.77200
x_4 = 0.5499(2.0) + 0.8352(0.5) = 1.51743

Notice the mean is pulling steadily toward zero (the schedule's implicit prior) while x_t is not monotonically decreasing in a simple way — at t=1 the noise term happens to push it slightly above x_0 before the shrinking mean dominates from t=2 onward. Now compute the score two independent ways and confirm they agree, which is the numerical content of the ε–score identity above. Way one, from the noise-prediction formula score_t = −ε / sqrt(1−ᾱ_t). Way two, by taking the actual gradient of the log-density log N(x_t; sqrt(ᾱ_t)x_0, 1−ᾱ_t) = −(x_t − sqrt(ᾱ_t)x_0)²/(2(1−ᾱ_t)) + const, whose derivative with respect to x_t is −(x_t − sqrt(ᾱ_t)x_0)/(1−ᾱ_t):

t   score (via ε)   score (via ∇ log N)
1   −1.581139       −1.581139
2   −0.944911       −0.944911
3   −0.709952       −0.709952
4   −0.598641       −0.598641

They match to full floating-point precision at every step, because they are algebraically the same quantity. Notice the score's magnitude shrinks monotonically as t increases — at t=1, x_t is still close to the sharp data point, so the gradient pointing back toward it is steep; by t=4, the sample has drifted far enough into the noise regime that the density is nearly flat and the corrective gradient is weaker in relative terms even though the absolute noise is larger. This is a direct arithmetic illustration of the signal-to-noise ratio, SNR_t = ᾱ_t / (1−ᾱ_t), which for this schedule runs 9.00, 2.57, 1.02, 0.43 across the four steps — dropping through 1.0 (the point where signal and noise contribute equally) between steps 3 and 4, which is exactly the region where a denoiser's job is hardest and most informative for training.

From a discrete chain to a continuous stochastic process

Song, Sohl-Dickstein, Kingma, Kumar, Ermon, and Poole's 2021 paper "Score-Based Generative Modeling through Stochastic Differential Equations" showed that as T → ∞ with steps shrinking accordingly, the DDPM forward chain becomes the discretization of a continuous-time stochastic differential equation, the variance-preserving SDE:

dx = −½β(t)·x dt + sqrt(β(t)) dw

where w is a standard Wiener process. Anderson's 1982 reverse-time result then gives an SDE that runs this process backward, provided the score is known at every noise level:

dx = [ −½β(t)·x − β(t)·∇x log p_t(x) ] dt + sqrt(β(t)) dw̄

with dw̄ a reverse-time Wiener process. Simulating this reverse SDE with a learned score network reproduces ancestral DDPM sampling in the continuous limit. But the same paper shows something more useful for deployment: there exists a companion ordinary differential equation, the probability-flow ODE, with no stochastic term at all, whose solution trajectories pass through exactly the same marginal densities p_t(x) at every time t as the stochastic version:

dx = [ −½β(t)·x − ½β(t)·∇x log p_t(x) ] dt

Because it is deterministic, this ODE can be solved with far larger step sizes than the noisy SDE tolerates, using standard numerical solvers. Song, Meng, and Ermon's 2021 "Denoising Diffusion Implicit Models" derives — independently, from a non-Markovian generalization of the forward process rather than from the SDE — a sampler with the identical deterministic update:

x̂_0 = ( x_t − sqrt(1−ᾱ_t)·ε_θ(x_t,t) ) / sqrt(ᾱ_t)
x_{t−1} = sqrt(ᾱ_{t−1}) · x̂_0 + sqrt(1−ᾱ_{t−1}) · ε_θ(x_t,t)

which the SDE paper shows is precisely an Euler discretization of the probability-flow ODE for the variance-preserving case. Applying this to the worked example: at t=4, x_4 = 1.51743 and, if the network's prediction were perfect (ε_θ = ε = 0.5), the recovered x̂_0 = (1.51743 − 0.8352·0.5)/0.5499 = 2.0000 — exactly the original data point, confirming the algebra is self-consistent. This is the mathematical reason DDIM sampling in Stable Diffusion can produce a full image in 20–50 steps instead of DDPM's original 1000: both are approximating the same underlying continuous trajectory, but the deterministic ODE view licenses a much coarser step size for a comparable accumulated error.

Visualizing the mechanism

Active recall

Attempt each question before reading its answer.

1. Using a two-step schedule β₁ = 0.1, β₂ = 0.2, x₀ = 3, and a fixed noise draw ε = −0.4, compute x₂ and the score ∇x log q(x₂|x₀) by hand.
2. Why does the forward process q(x_t | x_{t-1}) being a simple, known Gaussian not imply that the reverse conditional q(x_{t-1} | x_t) is also directly computable?
3. True or false, with justification: the noise-prediction network ε_θ used in DDPM and a score network s_θ used in score-based generative models are trained with unrelated objectives and learn unrelated functions.
4. Using the worked example's t=4 values (ᾱ₄ = 0.3024, x₄ = 1.51743) and assuming a perfectly trained ε_θ(x₄,4) = 0.5, compute the DDIM-predicted x̂₀.
5. Suppose you replace the linear schedule β_t (small at t=1, growing to β_T) with a cosine schedule in the style of Nichol and Dhariwal's 2021 "Improved Denoising Diffusion Probabilistic Models," which keeps ᾱ_t closer to 1 for longer before dropping off more gradually near the end. Trace the full ripple: what happens to (a) the ᾱ_t curve shape, (b) the SNR schedule ᾱ_t/(1−ᾱ_t), (c) which noise levels dominate training under the unweighted L_simple objective, and (d) the number of DDIM steps typically needed for comparable sample quality?

Answers.

1. α₁ = 0.9, α₂ = 0.8, so ᾱ₂ = 0.72, sqrt(ᾱ₂) = 0.84853, sqrt(1−ᾱ₂) = 0.52915. Then x₂ = 0.84853(3) + 0.52915(−0.4) = 2.54559 − 0.21166 = 2.33393. The score is −ε/sqrt(1−ᾱ₂) = −(−0.4)/0.52915 = 0.75593 — positive, because the negative noise draw pushed x₂ below the scaled mean 0.84853 × 3 = 2.54559, so the gradient of the log-density points back up toward that mean.

2. By Bayes' rule, q(x_{t-1}|x_t) ∝ q(x_t|x_{t-1})·q(x_{t-1}). The first factor is the known forward Gaussian, but q(x_{t-1}) is the marginal density of the entire data distribution at noise level t-1 — an integral over every possible image, which has no closed form. Knowing one conditional direction of a Markov step says nothing about the reverse conditional unless you also know the marginal densities involved, and those marginals are exactly the object a generative model exists to approximate. The reverse step only becomes tractable when additionally conditioned on the specific x₀ that generated the trajectory, which a sampler starting from noise doesn't have — it has to be supplied by a learned score/noise estimate instead.

3. False. They are related by the exact identity ε_θ(x_t,t) = −sqrt(1−ᾱ_t)·s_θ(x_t,t), derived by comparing DDPM's simplified loss to Vincent's denoising score-matching target ∇x̃ log q(x̃|x_0) = −ε/σ with σ = sqrt(1−ᾱ_t). A noise-prediction network is a score network up to a known, timestep-dependent linear rescaling — which is why the identical U-Net weights trained with the plain MSE-on-noise objective can be dropped into score-based samplers like the probability-flow ODE without retraining.

4. x̂₀ = (x₄ − sqrt(1−ᾱ₄)·ε_θ) / sqrt(ᾱ₄) = (1.51743 − 0.83522 × 0.5) / 0.54991 = (1.51743 − 0.41761) / 0.54991 = 1.09982 / 0.54991 = 2.00000. This recovers the original x₀ = 2.0 exactly, confirming that with a perfect noise estimate the DDIM update is self-consistent at any timestep, not just t=1.

5. (a) A cosine schedule keeps ᾱ_t near 1 for a larger fraction of early t and only drops sharply nearer t=T, instead of the linear schedule's roughly steady decay — visually, the curve stays flat longer then falls off, rather than descending at a constant rate. (b) The SNR curve ᾱ_t/(1−ᾱ_t) inherits this shape: under the linear schedule SNR crashes very early (most information is destroyed in the first few percent of steps, exactly as seen in the worked example where SNR fell from 9.0 to 0.43 in just four steps) and then changes very little near t=T, wasting many steps at a noise level that's already saturated; the cosine schedule spreads the SNR drop more evenly across the whole range of t. (c) Because L_simple averages ‖ε−ε_θ‖² uniformly over t with no explicit weighting, whichever schedule is used determines which actual noise levels those uniformly-sampled t indices land on — under the linear schedule, a large share of training compute is spent on t values that are already deep in the near-zero-SNR regime and contribute little useful gradient signal, whereas the cosine schedule allocates training more evenly across the perceptually important mid-SNR region where denoising is genuinely ambiguous and hardest to get right. (d) Because the SNR trajectory is smoother and has no near-vertical drop, the probability-flow ODE it corresponds to is less stiff, so a DDIM-style solver can take larger Euler steps for the same discretization error — in practice this is part of why cosine-scheduled models need noticeably fewer sampling steps than early linear-schedule DDPM to reach comparable sample fidelity.

Think About It

Think about this: How would you explain diffusion models: the mathematics of image generation to a friend who has never seen a computer? What real-world analogy would you use? Imagine you had to build a system using these concepts — what would be your first step? Try this: before moving on, write down three things you learned and one question you still have.

Key Takeaways — Summary and Recap

Let us recap what we covered: the core ideas behind diffusion models: the mathematics of image generation, how they connect to real-world applications, and why they matter for your journey in computer science. Remember these key points as you move forward. For competitive exam preparation (CBSE, JEE, BITSAT), focus on understanding the WHY behind each concept, not just the WHAT.