Diffusion Sampling Methods

A trained diffusion model is a denoiser: given a noisy image $x_t$ and its step $t$ , it predicts the noise $\epsilon_\theta(x_t, t)$ . A sampler turns it into images by choosing which noise levels to visit, how to move between them, and whether to inject fresh noise. Samplers trade the number of network evaluations (NFE), which dominates generation time, against sample quality, usually measured by the Fréchet Inception Distance (FID). All the samplers here reuse a network trained with the DDPM objective.

Problem

Given $\epsilon_\theta$ , trained with $T = 1000$ steps, draw images with as few evaluations as possible. In continuous time (Song et al., 2021) the target is fixed by the reverse SDE and the probability-flow ODE, which share the marginals $p_t$ and whose only unknown is the score $\nabla_x \log p_t(x) \approx -\epsilon_\theta(x_t, t)/\sqrt{1 - \bar\alpha_t}$ (see score matching). A sampler discretizes one of them; with an exact score, the error left is, in our reading, discretization error plus the small mismatch between $p_T$ and $\mathcal{N}(0, I)$ .

Inputs and Outputs

Input: $\epsilon_\theta$ (equivalently the clean-image estimate $\hat x_0$ or the score), its schedule, a step budget, and noise $x_T \sim \mathcal{N}(0, I)$ .
Output: an image $x_0$ ; deterministic samplers also give a map from $x_T$ to $x_0$ that can be run backward.

Intuition

Ancestral sampling is a random walk: each step denoises a little and adds new noise, which washes out errors but needs many small steps. Without the injected noise an ODE remains, whose smooth solution curves from noise to images can be followed with large steps.

Algorithm

All samplers run one loop:

Draw $x_T \sim \mathcal{N}(0, I)$ ; choose steps $\tau_S > \dots > \tau_1$ from $1, \dots, T$ , or continuous noise levels.
At each step, evaluate the network one or more times, form $\hat x_0$ and $\epsilon_\theta$ , and move to the next level, adding noise if the sampler is stochastic.
Return $\hat x_0$ from the last step.

Ancestral sampling (Ho et al., 2020, Algorithm 2; derived in the DDPM article) visits all $T$ steps and adds noise of variance $\sigma_t^2 = \beta_t$ or $\tilde\beta_t$ .

DDIM (Song, Meng & Ermon, 2021) replaces the update by

x_{\tau_{i-1}} = \sqrt{\bar\alpha_{\tau_{i-1}}}\, \hat x_0 + \sqrt{1 - \bar\alpha_{\tau_{i-1}} - \sigma_{\tau_i}^2}\; \epsilon_\theta(x_{\tau_i}, \tau_i) + \sigma_{\tau_i} z, \qquad \hat x_0 = \frac{x_{\tau_i} - \sqrt{1 - \bar\alpha_{\tau_i}}\, \epsilon_\theta(x_{\tau_i}, \tau_i)}{\sqrt{\bar\alpha_{\tau_i}}},

with $\bar\alpha_{\tau_0} = 1$ and $z \sim \mathcal{N}(0, I)$ (their Eq. 12; their $\alpha_t$ is this site’s $\bar\alpha_t$ ). The noise is set by one parameter $\eta \ge 0$ (Eq. 16):

\sigma_{\tau_i}(\eta) = \eta \sqrt{\frac{1 - \bar\alpha_{\tau_{i-1}}}{1 - \bar\alpha_{\tau_i}}} \sqrt{1 - \frac{\bar\alpha_{\tau_i}}{\bar\alpha_{\tau_{i-1}}}} .

With $\eta = 1$ this is DDPM’s ancestral step with $\sigma_t^2 = \tilde\beta_t$ on the chosen steps; with $\eta = 0$ no noise is added after $x_T$ and the sampler is deterministic.

Mathematical Formulation

Why DDIM needs no retraining

