Section C — Comparison table + code playbook

This section is a compact algorithm ↔ code playbook for three mainstream dynamics-based generative model families:

• Score-based diffusion (SDE): learn a time-dependent score \(s_\theta(x,t)\) and sample with a reverse-time SDE (often with a predictor–corrector scheme), or deterministically via the probability-flow ODE.
• Discrete diffusion (DDPM): learn a denoiser parameterization (commonly \(\varepsilon_\theta\)) and sample a reverse Markov chain.
• Flow matching (ODE): learn a velocity field \(v_\theta(x,t)\) by regression and sample by ODE integration.

To avoid duplication, we keep derivations elsewhere:

• Numerical ODE/SDE solvers: mathematics-foundation.md
• Continuity / Fokker–Planck and reverse-time constructions: mathematics-foundation.md and diffusion-models-intro.md
• Flow matching theory (CFM, Theorems 1–3, Gaussian paths): diffusion-models-intro.md and introduction-B-flow-matching.md

C0. Notation (shared)

We work in \(\mathbb{R}^D\).

• Data distribution: \(x_\text{data}\sim q(x)\) (we only have sample access).
• Simple prior: \(x_\text{prior}\sim p(x)\) (often \(p=\mathcal{N}(0,I)\)).
• Discrete time (DDPM): \(t\in\{1,\dots,T\}\).
• Continuous time (SDE/ODE): \(t\in[0,1]\).
• Model outputs:
• score: \(s_\theta(x,t)\approx \nabla_x\log p_t(x)\),
• noise: \(\varepsilon_\theta(x_t,t)\approx \varepsilon\) in the forward noising formula,
• velocity: \(v_\theta(x,t)\approx \dot\gamma_t\) along a chosen path \(\gamma_t\).

C1. Comparison table (what is modeled / learned / sampled)

Aspect	Score-based diffusion (SDE)	Discrete diffusion (DDPM)	Flow matching (ODE)
Dynamics type	Itô SDE \(dX_t=f(X_t,t)\,dt+g(t)\,dW_t\) (forward), sample via reverse-time SDE (or PF-ODE)	Forward Markov chain \(q(x_t\mid x_{t-1})\) with designed \(\beta_t\), sample via reverse chain	Deterministic ODE \(\frac{d}{dt}X_t=v_\theta(X_t,t)\)
What is designed	SDE coefficients \((f,g)\) / schedule (e.g., VP/VE)	discrete schedule \(\beta_t\) (equiv. \(\alpha_t,\bar\alpha_t\))	a probability path (bridge) \(\gamma_t\) and a time grid
What is learned	\(s_\theta(x,t)\) (score), or equivalent parameterization	typically \(\varepsilon_\theta(x_t,t)\) (noise), sometimes \(x_{0,\theta}\) or \(v_\theta\)	\(v_\theta(x,t)\) (velocity)
Training inputs	sample \((x_0,t)\), sample noise \(\varepsilon\), construct \(x_t\sim p_t(\cdot\mid x_0)\) in closed form	sample \((x_0,t)\), sample noise \(\varepsilon\), construct \(x_t=\sqrt{\bar\alpha_t}x_0+\sqrt{1-\bar\alpha_t}\,\varepsilon\)	sample endpoints (noise/data) and time \(t\), construct \(x_t=\gamma_t(x_0,x_1)\)
Supervision target	closed-form conditional score \(\nabla_{x_t}\log p(x_t\mid x_0)\) (DSM)	regression target \(\varepsilon\) in the forward corruption	target velocity \(\dot\gamma_t(x_0,x_1)\)
Sampling randomness	stochastic (reverse SDE), or deterministic (PF-ODE)	stochastic reverse chain (Gaussian transitions)	deterministic ODE (any ODE solver)
Typical solver	Euler–Maruyama + Langevin corrector (PC), or ODE solver for PF-ODE	ancestral time stepping	Euler / Heun / RK methods (see mathematics-foundation.md)

C2. Training pseudocode comparison (side-by-side)

C2.1 Score-based diffusion (denoising score matching)

Assume a forward SDE is chosen so that the conditional \(p(x_t\mid x_0)\) is tractable (often Gaussian). Then the DSM objective regresses the model score to the closed-form conditional score:

\[ \mathbb{E}_{t,x_0,x_t}\Big[\lambda(t)\,\big\|s_\theta(x_t,t)-\nabla_{x_t}\log p(x_t\mid x_0)\big\|^2\Big]. \]

Pseudocode:

repeat:
  sample x0 ~ p_data
  sample t  ~ Uniform(0, 1)
  sample ε  ~ N(0, I)
  construct xt ~ p(xt | x0)   # e.g. VE: x0 + σ(t) ε ; VP: √α̅(t) x0 + √(1-α̅(t)) ε
  set target score s*(xt,t) = ∇_{xt} log p(xt | x0)  # closed form
  loss = λ(t) · || sθ(xt,t) - s*(xt,t) ||^2
  take optimizer step on loss w.r.t. θ

C2.2 DDPM (noise-prediction loss)

DDPM training commonly regresses \(\varepsilon_\theta\) to the injected noise in the forward corruption:

repeat:
  sample x0 ~ p_data
  sample t  ~ Uniform({1, …, T})
  sample ε  ~ N(0, I)
  set xt = √α̅_t · x0 + √(1-α̅_t) · ε
  loss = || ε - εθ(xt, t) ||^2
  take optimizer step on loss w.r.t. θ

C2.3 Flow matching (velocity regression on a designed path)

Choose a bridge \(\gamma_t(x_0,x_1)\) connecting a prior endpoint and a data endpoint, and regress a velocity field to the path derivative \(\dot\gamma_t\):

\[ \mathbb{E}_{t,(x_0,x_1)}\Big[\|v_\theta(\gamma_t(x_0,x_1),t)-\dot\gamma_t(x_0,x_1)\|^2\Big]. \]

For the linear bridge with \(x_0\sim p\), \(x_1\sim q\),

\[ \gamma_t(x_0,x_1)=(1-t)x_0+t x_1,\qquad \dot\gamma_t(x_0,x_1)=x_1-x_0. \]

Pseudocode:

repeat:
  sample x1 ~ p_data
  sample x0 ~ N(0, I)
  sample t ~ Uniform(0, 1)
  xt = (1 - t) · x0 + t · x1
  target = x1 - x0
  loss = || vθ(xt, t) - target ||^2
  take optimizer step on loss w.r.t. θ

Markov vs continuous path sampling (practical point). In DDPM, training pairs \((x_t,t)\) are generated from a discrete Markov chain and indexed by integer \(t\). In score-based diffusion and flow matching, \((x_t,t)\) is sampled from a continuous-time path in one shot (via a tractable conditional perturbation or an explicit interpolation), without a Markov constraint between adjacent times.

C3. Sampling pseudocode comparison (side-by-side)

C3.1 Score-based diffusion: predictor–corrector reverse-SDE sampler

initialize x ~ N(0, I)  # terminal noise
for t = 1 → ε (descending):
  # Predictor: one reverse-SDE step (Euler–Maruyama)
  sample z ~ N(0, I)
  x ← x + [f(x,t) - g(t)^2 sθ(x,t)] Δt + g(t) √(|Δt|) · z

  # Corrector: K steps of Langevin at the same t (optional)
  repeat K times:
    sample z ~ N(0, I)
    x ← x + η_step · sθ(x,t) + √(2 η_step) · z
return x

The reverse-time formulas and the probability-flow ODE are summarized in diffusion-models-intro.md (Chapter 5.5) and in mathematics-foundation.md (Chapter 2).

C3.2 DDPM: ancestral reverse diffusion

initialize x_T ~ N(0, I)
for t = T, …, 1:
  sample z ~ N(0, I) if t > 1 else 0
  mean = 1/√α_t · (x_t - (1-α_t)/√(1-α̅_t) · εθ(x_t, t))
  x_{t-1} = mean + σ_t · z
return x_0

C3.3 Flow matching: ODE sampling by Heun (improved Euler)

set τ_0 = 0 < τ_1 < … < τ_S = 1
initialize x ~ N(0, I) at t = 0
for s = 0, …, S-1:
  t = τ_s,  t_next = τ_{s+1},  h = t_next - t
  k1 = vθ(x, t)
  x_pred = x + h · k1
  k2 = vθ(x_pred, t_next)
  x = x + (h/2) · (k1 + k2)
return x ≈ sample from p_data

C4. Minimal PyTorch code snippets (training + sampling)

