Why Image Generation Takes Fifty Steps

A diffusion model runs the same network dozens of times to make one picture. This post asks why it cannot run it once — and answers with an experiment where the network is provably perfect and the sampler still fails completely.
generative
diffusion
Published

August 1, 2026

Generate an image with a diffusion model and watch what the GPU does. It runs the same network, with the same weights, twenty to fifty times. Each run nudges the current picture. Only after the last one do you get an image.

That should feel wasteful. One network, one prompt — why fifty passes for one answer?

The tempting explanation is that the network isn’t good enough, and a better one would need fewer passes. That explanation is incomplete. Below we build a case where the model is provably exact and a single-pass sampler still fails completely. The experiment isolates a second error source that is easy to miss: a correct local prediction can be used in a bad long-range update.

TL;DR — The network answers a question about the state you have right now. The sampler must turn a sequence of those local answers into a complete trip from noise to image. At pure noise the correct answer is the dataset mean, so one giant step lands there even with a perfect model. Re-evaluating as the state changes lets the answer become mode-specific, then instance-specific. Sampling cost is the price of tracking a field that changes underneath you — and the error falls as \(1/n\) in the number of steps, which we measure.

Diffusion in one paragraph

Two processes, opposite directions.

Forward: take a clean sample and progressively mix in Gaussian noise until nothing of the original remains. No network, no learning — a fixed recipe you can run in one line.

Reverse: start from pure noise and undo that corruption to land on something that looks like data. This is the hard direction, and it is what a diffusion model is trained to enable.

The asymmetry is the whole subject. The forward process needs no model and no steps. The reverse needs a network, called repeatedly. This post is about why that repetition is structural rather than an implementation detail.

What the network actually computes

Diffusion and flow models use several noise parameterisations. To isolate the sampling problem we use the linear path

\[z_t = (1-t)\,x + t\,\epsilon, \qquad t \in [0,1] \tag{1}\]

with \(x\) a clean sample and \(\epsilon \sim \mathcal{N}(0, I)\). At \(t=0\) the state is data; at \(t=1\) it is pure noise. Models may predict \(\epsilon\), \(x\), a score, or a velocity; away from the endpoints these are algebraically interchangeable given the schedule. We use the clean-sample form and call it a denoiser.

Now the subtlety that turns out to be the whole story. Many clean samples could have produced the same noisy \(z_t\) — that is what noise does. The network is trained under squared error, and the value minimising squared error against a set of possibilities is their average. So the network does not pick one plausible clean sample. Its optimum is the conditional average:

\[D(z_t, t) = \mathbb{E}[x \mid z_t] \tag{2}\]

The output is an average, and how much it averages over depends on the noise level. At high noise many clean samples are consistent with \(z_t\), so the average is broad; at low noise few are, so it is sharp.

One word to fix now, because several later arguments depend on it. The prediction is local: valid at this state, at this noise level, and nowhere else. A sampler that treats it as valid over a long stretch of the journey is making an assumption the network never made.

A toy where the model cannot be wrong

On a real model a bad one-step image has two suspects — the network guessed wrong, or the sampler misused a correct guess — and no way to separate them.

So we use data simple enough that the exact denoiser is a formula. Any failure is then unambiguously the sampler’s. Our dataset is two Gaussian blobs at \((\pm 2, 0)\).

import torch, math
import matplotlib.pyplot as plt
torch.manual_seed(0)

A = 2.0                                   # distance of each blob from the origin
S = 0.35                                  # blob standard deviation
MU = torch.tensor([[-A, 0.0], [A, 0.0]])  # the two centres

def sample_data(n):
    which = torch.randint(0, 2, (n,))
    return MU[which] + S * torch.randn(n, 2)

Two blobs is the smallest dataset with a real generative choice: left mode or right mode. That is the one property of real images we need — a photo is of a dog or a cat, not a blend, and the central difficulty in generation is committing to one. It is deliberately not a stand-in for anything else about images. Two dimensions rather than one, because a curve in 1-D is a straight line by definition and we are going to measure curvature.

Deriving the exact denoiser

We want \(\mathbb{E}[x \mid z_t]\) in closed form. Two facts make that possible.

First, the forward path keeps Gaussians Gaussian. If \(x \sim \mathcal{N}(\mu_k, S^2 I)\) then \(z_t = (1-t)x + t\epsilon\) is Gaussian with

\[\text{mean } (1-t)\mu_k, \qquad \text{variance } \ \sigma_t^2 = (1-t)^2 S^2 + t^2 \tag{3}\]