Song, Meng and Ermon defined non-Markovian forward processes, Gaussian $q_\sigma(x_{t-1} \mid x_t, x_0)$ of mean $\sqrt{\bar\alpha_{t-1}}\, x_0 + \sqrt{1 - \bar\alpha_{t-1} - \sigma_t^2}\,(x_t - \sqrt{\bar\alpha_t}\, x_0)/\sqrt{1 - \bar\alpha_t}$ and variance $\sigma_t^2$ (Eq. 7), which keep every marginal $q(x_t \mid x_0) = \mathcal{N}(\sqrt{\bar\alpha_t}\, x_0, (1 - \bar\alpha_t) I)$ . The DDPM loss depends only on these marginals, so for $\sigma_t > 0$ each process’s variational objective is a weighted noise-prediction loss plus a constant (Theorem 1; $\eta = 0$ is the limit), and a network trained with Ho et al.’s simple loss serves all of them (exactly so if its parameters were not shared across $t$ ). DDIM substitutes $\hat x_0$ for $x_0$ in this posterior. The argument holds on any subsequence $\tau$ (§4.2, Appendix C.1), which licenses skipping steps; the authors spaced $\tau$ quadratically for CIFAR-10 and linearly elsewhere (Appendix D.2).

The ODE behind DDIM

Write $\nu_t = \sqrt{(1 - \bar\alpha_t)/\bar\alpha_t}$ , the inverse square root of the signal-to-noise ratio, and $\bar x = x_t / \sqrt{\bar\alpha_t} = x_0 + \nu_t\, \epsilon$ . The $\eta = 0$ update becomes $\bar x_{\tau_{i-1}} = \bar x_{\tau_i} + (\nu_{\tau_{i-1}} - \nu_{\tau_i})\, \epsilon_\theta$ (Eq. 13), an Euler step for

\frac{d\bar x}{d\nu} = \epsilon_\theta\Big(\frac{\bar x}{\sqrt{1 + \nu^2}},\ \nu\Big) = \frac{\bar x - \hat x_0}{\nu}

(Eq. 14), which Song, Meng and Ermon showed is equivalent to the probability-flow ODE of the variance-exploding SDE (Proposition 1). Karras et al. (2022, §3) reached the same form by setting the noise level equal to time, $\sigma(t) = t$ (their $\sigma$ is our $\nu$ ), without rescaling: $dx/dt = (x - D(x; t))/t$ with denoiser $D = \hat x_0$ , the choice DDIM makes. The ODE’s solution curves are nearly straight at high and low noise and bend in between, where the denoiser’s output changes between images; one Euler step to zero noise lands on $\hat x_0$ .

Euler, Heun and the step schedule

Euler’s method has local error $O(h^2)$ in the step size $h$ . Heun’s method (Karras et al., Algorithm 1) takes the Euler step, evaluates the slope at its end, and redoes the step with the average slope, reaching $O(h^3)$ for one extra evaluation; the last step, to zero noise, stays Euler (§3), so $N$ steps cost $2N - 1$ evaluations. Karras et al. placed the levels by

\nu_i = \Big(\nu_{\max}^{1/\rho} + \frac{i}{N - 1}\big(\nu_{\min}^{1/\rho} - \nu_{\max}^{1/\rho}\big)\Big)^{\rho}, \qquad i = 0, \dots, N - 1, \qquad \nu_N = 0

(Eq. 5, where they write $\sigma$ ). Large $\rho$ concentrates steps at low noise; $\rho = 3$ nearly equalized the truncation error per step, but 5 to 10 gave better images, and they chose $\rho = 7$ , with $\nu_{\min} = 0.002$ and $\nu_{\max} = 80$ for their models (Table 1).

Exponential integrators in log-SNR

Lu et al. (2022) solved the linear term of the probability-flow ODE exactly and changed variable to $\lambda = \tfrac{1}{2} \log \mathrm{SNR} = \log\big(\sqrt{\bar\alpha_t}/\sqrt{1 - \bar\alpha_t}\big)$ , which gives, in this site’s notation,

