Tutorial: Generative models as transport

This post is a companion to the hands-on tutorial at the 2026 IAIFI summer school; the notebooks, code, and setup instructions live in the tutorial repository.

A generative model in modern machine learning turns noise into structure. The goal of this tutorial is a hands-on understanding of how that works, from exact change of variables to diffusion models trainable at scale. From there, separate tracks explore the same picture under different constraints: regress the velocity instead of the score (flow matching), keep the exact density because the physics demands it (lattice field theory), or compress the whole transport into a single network call (distillation).

Part I The core: from maps to diffusion

How do we generate samples?

Physicists have been generating samples for many decades using Markov-chain Monte Carlo (MCMC): random walks whose stationary distribution is a known unnormalized density, often called energy in this context (Krauth, 1996). The same lineage produced energy-based models such as restricted Boltzmann machines, where an energy function is learned and then sampled via Markov chains. Modern generative modeling differs from classic MCMC in two important ways. First, what specifies the target: physics and inference problems give us explicit (unnormalized) log-densities, while in ML the target is usually known only through data. Second, and more radically, the energy frame of mind is often dropped altogether: instead of learning a log-density and then sampling it, learn the sampler itself, as a network carrying a latent variable $z$ to a data point. The generative adversarial network (GAN) takes this literally: $x = G(z)$ , a deterministic pushforward of noise trained adversarially, with no density we can evaluate (Goodfellow et al., 2014). The variational autoencoder (VAE) is subtly different. Its decoder maps $z$ not to a sample but to the parameters of a small output distribution (say, a Gaussian mean), and the sample is drawn around that; the model density is then an integral over all latents, intractable pointwise, so training goes through a variational bound (Kingma & Welling, 2014). This post follows a third route: transport a tractable base density into the target while keeping a usable density — the map GANs had, plus the density they lacked.

Generative modeling as transport

Fix a base density $\rho$ we can sample cheaply, e.g. a standard Gaussian, and a map $f_\theta$ . Define the model by construction: draw noise $z \sim \rho$ , output $x = f_\theta(z)$ . The distribution $q_\theta$ of the outputs is the pushforward of $\rho$ through $f_\theta$ , and all of generative modeling in this post is the task of shaping $f_\theta$ so that $q_\theta$ matches a target $p$ . Two capabilities must be kept apart, because methods make different trade-offs: sampling from $q_\theta$ (built in, by construction) and evaluating the density $q_\theta(x)$ at a given point (not automatic at all). The distinction matters because it decides which losses we can write down. Suppose we can express $q_\theta$ and sample from it — then what? Training means comparing $q_\theta$ to $p$ through some divergence or loss function, to be minimized. The classic choices are the forward or reverse varaints of the Kullback-Leibler (KL) divergence. Forward KL takes “ground truth” samples from $p$ , and reduces to maximum likelihood, which has the benefit that only $q_\theta(x)$ but not $p(x)$ must be evalauted. Reverse KL runs on draws from the model, scored under both the model and the target. Note that both need exactly the density evaluation of the model that is “not automatic”.

\text{Forward KL:} \quad D_{KL}(p \parallel q_\theta) = \int p(x) \log \frac{p(x)}{q_\theta(x)} \dd{x}

\text{Reverse KL:} \quad D_{KL}(q_\theta \parallel p) = \int q_\theta(x) \log \frac{q_\theta(x)}{p} \dd{x}

Exercise forward KL is maximum likelihood

Show that minimizing the forward KL over $\theta$ is equivalent to maximizing the likelihood of the data under $q_\theta$ , and explain why this never requires evaluating the data density $p(x)$ .

Hint

Split $\int p \log (p / q_\theta)$ into two terms; only one depends on $\theta$ .

objective	draw from	evaluate q(x)	evaluate p(x)
forward KL (max. likelihood)	p (target)	✓	— (const.)
reverse KL (variational)	q (model)	✓	✓ (up to Z)

Figure 1. The same target distribution two ways: a cloud of draws x ∼ p (left) and its density p(x) (right). Hover a sample to ring the same point in both panels — a draw is a location, the density is its value there. The table is why the choice matters: matching p by maximum likelihood (forward KL) only needs draws from p and the ability to score them under the model q; the reverse direction, used when the target is given as an (energy) density, instead needs draws from q, scored under both q and p.