The centre shrinks toward the origin and the spread grows. So the noisy distribution is a mixture of two known Gaussians — no approximation anywhere.

Second, the question splits in two: which blob did this come from? and given that blob, where in it?

For the first, Bayes’ rule over two equally likely blobs gives a posterior depending only on squared distance to each shifted centre:

\[w_k(z_t) = \operatorname{softmax}_k\!\left(-\frac{\lVert z_t - (1-t)\mu_k \rVert^2}{2\sigma_t^2}\right) \tag{4}\]

For the second, a Gaussian prior with a Gaussian likelihood has a Gaussian posterior whose mean is a precision-weighted blend of the prior centre and the rescaled observation — each pulled in proportion to how certain it is:

\[m_k(z_t) = \frac{\mu_k / S^2 \;+\; (1-t)\,z_t / t^2}{1/S^2 \;+\; (1-t)^2/t^2} \tag{5}\]

At small \(t\) the observation term dominates and \(m_k\) tracks \(z_t\); at large \(t\) the prior dominates and \(m_k\) collapses onto \(\mu_k\). Combining:

\[\mathbb{E}[x \mid z_t] = \sum_k w_k(z_t)\, m_k(z_t) \tag{6}\]

def denoiser(z, t):
    """Exact E[x | z_t] for two-blob data — Eq. (6). No network, no training."""
    t = torch.as_tensor(t, dtype=z.dtype)

    var_t = (1 - t)**2 * S**2 + t**2                                 # Eq. (3)
    sq = (z[:, None, :] - (1 - t) * MU).square().sum(-1)
    w = torch.softmax(-0.5 * sq / var_t, dim=1)                      # Eq. (4)

    precision = 1 / S**2 + (1 - t)**2 / t**2
    m = (MU / S**2 + (1 - t) * z[:, None, :] / t**2) / precision     # Eq. (5)

    return (w[..., None] * m).sum(1)                                 # Eq. (6)

Drawn as a map — each arrow runs from a point in space to the clean sample the denoiser predicts for it:

Code 
fig, axes = plt.subplots(1, 3, figsize=(11, 3.2))
gx, gy = torch.meshgrid(torch.linspace(-4, 4, 13), torch.linspace(-2, 2, 7), indexing="ij")
grid = torch.stack([gx.reshape(-1), gy.reshape(-1)], 1)

for ax, t in zip(axes, [0.99, 0.85, 0.40]):
    d = denoiser(grid, torch.tensor(t))
    ax.quiver(grid[:, 0], grid[:, 1], (d - grid)[:, 0], (d - grid)[:, 1],
              angles="xy", scale_units="xy", scale=1, width=0.004, alpha=0.75)
    ax.scatter(MU[:, 0], MU[:, 1], marker="x", s=70, c="crimson", zorder=3)
    ax.set_title(f"t = {t}"); ax.set_xlim(-4.5, 4.5); ax.set_ylim(-2.4, 2.4)
    ax.set_xticks([]); ax.set_yticks([])
plt.tight_layout(); plt.show()

The denoiser’s answer sharpens as noise falls. At t=0.99 every input maps to essentially one point, the dataset mean. At t=0.85 the answers split in two. At t=0.4 they are specific to the input.

Three regimes, worth naming because the rest of the post uses the terms:

noise level the denoiser’s answer is called
\(t \to 1\) the mean of the whole dataset global
middling \(t\) the centre of one blob mode-specific
\(t \to 0\) this particular point instance-specific

That changing specificity is the mechanism behind everything below.

Trajectories and the velocity field

To sample we must move, so we need a direction. Differentiate the path (1):

\[\frac{dz_t}{dt} = \epsilon - x \tag{7}\]

That is the velocity: the rate at which a state moves as noise increases. Sampling runs the other way, so we follow \(-dz_t/dt\).

We do not know the particular \(x\) and \(\epsilon\) behind a given state — that is exactly what the noise destroyed. So we substitute their conditional averages: \(\hat{x} = D(z_t,t)\) from Eq. (6), and the noise implied by the path, \(\hat{\epsilon} = (z_t - (1-t)\hat{x})/t\).

T_START, T_END = 1.0, 1e-3

def velocity(z, t):
    """Eq. (7) with conditional means substituted for the unknown x and ε."""
    t = torch.as_tensor(t, dtype=z.dtype).clamp(min=1e-3)
    x_hat = denoiser(z, t)
    eps_hat = (z - (1 - t) * x_hat) / t
    return eps_hat - x_hat