x_t = \sqrt{\frac{\bar\alpha_t}{\bar\alpha_s}}\, x_s - \sqrt{\bar\alpha_t} \int_{\lambda_s}^{\lambda_t} e^{-\lambda}\, \hat\epsilon_\theta(\hat x_\lambda, \lambda)\, d\lambda

(Proposition 3.1). DPM-Solver- $k$ Taylor-expands $\hat\epsilon_\theta$ in $\lambda$ to order $k - 1$ , costs $k$ evaluations per step, and has global error $O(h^k)$ in the largest $\lambda$ step (Theorem 3.2). The first-order solver is exactly DDIM (§4.1), which, Lu et al. argue, explains why DDIM beats Euler steps in $t$ . Since $\lambda = -\log \nu$ , $e^{-\lambda} d\lambda = -d\nu$ , and dividing by $\sqrt{\bar\alpha_t}$ gives $\bar x_t = \bar x_s + \int_{\nu_s}^{\nu_t} \epsilon_\theta\, d\nu$ : in our reading, DPM-Solver integrates the same ODE as DDIM and EDM and differs in where it evaluates the network, with steps uniform in $\lambda$ (§3.3).

Parameters

Step count. The main cost knob; Heun and DPM-Solver-2 or -3 spend 2 or 3 evaluations per step, so compare samplers at equal NFE.
Injected noise. $\eta \in [0, 1]$ for DDIM. EDM’s stochastic sampler (Algorithm 2) first raises the noise level by a factor $1 + \gamma_i$ , $\gamma_i = \min(S_{\text{churn}}/N, \sqrt{2} - 1)$ inside a window $[S_{t\min}, S_{t\max}]$ , by adding noise scaled by $S_{\text{noise}}$ , then takes a Heun step down.

Complexity

Time is NFE network passes; the solver’s arithmetic is negligible (Lu et al., 2022, Appendix E.7). Multistep solvers also store earlier evaluations.

Implementation

The NumPy script reuses the diffusion model article’s toy data: four $16 \times 16$ shapes (disk, square, plus sign, diagonal bar) with weights 0.4, 0.3, 0.2 and 0.1 and pixel noise of standard deviation 0.1. Every noisy marginal is a Gaussian mixture, so $\epsilon^\ast$ is exact, and all samplers start from the same 1,000 noise images, so they differ only by discretization. It prints the share of samples nearest each shape and the RMS distance to it, and for deterministic samplers the median distance from a 1,000-step Heun reference and the fraction whose nearest shape differs from the reference’s; then it tests determinism and inversion.

import sys
import numpy as np

seed = int(sys.argv[1]) if len(sys.argv) > 1 else 0

T = 1000
beta = np.linspace(1e-4, 0.02, T)               # linear schedule of Ho et al. (2020)
ab = np.cumprod(1 - beta)                       # ab[t-1] = alpha_bar_t
sig = np.sqrt((1 - ab) / ab)                    # nu_t of the text (EDM's sigma)

# Data of the diffusion model article: four 16x16 shapes (disk, square, plus, diagonal bar)
# in [-1, 1], x0 = mu_k + 0.1 * noise, with mixture weights 0.4, 0.3, 0.2, 0.1.
yy, xx = np.mgrid[:16, :16]
shapes = [(yy - 7.5) ** 2 + (xx - 7.5) ** 2 < 5.5 ** 2,
          (abs(yy - 7.5) < 5) & (abs(xx - 7.5) < 5),
          (abs(yy - 7.5) < 2) | (abs(xx - 7.5) < 2),
          abs(yy - xx) < 2.5]
mu = np.array([np.where(m, 1.0, -1.0).ravel() for m in shapes])   # (4, 256)
w = np.array([0.4, 0.3, 0.2, 0.1])
tau2 = 0.1 ** 2