So when is the density available? For a general map, evaluating $q_\theta(x)$ means marginalizing over every noise value that $f_\theta$ sends to $x$ — hopeless in general. When $f_\theta$ is bijective, however, every output $x$ has a unique preimage $z = f_\theta^{-1}(x)$ . If we can compute the inverse, we can then evaluate the probability of where we started from. We also have to account for how the map stretches volume (pulling nearby points in, or pushing them away), giving us the change-of-variables formula

\log q_\theta(x) = \log \rho(z) - \log \bigl| \det \partial_z f_\theta(z) \bigr| , \qquad z = f_\theta^{-1}(x) .\tag{1}

morph t1.00

γ (center)0.50

σ (width)0.80

λ (strength)3.60

Figure 2. A Gaussian (grey, dashed) deformed by a single cubic-rational map. Relabelling the prior gives ρ ∘ f⁻¹ (red, dashed); dividing by the local stretch |h′| gives the true density q (filled). The arrows under each sample point outward and red where the map pulls points apart and the density thins, and inward and blue where it compresses and the density piles up.

Building networks that are invertible and have a tractable Jacobian determinant is its own craft, yielding normalizing flows built from coupling or autoregressive layers; see e.g. the analytic bijections post for notes on that. In this tutorial we instead take the continuous route.

Continuous transport

Instead of trying to express $f_\theta$ as one big neural network, let a velocity field $v(x, t)$ generate the map by integration:

\dv{t} x(t) = v(x(t), t) ,\tag{2}

started at $x(0) = z \sim \rho$ and integrated to $t = 1$ ; the endpoint $x(1)$ is the model sample. Solving the ODE deforms the base density smoothly through a family of intermediate densities $p_t$ , and this evolution obeys two equivalent statements of the change of variables. We will need both later, so we write each with explicit arguments:

\pdv{t} p_t(x) + \nabla \cdot \bigl( p_t(x) \, v(x, t) \bigr) = 0 ,\tag{3}

the Eulerian (continuity-equation) form governing the density field at a fixed point, and

\dv{t} \log p_t(x(t)) = - (\nabla \cdot v)(x(t), t) ,\tag{4}

the Lagrangian form governing the log-density carried along a trajectory of (2). The Lagrangian form (4) is what a continuous normalizing flow integrates to report exact densities; the Eulerian form (3) is what the fixed noising paths of the next section satisfy.

Exercise continuity from an infinitesimal change of variables

Derive the Eulerian continuity equation (3) from the discrete change of variables (1) by taking one infinitesimal step $f(x) = x + \epsilon \, v(x)$ and expanding to first order in $\epsilon$ .

Hint

$\log \bigl| \det (I + \epsilon \, \partial_x v) \bigr| = \epsilon \, \nabla \cdot v + O(\epsilon^2)$ ; feed this into (1), divide by $\epsilon$ , and let $\epsilon \to 0$ .

Exercise the two forms are equivalent

Show directly that the Lagrangian form (4) and the Eulerian continuity equation (3) say the same thing.

Hint

Expand the total derivative along a trajectory: $\dv{t} \log p_t(x(t)) = \partial_t \log p_t + v \cdot \nabla \log p_t$ .

Both statements are the deterministic case of one more general object. Add to the transport an independent Brownian motion noise source, i.e. promote the ODE (2) to a stochastic differential equation (SDE):

\dd{x} = v(x, t) \, \dd{t} + g_t \, \dd{W}_t \,,\tag{5}

where $W_t$ is a standard Wiener process (Brownian motion), so $\dd{W}_t$ is a Gaussian increment of variance $\dd{t}$ and $g_t$ sets the noise scale. This is the equation of motion of a single sample; averaging over the ensemble of its solutions, the density $p_t$ obeys the Fokker–Planck equation

\pdv{t} p_t(x) = - \nabla \cdot \bigl( p_t(x) \, v(x, t) \bigr) + \tfrac{1}{2} g_t^2 \, \nabla^2 p_t(x) \,.\tag{6}

We will return to this stochastic form when we get to diffusion models. For now, consider the ODE case where $g_t \equiv 0$ .