Attach that arrow to every point in space and you have a velocity field. A trajectory is the path traced by following it: forward, a clean sample dissolving into noise; reverse, a noise sample resolving into data. Generation walks a reverse trajectory.

The field is local in both arguments — it depends on where you are and on the current noise level. Two consequences, and both are why one evaluation cannot suffice: the arrow at your current position tells you nothing about the arrow further along, and the whole field is different at every \(t\).

Euler steps, and why more of them help

Following a continuous path exactly would take infinitely many evaluations. Instead we take finite jumps. The Euler step is the simplest: read the velocity where you are, assume it holds for the whole interval, move.

\[z_{t+h} = z_t + h \cdot v(z_t, t) \tag{8}\]

The assumption is written into the formula, and it is false whenever the velocity changes along the interval. The resulting error is discretisation error — the cost of replacing a curve with a finite number of straight jumps. It has nothing to do with the model being wrong.

This also answers how local estimates ever produce a globally correct result. Individually they do not. Each step is wrong by however much the velocity changed across it. Halve the step and you roughly halve that change, so error per step halves — but you take twice as many steps, so the total halves rather than quartering. That is first-order convergence: error \(\propto 1/n\).

def sample(z, n_steps):
    """Walk from t=1 to t≈0 in n_steps equal Euler moves — Eq. (8)."""
    ts = torch.linspace(T_START, T_END, n_steps + 1)
    for i in range(n_steps):
        z = z + (ts[i + 1] - ts[i]) * velocity(z, ts[i])
    return z

start = torch.randn(3000, 2)
reference = sample(start.clone(), 8000)      # fine enough to count as exact

print(f"{'steps':>7}{'mean distance from reference':>32}{'improvement':>14}")
prev = None
for n in [1, 2, 4, 8, 16, 32, 64, 128]:
    err = (sample(start.clone(), n) - reference).norm(dim=-1).mean().item()
    print(f"{n:>7}{err:>32.4f}{'' if prev is None else f'{prev/err:>13.2f}x'}")
    prev = err
  steps    mean distance from reference   improvement
      1                          2.0188
      2                          0.6895         2.93x
      4                          0.2768         2.49x
      8                          0.1416         1.95x
     16                          0.0720         1.97x
     32                          0.0363         1.98x
     64                          0.0182         2.00x
    128                          0.0090         2.01x

From eight steps onward the improvement is a clean factor of two per doubling: exactly the \(1/n\) law. That last column is the honest answer to “why fifty” — each doubling buys a halving, so you stop once the residual falls below what you can see.

The one-step catastrophe

Now the extreme case: a single Euler step across the entire journey, with an exact denoiser, against a well-resolved run.

If the model were the only thing that mattered, one step should give a blurry but recognisable version of the data — points somewhere near each blob. Watch instead.

careful = sample(start.clone(), 4000)
single = sample(start.clone(), 1)

def near_a_blob(x):
    """Fraction within 3σ of either centre. In 2-D the radial distance is
    Rayleigh, so the ceiling is 1 − exp(−9/2) ≈ 98.9%, not the 1-D 99.7%."""
    return ((x[:, None, :] - MU).norm(dim=-1).min(1).values < 3 * S).float().mean()

print(f"{'sampler':>12}{'near a blob':>14}{'spread':>10}{'mean position':>20}")
for name, out in [("4000 steps", careful), ("1 step", single)]:
    print(f"{name:>12}{near_a_blob(out):>13.1%}{out.std(0).max():>10.4f}"
          f"{f'({out[:,0].mean():+.3f}, {out[:,1].mean():+.3f})':>20}")
     sampler   near a blob    spread       mean position
  4000 steps        98.7%    2.0173    (-0.082, -0.001)
      1 step         0.0%    0.0010    (-0.000, -0.000)
Code 
fig, axes = plt.subplots(1, 2, figsize=(10, 3.0), sharex=True, sharey=True)
cloud = sample_data(2000)
for ax, (name, out) in zip(axes, [("4000 steps", careful), ("1 step", single)]):
    ax.scatter(cloud[:, 0], cloud[:, 1], s=6, alpha=0.15, label="data")
    ax.scatter(out[:800, 0], out[:800, 1], s=6, alpha=0.5, c="crimson", label="samples")
    ax.set_title(name); ax.set_xlim(-4, 4); ax.set_ylim(-2, 2)
    ax.set_xticks([]); ax.set_yticks([])