These snippets are intentionally small and meant to make the variable ↔ formula correspondence explicit. For generic ODE/SDE solvers, see mathematics-foundation.md.

C4.1 Score-based diffusion (VE/VP schedules + DSM loss + PC sampler)

import torch


def ve_marginal_std(t, *, sigma_min=0.01, sigma_max=50.0):
    """
    VE schedule: choose a marginal std σ(t) (not the diffusion coefficient g(t)).
    t: shape (B,) in [0,1].
    """
    return sigma_min * (sigma_max / sigma_min) ** t


def vp_beta(t, *, beta_min=0.1, beta_max=20.0):
    """VP noise rate β(t), t in [0,1]."""
    return beta_min + t * (beta_max - beta_min)


def vp_alpha_bar(t, *, beta_min=0.1, beta_max=20.0):
    """
    α̅(t) = exp(-∫_0^t β(s) ds) for linear β(s).
    """
    integral = beta_min * t + 0.5 * (beta_max - beta_min) * t * t
    return torch.exp(-integral)


def dsm_loss_ve(score_model, x0, t, *, sigma_min=0.01, sigma_max=50.0, weight=None):
    """
    VE denoising score matching:
      x_t = x0 + σ(t) ε
      ∇_{x_t} log p(x_t|x0) = -(x_t-x0)/σ(t)^2 = -ε/σ(t)
    """
    eps = torch.randn_like(x0)
    sigma = ve_marginal_std(t, sigma_min=sigma_min, sigma_max=sigma_max)
    while sigma.ndim < x0.ndim:
        sigma = sigma[..., None]
    xt = x0 + sigma * eps
    target = -eps / sigma
    pred = score_model(xt, t)
    loss = (pred - target).pow(2)
    if weight is not None:
        w = weight(t)
        while w.ndim < loss.ndim:
            w = w[..., None]
        loss = w * loss
    return loss.mean()


@torch.no_grad()
def langevin_corrector(score_model, x, t, *, snr=0.16, n_steps=1):
    """
    K steps of overdamped Langevin dynamics at fixed t.
    Heuristic step size chosen to keep a target SNR.
    """
    for _ in range(n_steps):
        grad = score_model(x, t)
        noise = torch.randn_like(x)

        grad_norm = torch.norm(grad.reshape(grad.shape[0], -1), dim=-1).mean()
        noise_norm = torch.norm(noise.reshape(noise.shape[0], -1), dim=-1).mean()
        step_size = (snr * noise_norm / (grad_norm + 1e-12)) ** 2 * 2.0

        x = x + step_size * grad + torch.sqrt(2.0 * step_size) * noise
    return x


@torch.no_grad()
def pc_sampler(score_model, shape, *, n_steps=1000, sde="ve", device="cpu", corrector_steps=1, snr=0.16):
    """
    Predictor–corrector sampling.
    - VE predictor uses the identity g(t)^2 dt = d(σ(t)^2) for the chosen marginal σ(t).
    - VP predictor uses Euler–Maruyama with f(x,t)=-0.5β(t)x, g(t)=sqrt(β(t)).
    """
    x = torch.randn(shape, device=device)
    ts = torch.linspace(1.0, 1e-3, n_steps, device=device)  # reverse-time grid

    for i in range(len(ts) - 1):
        t = ts[i].expand(shape[0])
        t_next = ts[i + 1].expand(shape[0])
        dt = t_next - t  # negative

        z = torch.randn_like(x)

        if sde == "ve":
            sigma = ve_marginal_std(t)
            sigma_next = ve_marginal_std(t_next)
            while sigma.ndim < x.ndim:
                sigma = sigma[..., None]
                sigma_next = sigma_next[..., None]
            d_sigma2 = sigma_next**2 - sigma**2  # negative
            x = x - d_sigma2 * score_model(x, t) + torch.sqrt(torch.clamp(-d_sigma2, min=0.0)) * z

        elif sde == "vp":
            beta = vp_beta(t)
            while beta.ndim < x.ndim:
                beta = beta[..., None]
            f = -0.5 * beta * x
            g = torch.sqrt(beta)
            dt_b = dt
            while dt_b.ndim < x.ndim:
                dt_b = dt_b[..., None]
            x = x + (f - (g**2) * score_model(x, t)) * dt_b + g * torch.sqrt(-dt_b) * z

        else:
            raise ValueError("sde must be 've' or 'vp'")

        x = langevin_corrector(score_model, x, t_next, snr=snr, n_steps=corrector_steps)

    return x