Parametrizing $v_\theta$ by a neural network gives the continuous normalizing flow (CNF; Chen et al., 2018): sample by solving (2), get $\log q_\theta$ by integrating (4) alongside, train by maximum likelihood. Let’s consider the computational cost involved. For a generic network the exact divergence $\nabla \cdot v_\theta$ costs about $D$ backward passes per ODE step in $D$ dimensions; stochastic trace estimators cut this to $O(1)$ at the price of noise (FFJORD).¹ The divergence cost is also architecture-dependent: special constructions make it cheap by design, and the lattice-field-theory flows of track 2b compute it exactly at roughly one forward pass cost. That effort is worth it there because per-sample exact likelihood — a capability ML practice often happily gives up — is what makes lattice results exact. But the durable obstruction is not the divergence: it is that every likelihood or loss evaluation requires solving the ODE across all of $[0, 1]$ . The simulation is what the fixed-path trick of the next section removes.

Continue in Notebook 1 §1, where the CNF code is provided; we run it on the spiral target and query both capabilities.

Fix the path, learn the field

Now we come to the conceptual change that has unlocked the scaling of generative models to high resolution and perceptual quality. Step back and ask what any transport method must provide: a chain of distributions whose ends are a tractable base and (approximately) the target, and sample transformations that realize the chain efficiently. The CNF learns the whole chain, and pays for it with simulation inside the loss. But while constructing a path toward an unknown target is hard, it is easy to define one in the opposite direction, by progressively noising the data. So we fix a known path of distributions connecting the two ends in advance, and only learn the field that follows it. Concretely, choose a noising interpolation between a data sample $x$ and a noise sample $z \sim \mathcal{N}(0, I)$ :

x_t = \alpha_t \, x + \sigma_t \, z ,\tag{7}

with $\alpha_0 = 0, \sigma_0 = 1$ (pure noise) and $\alpha_1 = 1, \sigma_1 \approx 0$ (data). Read (7) conditionally on the data sample: given $x$ , the noisy point is exactly Gaussian,

x_t \mid x \; \sim \; \mathcal{N}\bigl( \alpha_t x, \; \sigma_t^2 I \bigr) .\tag{8}

Whatever the schedule, (8) is the marginal of a simple linear SDE.² Importantly, sampling any intermediate distribution needs no simulation at all: Draw $x$ from the dataset, draw $z$ , form $x_t$ . Averaging (8) over the dataset defines the intermediate marginals

p_t(x_t) = \int p_{\text{data}}(x) \, \mathcal{N}\bigl( x_t \mid \alpha_t x, \sigma_t^2 I \bigr) \, dx ,\tag{9}

the density of $x_t$ when we forget which data sample it came from. Read in the noising direction, this family is a solution of the Fokker–Planck equation (6) for a suitable linear drift and noise rate. Whatever object we decide to learn, the network can be trained on single $(x_t, t)$ points, one “time slice” at a time; the loss has become simulation-free.

Convention: why t = 0 is noise

Throughout, $t$ runs from $0$ (pure noise) to $1$ (data), matching the direction of generation. Many texts run it the other way, so their $t$ is our $1 - t$ . Depending on the perspective, placing noise at $t = 0$ can look unnatural. The noising process is an Ornstein–Uhlenbeck flow:

\dd{x} = -\, x \, \dd\tau + \sqrt{2} \, \dd{W}_\tau ,\tag{10}

which on its own time $\tau \in [0, \infty)$ forgets the data and approaches its stationary distribution $\mathcal{N}(0, I)$ as $\tau \to \infty$ . Any finite parametrization truncates that half-line. Our choice above compactifies it. For the variance-preserving path, $\alpha = e^{-\tau}$ maps $\tau \in [0, \infty)$ onto $t \in (0, 1]$ , landing the stationary end exactly at $t = 0$ , and (7) is the closed-form solution of (10) read through that map. Our generation time is more natural in the sense that the generative ODE is the only equation we integrate: the noising side is solved in closed-form via (8).

t · 0 noise → 1 data1.00

αₜ = 1.00 · σₜ = 0.00

Figure 3. The noising process: the two-armed swirl target (t = 1, right end of the slider) is joined to a standard Gaussian (t = 0) under the variance-preserving trig schedule αₜ = sin(πt/2), σₜ = cos(πt/2). Left: one fixed realization of the noising process — each point is kicked by fresh noise at every step (the exact Gaussian transition, no discretization), so the paths are rough, not straight lines. Right: the intermediate marginal pₜ. The marginals the cloud sweeps out are exactly the ones the fixed path xₜ = αₜ x + σₜ z samples in closed form with no simulation — draw a data point, draw noise, mix.

What should the network learn on this path? The central object is the score (or Stein score) $s(x, t) := \nabla_x \log p_t(x)$ , the gradient field of the log of the marginal density (9) (Figure 4). It determines everything about reversing the noising. In the next section we convert it directly into the velocity that transports noise to data. Thus, the natural loss is a direct regression onto it:

\mathcal{L}_{\text{SM}}(\theta) = \mathbb{E}_{t, \, x_t \sim p_t} \Bigl[ \, \bigl\| s_\theta(x_t, t) - \nabla_{x_t} \log p_t(x_t) \bigr\|^2 \, \Bigr] .\tag{11}

But the regression target is exactly what we cannot evaluate: the marginal (9) integrates over the whole dataset. What we can evaluate is its conditional counterpart. Fix the data sample $x$ : then $x_t$ is Gaussian by (8), and differentiating its log-density gives the conditional score in closed form,

\nabla_{x_t} \log p_t(x_t \mid x) = - \frac{x_t - \alpha_t x}{\sigma_t^2} = - \frac{z}{\sigma_t} .\tag{12}

Denoising score matching is the replacement of the intractable marginal target in (11) by this tractable conditional one:³

\mathcal{L}_{\text{DSM}}(\theta) = \mathbb{E}_{t, x, z} \Bigl[ \, \bigl\| s_\theta(x_t, t) + z / \sigma_t \bigr\|^2 \, \Bigr] .\tag{13}

The box below shows this leaves the objective untouched: (11) and (13) provably differ by a $\theta$ -independent constant, so they have the same gradient and the same minimizer (the marginal score). Still, $-z/\sigma_t$ equals the marginal score only in conditional mean. Once we use a finite number of samples to estimate the gradient (in SGD), we inherit the spread of $z$ , and the estimate becomes increasingly noisy near clean data ( $t=1$ ). In practice the network is often parametrized as a noise predictor, $\varepsilon_\theta(x_t, t) := -\sigma_t \, s_\theta(x_t, t)$ , which turns (13) into predicting the $O(1)$ target $z$ from $x_t$ by mean squared error.

Exercise score, noise, and denoised image are one target

The loss (13) can be read as learning the score, or the added noise $z$ , or the denoised sample $x$ . Show these three targets are equivalent, and give the map between them.

Hint

$x_t = \alpha_t x + \sigma_t z$ is one linear relation between $x_t$ , $x$ , and $z$ ; fix $x_t$ and any one of $\{s, z, x\}$ determines the rest.

t · 0 noise → 1 data0.50

αₜ = 0.71 · σₜ = 0.71

Figure 4. The central object of Part I as a picture: the exact score ∇log pₜ points up the log-density ridges of the intermediate marginal (faint underlay), and grows as σₜ shrinks — arrow lengths are saturated for display, directions are exact. Switch to the learned field: a network trained only on noisy points and their noise — never on a marginal quantity — carries the marginal score sθ = −εθ/σₜ, up to errors where it saw little data (between the arms, in the far tails).

Exercise why the learned score looks smoother at the edges

Toggle Figure 4 to the learned field. Near the domain edges it is smoother and more regular than the exact score. Why?

Hint

Two ingredients: what data did the network see out there, and why does low density not suppress the magnitude of $\nabla \log p_t$ ?

Why does regressing on a single noise draw work?

We regress $s_\theta$ on the conditional score (12), a target $-z/\sigma_t$ built from one $(x, z)$ pair — yet we want the marginal score, an average over all the data behind $x_t$ . Why does one draw recover the average? Two facts combine.

Squared-error regression returns the conditional mean. Let $\bar s := \mathbb{E}[-z/\sigma_t \mid x_t, t]$ be the target’s mean at fixed $(x_t, t)$ , and add and subtract it inside the loss:

\mathbb{E} \, \bigl\| s_\theta(x_t, t) + z/\sigma_t \bigr\|^2 = \underbrace{\mathbb{E} \, \bigl\| s_\theta(x_t, t) - \bar s \bigr\|^2}_{\theta\text{-dependent}} + \underbrace{\mathbb{E} \, \bigl\| {-z/\sigma_t} - \bar s \bigr\|^2}_{\text{constant}} .\tag{14}

The cross term vanishes: given $(x_t, t)$ the first factor is fixed and the second, $-z/\sigma_t - \bar s$ , averages to zero. Only the first term depends on $\theta$ , so the minimizer is $s_\theta(x_t, t) = \bar s$ . The single-draw noise shifts the loss by a constant, not the optimum!

The conditional mean is the marginal score. This is a matter of some rewriting and using Bayes’ theorem:⁴