axes[0].legend(loc="upper left", fontsize=8)
plt.tight_layout(); plt.show()

Left: 4000 steps recovers both modes. Right: one step puts every sample on the origin — a location the data never occupies.

Not degraded — annihilated. Every sample landed on essentially the same point, the average of the two blobs, a location the data never occupies. In image terms: every prompt, every seed, the same grey mush. And the denoiser was exact.

Why: at \(t=1\) the honest answer is the average of everything

At \(t=1\) the state is pure noise. Under the independent pairing used in training it carries no information about which clean sample it was paired with. Asked “which \(x\) produced this?”, every \(x\) is equally consistent, and the squared-error-optimal answer is the average of all of them.

guesses = denoiser(torch.randn(50_000, 2), torch.tensor(1.0))
print(f"spread of outputs over 50,000 different noise inputs at t=1: "
      f"{(guesses.max(0).values - guesses.min(0).values).max():.2e}")
print(f"the output: ({guesses[0,0]:.3f}, {guesses[0,1]:.3f})   [= the dataset mean]")
spread of outputs over 50,000 different noise inputs at t=1: 0.00e+00
the output: (0.000, 0.000)   [= the dataset mean]

Zero spread. The single step was therefore not a rough approximation of a good move. It was an accurate move in the one direction that was genuinely correct at that instant, extended across the entire journey.

The denoiser gives a locally valid answer. A sampling step converts it into a finite-distance move. The gap between those is the entire cost of generation.

What changes along the way

If the answer at \(t=1\) is useless and the answer at \(t=0\) is the image, something happens in between. Two clarifications before measuring it.

The sampler creates no information. It is deterministic — \(z_1\) already fixes the output. What changes is how much of that the denoiser can resolve at the current noise level.

A single correlation would mislead. The blobs sit four units apart, so any estimate that gets the left/right choice right correlates strongly with the endpoint regardless of whether it knows anything finer. So we track two things separately: whether the denoiser has picked the correct blob, and whether it knows the offset within that blob.

def trajectory(z, n_steps=600):
    """Heun (predict, re-evaluate, average) — accurate enough that these count
    as the true trajectories rather than an approximation of them."""
    ts = torch.linspace(T_START, T_END, n_steps + 1)
    path = [z.clone()]
    for i in range(n_steps):
        h = ts[i + 1] - ts[i]
        v = velocity(z, ts[i])
        z = z + h * 0.5 * (v + velocity(z + h * v, ts[i + 1]))
        path.append(z.clone())
    return torch.stack(path), ts

torch.manual_seed(0)
paths, ts = trajectory(torch.randn(2000, 2))
final_x = paths[-1][:, 0]
which_blob = torch.where(final_x >= 0, 1.0, -1.0)
within_blob = final_x - which_blob * A

print(f"{'t':>6}{'correct blob':>15}{'knows offset':>15}")
for t_target in [0.95, 0.90, 0.80, 0.60, 0.40, 0.20]:
    k = int((ts - t_target).abs().argmin())
    g = denoiser(paths[k], ts[k])[:, 0]
    acc = (torch.where(g >= 0, 1.0, -1.0) == which_blob).float().mean()
    corr = torch.corrcoef(torch.stack([g - which_blob * A, within_blob]))[0, 1]
    print(f"{t_target:>6.2f}{acc:>14.0%}{corr:>15.2f}")
     t   correct blob   knows offset
  0.95          100%           0.11
  0.90          100%           0.19
  0.80          100%           0.40
  0.60          100%           0.77
  0.40          100%           0.90
  0.20          100%           1.00

(Starting at \(t=0.95\): at exactly 1 the guess is identically zero, which makes both columns meaningless.)

Two completely different timescales. After 5% of the journey the denoiser already knows which blob every trajectory will end in — and has almost no idea where within it. That second kind of knowledge accumulates across the rest of the trip.

The category is settled almost immediately; the specifics take most of the journey. If that holds for real image models — a two-blob toy, so treat it as a hypothesis — it explains why nudging a seed tends to change details while leaving composition intact.

Why one large step fails: the velocity changes

Discretisation error is exactly the amount the velocity changes across a step. That change has two independent parts, and both cost you.

Direction — the path bends. Measure it against the straight chord between each trajectory’s own endpoints: project every point on the path onto that chord, then take the perpendicular residual. The largest such residual, as a fraction of the chord length, says how far the path bows from the shortest available route.

start_pt, end_pt = paths[0], paths[-1]
chord = end_pt - start_pt
along = (((paths - start_pt[None]) * chord[None]).sum(-1)
         / chord.square().sum(-1).clamp_min(1e-8)[None])
