Score-Based Diffusion Models: Denoising and Generative Modeling via Score Functions

The Question Behind Every AI-Generated Image

Open Stable Diffusion or DALL-E 2 and type "a tiger sitting on the steps of the Gateway of India." A few seconds later, a coherent image appears — not retrieved from a database, but built up from a canvas of pure random static. Strip away the marketing language and one exact, answerable mathematical question sits underneath: if you are handed a grid of random numbers (static) and told "this should eventually look like a photograph," which direction should each pixel move, right now, to look a little more like a real photograph and a little less like noise?

That question has a precise mathematical answer, and this chapter derives it from first principles. The object that answers it is called the score function. Once you understand what it is, how to compute it, how to estimate it when you cannot compute it directly, and how to use it to walk a random point from pure noise into a realistic sample, you understand the mathematical core of every major diffusion-based image, audio, and video generator in use today.

The Score Function: A Compass on the Probability Landscape

Let p(x) be a probability density function describing "how likely is this point x" — for now, imagine x is a single real number, and later we will let it be a vector of a million pixel values. High p(x) means x looks like real data; low p(x) means it looks like noise.

Throughout this chapter, "log" means the natural logarithm, as is standard in machine learning. Define the score function as

s(x) = d/dx log p(x)

— the derivative of the log-density with respect to x (a gradient, once x becomes a vector). Notice two things immediately. First, this is a vector (or, in one dimension, a signed number), not a probability — it tells you a direction and a strength, never a "how likely." Second, because log is an increasing function, wherever p(x) increases, log p(x) increases too, and its derivative points the same way p(x) is rising. So s(x) always points toward higher density — it is a compass needle that tells a lost point which way "more realistic" lies.

If you have studied electrostatics for JEE Physics, you have already used this exact mathematical object. The electric field is E = −∇V, the negative gradient of the electric potential; it points from high potential to low potential and its magnitude tells you the strength of the pull. The score is the same construction applied to log-probability instead of potential, and (with the sign flipped) it points from low density to high density. Keep this analogy — it will make the sampling algorithm in a few sections feel like a familiar physics problem rather than a new one.

One warning about vocabulary: in classical statistics, "score function" sometimes means the gradient of the log-likelihood with respect to a parameter θ (used in maximum-likelihood estimation and the Cramér–Rao bound). The score function in this chapter is the gradient with respect to the data point x, a completely different object introduced for generative modeling by Aapo Hyvärinen in 2005. Both are real, standard terms — just don't let a search engine mix them up for you.

Worked Example: The Score of a Two-Mode Distribution

Take a concrete data distribution: a 50/50 mixture of two Gaussians, one centered at μ₁ = −2 and one at μ₂ = 3, each with variance σ² = 0.25. This is a reasonable toy stand-in for "real data lives in a few clusters, with empty space between them" — exactly the shape of a real high-dimensional image distribution, just squashed into one dimension so we can compute everything by hand.

For a single Gaussian N(μ, σ²), log p(x) = −(x − μ)² / (2σ²) + constant, so its score is a one-line derivative:

s(x) = d/dx log p(x) = −(x − μ) / σ²

This is exactly Hooke's law in disguise: a linear restoring force pulling x back toward μ, growing without bound the farther you stray. For our mixture p(x) = 0.5·p₁(x) + 0.5·p₂(x), the score is not simply the average of the two component scores — you have to differentiate log(p₁ + p₂), which by the chain rule gives a weighted average where the weights are the posterior responsibilities:

s(x) = γ₁(x)·s₁(x) + γ₂(x)·s₂(x), where γₖ(x) = pₖ(x) / (p₁(x) + p₂(x))

γₖ(x) is "how much of the density at this exact point comes from component k" — it shifts smoothly from 1 to 0 as x moves from one mode's territory into the other's. This responsibility-weighted average is the general rule for the score of any mixture, and it is the same weighting rule used in the E-step of the EM algorithm for Gaussian Mixture Models.

Plugging in numbers (s₁(x) = −4(x + 2), s₂(x) = 4(3 − x) since 1/σ² = 4) gives exact values you can check by hand:

• x = −5 (far left of both modes): γ₁ ≈ 1, so s(x) ≈ s₁(−5) = −4(−3) = +12 — a strong pull rightward toward the near mode.
• x = −2 (exactly at mode 1): s(x) = 0. The peak is a stationary point of the density.
• x = −1 (just past mode 1): s(x) = −4(1) = −4 — pulled back toward the peak, like a restoring spring.
• x = 0.5, the exact midpoint between −2 and 3: by symmetry, γ₁ = γ₂ = 0.5 and s₁(0.5) = −10, s₂(0.5) = +10, so they cancel exactly: s(x) = 0.
• x = 3 (mode 2): s(x) = 0 again.
• x = 6 (far right): s(x) ≈ s₂(6) = 4(3 − 6) = −12.

That midpoint result deserves a second look: the score is zero at both peaks and at the valley between them. A single Gaussian's score is an unbounded, ever-strengthening restoring force; a mixture's score is far richer, because the responsibility weights γₖ(x) can swing the balance to exactly zero in low-density territory too. This single fact — that "score equals zero" does not mean "you are at the data" — is the most commonly misunderstood point in the whole subject, and the figure below marks it in red for exactly that reason.

Why the Score, and Not the Density Itself?

A natural question: why not just estimate p(x) directly and use that? The answer is a specific algebraic escape hatch, and it is the entire reason score-based methods exist. Any probability density that comes from a learned scoring function f(x) — as in classical energy-based models — must be written as

p(x) = f(x) / Z, where Z = ∫ f(x) dx

Z is the normalizing constant that forces the total probability to equal 1. For any interesting high-dimensional model (a million-pixel image space), that integral is completely intractable — there is no way to compute it, and every step of ordinary maximum-likelihood training needs its exact value. This is the single biggest obstacle in classical energy-based modeling.

Now take the score of that same p(x):

s(x) = d/dx log p(x) = d/dx [log f(x) − log Z] = d/dx log f(x) − d/dx log Z

Because Z does not depend on x at all, log Z is just a constant, and the derivative of a constant is zero. The intractable term vanishes:

s(x) = d/dx log f(x)

The score only ever needs f(x), never Z. This is why Hyvärinen's original 2005 score matching framework, and every diffusion model built on it since, can train on unnormalized models without ever touching the partition function. You give up knowing "exactly how likely is x" and, in exchange, get "exactly which way is x should move" — for a generative model that only needs to produce samples, that is a trade you are happy to make.

Sampling With the Score: Langevin Dynamics