\begin{align} \nabla \log p_t(x_t) &= \frac{\nabla p_t(x_t)}{p_t(x_t)} = \frac{\int p_{\text{data}}(x) \, \nabla p_t(x_t \mid x) \, \dd{x}}{p_t(x_t)} \nonumber \\ &= \mathbb{E}_{x \mid x_t}\!\bigl[ \nabla \log p_t(x_t \mid x) \bigr] = \bar s \vphantom{\Huge|}. \end{align}

So $s_\theta$ trained on single noisy pairs lands on the marginal score, $s_\theta = \nabla \log p_t$ . The same two steps justify conditional flow matching in track 2a (Vincent, 2011).

Exercise the score is the optimal denoiser (Tweedie)

The best guess of the clean sample given a noised one is the posterior mean $\hat x = \mathbb{E}[x \mid x_t]$ . Show it is fixed entirely by the score,

\hat x = \frac{1}{\alpha_t} \bigl( x_t + \sigma_t^2 \, \nabla \log p_t(x_t) \bigr) ,

so denoising and score estimation are the same problem (Tweedie’s formula).

Hint

Average the conditional score (12) over the posterior, exactly as in the box above: $\nabla \log p_t(x_t) = \mathbb{E}_{x \mid x_t}[ -(x_t - \alpha_t x)/\sigma_t^2 ]$ , then solve for $\mathbb{E}[x \mid x_t]$ .

Continue in Notebook 1 §§2–3, where we implement the noising and the loss.

Diffusion models

Diffusion models are the canonical instance of the previous section: the Gaussian path (7) plus the noise-prediction regression. Our canonical schedule is the trigonometric choice $\alpha_t = \sin(\pi t / 2)$ , $\sigma_t = \cos(\pi t / 2)$ , an instance of what is called a variance-preserving⁵ schedule because $\alpha_t^2 + \sigma_t^2 = 1$ . The training objective is the regression of the previous section with a schedule-dependent weighting $w(t)$ :

\mathcal{L}(\theta) = \mathbb{E}_{x, z, t} \Bigl[ \, w(t) \, \bigl\| \varepsilon_\theta(x_t, t) - z \bigr\|^2 \, \Bigr] ,\tag{15}

which for suitable $w(t)$ recovers the DDPM loss (Ho et al., 2020); we take $w(t) = 1$ .

Sampling now means transporting noise back along the same path but in the opposite direction, towards the data distribution. Differentiating the conditional path gives the conditional velocity $\dot{x}_t = \dot\alpha_t x + \dot\sigma_t z$ . Taking its conditional expectation given $x_t$ (the same trick as above) yields a marginal velocity field expressed entirely through the trained network:

v_\theta(x, t) = \frac{\dot\alpha_t}{\alpha_t} \, x + \Bigl( \dot\sigma_t - \sigma_t \, \frac{\dot\alpha_t}{\alpha_t} \Bigr) \, \varepsilon_\theta(x, t) .\tag{16}

The PF-ODE is not the only sampler. Injected noise and transport along the score are exchangeable, so the PF-ODE is the zero-noise member of a one-parameter family of stochastic samplers that share the same marginals $p_t$ (Song et al., 2021); all reach the same distribution along different trajectories (Figure 5).

One network, many samplers

Why are noise and score-transport exchangeable? Injected noise diffuses samples down the density; adding a score drift $\lambda_t \, s_\theta$ transports them back up at the matching rate, so the ensemble $p_t$ is untouched — only individual trajectories jitter.⁶ Injecting noise at rate $2\lambda_t$ with that compensating drift turns the PF-ODE into a one-parameter family (Karras et al., 2022),

\dd{x} = \bigl( v_\theta(x, t) + \lambda_t \, s_\theta(x, t) \bigr) \dd{t} + \sqrt{2 \lambda_t} \, \dd{W}_t , \qquad \lambda_t \ge 0 ,\tag{17}

where $\lambda_t = 0$ gives the PF-ODE and $\lambda_t$ at half the noising rate gives Song et al.’s reverse-time SDE — Notebook 1’s stochastic sampler. With the exact score every member shares the marginals $p_t$ ; with a learned $s_\theta$ (and finite steps) they drift apart, and a little noise can even mask score errors — one reason stochastic samplers sometimes win.

t · 0 noise → 1 data1.00

αₜ = 1.00 · σₜ = 0.00