perpendicular = (paths - (start_pt[None] + along[..., None] * chord[None])).norm(dim=-1)
bend = perpendicular.max(0).values / chord.norm(dim=-1)

speed = torch.stack([velocity(paths[k], ts[k]).norm(dim=-1).mean()
                     for k in range(0, 601, 10)])
print(f"median bend, as a fraction of the direct distance : {bend.median():.2f}")
print(f"speed varies along the path by a factor of        : {speed.max()/speed.min():.2f}x")
median bend, as a fraction of the direct distance : 0.17
speed varies along the path by a factor of        : 2.69x

The path bows by roughly a sixth of the endpoint distance, toward the mode it eventually commits to. A single step follows only the initial velocity, which points at the dataset average instead.

Code 
torch.manual_seed(3)
z0 = torch.tensor([[1.1, 0.9], [-0.7, -1.0], [0.4, 1.3]])
true_path, _ = trajectory(z0.clone())

fig, axes = plt.subplots(1, 2, figsize=(10, 3.5), sharex=True, sharey=True)
for ax, n in zip(axes, [1, 4]):
    tt = torch.linspace(T_START, T_END, n + 1)
    z, pts = z0.clone(), [z0.clone()]
    for i in range(n):
        z = z + (tt[i + 1] - tt[i]) * velocity(z, tt[i]); pts.append(z.clone())
    E = torch.stack(pts)
    for j in range(z0.shape[0]):
        ax.plot(true_path[:, j, 0], true_path[:, j, 1], lw=1.6, alpha=0.85, color="#2C4F6B")
        for i in range(n):
            d = E[i + 1, j] - E[i, j]
            ax.arrow(float(E[i, j, 0]), float(E[i, j, 1]), float(d[0]), float(d[1]),
                     head_width=0.09, length_includes_head=True, color="#B23A3A", alpha=0.9)
    ax.scatter(MU[:, 0], MU[:, 1], marker="x", s=70, c="k", zorder=3)
    ax.set_title(f"{n} Euler step{'s' if n > 1 else ''}")
    ax.set_xlim(-3.6, 3.6); ax.set_ylim(-1.9, 1.9); ax.set_xticks([]); ax.set_yticks([])
plt.tight_layout(); plt.show()

Three trajectories (blue) with Euler steps overlaid (red). One step leaves the data entirely and converges on the origin; four steps begin to track the curve.

Speed — the path is not covered at a uniform rate. This part is less intuitive, so worth being explicit: even a perfectly straight path defeats a single Euler step if it is traversed at a varying rate, because the step assumes one constant velocity and will overshoot or undershoot. Here the speed varies by a factor of about 2.7.

The precise condition is that one Euler step is exact only when the trajectory is affine in time — a straight line covered at constant speed. Geometric straightness alone is not enough. In 2-D that is easy to picture; for images the same statement holds in pixel or latent space, where “straight” means the intermediate states are a constant-rate blend of the endpoints.

What to carry away

  • The network answers a narrow question: given this noisy input, what is the average clean sample consistent with it. Not a plausible sample — the average.
  • How much it averages over depends on noise. At \(t=1\) the answer is the dataset mean, with zero spread across 50,000 inputs; then mode-specific; then instance-specific.
  • The one-step collapse is not model error. With a perfect denoiser one step still lands on the mean, because a locally correct direction was extended across the whole journey. That is discretisation error — a numerical problem with numerical answers.
  • Error falls as \(1/n\). Measured: a clean factor of two per doubling of steps. That is what makes “fifty” a budget rather than a magic number.
  • Sampling creates no information. The seed fixes the output; what changes is how much the denoiser can resolve.
  • Category first, detail later. The correct blob is decided within 5% of the trajectory; position within it takes most of the rest.
  • Two things vary along the way. The path bends by about a sixth of the endpoint distance and its speed changes by about 2.7×. One step is exact only if neither happens.

If one sentence survives: the model tells you how to move from the state you have now, and the sampler pays repeated network calls to turn those local answers into an accurate finite journey.

Where this goes next

The problem is now stated precisely: a changing velocity forces repeated evaluation, and repeated evaluation costs network calls. That makes the fast-sampling literature legible, because there are only three responses — follow the curve more cleverly, make the curve straighter, or learn a shortcut past it. The next post works through all three, including a result that catches people out: a mathematically better solver can be worse than a crude one at very low step counts.

References