Suppose you have s(x) for every x. How do you actually turn that into random samples from p(x), rather than just walking uphill to the single highest peak (which would collapse every sample onto one mode)? You need to climb the gradient while adding just enough randomness to keep exploring. The tool for this is Langevin dynamics, named after the same physical process (Paul Langevin's 1908 description of Brownian motion) that governs a pollen grain jostled by water molecules while a concentration gradient nudges it toward denser regions. The discrete update rule, at step size ε, is:

xₜ₊₁ = xₜ + (ε/2)·s(xₜ) + √ε·zₜ, where zₜ ~ N(0, 1) is fresh random noise each step

The first term after xₜ is the "climb toward higher density" drift, using exactly the score you computed. The second term is pure noise, which is what stops every sample from marching straight to a peak and stopping there — it keeps a cloud of samples spread out with a density that (as ε → 0 and the number of steps → ∞) provably converges to p(x) itself. We use this fact — that this particular drift-plus-noise recipe has p(x) as its long-run stationary distribution — without deriving the full proof, which requires the Fokker–Planck equation from stochastic calculus; the physical intuition (uphill pull balanced against random spreading, exactly like heat diffusing against a temperature gradient) is enough to trust and use the rule correctly.

Worked Example: Two Steps of Langevin Dynamics by Hand

Using the same mixture as before (μ₁ = −2, μ₂ = 3, σ² = 0.25) and ε = 0.1, start at x₀ = 0 and trace two steps by hand, using illustrative random draws z₁ = 0.3 and z₂ = −0.5 (chosen here for a clean worked example, not from an actual random number generator).

Step 1. At x₀ = 0, compute the score: since x₀ = 0 is 4 standard deviations from μ₁ but far more standard deviations from μ₂, the responsibility γ₁(0) ≈ 0.99995 ≈ 1, so s(0) ≈ s₁(0) = −4(0 + 2) = −8. Update:

x₁ = 0 + 0.5(0.1)(−8) + √0.1 · (0.3) = 0 − 0.4 + (0.316)(0.3) = −0.4 + 0.095 = −0.305

Step 2. At x₁ = −0.305, again γ₁ ≈ 1 (p₁ at this point is about 0.0032, while p₂ is on the order of 10⁻¹⁰), so s(x₁) ≈ s₁(−0.305) = −4(−0.305 + 2) = −4(1.695) = −6.78. Update:

x₂ = −0.305 + 0.5(0.1)(−6.78) + √0.1 · (−0.5) = −0.305 − 0.339 − 0.158 = −0.802

Two steps, and despite the random pushes in both directions, x has already drifted from 0 to −0.802, heading toward the mode at −2 exactly because the score consistently pulls that way. Run this for a thousand steps with a small ε, and the point spends most of its time hovering near −2 or 3 in proportion to how much probability mass each mode holds — which is precisely what "sampling from p(x)" means.

import numpy as np

def sample_data(n):
    choice = np.random.rand(n) < 0.5
    means = np.where(choice, -2.0, 3.0)
    return means + 0.5 * np.random.randn(n)

def gaussian_score(x, mu, var):
    return -(x - mu) / var

def true_score(x):
    var = 0.25
    p1 = np.exp(-(x + 2.0)**2 / (2*var))
    p2 = np.exp(-(x - 3.0)**2 / (2*var))
    s1 = gaussian_score(x, -2.0, var)
    s2 = gaussian_score(x, 3.0, var)
    w1 = p1 / (p1 + p2)
    w2 = p2 / (p1 + p2)
    return w1 * s1 + w2 * s2

def langevin_sample(steps=1000, eps=0.01, x0=0.0):
    x = x0
    for t in range(steps):
        noise = np.random.randn()
        x = x + 0.5 * eps * true_score(x) + np.sqrt(eps) * noise
    return x

Trace it: true_score(-5) computes p1 ≈ e⁻¹⁸, p2 ≈ 0, so w1 ≈ 1 and it returns s1 = −4(−5+2) = 12 — matching the +12 arrow in the figure exactly. true_score(0.5) gives p1 = p2 = e⁻¹²·⁵ (by symmetry), so w1 = w2 = 0.5, s1 = −10, s2 = +10, and the function returns 0.5(−10) + 0.5(10) = 0 — the "valley" zero from the figure, reproduced exactly by running the code.

Diffusion: Corrupting Data on Purpose

Everything so far assumed we already had s(x) for the real data distribution. We do not — we only have a finite set of training images, not a formula for p(x). Diffusion models solve this with a clever trick: deliberately corrupt the data with known, controllable noise, because that noising process has a score we can write down exactly, and — as the next section shows — that gives us a target we can train a neural network to hit.

Define a sequence of increasingly noisy versions of a real data point x₀, using a fixed, non-learned schedule of noise levels σ₁ < σ₂ < ... < σ_L:

xₜ = x₀ + σₜ · z, where z ~ N(0, I) is fresh standard Gaussian noise

At small σ, xₜ is barely disturbed. At large σ, xₜ is almost entirely noise, and its distribution approaches N(0, σ²I) regardless of what x₀ was — the data's structure gets washed out completely. This is the forward process, and it is deliberately simple: it needs no neural network, no training, nothing learned — you can sample from it in one line, as shown in the figure's top-row arrow. Yang Song and Stefano Ermon (Stanford, NeurIPS 2019) formalized this as Noise Conditional Score Networks; Jonathan Ho, Ajay Jain, and Pieter Abbeel (UC Berkeley, NeurIPS 2020) formalized a closely related variant, variance-preserving diffusion, in the paper that coined the now-standard name Denoising Diffusion Probabilistic Models (DDPM).

Why One Noise Level Is Not Enough: The Manifold Problem

Here is the obstacle that makes multiple noise levels necessary rather than a mere implementation detail. Real data — natural images, in particular — does not fill up the entire space of possible pixel values; it clusters on a thin, low-dimensional manifold inside a vastly larger space (a random grid of pixel noise looks nothing like a photograph, and the "photograph-like" region is a tiny sliver of everything a pixel grid could be). Far away from that manifold, p(x) is close to zero everywhere, log p(x) plunges toward negative infinity, and any score estimate learned from finite real samples in that region is either or wildly unreliable — there is essentially no training data anywhere near a random starting point.

If you tried Langevin dynamics using only the score of the clean, unnoised data distribution, a sample starting at random noise would have almost no useful gradient signal to follow and would wander for an impractically long time, if it found the manifold at all. Song and Ermon's fix is exactly the multi-scale noising you just saw: adding noise at a large σ smears the tight data manifold out until it fills the whole space, making the score well-defined and learnable everywhere. Sampling then starts at a large σ, where the smoothed density is easy to climb, and gradually anneals σ down toward zero, tightening the target step by step until the sample lands on the real data manifold. This procedure — called annealed Langevin dynamics — is the actual sampling algorithm inside modern diffusion models, and the bottom panel of the figure is its cartoon: reading right to left, a noise-filled patch is nudged, at progressively finer noise scales, back onto the shape of the original image.

Learning the Score You Cannot Compute: Denoising Score Matching

We still need s(xₜ), the score of the noised marginal distribution, and we do not have a formula for it — we only know the forward corruption rule. The key result, due to Pascal Vincent (2011), is that we can get an exact, tractable training target out of the one thing we do control: the noise we ourselves added.

For a single clean point x₀, the conditional distribution of the noised version is Gaussian, p(x | x₀) = N(x; x₀, σ²), so exactly as in the earlier worked example:

∇ₓ log p(x | x₀) = −(x − x₀) / σ² = −z / σ, using x − x₀ = σz

That is a formula for the score of one noised point given its clean origin. The result we actually need is the score of the full noised marginal p(x) (averaged over all possible x₀), and Vincent's identity connects the two using Bayes' rule and the log-derivative identity ∇f = f∇(log f):

∇ₓ p(x) = ∫ p(x₀) ∇ₓ p(x|x₀) dx₀ = ∫ p(x₀) p(x|x₀) ∇ₓ log p(x|x₀) dx₀

Dividing both sides by p(x), and noting that p(x₀)p(x|x₀)/p(x) = p(x₀|x) is exactly the posterior probability of "which clean image produced this noisy one" by Bayes' rule, this becomes an expectation:

∇ₓ log p(x) = E[∇ₓ log p(x|x₀) | x] = −E[z | x] / σ

In words: the true score of the noised data equals minus the average noise that was added, given only the noisy result — divided by σ. This is Tweedie's formula, and it converts an unreachable target (the marginal score) into a supervised learning problem with a target we know exactly during training: the noise z we drew ourselves. Train a network εθ(x, σ) to predict that noise from the noisy input, using ordinary mean-squared error, and the trained network gives you the score for free: s(x) ≈ −εθ(x, σ) / σ.

import numpy as np

def training_loss(eps_theta, x0, sigma):
    z = np.random.randn(*x0.shape)
    x_noisy = x0 + sigma * z
    z_pred = eps_theta(x_noisy, sigma)
    return np.mean((z_pred - z) ** 2)

def train_step(eps_theta, optimizer, data_batch, sigma_min=0.01, sigma_max=1.0):
    sigma = np.exp(np.random.uniform(np.log(sigma_min), np.log(sigma_max)))
    loss = training_loss(eps_theta, data_batch, sigma)
    optimizer.step(loss)
    return loss

This is simplified pseudocode — eps_theta stands in for a real neural network with learnable weights and optimizer.step(loss) stands in for a real backpropagation call — but the loss itself, mean((z_pred - z)**2), is the exact training objective used in real diffusion models: sample a noise level, corrupt a real example by that much, ask the network to guess the exact noise you added, and penalize the squared error. Nothing about "generating images" appears anywhere in this loss; the network only ever learns to denoise.

From Trained Network to Generated Image

Put the two pieces together and you have the complete recipe used, with engineering variations, by Stable Diffusion (Rombach et al., CompVis/Stability AI/Runway, 2022, operating in a compressed latent space rather than raw pixels), DALL-E 2 (OpenAI, 2022), and Imagen (Google Research, 2022): start from pure noise x_L ~ N(0, σ_L²I); at each noise level from large σ down to small σ, run several steps of Langevin dynamics using the score −εθ(x, σ)/σ recovered from the trained noise-predicting network; anneal σ down to (near) zero; output the final x₀. Text-to-image conditioning adds a text embedding as an extra input to εθ, so the same denoising machinery is steered toward images consistent with the prompt — the score-and-Langevin core described in this chapter is unchanged.

Common Misconceptions, Corrected

• "Score = 0 means you have found real data." False, and this chapter's worked example proves it directly: the valley at x = 0.5 between the two modes also has score exactly zero, yet it is the least likely region on the whole curve. Score zero only means you are at a stationary point of log p(x) along that direction — it could be a peak (stable, real data) or a saddle/valley (unstable, not real data). Distinguishing the two requires checking the second derivative (curvature), not just the score.
• "Diffusion models generate the image in one evaluation of the network." False. A single forward pass only predicts noise at one specific noise level; a full generation requires many Langevin steps across an annealed schedule of decreasing σ, as the manifold-problem discussion showed is mathematically necessary, not just a convenient design choice.
• "The score function tells you the probability of x." False — it is a gradient (a direction and a magnitude), not a probability value; you can have identical score magnitude at wildly different densities, and score alone never tells you which.

Where This Fits in Your Exams

The calculus here — differentiating log p(x), the chain rule used to get the mixture score, and completing the square inside a Gaussian exponent — is squarely within the differentiation syllabus tested in CBSE Class 12 and JEE Main/Advanced. The E = −∇V analogy connects directly to electrostatics in JEE Physics, and recognizing that structure is a genuinely transferable exam skill: whenever a problem gives you a scalar potential-like quantity and asks for the "force" or "pull," differentiating it is the method, whether the potential is electrical, gravitational, or (as here) a log-probability. The full multivariable gradient and the Bayes'-rule argument behind Tweedie's formula go beyond the JEE syllabus into first-year engineering mathematics and machine learning theory covered at the GATE and Olympiad-preparation level — exactly the kind of research-adjacent material that rewards a strong Class 11–12 calculus foundation.

Check Your Understanding

• Q1. For a single Gaussian N(μ=5, σ²=2), compute the score at x=7. Answer: s(7) = −(7−5)/2 = −1.
• Q2. A two-component mixture has equal weights, shared variance 0.25, and means −1 and 4. Besides the two means, find another x where the score is guaranteed to be exactly zero. Answer: the symmetric midpoint x = 1.5, by the same cancellation argument as the x = 0.5 valley worked out in this chapter.
• Q3. Why can score-based training skip the normalizing constant Z, while ordinary maximum-likelihood training of the same model cannot? Answer: the score is a derivative with respect to x, and log Z does not depend on x, so it differentiates to zero; maximum likelihood instead needs the actual numeric value of log p(x), which includes log Z and cannot drop it.
• Q4. In the forward process x = x₀ + σz with σ = 2, if a training example has x₀ = 6 and the noised version is x = 10, find the noise z that was added, and the exact conditional score ∇ₓ log p(x|x₀) at that point. Answer: z = (x−x₀)/σ = 4/2 = 2. Score = −z/σ = −2/2 = −1, which matches the direct formula −(x−x₀)/σ² = −4/4 = −1.
• Q5. True or false, with correction: "A point where the score is zero is always a point of high data density." Answer: False — the valley worked example in this chapter is a counterexample: score is zero there too, despite being the lowest-density region shown.
• Q6. Why do diffusion models need many noise levels σ₁ > σ₂ > ... > σ_L rather than training the score once on clean data? Answer: real data occupies a thin manifold in a much larger space; far from that manifold the clean-data score is or unreliable, so a single random starting point has no useful gradient to follow. Noising at a large σ smears the manifold to fill the space, making the score learnable everywhere, and annealing σ down during sampling guides the point from that smoothed landscape onto the true manifold step by step.

Summary

The score function s(x) = ∇ log p(x) is a vector field that always points toward higher probability density, and — crucially — it can be computed from an unnormalized model without ever touching the intractable normalizing constant Z, because differentiating with respect to x kills the constant log Z term. Langevin dynamics turns this vector field into a sampler by alternating a drift along the score with injected random noise, converging to samples from p(x) itself. Diffusion models supply a score we can actually estimate by deliberately corrupting data with known Gaussian noise at many scales; Vincent's denoising-score-matching identity shows that training a network to predict the added noise, via a simple mean-squared-error loss, recovers the true score exactly, via s(x) ≈ −εθ(x,σ)/σ. Generation then runs this whole pipeline backward: start from pure noise, and use annealed Langevin dynamics — many small score-guided steps across decreasing noise levels — to walk that noise into a sample that looks like real data. That entire chain, from a one-line derivative to a trained noise-predicting network to annealed sampling, is the complete mathematical machinery inside Stable Diffusion, DALL-E 2, and Imagen.

Practice Exercises

Now it is time to practice! Complete these challenges to solidify your understanding:

• Exercise 1: Write a short program that demonstrates the core concept from this chapter. Test it with at least 3 different inputs.
• Exercise 2: Find a real-world example where score-based diffusion models: denoising and generative modeling via score functions is used in an Indian company (like TCS, Infosys, Flipkart, or ISRO). Write a paragraph explaining the connection.
• Exercise 3: Create a mind-map connecting score-based diffusion models: denoising and generative modeling via score functions to at least 3 other topics you have studied.