Figure 5. One trained network, two samplers: the stochastic SDE (left) and the deterministic PF-ODE (right) share the same intermediate marginals and reach the same destination distribution along very different trajectories — the PF-ODE is the zero-noise member of the family. Both panels integrate the same 14 fixed noise draws with 96 uniform steps; the slider moves along the trajectories, each shown as a head dot with a short fading tail. The underlay is the intermediate marginal pₜ — identical in both panels at every t — so both clouds visibly ride the same evolving density; at t = 1 (the default) it is the target itself.

The conceptual point of Part I is that, given data, we can empirically define a path of distributions from noise to data such that the object we must learn has a known, simulation-free loss. Everything downstream — flow matching, distillation — are variations of which object that is.

Continue in Notebook 1 §§4–5, where we train the spiral model and write the PF-ODE Euler sampler.

Part II The tracks

Track 2a — flow matching

Diffusion learned the score and converted it to a velocity at sampling time via (16); flow matching skips the detour and regresses the velocity field $v_\theta(x_t, t)$ directly (Lipman et al., 2023). The simplest path is the linear interpolant $x_t = t \, x + (1 - t) \, z$ , whose conditional velocity is constant along each path: $u = x - z$ . Why regressing this conditional velocity finds the marginal field is exactly the argument in the box of Fix the path, learn the field: the conditional-target MSE differs from the intractable marginal-target MSE by a $\theta$ -independent constant. The two methods are members of one family:

	path $(\alpha_t, \sigma_t)$	network	regression target
diffusion (Notebook 1)	$\sin(\pi t/2), \; \cos(\pi t/2)$	$\varepsilon_\theta(x_t, t)$	noise $z$
flow matching (2a)	$t, \; 1 - t$	$v_\theta(x_t, t)$	velocity $x - z$

While much of the appeal of flow matching lies in finding simpler “straighter” paths, we must be precise in what this means. The conditional paths are straight lines by construction. The marginal ODE trajectories are only straighter than the diffusion parametrization’s, not necessarily straight. Multimodal target forces curvature, and methods have been developed to actually straighten them (Liu et al., 2022). The important point is that straighter paths tolerate coarser integration: fewer Euler steps at equal quality, and thus cheaper inference.

Euler steps per samplek = 8

Figure 6. The same fixed batch of 400 noise draws, integrated with k Euler steps of each learned ODE (NFE: number of network forward evaluations per generated sample). The points move as k changes — they are never re-randomized — and the exact target density sits faintly behind them. The flow-matching model stays recognizable down to very few steps; the diffusion samples fall apart earlier. The gap mixes two causes: the linear path gives straighter marginal trajectories, and the ε-parametrization’s PF-ODE coefficients blow up like 1/αₜ near the noise end, which a coarse Euler grid resolves badly.

Continue in Notebook 2a, where we train the flow-matching model and run the step-count study.

Track 2b — flows for lattice field theory

Here, we turn to an application in science. In lattice field theory the target density is known up to the normalization constant (called partition function $Z$ in this context, the estimation of which is a goal in itself): $p(\phi) \propto e^{-S(\phi)}$ for an action $S$ we can write down. But there are no (or few) samples: producing them is the task, as physics prediction derive from expectation values $\langle O \rangle$ under $p$ . In addition (in contrast to some other common ML applications), the goal is faithful ensemble statistics and not individually impressive samples. With no dataset, the simulation-free losses of Part I do not apply, and the CNF’s exact density — the expensive capability Part I side-steps — becomes central. Training uses the reverse KL divergence: minimize $\mathbb{E}_{\phi \sim q_\theta}[\log q_\theta(\phi) + S(\phi)]$ , which needs only model samples and the action. The unknown normalization of $p$ shifts the loss by a constant and drops out of the gradient. Exact model density gets us another imporant feature, namely exactness of the physics predictions: importance reweighting corrects the residual model error, giving unbiased results in the infinite-sample limit. This can be interpreted in the context of classical Monte Carlo as learning a proposal distribution, where historically it was often painstakingly handcrafted.

Reverse KL is mode-seeking, and a flow that silently drops one mode of a symmetric target still reweights to perfect-looking values for every symmetry-even observable. We reproduce this failure mode in the notebook, catch it with symmetry-odd observables, and fixes it in the architecture. The appendix carries the same transport idea onto the SU(3) group manifold as a first step toward lattice gauge theory.