# In the variable xbar = x_t / sqrt(ab_t) = x0 + sigma * eps the noisy density is a
# mixture of N(mu_k, (tau2 + sigma^2) I), so the optimal noise prediction
# eps*(xbar; sigma) = E[eps | xbar] = (xbar - E[x0 | xbar]) / sigma is exact.
def eps_star(xb, s):
    v = tau2 + s ** 2
    d2 = (xb ** 2).sum(1, keepdims=True) - 2 * xb @ mu.T + (mu ** 2).sum(1)
    lg = np.log(w) - d2 / (2 * v)
    r = np.exp(lg - lg.max(1, keepdims=True))
    r /= r.sum(1, keepdims=True)
    return (xb - r @ mu) * (s / v)              # well defined at sigma = 0, where it is 0

def grid(S):                                    # S evenly spaced steps of 1..T, ending at T
    return np.round(np.linspace(1, T, S)).astype(int)

def ddim(xT, S, eta, rng=None, ancestral=False):
    """DDIM update (Song, Meng & Ermon, 2021, Eq. 12) on S grid steps. eta = 1 is DDPM
    with sigma_t^2 = beta~_t; ancestral=True adds noise of variance beta_t (Ho et al.)."""
    tau, x = grid(S), xT.copy()
    for i in range(S - 1, -1, -1):
        a = ab[tau[i] - 1]
        a_prev = ab[tau[i - 1] - 1] if i > 0 else 1.0
        e = eps_star(x / np.sqrt(a), sig[tau[i] - 1])
        x0_hat = (x - np.sqrt(1 - a) * e) / np.sqrt(a)
        s_eta = eta * np.sqrt((1 - a_prev) / (1 - a) * (1 - a / a_prev))
        x = np.sqrt(a_prev) * x0_hat + np.sqrt(1 - a_prev - s_eta ** 2) * e
        noise_sd = np.sqrt(1 - a / a_prev) if (ancestral and i > 0) else s_eta
        if noise_sd > 0:                        # eta = 0 uses no randomness after x_T
            x = x + noise_sd * rng.standard_normal(x.shape)
    return x

def edm_sigmas(N, rho=7.0):                     # Karras et al. (2022), Eq. 5, then sigma = 0
    smax, smin, i = sig[-1], sig[0], np.arange(N)
    s = (smax ** (1 / rho) + i / (N - 1) * (smin ** (1 / rho) - smax ** (1 / rho))) ** rho
    return np.append(s, 0.0)

def ode(xT, sigmas, heun):
    """Euler or Heun (Karras et al., Algorithm 1 with sigma(t) = t, s(t) = 1) on the
    probability-flow ODE d xbar / d sigma = eps*(xbar; sigma)."""
    x = xT / np.sqrt(ab[-1])                    # the same starting noise, as xbar
    for s, s_next in zip(sigmas[:-1], sigmas[1:]):
        d = eps_star(x, s)
        x_next = x + (s_next - s) * d
        if heun and s_next > 0:                 # 2nd-order correction, skipped at sigma = 0
            x_next = x + (s_next - s) * 0.5 * (d + eps_star(x_next, s_next))
        x = x_next
    return x                                    # at sigma = 0, xbar = x

def nearest(x):
    d2 = ((x[:, None, :] - mu[None]) ** 2).mean(-1)
    return d2.argmin(1), np.sqrt(d2.min(1))     # nearest shape, RMS distance to it

def dist(a, b):                                 # per-image RMS distance
    return np.sqrt(((a - b) ** 2).mean(1))

n = 1000
xT = np.random.default_rng(seed).standard_normal((n, 256))     # one x_T for every sampler
ref = ode(xT, edm_sigmas(1000), heun=True)                     # 1,000-step Heun reference
k_ref = nearest(ref)[0]
print(f"seed {seed}. 1000-step DDIM eta=0 vs reference: median distance "
      f"{np.median(dist(ddim(xT, 1000, 0.0), ref)):.1e}")
