How to read these notes↩︎

The intended audience is beginning graduate students who have seen probability, linear algebra, and multivariable calculus, but may not have prior exposure to stochastic differential equations or diffusion models. The notes therefore separate three levels of material.

• Core definitions and identities are proved in detail. These run through every section and include the Fokker–Planck equation and entropy dissipation, Tweedie’s identity and the reverse SDE, the probability flow ODE, the exact and Gaussian DDPM reverse kernels, the denoising-score equivalence, the finite-state reverse kernel of discrete diffusion, and the reward-tilting identity.

• Representative estimates are proved under simplifying assumptions, to show a mechanism in its cleanest form. For example, we prove KL contraction under a log-Sobolev inequality, compute the ULA bias exactly for a Gaussian target, and bound the one-step discretization and score-error contributions to the diffusion-model KL.

• Research-level theorems are stated in simplified form with a proof roadmap. These include the ULA and MALA convergence guarantees, high-accuracy diffusion sampling via first-order rejection sampling and the discrete-diffusion error analysis. The goal is not to reproduce every technical lemma of the original papers, but to make clear what each theorem is saying, why the assumptions appear, and how the proof is organized.

Two appendices collect the analytic tools used repeatedly: Itô calculus and the Girsanov change-of-measure formula in Appendix 9, and the Gaussian identities in Appendix 10. A reader comfortable with these may skip them and refer back as needed.

These notes are self-contained but deliberately selective, and they overlap with several excellent treatments that a reader may wish to consult in parallel. Yuansi Chen’s lecture notes for the ETH course Computational and Statistical Aspects of Diffusion Models [1] cover much of the same ground with a stronger emphasis on convergence proofs, guidance, and discrete diffusion, and are a useful companion to the later sections here. For the Langevin and log-concave sampling side, Chewi’s book [2] is the standard reference. For stochastic calculus and numerical analysis of SDEs and CTMCs, we refer to the book [3] by E, Li, and Vanden-Eijnden. Throughout, when a result first appeared in a research paper we cite the original; the lecture notes above are flagged here once, as a unified entry point to several of the topics that follow.

These notes do not try to survey the full modern generative-modeling landscape. We focus on diffusion and Langevin mechanisms for which reverse-time SDEs, probability-flow ODEs, score-estimation identities, and inference-time control are the central actors. The notes treat learned-score error as an input to sampler guarantees, but do not discuss statistical learning theory for estimating scores from finite data. We also do not give detailed treatments of likelihood-based normalizing flows, adversarial and variational models, autoregressive architectures, large-scale latent-diffusion engineering, or several newer continuous-time alternatives. For entry points to the latter, see flow matching [4], rectified flows [5], stochastic interpolants [6], consistency models [7], and recent one-step formulations such as mean flow [8]. These methods share much of the same language of probability paths, transport equations, denoisers, and ODE/SDE samplers.

Minimal probability background↩︎

We use the following conventions repeatedly. If \(X\) is a random variable, then \(\mathop{\mathrm{Law}}(X)\) denotes its distribution. If \(X\) has density \(p\), expectations are written either as \(\mathbb{E}[f(X)]\) or \(\int f(x)p(x)\,\mathrm{d}x\). A conditional expectation \(\mathbb{E}[Y\mid X=x]\) is best thought of as the best prediction of \(Y\) as a function of the observed value \(x\). For a Markov chain, \(P(x,\,\mathrm{d}y)\) denotes the distribution of the next state given the current state \(x\). For an SDE, all formal differentiations can be justified under smoothness and decay assumptions; the notes use this formal calculus as a way to keep the main ideas visible.

We often use the same symbol for a probability law and its density when the meaning is clear from context. For example, \(p_{\mathsf{data}}\) may denote the data law, the density of that law, or the measure \(p_{\mathsf{data}}(\,\mathrm{d}x)\). Similarly, \(p_t\) may denote either \(\mathop{\mathrm{Law}}(X_t)\) or its density. This common abuse of notation keeps the formulas readable, but when a density value is needed we write expressions such as \(\nabla\log p_t(x)\).

Recurring notation↩︎

The following symbols are used across several sections. More local notation is introduced where it is needed.

Algorithm acronyms↩︎

The following acronyms name algorithms or algorithmic discretizations used in several sections.

AI usage disclosure↩︎

Large-language-model tools were used in the preparation of these notes to help with drafting, editing, reorganization, and consistency checks. The mathematical content, exposition choices, references, and any remaining errors are the responsibility of the author.

1 Introduction to Sampling and Langevin Dynamics↩︎

Sampling is the problem of producing random points that look as if they were drawn from a given distribution. This is easy to state but often hard to do: the distribution may be known only up to a normalizing constant, or only through data, and drawing from it directly may be infeasible. The strategy running through these notes is indirect. Rather than sample the hard target in one shot, we build a random process that is easy to simulate and whose distribution gradually drifts toward the target, and we read off a sample once the process has run long enough. Making this idea precise calls for one habit above all: instead of following a single trajectory, track the whole distribution as it moves through the algorithm. This first section deliberately goes slowly, fixing that habit and the language we use to measure how far one distribution is from another.

1.1 The sampling problem↩︎

Let \(\pi\) be a probability distribution on \(\mathbb{R}^d\). In most applications, \(\pi\) is not available through exact samples. Instead we may have access to:

• the unnormalized density \(\pi(x)\propto e^{-U(x)}\);

• the gradient \(\nabla\log\pi(x)=-\nabla U(x)\);

• samples from a noisy distribution related to \(\pi\);

The goal is to construct a random variable \(\widehat{X}\) with law \(\widehat{\pi}\) such that \(D(\widehat{\pi},\pi)\leq \varepsilon\), where \(D\) is an appropriate discrepancy. There is no single best choice of \(D\). Let \(\mu\) and \(\nu\) be two probability laws on \(\mathbb{R}^d\). We recall three basic discrepancies.

1. Total variation asks whether every event has nearly the right probability: \[\mathop{\mathrm{D_{\mathsf{TV}}}}(\mu,\nu)=\sup_A\left\lvert \mu(A)-\nu(A) \right\rvert =\frac{1}{2}\int\left\lvert \frac{\,\mathrm{d}\mu}{\,\mathrm{d}\lambda} -\frac{\,\mathrm{d}\nu}{\,\mathrm{d}\lambda} \right\rvert\,\mathrm{d}\lambda.\] Here \(\lambda\) is any measure such that both \(\mu\) and \(\nu\) are absolutely continuous with respect to \(\lambda\), written \(\mu,\nu\ll\lambda\); for example, we may take \(\lambda=\mu+\nu\). The notation \(\,\mathrm{d}\mu/\,\mathrm{d}\lambda\) denotes the Radon–Nikodym derivative, the density of \(\mu\) with respect to \(\lambda\), characterized by \[\mu(A)=\int_A \frac{\,\mathrm{d}\mu}{\,\mathrm{d}\lambda}\,\,\mathrm{d}\lambda\] for measurable sets \(A\). The same applies to \(\,\mathrm{d}\nu/\,\mathrm{d}\lambda\), and the value of the total-variation distance does not depend on which such reference measure \(\lambda\) is chosen.

2. Wasserstein-\(2\) distance asks whether probability mass can be transported a short geometric distance: \[W_2^2(\mu,\nu)=\inf_{\gamma\in\Pi(\mu,\nu)} \int \left\lVert x-y \right\rVert^2\,\gamma(\,\mathrm{d}x,\,\mathrm{d}y).\] Here \(\Pi(\mu,\nu)\) is the set of couplings of \(\mu\) and \(\nu\), meaning probability measures \(\gamma\) on \(\mathbb{R}^d\times\mathbb{R}^d\) whose first marginal is \(\mu\) and whose second marginal is \(\nu\).

3. KL divergence asks whether \(\mu\) can be encoded efficiently using \(\nu\) as a reference: \[\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu) = \int \log\!\left(\frac{\,\mathrm{d}\mu}{\,\mathrm{d}\nu}\right)\,\mathrm{d}\mu,\] when \(\mu\) is absolutely continuous with respect to \(\nu\), and \(+\infty\) otherwise. Note that the KL divergence is asymmetric: \(\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu) \neq \mathop{\mathrm{D_{\mathsf{KL}}}}(\nu\|\mu)\). In particular, KL is not a metric.

KL is often the main bookkeeping divergence we will use in this course. The other metrics translate this information into statements about events, geometric displacement, or weak test functions. The basic comparison begins with the Csiszár–Kullback–Pinsker inequality; see, for example, Bakry, Gentil, and Ledoux [9]: \[\mathop{\mathrm{D_{\mathsf{TV}}}}(\mu,\nu) \leq \sqrt{\frac{1}{2}\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu)}.\] Thus a KL guarantee immediately gives a TV guarantee, and hence also a weak test-function guarantee. Another recurring rule is data processing: applying the same observation map to two random objects cannot increase KL. If \(T\) is measurable and \(T_\#\mu\) denotes the law of \(T(X)\) for \(X\sim\mu\), then \[\mathop{\mathrm{D_{\mathsf{KL}}}}(T_\#\mu\|T_\#\nu)\leq \mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu).\] The same monotonicity holds for total variation. We use this principle repeatedly to pass from path-space comparisons to endpoint laws; a proof is given in Lemma 2.

Wasserstein distance is more geometric. If \(\mu\) and \(\nu\) are supported on a set of diameter \(R\), then a maximal coupling gives \[W_2^2(\mu,\nu) \leq R^2\mathop{\mathrm{D_{\mathsf{TV}}}}(\mu,\nu) \leq R^2\sqrt{\frac{1}{2}\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu)}.\] Without a bounded-diameter or moment/functional-inequality assumption, there is no universal comparison in the other direction: small \(W_2\) need not imply small TV or KL, and small TV or KL need not control \(W_2\) if a tiny amount of mass can move very far away.

A more powerful comparison is available relative to a fixed reference law. A probability law \(\nu\) satisfies a Talagrand \(T_2\) transport-entropy inequality with constant \(C_T\) if \[W_2^2(\mu,\nu)\leq 2C_T\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu) \qquad\text{for all }\mu\ll\nu.\] One standard route to this inequality is through log-Sobolev. In the convention of these notes, if \(\nu\) satisfies \[\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu) \leq \frac{C_{\mathsf{LSI}}}{2}\mathop{\mathrm{FI}}(\mu\|\nu),\] where \[\mathop{\mathrm{FI}}(\mu\|\nu) = \int \left\lVert \nabla\log\frac{\,\mathrm{d}\mu}{\,\mathrm{d}\nu} \right\rVert^2\,\mathrm{d}\mu\] when the derivative is well defined, then the Otto–Villani theorem gives the transport bound \[W_2^2(\mu,\nu)\leq 2C_{\mathsf{LSI}}\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu).\] This is one reason log-Sobolev and transport inequalities appear naturally in sampling theory; see Villani [10] and Bakry, Gentil, and Ledoux [9].

1.2 Markov chain Monte Carlo↩︎

Markov chain Monte Carlo (MCMC) turns sampling into the problem of running a Markov chain. Instead of drawing directly from \(\pi\), we choose a transition rule whose repeated application should have \(\pi\) as its equilibrium law. A Markov kernel \(P\) maps a current state \(x\) to a distribution \(P(x,\cdot)\). If \(\pi P=\pi\), then \(\pi\) is invariant, and the ideal MCMC chain runs \[X_{n+1}\sim P(X_n,\cdot),\qquad n=0,1,\ldots,\] with the hope that \(\mathop{\mathrm{Law}}(X_n)\) approaches \(\pi\) as \(n\to\infty\).