C4.2 DDPM (noise-prediction training + ancestral sampler)

import torch


def ddpm_forward_sample(x0, t, alpha_bar):
    """
    x_t = sqrt(alpha_bar[t]) x0 + sqrt(1-alpha_bar[t]) eps
    alpha_bar: shape (T+1,) with alpha_bar[0]=1.
    t: int64 tensor of shape (B,) with values in {1,...,T}.
    """
    eps = torch.randn_like(x0)
    a_bar = alpha_bar[t]
    while a_bar.ndim < x0.ndim:
        a_bar = a_bar[..., None]
    xt = torch.sqrt(a_bar) * x0 + torch.sqrt(1.0 - a_bar) * eps
    return xt, eps


def ddpm_loss_eps(eps_model, x0, t, alpha_bar):
    xt, eps = ddpm_forward_sample(x0, t, alpha_bar)
    pred_eps = eps_model(xt, t)
    return (pred_eps - eps).pow(2).mean()


@torch.no_grad()
def ddpm_sample(eps_model, shape, *, alpha, alpha_bar, sigma):
    """
    Ancestral sampling for a fixed-variance DDPM parameterization.
    Arrays alpha/alpha_bar/sigma are shape (T+1,) with index 0 unused for sampling.
    """
    device = alpha.device
    T = alpha.shape[0] - 1
    x = torch.randn(shape, device=device)

    for t in range(T, 0, -1):
        t_vec = torch.full((shape[0],), t, device=device, dtype=torch.long)
        eps = eps_model(x, t_vec)

        mean = (1.0 / torch.sqrt(alpha[t])) * (x - ((1.0 - alpha[t]) / torch.sqrt(1.0 - alpha_bar[t])) * eps)

        if t > 1:
            x = mean + sigma[t] * torch.randn_like(x)
        else:
            x = mean

    return x

C4.3 Flow matching (linear bridge + ODE sampling)

import torch


def flow_matching_loss(v_model, x1, t):
    """
    Linear bridge between x0 ~ N(0,I) and x1 ~ data:
      x_t = (1-t) x0 + t x1
      target = d/dt x_t = x1 - x0
    """
    x0 = torch.randn_like(x1)
    t_for_path = t
    while t_for_path.ndim < x1.ndim:
        t_for_path = t_for_path[..., None]
    xt = (1.0 - t_for_path) * x0 + t_for_path * x1
    target = x1 - x0
    pred = v_model(xt, t)  # convention: model expects t with shape (B,)
    return (pred - target).pow(2).mean()


@torch.no_grad()
def flow_sample(v_model, shape, *, n_steps=50, method="heun", device="cpu"):
    """Sample by integrating dx/dt = v_model(x,t) from t=0 to t=1."""
    x = torch.randn(shape, device=device)
    ts = torch.linspace(0.0, 1.0, n_steps + 1, device=device)

    for i in range(n_steps):
        t = float(ts[i].item())
        t_next = float(ts[i + 1].item())
        h = t_next - t

        t_vec = torch.full((shape[0],), t, device=device)
        if method == "euler":
            x = x + h * v_model(x, t_vec)

        elif method == "heun":
            k1 = v_model(x, t_vec)
            x_pred = x + h * k1
            t_next_vec = torch.full((shape[0],), t_next, device=device)
            k2 = v_model(x_pred, t_next_vec)
            x = x + 0.5 * h * (k1 + k2)

        else:
            raise ValueError("method must be 'euler' or 'heun'")

    return x

C5. ELBO policy (keep it light)

DDPM was originally introduced via a variational bound (an ELBO over a latent diffusion chain), but most practical training uses a simplified regression loss (e.g., noise prediction MSE) that can be derived as a special case of that ELBO under common parameterizations and fixed variances.

For details, see introduction-beta-CFM-and-ELBO.md and diffusion-models-intro.md (Chapter 4).

C6. Section close template

1. We model: pick a view (score / noise / velocity).
2. We learn: the corresponding network output (\(s_\theta\) / \(\varepsilon_\theta\) / \(v_\theta\)).
3. We sample: the matching dynamics + numerical method (reverse SDE / reverse chain / ODE solver).