same = np.append(sig[grid(20)[::-1] - 1], 0.0)                 # DDIM's 20 noise levels
print(f"max |DDIM eta=0 - Euler on the same 20 noise levels| = "
      f"{np.abs(ddim(xT, 20, 0.0) - ode(xT, same, heun=False)).max():.1e}\n")

print(f"{'sampler':15s}{'NFE':>5s}{'disk':>7s}{'square':>7s}{'plus':>7s}{'bar':>7s}"
      f"{'RMS':>7s}{'to ref':>8s}{'shape':>7s}")
print(f"{'data':15s}{'':5s}" + "".join(f"{v:7.3f}" for v in w) + f"{0.1:7.3f}")
rng = lambda: np.random.default_rng(seed + 1)
for S in [10, 20, 50, 100]:
    for name, nfe, f in [("DDPM ancestral", S, lambda: ddim(xT, S, 1.0, rng(), ancestral=True)),
                         ("DDIM eta=1", S, lambda: ddim(xT, S, 1.0, rng())),
                         ("DDIM eta=0", S, lambda: ddim(xT, S, 0.0)),
                         ("Euler, rho=7", S, lambda: ode(xT, edm_sigmas(S), heun=False)),
                         ("Heun, rho=7", S - 1, lambda: ode(xT, edm_sigmas(S // 2), heun=True))]:
        x = f()
        k, d = nearest(x)
        line = (f"{name:15s}{nfe:5d}" + "".join(f"{v:7.3f}" for v in np.bincount(k, minlength=4) / n)
                + f"{np.sqrt((d ** 2).mean()):7.3f}")
        if name not in ("DDPM ancestral", "DDIM eta=1"):        # deterministic samplers
            line += f"{np.median(dist(x, ref)):8.4f}{np.mean(k != k_ref):7.1%}"
        print(line)
k, d = nearest(ddim(xT, 1000, 1.0, rng(), ancestral=True))
print(f"{'DDPM ancestral':15s}{1000:5d}" + "".join(f"{v:7.3f}" for v in np.bincount(k, minlength=4) / n)
      + f"{np.sqrt((d ** 2).mean()):7.3f}")

# Determinism: two runs from the same x_T with different noise generators.
print()
for eta in [0.0, 1.0]:
    x1 = ddim(xT, 50, eta, np.random.default_rng(100))
    x2 = ddim(xT, 50, eta, np.random.default_rng(200))
    print(f"DDIM eta={eta:.0f}, 50 steps, same x_T, two noise draws: max difference "
          f"{np.abs(x1 - x2).max():.2f}, nearest shape differs for "
          f"{np.mean(nearest(x1)[0] != nearest(x2)[0]):.1%}")

# Inversion: run the eta = 0 update upward in sigma (an Euler step with eps at the
# current point) from data images to x_T, then generate from that x_T with the same S.
def invert(x0, S):
    sigmas = np.append(0.0, sig[grid(S) - 1])
    x = x0.copy()
    for s, s_next in zip(sigmas[:-1], sigmas[1:]):
        x = x + (s_next - s) * eps_star(x, s)
    return x * np.sqrt(ab[-1])                  # x_T = sqrt(ab_T) * xbar_T

r = np.random.default_rng(seed + 2)
x0 = mu[r.choice(4, n, p=w)] + 0.1 * r.standard_normal((n, 256))
print("\ninversion of 1,000 data images, x0 -> x_T -> x0' with S steps each way")
print(f"{'S':>5s}{'median RMS(x0 - x0prime)':>26s}{'shape changed':>15s}{'std of x_T':>12s}")
for S in [10, 20, 50, 100, 1000]:
    z = invert(x0, S)
    x0p = ddim(z, S, 0.0)
    print(f"{S:5d}{np.median(dist(x0, x0p)):26.4f}"
          f"{np.mean(nearest(x0)[0] != nearest(x0p)[0]):15.1%}{z.std():12.3f}")

Output:

seed 0. 1000-step DDIM eta=0 vs reference: median distance 1.2e-03
max |DDIM eta=0 - Euler on the same 20 noise levels| = 4.7e-14

sampler          NFE   disk square   plus    bar    RMS  to ref  shape
data                  0.400  0.300  0.200  0.100  0.100
DDPM ancestral    10  0.461  0.304  0.158  0.077  0.353
DDIM eta=1        10  0.457  0.309  0.160  0.074  0.027
DDIM eta=0        10  0.424  0.312  0.169  0.095  0.035  0.0655   5.6%
Euler, rho=7      10  0.434  0.312  0.166  0.088  0.068  0.0337   7.2%
Heun, rho=7        9  0.385  0.316  0.189  0.110  0.429  0.2604   2.6%
DDPM ancestral    20  0.408  0.320  0.172  0.100  0.196
DDIM eta=1        20  0.405  0.327  0.170  0.098  0.046
DDIM eta=0        20  0.402  0.310  0.182  0.106  0.054  0.0457   2.6%
Euler, rho=7      20  0.408  0.307  0.181  0.104  0.082  0.0179   3.3%
Heun, rho=7       19  0.386  0.306  0.195  0.113  0.143  0.0433   0.3%
DDPM ancestral    50  0.398  0.304  0.195  0.103  0.129
DDIM eta=1        50  0.398  0.305  0.195  0.102  0.068
DDIM eta=0        50  0.392  0.308  0.190  0.110  0.078  0.0214   1.0%
Euler, rho=7      50  0.394  0.306  0.190  0.110  0.092  0.0074   1.2%
Heun, rho=7       49  0.387  0.304  0.195  0.114  0.104  0.0049   0.0%
DDPM ancestral   100  0.383  0.312  0.205  0.100  0.113
DDIM eta=1       100  0.380  0.315  0.205  0.100  0.079
DDIM eta=0       100  0.390  0.306  0.192  0.112  0.088  0.0112   0.5%
Euler, rho=7     100  0.390  0.307  0.192  0.111  0.096  0.0037   0.6%
Heun, rho=7       99  0.387  0.304  0.195  0.114  0.101  0.0011   0.0%
DDPM ancestral  1000  0.384  0.322  0.200  0.094  0.101

DDIM eta=0, 50 steps, same x_T, two noise draws: max difference 0.00, nearest shape differs for 0.0%
DDIM eta=1, 50 steps, same x_T, two noise draws: max difference 2.38, nearest shape differs for 69.9%

inversion of 1,000 data images, x0 -> x_T -> x0' with S steps each way
    S  median RMS(x0 - x0prime)  shape changed  std of x_T
   10                    0.0886          11.2%       0.346
   20                    0.0709           2.5%       0.541
   50                    0.0388           1.1%       0.782
  100                    0.0220           0.9%       0.883
 1000                    0.0033           0.2%       0.983

The second line confirms that DDIM with $\eta = 0$ is Euler’s method. At 10 evaluations the samplers fail differently: ancestral DDPM leaves samples noisy (RMS 0.353 against the data’s 0.1), DDIM over-smooths them (0.027 for $\eta = 1$ , 0.035 for $\eta = 0$ ), and both stochastic samplers over-produce the common disk and under-produce the rare bar. Heun with 5 steps has not converged, but from about 50 evaluations it leads: its median distance from the reference falls fourfold from 49 to 99 NFE, the second-order rate, to about that of 1,000 DDIM steps, while DDIM’s error at best halves when its steps double. Euler on $\rho = 7$ levels beats DDIM’s uniform steps in median distance but changes more shapes, in our reading because its steps are longer at high noise, where the shape is decided.

From one $x_T$ , $\eta = 0$ gives identical images, while $\eta = 1$ picks a different shape for 69.9%, near the $1 - \sum_k w_k^2 = 0.70$ expected for independent runs. Inversion’s median error is 0.0033 with 1,000 steps but 0.0886 with 10, whose latents have standard deviation 0.346, unlike Gaussian noise. Over seeds 0 to 4 (the script’s argument), Heun’s median distance at 99 NFE was always 0.0011, the deterministic samplers at 99 or 100 NFE changed at most 0.7% of the reference’s shapes, and inversion changed the shape of 9.3% to 11.2% of images at 10 steps and 0.0% to 0.2% at 1,000.

Properties and Behavior

NFE versus FID. On CIFAR-10, DDIM’s FID was 13.36, 6.84, 4.67, 4.16 and 4.04 with 10, 20, 50, 100 and 1,000 steps; with $\eta = 1$ it was 41.07, 18.36, 8.01, 5.78 and 4.73, and ancestral sampling with $\sigma_t^2 = \beta_t$ gave 367.43 at 10 steps but the best value, 3.17, at 1,000 (Table 1, identical in arXiv v1, v2 and v4). The abstract claims wall-clock speed-ups of 10 to 50 times, §5 10 to 100. Karras et al.’s changes cut the NFE for near-best FID by 7.3, 300 and 3.2 times on three pretrained models (§3), and their retrained networks reached CIFAR-10 FIDs of 1.79 (class-conditional) and 1.97 (unconditional) at 35 NFE (Table 2, arXiv v2). DPM-Solver reported 4.70 at 10 NFE and 2.87 at 20 on CIFAR-10, the two with different end times, $10^{-3}$ and $10^{-4}$ (Table 3, arXiv v3).
Stochastic versus deterministic. Both reach the same distribution in theory, but Karras et al. argue that injected noise corrects earlier errors: on ImageNet $64 \times 64$ their stochastic sampler reached 1.55 at 511 NFE against 2.64 at 79 deterministically (Tables 3–4, arXiv v2). Their text starts from ADM’s published 2.07; Table 4 lists 2.01 for the original sampler. Too much churn lost detail and oversaturated colors, and the settings needed a per-model grid search (§4).
Determinism and inversion. A deterministic sampler draws all its diversity from $x_T$ . Trajectories of different lengths from one $x_T$ kept the same high-level features, making interpolation in $x_T$ meaningful, and encoding then decoding CIFAR-10 test images gave a per-dimension mean squared error (pixels in $[0, 1]$ ) of 0.014 with 10 steps and 0.0001 with 1,000 (Song, Meng & Ermon, §5.2–5.4, Table 2 of arXiv v2 and later). Mokady et al. (2023) note that inversion works for unconditional models, but the large classifier-free guidance scales used for editing with Stable Diffusion amplify its error; they take the unguided inversion as a pivot and optimize the unconditional text embedding.

Limitations

Strong guidance destabilizes high-order solvers; Lu et al. (2025) report earlier ones becoming slower than DDIM at large guidance scales.
Step placement, churn and the last step are tuned per model, so a reported FID reflects the sampler as well as the network.
Inversion is exact only in the small-step limit; with large steps the toy example regenerates some images as other shapes.

Variants

Predictor–corrector samplers (Song et al., 2021, §4.2) follow each reverse-SDE step with Langevin steps toward $p_t$ ; on CIFAR-10, one corrector step per predictor step always beat the predictor alone.

DPM-Solver++ (Lu et al., 2025, read in arXiv v3) integrates the $\hat x_0$ prediction, adds thresholding and a multistep variant, and reports good guided samples in 15 to 20 steps.

UniPC (Zhao et al., 2023) adds a corrector that raises the order of accuracy without extra evaluations and reduces to DDIM at first order (§3.3).

Distillation trains a new network: progressive distillation reached a CIFAR-10 FID of 3.00 in 4 steps (Salimans & Ho, 2022, Table 2, arXiv v2), and consistency models distilled from a diffusion model reached one-step FIDs of 3.55 on CIFAR-10 and 6.20 on ImageNet $64 \times 64$ (Song et al., 2023). Latent diffusion instead makes each evaluation cheaper.