Continue in Notebook 2b, where we train a $\phi^4$ flow by reverse KL, reweight observables, and break the flow on purpose.

Track 2c — few-step distillation

Sampling in Notebook 1 costs on the order of a hundred network evaluations per sample (NFE: the number of neural-network forward evaluations per generated sample), and inference cost is the primary bottleneck at deployment. The key is that the map $z \mapsto x$ defined by integrating the PF-ODE (16) is deterministic. It is a single-valued function, not a distribution. Thus, we can try to learn this function: run a trained diffusion model teacher offline to generate $(z, x)$ pairs, then train a student $f_\theta(z)$ to jump straight to the endpoint. The notebook trains exactly this one-step student, dissects its characteristic artifacts (the teacher’s map is nearly discontinuous where neighboring $z$ land on different spiral arms, and a smooth student must interpolate across the gap), and explores various modern families of approaches. The step-count figure (6) above shows the exact quality axis distillation moves along. Network calls are not the only cost axis. For example, quantization and caching aim to reduce the price of each call. See the TinyML course for that direction.

Continue in Notebook 2c, where we distill the provided teacher into a one-step student and study its failure modes.

Stepping back

One idea carried through the whole story: a generative model transports probability mass from a simple base distribution to the generated data. The different method we encountered are different ways to realize that transport, and make different trade-offs: how much to fix in advance and how much to learn, and trading exact likelihood against the cost of training and sampling.

Several places to go from here: the tutorial repository with the notebooks, the MIT flow-matching and diffusion course notes are a good student-level treatment of exactly this material, and the Flow Matching Guide and Code goes deeper on track 2a.

The Hutchinson estimator replaces $\nabla \cdot v = \operatorname{tr}(\partial_x v)$ by $\mathbb{E}_\eta[\eta^\top (\partial_x v) \, \eta]$ with random probe vectors $\eta$ — one Jacobian-vector product per probe instead of $D$ . One caveat matters for track 2b: the estimate is unbiased for $\log q_\theta$ but not for $q_\theta$ itself — exponentiating a noisy estimate biases it upward (Jensen) — so stochastic divergences are disqualified where the density feeds an exactness argument, on correctness grounds and not just noise. ↩
Explicitly $\dd{x_t} = \frac{\dot\alpha_t}{\alpha_t} \, x_t \, \dd{t} + \sqrt{2 \sigma_t \, ( \tfrac{\dot\alpha_t}{\alpha_t} \sigma_t - \dot\sigma_t )} \; \dd{W}_t$ , run in the noising direction (decreasing $t$ ); the drift is just the mean law $\mathbb{E}[x_t \mid x] = \alpha_t x$ . ↩
Denoising is not the only escape from the intractable target: integration by parts removes it entirely at the cost of a Jacobian term (implicit score matching, Hyvärinen, 2005), and random projections make that term scalable (sliced score matching, Song et al., 2019). The denoising route is the one many modern diffusion models use, with different choices of time-dependent weighting. ↩
Two things easy to miss. First, the gradient is with respect to $x_t$ , not $x$ and thus $p_{\text{data}}(x)$ is like a constant. Second, the move to the expectation simultaneously applies Bayes’ theorem and the log-derivative identity $\nabla p_t(x_t \mid x) = p_t(x_t \mid x) \, \nabla \log p_t(x_t \mid x)$ . ↩
“Variance-preserving” is a slight misnomer: the name says that the path ends at (and, run as a process, converges to) a unit Gaussian, distinguishing it from variance-exploding parametrizations — it does not preserve the variance of your data distribution unless that variance is already one. ↩
One line of Fokker–Planck. With $g_t^2 = 2\lambda_t$ the diffusion term of (6) is $\lambda_t \nabla^2 p_t$ , and since $\nabla p = p \, \nabla\!\log p$ it is a transport in disguise, $\lambda_t \nabla^2 p_t = \nabla \cdot \bigl( p_t \, \lambda_t \, \nabla \log p_t \bigr)$ ; this exactly offsets the compensating $\lambda_t \, s_\theta$ drift, leaving effective velocity $v_\theta$ for every $\lambda_t \ge 0$ . Run it the other way — trade all of the noising process’s diffusion for score transport — and the same identity re-derives (16). ↩

How do we generate samples?

Generative modeling as transport

Continuous transport

Fix the path, learn the field

Diffusion models

Track 2a — flow matching

Track 2b — flows for lattice field theory

Track 2c — few-step distillation

Stepping back

Footnotes