This viewpoint immediately separates two questions. First, is the transition kernel designed so that the right law is invariant? Second, how long must the chain run before it is close to equilibrium? For a genuine metric \(D\), such as TV or \(W_2\), a typical algorithmic decomposition is \[D(\mathop{\mathrm{Law}}(X_n),\pi) \leq \underbrace{D(\mathop{\mathrm{Law}}(X_n),\pi_h)}_{\text{mixing to the algorithm's invariant law}} + \underbrace{D(\pi_h,\pi)}_{\text{bias from approximation}},\] where \(\pi_h\) is the invariant distribution of the implementable algorithm. This is a triangle-inequality argument, so it applies to genuine metrics but not to KL. In general there is no inequality of the form \(\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\rho)\leq\mathop{\mathrm{D_{\mathsf{KL}}}}(\mu\|\nu)+\mathop{\mathrm{D_{\mathsf{KL}}}}(\nu\|\rho)\).

This decomposition separates the error caused by not running the algorithm long enough from the error caused by using a transition rule whose invariant law is not exactly \(\pi\). For Langevin sampling, that second error appears when a continuous diffusion is replaced by a time-stepping rule. The discretization may lead to an implementable chain with \(\pi_h\neq\pi\). In exact MCMC methods, such as Metropolis-adjusted Langevin, the transition is constructed so that \(\pi_h=\pi\).

1.2.0.1 Caution: invariance is not convergence.

Invariance is only a fixed-point statement: if \(X_0\sim\pi\), then one step of the chain still has law \(\pi\). It does not by itself say that a chain started from another law will approach \(\pi\). For example, the identity kernel \(P(x,\cdot)=\delta_x\) leaves every distribution invariant, but it never mixes. Even uniqueness of the invariant law is not quite enough without excluding periodic behavior: on the two-point space \(\{0,1\}\), the deterministic flip \(0\mapsto1\), \(1\mapsto0\) has the uniform distribution as an invariant reversible law, but a chain started from \(0\) oscillates forever.

Thus detailed balance is a convenient way to verify invariance, not a convergence theorem. To justify MCMC, one also needs an ergodicity mechanism, such as irreducibility and aperiodicity in finite state spaces, or the appropriate Harris recurrence and minorization conditions in general state spaces. Quantitative sampling bounds then add still more structure: spectral gaps, conductance, contraction, or functional inequalities. See Meyn and Tweedie [11] for a classical treatment of these Markov-chain conditions.

A common design principle for MCMC is therefore to start one level above the implementable chain. First design a continuous-time Markov process whose invariant law is \(\pi\) and whose convergence mechanism can be analyzed. Then discretize that process to obtain a computable Markov chain, possibly adding a Metropolis correction if one wants to remove discretization bias. Langevin dynamics is the canonical example of this principle: the continuous diffusion has \(\pi\) as its invariant law, while ULA and MALA are two different ways to turn that diffusion into an algorithm.

1.3 Overdamped Langevin diffusion↩︎

Assume \(\pi(x)\propto e^{-U(x)}\) with \(U:\mathbb{R}^d\to\mathbb{R}\) smooth. The overdamped Langevin diffusion is the continuous-time Markov process \[\label{eq:langevin} \,\mathrm{d}X_t = -\nabla U(X_t)\,\mathrm{d}t+\sqrt{2}\,\,\mathrm{d}B_t = \nabla\log\pi(X_t)\,\mathrm{d}t+\sqrt{2}\,\,\mathrm{d}B_t.\tag{1}\] The drift \(-\nabla U=\nabla\log\pi\) points toward regions where \(\pi\) is larger, much as gradient descent moves toward lower potential. The diffusion coefficient in the Brownian term is tuned so that the ensemble keeps precisely the spread prescribed by \(\pi\).

This is the basic difference between optimization and sampling. Optimization asks for a point \(x_\star\) with small objective value, whereas sampling asks for a law. If \(\pi\) is bimodal, with half of its mass near each of two separated modes, an optimizer may correctly return one mode, but a sampler must visit both modes with the correct frequencies.

There are two complementary ways to read 1 . Pathwise, each particle is pulled down the potential landscape by \(-\nabla U\) while Brownian motion keeps the ensemble spread out. Distributionally, the density \(q_t\) of \(X_t\) evolves by a deterministic PDE. The corresponding Fokker–Planck equation is \[\label{eq:fokker-planck} \partial_t q_t = \nabla\cdot(q_t\nabla U)+\Delta q_t = \nabla\cdot\!\left(q_t\nabla\log\frac{q_t}{\pi}\right).\tag{2}\] At this stage the point of 2 is simply that the stochastic particle system has a deterministic law-level description. The SDE picture is good for intuition and simulation. The PDE picture is good for proving invariance and entropy decay.

To derive the Fokker–Planck equation, we recall Itô calculus in the form needed for our calculation, more details can be found in Appendix 9. If \[\,\mathrm{d}X_t=b(X_t)\,\mathrm{d}t+\sigma\,\,\mathrm{d}B_t\] with constant diffusion matrix \(\sigma\), then for a smooth test function \(\varphi\), \[\,\mathrm{d}\varphi(X_t) = \left\langle \nabla\varphi(X_t),b(X_t) \right\rangle\,\mathrm{d}t +\frac{1}{2}\mathop{\mathrm{Tr}}\!\left(\sigma\sigma^\top\nabla^2\varphi(X_t)\right)\,\mathrm{d}t +\left\langle \nabla\varphi(X_t),\sigma\,\,\mathrm{d}B_t \right\rangle.\] The last term is a martingale increment, so it disappears after taking expectations. For 1 , \(b=-\nabla U\) and \(\sigma=\sqrt2\,I\). Applying the preceding formula to 1 gives \[\frac{\,\mathrm{d}}{\,\mathrm{d}t}\mathbb{E}[\varphi(X_t)] = \mathbb{E}[-\left\langle \nabla U(X_t),\nabla\varphi(X_t) \right\rangle+\Delta\varphi(X_t)].\] If \(X_t\) has density \(q_t\), then \[\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int \varphi q_t\,\mathrm{d}x = \int\left[-\left\langle \nabla U,\nabla\varphi \right\rangle+\Delta\varphi\right]q_t\,\mathrm{d}x.\] Integrating by parts, assuming boundary terms vanish, \[\int -q_t\left\langle \nabla U,\nabla\varphi \right\rangle\,\mathrm{d}x = \int \varphi\,\nabla\cdot(q_t\nabla U)\,\mathrm{d}x, \qquad \int q_t\Delta\varphi\,\mathrm{d}x = \int \varphi\,\Delta q_t\,\mathrm{d}x.\] Since this holds for all test functions \(\varphi\), we obtain \[\partial_tq_t=\nabla\cdot(q_t\nabla U)+\Delta q_t.\] Finally, \(\pi\propto e^{-U}\) implies \(\nabla U=-\nabla\log\pi\), and \[\nabla\cdot(q_t\nabla U)+\Delta q_t = \nabla\cdot\left(q_t\nabla\log\frac{q_t}{\pi}\right).\]

The differential operator that appeared in the test-function calculation, \[\mathcal{L}\varphi = -\left\langle \nabla U,\nabla\varphi \right\rangle+\Delta\varphi,\] is called the infinitesimal generator of the Langevin diffusion. The Fokker–Planck equation is the adjoint equation on densities. In weak form, the calculation above says \[\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int \varphi q_t\,\mathrm{d}x = \int (\mathcal{L}\varphi) q_t\,\mathrm{d}x = \int \varphi(\mathcal{L}^\ast q_t)\,\mathrm{d}x,\] so \(\partial_t q_t=\mathcal{L}^\ast q_t\), where \[\mathcal{L}^\ast q = \nabla\cdot(q\nabla U)+\Delta q.\]

1.4 Entropy dissipation↩︎

Differentiating KL along the flow gives \[\label{eq:entropy-dissipation} \frac{\,\mathrm{d}}{\,\mathrm{d}t}\mathop{\mathrm{D_{\mathsf{KL}}}}(q_t\|\pi) = -\int \left\lVert \nabla\log\frac{q_t}{\pi} \right\rVert^2 q_t\,\,\mathrm{d}x =-\mathop{\mathrm{FI}}(q_t\|\pi),\tag{3}\] where \(\mathop{\mathrm{FI}}(q\|\pi)\) is the relative Fisher information.

The identity just proved becomes a quantitative convergence estimate once it is combined with a functional inequality. If \(\pi\) satisfies a log-Sobolev inequality \[\mathop{\mathrm{D_{\mathsf{KL}}}}(q\|\pi)\leq \frac{C_{\mathsf{LSI}}}{2}\mathop{\mathrm{FI}}(q\|\pi),\] then 3 implies \[\frac{\,\mathrm{d}}{\,\mathrm{d}t}\mathop{\mathrm{D_{\mathsf{KL}}}}(q_t\|\pi) \leq -\frac{2}{C_{\mathsf{LSI}}}\mathop{\mathrm{D_{\mathsf{KL}}}}(q_t\|\pi).\] By Grönwall’s inequality, \[\mathop{\mathrm{D_{\mathsf{KL}}}}(q_t\|\pi) \leq e^{-2t/C_{\mathsf{LSI}}}\mathop{\mathrm{D_{\mathsf{KL}}}}(q_0\|\pi).\] Thus the continuous-time process reduces KL error exponentially.

In particular, if the potential \(U\) is \(m\)-strongly convex, that is \(\nabla^2U\succeq mI\) with \(m>0\), then by the Bakry–Émery criterion [9] the target satisfies the log-Sobolev inequality with \(C_{\mathsf{LSI}}\leq 1/m\), and the estimate above becomes the clean exponential rate \[\mathop{\mathrm{D_{\mathsf{KL}}}}(q_t\|\pi)\leq e^{-2mt}\mathop{\mathrm{D_{\mathsf{KL}}}}(q_0\|\pi).\]

The main lesson is the mechanism: KL divergence is a Lyapunov function, the relative Fisher information is its dissipation rate, and a structural inequality for the target converts dissipation into convergence. Sampling arguments often follow this pattern: choose a discrepancy, compute its derivative along the dynamics, and use an inequality for the target law to turn the derivative into a rate.

The same identity also has a geometric reading. In optimal-transport language, overdamped Langevin is the \(W_2\)-gradient flow of \(\mathop{\mathrm{D_{\mathsf{KL}}}}(\cdot\|\pi)\), and 3 is its energy-dissipation identity; see Ambrosio, Gigli, and Savaré [12]. For background on diffusion semigroups, entropy dissipation, and log-Sobolev inequalities, see Bakry, Gentil, and Ledoux [9]; for the transport viewpoint behind Wasserstein geometry, see Villani [10].

2 Convergence of Langevin Diffusion and ULA↩︎

The previous section established continuous-time convergence through the entropy identity and log-Sobolev inequality. We now ask what remains of that story after replacing the diffusion by an implementable Markov chain. The main new issue is numerical error: a time-stepping rule may converge as a Markov chain, but its invariant law need not be the target law. For a modern treatment of log-concave sampling algorithms, including Langevin algorithms and their Metropolis-adjusted variants, see the excellent book by Chewi [2].

2.1 Unadjusted Langevin algorithm↩︎

Euler–Maruyama with step size \(h>0\) gives the unadjusted Langevin algorithm (ULA) \[\label{eq:ula} X_{n+1}=X_n-h\nabla U(X_n)+\sqrt{2h}\,\xi_n, \qquad \xi_n\sim\mathcal{N}\!\left(0,I\right).\tag{4}\] ULA is often the first sampler one considers because it is simple, gradient-based, and cheap. Its limitation is equally important: 4 is not an exact transition of 1 , and in general does not preserve \(\pi\).

It is tempting to view ULA as “almost Langevin” and therefore assume that it must have almost the right stationary law. Already for a Gaussian target one can see exactly what “almost” means. For small step size the bias is small, but it does not vanish at fixed \(h\) no matter how long the chain is run. This is the difference between mixing error, which decreases with more iterations, and discretization bias, which is built into the transition rule itself.

2.2 One-step KL discretization error↩︎

Let \(P_h(x,\cdot)\) be the exact Langevin transition for time \(h\), and let \(\widehat{P}_h(x,\cdot)\) be the ULA transition. To compare them in KL, interpolate one ULA step over the interval \([0,h]\) by \[\,\mathrm{d}\widehat{X}_s=-\nabla U(x)\,\mathrm{d}s+\sqrt2\,\,\mathrm{d}B_s,\qquad \widehat{X}_0=x.\] The exact Langevin law started from the same point has drift field \(y\mapsto-\nabla U(y)\). Since the KL below is \(\mathop{\mathrm{Law}}(\widehat{X}_{[0,h]})\) relative to \(\mathop{\mathrm{Law}}(X_{[0,h]})\), the Girsanov formula evaluates both drift fields along the first path, namely at \(\widehat{X}_s\). Thus the formula reviewed in Appendix 9 gives the path-space comparison \[\mathop{\mathrm{D_{\mathsf{KL}}}}\!\left(\mathop{\mathrm{Law}}(\widehat{X}_{[0,h]})\middle\|\mathop{\mathrm{Law}}(X_{[0,h]})\right) = \frac{1}{4} \int_0^h \mathbb{E}\left\lVert \nabla U(\widehat{X}_s)-\nabla U(x) \right\rVert^2\,\mathrm{d}s.\] By the data-processing inequality reviewed in Lemma 2, the KL between the endpoint laws is no larger. If \(\nabla U\) is \(L\)-Lipschitz, then \[\mathop{\mathrm{D_{\mathsf{KL}}}}\!\left(\widehat{P}_h(x,\cdot)\middle\|P_h(x,\cdot)\right) \leq \frac{L^2}{4} \int_0^h \mathbb{E}\left\lVert \widehat{X}_s-x \right\rVert^2\,\mathrm{d}s.\] Since \(\widehat{X}_s=x-s\nabla U(x)+\sqrt2\,B_s\), \[\mathbb{E}\left\lVert \widehat{X}_s-x \right\rVert^2=s^2\left\lVert \nabla U(x) \right\rVert^2+2ds,\] and therefore \[\label{eq:ula-local-kl} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left(\widehat{P}_h(x,\cdot)\middle\|P_h(x,\cdot)\right) \leq \frac{L^2}{4} \left(dh^2+\frac{h^3}{3}\left\lVert \nabla U(x) \right\rVert^2\right).\tag{5}\] This is the KL analogue of a local truncation error: freezing the drift for one short interval costs order \(dh^2\), plus a term depending on the local drift size.

2.3 Direct KL recursion for ULA↩︎

The one-step estimate explains the scale of the discretization error, but it is not by itself a convergence theorem to \(\pi\). Indeed, a comparison such as \(\mathop{\mathrm{D_{\mathsf{KL}}}}(\widehat{q}_k\|q_{kh})\) cannot simply be added to \(\mathop{\mathrm{D_{\mathsf{KL}}}}(q_{kh}\|\pi)\), because KL has no triangle inequality. The modern KL-based analysis instead tracks the target-relative quantity \(\mathop{\mathrm{D_{\mathsf{KL}}}}(\widehat{q}_k\|\pi)\) directly.

The following result is a translation of the KL theorem of Vempala and Wibisono [13] into our notation.

No convexity is assumed in this theorem. The log-Sobolev inequality supplies the global mixing mechanism, while the bounded Hessian assumption supplies the local control needed to discretize the diffusion.

Thus, to reach KL target \(\mathop{\mathrm{D_{\mathsf{KL}}}}(\widehat{q}_k\|\pi)\leq \varepsilon^2\), it suffices to choose \[h\lesssim \min\left\{ \frac{1}{C_{\mathsf{LSI}}L^2}, \frac{\varepsilon^2}{C_{\mathsf{LSI}}L^2d} \right\}\] and run for \[k \gtrsim \frac{C_{\mathsf{LSI}}}{h} \log\frac{\mathop{\mathrm{D_{\mathsf{KL}}}}(\widehat{q}_0\|\pi)}{\varepsilon^2}\] iterations. Equivalently, the iteration complexity is \[k=\widetilde{O}\!\left( \frac{C_{\mathsf{LSI}}^2L^2d}{\varepsilon^2} \right).\] The dependence on the ambient dimension is linear in this KL analysis. By the Bakry–Émery criterion for log-Sobolev inequalities (see Bakry, Gentil, and Ledoux [9]), under \(m\)-strong convexity one has \(C_{\mathsf{LSI}}\leq 1/m\), so the bound becomes \(\widetilde{O}(\kappa^2d/\varepsilon^2)\) complexity, where \(\kappa=L/m\) is the condition number.

Let us take a look at the result of the Theorem. The dependence on dimension is linear, but the dependence on accuracy is polynomial: to make the KL error of order \(\varepsilon^2\), ULA uses \(\widetilde{O}(C_{\mathsf{LSI}}^2L^2d/\varepsilon^2)\) steps. This happens because the Euler discretization has a stationary bias of order \(h\), so high accuracy forces a small step size. It is therefore natural to ask whether one can keep the Langevin proposal but remove the discretization bias by an exact correction. Classical Metropolis adjustment does precisely this when density ratios are available.

2.4 Metropolis-adjusted Langevin algorithm↩︎

The Metropolis-adjusted Langevin algorithm (MALA) corrects the ULA proposal with an accept-reject step. The distinction from the gradient-only setting is central: density evaluations enable the classical high-accuracy rejection ideas, whereas gradient information alone does not provide the density ratios those ideas require.

Let \(\pi(\,\mathrm{d}x)\propto e^{-U(x)}\,\mathrm{d}x\). From a current state \(x\), MALA first draws the ULA proposal \[y\sim q_h(x,\cdot) = \mathcal{N}\!\left(x-h\nabla U(x),2hI\right).\] It then accepts \(y\) with probability \[\label{eq:mala-acceptance} a_h(x,y) = 1\wedge \frac{\pi(y)q_h(y,x)}{\pi(x)q_h(x,y)}.\tag{7}\] If the proposal is rejected, the chain stays at \(x\). The Metropolis ratio 7 enforces detailed balance: \[\pi(\,\mathrm{d}x)P_h^{\mathsf{MALA}}(x,\,\mathrm{d}y) = \pi(\,\mathrm{d}y)P_h^{\mathsf{MALA}}(y,\,\mathrm{d}x).\] Consequently, \(\pi\) is exactly stationary for MALA. This is the basic advantage over ULA: the discretization bias is removed, at the price of needing the density ratio \(\pi(y)/\pi(x)\).

MALA is the cleanest example of a correction mechanism. The proposal uses only gradient information, just like ULA, but the accept-reject step asks whether the proposed move has the right probability under the target density. This single density-ratio check changes the invariant distribution from an approximation \(\pi_h\) back to the exact target \(\pi\).

Exact stationarity, however, is only the invariance part of the story. To turn MALA into a quantitative sampling algorithm, one still has to bound how many accepted-or-rejected proposals are needed before the chain is close to \(\pi\). The modern high-accuracy theory separates this analysis into two tasks: sampling efficiently once a good initialization is available, and producing that initialization in the first place. Both are phrased relative to a warm start. A standard warm-start condition is the following: an initial law \(\mu_0\) is \(M\)-warm with respect to \(\pi\) if \[\mu_0(A)\leq M\pi(A) \qquad\text{for every measurable set }A,\] or equivalently \[\label{eq:warmstartinfty} \left\|\frac{\,\mathrm{d}\mu_0}{\,\mathrm{d}\pi}\right\|_{L^\infty(\pi)} \leq M.\tag{8}\] Some warm-start preparation results use the finite-order Rényi version: \[D_q(\mu_0\|\pi) = \frac{1}{q-1} \log \int \left(\frac{\,\mathrm{d}\mu_0}{\,\mathrm{d}\pi}\right)^q \,\mathrm{d}\pi =O(1)\] for some \(q>1\). Note that in the Rényi notation, 8 can be written as \(D_\infty(\mu_0\|\pi)\leq \log M\).

For the first task, fix a target total-variation accuracy \(\varepsilon\in(0,1)\). Assuming an \(M\)-warm start, Wu, Schmidler, and Chen [14] prove that, for \(m\)-strongly log-concave and \(L\)-smooth targets in \(\mathbb{R}^d\), MALA mixes in \[O\!\left( \kappa\sqrt d\, \log^3\!\max\left\{\kappa,d,\frac{M}{\varepsilon}\right\} \right), \qquad \kappa=L/m,\] iterations, up to universal constants. Chen and Gatmiry [15] extend this warm-start MALA picture under smoothness and isoperimetry, recovering the same leading \(\kappa\sqrt d\) behavior in the strongly log-concave case, again with logarithmic dependence on \(M/\varepsilon\).

For the second task, Altschuler and Chewi [16] show that kinetic Langevin, also called underdamped Langevin, can produce the required finite-order Rényi warm start with the same \(\widetilde{O}(\sqrt d)\) dimension dependence. At the level of dimension dependence, the high-accuracy picture is: use kinetic Langevin to prepare a warm start, then use MALA as the exact Metropolis-corrected sampler.

The discussion above also shows the limitation of classical adjusted MCMC. MALA removes discretization bias when density evaluations are available, but its sharp convergence guarantees require structural assumptions such as strong log-concavity, log-Sobolev or isoperimetric inequalities, smoothness, and warm starts. A general data distribution may be multimodal, singular, or available only through samples, so these assumptions are not a natural starting point. This motivates a different question: if we are given data rather than an explicit target density, how should we design an algorithm to generate new samples?

3 Score-Based Diffusion Models↩︎

The previous section treated sampling from an explicit density: the algorithm could use \(\nabla U\), and MALA could even use density ratios to remove discretization bias. Diffusion models begin from a different premise. The target law is represented by data, not by a tractable formula, so we do not try to run Langevin dynamics directly on the data distribution. Instead we add noise to the data, forming a path of smoother laws. At each positive noise level, the score of the noised law is the local vector field that guides a reverse dynamics back toward the data.

The denoising-score connection underlying this approach goes back in part to the DDPM formulation of Ho, Jain, and Abbeel [17]. The continuous-time SDE and probability-flow formulation used here follows Yang Song, Sohl-Dickstein, Kingma, Kumar, Ermon, and Poole [18]. The basic picture is simple: at time \(t\), the noised distribution \(p_t\) is a blurred version of the data, and \(\nabla\log p_t(x)\) points toward nearby regions where this blurred density is larger. Diffusion models learn this time-indexed field of local denoising directions.

3.1 Continuous-time forward noising↩︎

The cleanest conceptual starting point is continuous time. Let \(X_0\sim p_{\mathsf{data}}\) and let \((X_t)_{t\in[0,T]}\) solve a forward noising SDE \[\label{eq:forward-sde} \,\mathrm{d}X_t=f_t(X_t)\,\mathrm{d}t+g_t\,\mathrm{d}B_t, \qquad X_0\sim p_{\mathsf{data}}.\tag{9}\] Here \(B_t\) is Brownian motion, \(f_t\) is a drift field, and \(g_t\) is a scalar diffusion coefficient. We write \(p_t\) for the density of \(X_t\). Most diffusion models use linear Gaussian forward processes for which, for deterministic functions \(a_t\) and \(\sigma_t\), \[\label{eq:continuous-marginal} X_t\mid X_0\sim\mathcal{N}\!\left(a_tX_0,\sigma_t^2I\right).\tag{10}\]

The two scalar functions \(a_t\) and \(\sigma_t\) summarize the signal-to-noise ratio. The coefficient \(a_t\) says how much of the original sample remains visible in the conditional mean, while \(\sigma_t\) says how much independent Gaussian uncertainty has been added. Early in the forward process, \(\sigma_t\) is small and the score may be complicated because it reflects the fine structure of the data. At large noise, \(p_t\) is smoother and closer to a simple reference law, so the reverse sampler has an easier starting point.

3.2 Continuous-time Tweedie identity↩︎

The most natural way to undo the forward noising is to denoise: given a noisy observation \(X_t=x\), estimate the clean sample that produced it, namely the posterior mean \(\mathbb{E}[X_0\mid X_t=x]\). The object that encodes this denoising is the score of the noised law, \[\label{eq:continuous-score} \mathsf{s}^{\star}_t(x)=\nabla\log p_t(x),\tag{11}\] the true score at time \(t\), approximated in practice by a learned model \(\mathsf{s}_t(x)\approx\mathsf{s}^{\star}_t(x)\). The Gaussian marginal formula 10 makes the link between the score and denoising precise through a continuous-time form of Tweedie’s identity: \[\label{eq:continuous-tweedie} \mathsf{s}^{\star}_t(x) = \frac{1}{\sigma_t^2} \mathbb{E}[a_tX_0-X_t\mid X_t=x].\tag{12}\] Equivalently, whenever \(a_t\neq0\), define the optimal denoiser \[\mathsf{D}^{\star}_t(x) := \mathbb{E}[X_0\mid X_t=x] = a_t^{-1}\bigl(x+\sigma_t^2\mathsf{s}^{\star}_t(x)\bigr).\] A learned score \(\mathsf{s}_t\) similarly determines a learned denoiser \[\mathsf{D}_t(x) := a_t^{-1}\bigl(x+\sigma_t^2\mathsf{s}_t(x)\bigr).\] Thus denoising and score estimation are two ways of describing the same posterior information.

This Tweedie identity explains why the score can be learned from empirical data. Fix \(t>0\). Under the Gaussian corruption \(X_t=a_tX_0+\sigma_t Z\), with \(Z\sim\mathcal{N}\!\left(0,I\right)\) independent of \(X_0\), the random vector inside this conditional expectation becomes \[\frac{a_tX_0-X_t}{\sigma_t^2} = -\frac{Z}{\sigma_t}.\] Thus the conditional mean of \(-Z/\sigma_t\) given \(a_tX_0+\sigma_t Z=x\) is exactly \(\mathsf{s}^{\star}_t(x)\), which motivates the regression population loss: \[\mathbb{E}\left\lVert \mathsf{s}_t(a_tX_0+\sigma_t Z)+\frac{Z}{\sigma_t} \right\rVert^2.\] Replacing the expectation over \(X_0\sim p_{\mathsf{data}}\) by an empirical average over data points, while resampling the Gaussian noise, gives the basic denoising score-matching objective.

This is why the learned field has a denoising interpretation rather than being an arbitrary vector field. The regression target is noisy for each individual corruption, but its conditional average points from the observation toward the posterior mean \(\mathsf{D}^{\star}_t(x)\) of the clean sample. Equivalently, the score records the Bayes correction encoded by the noised data distribution.

3.3 Continuous-time reverse SDE↩︎

The forward SDE 9 is designed to move data toward a simple law. For sampling, the object we need is not a pathwise inverse of the Brownian motion. It is enough to construct a Markov process whose one-time marginals run through the same densities in the opposite order: if \((p_t)_{0\leq t\leq T}\) are the forward noising marginals, then the reverse sampler \((Y^{\leftarrow}_s)_{0\leq s\leq T}\) should satisfy \[\mathop{\mathrm{Law}}(Y^{\leftarrow}_s)=p_{T-s} \qquad\text{for every }s\in[0,T].\] This requirement is deliberately distributional. It says that the reverse sampler has the right snapshots, not that it retraces individual forward Brownian paths. We therefore first match the evolution of densities. Set \(q_s=p_{T-s}\) and write \(t=T-s\). The forward density satisfies \[\partial_t p_t = -\nabla\cdot(f_tp_t) +\frac{1}{2}g_t^2\Delta p_t.\] Therefore, with \(t=T-s\), \[\partial_s q_s = \nabla\cdot(f_t p_t) -\frac{1}{2}g_t^2\Delta p_t.\] The key step uses the identity \(\Delta p_t=\nabla\cdot(p_t\nabla\log p_t)\), which brings in precisely the score \(\mathsf{s}^{\star}_t=\nabla\log p_t\) from 11 . This rewrites the density evolution as \[\partial_s q_s = -\nabla\cdot\left[ \left(-f_t+g_t^2\mathsf{s}^{\star}_t\right)p_t \right] +\frac{1}{2}g_t^2\Delta p_t.\] This is the Fokker–Planck equation for the reverse-time diffusion \[\label{eq:reverse-sde} \,\mathrm{d}Y^{\leftarrow}_s = \left[-f_{T-s}(Y^{\leftarrow}_s) +g_{T-s}^2\,\mathsf{s}^{\star}_{T-s}(Y^{\leftarrow}_s)\right]\,\mathrm{d}s +g_{T-s}\,\mathrm{d}B^{\leftarrow}_s, \qquad Y^{\leftarrow}_0\sim p_T.\tag{13}\] Here \(B^{\leftarrow}_s\) is Brownian motion in the reverse sampling time. By construction, if \(Y^{\leftarrow}_0\sim p_T\), then \(Y^{\leftarrow}_s\sim p_{T-s}\) for every \(s\in[0,T]\).

The drift in the reverse SDE has two pieces. The term \(-f_{T-s}\) reverses the deterministic transport part of the forward dynamics. The extra term \(g_{T-s}^2\,\mathsf{s}^{\star}_{T-s}\) is exactly the score introduced in Subsection 3.2 as the denoising direction: it corrects the density evolution caused by diffusion by nudging each sample back toward regions of higher noised density, i.e.toward the posterior mean of the clean sample. This is the precise sense in which the denoiser of the previous subsection is the engine of the reverse dynamics. Replacing the true score \(\mathsf{s}^{\star}_t\) by a learned score \(\mathsf{s}_t\) yields the score-based reverse SDE sampler, the continuous object that DDPM-type algorithms discretize.

It is instructive to specialize 13 to the two schedules of the running example; we will revisit the same two processes from the probability-flow viewpoint in the next subsection.

In practice the reverse SDE is run in discrete time, and this discretization is exactly DDPM. One fixes a grid in the sampling time \(s\), replaces the true score \(\mathsf{s}^{\star}_{T-s}\) by the learned model \(\mathsf{s}_{T-s}\), and takes Euler–Maruyama steps of 13 : each step nudges the current sample along the denoising drift and then adds fresh Gaussian noise of the appropriate variance. This is the DDPM sampler of Ho, Jain, and Abbeel [17], the stochastic counterpart of the deterministic DDIM update we will read off from the probability-flow ODE in the next subsection. The detailed Gaussian transition kernels and the resulting discretization error are taken up in Section 5.

3.4 Probability flow ODE↩︎

There is a deterministic ODE whose one-time marginals are the same as those of the forward SDE 9 . This ODE is called the probability flow ODE. The probability flow ODE replaces Brownian randomness by a transport field that has the same effect on marginal densities. Define the velocity field \[\label{eq:pf-velocity} v_t(x) = f_t(x)-\frac{1}{2} g_t^2\mathsf{s}^{\star}_t(x) = f_t(x)-\frac{1}{2} g_t^2\nabla\log p_t(x).\tag{14}\] The probability flow ODE is \[\label{eq:pf-ode} \frac{\,\mathrm{d}X_t}{\,\mathrm{d}t}=v_t(X_t).\tag{15}\] If \(X_0\sim p_0\), then the solution of 15 has marginal density \(p_t\) at every time \(t\) for which the ODE is well posed.

The reverse SDE and the probability flow ODE are two different machines that produce the same snapshots when the score is exact. The SDE machine keeps injecting randomness; the ODE machine deterministically transports particles (starting from random initial data). Since many numerical and control questions depend on paths, not only snapshots, one should not freely replace one machine by the other without checking what quantity is being analyzed.

For sampling, we use the same reverse-time convention as in 13 : set \(s=T-t\), start from a high-noise draw with law \(p_T\), and integrate from \(s=0\) to \(s=T\). The reverse-time probability-flow sampler is the ODE \[\frac{\,\mathrm{d}Z^{\leftarrow}_s}{\,\mathrm{d}s} = -v_{T-s}(Z^{\leftarrow}_s) = -f_{T-s}(Z^{\leftarrow}_s) +\frac{1}{2} g_{T-s}^2\mathsf{s}^{\star}_{T-s}(Z^{\leftarrow}_s), \qquad Z^{\leftarrow}_0\sim p_T.\] By the same continuity-equation calculation, \(Z^{\leftarrow}_s\) has law \(p_{T-s}\) when the score is exact. Compare this with the reverse SDE drift in 13 : \(-f_{T-s}+g_{T-s}^2\mathsf{s}^{\star}_{T-s}\). The DDIM sampler of Song, Meng, and Ermon [19] is the standard sampler on this probability-flow side: it uses the same trained denoising model as DDPM, but follows the deterministic update suggested by the ODE viewpoint.

The deterministic nature of the ODE has an important consequence for likelihoods. Along a solution of \(\dot{X}_t=v_t(X_t)\), the continuity equation 17 implies \[\label{eq:instantaneous-change} \frac{\,\mathrm{d}}{\,\mathrm{d}t}\log p_t(X_t) = -\nabla\cdot v_t(X_t).\tag{18}\] Indeed, \[\frac{\,\mathrm{d}}{\,\mathrm{d}t}\log p_t(X_t) = \partial_t\log p_t(X_t)+ \left\langle v_t(X_t),\nabla\log p_t(X_t) \right\rangle,\] and dividing 17 by \(p_t\) gives \[\partial_t\log p_t+\left\langle v_t,\nabla\log p_t \right\rangle = -\nabla\cdot v_t.\] Therefore, if the map \(X_0\mapsto X_T\) is obtained by integrating the probability flow ODE forward, \[\log p_0(X_0) = \log p_T(X_T) + \int_0^T \nabla\cdot v_t(X_t)\,\mathrm{d}t.\] This is the continuous normalizing-flow identity used to compute likelihoods from probability flow ODEs.

4 Stochastic Localization and Polchinski Flow↩︎

The previous section introduced diffusion models through Gaussian noising. We started from data, followed the noised marginals \(p_t\), learned the score \(\nabla\log p_t\), and used that score to run either a reverse SDE or a probability-flow ODE. This section looks at the same Gaussian channel from two complementary angles: stochastic localization and Polchinski flow. All three viewpoints start from a clean random variable and the noisy observation (here and below we take the variance-exploding one) \[X_t=X_\star+\sqrt{t}Z,\] but they organize the information in different ways. Diffusion models emphasize reverse-time sampling. Stochastic localization asks a Bayesian question: as the noisy observation becomes more precise, how does the posterior law of the hidden signal \(X_\star\) evolve? Polchinski flow asks a density-level question: as the noise level changes, how do the smoothed density \(p_t\) and the effective potential \(U_t=-\log p_t\) evolve?

These two changes of perspective are useful because they expose structure that is less visible from the reverse SDE alone. In stochastic localization, the score appears as a posterior denoising correction. In Polchinski flow, the score is the negative gradient of a coarse-grained energy landscape. Both viewpoints lead naturally to quantitative tools, especially covariance identities, that will be useful later.

We begin by treating stochastic localization as an observation model for the same Gaussian channel. This lets us compare diffusion time with localization precision, derive the posterior law, and connect its mean to the score. We then record the martingale structure that makes localization useful. Finally, we return to the density \(p_t\) itself and write the corresponding Polchinski flow for the effective potential \(U_t\).

4.1 Stochastic localization↩︎

Let \(X_\star\sim\mu\) be the signal, or clean data point. Here we write the law as \(\mu\) rather than \(p_{\mathsf{data}}\), to emphasize that the construction applies to a general distribution and not only to the data law. Start with the Gaussian smoothing channel \[X_t=X_\star+\sqrt{t}Z,\qquad Z\sim\mathcal{N}\!\left(0,I\right).\] Here \(t\) is the noise variance, and \(p_t=\mathop{\mathrm{Law}}(X_t)\) is the Gaussian convolution of the data law with covariance \(tI\). As \(t\) increases, the density becomes smoother; as \(t\) decreases, the channel approaches the original data distribution.

Stochastic localization uses the same Gaussian channel, but parametrizes it by precision rather than variance. Observe the same unknown \(X_\star\) through the continuous-time Gaussian observation process \[\label{eq:sl-observation} \,\mathrm{d}Y_u=X_\star\,\mathrm{d}u+\,\mathrm{d}W_u,\qquad Y_0=0.\tag{19}\] Here \(W_u\) is Brownian motion independent of \(X_\star\). Conditional on \(X_\star=x\), we have \(Y_u\sim\mathcal{N}\!\left(ux,uI\right)\), and thus after a rescaling: \[\bar Y_u:=u^{-1}Y_u\sim\mathcal{N}\!\left(x,u^{-1}I\right).\] Hence the localization time \(u\) corresponds to the diffusion noise variance \(t=u^{-1}\): \[\bar Y_u=X_\star+\sqrt{t}Z, \qquad t=\frac{1}{u}.\] The same noised observation can therefore be indexed either by variance \(t\) or by precision \(u\). Diffusion notation emphasizes smoothing as \(t\) increases. Localization notation emphasizes Bayesian inference as \(u\) increases and the observation becomes more informative. This precision parametrization is useful because the posterior mean and covariance then obey simple martingale identities, as we will discuss below.

We now write this posterior law explicitly. At precision \(u\), the observation is \(Y_u\), and the basic object is the conditional law of the hidden signal \(X_\star\). Let \(\mu_u(\cdot \mid y)=\mathop{\mathrm{Law}}(X_\star\mid Y_u = y)\) be the posterior law after observing \(Y_u = y\). Bayes’ rule gives \[\label{eq:sl-posterior} \mu_{u}(\,\mathrm{d}x \mid Y_u = y) = \frac{1}{Z_u(y)} \exp\left\{\left\langle y,x \right\rangle-\frac{u}{2}\left\lVert x \right\rVert^2\right\}\mu(\,\mathrm{d}x).\tag{20}\] The normalizing constant is \[Z_u(y)= \int \exp\left\{\left\langle y,x \right\rangle-\frac{u}{2}\left\lVert x \right\rVert^2\right\}\mu(\,\mathrm{d}x).\] This is why the method is called localization: the posterior is the original measure tilted by a random linear field and penalized by a growing quadratic function. As \(u\) grows, the posterior becomes increasingly concentrated near the hidden signal.

Now put the random observation back in and write \[\mu_u=\mu_u(\cdot\mid Y_u), \qquad M_u=\int x\,\mu_u(\,\mathrm{d}x).\] The original stochastic-localization viewpoint is that the random measure \((\mu_u)_{u\geq0}\) itself evolves by a measure-valued SDE, which we now derive. Let \(\mathcal{F}_u^Y=\sigma(Y_r:0\leq r\leq u)\) and define the innovation process \[I_u=Y_u-\int_0^u M_s\,\mathrm{d}s, \qquad \,\mathrm{d}I_u=\,\mathrm{d}Y_u-M_u\,\mathrm{d}u.\] This subtracts the part of the next observation that is already predictable from the current posterior. Intuitively, \(I_u\) records only the “surprise” in the observation stream: after conditioning on \(\mathcal{F}_u^Y\), the drift of \(\,\mathrm{d}Y_u\) is \(M_u\,\mathrm{d}u\), so subtracting it leaves an increment with no predictable component and the same accumulated covariance as the original Brownian noise. More formally, one can check that \(I_u\) is a continuous \(\mathcal{F}_u^Y\)-martingale whose accumulated covariance over \([0,u]\) is \(uI\); Levy’s characterization then identifies it as a Brownian motion with respect to \(\mathcal{F}_u^Y\). The posterior law satisfies \[\label{eq:sl-measure-sde} \,\mathrm{d}\mu_u(\,\mathrm{d}x) = \left\langle x-M_u,\,\mathrm{d}I_u \right\rangle\,\mu_u(\,\mathrm{d}x).\tag{21}\] Equivalently, for every bounded test function \(f\), \[\,\mathrm{d}\int f(x)\,\mu_u(\,\mathrm{d}x) = \left\langle \int f(x)(x-M_u)\,\mu_u(\,\mathrm{d}x),\,\mathrm{d}I_u \right\rangle.\] The same identity applies componentwise to vector-valued test functions when the relevant moments are finite. In particular, for every Borel set \(A\), the process \(u\mapsto\mu_u(A)=\Pr(X_\star\in A\mid\mathcal{F}_u^Y)\) is a martingale: \[\mathbb{E}[\mu_v(A)\mid\mathcal{F}_u^Y]=\mu_u(A), \qquad v\geq u.\] Thus, in this formulation, stochastic localization is a measure-valued martingale whose random tilt concentrates the prior while preserving conditional expectations.

4.2 Posterior means and martingale structure↩︎

The previous subsection constructed the posterior law and its measure-valued martingale equation. We now pass to its first two moments. First we freeze the observation and connect the posterior mean to the diffusion-model score; then we put the random observation back in and record the martingale dynamics.

For a deterministic observation \(y\), define the posterior mean and covariance \[m_u(y)=\mathbb{E}[X_\star\mid Y_u=y], \qquad \Sigma_u(y)=\mathop{\mathrm{Cov}}(X_\star\mid Y_u=y).\] The posterior mean is exactly a denoiser. Indeed, since \(\bar Y_u=Y_u/u=X_\star+u^{-1/2}Z\), a deterministic observation \(y\) corresponds to the noised sample \(y/u\), and Tweedie’s identity 12 gives \[\label{eq:sl-tweedie} m_u(y) = \frac{y}{u} + \frac{1}{u} \nabla\log p_{1/u}\!\left(\frac{y}{u}\right),\tag{22}\] where \(p_t\) is the density of \(X_\star+\sqrt{t}Z\). Equivalently, \[\nabla\log p_{1/u}(y/u) = u\bigl(m_u(y)-y/u\bigr).\] Thus learning a score is the same problem as learning the localization posterior mean.

The covariance, in turn, is the linear response of this denoiser: differentiating the posterior tilt 20 with respect to the observation gives \[\label{eq:posterior-mean-jacobian} \nabla_y m_u(y)=\Sigma_u(y).\tag{23}\] Combining this with Tweedie’s identity 22 , and writing \(t=1/u\), gives the score-Jacobian identity \[\label{eq:score-hessian-covariance} \nabla \mathsf{s}^{\star}_t(x) = t^{-2}\Sigma_{1/t}(x/t) - t^{-1}I.\tag{24}\] Up to this time change, then, the Hessian of the noised log-density is precisely the posterior covariance of the localization process.

Now return to the random localization path. We write \(M_u\) and \(\Sigma_u\) below for the random evaluations \(m_u(Y_u)\) and \(\Sigma_u(Y_u)\). The measure-valued equation 21 then turns the fixed-observation posterior moments into stochastic dynamics.

Start with the mean. Applying the weak form of 21 componentwise to the test function \(f(x)=x\) gives, under standard integrability assumptions, \[\label{eq:sl-mean-sde} \,\mathrm{d}M_u = \left(\int x(x-M_u)^\top \mu_u(\,\mathrm{d}x)\right)\,\mathrm{d}I_u = \Sigma_u\,\,\mathrm{d}I_u.\tag{25}\] Thus the mean moves only in response to the innovation Brownian motion, and the size and direction of that motion are governed by the current posterior covariance. In particular, the posterior mean is a martingale: before seeing the next small piece of data, the current posterior mean is the best prediction of the next posterior mean.

The covariance, which set the size of these mean fluctuations, has an evolution of its own. Unlike the mean, it carries a drift, and that drift is strictly dissipative: in one dimension, \[\frac{\,\mathrm{d}}{\,\mathrm{d}u}\mathbb{E}[\Sigma_u]=-\mathbb{E}[\Sigma_u^2]\leq0,\] and in multiple dimensions the trace satisfies \[\label{eq:localization-covariance-dissipation} \frac{\,\mathrm{d}}{\,\mathrm{d}u}\mathbb{E}[\mathop{\mathrm{Tr}}\Sigma_u] = -\mathbb{E}[\mathop{\mathrm{Tr}}(\Sigma_u^2)]\leq0.\tag{26}\] Integrate 26 over \(u\in[0,\infty)\): the fundamental theorem of calculus gives \[\int_0^\infty\mathbb{E}\mathop{\mathrm{Tr}}(\Sigma_u^2)\,\,\mathrm{d}u = \mathbb{E}\mathop{\mathrm{Tr}}\Sigma_0-\lim_{u\to\infty}\mathbb{E}\mathop{\mathrm{Tr}}\Sigma_u \leq \mathbb{E}\mathop{\mathrm{Tr}}\Sigma_0,\] the inequality because \(\mathop{\mathrm{Tr}}\Sigma_u\geq0\). At \(u=0\) no observation has yet been made, so \(\Sigma_0=\mathop{\mathrm{Cov}}(X_\star)\); hence if \(\mathop{\mathrm{Cov}}(X_\star)\preceq I\), \[\label{eq:localization-hessian-budget} \int_0^\infty \mathbb{E}\mathop{\mathrm{Tr}}(\Sigma_u^2)\,\,\mathrm{d}u \leq \mathbb{E}\mathop{\mathrm{Tr}}\Sigma_0 \leq d,\tag{27}\] so the total squared covariance accumulated along the entire path is bounded by the dimension. This is the geometric content of localization: the random posterior measure steadily concentrates, while its mean drifts by martingale increments whose size is set by the uncertainty that remains.

This martingale-and-covariance structure is what makes stochastic localization useful as an analytical tool. Rather than bounding the original high-dimensional measure directly, one embeds it in a random path of progressively tilted posteriors and argues along the path, where the mean is a martingale and the covariance dissipates. Many static questions thereby reduce to estimating quantities averaged over the localization path, with the covariance identities 26 –27 supplying the bookkeeping.

The main early applications of stochastic localization [20] were in high-dimensional geometry. For an isotropic log-concave measure, the Kannan–Lovász–Simonovits (KLS) conjecture asks, in one formulation, for a dimension-free lower bound on the Cheeger isoperimetric constant. This problem is closely tied to spectral gaps, Poincaré inequalities, thin-shell estimates, and concentration of measure. Eldan introduced stochastic localization to relate thin-shell and spectral-gap/KLS bounds up to logarithmic factors [20]. Lee and Vempala developed the method further for isoperimetry, concentration, and mixing, improving the known KLS-type bounds for isotropic log-concave measures [21]. Later work of Chen obtained an almost constant lower bound for the KLS isoperimetric coefficient using related localization techniques [22]. Chen and Eldan also made the sampling connection explicit: their localization-schemes framework associates Markov chains to localized martingales of measures and proves mixing bounds by analyzing the localization process [23]. For our purposes, the lesson is that stochastic localization converts global geometric or sampling questions into estimates on the random posterior covariance process. The same covariance-budget viewpoint is what later reappears in analyses of diffusion models [24], [25].

4.3 Polchinski flow for the effective potential↩︎

The localization picture describes the Gaussian channel through posterior quantities. Polchinski flow describes the same channel at the level of the noised marginals. In diffusion-model language, this is the variance-exploding forward process from the previous section: \[X_t=X_0+\sqrt{t}Z,\] whose marginal density is the Gaussian-smoothed density \[p_t=p_0*\mathcal{N}\!\left(0,tI\right).\] It solves the heat equation \[\label{eq:heat-flow-density} \partial_t p_t=\frac{1}{2}\Delta p_t.\tag{29}\] Define the effective potential \[U_t(x)=-\log p_t(x).\] The diffusion-model score at noise level \(t\) is therefore \(\mathsf{s}^{\star}_t=\nabla\log p_t=-\nabla U_t\). Thus Polchinski flow is not a new sampling process; it is a way to track how the same score field evolves as the forward diffusion smooths the data distribution. The heat equation becomes the nonlinear PDE \[\label{eq:polchinski-finite} \partial_t U_t = \frac{1}{2}\Delta U_t - \frac{1}{2}\left\lVert \nabla U_t \right\rVert^2.\tag{30}\] This is the finite-dimensional Polchinski equation. It is the isotropic Gaussian-convolution version of the Polchinski renormalization flow [26], [27]. The density formulation is linear, but the effective potential formulation is nonlinear because logarithms turn smoothing into a viscous Hamilton–Jacobi equation.

Because the score is the negative gradient of the effective potential, \[\mathsf{s}^{\star}_t(x)=\nabla\log p_t(x)=-\nabla U_t(x),\] the score satisfies the viscous Burgers equation \[\label{eq:burgers-score} \partial_t \mathsf{s}^{\star}_t = \frac{1}{2}\Delta \mathsf{s}^{\star}_t + \nabla\left(\frac{1}{2}\left\lVert \mathsf{s}^{\star}_t \right\rVert^2\right).\tag{31}\] Diffusion models are usually trained by score matching, rather than solving the Burgers equation. The equation nevertheless explains why regularity of the score improves at positive noise levels and why small-noise estimates are delicate: as \(t\downarrow0\), the smoothing that controls derivatives disappears.

The Polchinski equation 30 is useful because it shows what forward noising does to the energy landscape. Sharp wells, ridges, and nonsmooth features of \(U_0=-\log p_0\) are averaged out by convolution. At positive noise, \(U_t\) is smoother and its gradient is better behaved. The reverse model is trying to follow this family of smoothed landscapes backward, where the geometry becomes progressively sharper.

This is the renormalization interpretation of diffusion models. Historically, renormalization entered statistical physics through the idea that one should describe a system at a chosen resolution and then track how the effective description changes when microscopic degrees of freedom are averaged out. Kadanoff’s block-spin scaling picture [28] gave the physical image of coarse-graining near criticality. Wilson’s renormalization group made this image into a systematic flow of effective theories and explained universality at critical points [29], [30]; the review of Wilson and Kogut [31] is the classical reference. Polchinski’s exact renormalization group equation [26] then expressed this flow directly at the level of effective Lagrangians. The finite-dimensional equation 30 is the diffusion-model shadow of the same idea.

Renormalization group language says that forward noising integrates out microscopic information. In image models this statement is metaphorical but useful: the data distribution is gradually blurred, high-frequency details become harder to distinguish, and the effective potential \(U_t\) describes the log-density of the coarse-grained law. In field-theoretic applications the same idea can be literal: one chooses a covariance schedule that removes degrees of freedom by scale, and the corresponding effective action follows a Polchinski-type equation. In the work of Bauerschmidt and collaborators, this flow is used precisely in this rigorous sense: it is an evolution equation for scale-dependent effective potentials, so estimates along the flow turn renormalization into analytic control of stochastic dynamics and Gibbs measures [27].

4.4 A dictionary↩︎

The same Gaussian smoothing operation can now be summarized in three languages: diffusion dynamics, posterior localization, and renormalization by Gaussian coarse-graining.

5 Discretizing Continuous Diffusion Models↩︎

We now pass from ideal continuous dynamics to algorithms. A sampler must choose a finite time grid and replace each exact reverse transition by a computable transition kernel. This section sets up the formulas needed for the numerical analysis.

A useful way to keep the presentation organized is to separate the three implementation questions. First, the reverse chain must be initialized from a simple high-noise law rather than the exact terminal law. Second, the reverse dynamics use learned scores or denoisers rather than the true ones. Third, the continuous reverse dynamics, or equivalently the exact Bayes reverse kernels on a grid, must be discretized. The analysis is built around exactly this trichotomy: initialization error, statistical score error, and discretization error.

5.1 The DDPM grid, scores, and exact reverse kernels↩︎

Let \(X_0\sim p_{\mathsf{data}}\) on \(\mathbb{R}^d\). Choose a time grid \[0=t_0<t_1<\cdots<t_K=T.\] The goal is to replace the continuous noising process by a finite Markov chain with the same prescribed one-time marginals. Recall from 10 that these marginals are \[X_t\mid X_0\sim\mathcal{N}\!\left(a_tX_0,\sigma_t^2I\right).\] On the grid, we write \[X_k = X_{t_k}, \qquad a_k=a_{t_k}, \qquad \sigma_k=\sigma_{t_k}.\] The grid schedule determines the one-step noising parameters. We choose a Gaussian forward transition \[\label{eq:one-step-noising} X_{k+1}\mid X_k \sim \mathcal{N}\!\left(\alpha_k X_k,\alpha_k^2\eta_kI\right).\tag{32}\] If \[X_k=a_kX_0+\sigma_k Z_k,\qquad Z_k\sim\mathcal{N}\!\left(0,I\right),\] then 32 is equivalently \[X_{k+1}=\alpha_kX_k+\alpha_k\sqrt{\eta_k}\xi_k, \qquad \xi_k\sim\mathcal{N}\!\left(0,I\right)\] and hence \[X_{k+1}\mid X_0 \sim \mathcal{N}\!\left( \alpha_k a_kX_0 , \alpha_k^2(\sigma_k^2+\eta_k)I \right).\] Matching this with the target marginal \(X_{k+1}\mid X_0\sim\mathcal{N}\!\left(a_{k+1}X_0,\sigma_{k+1}^2I\right)\) gives \[\alpha_k=\frac{a_{k+1}}{a_k}, \qquad \eta_k=\frac{\sigma_{k+1}^2}{\alpha_k^2}-\sigma_k^2.\] The variance-preserving case has \(a_k^2+\sigma_k^2=1\). The variance-exploding case has \(a_k\equiv1\) and an increasing noise scale \(\sigma_k^2=t_k\).

The reverse sampler has to invert one grid step: given a noisy point at level \(k+1\), it needs the conditional law of the previous point at level \(k\). This is where denoising and score information enter. Let \(p_k\) denote the density of \(X_k\). The true score at time \(k\) is \[\mathsf{s}^{\star}_k(x)=\nabla\log p_k(x).\]

The score determines the posterior mean of the clean sample. This is the grid-time specialization of the continuous-time Tweedie identity 12 , with \(a_t,\sigma_t\) read off on the grid: \[\label{eq:tweedie} \mathsf{s}^{\star}_k(x) = \frac{1}{\sigma_k^2} \mathbb{E}[a_kX_0-X_k\mid X_k=x].\tag{33}\] Equivalently, the grid-time version of the optimal denoiser from Subsection 3.2 is \[\mathsf{D}^{\star}_k(x) := \mathbb{E}[X_0\mid X_k=x] = a_k^{-1}\bigl(x+\sigma_k^2\mathsf{s}^{\star}_k(x)\bigr).\] This is the grid-time notation corresponding to the posterior mean \(m_u(y)\) from stochastic localization: in the variance-exploding normalization \(\sigma_t^2=t\), one has \(\mathsf{D}^{\star}_t(x)=x+t\mathsf{s}^{\star}_t(x)=m_{1/t}(x/t)\). A score model \(\mathsf{s}_k\) is equivalently a denoiser model \[\mathsf{D}_k(x) := a_k^{-1}\bigl(x+\sigma_k^2\mathsf{s}_k(x)\bigr), \qquad \mathsf{D}_k(x)-\mathsf{D}^{\star}_k(x) = a_k^{-1}\sigma_k^2\bigl(\mathsf{s}_k(x)-\mathsf{s}^{\star}_k(x)\bigr).\] In the variance-exploding normalization \(a_t\equiv1\) and \(\sigma_t^2=t\), the reverse drift can be written as \[\mathsf{s}^{\star}_t(x) = \frac{\mathsf{D}^{\star}_t(x)-x}{t}.\] Thus one may either think of the sampler as using a score field, or as using an optimal denoiser whose value is frozen over each discrete reverse step.

Before approximating the reverse step with a score model, let us write the exact Bayes kernel in terms of the marginal density \(p_k\). Let \(P_k^{\leftarrow}(x' \to \cdot)\) be the density of the backward transition \(\mathop{\mathrm{Law}}(X_k\mid X_{k+1}=x')\), and let \(P_k^{\to}(x \to x')\) denote the forward transition density in 32 . For fixed \(x\), \[P_k^{\to}(x \to x') = (2\pi\alpha_k^2\eta_k)^{-d/2} \exp\left( -\frac{\left\lVert x'-\alpha_kx \right\rVert^2}{2\alpha_k^2\eta_k} \right).\] Bayes’ rule gives \[P_k^{\leftarrow}(x' \to x) = \frac{p_k(x)P_k^{\to}(x \to x')}{p_{k+1}(x')}, \qquad p_{k+1}(x')=\int p_k(y)P_k^{\to}(y \to x')\,\,\mathrm{d}y.\] So the exact reverse kernel is \[\label{eq:backward-kernel} P_k^{\leftarrow}(x' \to x) \propto p_k(x) \exp\!\left( -\frac{\left\lVert x-\alpha_k^{-1}x' \right\rVert^2}{2\eta_k} \right).\tag{34}\] Thus each reverse step is a Bayes update: \(p_k\) is the prior at the previous noise level, and the Gaussian factor is the likelihood of the rescaled observation \(\alpha_k^{-1}x'\) under the channel \(\alpha_k^{-1}X_{k+1}=X_k+\sqrt{\eta_k}\xi_k\). If \(p_k\) were known as a density, 34 could be used directly. In diffusion models, the learned object is only the score \(\nabla\log p_k\), so the next step is to approximate this Bayes update locally.

5.2 The Gaussian reverse update↩︎

A common sampler replaces 34 by a Gaussian transition. With the notation of this section, \[\label{eq:ddpm-approx} \widehat{P}_k^{\leftarrow}(x' \to \cdot) = \mathcal{N}\!\left(\alpha_k^{-1}x'+\eta_k\alpha_k\mathsf{s}_{k+1}(x'),\eta_kI\right),\tag{35}\] where we use the same step-size parameter \(\eta_k\) for the Gaussian variance. This should be contrasted with the exact reverse kernel 34 , which need not be Gaussian because of the prior factor \(p_k(x)\).

To see why the mean in 35 has this form, apply Tweedie’s identity to the single forward step itself. After rescaling the next observation, \[\alpha_k^{-1}X_{k+1} = X_k+\sqrt{\eta_k}\xi_k.\] Tweedie’s identity for this channel gives \[\begin{align} \mathbb{E}[X_k\mid X_{k+1}=x'] & = \alpha_k^{-1}x' + \alpha_k \eta_k\nabla\log\bigl( p_{k+1}(x')\bigr) \\ & = \alpha_k^{-1}x' + \alpha_k \eta_k \mathsf{s}^{\star}_{k+1}(x'). \end{align}\] Thus the center of 35 is the exact posterior mean with \(\mathsf{s}^{\star}_{k+1}\) replaced by the learned score \(\mathsf{s}_{k+1}\). The remaining approximation is to replace the exact posterior covariance by the local Gaussian scale \(\eta_kI\): when the prior density \(p_k\) is smooth, the likelihood in 34 confines \(X_k\) to a ball of radius \(\sqrt{\eta_k}\) around the observation, up to leading order.

Let us also present a complementary derivation from the numerical SDE viewpoint. Choose a forward noising SDE on \([0,T]\) whose marginals agree with the channel \(X_t\mid X_0\sim\mathcal{N}\!\left(a_tX_0,\sigma_t^2I\right)\). Assume the schedule is smooth, \(a_t>0\), with the usual normalization \(a_0=1\) and \(\sigma_0=0\), and define \[\lambda_t=\frac{\dot{a}_t}{a_t}, \qquad f_t(x)=\lambda_t x, \qquad g_t^2=\frac{\,\mathrm{d}}{\,\mathrm{d}t}\sigma_t^2-2\lambda_t\sigma_t^2.\] Assume the last quantity is non-negative, so that \(g_t\) is well defined, then the linear SDE \[\,\mathrm{d}X_t=\lambda_tX_t\,\mathrm{d}t+g_t\,\mathrm{d}B_t\] has conditional mean \(a_tX_0\) and conditional variance \(\sigma_t^2I\): the variance solves \(\dot{v}_t=2\lambda_tv_t+g_t^2\) with \(v_0=0\), hence \(v_t=\sigma_t^2\). Applying 13 to this particular choice of \(f_t\) and \(g_t\) gives the reverse SDE \[\,\mathrm{d}Y^{\leftarrow}_s = \left\{ -\lambda_{T-s}Y^{\leftarrow}_s + g_{T-s}^2\nabla\log p_{T-s}(Y^{\leftarrow}_s) \right\}\,\mathrm{d}s + g_{T-s}\,\mathrm{d}B^{\leftarrow}_s.\] Now apply Euler–Maruyama to this reverse SDE on the grid. For the reverse step from time \(t_{k+1}\) to time \(t_k\), write \(h_k=t_{k+1}-t_k\) and start from \(x'\) at noise level \(t_{k+1}\). Freezing the coefficients and the score at the beginning of this reverse step, namely at \(t_{k+1}\) in forward time, and using the grid shorthand \(\lambda_{k+1}:=\lambda_{t_{k+1}}\) and \(g_{k+1}:=g_{t_{k+1}}\), gives \[x' -h_k\lambda_{k+1}x' +h_kg_{k+1}^2\mathsf{s}_{k+1}(x') +g_{k+1}\sqrt{h_k}\,\xi.\] The discrete parameters are exactly the grid versions of the same schedule: \[\alpha_k=\frac{a_{k+1}}{a_k}, \qquad \eta_k=\frac{\sigma_{k+1}^2}{\alpha_k^2}-\sigma_k^2.\] A Taylor expansion at \(t_{k+1}\) gives \[\alpha_k^{-1} = 1-h_k\lambda_{k+1}+O(h_k^2), \qquad \alpha_k\eta_k = h_kg_{k+1}^2+O(h_k^2),\] which agrees with the Gaussian reverse update discussed above. It is worth being explicit about the two step parameters here, since both recur in the error analysis: \(h_k=t_{k+1}-t_k\) is the time increment of the reverse step, while \(\eta_k\) is the variance of the one-step Gaussian kernel 32 . In a general schedule they are different quantities, related by \(\alpha_k\eta_k=h_kg_{k+1}^2+O(h_k^2)\), but they coincide, \(\eta_k=h_k\), in the variance-exploding normalization (\(\alpha_k=1\), \(g_t\equiv1\)) used in Section 6. The update 35 can therefore be viewed either as a local approximation to the exact Bayes kernel or as an Euler–Maruyama discretization of the reverse SDE associated with the same Gaussian noising channel. The discretization analysis asks how much error is introduced by making this local replacement at every reverse step.

6 Error Analysis for Diffusion Models↩︎

The previous section isolated the local numerical question: each step has an exact Bayes kernel 34 , while a practical sampler uses an implementable approximation such as 35 . We now turn this local comparison into a global guarantee. Beyond the early-stopping gap between the smoothed target and the raw data law, the reverse-chain error splits into the same three sources isolated in the previous section: how the chain is initialized, how accurately the score has been learned, and how accurately each reverse kernel is discretized.

Before getting into the details, it is worth laying out the roadmap of the whole analysis; each piece is made precise in the subsections that follow. Throughout, \(P_k^{\leftarrow}(x_{k+1} \to \cdot)\) denotes the exact reverse kernels 34 , and \(\widehat{P}_k^{\leftarrow}(x_{k+1} \to \cdot)\) the implemented reverse kernels of a sampler. For each time \(k\), the \(L^2(p_k)\) score error is \[\label{eq:score-error} \varepsilon_{k,\mathsf{score}}^2 = \mathbb{E}_{X_k\sim p_k} \left\lVert \mathsf{s}_k(X_k)-\mathsf{s}^{\star}_k(X_k) \right\rVert^2.\tag{36}\] The implemented reverse chain is initialized at the high-noise endpoint with error \[\mathop{\mathrm{D_{\mathsf{KL}}}}(p_K\|\widehat{p}_K)\leq \varepsilon_{\mathsf{init}}^2,\] and the accumulated one-step kernel error splits into a pure discretization part and a score-error part, \[\sum_{k=1}^{K-1} \mathbb{E}_{X_{k+1}\sim p_{k+1}} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left( P_k^{\leftarrow}(X_{k+1} \to \cdot) \middle\| \widehat{P}_k^{\leftarrow}(X_{k+1} \to \cdot) \right) \leq \varepsilon_{\mathsf{disc}}^2 + c\sum_{k=1}^{K-1}\eta_k\varepsilon_{k+1,\mathsf{score}}^2.\] Putting these together, the error decomposition we will establish reads \[\label{eq:ddpm-error-decomposition} \mathop{\mathrm{D_{\mathsf{KL}}}}(p_1\|\widehat{p}_1) \leq \varepsilon_{\mathsf{init}}^2 + \varepsilon_{\mathsf{disc}}^2 + c\sum_{k=1}^{K-1}\eta_k\varepsilon_{k+1,\mathsf{score}}^2,\tag{37}\] with the early-stopping comparison between \(p_1\) and \(p_{\mathsf{data}}\) handled separately as a geometric smoothing estimate, naturally controlled in \(W_2\).

This decomposition is mainly a bookkeeping statement. Instead of comparing only the final samples, we compare the exact and implemented reverse chains as full paths. Data processing then says that the KL between the final marginals is no larger than the path-space KL. The path-space chain rule breaks this latter quantity into two pieces: the mismatch between the two initial laws at the high-noise endpoint, and the sum of the one-step KL errors accumulated along the reverse chain.

The factor multiplying the score error has a simple origin. In Gaussian DDPM or Euler–Maruyama updates, the learned score changes the mean of the one-step Gaussian kernel by about \(\eta_k(\mathsf{s}_{k+1}-\mathsf{s}^{\star}_{k+1})\). KL between Gaussians with the same covariance \(\eta_kI\) divides the squared mean error by the variance, leaving a contribution of order \(\eta_k\varepsilon_{k+1,\mathsf{score}}^2\). Thus score errors at times with larger reverse steps matter more, while errors on very small steps are naturally downweighted. Here \(\eta_k\) is the one-step variance parameter of Section 5.1; in the variance-exploding chain analyzed in the rest of this section it is exactly the reverse-time step \(h_k=t_{k+1}-t_k\), which is why these weights act as step sizes.

The term \(\varepsilon_{\mathsf{disc}}^2\) collects the error that would remain even with the exact score. For Euler–Maruyama this is the cost of freezing the score or denoiser during a step. The stochastic-localization argument controls that cost through posterior covariance budgets, and the high-accuracy correction later reduces it by sampling the local Gaussian tilt more faithfully. The subsections below unpack the pieces in this order: early stopping, KL telescoping, score-error contribution, Euler–Maruyama discretization and its localization refinement, and finally the high-accuracy correction.

6.1 Early stopping↩︎

The first issue is the endpoint of the reverse process. We stop at a small positive noise level and target the smoothed law \(p_1\). This early-stopped law is the right object for the KL analysis: path-space chain rules and Girsanov-type estimates give direct control of \(\mathop{\mathrm{D_{\mathsf{KL}}}}(p_1\|\widehat{p}_1)\) for the sampler output \(\widehat{p}_1\). The raw data distribution may be singular, supported on a lower-dimensional set, or otherwise mutually singular with any smooth sampler output, so it cannot generally be compared to \(\widehat{p}_1\) in KL.

Thus we keep two error measurements separate. The early-stopping error is a geometric smoothing error and is controlled in \(W_2\). The algorithmic error between \(p_1\) and \(\widehat{p}_1\) is controlled later in KL.

The \(W_2\) estimate is immediate from the forward noising coupling. In the grid notation of the previous section, \(X_1=a_1X_0+\sigma_1Z\), so coupling \(X_1\) with the same clean sample \(X_0\) gives \[\label{eq:early-stop-w2} W_2^2(p_{\mathsf{data}},p_1) \leq (1-a_1)^2\mathbb{E}\left\lVert X_0 \right\rVert^2+\sigma_1^2d.\tag{38}\] Thus early stopping is harmless when \(\sigma_1\) and \(1-a_1\) are small relative to the scale of the data.

6.2 KL telescoping transition error↩︎

The algorithmic comparison in this section is \[\mathop{\mathrm{D_{\mathsf{KL}}}}(p_1\|\widehat{p}_1),\] where \(p_1\) is the law produced by the exact reverse chain and \(\widehat{p}_1\) is the law produced by the implemented chain, both stopped at the first positive noise level. We control this final KL by first comparing the two reverse chains as full paths and then projecting those paths to their endpoint.

The reason this helps is simple: a Markov chain path density is a product of an initial density and transition densities. Taking a logarithm turns that product into a sum, and taking expectation gives an additive KL identity. This is the basic telescoping mechanism behind the error decomposition: a long reverse sampling procedure is reduced to many one-step reverse-kernel comparisons.

This proposition is the path-space part of the error decomposition 37 . The remaining work is to estimate the one-step terms in the sum.

6.3 How score error contributes to KL↩︎

Proposition 4 reduces the global KL error to local KL comparisons between the exact reverse kernel \(P_k^{\leftarrow}\) and the implemented kernel \(\widehat{P}_k^{\leftarrow}\). Each local comparison can contain several sources of error. One source is numerical: even with the exact score, the implemented transition may only approximate the Bayes reverse kernel. Another source is statistical: the sampler uses a learned score \(\mathsf{s}\) in place of the true score \(\mathsf{s}^{\star}\).

In the rest of this section we specialize the DDPM grid of Section 5.1 to the variance-exploding case (\(\alpha_k=1\), \(\sigma_k^2=t_k\)), now run backward from a large noise level \(T\) down to an early-stopping level \(\delta>0\): \[\delta=t_1<t_2<\cdots<t_K=T,\qquad h_k=t_{k+1}-t_k.\] The reverse step indexed by \(k\) moves from noise level \(t_{k+1}\) to \(t_k\), for \(k=1,\ldots,K-1\), so the chain ends at \(t_1=\delta\); the final law written \(p_1\) in the earlier subsections is thus the early-stopped law \(p_\delta\). The one-step variance parameter then reduces to the time increment, \(\eta_k=h_k\), and we write \(h_k\) from here on to emphasize that the reverse step is a genuine time step of an SDE discretization.

This subsection isolates the statistical part. The implemented reverse step is the Gaussian update 35 – equivalently, an Euler–Maruyama discretization of the reverse SDE – in which the score enters only through the mean of the kernel. A score error therefore becomes a mean error, which we now quantify.

To keep the notation separate, first fix the Gaussian discretization. Let \[\label{eq:em-exact-score} P_k^{\mathsf{EM}}(x' \to \cdot) = \mathcal{N}\!\left(x'+h_k\mathsf{s}^{\star}_{k+1}(x'),h_kI\right)\tag{39}\] be the exact-score Euler kernel, and let \[\label{eq:em-learned-score} \widehat{P}_k^{\mathsf{EM}}(x' \to \cdot) = \mathcal{N}\!\left(x'+h_k\mathsf{s}_{k+1}(x'),h_kI\right)\tag{40}\] be the learned-score Euler kernel. These are not the exact Bayes reverse kernel \(P_k^{\leftarrow}\); they are the true-score and learned-score versions of the same Gaussian approximation 35 . The following KL therefore isolates the statistical score error after the discretization has been fixed. For fixed \(x'\), applying Lemma 4 gives \[\mathop{\mathrm{D_{\mathsf{KL}}}}\!\left( P_k^{\mathsf{EM}}(x' \to \cdot) \middle\| \widehat{P}_k^{\mathsf{EM}}(x' \to \cdot) \right) = \frac{1}{2h_k} \left\lVert h_k \bigl(\mathsf{s}_{k+1}(x')-\mathsf{s}^{\star}_{k+1}(x')\bigr) \right\rVert^2.\] Averaging over \(X_{k+1}\sim p_{k+1}\) and using the definition of \(\varepsilon_{k+1,\mathsf{score}}^2\) gives \[\label{eq:gaussian-score-step} \mathbb{E}_{X_{k+1}\sim p_{k+1}} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left( P_k^{\mathsf{EM}}(X_{k+1} \to \cdot) \middle\| \widehat{P}_k^{\mathsf{EM}}(X_{k+1} \to \cdot) \right) = \frac{h_k}{2} \varepsilon_{k+1,\mathsf{score}}^2,\tag{41}\] which is the score-estimation contribution in the error decomposition 37 , after the discretization choice has been separated, with \(\eta_k=h_k\) and \(c=\tfrac12\).

Summing these one-step terms over the reverse chain, the total score contribution has the form \[\label{eq:score-sum} \sum_{k=1}^{K-1}h_k \varepsilon_{k+1,\mathsf{score}}^2,\tag{42}\] the discrete analogue of the path-space estimate \[\int_\delta^T \mathbb{E}_{X_t\sim p_t} \left\lVert \mathsf{s}_t(X_t)-\mathsf{s}^{\star}_t(X_t) \right\rVert^2 \,\,\mathrm{d}t.\] The weight \(h_k\) is what makes this estimate informative: a score model need not be equally accurate at all times, since an error matters less when the sampler takes only a small step, while an error during a large reverse step has a larger effect.

6.4 Euler–Maruyama↩︎

We now turn to the numerical error of the simplest reverse-time integrator, the Euler–Maruyama scheme, in the variance-exploding setting fixed at the start of this section. Here the initialization term in the error decomposition 37 is explicit. If the implemented reverse chain is initialized from \(\mathcal{N}\!\left(0,TI\right)\), then \[\label{eq:ve-init-kl} \varepsilon_{\mathsf{init}}^2 := \mathop{\mathrm{D_{\mathsf{KL}}}}(p_T\|\mathcal{N}\!\left(0,TI\right)) \leq \int \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left(\mathcal{N}\!\left(x,TI\right)\middle\|\mathcal{N}\!\left(0,TI\right)\right) p_{\mathsf{data}}(\,\mathrm{d}x) = \frac{\mathbb{E}\left\lVert X_0 \right\rVert^2}{2T}.\tag{43}\] Here we used the representation \(p_T=\int\mathcal{N}\!\left(x,TI\right)\,p_{\mathsf{data}}(\,\mathrm{d}x)\) and convexity of KL in its first argument.

The other, and more substantial, contribution is the discretization error accumulated over the reverse steps, whose size is governed by the step lengths \(h_k\). Choosing the grid well is therefore the central design question, and the relevant local scale against which to measure \(h_k\) is the current noise level. Indeed, the reverse drift contains \(\mathsf{s}^{\star}_t=(\mathsf{D}^{\star}_t-x)/t\), so the same absolute step size is much more aggressive near the early-stopping level than it is at large noise. A uniform grid would therefore force all steps to be as small as the lowest noise scale \(\delta\), leading to a poor dependence on \(T/\delta\). The standard non-uniform-grid remedy is to use a relative control, as in Chen, Lee, and Lu [32]: \[\label{eq:relative-mesh} h_k\leq \bar h\,t_{k+1}, \qquad \bar h\leq \frac{1}{2}.\tag{44}\] The canonical example is the geometric grid \[\label{eq:geometric-ve-grid} t_k=(1-\bar h)^{K-k}T, \qquad \delta=t_1=(1-\bar h)^{K-1}T.\tag{45}\] For this grid, \[K\asymp \bar h^{-1}\log\frac{T}{\delta}, \qquad \bar h\asymp \frac{\log(T/\delta)}{K}.\]

With this convention, the exact Bayes reverse kernel is \[\label{eq:ve-backward-kernel} P_k^{\leftarrow}(x' \to \,\mathrm{d}x) \propto_x p_k(x) \exp\!\left(-\frac{\left\lVert x-x' \right\rVert^2}{2h_k}\right)\,\mathrm{d}x,\tag{46}\] where \(x'\) is the state at noise level \(t_{k+1}\). A first-order Gaussian approximation freezes the score at the beginning of the reverse step. This is the exact-score Euler kernel \(P_k^{\mathsf{EM}}\) from 39 ; the learned sampler uses \(\widehat{P}_k^{\mathsf{EM}}\) from 40 , replacing \(\mathsf{s}^{\star}_{k+1}\) by \(\mathsf{s}_{k+1}\). The superscript \(\mathsf{EM}\) is just a reminder that the score has been frozen during the step. In the general DDPM notation of 35 , these formulas are obtained by setting \(\alpha_k=1\) and \(\eta_k=h_k\).

There are two local error sources. The score-estimation part has already been isolated in the previous subsection: in the present variance-exploding normalization, 41 contributes \(h_k\varepsilon_{k+1,\mathsf{score}}^2\) up to constants. The new issue in this subsection is the numerical error caused by Euler–Maruyama itself. Thus we first compare the exact reverse kernel \(P_k^{\leftarrow}\) with the exact-score Euler kernel \(P_k^{\mathsf{EM}}\). This is the error made by freezing the true reverse drift during one time step.

Proposition 5 is only the exact-score discretization estimate. When the sampler uses \(\mathsf{s}_{k+1}\) instead of \(\mathsf{s}^{\star}_{k+1}\), the extra contribution is the statistical Gaussian mean-error term computed in 41 , which becomes \(h_k\varepsilon_{k+1,\mathsf{score}}^2\) in the variance-exploding normalization.

Combining Proposition 5, the score-error contribution 41 , and Proposition 4 gives the standard Euler–Maruyama diffusion-sampling bound: \[\label{eq:em-global-bound} \mathop{\mathrm{D_{\mathsf{KL}}}}(p_\delta\|\widehat{p}_\delta) \lesssim \varepsilon_{\mathsf{init}}^2 + \sum_{k=1}^{K-1} d\,L_k^2h_k^2 + \sum_{k=1}^{K-1} h_k\varepsilon_{k+1,\mathsf{score}}^2.\tag{47}\] The three terms are initialization at the high-noise endpoint, Euler discretization, and score estimation. This is the KL telescoping estimate of Proposition 4 with the local Gaussian calculation inserted, so it is the Euler–Maruyama specialization of the error decomposition 37 .

In fact, the discretization error in estimate 47 can be sharpened, as the next subsection shows. Once initialization and score estimation have been separated, the exact-score Euler sampler pays only for freezing the true score during each reverse step. Let \(t^+\) denote the next grid point at or above the current noise level, so \(t^+=t_{k+1}\) for \(t\in[t_k,t_{k+1}]\). The discretization cost is measured by the pathwise quantity \[\label{eq:cll-bottleneck} \int_\delta^T \mathbb{E}\left\lVert \mathsf{s}^{\star}_t(X_t)-\mathsf{s}^{\star}_{t^+}(X_{t^+}) \right\rVert^2\,\mathrm{d}t,\tag{48}\] up to deterministic schedule factors in more general schedules. This is the pathwise version of the finite-grid \(\sum_k dL_k^2h_k^2\) term: Euler–Maruyama freezes the score at the high-noise endpoint \(t^+\) of each reverse step, and the discretization error is governed by the change in the true score along the ideal reverse path as the noise time decreases from \(t^+\) to \(t\).

In the error analysis above, the passage from this score difference to the bound \(dL_k^2h_k^2\) used a worst-case control: the score was bounded by a Lipschitz constant \(L_k\), or equivalently by a pointwise Hessian bound for \(\log p_t\). The factor \(d\) comes from the Brownian displacement during one step; the potential loss comes from multiplying it by a uniform operator-norm bound on the Hessian. Under weak assumptions, this worst-case curvature can be much larger than the average curvature seen by the reverse process, especially near the early-stopping time. The next subsection replaces this pointwise control by an averaged posterior-covariance quantity, leading to the nearly linear-in-\(d\) discretization bound.

6.5 Hessian control↩︎

The difficulty left by the previous subsection is the pathwise score difference in 48 . The aim of this subsection is to replace the worst-case Hessian estimate by the more averaged bound \[\text{score-freezing cost} \;\lesssim\; \bar h \times \text{posterior-covariance budget}.\] We first derive the local score-freezing estimate where the factor \(\bar h\) appears, and then rewrite the remaining weighted quadratic variation in terms of posterior covariance.

For \(t\in[t_k,t_{k+1}]\), write \(t^+=t_{k+1}\). The Girsanov score-freezing cost in 48 is, up to the constant \(\frac{1}{2}\) and schedule factors, the sum over grid intervals of \[\int_{t_k}^{t_{k+1}} \mathbb{E}\left\lVert \mathsf{s}^{\star}_t(X_t)-\mathsf{s}^{\star}_{t^+}(X_{t^+}) \right\rVert^2\,\,\mathrm{d}t.\]

The useful variable is not the score itself but the optimal denoiser. In the variance-exploding normalization, \[\mathsf{D}^{\star}_t(x):=\mathbb{E}[X_0\mid X_t=x]=x+t\mathsf{s}^{\star}_t(x), \qquad \mathsf{s}^{\star}_t(x)=\frac{\mathsf{D}^{\star}_t(x)-x}{t}.\] The stochastic-localization martingale identity 25 , rewritten in the noise-time coordinate, says that for \(t\leq t^+\), \[\mathbb{E}\left\lVert \mathsf{D}^{\star}_t(X_t)-\mathsf{D}^{\star}_{t^+}(X_{t^+}) \right\rVert^2 = \int_t^{t^+} \mathbb{E}\left\lVert \nabla \mathsf{D}^{\star}_r(X_r) \right\rVert_{\mathrm F}^2\,\,\mathrm{d}r.\] This is the replacement for a pointwise Hessian bound: it measures the actual quadratic variation of the denoiser along the noising path.

Now estimate one grid interval. Using \(\mathsf{s}^{\star}_t=(\mathsf{D}^{\star}_t-x)/t\) and keeping only the curvature term explicitly, the score increment is bounded by the denoiser increment with the natural \(t^{-2}\) weight: \[\mathbb{E}\left\lVert \mathsf{s}^{\star}_t(X_t)-\mathsf{s}^{\star}_{t^+}(X_{t^+}) \right\rVert^2 \lesssim \frac{1}{t^2} \mathbb{E}\left\lVert \mathsf{D}^{\star}_t(X_t)-\mathsf{D}^{\star}_{t^+}(X_{t^+}) \right\rVert^2.\] Here we have omitted some error terms coming from the explicit factor \(x/t\), which can be handled easily. Substituting the martingale identity and applying Fubini gives \[\begin{align} &\int_{t_k}^{t_{k+1}} \mathbb{E}\left\lVert \mathsf{s}^{\star}_t(X_t)-\mathsf{s}^{\star}_{t^+}(X_{t^+}) \right\rVert^2\,\,\mathrm{d}t \\ &\qquad\lesssim \int_{t_k}^{t_{k+1}} \frac{1}{t^2} \int_t^{t_{k+1}} \mathbb{E}\left\lVert \nabla \mathsf{D}^{\star}_r(X_r) \right\rVert_{\mathrm F}^2\,\,\mathrm{d}r\,\,\mathrm{d}t \\ &\qquad= \int_{t_k}^{t_{k+1}} \mathbb{E}\left\lVert \nabla \mathsf{D}^{\star}_r(X_r) \right\rVert_{\mathrm F}^2 \left(\int_{t_k}^r\frac{\,\mathrm{d}t}{t^2}\right)\,\mathrm{d}r \\ &\qquad= \int_{t_k}^{t_{k+1}} \mathbb{E}\left\lVert \nabla \mathsf{D}^{\star}_r(X_r) \right\rVert_{\mathrm F}^2 \frac{r-t_k}{r\,t_k}\,\,\mathrm{d}r. \end{align}\] Since \(r-t_k\leq h_k\) and the mesh condition 44 with \(\bar h\leq1/2\) gives \(t_k\geq t_{k+1}/2\), we get \[\int_{t_k}^{t_{k+1}} \mathbb{E}\left\lVert \mathsf{s}^{\star}_t(X_t)-\mathsf{s}^{\star}_{t^+}(X_{t^+}) \right\rVert^2\,\,\mathrm{d}t \lesssim \frac{h_k}{t_{k+1}} \int_{t_k}^{t_{k+1}} \frac{\mathbb{E}\left\lVert \nabla \mathsf{D}^{\star}_r(X_r) \right\rVert_{\mathrm F}^2}{r}\,\,\mathrm{d}r.\] This is the point at which \(\bar h\) enters: the ratio \(h_k/t_{k+1}\) is the relative step size on the interval, and 44 says it is at most \(\bar h\). Summing over \(k\) therefore gives \[\sum_k \int_{t_k}^{t_{k+1}} \mathbb{E}\left\lVert \mathsf{s}^{\star}_t(X_t)-\mathsf{s}^{\star}_{t^+}(X_{t^+}) \right\rVert^2\,\,\mathrm{d}t \lesssim \bar h \int_\delta^T \frac{\mathbb{E}\left\lVert \nabla \mathsf{D}^{\star}_t(X_t) \right\rVert_{\mathrm F}^2}{t}\,\,\mathrm{d}t.\]

It remains to express the weighted quadratic variation in a more intrinsic form. Let \[\Sigma_t(x)=\mathop{\mathrm{Cov}}(X_0\mid X_t=x)\] be the posterior covariance in the variance-exploding noising channel. The linear-response identity 23 and the score-Jacobian identity 24 , translated to this notation, give \[\nabla \mathsf{D}^{\star}_t(x) = I+t\nabla\mathsf{s}^{\star}_t(x) = t^{-1}\Sigma_t(x).\] The covariance-dissipation identity 26 , \(\frac{\,\mathrm{d}}{\,\mathrm{d}u}\mathbb{E}[\mathop{\mathrm{Tr}}\Sigma_u] = -\mathbb{E}[\mathop{\mathrm{Tr}}(\Sigma_u^2)]\), after the change of variables \(t=1/u\), gives \[\mathbb{E}\left\lVert \nabla \mathsf{D}^{\star}_t(X_t) \right\rVert_{\mathrm F}^2 = \partial_t\,\mathbb{E}\mathop{\mathrm{Tr}}\Sigma_t(X_t).\] Thus \[\int_\delta^T \frac{\mathbb{E}\left\lVert \nabla \mathsf{D}^{\star}_t(X_t) \right\rVert_{\mathrm F}^2}{t} \,\,\mathrm{d}t = \int_\delta^T \frac{\partial_t\,\mathbb{E}\mathop{\mathrm{Tr}}\Sigma_t(X_t)}{t}\,\,\mathrm{d}t.\] Integrating by parts, \[\int_\delta^T \frac{\partial_t\,\mathbb{E}\mathop{\mathrm{Tr}}\Sigma_t(X_t)}{t}\,\,\mathrm{d}t = \frac{\mathbb{E}\mathop{\mathrm{Tr}}\Sigma_T(X_T)}{T} - \frac{\mathbb{E}\mathop{\mathrm{Tr}}\Sigma_\delta(X_\delta)}{\delta} + \int_\delta^T \frac{\mathbb{E}\mathop{\mathrm{Tr}}\Sigma_t(X_t)}{t^2}\,\,\mathrm{d}t.\] The middle term is nonpositive, so the first and third terms give an upper bound. Following the notation of Chewi [2], define the covariance budget \[\label{eq:D-delta-T} \mathfrak D_{\delta,T}(p_{\mathsf{data}}) := \frac{\mathbb{E}_{X_T\sim p_T}\mathop{\mathrm{Tr}}\Sigma_T(X_T)}{T} + \int_\delta^T \frac{\mathbb{E}_{X_t\sim p_t}\mathop{\mathrm{Tr}}\Sigma_t(X_t)}{t^2}\,\,\mathrm{d}t.\tag{49}\] Combining the previous displays gives the clean summary \[\sum_k \int_{t_k}^{t_{k+1}} \mathbb{E}\left\lVert \mathsf{s}^{\star}_t(X_t)-\mathsf{s}^{\star}_{t^+}(X_{t^+}) \right\rVert^2\,\,\mathrm{d}t \lesssim \bar h\,\mathfrak D_{\delta,T}(p_{\mathsf{data}}).\] This is the key identification: the discretization error is governed by the weighted quadratic variation of the optimal denoiser, and \(\mathfrak D_{\delta,T}\) packages that variation after the DDPM time weights have been accounted for. A careful version of this estimate, together with the score-error term, gives \[\label{eq:chewi-ddpm-D} \mathop{\mathrm{D_{\mathsf{KL}}}}(p_\delta\|\widehat{p}_\delta) \leq \varepsilon_{\mathsf{init}}^2 + 4\sum_{k=1}^{K-1} h_k\, \mathbb{E}_{p_{k+1}} \left\lVert \mathsf{s}_{k+1}-\mathsf{s}^{\star}_{k+1} \right\rVert^2 + 2\bar h\,\mathfrak D_{\delta,T}(p_{\mathsf{data}}),\tag{50}\] where the middle term is the time-weighted score error. The important lesson is that \(\mathfrak D_{\delta,T}(p_{\mathsf{data}})\) replaces a global Lipschitz constant for the score as the data-dependent factor in the discretization complexity.

The remaining question is how large this covariance budget can be. A useful ambient-dimensional estimate comes directly from its posterior-variance form: for any distribution with finite second moment, \[\mathbb{E}\mathop{\mathrm{Tr}}\Sigma_t(X_t) = \mathbb{E}\left\lVert X_0-\mathbb{E}[X_0\mid X_t] \right\rVert^2 \leq \mathbb{E}\left\lVert X_0-X_t \right\rVert^2 = dt,\] because \(X_t\) itself is an admissible estimator of \(X_0\). Substituting this into 49 gives \[\mathfrak D_{\delta,T}(p_{\mathsf{data}}) \leq d\left(1+\log\frac{T}{\delta}\right).\] Combining this estimate with 50 and the geometric-grid relation 45 , the exact-score discretization term is controlled by \[O\!\left( \frac{d\log^2(T/\delta)}{K} \right).\] Thus, when the initialization and score-estimation errors are also \(O(\varepsilon^2)\), it is enough to take \[K=\widetilde{O}\!\left( \frac{d\log^2(T/\delta)}{\varepsilon^2} \right)\] reverse steps. This is the nearly \(d\)-linear dependence proved by Benton, De Bortoli, Doucet, and Deligiannidis [25] under a finite second-moment assumption. Related KL guarantees under finite Fisher-information assumptions were obtained by Conforti, Durmus, and Gentiloni Silveri [33]. The covariance budget can be bounded more sharply when the data are intrinsically low-dimensional: entropy or covering-number estimates for the smoothed law \(p_\delta\) replace the ambient dimension by an intrinsic dimension at scale \(\sqrt\delta\).

6.6 First-order rejection sampling and high accuracy↩︎

The Hessian-control bound explains why Euler–Maruyama can be analyzed under weak assumptions, but it also exposes a limitation of the sampler itself. Euler–Maruyama does not try to sample the exact one-step Bayes kernel; it linearizes the log density and relies on the step size being small enough that the nonlinear remainder is negligible. Thus, as for ULA before the Metropolis correction, high accuracy is bought by a fine grid and hence by a polynomial dependence on the target accuracy. The natural MALA-like question is whether this dependence can be reduced to logarithmic, or at least polylogarithmic, dependence on \(1/\varepsilon\).

Here the analogy with MALA also reveals the main obstruction. MALA obtains its correction from a log-density ratio. In a diffusion model, the one-step target involves the noised density \(p_k\), but \(p_k(x)\) and \(\log p_k(x)\) are not available as numerical quantities; the learned object is only the score \(\mathsf{s}^{\star}_k=\nabla\log p_k\). Thus the question is sharper than simply applying a Metropolis step: can a score-only correction emulate the missing log-density-ratio test and correct the local Gaussian proposal toward the exact Bayes kernel?

For the variance-exploding chain, this question has a useful local structure. The exact backward kernel 46 is not merely a Gaussian with a shifted mean; it is a Gaussian tilt of the form \[\label{eq:gaussian-tilt} p_k(x)\exp\!\left(-\frac{\left\lVert x-x' \right\rVert^2}{2h_k}\right),\tag{51}\] where \(x'\) is the state at noise level \(t_{k+1}\). Euler–Maruyama corresponds to replacing this tilt by the Gaussian obtained from a first-order expansion of \(\log p_k\). The high-accuracy idea is to sample the tilted law 51 more faithfully while still using only first-order information.

This is where the algorithmic innovation enters. First-order rejection sampling, denoted \(\textsf{FORS}\) [34], replaces direct log-density-ratio evaluation by randomized first-order estimates. Although \(\log p_k(x)\) and \(\log p_k(y)\) are unavailable separately, their difference can be written as a path integral involving only the score: for any smooth path \(\gamma\) with \(\gamma_0=y\) and \(\gamma_1=x\), \[\log p_k(x)-\log p_k(y) = \int_0^1 \left\langle \dot{\gamma}_r,\mathsf{s}^{\star}_k(\gamma_r) \right\rangle\,\,\mathrm{d}r.\] Sampling a random point along the path turns this integral into an unbiased estimate of the log-density difference. Thus, even though the endpoint log densities themselves are unavailable, score queries can produce a randomized estimate of exactly the log tilt that a rejection correction needs.

There is still one conversion to make. Rejection sampling needs an acceptance weight proportional to the exponential of the log tilt, not just an estimate of the log tilt itself. Exponentiating an unbiased estimate would not usually give the right average, so \(\textsf{FORS}\) uses a Poisson product identity to build the exponential tilt from first-order random estimates. In abstract form, suppose a proposal is drawn from a law \(Q\) and the desired law is the exponential tilt of \(Q\) with Radon–Nikodym derivative proportional to \(e^{w(x)}\). Assume that, given the proposal \(x\), one can sample bounded random variables \(W_1,W_2,\ldots\) with \(\mathbb{E}[W_1\mid x]=w(x)\). In the simplest bounded form, assume \(|W_j|\leq R\) almost surely for a deterministic envelope \(R\). Take \(J\sim\mathsf{Poisson}(2R)\) and accept a proposal \(x\) with probability \[\prod_{j=1}^{J}\frac{R+W_j}{2R}.\] The probability generating function of a Poisson random variable makes the average acceptance factor proportional to \(e^{w(x)}\). Hence the accepted proposals follow the exponentially tilted law. If \(R=\Theta(1)\), the number of first-order queries is constant in expectation and logarithmic with high probability.

Thus \(\textsf{FORS}\) is the score-only analogue of a Metropolis correction. MALA uses density ratios to correct a Gaussian proposal; \(\textsf{FORS}\) uses randomized first-order estimates to correct a Gaussian tilt without evaluating the density. Instead of accepting the Gaussian linearization as the sampler, the algorithm corrects that local proposal toward the true Gaussian tilt.

In the notation of the error decomposition 37 , \(\textsf{FORS}\) changes the size of the exact-score algorithmic term \(\varepsilon_{\mathsf{disc}}^2\). By sampling each local Gaussian tilt accurately enough, the accumulated local algorithmic error can be made \(O(\varepsilon^2)\) using only polylogarithmically many reverse steps in the target accuracy \(\varepsilon\). In a simplified ambient-dimension form, the result of Chen, Chewi, Daskalakis, and Rakhlin [34] says that, suppressing schedule constants and taking the initialization and local algorithmic budgets of order \(\varepsilon^2\), there is a reverse sampler using score evaluations for which \[\label{eq:high-acc-prototype} \mathop{\mathrm{D_{\mathsf{KL}}}}(p_\delta\|\widehat{p}_\delta) \lesssim \varepsilon^2 + \sum_{k=1}^{K-1} h_k\varepsilon_{k+1,\mathsf{score}}^2\tag{52}\] with \[K = \widetilde{O}\!\left(d\,\mathop{\mathrm{polylog}}(1/\varepsilon)\right)\] in a basic ambient-dimension reading. The full theorem contains sharper dimension measures and refinements under additional structure, but the key message here is simpler: the expensive dependence on the accuracy parameter is replaced by a polylogarithmic one. Together with the variance-exploding early-stopping estimate \(W_2^2(p_{\mathsf{data}},p_\delta)\leq \delta d\), this separates the two error measurements: the initial smoothing is controlled in \(W_2\), while the implemented reverse chain is controlled in KL. Note that the same time-integrated score-error term appears in 52 : \(\textsf{FORS}\) improves the local sampling step, but it does not remove the need for an accurate learned score.

7 Discrete Diffusion Models↩︎

The development so far built and analyzed diffusion models on continuous state spaces. We now turn to a different modeling question: what if the state space itself is not continuous? Reverse-time Bayes rules, denoising objectives, and path-space error decompositions survive on finite state spaces; Euclidean gradients and Brownian motion do not. Discrete diffusion models therefore force us to state the sampling ideas in a form that does not rely on calculus.

7.1 Forward and reverse kernels on a finite state space↩︎

Images in pixel space can often be treated as continuous, but language, symbolic sequences, molecules, graphs, and many scientific configurations are inherently discrete. If \(x\) is a token sequence, an expression such as \(x+\sqrt t Z\) is meaningless. There is no small Gaussian perturbation of the word “cat” inside the vocabulary. Instead of adding Brownian noise, a discrete diffusion model applies a sequence of Markov kernels on a finite state space.

The probabilistic skeleton is exactly the one we have used throughout: a forward chain that degrades the data toward a reference law, and a reverse chain recovered by Bayes’ rule. We therefore keep the same symbols as in the continuous case. The forward transition \(P_k^{\to}\) and the backward transition \(P_k^{\leftarrow}\) of Subsection 5.1 are now read as stochastic matrices on the finite set, rather than as transition densities on \(\mathbb{R}^d\). The only genuine loss is the calculus: there is no gradient \(\nabla\log p_k\), so the score must be replaced by a finite-difference object.

Let \(\mathcal{X}\) be a finite state space. The data law \(p_{\mathsf{data}}\) may be supported on a subset \(\mathcal{X}_0\subseteq\mathcal{X}\); we view it as a law on \(\mathcal{X}\) by assigning zero mass outside \(\mathcal{X}_0\). For ordinary categorical sequences one can take \(\mathcal{X}_0=\mathcal{X}=\mathcal{V}^L\), while for the masked diffusions treated below the noisy state space \(\mathcal{X}\) also contains mask symbols. A discrete forward diffusion is a Markov chain \[\label{eq:discrete-forward} X_0\sim p_{\mathsf{data}},\qquad X_{k+1}\sim P_k^{\to}(X_k \to \cdot), \qquad k=0,\ldots,K-1,\tag{53}\] where each \(P_k^{\to}\) is a stochastic matrix on \(\mathcal{X}\): \[P_k^{\to}(x \to y)\geq0,\qquad \sum_{y\in\mathcal{X}}P_k^{\to}(x \to y)=1.\] Let \(p_k\) be the law of \(X_k\). In row-vector notation, \[p_{k+1}=p_kP_k^{\to},\qquad p_k=p_{\mathsf{data}}P_0^{\to}P_1^{\to}\cdots P_{k-1}^{\to}.\] The forward kernels are chosen so that \(p_K\) is close to a simple reference distribution \(\pi_{\mathsf{ref}}\), such as all-mask tokens, independent uniform tokens, or a product distribution.

The design principle is the same as in continuous diffusion. The forward process should be simple enough that we can sample from it and know its transition probabilities, and destructive enough that after many steps the distribution is close to a reference law. The reverse process is then the learned part. What changes is the nature of corruption: tokens are replaced, masked, or resampled rather than moved by small Gaussian increments.

7.2 Exact reverse kernels and ratio scores↩︎

The exact reverse kernel is again obtained from Bayes’ rule. For fixed \(x'\in\mathcal{X}\), \[\label{eq:discrete-reverse} P_k^{\leftarrow}(x' \to x) = \mathbb{P}(X_k=x\mid X_{k+1}=x') = \frac{p_k(x)P_k^{\to}(x \to x')}{p_{k+1}(x')}.\tag{54}\] This is the discrete analogue of 34 , in the same \(P_k^{\leftarrow}\) notation: there the continuous backward kernel equalled \(p_k\) times a Gaussian likelihood, renormalized; here it equals \(p_k\) times a transition likelihood \(P_k^{\to}(x \to x')\), renormalized.

If we knew \(p_k\) exactly, 54 would give an exact reverse sampler. In learning, we approximate either the reverse kernel \(P_k^{\leftarrow}(x' \to \cdot)\) directly or an object from which it can be computed, such as the posterior distribution of \(X_0\) given \(X_k\).

The reverse-kernel identity 54 shows exactly what must be learned. The forward transition \(P_k^{\to}\) is known by design, but the intermediate law \(p_k\) contains information about the data distribution and is unknown. Learning a discrete diffusion model means learning enough about these intermediate laws to simulate the reverse Bayes kernels. One way to represent this missing information is through ratios of the unknown law. Indeed, for two candidate predecessors \(x\) and \(y\) of the same current state \(x'\), the normalizing factor \(p_{k+1}(x')\) in 54 cancels: \[\frac{P_k^{\leftarrow}(x' \to y)}{P_k^{\leftarrow}(x' \to x)} = \frac{p_k(y)P_k^{\to}(y \to x')}{p_k(x)P_k^{\to}(x \to x')}.\] Since \(P_k^{\to}\) is known, the unknown part is the same-time density ratio \[\label{eq:discrete-ratio-score} \mathsf{s}^{\star}_k(x,y)=\frac{p_k(y)}{p_k(x)},\qquad x,y\in\mathcal{X}.\tag{55}\] In continuous space, the score tells us the infinitesimal log-density change when we move from \(x\) to \(x+\,\mathrm{d}x\). On a graph or finite set, the ratio \(p_k(y)/p_k(x)\) tells us the finite log-density change when we move from \(x\) to a neighboring state \(y\). If the forward kernel permits only local changes, these ratios are enough to specify local reverse transition probabilities.

For many modern discrete diffusion models, however, the denoising posterior is the more convenient target: \[\mathsf{D}^{\star}_k(x_0\mid x)=\mathbb{P}(X_0=x_0\mid X_k=x), \qquad x_0\in\mathcal{X}_0,\;x\in\mathcal{X}.\] This is the discrete counterpart of the continuous posterior mean \(\mathbb{E}[X_0\mid X_k=x]\) in Tweedie’s identity.

The ratio-score and denoiser views are closely related by Bayes’ rule: \[\mathsf{D}^{\star}_k(x_0\mid x) = \frac{p_{\mathsf{data}}(x_0)\mathbb{P}(X_k=x\mid X_0=x_0)}{p_k(x)}.\] Thus comparing the posterior probabilities at two noisy states recovers the density ratio \(p_k(y)/p_k(x)\) up to the known forward-likelihood ratio. The denoiser therefore contains the ratio-score information, but in a form that is often easier to learn by supervised clean-token prediction.

7.3 Masked diffusion↩︎

The cleanest example is masked diffusion for token sequences. Let the vocabulary be \(\mathcal{V}\), and add a special mask symbol \([\mathsf M]\). Here \(\mathcal{X}_0=\mathcal{V}^L\), while the noisy state space is the augmented sequence space \[\mathcal{X}_{\mathsf M}:=(\mathcal{V}\cup\{[\mathsf M]\})^L.\] For \(x,y\in\mathcal{X}_{\mathsf M}\), write \(x_i\) and \(y_i\) for their \(i\)th coordinates. Define the forward step by \[\label{eq:mask-forward} P_k^{\to}(x\to y) = \prod_{i=1}^L \begin{cases} 1-\beta_k, & y_i=x_i,\;x_i\in\mathcal{V},\\ \beta_k, & y_i=[\mathsf M],\;x_i\in\mathcal{V},\\ 1, & x_i=y_i=[\mathsf M],\\ 0, & \text{otherwise}. \end{cases}\tag{56}\] Thus an unmasked token either stays unchanged or becomes masked, and once a token is masked it remains masked. Variants may use coordinate-dependent masking rates, but the same sequence-level form applies.

Let \[\bar\alpha_k=\prod_{s=0}^{k-1}(1-\beta_s).\] For each coordinate \(i\), conditional on \(X_{0,i}=a\in\mathcal{V}\), \[\mathbb{P}(X_{k,i}=a\mid X_{0,i}=a)=\bar\alpha_k,\qquad \mathbb{P}(X_{k,i}=[\mathsf M]\mid X_{0,i}=a)=1-\bar\alpha_k.\] This is analogous to the Gaussian formula \(X_k\mid X_0\sim\mathcal{N}\!\left(a_kX_0,\sigma_k^2I\right)\): the scalar attenuation measures how much information about \(X_0\) remains.

Masked diffusion is especially transparent because corruption is irreversible in the forward direction. Once a token is masked, the forward process has forgotten its identity. The reverse model must therefore use the surrounding context and the learned data distribution to infer plausible clean tokens. This makes the Bayesian nature of denoising visible without any calculus. At the coordinate level, the reverse step has only two cases.

• If \(X_{k+1,i}=a\in\mathcal{V}\) is unmasked, then necessarily \(X_{k,i}=a\). There is nothing to sample at that coordinate.

• If \(X_{k+1,i}=[\mathsf M]\), then \(X_{k,i}\) may either already be masked or may be the original clean token. The reverse sampler must decide whether to unmask, and if so which token to place.

Specializing the reverse kernel 54 to 56 gives the sequence-level formula \[\label{eq:mask-reverse-sequence} P_k^{\leftarrow}(x'\to x) = \frac{p_k(x)}{p_{k+1}(x')} P_k^{\to}(x\to x').\tag{57}\] This is still a coordinatewise masking rule through the known factor \(P_k^{\to}(x\to x')\), but the unknown data-dependent weight is the full sequence law \(p_k(x)\). Thus the distribution of an unmasked token generally depends on the entire partially masked sequence.

7.4 Training objective for masked diffusion↩︎

The model is usually trained to predict masked clean tokens from a corrupted sequence. Let \(M_k\subseteq\{1,\ldots,L\}\) be the set of masked coordinates at time \(k\), and let \(x_k\) be the masked sequence. A denoising model \(\mathsf{D}_k(i,\cdot\mid x_k)\) outputs a distribution over vocabulary tokens for coordinate \(i\). The population loss is \[\label{eq:masked-loss} \mathsf L_{\mathsf{mask}}(\mathsf{D}) = \mathbb{E}\left[ \sum_{i\in M_k} -\log \mathsf{D}_k(i,X_{0,i}\mid X_k) \right].\tag{58}\] This is cross-entropy with the clean token as the label. The population minimizer is the coordinate marginal of the posterior, \[\mathsf{D}^{\star}_k(i,a\mid x_k) = \mathbb{P}(X_{0,i}=a\mid X_k=x_k).\]

Thus the cross-entropy objective is not an arbitrary language-modeling loss. It is the discrete analogue of denoising score matching. In the Gaussian case, the optimal prediction is a posterior mean and can be converted into a score by Tweedie’s identity. In the masked-token case, the optimal prediction is the full posterior distribution of the hidden clean token.

7.5 Sampling from a masked model↩︎

A trained denoiser turns the posterior predictions above into a reverse-time update rule. In the absorbing case the terminal state is the all-mask sequence. Starting from that state, a sampler repeatedly chooses a set \(A\subseteq M_k\) of currently masked coordinates and fills those coordinates by drawing \[\widehat{x}_i\sim \mathsf{D}_k(i,\cdot\mid x_k), \qquad i\in A.\] Coordinates outside \(A\) are left unchanged, unless the chosen algorithm deliberately remasks or resamples them. Thus the schedule is the rule that chooses \(A\) and the time label \(k\) at each stage of the reverse process.

This viewpoint makes the any-order character of masked diffusion explicit. The model is not tied to a left-to-right factorization: it is trained to denoise from partially masked contexts, so at inference time the sampler may reveal coordinates in any order. This order-agnostic view goes back to the order-agnostic training of NADE [35], and was tied directly to absorbing/masked diffusion by the autoregressive diffusion models of Hoogeboom et al. [36]. For an ordering or sequence of blocks \(A_1,A_2,\ldots,A_m\subseteq\{1,\ldots,L\}\), the sampler updates the still-masked coordinates in \(A_j\) using posterior predictions conditioned on the current partial sequence. Singleton blocks give an autoregressive-like sampler with a chosen order; larger blocks give parallel decoding; adaptive blocks chosen from model confidence give a coarse-to-fine or easy-to-hard reveal schedule.

The reason this freedom is legitimate is that training never fixes an order in the first place. The objective 58 masks a random subset of coordinates and asks the model to predict the clean tokens at the masked positions from the unmasked ones. Across training examples, noise levels, and random masks, the masked set \(M_k\) ranges over essentially all subsets of \(\{1,\ldots,L\}\), so the model is implicitly trained to approximate the entire family of conditionals \[\mathbb{P}\!\left(X_{0,i}=a \mid X_{k,M_k^c}=x_{M_k^c}\right), \qquad M_k\subseteq\{1,\ldots,L\},\;i\in M_k,\] of a clean coordinate \(i\) given the visible entries of the partially masked sequence. Two features of masked diffusion make this exact. First, the forward process only masks tokens, never alters them, so \(X_{k,M_k^c}=X_{0,M_k^c}\): conditioning on the visible part of \(X_k\) is the same as conditioning on the corresponding clean values. Second, whether a coordinate is masked is decided independently of the token values, so the masking pattern carries no information beyond which coordinates are observed. Hence the population optimum \(\mathsf{D}^{\star}_k(i,\cdot\mid x_k)\) of 58 is exactly the above conditional.

An inference order is therefore just a rule for choosing which conditional to query next. Revealing coordinate \(i\) from the currently observed set \(M_k^c\) uses \(\mathbb{P}(X_{0,i}\mid X_{k,M_k^c}=x_{M_k^c})\), a member of the family the model already learned. Every order factorizes the joint into conditionals drawn from this one family, which is why any order is supported by the same trained model. The choice of block size, however, is still a sampling choice. If a block contains several coordinates and the sampler draws them independently from the same context, it is using a factorized approximation to the joint conditional for that block. The next two subsections make this distinction precise: first by writing the general KL accounting for approximate reverse kernels, and then by isolating the extra error created by parallel block updates.

7.6 Error analysis parallel to the continuous case↩︎

A schedule together with a denoiser defines approximate reverse kernels \(\widehat{P}_k^{\leftarrow}(x' \to \cdot)\), while the exact kernels are \(P_k^{\leftarrow}(x' \to \cdot)\) from 54 . Let \(\widehat{p}_0\) be the final distribution generated by the approximate reverse chain, and let \(\widehat{p}_K\) be the law used to initialize that chain at time \(K\). A KL path-space argument gives \[\label{eq:discrete-error} \mathop{\mathrm{D_{\mathsf{KL}}}}(p_0\|\widehat{p}_0) \leq \mathop{\mathrm{D_{\mathsf{KL}}}}(p_K\|\widehat{p}_K) + \sum_{k=0}^{K-1} \mathbb{E}_{X_{k+1}\sim p_{k+1}} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left( P_k^{\leftarrow}(X_{k+1} \to \cdot) \middle\| \widehat{P}_k^{\leftarrow}(X_{k+1} \to \cdot) \right).\tag{59}\] This is the exact analogue of the continuous diffusion decomposition: \[\begin{array}{c|c|c} \text{continuous diffusion} & \text{discrete diffusion} & \text{meaning}\\ \hline p_K\approx\mathcal{N}\!\left(0,I\right) & p_K\approx\pi_{\mathsf{ref}} & \text{initialization at noise}\\ P_k^{\leftarrow}\propto p_k\times\text{Gaussian likelihood} & P_k^{\leftarrow}\propto p_k\times P_k^{\to} & \text{exact reverse kernel}\\ \mathsf{s}_k\approx\mathsf{s}^{\star}_k & \widehat{P}_k^{\leftarrow}\approx P_k^{\leftarrow} \text{ or }\mathsf{D}_k\approx\mathsf{D}^{\star}_k & \text{learned denoising information}\\ \sum_k\eta_k\varepsilon_{k,\mathsf{score}}^2 & \sum_k\mathbb{E}\,\mathop{\mathrm{D_{\mathsf{KL}}}}(P_k^{\leftarrow}\|\widehat{P}_k^{\leftarrow}) & \text{statistical/model error} \end{array}\]

The justification is verbatim the KL telescoping argument of Section 6: apply Lemma 1 to the exact and approximate reverse path laws, then project the paths to their final state using data processing.

The decomposition 59 isolates the same two error sources as in the continuous case: the initialization gap \(\mathop{\mathrm{D_{\mathsf{KL}}}}(p_K\|\widehat{p}_K)\) and the per-step reverse-kernel KL summed over the reverse steps. In absorbing masked diffusion, finiteness of these KL terms is already informative. If the sampler starts from the deterministic all-mask state, then \(\widehat{p}_K=\delta_{[\mathsf M]^L}\), so the initialization term is finite only when the forward terminal law \(p_K\) is also supported on the all-mask state. For the independent masking chain 56 , this requires \(\bar\alpha_K=0\), for instance through a terminal step with \(\beta_k=1\).

Once the supports are compatible so that the KL terms are finite, the second term in 59 is tied to the quantity the model actually trains on. For masked diffusion the model does not learn the reverse kernel directly; it learns posterior clean-token distributions through the cross-entropy objective 58 . Schematically, we have the bound \[\mathbb{E}\,\mathop{\mathrm{D_{\mathsf{KL}}}}(P_k^{\leftarrow}(X_{k+1} \to \cdot)\|\widehat{P}_k^{\leftarrow}(X_{k+1} \to \cdot)) \lesssim \mathbb{E}\sum_{i\in M_{k+1}} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left( \mathsf{D}^{\star}_{k+1}(i,\cdot\mid X_{k+1}) \middle\| \mathsf{D}_{k+1}(i,\cdot\mid X_{k+1}) \right).\] This bound is the cross-entropy regret of the posterior predictions, the part that training controls: it vanishes as the learned posterior \(\mathsf{D}_{k+1}\) approaches the true posterior \(\mathsf{D}^{\star}_{k+1}\). It is the entire per-step error when the sampler unmasks one coordinate at a time, because each single-coordinate reverse step queries a true conditional (as explained in Subsection 7.5), so the approximate reverse chain reproduces the exact one. With a perfect denoiser, single-coordinate decoding is therefore exact in any order. Its drawback is cost—revealing a single coordinate per step takes as many network evaluations as there are coordinates—which raises the question of the next subsection: can we instead reveal several coordinates at once, and what accuracy does that sacrifice?

7.7 Parallel masked inference↩︎

To avoid spending one network evaluation per coordinate, practical masked samplers often reveal a whole block of coordinates in a single reverse step. This can mean two different things. The simplest implementation is factorized block decoding: if \(A\) is the set of coordinates to reveal from the current partially observed sequence \(x_k\), it samples \[a_A\sim \prod_{i\in A}\mathsf{D}_k(i,a_i\mid x_k).\] The exact object, however, is the joint conditional law \[\mathbb{P}(X_{0,A}=a_A\mid X_k=x_k).\] Equivalently, for any ordering \(A=\{i_1,\ldots,i_m\}\), the chain rule writes this joint conditional as \[\prod_{\ell=1}^{m} \mathbb{P}\!\left( X_{0,i_\ell}=a_{i_\ell} \mid X_k=x_k,\, X_{0,i_1}=a_{i_1},\ldots,X_{0,i_{\ell-1}}=a_{i_{\ell-1}} \right).\] Thus a block can be sampled exactly by sequentially querying updated conditionals, but the product of one-coordinate posteriors from the same context is exact only under conditional independence. The gap between these two laws is the factorization error. The inference schedule therefore sets a speed–accuracy tradeoff: larger blocks shorten the reverse chain but introduce conditional-independence bias, while singleton blocks remove the bias at the cost of as many reverse evaluations as there are coordinates. Lavenant and Zanella [37] analyze this error directly for masked diffusions with factorized approximations, decomposing the relative-entropy error into learning and factorization terms and relating the optimal block-size schedule to an information profile of the data distribution.

This is not the only possible approach of parallel inference. The block sampler above is a factorized approximation to the joint conditional; a different question is whether one can parallelize the exact sequential conditional-sampling procedure itself. In an oracle model with access to exact conditional marginals of a target law on \(\mathcal{V}^L\), Anari, Gao, and Rubinstein [38] show how to organize the queries to sample arbitrary product-space laws in \(\widetilde{O}(L^{2/3})\) parallel time. Anari et al. [39] improve the parallel time to \(\widetilde{O}(L^{1/2})\) using autospeculative rejection sampling. These results are best read as a corrected-parallel counterpart to the simple block update: they do not replace a joint conditional by independent marginals, but instead use parallel queries to simulate the sequential conditionals with a correction.

7.8 CTMC setup and learning↩︎

The discussion so far used a discrete-time masked chain because it makes the Bayes rule and the denoising posterior easy to see. For analyzing practical samplers, however, the sharper language is continuous time: it lets us describe the exact reverse dynamics as one coordinate jump at a time, and then view parallel updates as a numerical approximation to those dynamics.

This is one reason that most modern discrete diffusion models are phrased not with the stochastic matrices \(P_k^{\to}\) but as a continuous-time Markov chain (CTMC), the finite-state analogue of the forward SDE. The forward corruption is run as a CTMC whose rate matrix factorizes over coordinates, so each coordinate is noised independently. This has a decisive consequence for the reverse process: at any instant only one coordinate changes, because the probability that two independent coordinates jump in the same infinitesimal interval is \(o(\,\mathrm{d}t)\). The exact reverse dynamics are therefore unambiguous and update a single coordinate at a time, and the parallel-update bias described above reappears in a controlled form, as the error of approximating these dynamics on a finite time grid (\(\tau\)-leaping). Continuous time thus separates the posterior/score error from the discretization error cleanly, produces the cleanest analogue of the score, and admits a flexible family of reverse samplers [40], [41].

A CTMC on a finite set \(\mathcal{X}\) is generated by a time-dependent forward rate matrix \(R_t^{\to}\) with \[R_t^{\to}(x,y)\geq0\;(y\neq x), \qquad R_t^{\to}(x,x)=-\sum_{y\neq x}R_t^{\to}(x,y).\] Thus the off-diagonal entries are jump rates, and the diagonal entry makes each row sum to zero. Over an infinitesimal step, \[\mathbb{P}(X_{t+\,\mathrm{d}t}=y\mid X_t=x) = \begin{cases} R_t^{\to}(x,y)\,\mathrm{d}t+o(\,\mathrm{d}t), & y\neq x,\\ 1+R_t^{\to}(x,x)\,\mathrm{d}t+o(\,\mathrm{d}t), & y=x. \end{cases}\] Equivalently, \[\mathbb{P}(X_{t+\,\mathrm{d}t}=y\mid X_t=x) = \delta_{xy}+R_t^{\to}(x,y)\,\mathrm{d}t+o(\,\mathrm{d}t).\] If \(p_t(x)=\mathbb{P}(X_t=x)\) is written as a row vector, then its evolution is \[\label{eq:ctmc-forward-equation} \partial_t p_t=p_tR_t^{\to}, \qquad \partial_t p_t(y)=\sum_{x\in\mathcal{X}}p_t(x)R_t^{\to}(x,y).\tag{60}\] This is the finite-state counterpart of the Fokker–Planck equation 2 : a rate matrix replaces the drift-and-diffusion operator acting on continuous densities.

The rate matrix is also where the corruption family is specified. An absorbing (masked) generator sends each token to \([\mathsf M]\) at a schedule-dependent rate and is the continuous-time version of the masked diffusion of the previous subsections; a uniform generator pushes each token toward the uniform distribution over \(\mathcal{V}\); other choices encode structured token graphs. Masked diffusion is therefore a special case of this CTMC picture.

The reverse process is again a CTMC, and its rates are fixed by a Bayes-ratio formula that is the time-continuous version of the reverse kernel 54 . For \(y\neq x\), \[\label{eq:ctmc-reverse-rate} R_t^{\leftarrow}(y,x) = R_t^{\to}(x,y)\,\frac{p_t(x)}{p_t(y)}.\tag{61}\] With the diagonal chosen so that rows sum to zero, these rates are exactly the ones that track the same marginal in reverse time. Indeed, for each \(y\) with \(p_t(y)>0\), \[\begin{align} (p_tR_t^{\leftarrow})(y) &= \sum_{x\neq y}p_t(x)R_t^{\leftarrow}(x,y) - p_t(y)\sum_{x\neq y}R_t^{\leftarrow}(y,x) \\ &= \sum_{x\neq y}p_t(x)R_t^{\to}(y,x)\frac{p_t(y)}{p_t(x)} - p_t(y)\sum_{x\neq y}R_t^{\to}(x,y)\frac{p_t(x)}{p_t(y)} \\ &= p_t(y)\sum_{x\neq y}R_t^{\to}(y,x) - \sum_{x\neq y}p_t(x)R_t^{\to}(x,y) \\ &= -\partial_t p_t(y), \end{align}\] where the last equality is the forward equation 60 . Thus \(-\partial_t p_t=p_tR_t^{\leftarrow}\), which is the forward equation for the chain when it is simulated from time \(T\) down to time \(0\). Starting the reverse CTMC from \(p_T\) therefore gives marginal \(p_t\) at every intermediate time. As in the continuous case, the forward rates \(R_t^{\to}\) are known by design and the only unknown ingredient is the collection of same-time density ratios over neighboring states. We write the exact score as \[\mathsf{s}^{\star}_t(x,y):=\frac{p_t(y)}{p_t(x)}, \qquad y\neq x,\;R_t^{\to}(y,x)>0.\] This is the continuous-time form of the ratio score 55 , the discrete analogue of \(\nabla\log p_t\) that records the finite relative change of the noised law along each admissible move. Let \(\mathsf{s}_t(x,y)>0\) be the learned prediction, \(\mathsf{s}_t(x,y)\approx\mathsf{s}^{\star}_t(x,y)\). The orientation is chosen for reverse sampling: when the reverse chain is at \(x\), a candidate predecessor \(y\) is admissible exactly when the forward CTMC can jump from \(y\) to \(x\). With this convention the learned reverse rate is \[\widehat{R}_t^{\leftarrow}(x,y) = R_t^{\to}(y,x)\,\mathsf{s}_t(x,y), \qquad y\neq x.\] At first the training target looks inaccessible, since the ratio contains the unknown marginal law \(p_t\). The key denoising identity is that this unknown ratio can be written as a posterior average of known forward likelihood ratios. \[\label{eq:discrete-ratio-denoising-identity} \mathsf{s}^{\star}_t(x,y) = \mathbb{E}\!\left[ \frac{\mathbb{P}(X_t=y\mid X_0)}{\mathbb{P}(X_t=x\mid X_0)} \,\middle|\, X_t=x \right].\tag{62}\] Thus training can sample clean data \(X_0\), corrupt it to \(X_t=x\), and use the known forward likelihood ratio inside the loss.

Score-entropy training [41] turns this identity into a supervised loss for positive ratios by first choosing a scalar discrepancy for one reverse edge. Let \(\mathsf{s}^{\star}>0\) denote the exact ratio, let \(\mathsf{s}>0\) be the learned prediction, and set \(F(u)=-\log u\) on \(u>0\). Score entropy uses the scaled Bregman divergence \[\mathsf{s}^{\star}D_F(\mathsf{s},\mathsf{s}^{\star}) = \mathsf{s}^{\star}\bigl(F(\mathsf{s})-F(\mathsf{s}^{\star})-F'(\mathsf{s}^{\star})(\mathsf{s}-\mathsf{s}^{\star})\bigr) = \mathsf{s}-\mathsf{s}^{\star}\log\mathsf{s}+\mathsf{s}^{\star}\log\mathsf{s}^{\star}-\mathsf{s}^{\star}.\] This quantity is nonnegative, has derivative \(1-\mathsf{s}^{\star}/\mathsf{s}\), and is minimized at \(\mathsf{s}=\mathsf{s}^{\star}\); note that it is defined only for positive ratios, so the learned \(\mathsf{s}_t(x,y)\) must be constrained to remain positive. In the learning objective the terms \(\mathsf{s}^{\star}\log\mathsf{s}^{\star}-\mathsf{s}^{\star}\) are independent of the learned score. Thus the model-dependent part of one edge is \(\mathsf{s}_t(x,y)-\mathsf{s}^{\star}_t(x,y)\log \mathsf{s}_t(x,y)\). In training, the unknown ratio \(\mathsf{s}^{\star}_t(x,y)\) is replaced by the computable forward likelihood ratio from the noising process. A schematic denoising score-entropy loss is \[\label{eq:score-entropy-loss} \mathbb{E}_{t,X_0} \sum_{x\in\mathcal{X}}\mathbb{P}(X_t=x\mid X_0) \sum_{y:\,R_t^{\to}(y,x)>0} w_t(x,y) \left[ \mathsf{s}_t(x,y) - \frac{\mathbb{P}(X_t=y\mid X_0)}{\mathbb{P}(X_t=x\mid X_0)} \log \mathsf{s}_t(x,y) \right],\tag{63}\] up to terms independent of the learned score. Here \(w_t(x,y)\geq0\) is a chosen weight, often the incoming forward rate \(R_t^{\to}(y,x)\) or a multiple of it.

The reason this objective has the correct population target is exactly the denoising identity 62 above. Since the displayed loss is affine in the likelihood ratio appearing in 62 , averaging over the unknown clean data simply replaces it by its conditional mean. The resulting edgewise risk is, up to constants, \(\mathsf{s}_t(x,y) - \mathsf{s}^{\star}_t(x,y)\log \mathsf{s}_t(x,y),\) and is minimized at the exact score \(\mathsf{s}_t(x,y)=\mathsf{s}^{\star}_t(x,y)\) on each reverse-admissible edge.

The special choice \(w_t(x,y)=R_t^{\to}(y,x)\) gives this denoising loss its sampling interpretation. On the edge from the current state \(x\) back to a possible predecessor \(y\), the true reverse rate is \(R_t^{\to}(y,x)\mathsf{s}^{\star}_t(x,y)\), while the learned reverse rate is \(R_t^{\to}(y,x)\mathsf{s}_t(x,y)\). The corresponding per-edge contribution to the CTMC path-space relative entropy is (as will be shown in 64 ) \[R_t^{\to}(y,x)\left[ \mathsf{s}_t(x,y) - \mathsf{s}^{\star}_t(x,y) + \mathsf{s}^{\star}_t(x,y) \log \frac{\mathsf{s}^{\star}_t(x,y)}{\mathsf{s}_t(x,y)} \right].\] After dropping terms independent of the learned score, this is exactly the objective in 63 with \(w_t(x,y)=R_t^{\to}(y,x)\). Thus the rate-weighted choice matches the path-space error quantity used in the CTMC sampling analysis below.

7.9 CTMC sampling↩︎

Once the learned ratios have specified the reverse rates, the remaining problem is numerical simulation of the resulting jump process. \(\tau\)-leaping is a standard acceleration of the stochastic simulation algorithm for chemical reaction networks, introduced by Gillespie [42]. The idea is to group many small jumps over a short interval, while pretending that the jump rates are essentially constant during that interval.

For the learned reverse CTMC, keep the same forward-time grid convention as before, \[0=t_0<t_1<\cdots<t_K=T, \qquad h_k=t_{k+1}-t_k.\] The reverse sampler moves from \(t_{k+1}\) down to \(t_k\). If the current state is \(x\) at time \(t_{k+1}\), the simplest tau-leap freezes the learned reverse rates at \((t_{k+1},x)\) and uses the first-order kernel \[\mathbb{P}(\widehat{X}_{t_k}=z\mid \widehat{X}_{t_{k+1}}=x) = h_k\,\widehat{R}_{t_{k+1}}^{\leftarrow}(x,z) +O(h_k^2), \qquad z\neq x,\] with the remaining probability assigned to staying at \(x\). Equivalently, for product spaces one draws independent Poisson clocks \(N_z\sim\mathsf{Poisson}(h_k\widehat{R}_{t_{k+1}}^{\leftarrow}(x,z))\) for the admissible local moves and applies the proposed coordinate changes in parallel, with a fixed tie-breaking convention if several incompatible clocks ring. This is the discrete analogue of an Euler–Maruyama step: it replaces the time-varying reverse generator by a frozen one and, in parallel implementations, introduces a local independence error by allowing several coordinates to update during one step. We refer to [3] for a more detailed discussion of \(\tau\)-leaping and related algorithms.

Tau-leaping is convenient, but it is still an approximation because the rates are frozen over a whole interval. A closely related exact construction is uniformization, originally introduced by Grassmann [43], applied to chemical reaction simulation by Beentjes and Baker [44], and applied to discrete diffusion models by Chen and Ying [45].

Let us describe how uniformization works. On the interval \([t_k,t_{k+1}]\), choose a clock rate \[\Lambda_k \geq \sup_{\substack{s\in[t_k,t_{k+1}]\\ x\in\mathcal{X}}} \sum_{z\neq x}\widehat{R}_s^{\leftarrow}(x,z).\] Starting from the state at time \(t_{k+1}\), draw Poisson event times in \([t_k,t_{k+1}]\) with rate \(\Lambda_k\) and process them in decreasing time. At an event time \(s\), if the current state is \(x\), jump to \(z\neq x\) with probability \(\widehat{R}_s^{\leftarrow}(x,z)/\Lambda_k\) and otherwise stay at \(x\). The stay-put events are virtual jumps. Since the clock dominates all outgoing rates, this thinning construction has exactly the jump law of the learned reverse CTMC on the interval, without a frozen-rate approximation.

7.10 Error analysis for CTMC samplers↩︎

Having described the two CTMC simulation schemes, we now ask how their output laws differ from the desired data law. The numerical analysis can be organized in the same way as in Section 6. There are three contributions. First, the exact reverse process should start from \(p_T\), whereas the sampler is initialized from some implementable law \(\widehat{p}_T\). This gives the common initialization gap. Define \[\varepsilon_{\mathsf{init}}^2:=\mathop{\mathrm{D_{\mathsf{KL}}}}(p_T\|\widehat{p}_T).\] Two common choices of \(\widehat{p}_T\) should be interpreted differently. If \(\widehat{p}_T\) is a point mass, then \(\varepsilon_{\mathsf{init}}^2\) is mainly a support condition. In absorbing masked diffusion, for instance, one often takes \(\widehat{p}_T=\delta_{[\mathsf M]^L}\); the KL is finite only if the forward terminal law \(p_T\) is also supported on the all-mask state. With a finite integrated masking rate this may fail, so one should either force exact terminal absorption or measure the endpoint error in a weaker metric.

If \(\widehat{p}_T\) is full-support, the same term is an ordinary initialization error. For example, if \(\widehat{p}_T\) is uniform and the forward CTMC contracts KL to \(\widehat{p}_T\) at rate \(\rho\), then \[\varepsilon_{\mathsf{init}}^2 \leq e^{-\rho T}\mathop{\mathrm{D_{\mathsf{KL}}}}(p_0\|\widehat{p}_T) \leq e^{-\rho T}\log|\mathcal{X}|.\] On product token spaces \(\mathcal{X}=\mathcal{V}^L\), this term is of order \(L(\log|\mathcal{V}|)e^{-T}\) for the usual independent-coordinate corruption.

After initialization, fix a simulator and first imagine running that simulator with the true reverse rates, and only afterwards replace the true ratio score by the learned ratio score: \[p_0 \quad\longrightarrow\quad p_0^{\mathsf{num},*} \quad\longrightarrow\quad \widehat{p}_0^{\mathsf{num}}.\] Here \(p_0^{\mathsf{num},*}\) denotes the output law of the chosen numerical simulator when it is run with the true reverse rates. The first comparison isolates numerical error, while the second comparison isolates score-estimation error. Strictly speaking, KL has no triangle inequality, so the rigorous proof usually performs this split on path space, or directly in the log-intensity integrand, rather than by adding two marginal KL divergences. Nevertheless, this bookkeeping is useful for organizing the error analysis.

The score-estimation term has a path-space form analogous to Girsanov’s theorem for diffusions. We write it in the same epsilon-squared convention as the continuous score error: \[\begin{align} \label{eq:ctmc-score-entropy-error} \varepsilon_{\mathsf{SE}}^2 := \int_0^T \mathbb{E}_{X_t\sim p_t} \sum_{y:\,R_t^{\to}(y,X_t)>0} R_t^{\to}(y,X_t) \left[ \mathsf{s}_t(X_t,y) -\mathsf{s}^{\star}_t(X_t,y) +\mathsf{s}^{\star}_t(X_t,y) \log \frac{\mathsf{s}^{\star}_t(X_t,y)}{\mathsf{s}_t(X_t,y)} \right]\,\mathrm{d}t. \end{align}\tag{64}\] This is the same score-entropy Bregman loss that appears in the learning objective 63 , which explains why it is a natural training loss.

For \(\tau\)-leaping, the exact-score numerical law is not the true reverse CTMC law. On each interval \([t_k,t_{k+1}]\) it replaces the time-dependent true reverse rates by the frozen rates \(R_{t_{k+1}}^{\leftarrow}\). Replacing these frozen true rates by \(\widehat{R}_{t_{k+1}}^{\leftarrow}\) then adds the score-estimation part. Let \(\bar h\) denote the mesh parameter controlling the reverse step sizes, and let \(\bar R\) be an upper bound on the total outgoing reverse rate, \(\sum_{z\neq x}R_s^{\leftarrow}(x,z)\leq\bar R\) on \([0,T]\). A stochastic-integral analysis of the combined comparison gives, under bounded-rate and regularity assumptions, a bound of the schematic form \[\label{eq:tau-leaping-kl-general} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left(p_0\middle\|\widehat{p}_0^\tau\right) \lesssim \varepsilon_{\mathsf{init}}^2 +\varepsilon_{\mathsf{SE}}^2 +\bar R^{\,2}\bar h T.\tag{65}\] The three terms are the initialization error, the score-estimation error, and the exact-score tau-leaping discretization error.

For product token spaces \(\mathcal{X}=\mathcal{V}^L\), more recent work gives a sharper version of the same message with explicit dependence on the vocabulary size. If the learned ratios are clipped to \([B^{-1},B]\) and the reverse grid is controlled by mesh parameter \(\bar h\), then the standard tau-leaping sampler obeys, suppressing logarithmic factors and endpoint regularity constants, \[\label{eq:tau-leaping-kl-product} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left(p_0\middle\|\widehat{p}_0^\tau\right) \lesssim \varepsilon_{\mathsf{init}}^2 +\varepsilon_{\mathsf{SE}}^2 +\bar h\, L^2|\mathcal{V}|\,T,\tag{66}\] as in the analysis of Liang et al. [46]. For the usual independent-coordinate corruption one may substitute \(\varepsilon_{\mathsf{init}}^2\lesssim L(\log|\mathcal{V}|)e^{-T}\). Thus making the tau-leaping discretization error of order \(\varepsilon\) requires \(\bar h\) of order \(\widetilde{O}(\varepsilon/(L^2|\mathcal{V}|T))\), and hence a deterministic grid with roughly \(\widetilde{O}(L^2|\mathcal{V}|/\varepsilon)\) reverse steps in this analysis.

We can achieve high-accuracy sampling using uniformization. With the true ratio score, it simulates the true reverse CTMC exactly, so \(p_0^{\mathsf{uni},*}=p_0\) if the reverse process is initialized from \(p_T\). Thus, for the learned score, the error estimate of uniformization reduces to \[\label{eq:ctmc-exact-kl} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left(p_0\middle\|\widehat{p}_0^{\rm uni}\right) \leq \varepsilon_{\mathsf{init}}^2 +\varepsilon_{\mathsf{SE}}^2.\tag{67}\] On product token spaces \(\mathcal{X}=\mathcal{V}^L\) with the usual independent-coordinate corruption, one may substitute \(\varepsilon_{\mathsf{init}}^2\lesssim L(\log|\mathcal{V}|)e^{-T}\) in 67 , giving \(\mathop{\mathrm{D_{\mathsf{KL}}}}(p_0\|\widehat{p}_0^{\rm uni})\lesssim L(\log|\mathcal{V}|)e^{-T}+\varepsilon_{\mathsf{SE}}^2\). Choosing \(T\asymp\log(L\log|\mathcal{V}|/\varepsilon^2)\) and learning the ratio score so that \(\varepsilon_{\mathsf{SE}}^2\lesssim\varepsilon^2\) gives \(\mathop{\mathrm{D_{\mathsf{KL}}}}(p_0\|\widehat{p}_0^{\rm uni})\lesssim\varepsilon^2\). Chen and Ying [45] further show that, with adaptive dominating rates, the expected number of Poisson events is nearly linear in the hypercube dimension. Their theorem is stated for \(\{0,1\}^d\); for an alphabet of size \(|\mathcal{V}|\), a fixed binary encoding uses \(q=\lceil\log_2|\mathcal{V}|\rceil\) bits per token, so \(d=Lq\). In token notation this gives an expected event count \(\widetilde{O}\!\left(L\log|\mathcal{V}|\right)\) for \(\varepsilon^2\) KL error.

8 Guidance, Reward Tilting, and Inference-Time RL↩︎

So far we have built diffusion samplers, continuous and discrete, that aim to reproduce the data distribution, or a slightly smoothed version of it. In many applications one instead wants a controlled modification of that base sampler: generate samples satisfying a condition, prefer high-reward outputs, or adapt a model at inference time. Guidance, reward tilting, and inference-time RL are all mechanisms for biasing the base model toward preferred outputs while staying close to the pretrained distribution, and their common mathematical language is the KL-regularized change of measure developed in this section.

8.1 Reward-tilted targets and KL-regularized optimization↩︎

Let \(p_0\) denote the clean-output distribution produced by the pretrained sampler; we keep the diffusion convention that clean samples are written as \(x_0\). At this point we specify only the desired law on clean outputs; the corresponding reverse-time dynamics will be derived below. Given a reward \(r:\mathbb{R}^d\to\mathbb{R}\) and inverse temperature \(\beta\geq0\), define the exponentially tilted target \[\label{eq:reward-tilt} p_0^\beta(x_0) = \frac{1}{Z_\beta}p_0(x_0)e^{\beta r(x_0)}.\tag{68}\] Here \[Z_\beta=\int p_0(x_0)e^{\beta r(x_0)}\,\mathrm{d}x_0 = \mathbb{E}_{X_0\sim p_0}e^{\beta r(X_0)}\] is the normalizing constant, assumed finite. Equivalently, for any test function \(f\), \[\mathbb{E}_{p_0^\beta}f(X_0) = \frac{\mathbb{E}_{p_0}\!\left[f(X_0)e^{\beta r(X_0)}\right]}{\mathbb{E}_{p_0} e^{\beta r(X_0)}}.\] Thus the target is obtained from the base model by reweighting samples according to their clean-output reward. The case \(\beta=0\) recovers \(p_0\), while larger \(\beta\) places more mass on high-reward regions and therefore trades diversity under the base model for reward improvement.

Conditional generation fits the same form. If \(y\) is an observation, class label, or prompt and \(p(y\mid x_0)\) is the corresponding likelihood or compatibility model, then Bayes’ rule gives \[p_0(x_0\mid y)\propto p_0(x_0)p(y\mid x_0).\] This is an exponential tilt with reward function \(x_0\mapsto\log p(y\mid x_0)\) and \(\beta=1\); a classifier score or learned preference model plays the same role when an explicit likelihood is unavailable. The inference-time problem is to modify the reverse sampler so that its clean-output law approximates such a tilted or conditional target while still reusing the pretrained base dynamics.

The same tilted law has a useful variational meaning. Among all distributions \(q\) over final samples, it is the optimizer of \[\label{eq:gibbs-variational} \sup_q \left\{ \beta\mathbb{E}_q r-\mathop{\mathrm{D_{\mathsf{KL}}}}(q\|p_0) \right\}.\tag{69}\] Equivalently, if \(\beta>0\), \[p_0^\beta = \mathop{\mathrm{arg\,max}}_q \left\{ \mathbb{E}_q r-\frac{1}{\beta}\mathop{\mathrm{D_{\mathsf{KL}}}}(q\|p_0) \right\}.\] Thus \(\beta\) controls the tradeoff between reward improvement and staying close to the base model. Large \(\beta\) pushes hard toward high reward and risks mode collapse or reward hacking; small \(\beta\) preserves the base distribution but gives weaker alignment.

8.2 Guidance as score tilting↩︎

To run a diffusion sampler for the tilted target \(p_0^\beta\) we need the score of its noised marginals, so the question is how that score differs from the base score we have already learned. Let \(p_t\) and \(p_t^\beta\) be the noised marginals of \(p_0\) and \(p_0^\beta\). Write \(P_{0,t}^{\to}(x_0\to x)\) for the density of the clean-to-noisy marginal forward kernel \(\mathop{\mathrm{Law}}(X_t\mid X_0=x_0)\). Because the reward acts only on the clean sample \(x_0\), noising the tilted law gives \[p_t^\beta(x) = \frac{1}{Z_\beta}\int P_{0,t}^{\to}(x_0\to x)\, e^{\beta r(x_0)}p_0(x_0)\,\mathrm{d}x_0 = \frac{p_t(x)}{Z_\beta}\,h_t^\beta(x), \qquad h_t^\beta(x):=\mathbb{E}\!\left[e^{\beta r(X_0)}\mid X_t=x\right].\] The noised tilted law is thus the base noised law reweighted by the posterior tilt factor \(h_t^\beta\), which reads off, from the current noisy state, the expected exponential reward of the clean sample it will denoise to. Because this marginal is a product, its log-gradient splits into the base score plus a correction, \[\nabla\log p_t^\beta(x) = \nabla\log p_t(x)+\nabla\log h_t^\beta(x),\] so exact guidance merely adds the tilt gradient \(\nabla\log h_t^\beta\) to the base score, and practical methods differ only in how they approximate it.

For conditional generation, the reward function is \(x_0\mapsto\log p(y\mid x_0)\). The same posterior-tilt notation gives the noisy likelihood \[h_t(x\mid y) := \mathbb{E}[p(y\mid X_0)\mid X_t=x].\] This is the quantity a noisy classifier estimates. It is not the noised conditional density itself; rather, it tilts the base noisy density. If the condition is fixed, the noised conditional density is \[p_t(x\mid y)\propto p_t(x)h_t(x\mid y),\] where the missing normalizing constant is independent of \(x\). Taking a log-gradient at noise level \(t\) gives \[\label{eq:classifier-guidance} \nabla\log p_t(x\mid y) = \nabla\log p_t(x) + \nabla\log h_t(x\mid y).\tag{70}\] This is classifier guidance [47]: train or use a classifier on noisy inputs, then add its gradient to the unconditional score. Classifier-free guidance [48] estimates the same increment without a separate classifier. We use the already established score notation in the form \[\mathsf{s}_t(x)\approx\nabla\log p_t(x), \qquad \mathsf{s}_t(x\mid y)\approx\nabla\log p_t(x\mid y).\] Thus 70 suggests \[\mathsf{s}_t(x\mid y)-\mathsf{s}_t(x) \approx \nabla\log h_t(x\mid y).\] With guidance strength \(\beta\geq0\), classifier-free guidance replaces the unconditional score by \[\label{eq:cfg} \mathsf{s}_t(x) \quad\longmapsto\quad \mathsf{s}_t(x) + \beta\bigl(\mathsf{s}_t(x\mid y)-\mathsf{s}_t(x)\bigr).\tag{71}\] The choice \(\beta=1\) gives the conditional score estimate \(\mathsf{s}_t(x\mid y)\), while \(\beta>1\) extrapolates the conditional-score increment; empirically this often improves condition satisfaction, at the cost of moving farther from the base distribution. If the two learned scores were exact and \(\beta=1\), the approximation 71 would recover the conditional score in 70 . For \(\beta\neq1\), it instead gives the score of the noise-level power tilt \(p_t(x)h_t(x\mid y)^\beta\), whereas the clean reward tilt 68 would involve the posterior factor \(h_t^\beta(x)=\mathbb{E}[p(y\mid X_0)^\beta\mid X_t=x]\). Powering the noisy likelihood is not the same as tilting by the noisy expectation of the powered clean likelihood, so the two coincide only at \(\beta=1\).

For small tilts, the posterior tilt factor has a linear approximation. Since \(h_t^\beta(x)=1+\beta V_t(x)+O(\beta^2)\) with \(V_t(x)=\mathbb{E}[r(X_0)\mid X_t=x]\), we have \[\label{eq:value-guidance} \nabla\log p_t^\beta(x) \approx \nabla\log p_t(x) + \beta\nabla V_t(x), \qquad V_t(x)=\mathbb{E}[r(X_0)\mid X_t=x].\tag{72}\] In the small-tilt regime, then, guidance perturbs the score by the gradient of the posterior expected reward \(V_t\), so a reward model need only capture this first-order landscape. The simplification is that \(V_t\) averages the reward \(r\) itself, whereas the exact factor \(h_t^\beta\) averages \(e^{\beta r}\); the two agree only to first order in \(\beta\). Exact guidance uses the full \(h_t^\beta\), and the small-tilt form trades it for a cheaper reward-gradient correction.

8.3 Reward tilting as a Polchinski flow↩︎

Stepping back from the practical methods, the tilt factor \(h_t^\beta\) ties guidance to the renormalization viewpoint of Section 4. Recall the variance-exploding normalization \(X_t=X_0+\sqrt t\,Z\) used there, in which \(p_t=p_0*\mathcal{N}\!\left(0,tI\right)\) and the effective potential \(U_t=-\log p_t\) solves the finite-dimensional Polchinski equation 30 . Reading the integral that defines \(h_t^\beta\) the other way, the unnormalized tilted marginal \(h_t^\beta\,p_t=(e^{\beta r}p_0)*\mathcal{N}\!\left(0,tI\right)\) is the same forward channel applied to the reward-reweighted data measure \(e^{\beta r}p_0\)—noising the tilted data and tilting the noised data by \(h_t^\beta\) coincide—so it is a heat flow and solves 29 just as \(p_t\) does. Dividing out \(p_t\) leaves a forward Kolmogorov equation for the tilt factor, \[\label{eq:tilt-pde} \partial_t h_t^\beta = \tfrac12\Delta h_t^\beta + \mathsf{s}^{\star}_t\cdot\nabla h_t^\beta,\tag{73}\] so \(h_t^\beta\) is transported by the score-driven generator \(\tfrac12\Delta+\mathsf{s}^{\star}_t\cdot\nabla\). Equivalently, the tilted effective potential \(U_t^\beta:=-\log(h_t^\beta\,p_t)\) obeys the same Polchinski equation, \[\label{eq:polchinski-tilt} \partial_t U_t^\beta = \tfrac12\Delta U_t^\beta - \tfrac12\left\lVert \nabla U_t^\beta \right\rVert^2,\tag{74}\] only with the bare potential shifted from \(U_0=-\log p_0\) to \(U_0^{\beta} = -\log p_0-\beta r\). Reward tilting therefore runs the renormalization flow from a reward-shifted initial potential, and the guided score \(-\nabla U_t^\beta=\mathsf{s}^{\star}_t+\nabla\log h_t^\beta\) is exactly the additive guidance correction of the previous subsection.

It should be emphasized that this is a reformulation, not a recipe. Solving the tilt equation 73 or the Polchinski equation 74 is no easier than evaluating the posterior expectation \(h_t^\beta(x)=\mathbb{E}[e^{\beta r(X_0)}\mid X_t=x]\) that defines it: the PDE carries the same computational difficulty as the conditional expectation.

8.4 Path-space control and the Doob transform↩︎

The score-tilting, classifier-guidance, and Polchinski-flow descriptions of the previous subsections are all exact, but they all rest on the same intractable object, the posterior tilt \(h_t^\beta\). To see how it is approximated in practice it helps to pass from marginals to whole reverse trajectories: on path space the tilt becomes a stochastic control problem, and its exact solution—the Doob transform of the base reverse chain—exposes precisely which conditional expectations must be estimated.

Let us write a continuous reverse sampler in generation time as \[\,\mathrm{d}Y^{\leftarrow}_s=b_s(Y^{\leftarrow}_s)\,\mathrm{d}s+\sigma_s\,\mathrm{d}B^{\leftarrow}_s, \qquad 0\leq s\leq T,\] and let \(\mathbb{P}^0\) be its path law on full reverse trajectories \((Y^{\leftarrow}_s)_{0\leq s\leq T}\), with clean output \(Y^{\leftarrow}_T\).

Reward-tilted sampling trades terminal reward against the cost of departing from this base law. At the level of path laws it is the KL-regularized optimization \[\label{eq:path-gibbs} \sup_{\mathbb{Q}} \left\{ \beta\mathbb{E}_{\mathbb{Q}}r(Y^{\leftarrow}_T)-\mathop{\mathrm{D_{\mathsf{KL}}}}(\mathbb{Q}\|\mathbb{P}^0) \right\} = \log\mathbb{E}_{\mathbb{P}^0}e^{\beta r(Y^{\leftarrow}_T)},\tag{75}\] the path-space version of the Gibbs variational principle 69 . Its optimizer is the Gibbs path law \[\label{eq:path-gibbs-optimizer} \frac{\,\mathrm{d}\mathbb{Q}^\star}{\,\mathrm{d}\mathbb{P}^0} \!\left((Y^{\leftarrow}_s)_{0\leq s\leq T}\right) = \frac{e^{\beta r(Y^{\leftarrow}_T)}}{\mathbb{E}_{\mathbb{P}^0}e^{\beta r(Y^{\leftarrow}_T)}}.\tag{76}\]

To turn this into a sampler we express it as a control problem on the drift. A controlled sampler changes only the drift, \[\,\mathrm{d}Y^{\leftarrow}_s=\bigl(b_s(Y^{\leftarrow}_s)+\sigma_s u_s(Y^{\leftarrow}_s)\bigr)\,\mathrm{d}s+\sigma_s\,\mathrm{d}B^{\leftarrow}_s,\] with path law \(\mathbb{P}^u\); using the same diffusion coefficient means the two path laws differ only through the drift. Under the usual absolute-continuity and integrability assumptions, the Girsanov KL formula in Appendix 9 gives \[\label{eq:girsanov-control-cost} \mathop{\mathrm{D_{\mathsf{KL}}}}(\mathbb{P}^u\|\mathbb{P}^0) = \frac{1}{2} \mathbb{E}_{\mathbb{P}^u}\int_0^T\left\lVert u_s(Y^{\leftarrow}_s) \right\rVert^2\,\mathrm{d}s,\tag{77}\] with \(u_s(Y^{\leftarrow}_s)\) a progressively measurable control, so the quadratic control energy is exactly the KL cost of steering the sampler away from its base law. Restricting \(\mathbb{Q}\) to controlled diffusions \(\mathbb{P}^u\) and substituting this cost into 75 makes the Gibbs problem the equivalent stochastic control problem over drifts \[\label{eq:control-objective} \sup_u \mathbb{E}_{\mathbb{P}^u}\left[ \beta r(Y^{\leftarrow}_T) - \frac{1}{2}\int_0^T\left\lVert u_s(Y^{\leftarrow}_s) \right\rVert^2\,\mathrm{d}s \right].\tag{78}\] The two are the same optimization stated at the level of measures and of drifts.

Reading the optimizer 76 back at the level of drifts recovers the slogan “guidance is control”: the optimal control’s drift increment is exactly the value-gradient term of the ideal guided reverse dynamics of Subsection 8.2, now in this generation-time parametrization. The KL cost 77 , which measured discretization error in Section 6, is here a deliberately chosen guidance budget, spent where it raises expected terminal reward most efficiently.

From here to the end of the section we adopt the discrete-time point of view, working with the time-discretized reverse chain \((X_K,\ldots,X_0)\) of Section 5 rather than the continuous SDE. The reason is that the objects we now build—the Doob transform below, the sequential Monte Carlo weights of Subsection 8.5, and the policies of Subsection 8.6—all act step by step, and a finite chain of reverse kernels states them without stochastic-calculus bookkeeping; the continuous formulas above are recovered in the small-step limit. In this notation the clean output is indexed by \(0\), so a terminal reward is written as \(r(X_0)\). The path-space optimizer has an explicit Markov-kernel form: each base reverse proposal is reweighted by the tilt factor of its continuations and normalized by the tilt factor at the current state. This kernel reweighting is the Doob transform of the base reverse chain.

The value function \(h_k\) is intractable, and the two ways of coping with this organize the rest of the section. One route uses an approximate value only as a proposal and removes its bias by reweighting, giving an exact sampler in the limit (Subsection 8.5); the other approximates the value and accepts the resulting bias, the inference-time-RL view (Subsection 8.6).

8.5 Feynman–Kac correction and sequential Monte Carlo↩︎

Theorem 6 gives the exact guided sampler, but its kernels require the value functions \(h_k(x_k)=\mathbb{E}_{\mathbb{P}^0}[e^{\beta r(X_0)}\mid X_k=x_k]\), which are exactly as unavailable as the posterior reward-to-go in the score picture. This subsection takes the first route of the fork above: use an approximate value only as a proposal and remove its bias by reweighting. That is the Feynman–Kac / sequential Monte Carlo route, and it is the mechanism behind inference-time scaling of guided samplers.

The name comes from the value function. Along the base reverse process, \(h\) is a conditional expectation of a terminal functional, hence a martingale, \[h_{k+1}(x_{k+1})=\int h_k(y)\, P_k^{0,\leftarrow}(x_{k+1}\to \,\mathrm{d}y),\] the same identity that appeared in Theorem 6. In continuous time it solves the backward equation \(\partial_s h_s+\mathcal{L}_s h_s=0\) with terminal data \(h=e^{\beta r}\), where \(\mathcal{L}_s\) is the generator of the base reverse diffusion; this is the reverse-time companion of the noising-time equation 73 for the same tilt factor in Subsection 8.3. At the clean end \(k=0\) the state is \(X_0\) itself, so the terminal value is exactly \(h_0(x_0)=e^{\beta r(x_0)}\).

The key fact is that the entire telescoping of the Doob transform collapses to a single terminal weight.

This suggests the simplest exact algorithm: draw \(N\) trajectories from the base model, weight each by \(e^{\beta r(x_0)}\), and resample. As \(N\to\infty\) the weighted empirical law converges to \(\mathbb{Q}^\star\), so the clean marginal converges to the tilted target \(p_0^\beta\). The catch is variance: if the base model rarely produces high-reward outputs, almost all weight lands on a few trajectories and the effective sample size collapses.

Sequential Monte Carlo fixes this by twisting the proposal with the approximate value and resampling along the way. Replace the base kernel by a proposal that follows an approximate Doob transform, \[\label{eq:twisted-proposal} \widehat{P}_k^{\leftarrow}(x_{k+1}\to \,\mathrm{d}x_k) = \frac{P_k^{0,\leftarrow}(x_{k+1}\to \,\mathrm{d}x_k)\,\widehat{h}_k(x_k)}{\widehat{g}_{k+1}(x_{k+1})}, \qquad \widehat{g}_{k+1}(x_{k+1})=\int \widehat{h}_k(y)\, P_k^{0,\leftarrow}(x_{k+1}\to \,\mathrm{d}y),\tag{80}\] where \(\widehat{h}_k\approx h_k\) is any tractable value estimate (a classifier, a reward model, a learned value, or the linearized small-tilt value of 72 ) with the exact terminal value \(\widehat{h}_0=e^{\beta r}\). If the particles are initialized from the base noisy law \(p_K^0\), they also carry the initial weight \(\widehat{h}_K(X_K)\); if one can instead sample from the twisted noisy law proportional to \(p_K^0\widehat{h}_K\), this initial weight is constant. Running \(N\) particles through \(\widehat{P}_k^{\leftarrow}\) and carrying an incremental importance weight at each reverse step, then resampling when the effective sample size drops, gives a consistent estimator of \(\mathbb{Q}^\star\).

At each reverse step the particle picks up an incremental importance weight: the ratio of the base transition weighted by the current value \(\widehat{h}_k(x_k)\) to the twisted proposal weighted by the next value \(\widehat{h}_{k+1}(x_{k+1})\). Substituting the proposal 80 , this collapses to \[\frac{ P_k^{0,\leftarrow}(x_{k+1}\to \,\mathrm{d}x_k)\,\widehat{h}_k(x_k)}{ \widehat{P}_k^{\leftarrow}(x_{k+1}\to \,\mathrm{d}x_k)\,\widehat{h}_{k+1}(x_{k+1})} = \frac{\widehat{g}_{k+1}(x_{k+1})}{\widehat{h}_{k+1}(x_{k+1})}.\] Thus the incremental weight is \[\label{eq:smc-incremental-weight} w_k(x_{k+1}) = \frac{\widehat{g}_{k+1}(x_{k+1})}{\widehat{h}_{k+1}(x_{k+1})}.\tag{81}\] It measures the one-step mismatch between the value \(\widehat{h}_{k+1}\) assumed at the current state and the value \(\widehat{g}_{k+1}\) obtained by propagating \(\widehat{h}_k\) one base step. If the twist is exact, \(\widehat{h}=h\), then \(\widehat{g}_{k+1}=h_{k+1}\) by the martingale identity and all incremental weights are one; the only nonconstant weight left is the initial twist weight if the sampler starts from \(p_K^0\). All of the approximation error is thus pushed into nonuniform weights, which resampling removes in the large-\(N\) limit; a good \(\widehat{h}\) keeps the weights near uniform and the variance low.

For conditional generation, take the terminal reward to be \(\log p(y\mid x_0)\). Then \(h_k(x_k)=\mathbb{P}^0(y\mid X_k=x_k)\), and a noisy classifier or conditional likelihood model gives a tractable \(\widehat{h}_k\). When this twist is realized as a guidance gradient \(\nabla\log\widehat{h}_k\) that only moves the particles, the proposal is exactly a classifier-guided sampler, and the weights above are the SMC correction that restores the conditional target in the large-particle limit.

8.6 Inference-time reinforcement learning↩︎

The previous subsections describe the ideal sampler for a reward tilt: if the posterior value functions \(h_k\) were known exactly, the reverse kernels would be the Doob transforms in Theorem 6. The setting here is genuinely inference-time: the pretrained reverse sampler stays fixed as the reference law, and a policy changes only how that sampler is run at test time. The RL problem is therefore over sampler choices during generation, not over retraining the base score model.

In a finite-horizon formulation, the reverse sampler takes actions during sampling. Let \(\pi=(\pi_k)\) be a possibly randomized policy: after observing \(X_{k+1}\) at step \(k\), it chooses an action \(a_k\sim\pi_k(\cdot\mid X_{k+1})\), and the reverse step uses the corresponding kernel, \[a_k\sim\pi_k(\cdot\mid X_{k+1}), \qquad X_k\sim P_k^{a_k,\leftarrow}(X_{k+1}\to\cdot), \qquad a_k\in\mathcal{A}.\] The action may be a guidance scale, a timestep choice, an injected-noise level, a rejection or resampling decision, or an additive drift correction. In a masked language diffusion sampler, it may choose how many tokens to unmask, which positions to update, or how sharply to sample from a token posterior. The terminal reward is measured only after the final object \(X_0\) is produced.

The objective is the same KL-regularized change of measure, restricted to the family of sampler modifications available at inference time: \[\label{eq:rl-kl-objective} J(\pi) = \mathbb{E}_{\pi} \left[ r(X_0) - \frac{1}{\beta}\sum_k \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left( P_k^{a_k,\leftarrow}(X_{k+1}\to\cdot) \middle\| P_k^{0,\leftarrow}(X_{k+1}\to\cdot) \right) \right].\tag{82}\] Here \(\beta>0\) is the same inverse temperature as in the variational principle 69 ; larger \(\beta\) means weaker KL regularization and stronger reward seeking. The KL term keeps the policy-induced sampler close to the pretrained sampler at each reverse step. This mirrors the variational identity: reward improvement is meaningful only relative to a reference distribution. Without this reference cost, policy optimization can exploit flaws in the reward model or collapse to a narrow set of high-scoring samples.

If the action class is rich enough to choose the whole next-step kernel, this optimization does not produce a new sampler: it recovers the Doob-transformed kernel of Theorem 6, now reached by dynamic programming rather than by tilting the path law. To see the same object appear, define the regularized value-to-go from \(x_k\) by \[v_k(x_k) = \sup_{\pi_0,\ldots,\pi_{k-1}} \mathbb{E}_{\pi} \left[ r(X_0) - \frac{1}{\beta}\sum_{j=0}^{k-1} \mathop{\mathrm{D_{\mathsf{KL}}}}\!\left( P_j^{a_j,\leftarrow}(X_{j+1}\to\cdot) \middle\| P_j^{0,\leftarrow}(X_{j+1}\to\cdot) \right) \,\middle|\, X_k=x_k \right],\] with the empty-sum convention \(v_0(x_0)=r(x_0)\). The local Gibbs variational principle gives \[P_k^{\star,\leftarrow}(x_{k+1}\to \,\mathrm{d}x_k) \propto P_k^{0,\leftarrow}(x_{k+1}\to \,\mathrm{d}x_k) \exp\left(\beta v_k(x_k)\right).\] Writing \(h_k(x_k)=\exp(\beta v_k(x_k))\) makes this exactly the normalization in 79 : the same reweighting we already have, only now the value is presented as a Bellman value-to-go rather than a posterior tilt factor. The two letters are the two coordinates of one object: the multiplicative tilt factor \(h_k\) that enters the Doob kernel and the SMC weights, and its log-domain soft value \(v_k=\beta^{-1}\log h_k\), so that \(\beta\) times its gradient (in the continuous picture) is added to the score, matching \(\nabla\log h_t^\beta=\beta\nabla v_t\). The posterior expected reward \(V_t\) of 72 is the small-\(\beta\) linearization of this soft value. Soft dynamic programming is therefore the finite-state computational form of the path-space tilt; what is genuinely new here is not the optimal kernel but how its value is obtained, since guidance methods approximate this value information while RL methods estimate it from sampled rewards.

For a parameterized policy \(\pi_\theta\) with reverse-path density \(p_\theta(x_K,\ldots,x_0)\), sampled rewards can be turned into updates by the likelihood-ratio identity \[\nabla_\theta\mathbb{E}_{(X_K,\ldots,X_0)\sim p_\theta}[r(X_0)] = \mathbb{E}_{(X_K,\ldots,X_0)\sim p_\theta} \left[ r(X_0)\nabla_\theta\log p_\theta(X_K,\ldots,X_0) \right],\] In implementations one often subtracts a baseline from the reward to reduce variance, but that is not part of the change-of-measure identity. Along a reverse diffusion chain, \[\log p_\theta(x_K,\ldots,x_0) = \log p_{\theta,K}(x_K) + \sum_k \log p_\theta(x_k\mid x_{k+1}),\] with the initial-law term disappearing from the gradient when the noisy initialization is fixed. Thus the update decomposes across the guidance decisions made along the trajectory. For the KL-regularized objective 82 , the return is replaced by the regularized return inside the brackets, with any explicit derivative of the KL penalty added when the transition family is differentiable.

This perspective is also how to read current inference-time RL methods. One can combine the base sampler with a reward-aligned proposal, spend particles on resampling and trajectory correction, or estimate a drift/Doob correction more directly. These are different approximations to the same ideal controlled path law. The useful message for these notes is not that one particular recent algorithm has replaced the others; it is that reward tilting, KL-regularized control, and Doob transforms provide the common language for comparing ways of changing the sampler at inference time.

8.7 Guidance and reward tilting in the discrete case↩︎

Everything in this section carries over almost unchanged to the discrete diffusion models of Section 7, because the target and objective never used continuous structure. The reward-tilted target 68 and the KL-regularized variational principle 69 are defined on any space; the Feynman–Kac/twisted-SMC correction of Subsection 8.5 reweights particles by the same potentials; and inference-time RL 82 differentiates the trajectory likelihood—for a discrete chain a sum \(\sum_k\log P_{k,\theta}^{\leftarrow}(x_{k+1}\to x_k)\) of log transition probabilities—with respect to \(\theta\) rather than the state, so it transfers verbatim.

The one step that changes form is guidance. The optimal tilt is again the Doob \(h\)-transform of Theorem 6, now of a jump process: with the soft value \(v_t(x)=\beta^{-1}\log\mathbb{E}[\exp(\beta r(X_0))\mid X_t=x]\), the guided reverse rates are reweighted by value differences, \[R_t^{\leftarrow,\mathrm{guided}}(y\to x) =R_t^{\leftarrow}(y\to x)\,\exp\!\big(\beta(v_t(x)-v_t(y))\big),\] the exact analogue of adding \(\beta\nabla v_t\) to the reverse drift, with the gradient replaced by a finite difference of \(v_t\) across admissible moves—the same score-to-ratio substitution as in 55 . As before \(v_t\) is intractable and is approximated, by a learned predictor of the reward, a predictor-free interpolation, or a Taylor expansion around the predicted clean token [50].

9 Itô Calculus and Girsanov Theorem↩︎

This appendix collects the stochastic-calculus facts used in the notes. The statements below are written in the smooth, non-explosive setting in which the formal calculations in the text are valid. For rigorous hypotheses and proofs, see the texts of E, Li and Vanden-Eijnden [3], Karatzas and Shreve [51] or Øksendal [52].

Itô formula, generator, and adjoint↩︎

Let \(X_t\) solve \[\,\mathrm{d}X_t=b_t(X_t)\,\mathrm{d}t+\sigma_t(X_t)\,\mathrm{d}B_t,\] where \(B_t\) is an \(m\)-dimensional Brownian motion and \(\sigma_t(x)\in\mathbb{R}^{d\times m}\). Write \(a_t(x)=\sigma_t(x)\sigma_t(x)^\top\). For a smooth test function \(\varphi\), Itô’s formula gives \[\,\mathrm{d}\varphi(X_t) = \left\langle \nabla\varphi(X_t),b_t(X_t) \right\rangle\,\mathrm{d}t +\frac{1}{2}\mathop{\mathrm{Tr}}\!\left(a_t(X_t)\nabla^2\varphi(X_t)\right)\,\mathrm{d}t +\left\langle \nabla\varphi(X_t),\sigma_t(X_t)\,\mathrm{d}B_t \right\rangle.\] The last term is a martingale increment and has mean zero under the usual integrability assumptions. Thus \[\frac{\,\mathrm{d}}{\,\mathrm{d}t}\mathbb{E}[\varphi(X_t)] = \mathbb{E}[(\mathcal{L}_t\varphi)(X_t)],\] where the infinitesimal generator is \[\mathcal{L}_t\varphi = \left\langle b_t,\nabla\varphi \right\rangle +\frac{1}{2}\mathop{\mathrm{Tr}}(a_t\nabla^2\varphi).\] If \(X_t\) has density \(q_t\), then the Fokker–Planck equation is the adjoint equation \[\partial_tq_t=\mathcal{L}_t^\ast q_t,\] where, in coordinates, \[\mathcal{L}_t^\ast q = -\nabla\cdot(b_tq) +\frac{1}{2}\sum_{i,j=1}^d \partial_i\partial_j\!\left((a_t)_{ij}q\right).\] For overdamped Langevin, \(b=-\nabla U\) and \(a=2I\), so this reduces to \[\partial_tq_t=\nabla\cdot(q_t\nabla U)+\Delta q_t.\]

Girsanov’s theorem↩︎

Girsanov’s theorem compares two diffusions that share the same initial law and the same diffusion coefficient but may carry different drifts. We fix the coefficient \(\sqrt2\,I\), the Langevin convention used in the notes. Let \(\mathbb{P}\) and \(\mathbb{Q}\) be the laws on path space \(C([0,T];\mathbb{R}^d)\) of \[\,\mathrm{d}X_t=b^\mathbb{P}_t(X_t)\,\mathrm{d}t+\sqrt2\,\,\mathrm{d}B_t \qquad\text{and}\qquad \,\mathrm{d}X_t=b^{\mathbb{Q}}_t(X_t)\,\mathrm{d}t+\sqrt2\,\,\mathrm{d}B_t.\] Under the usual absolute-continuity and integrability assumptions, the two laws are equivalent on \(\mathcal{F}_T\), the terminal \(\sigma\)-algebra of the Brownian filtration, with Radon–Nikodym density \[\label{eq:girsanov-density} \frac{\,\mathrm{d}\mathbb{P}}{\,\mathrm{d}\mathbb{Q}}\bigg|_{\mathcal{F}_T} = \exp\!\left( \frac{1}{2}\int_0^T \left\langle b^\mathbb{P}_t-b^{\mathbb{Q}}_t,\,\mathrm{d}X_t-b^{\mathbb{Q}}_t\,\,\mathrm{d}t \right\rangle - \frac{1}{4}\int_0^T \left\lVert b^\mathbb{P}_t-b^{\mathbb{Q}}_t \right\rVert^2\,\,\mathrm{d}t \right),\tag{83}\] where all integrands are evaluated at \(X_t\).

To obtain the KL divergence, take logarithms in 83 and average under \(\mathbb{P}\): \[\mathop{\mathrm{D_{\mathsf{KL}}}}(\mathbb{P}\|\mathbb{Q}) = \mathbb{E}_\mathbb{P}\!\left[\log\frac{\,\mathrm{d}\mathbb{P}}{\,\mathrm{d}\mathbb{Q}}\right] = \frac{1}{2}\mathbb{E}_\mathbb{P}\int_0^T \left\langle b^\mathbb{P}_t-b^{\mathbb{Q}}_t,\,\mathrm{d}X_t-b^{\mathbb{Q}}_t\,\,\mathrm{d}t \right\rangle - \frac{1}{4}\mathbb{E}_\mathbb{P}\int_0^T\left\lVert b^\mathbb{P}_t-b^{\mathbb{Q}}_t \right\rVert^2\,\mathrm{d}t.\] The second term is already a plain expectation of a path integral, so it is left as is. For the first term, recall that under \(\mathbb{P}\), we have \(\,\mathrm{d}X_t-b^{\mathbb{Q}}_t\,\,\mathrm{d}t=(b^\mathbb{P}_t-b^{\mathbb{Q}}_t)\,\mathrm{d}t+\sqrt2\,\,\mathrm{d}B_t\), and \[\frac{1}{2}\mathbb{E}_\mathbb{P}\int_0^T \left\langle b^\mathbb{P}_t-b^{\mathbb{Q}}_t,\,\mathrm{d}X_t-b^{\mathbb{Q}}_t\,\,\mathrm{d}t \right\rangle = \frac{1}{2}\mathbb{E}_\mathbb{P}\int_0^T\left\lVert b^\mathbb{P}_t-b^{\mathbb{Q}}_t \right\rVert^2\,\mathrm{d}t + \frac{1}{\sqrt2}\,\mathbb{E}_\mathbb{P}\int_0^T\left\langle b^\mathbb{P}_t-b^{\mathbb{Q}}_t,\,\mathrm{d}B_t \right\rangle.\] The expectation of the final integral vanishes, and we arrive at \[\label{eq:girsanov-kl} \mathop{\mathrm{D_{\mathsf{KL}}}}(\mathbb{P}\|\mathbb{Q}) = \frac{1}{4} \mathbb{E}_{\mathbb{P}}\int_0^T \left\lVert b^\mathbb{P}_t(X_t)-b^{\mathbb{Q}}_t(X_t) \right\rVert^2\,\mathrm{d}t.\tag{84}\] So far the two processes share the same initial law. If instead they start from different initial laws, the density 83 acquires the extra factor \(\,\mathrm{d}\mathop{\mathrm{Law}}_\mathbb{P}(X_0)/\,\mathrm{d}\mathop{\mathrm{Law}}_{\mathbb{Q}}(X_0)\) at \(t=0\), and the KL gains the corresponding initial term: \[\mathop{\mathrm{D_{\mathsf{KL}}}}(\mathbb{P}\|\mathbb{Q}) = \mathop{\mathrm{D_{\mathsf{KL}}}}\bigl(\mathop{\mathrm{Law}}_\mathbb{P}(X_0)\,\big\|\,\mathop{\mathrm{Law}}_{\mathbb{Q}}(X_0)\bigr) + \frac{1}{4}\mathbb{E}_\mathbb{P}\int_0^T\left\lVert b^\mathbb{P}_t(X_t)-b^{\mathbb{Q}}_t(X_t) \right\rVert^2\,\mathrm{d}t.\]

Thus, for any measurable map \(F\) on path space, \[\mathop{\mathrm{D_{\mathsf{KL}}}}(F_\#\mathbb{P}^u\|F_\#\mathbb{P}^0)\leq \mathop{\mathrm{D_{\mathsf{KL}}}}(\mathbb{P}^u\|\mathbb{P}^0).\] Taking \(F(\omega)=\omega_T\) gives the endpoint KL bound used in the one-step ULA calculation.

10 Gaussian Toolbox↩︎

This short toolbox collects facts used throughout the notes. These are good warm-up exercises for students who have not recently worked with Gaussian densities.

References↩︎