The DDPM note describes how a diffusion model turns a finite set of samples into a sampler: one learns the drift of a reversed noising process, and this drift is fitted by least-squares regression (Ho, Jain, and Abbeel 2020; Song et al. 2020). The usual Monte Carlo situation is the mirror image. We are given a target distribution \(\pi(x) = \widetilde{\pi}(x) / \mathcal{Z}_\pi\) on \(\mathbb{R}^D\) with \(\widetilde{\pi}(x) = e^{-E(x)}\) and an unknown normalizing constant \(\mathcal{Z}_\pi\). We can evaluate the energy \(E\) and its gradient, and we do not have a single sample. MCMC and SMC methods exploit the energy directly, but each new batch of samples requires a new run and nothing reusable is learned along the way. Diffusion models, on the other hand, learn a transport from noise to data that can be reused at no cost, but they need data to do so.
Can one train the diffusion model of \(\pi\) from its energy alone? A recent line of work does exactly this, under names such as adjoint matching and adjoint sampling (Domingo-Enrich et al. 2024; Havens et al. 2025; G.-H. Liu et al. 2025). These notes are my attempt to understand the simple structure underneath, and to phrase it in the language of Monte Carlo. Three facts do most of the work. First, the drift of the diffusion model of \(\pi\) is a conditional expectation of a functional of the future path, what one would call a delta in option pricing. Second, this delta has a pathwise and a likelihood-ratio representation, related by an integration by parts, and it is this integration by parts that trades samples of \(\pi\) for the score of \(\pi\). Third, the conditional expectation can be estimated by least-squares regression on paths that we simulate ourselves, so that the sampler bootstraps from its own output. There is one subtlety: on a finite time horizon, the start of the process remembers its end, and even a perfect drift leaves non-trivial importance weights. Most of the variety in the recent literature comes from the different ways of dealing with this.
A noise process that ends at \(\pi\)
The reference process
Consider the Ornstein–Uhlenbeck (OU) process of the DDPM note, started from its stationary distribution,
\[ dX_t = -\tfrac12 X_t \, dt + dW_t, \qquad X_0 \sim \pi_{\mathrm{ref}} = \mathcal{N}(0, I), \tag{1}\]
on the time interval \([0,T]\), and denote by \(\mathbb{P}\) its law on path-space. Contrary to the DDPM note, time runs here from noise at \(t=0\) to the target at \(t=T\). With \(\alpha_s = e^{-s/2}\) and \(\sigma_s^2 = 1 - \alpha_s^2\), the transitions read
\[ X_{t+s} = \alpha_s \, X_t + \sigma_s \, \mathbf{n} \tag{2}\]
for a standard Gaussian \(\mathbf{n}\sim \mathcal{N}(0,I)\) independent of \(X_t\). In other words, the transition density is \(p_s(y \mid x) = \mathcal{N}(y; \alpha_s x, \sigma_s^2 I)\): shrink, then blur. When the current time \(t\) is fixed, I will write \(\alpha \equiv \alpha_{T-t}\) and \(\sigma \equiv \sigma_{T-t}\) for the shrinkage factor and noise level over the remaining time \(T-t\).
Three properties of Equation 1 are used over and over. First, nudges fade. By Equation 2 we have \(X_T = \alpha X_t + \sigma \mathbf{n}'\) with \(\mathbf{n}'\) independent of \(X_t\); holding \(\mathbf{n}'\) fixed and moving \(X_t\) by \(\delta\) moves \(X_T\) by \(\alpha \, \delta\). Second, the process is reversible: its path distribution is invariant under time reversal, so that, conditionally on the future \(X_{t+s} = y\), the past \(X_t\) is distributed as \(\mathcal{N}(\alpha_s y, \sigma_s^2 I)\). Looking backward, one again shrinks and blurs. Third, bridges are Gaussian: conditionally on both endpoints,
\[ X_t \mid (X_0, X_T) = (x_0,x_T) \;\sim\; \mathcal{N} {\left( \alpha_t \, x_0 + \frac{\alpha_{T-t} \, \sigma_t^2}{\sigma_T^2} \, (x_T - \alpha_T \, x_0), \; \frac{\sigma_t^2 \, \sigma_{T-t}^2}{\sigma_T^2} \, I \right)} . \tag{3}\]
Sampling an intermediate point given the two endpoints costs a single Gaussian draw; no SDE needs to be simulated. Two special cases will be useful: when the process is started at \(\pi_{\mathrm{ref}}\) and only the endpoint is given, reversibility gives \(X_t \mid X_T \sim \mathcal{N}(\alpha_{T-t} X_T, \sigma_{T-t}^2 I)\); when the process is started at \(X_0 = 0\), Equation 3 gives \(X_t \mid X_T \sim \mathcal{N}(\alpha_{T-t}\sigma_t^2 X_T / \sigma_T^2, \, \sigma_t^2 \sigma_{T-t}^2 I / \sigma_T^2)\).
Some details:
The transitions follow from the variation of constants formula \(X_{t+s} = e^{-s/2} X_t + \int_t^{t+s} e^{-(t+s-r)/2} \, dW_r\); the stochastic integral is a centred Gaussian with covariance \(\int_0^s e^{-(s-r)} \, dr \, I = (1-e^{-s}) \, I\). For reversibility, the joint density of \((X_t, X_{t+s})\) at stationarity is \[ \pi_{\mathrm{ref}}(x) \, p_s(y \mid x) \; \propto \; \exp {\left\{ -\frac{\sigma_s^2 \|x\|^2 + \|y - \alpha_s x\|^2}{2 \sigma_s^2} \right\}} = \exp {\left\{ -\frac{\|x\|^2 + \|y\|^2 - 2 \alpha_s \left< x,y \right>}{2\sigma_s^2} \right\}} , \] where we have used that \(\alpha_s^2 + \sigma_s^2 = 1\); this expression is symmetric in \((x,y)\). For the bridge, conditionally on \(X_0 = x_0\) we have \(X_t = \alpha_t x_0 + \sigma_t \mathbf{n}_1\) and \(X_T = \alpha_{T-t} X_t + \sigma_{T-t} \mathbf{n}_2\). The pair \((X_t, X_T)\) is Gaussian with \(\mathbb{E}[X_t] = \alpha_t x_0\), \(\mathop{\mathrm{Var}}(X_t) = \sigma_t^2 I\), \(\mathop{\mathrm{Cov}}(X_t, X_T) = \alpha_{T-t} \sigma_t^2 I\) and \(\mathop{\mathrm{Var}}(X_T) = (\alpha_{T-t}^2 \sigma_t^2 + \sigma_{T-t}^2) I = \sigma_T^2 I\). Gaussian conditioning gives the mean in Equation 3, and the variance \(\sigma_t^2 - \alpha_{T-t}^2 \sigma_t^4 / \sigma_T^2 = \sigma_t^2 \sigma_{T-t}^2 / \sigma_T^2\) since \(\sigma_T^2 - \alpha_{T-t}^2 \sigma_t^2 = \sigma_{T-t}^2\).
Several of the papers cited below use a Brownian motion as reference instead of Equation 1. This changes nothing conceptually: the OU process is a Brownian motion in moving coordinates, since \(Y_t = e^{t/2} X_t\) satisfies \(dY_t = e^{t/2} \, dW_t\). All the formulas below translate to the Brownian case by replacing \(\alpha_s\) by \(1\) and \(\sigma_s^2\) by \(s\). I use the OU process because it comes with a natural reference distribution \(\pi_{\mathrm{ref}}\), because it is reversible, and because it is the noising process of diffusion models, which makes the identification of the next section exact.
Reweighting the end
As in these notes on Doob h-transforms, reweight each path of the reference process by a function of its endpoint,
\[ \frac{d\mathbb{Q}}{d\mathbb{P}}(X) = \frac{1}{\mathcal{Z}_\pi} \, \frac{\widetilde{\pi}(X_T)}{\pi_{\mathrm{ref}}(X_T)} = e^{g(X_T)} \qquad \text{with} \qquad g(x) = \log \frac{\pi(x)}{\pi_{\mathrm{ref}}(x)} = -E(x) + \tfrac12 \|x\|^2 + \text{const}. \tag{4}\]
The normalizing constant is exactly the unknown \(\mathcal{Z}_\pi\) since \(\mathbb{E}_{\mathbb{P}}[\widetilde{\pi}(X_T)/\pi_{\mathrm{ref}}(X_T)] = \int \widetilde{\pi}(x) \, dx\), the endpoint \(X_T\) being distributed as \(\pi_{\mathrm{ref}}\) under \(\mathbb{P}\). It shifts \(g\) by a constant and will never enter any of the drifts below. The endpoint of the reweighted process is distributed as \(\pi_{\mathrm{ref}} \times \pi / \pi_{\mathrm{ref}} = \pi\): simulating \(\mathbb{Q}\) is sampling from \(\pi\).
What does \(\mathbb{Q}\) look like? A weight that only depends on the endpoint does not modify the distribution of a path given its endpoint. Indeed, Bayes’ rule gives \(\mathbb{Q}(dX \mid X_T = y) = \mathbb{P}(dX) \, e^{g(y)} / [\mathbb{P}(X_T \in dy) \, e^{g(y)}] = \mathbb{P}(dX \mid X_T = y)\). By reversibility, \(\mathbb{P}(\cdot \mid X_T = y)\) is the OU process started at \(y\) and read backward in time. This means that \(\mathbb{Q}\) is obtained by drawing \(X_T \sim \pi\) and running the noising process backward from it: \(\mathbb{Q}\) is the diffusion model of \(\pi\), as described in the DDPM note, with noise at \(t=0\) and the target at \(t=T\). In particular, the marginal \(\mathbb{Q}_t\) at time \(t\) is the distribution \(\pi\) noised for a duration \(T-t\), i.e. the law of \(\alpha_{T-t} Y + \sigma_{T-t} \mathbf{n}\) with \(Y \sim \pi\).
The Doob h-transform describes the dynamics of \(\mathbb{Q}\) as a diffusion with an additional drift,
\[ dX_t = -\tfrac12 X_t \, dt + \textcolor{blue}{\nabla \log h(t, X_t)} \, dt + dW_t \qquad \text{with} \qquad h(t,x) = \mathbb{E}_{\mathbb{P}} {\left[ e^{g(X_T)} \mid X_t = x \right]} . \tag{5}\]
In the language of the HJB note, the optimal control is \(u^\star = \nabla \log h = \nabla V\) with value function \(V = \log h\). The rest of these notes is about learning this single vector field. As a sanity check, integrating Equation 4 over everything but \(X_t\) gives the marginal
\[ \mathbb{Q}_t(x) = \pi_{\mathrm{ref}}(x) \, h(t,x), \tag{6}\]
so that \(\nabla \log h(t,x) = \nabla \log \mathbb{Q}_t(x) + x\) and the total drift in Equation 5 equals \(\tfrac12 x + \nabla \log \mathbb{Q}_t(x)\). This is exactly the drift of the reverse diffusion in the DDPM note, with the noised data replaced by the noised target.
Some details:
Integrating out everything but \(X_t\) gives \(\mathbb{Q}_t(x) = \int \pi_{\mathrm{ref}}(x) \, p_{T-t}(y \mid x) \, e^{g(y)} \, dy = \pi_{\mathrm{ref}}(x) \, h(t,x)\). By reversibility \(\pi_{\mathrm{ref}}(x) \, p_{T-t}(y \mid x) = \pi_{\mathrm{ref}}(y) \, p_{T-t}(x \mid y)\), and since \(\pi_{\mathrm{ref}} \, e^{g} = \pi\) this also reads \(\mathbb{Q}_t(x) = \int \pi(y) \, p_{T-t}(x \mid y) \, dy\), the distribution \(\pi\) noised for a duration \(T-t\).
The look-ahead function
The function \(h(t,x)\) measures how promising the position \(x\) is at time \(t\), and the drift climbs its logarithm. By Equation 2,
\[ h(t,x) = \int \mathcal{N}(y; \alpha x, \sigma^2 I) \, e^{g(y)} \, dy : \tag{7}\]
the weight \(e^{g} = \pi / \pi_{\mathrm{ref}}\) blurred by a Gaussian kernel of width \(\sigma\) and evaluated at the shrunk point \(\alpha x\). Far from the end (\(\alpha \approx 0\), \(\sigma \approx 1\)), the drift sees a heavily blurred version of the target; close to the end (\(\alpha \approx 1\), \(\sigma \approx 0\)), it sees its fine details. This can be made precise. Writing \(m_\pi\) and \(C_\pi\) for the mean and covariance of \(\pi\), an expansion in powers of \(\alpha = e^{-(T-t)/2}\) gives
\[ \nabla \log h(t,x) = \alpha \, m_\pi + \alpha^2 \, (C_\pi - I) \, x + O(\alpha^3). \tag{8}\]
Far from the end, the drift only knows the mean of \(\pi\); one order later it knows its covariance, and finer features of \(\pi\) only appear at higher orders.
Some details:
Completing the square and using \(\sigma^2 = 1 - \alpha^2\) gives \(\mathcal{N}(y; \alpha x, \sigma^2 I) / \pi_{\mathrm{ref}}(y) = \sigma^{-D} \exp {\left\{ (2\alpha \left< x,y \right> - \alpha^2 \|x\|^2 - \alpha^2 \|y\|^2) / (2\sigma^2) \right\}} \), so that \[ h(t,x) = \sigma^{-D} \, e^{-\alpha^2 \|x\|^2 / (2\sigma^2)} \; \mathbb{E}_{\pi} {\left[ \exp {\left\{ \frac{\alpha \left< x,Y \right>}{\sigma^2} - \frac{\alpha^2 \|Y\|^2}{2\sigma^2} \right\}} \right]} . \] Differentiating gives \(\nabla \log h(t,x) = \frac{\alpha}{\sigma^2} {\left( \widehat{x} - \alpha x \right)} \) where \(\widehat{x}\) is the mean of \(\pi\) tilted by the exponential factor above; this tilted distribution is the conditional distribution of \(X_T\) given \(X_t = x\) under \(\mathbb{Q}\), and this formula will reappear in the next section. For small \(\alpha\), the tilted mean is \(m_\pi + \alpha \, C_\pi \, x + O(\alpha^2)\) and \(\sigma^2 = 1 + O(\alpha^2)\), which gives Equation 8.
The marginals \(\mathbb{Q}_t\) form a path of distributions from (nearly) \(\pi_{\mathrm{ref}}\) to \(\pi\), just as the tempering path \(\pi_\beta \propto \pi_{\mathrm{ref}}^{1-\beta} \, \pi^{\beta}\) used in annealed importance sampling and SMC (see these notes on Jarzynski and Crooks). The difference is that the diffusion path blurs the target instead of flattening it. Since \(\mathbb{Q}_t\) is the law of \(\alpha Y + \sigma \mathbf{n}\) with \(Y \sim \pi\), its mean is \(\alpha \, m_\pi\), its covariance is \(I + \alpha^2 (C_\pi - I)\), and its cumulants of order \(k \geq 3\) scale as \(\alpha^k\). The modes of \(\pi\) are merged into a single Gaussian-looking blob early on and only separate close to the end, whereas tempering keeps the modes where they are and lowers the barriers between them. A diffusion sampler decides “which mode” early, when the modes are merged, and “where in the mode” late.
A perfect sampler is a perfect proposal
For a Monte Carlo reader, \(\mathbb{Q}\) is also the zero-variance importance sampling proposal on path-space for estimating \(\mathcal{Z}_\pi = \mathbb{E}_{\mathbb{P}}[\widetilde{\pi}(X_T) / \pi_{\mathrm{ref}}(X_T)]\), a remark already made at the end of the Girsanov note. Any drift \(u(t,x)\) defines a proposal \(\mathbb{P}^u\), the law of \(dX_t = -\tfrac12 X_t \, dt + u(t,X_t) \, dt + dW_t\) started at \(X_0 \sim \pi_{\mathrm{ref}}\), and the Girsanov theorem gives computable importance weights,
\[ w(X) = \frac{\widetilde{\pi}(X_T)}{\pi_{\mathrm{ref}}(X_T)} \, \exp {\left\{ -\int_0^T u(t,X_t)^\top \, dW_t - \frac 12 \int_0^T \|u(t,X_t)\|^2 \, dt \right\}} \qquad \text{with} \qquad \mathbb{E}_{\mathbb{P}^u}[w] = \mathcal{Z}_\pi, \tag{9}\]
where \(dW_t\) denotes the noise increments used to simulate the path. Self-normalized importance sampling with these weights removes the bias of any learned sampler and estimates \(\mathcal{Z}_\pi\); resampling along the way gives an SMC sampler (Phillips et al. 2024). The weights are constant, i.e. the proposal is perfect, exactly when \(\mathbb{P}^u = \mathbb{Q}\), start included. The effective sample size of the weights Equation 9 is the natural diagnostic of a learned sampler, and “start included” is the catch discussed at length further below.
The drift is an average
Under \(\mathbb{Q}\) and conditionally on \(X_t = x\), the endpoint \(X_T\) is distributed as
\[ \mathbb{Q}(X_T \in dy \mid X_t = x) = \frac{p_{T-t}(y \mid x) \, e^{g(y)}}{h(t,x)} \, dy . \tag{10}\]
This is the Doob recipe: condition the reference on \(X_t = x\), then reweight the future by the endpoint weight \(e^{g}\), the look-ahead function \(h(t,x)\) being the normalizing constant. Equivalently, by reversibility, \(\mathbb{Q}(X_T \in dy \mid X_t = x) \propto \pi(y) \, \mathcal{N}(x; \alpha y, \sigma^2 I)\): the posterior of a clean sample \(y\) given its noisy version \(x\), exactly as in the DDPM note. The drift \(\nabla \log h\) can be written as an average under Equation 10 in (at least) two ways.
Push the endpoint
By Equation 2, the look-ahead function reads \(h(t,x) = \mathbb{E}[ e^{g(\alpha x + \sigma \mathbf{n})} ]\) for a standard Gaussian \(\mathbf{n}\). Differentiating under the expectation while holding the noise \(\mathbf{n}\) fixed gives \(\nabla h(t,x) = \alpha \, \mathbb{E}[ \nabla g(\alpha x + \sigma \mathbf{n}) \, e^{g(\alpha x + \sigma \mathbf{n})} ] = \alpha \int p_{T-t}(y \mid x) \, e^{g(y)} \, \nabla g(y) \, dy\). Dividing by \(h(t,x)\) and recognizing Equation 10,
\[ \nabla \log h(t,x) = \alpha \; \mathbb{E}_{\mathbb{Q}} {\left[ \nabla g(X_T) \mid X_t = x \right]} = \alpha \; \mathbb{E}_{\mathbb{Q}} {\left[ \textcolor{blue}{\nabla \log \pi(X_T) + X_T} \mid X_t = x \right]} , \tag{11}\]
where the second expression uses \(\nabla g(y) = \nabla \log \pi(y) + y\). Only the gradient of the energy is needed, since \(\nabla \log \pi = -\nabla E\) does not involve \(\mathcal{Z}_\pi\). The interpretation is intuitive: the drift at time \(t\) is the average force \(\nabla g\) felt at the endpoint, transported back to the present by the factor \(\alpha = \alpha_{T-t}\) by which a nudge at time \(t\) survives until time \(T\).
The reader of these notes has in fact already seen Equation 11: it is the last formula of the note on Tweedie’s formula, the generalized Tweedie formula for a reversible AR(1) model \(Y = \alpha X + Z\). The dictionary is as follows: their \(X\) is our endpoint \(X_T \sim \pi\), their \(Y\) is our \(X_t\) (by reversibility), their reference distribution \(\mathcal{N}(0, \sigma^2/(1-\alpha^2))\) is our \(\pi_{\mathrm{ref}}\) since \(\sigma^2 = 1 - \alpha^2\), and their \(p_Y / \rho\) is our \(\mathbb{Q}_t / \pi_{\mathrm{ref}} = h\) by Equation 6. The identity is a version of the target score identity of (De Bortoli et al. 2024). In the control literature, the pathwise sensitivity \(\alpha \, \nabla g(X_T)\) is an adjoint (or costate), and the methods that regress onto it are called adjoint matching and adjoint sampling (Domingo-Enrich et al. 2024; Havens et al. 2025); we will come back to the name later.
Pull toward the endpoint
One can also differentiate the Gaussian kernel instead of the path: \(\nabla h(t,x) = \int \nabla_x p_{T-t}(y \mid x) \, e^{g(y)} \, dy\) with \(\nabla_x \log p_{T-t}(y \mid x) = \alpha (y - \alpha x) / \sigma^2\). Dividing by \(h(t,x)\) gives
\[ \nabla \log h(t,x) = \mathbb{E}_{\mathbb{Q}} {\left[ \textcolor{blue}{\frac{\alpha \, (X_T - \alpha x)}{\sigma^2}} \;\Big|\; X_t = x \right]} = \frac{\alpha}{\sigma^2} {\left( \widehat{x}(t,x) - \alpha x \right)} \qquad \text{with} \qquad \widehat{x}(t,x) = \mathbb{E}_{\mathbb{Q}}[X_T \mid X_t = x]. \tag{12}\]
The quantity \(\alpha (y - \alpha x)/\sigma^2\) is the additional drift of the OU bridge pinned at \(X_T = y\), the OU analogue of the drift \(-(x - x_\star)/(T-t)\) of the Brownian bridge in the Doob note. Each endpoint that the process might reach pulls it toward itself, and the drift is the average pull. Equivalently, the process guesses its endpoint with the denoiser \(\widehat{x}(t,x)\), the time-reversed version of the denoiser \(\widehat{x}_0\) of the DDPM note, and moves toward that guess: this is Tweedie’s formula read as a drift.
Why the two agree
The Gaussian kernel \(p_{T-t}(y \mid x)\) depends on \((x,y)\) only through \(y - \alpha x\), so that \(\nabla_x p_{T-t}(y \mid x) = -\alpha \, \nabla_y p_{T-t}(y \mid x)\). An integration by parts in \(y\) then moves the derivative from the kernel onto the weight,
\[ \int \nabla_x p_{T-t}(y \mid x) \, e^{g(y)} \, dy = -\alpha \int \nabla_y p_{T-t}(y \mid x) \, e^{g(y)} \, dy = \alpha \int p_{T-t}(y \mid x) \, \nabla_y e^{g(y)} \, dy, \tag{13}\]
and Equation 12 becomes Equation 11. In Monte Carlo terms, \(h(t,x)\) is the price of the payoff \(e^{g(X_T)}\) from state \(x\) and \(\nabla \log h\) is a normalized delta: Equation 11 is the pathwise delta estimator and Equation 12 is the likelihood-ratio (score function) estimator (Glasserman 2003). For general diffusions, the link between the two is the Malliavin integration by parts used to compute Greeks (Fournié et al. 1999); with Gaussian transitions it is the one line Equation 13.
Integrating by parts “the other way”, on the marginal \(\mathbb{Q}_t(x) = \int \pi(y) \, p_{T-t}(x \mid y) \, dy\), gives a third representation. The kernel \(p_{T-t}(x \mid y)\) depends on \(x - \alpha y\), so that \(\nabla_x p_{T-t}(x \mid y) = -\alpha^{-1} \nabla_y p_{T-t}(x \mid y)\) and \(\nabla \log \mathbb{Q}_t(x) = \alpha^{-1} \, \mathbb{E}_{\mathbb{Q}}[ \nabla \log \pi(X_T) \mid X_t = x]\), which is the target score identity as quoted in the Tweedie note. With Equation 6,
\[ \nabla \log h(t,x) = \mathbb{E}_{\mathbb{Q}} {\left[ \textcolor{blue}{\frac{\nabla \log \pi(X_T)}{\alpha} + x} \;\Big|\; X_t = x \right]} . \tag{14}\]
These are three different integrands with the same conditional expectation; they differ in their variance, as discussed below.
Samples or scores
A conditional expectation is a least-squares projection, as already exploited in the DDPM note. If \(\xi_t\) is any random variable such that \(\mathbb{E}[\xi_t \mid X_t] = \nabla \log h(t, X_t)\), and \(\omega(t) > 0\) is any time weighting, then
\[ \nabla \log h = \mathop{\mathrm{argmin}}_{u} \; \int_0^T \omega(t) \; \mathbb{E}_{\mathbb{Q}} \big\| u(t, X_t) - \xi_t \big\|^2 \, dt . \tag{15}\]
Some details:
For any function \(u\), we have \(\mathbb{E}\|u(X) - \xi\|^2 = \mathbb{E}\|u(X) - \mathbb{E}[\xi \mid X]\|^2 + \mathbb{E}\|\mathbb{E}[\xi \mid X] - \xi\|^2\). Only the first term depends on \(u\), and it vanishes for \(u(X) = \mathbb{E}[\xi \mid X]\). The time weighting \(\omega(t)\) changes the variance of stochastic gradient estimates of the loss, never its minimizer.
The two representations Equation 11 and Equation 12 give two regression problems with the same solution, with integrands
\[ \xi^{\text{push}}_t = \textcolor{blue}{\alpha_{T-t} \, {\left( \nabla \log \pi(X_T) + X_T \right)} } \qquad \text{and} \qquad \xi^{\text{pull}}_t = \textcolor{blue}{\frac{\alpha_{T-t} \, (X_T - \alpha_{T-t} X_t)}{\sigma_{T-t}^2}}. \tag{16}\]
Training pairs under \(\mathbb{Q}\) are cheap to generate: draw \(X_T \sim \pi\) and set \(X_t = \alpha_{T-t} X_T + \sigma_{T-t} \mathbf{n}\), which is correct by reversibility. The pull regression needs samples of \(\pi\) and nothing else: it is denoising score matching (Vincent 2011), i.e. the way diffusion models are trained. The push regression needs no sample of \(\pi\) at the endpoints (we will see below where the endpoints come from), but it needs the gradient \(\nabla \log \pi\) there: it is the regression behind adjoint matching and adjoint sampling (Domingo-Enrich et al. 2024; Havens et al. 2025). Samples tell us where \(\pi\) is, the score tells us which way is uphill, and an integration by parts trades one for the other. A diffusion sampler of this type is a diffusion model trained with \(\nabla \log \pi\) instead of data.
In Monte Carlo language, Equation 15 is least-squares Monte Carlo, the regression step of Longstaff and Schwartz (Longstaff and Schwartz 2001): simulate paths, compute a noisy functional of each path, regress it on the current state, and read off the conditional expectation. The neural network acts as a Rao-Blackwellizer that averages many one-sample estimates of the drift into a smooth vector field. Particle methods estimate the same conditional expectations state by state, with importance sampling or SMC (Phillips et al. 2024); the regression amortizes this over all states at once.
The integrands \(\xi^{\text{push}}\), \(\xi^{\text{pull}}\) and the integrand of Equation 14 have the same conditional mean but very different fluctuations. The push integrand is of order \(\alpha\): it is small far from the end, stays bounded close to the end, and inherits the roughness of \(\nabla \log \pi\) (stiff molecular energies require gradient clipping in practice (Havens et al. 2025)). The fluctuations of the pull integrand equal \(\alpha / \sigma^2\) times the spread of \(X_T\) given \(X_t\); far from the end this spread is of order one and the variance is of order \(\alpha^2\), but close to the end the spread is about \(\sigma / \alpha\) and the variance blows up as \(1/\sigma^2\). The integrand of Equation 14 is of order \(1/\alpha\), so its variance grows like \(e^{T-t}\) far from the end. As an illustration, for a one-dimensional bimodal target with \(T=3\) and \(x = 1/2\), the conditional variances of the push, pull and Equation 14 integrands are \(0.21\), \(0.16\) and \(47.7\) at time-to-go \(2.8\); \(0.62\), \(0.53\) and \(11.4\) at time-to-go \(1.5\); and \(0.17\), \(8.3\) and \(0.57\) at time-to-go \(0.1\). This suggests pulling early and pushing late. Any mixture \(\lambda_t \, \xi^{\text{push}}_t + (1-\lambda_t) \, \xi^{\text{pull}}_t\) is also a valid integrand: the difference of two unbiased integrands is a control variate in the most literal sense, and the variance-minimizing weight is \(\lambda_t = [\mathop{\mathrm{Var}}\xi^{\text{pull}} - \mathop{\mathrm{Cov}}(\xi^{\text{push}}, \xi^{\text{pull}})] / \mathop{\mathrm{Var}}(\xi^{\text{push}} - \xi^{\text{pull}})\), with all moments conditional on \(X_t\). The family of integrands indexed by a function \(c(t)\) in (Blessing et al. 2026) is this freedom. In the sampling setting, the pull integrand cannot be used on its own, for a reason explained in the next section.
Pulling yourself up by your own samples
Both regressions average over \(\mathbb{Q}\), and therefore require endpoints \(X_T \sim \pi\), which is precisely what we do not have. Two observations get around this.
First, only the endpoints matter. Given its endpoint, a path of \(\mathbb{Q}\) is a reference path conditioned on that endpoint, so a training pair \((X_t, X_T)\) is a stored endpoint plus a single Gaussian draw, \(X_t = \alpha_{T-t} X_T + \sigma_{T-t} \mathbf{n}\) (or the bridge Equation 3 when the reference is started from a point). This is used in (Havens et al. 2025) to decouple expensive energy evaluations from cheap gradient steps.
Second, one can use the endpoints produced by the current model. Simulate the current model \(dX_t = -\tfrac12 X_t \, dt + u_\theta(t,X_t) \, dt + dW_t\) from \(X_0 \sim \pi_{\mathrm{ref}}\), store the endpoints \(X_T^{(k)}\) together with the gradients \(\nabla \log \pi(X_T^{(k)})\) in a replay buffer, perform many regression steps on pairs built from the buffer, and refresh the buffer from time to time.
The algorithm in a few lines:
repeat
simulate K paths of dX = (-X/2 + u_theta(t, X)) dt + dW from X_0 ~ N(0, I) # no gradient needed
add the pairs (X_T, grad log pi(X_T)) to the buffer # K energy gradients
repeat M times # M >> K
draw (X_T, G) from the buffer, t ~ Unif[0, T] and n ~ N(0, I)
set X_t = alpha_{T-t} X_T + sigma_{T-t} n
take a gradient step on omega(t) * | u_theta(t, X_t) - alpha_{T-t} (G + X_T) |^2What does one round of this procedure compute? Let \(\rho\) denote the distribution of the endpoints produced by the current model, and let \(\mathbb{Q}^{\rho}\) be the reference reweighted by \((\rho / \pi_{\mathrm{ref}})(X_T)\), i.e. the diffusion model of \(\rho\), with look-ahead function \(h^{\rho}\). The training pairs are exactly pairs of \(\mathbb{Q}^{\rho}\), so the regression returns \(u_{\text{new}}(t,x) = \alpha \, \mathbb{E}_{\mathbb{Q}^\rho}[\nabla \log \pi(X_T) + X_T \mid X_t = x]\). Writing \(\log(\pi / \pi_{\mathrm{ref}}) = \log(\rho / \pi_{\mathrm{ref}}) + \log(\pi / \rho)\) and applying Equation 11 to \(\rho\) gives
\[ u_{\text{new}}(t,x) = \nabla \log h^{\rho}(t,x) \; + \; \textcolor{blue}{\alpha \, \mathbb{E}_{\mathbb{Q}^{\rho}} {\left[ \nabla \log \frac{\pi}{\rho}(X_T) \;\Big|\; X_t = x \right]} }. \tag{17}\]
The update is the diffusion model of what the current model produces, plus a correction that points from \(\rho\) toward \(\pi\), averaged over where the path is heading. If \(\rho = \pi\) the correction vanishes: the target is a fixed point. Conversely, suppose that the model is already the diffusion model of its own output, \(u = \nabla \log h^{\rho}\); it is a fixed point only if \(\mathbb{E}_{\mathbb{Q}^\rho}[\nabla \log(\pi/\rho)(X_T) \mid X_t = x] = 0\) for all \((t,x)\). Letting \(t \to T\), the conditional distribution of \(X_T\) given \(X_t = x\) concentrates at \(x\), so that \(\nabla \log(\pi / \rho) = 0\) and \(\rho = \pi\). This is the same condition under which uniqueness is established in (Havens et al. 2025); I do not know of a general argument without it. In a Gaussian toy example where everything is exact (Brownian reference started at zero, Gaussian target, affine drifts), the iteration converges to the target from every initialization I tried, and a grid search finds no other fixed point.
Equation Equation 17 also explains why the score is indispensable. Running the same procedure with the pull integrand \(\xi^{\text{pull}}\) returns \(\nabla \log h^{\rho}\) exactly, by Equation 12 applied to \(\rho\): the diffusion model of the current output, which reproduces the current output. Every \(\rho\) is then a fixed point and \(\pi\) never enters the computation. With self-generated samples, the data-driven way of training a diffusion model learns nothing new; the gradient \(\nabla \log \pi\) is the only channel through which the target enters. Mixtures of the pull and push integrands remain legitimate, but they must keep some push.
In Monte Carlo terms, this is adaptive importance sampling in the spirit of the cross-entropy method (Rubinstein and Kroese 2004): sample from the current proposal, refit the proposal, repeat. Here the refit is a regression of the drift instead of a weighted maximum likelihood fit, and no importance weights are used; the replay buffer is deliberately off-policy, mixing endpoints from successive models.
The price of this bootstrap is that one only learns where one looks. In Equation 17, the target enters only through \(\nabla \log \pi\) evaluated at the model’s own endpoints, and a mode of \(\pi\) never visited by \(\rho\) contributes nothing. Far from the end, the conditional distribution of \(X_T\) given \(X_t\) under \(\mathbb{Q}^\rho\) is close to \(\rho\) itself, so that
\[ u_{\text{new}}(t,x) \approx \alpha \, {\left( \mathbb{E}_{\rho}[\nabla \log \pi(Y)] + \mathbb{E}_{\rho}[Y] \right)} . \]
Everything the early “which mode” stage of the dynamics learns about \(\pi\) is the Stein-type moment \(\mathbb{E}_\rho[\nabla \log \pi]\), which vanishes for \(\rho = \pi\) by Stein’s identity and can be small when \(\rho\) misses a mode altogether. This is the mode-seeking behaviour familiar from adaptive importance sampling and variational methods, and the usual remedies apply: broad reference distributions, long horizons (more blurring), and the importance weights Equation 9. It also helps to move slowly. Replacing the model by \((1-\lambda) \, u_{\text{old}} + \lambda \, u_{\text{new}}\), a Robbins–Monro step on the fixed-point map, is equivalent to adding a penalty \(\eta \, \mathbb{E}\|u - u_{\text{old}}\|^2\) to the regression with \(\lambda = 1/(1+\eta)\), and (Blessing et al. 2026) report that this damping helps to preserve modes.
Some details:
Adding \(\eta \, \mathbb{E}\|u(t,X_t) - u_{\text{old}}(t,X_t)\|^2\) to the regression, the minimizer satisfies the pointwise condition \((u - u_{\text{new}}) + \eta \, (u - u_{\text{old}}) = 0\), i.e. \(u = (u_{\text{new}} + \eta \, u_{\text{old}}) / (1 + \eta)\).
The computational structure is attractive. Each endpoint costs one gradient of the energy and is reused for many pairs \((t, X_t)\); no backpropagation through the SDE solver is required, and the SDE is only simulated to refresh the buffer. The number of gradient steps per energy evaluation can be in the hundreds, which is what makes these methods applicable to large molecular systems (Havens et al. 2025).
A recipe with three dials
Everything so far fits in one sentence: take a reference noise process, tie its start and its end to a prior and to \(\pi\) by reweighting, write the resulting drift as an average over where the path is heading, and learn it by regression on endpoints that one simulates oneself. The recipe has three dials. The first is the estimator of the average: push the endpoint back to the present, pull toward it, or a mixture of the two. The second is the way the start and the end are tied together: only the end reweighted, both ends reweighted, or the two ends paired independently. The third is the reference process itself: a Brownian motion, an OU process, a linear process on a space of paths, a pretrained generative model, or a jump process. Most recent diffusion samplers are settings of these three dials. The second dial is where most of the variety comes from, because of a subtlety ignored so far.
The start remembers the end
The leak
On a finite horizon, the reference process still remembers where it started: \(\mathop{\mathrm{Cov}}(X_0, X_T) = \alpha_T \, I\) with \(\alpha_T = e^{-T/2}\). Reweighting the end of the path therefore also reweights its start, as already noted in the Doob note where the conditioned process starts from \(\mu_0 \, h(0,\cdot)/\mathcal{Z}\). Here, by Equation 6 at \(t=0\),
\[ \mathbb{Q}_0(x_0) = \pi_{\mathrm{ref}}(x_0) \, h(0, x_0), \tag{18}\]
the distribution \(\pi\) noised for a duration \(T\), which is close to but not equal to \(\pi_{\mathrm{ref}}\). Good endpoints are easier to reach from some starting points than from others, and the reweighting leaks backward in time. The algorithm of the previous section starts all its paths from \(\pi_{\mathrm{ref}}\) instead of \(\mathbb{Q}_0\). Even with the exact drift, it therefore simulates the path distribution \(R(dX) = \pi_{\mathrm{ref}}(dX_0) \, \mathbb{Q}(dX \mid X_0)\), and since \(R\) and \(\mathbb{Q}\) share the same dynamics given the starting point,
\[ \frac{d\mathbb{Q}}{dR}(X) = \frac{\mathbb{Q}_0(X_0)}{\pi_{\mathrm{ref}}(X_0)} = h(0, X_0). \tag{19}\]
The leak is the part of the importance weights that no drift can remove. With the weights Equation 9, it disappears: the weighted endpoint distribution is exactly \(\pi\), whatever the horizon \(T\). Without weights, it is a bias. A Monte Carlo simulation of the exact drift in the one-dimensional bimodal example with \(T=1\) confirms Equation 19: the path weights have correlation \(0.997\) with \(h(0,X_0)\), and both have an effective sample size of \(0.725\).
How large is the bias? The endpoint distribution of \(R\) is
\[ R_T(y) = \pi(y) \; \mathbb{E} {\left[ \frac{1}{h(0, \alpha_T \, y + \sigma_T \, \mathbf{n})} \right]} = \pi(y) \, \Big( 1 - e^{-T} \, \left< m_\pi, y - m_\pi \right> + O(e^{-3T/2}) \Big). \tag{20}\]
The unweighted sampler is biased toward the mean of the reference: its mean is \(m_\pi - e^{-T} \, C_\pi \, m_\pi\) to leading order, and the residual weights \(h(0, X_0)\) have an effective sample size \(1 - e^{-T}\|m_\pi\|^2\) to leading order. In the one-dimensional example with \(T = 7\), the measured shift of the mean is \(-1.191 \times 10^{-3}\) against a prediction of \(-1.195 \times 10^{-3}\), and the measured loss of effective sample size is \(1.865 \times 10^{-4}\) against \(1.847 \times 10^{-4}\).
Some details:
The endpoint distribution is \(R_T(y) = \int \pi_{\mathrm{ref}}(x_0) \, p_T(y \mid x_0) \, e^{g(y)} / h(0,x_0) \, dx_0\). By reversibility \(\pi_{\mathrm{ref}}(x_0) \, p_T(y \mid x_0) = \pi_{\mathrm{ref}}(y) \, p_T(x_0 \mid y)\), and \(\pi_{\mathrm{ref}} \, e^{g} = \pi\), which gives the first expression in Equation 20. For the expansion, the computation leading to Equation 8 gives \(h(0,x_0) = 1 + \alpha_T \left< x_0, m_\pi \right> + \alpha_T^2 \, q(x_0) + O(\alpha_T^3)\) for a quadratic function \(q\) with \(\mathbb{E}_{\pi_{\mathrm{ref}}}[q] = 0\); the latter holds because \(\mathbb{E}_{\pi_{\mathrm{ref}}}[h(0,X_0)] = 1\) for every \(T\). Inverting and averaging over \(x_0 \sim \mathcal{N}(\alpha_T y, \sigma_T^2 I)\), the first-order term contributes \(-\alpha_T^2 \left< m_\pi, y \right>\), the square of the first-order term contributes \(\alpha_T^2 \|m_\pi\|^2\), and \(q\) only contributes at higher order. Finally \(\alpha_T^2 = e^{-T}\).
The leading term is proportional to the mean of \(\pi\) in the coordinates of the reference. Shifting coordinates so that \(m_\pi \approx 0\), using any rough estimate of the mean, removes it and the leak becomes of order \(e^{-2T}\). In the one-dimensional example, the error of the endpoint distribution at \(T = 6\) drops from \(1.6 \times 10^{-3}\) to \(1.3 \times 10^{-5}\) after centering.
Diffusion models live with the same mismatch when they start their reverse process from \(\mathcal{N}(0,I)\) instead of the exactly noised data distribution. In the control literature it is called the initial value function bias (Domingo-Enrich et al. 2024). In the HJB note, it appears in the joint optimization over the initial distribution \(q\) and the control \(u\): the optimal initial distribution is \(p_0 \, h(0,\cdot)/\mathcal{Z}\), and fixing the initial distribution to \(p_0\) is suboptimal. There are three ways out.
Make the reference forget
If the start and the end of the reference are independent, then \(h(0,\cdot) = \mathbb{E}_{\mathbb{P}}[e^{g(X_T)}]\) is constant and nothing leaks. The first option is a long horizon, \(T = \infty\) in theory and \(T\) large in practice, which is what diffusion models do. The second option is to start the reference from a point, \(X_0 = 0\), as in (Havens et al. 2025). The endpoint of the reference is then distributed as \(p_T(\cdot \mid 0) = \mathcal{N}(0, \sigma_T^2 I)\), one reweights by \(\pi / p_T(\cdot \mid 0)\), and \(\mathbb{Q}_0 = \delta_0\): the construction is exact for any horizon. The push representation Equation 11 only used the transitions Equation 2, not stationarity, and the integrand becomes
\[ \xi^{\text{push}}_t = \alpha_{T-t} {\left( \nabla \log \pi(X_T) + \frac{X_T}{\sigma_T^2} \right)} , \tag{21}\]
with training pairs obtained from the bridge Equation 3 started at \(x_0 = 0\). With a Brownian reference, this reads \(\nabla \log \pi(X_T) + X_T / T\) with \(X_t \mid X_T \sim \mathcal{N}(t X_T / T, \, t(T-t)/T \, I)\); for \(T = 1\), \(\mathbb{Q}\) is the Föllmer process, a Brownian motion conditioned to be distributed as \(\pi\) at time one.
When the reference is a pretrained generative model, one cannot start from a point, but one can choose how much noise the model injects. The noise schedule that makes the start of such a model independent of its end, the memoryless schedule of (Domingo-Enrich et al. 2024), has a simple description: it is an OU process run for an infinite time on a squeezed clock. Consider a generative model built on an interpolation \(X_t = a_t \, Y + b_t \, \mathbf{n}\) between noise (\(a_0 = 0\), \(b_0 = 1\)) and data \(Y\) (\(a_1 = 1\), \(b_1 = 0\)). After rescaling, \(Z_t = X_t / \sqrt{a_t^2 + b_t^2}\) is distributed given \(Y\) as \(\mathcal{N}(\kappa_t Y, (1 - \kappa_t^2) I)\) with \(\kappa_t = a_t / \sqrt{a_t^2 + b_t^2}\): this is the OU transition Equation 2 from \(Y\) after an OU time \(\tau_t = \log[(a_t^2 + b_t^2)/a_t^2]\). The OU time \(\tau_t\) diverges as \(t \to 0\), so the noising process is an OU process started at the data and run forever, squeezed into the interval \([0,1]\), and its final value is independent of the data. Its time reversal is a generative SDE with the same diffusion coefficient, which in the original time variable is \(\sigma_t^2 = 2 \eta_t\) with \(\eta_t = b_t {\left( \frac{\dot a_t}{a_t} b_t - \dot b_t \right)} \): the memoryless schedule. For flow matching, \(a_t = t\) and \(b_t = 1 - t\) give \(\sigma_t^2 = 2(1-t)/t\), and the blow-up at \(t = 0\) is the squeezing of an infinite OU time into a finite interval.
Some details:
Let \(t < t'\). The noising process runs from data at \(t=1\) toward noise at \(t=0\), and \(X_t \mid X_{t'} \sim \mathcal{N} {\left( \frac{a_t}{a_{t'}} X_{t'}, {\left( b_t^2 - \frac{a_t^2}{a_{t'}^2} b_{t'}^2 \right)} I \right)} \). After rescaling, \(Z_t \mid Z_{t'}\) is Gaussian with mean coefficient \(\frac{a_t}{a_{t'}} \frac{\sqrt{a_{t'}^2 + b_{t'}^2}}{\sqrt{a_t^2 + b_t^2}} = e^{-(\tau_t - \tau_{t'})/2}\) and variance \(1 - e^{-(\tau_t - \tau_{t'})}\), the OU transition over an OU time \(\tau_t - \tau_{t'}\). Writing the noising process as an SDE in the time variable \(t\) and matching the rate of growth of its conditional variance gives the squared diffusion coefficient \(-\frac{d}{dt} b_t^2 + 2 \frac{\dot a_t}{a_t} b_t^2 = 2 \eta_t\).
Tilt the start
Keep any prior \(p_0\) (a harmonic prior for molecules, or \(\pi_{\mathrm{ref}}\) itself at finite \(T\)), start the reference from it, and reweight both ends, as in the note on Schrödinger bridges:
\[ \frac{d\mathbb{Q}}{d\mathbb{P}}(X) \propto e^{g_0(X_0) + g(X_T)} \qquad \text{with} \qquad \mathbb{Q}_0 = p_0 \quad \text{and} \quad \mathbb{Q}_T = \pi. \tag{22}\]
A factor that only depends on \(X_0\) changes the distribution of the start, not the dynamics given the start. The drift is therefore still \(\nabla \log h\) with \(h(t,x) = \mathbb{E}_{\mathbb{P}}[e^{g(X_T)} \mid X_t = x]\), and everything in the previous sections applies verbatim with the push integrand \(\alpha_{T-t} \nabla g(X_T)\). The two constraints determine the two weights. At the start, \(\mathbb{Q}_0 \propto p_0 \, e^{g_0} \, h(0,\cdot) = p_0\) forces \(e^{g_0} = 1/h(0,\cdot)\): the start weight cancels the leak. At the end, \(\mathbb{Q}_T \propto e^{g} \, \psi = \pi\) with
\[ \psi(y) = \int p_0(x_0) \, e^{g_0(x_0)} \, p_T(y \mid x_0) \, dx_0, \tag{23}\]
the distribution produced at time \(T\) by the reference once its start has been tilted, so that \(e^{g} = \pi / \psi\) and the push integrand is \(\alpha_{T-t} [\nabla \log \pi(X_T) - \nabla \log \psi(X_T)]\). These are the Schrödinger equations of the SB note in disguise: \(e^{g_0}\) and \(e^{g}\) are its potentials, and \(\psi = p^{\mathrm{ref}}_T \, \widehat{\varphi}_T\). The previous subsection is the special case where \(\psi\) is the endpoint distribution of the reference. The function \(\psi\) is unknown, but its score is again an average,
\[ \nabla \log \psi(y) = \mathbb{E}_{\mathbb{Q}} {\left[ \nabla_y \log p_T(y \mid X_0) \mid X_T = y \right]} = \mathbb{E}_{\mathbb{Q}} {\left[ -\frac{y - \alpha_T X_0}{\sigma_T^2} \;\Big|\; X_T = y \right]} , \tag{24}\]
since Equation 23 is a mixture of Gaussians in \(y\) and the score of a mixture is the average score of its components, weighted by the posterior of the component (see the DDPM note); here this posterior is \(\mathbb{Q}(X_0 \in dx_0 \mid X_T = y) \propto p_0(x_0) \, e^{g_0(x_0)} \, p_T(y \mid x_0) \, dx_0\). This is Tweedie’s formula from the start to the end, and \(\nabla \log \psi\) can be learned by regressing \(-(X_T - \alpha_T X_0)/\sigma_T^2\) on \(X_T\) over endpoint pairs.
This suggests alternating two regressions. For the current \(\psi\), fit the drift with the push integrand \(\alpha_{T-t}[\nabla \log \pi(X_T) - \nabla \log \psi(X_T)]\), drawing \(X_t\) from the bridge Equation 3 given both endpoints. Since the model is started from \(p_0\), its path distribution is automatically the reference tilted by \(e^{g_0} = 1/h(0,\cdot)\) at the start and by \(e^{g}\) at the end, so that the start constraint is enforced for free; its endpoint pairs then give a new \(\psi\) through Equation 24, which updates the end potential \(g = \log(\pi / \psi)\). Alternately enforcing the two marginal constraints is the Sinkhorn algorithm of the SB note, on path-space. A natural initialization is \(\psi\) equal to the endpoint distribution of the reference, i.e. the previous subsection. This is the construction of (G.-H. Liu et al. 2025), with a Brownian reference and a harmonic prior, and of its discrete counterpart (Guo et al. 2026); with samples instead of scores, the same alternation is the diffusion Schrödinger bridge of (De Bortoli et al. 2021). Among all path distributions with the prescribed marginals at both ends, \(\mathbb{Q}\) is the closest to the reference in Kullback–Leibler divergence.
Some details:
A model with drift \(\nabla \log h_g\) started from \(p_0\) has path distribution \(p_0(x_0) \, \mathbb{P}(dX \mid X_0 = x_0) \, e^{g(X_T)} / h_g(0, x_0)\), which is \(\mathbb{P}\) tilted by \(e^{g_0(X_0) + g(X_T)}\) with \(e^{g_0} = 1/h_g(0,\cdot)\). Its endpoint pairs have joint density proportional to \(p_0(x_0) \, e^{g_0(x_0)} \, p_T(y \mid x_0) \, e^{g(y)}\), so the conditional distribution of \(X_0\) given \(X_T = y\) is the posterior appearing in Equation 24.
With an OU reference, the correction fades with the memory: \(\nabla \log \psi(y) \to -y\) as \(T\) grows. In a one-dimensional example with a non-Gaussian prior, \(\max_y |\nabla \log \psi(y) + y|\) equals \(0.54\), \(0.069\) and \(0.0088\) for \(T = 1, 3, 6\). A moderately long horizon makes the second regression almost redundant.
Pair the ends independently
A different option is to give up on reweighting the reference altogether. Draw \(X_0 \sim p_0\) and \(X_T \sim \pi\) independently and connect them with OU bridges Equation 3. The resulting path distribution \(\Pi\) is a mixture of bridges; it is not a reweighting of the reference and it is not Markovian. Nonetheless, the Markov diffusion with drift
\[ b(t,x) = -\tfrac12 x + \mathbb{E}_{\Pi} {\left[ \nabla_x \log p_{T-t}(X_T \mid x) \mid X_t = x \right]} , \tag{25}\]
the average of the drifts of the bridges passing through \(x\) at time \(t\), has the same marginals as \(\Pi\) at all times; started from \(p_0\), it ends at \(\pi\) (Brunick and Shreve 2013).
Some details:
By the tower property, \(\mathbb{E}[X_{t+\delta} - X_t \mid X_t = x] = \mathbb{E}\big[ \mathbb{E}[X_{t+\delta} - X_t \mid X_t, X_0, X_T] \mid X_t = x \big]\), and given \((X_t, X_T)\) the bridge moves with drift \(-\tfrac12 x + \nabla_x \log p_{T-t}(X_T \mid x)\), the Doob transform of the OU process with \(h = p_{T-t}(X_T \mid \cdot)\). For the marginals, each bridge has marginals solving its own Fokker–Planck equation; averaging these equations over the endpoints gives a Fokker–Planck equation for the marginals \(\Pi_t\) with drift equal to the conditional average Equation 25.
The same integration by parts as in Equation 13 gives an explicit integrand,
\[ b(t,x) + \tfrac12 x = \alpha_{T-t} \; \mathbb{E}_{\Pi} {\left[ \textcolor{blue}{\nabla \log \pi(X_T) + \frac{X_T - \alpha_T X_0}{\sigma_T^2}} \;\Big|\; X_t = x \right]} , \tag{26}\]
a special case of the general formula of the last subsection. Everything is explicit: there is no second regression and no alternation. Where the Schrödinger construction learns \(\nabla \log \psi(X_T)\), the explicit transition score appears inside the average, because the path distribution is never required to be a reweighting of the reference. The algorithm is the one of the previous section, except that current start samples and current end samples are paired at random and \(X_t\) is drawn from the bridge given both. This is the route of (Blessing et al. 2026) (with a Brownian reference, and the damping mentioned earlier); with samples instead of scores, it is bridge matching, or flow matching with independent coupling (Lipman et al. 2022; Albergo, Boffi, and Vanden-Eijnden 2023; X. Liu et al. 2022; Peluchetti 2023). The price is that \(\Pi\) is not the path distribution closest to the reference, and the learned process only matches it one time at a time. The flow matching literature suggests that couplings closer to that of the reference (Schrödinger, or optimal transport) lead to straighter and easier-to-learn dynamics; the independent pairing is the simplest and the least economical.
With a long horizon, all three agree
As \(T \to \infty\), \(\alpha_T \to 0\) and \(\sigma_T \to 1\). The point-start integrand Equation 21 tends to \(\alpha_{T-t}(\nabla \log \pi(X_T) + X_T)\), the function \(\psi\) tends to a multiple of \(\pi_{\mathrm{ref}}\) so that the Schrödinger integrand has the same limit, the independent integrand Equation 26 as well, and the leak Equation 20 vanishes. All three constructions become the diffusion model of \(\pi\). In the one-dimensional example with prior \(\pi_{\mathrm{ref}}\) and a fixed time-to-go of \(1/2\), the three drifts differ from the drift of the diffusion model by at most \(0.6\), \(0.14\), \(0.017\) and \(3 \times 10^{-4}\) for \(T = 1, 2, 4, 8\). Forgetting, tilting and independent pairing are three ways of paying for a finite horizon: the first is the simplest, the second stays closest to the reference, and the third is fully explicit.
One formula for all ties
The three constructions are instances of a single formula. For any endpoint distribution \(\Gamma(x_0, y)\) whose first marginal is the prior \(p_0\), connected with reference bridges, the Markov drift is \(-\tfrac12 x + u(t,x)\) with
\[ u(t,x) = \alpha_{T-t} \; \mathbb{E} {\left[ \nabla_y \log w(X_0, y) \big|_{y = X_T} \;\Big|\; X_t = x \right]} \qquad \text{where} \qquad w(x_0, y) = \frac{\Gamma(x_0, y)}{p_0(x_0) \, p_T(y \mid x_0)} \tag{27}\]
is the endpoint weight relative to the reference. Choose how the ends are tied, compare with how the reference ties them, transport the comparison back to the present, and average. Tilting both ends gives \(w = e^{g_0(x_0) + g(y)}\) and the integrand \(\alpha_{T-t} \nabla g(X_T)\); a point start gives \(w = \pi(y) / p_T(y \mid 0)\) and Equation 21; independent pairing gives \(w = \pi(y) / p_T(y \mid x_0)\) and Equation 26.
Some details:
The density of a bridge at time \(t\) is \(p_t(x \mid x_0) \, p_{T-t}(y \mid x) / p_T(y \mid x_0)\), so that the density of \((X_t, X_T) = (x,y)\) is \(p_{T-t}(y \mid x) \, W(x,y)\) with \(W(x,y) = \int \Gamma(x_0, y) \, p_t(x \mid x_0) / p_T(y \mid x_0) \, dx_0\). The average bridge drift equals \(\int \nabla_x p_{T-t}(y \mid x) \, W(x,y) \, dy / \Pi_t(x)\), and the integration by parts Equation 13 in the variable \(y\) turns it into \(\alpha_{T-t} \, \mathbb{E}[\nabla_y \log W(x, X_T) \mid X_t = x]\). Finally, \(\nabla_y \log W(x,y)\) is the average of \(\nabla_y \log[\Gamma(X_0, y) / p_T(y \mid X_0)]\) under the conditional distribution of \(X_0\) given \((X_t, X_T) = (x,y)\), the tower property removes this inner conditioning, and the factor \(p_0(X_0)\) does not depend on \(y\).
How nudges travel
The reference process enters the push representation Equation 11 only through the way a nudge at time \(t\) propagates to the endpoint, and it enters the algorithm through two conveniences: explicit bridges and an explicit endpoint distribution. Changing the reference changes these three things and nothing else.
Linear references, and paths as states
For a linear reference \(dX_t = -A X_t \, dt + C \, dW_t\) with a control entering as \(C \, u\), the transitions read \(X_T = e^{-A(T-t)} X_t + \xi\) with a Gaussian \(\xi\) independent of \(X_t\), and a nudge travels through the matrix \(e^{-A(T-t)}\). The computation leading to Equation 11 gives
\[ \nabla \log h(t,x) = e^{-A^\top (T-t)} \; \mathbb{E}_{\mathbb{Q}} {\left[ \nabla g(X_T) \mid X_t = x \right]} , \tag{28}\]
and the Doob control is \(u^\star = C^\top \nabla \log h\). In the eigenbasis of a symmetric matrix \(A\), each mode forgets at its own rate: a nudge along the \(k\)-th eigenvector is multiplied by \(e^{-\lambda_k (T-t)}\). This is useful when the unknown is itself a path, for example a transition path of a molecule between two metastable states. Discretize the path, \(x = (x_1, \ldots, x_K)\), with its two ends fixed, and take \(A\) a discrete Laplacian with Dirichlet boundary conditions: the reference then favours smooth paths and its high-frequency modes forget fastest. Writing \(x = \bar{x} + \sum_k c_k \, \varphi_k\), with \(\bar{x}\) any fixed path between the two states and \(\varphi_k\) the Dirichlet eigenvectors, the boundary conditions are satisfied by construction and the reference acts on the coefficients \(c_k\) only. With a point start, this is the construction of (Park, Lee, and Liu 2025), whose infinite-dimensional formulation is the \(K \to \infty\) limit of this picture; bridges are Gaussian mode by mode and everything else is unchanged.
A pretrained model as reference
For fine-tuning a pretrained diffusion or flow model \(dX_t = b(t, X_t) \, dt + \sigma_t \, dW_t\) toward a reward \(r\), reweight its paths by \(e^{r(X_T)}\), the model being run with a memoryless noise schedule so that nothing leaks. The endpoint distribution of \(\mathbb{Q}\) is the output of the model tilted by \(e^{r}\). The endpoint distribution of the reference cancels from the weight; sampling from a given target \(\pi\) with a nonlinear reference would instead require the weight \(\pi / p_T\), and \(p_T\) is unknown. This is the setting of (Domingo-Enrich et al. 2024).
Holding the noise fixed, \(X_T\) is a function of \(X_t\) with Jacobian \(J_{t \to T} = \partial X_T / \partial X_t\), which solves \(\frac{d}{ds} J_{t \to s} = \nabla b(s, X_s) \, J_{t \to s}\) with \(J_{t \to t} = I\); there is no noise term since the noise is additive. Nudges travel through this Jacobian, and the computation leading to Equation 11 gives
\[ \nabla \log h(t,x) = \mathbb{E}_{\mathbb{Q}} {\left[ a(t) \mid X_t = x \right]} \qquad \text{with} \qquad a(t) = J_{t \to T}^\top \, \nabla r(X_T), \quad \dot{a}(t) = -\nabla b(t, X_t)^\top a(t), \quad a(T) = \nabla r(X_T). \tag{29}\]
The backward equation for \(a(t)\) follows from \(J_{t \to T} = J_{s \to T} \, J_{t \to s}\), which gives \(\frac{d}{dt} J_{t \to T} = -J_{t \to T} \, \nabla b(t, X_t)\). It is the adjoint equation of the adjoint method note, and \(a(t)\) is what (Domingo-Enrich et al. 2024) call the lean adjoint; regressing \(u\) onto \(\sigma_t \, a(t)\) is their adjoint matching loss. The identity Equation 29 is exact: the derivative is taken along the reference flow, inside an average under \(\mathbb{Q}\), so derivatives of the control never appear, and there is no need to start from the full adjoint of the controlled flow and drop terms that vanish at the optimum. For the OU reference, \(\nabla b = -\tfrac12 I\) and \(a(t) = \alpha_{T-t} \nabla r(X_T)\), which is Equation 11. The price of a nonlinear reference is the absence of explicit bridges: whole trajectories must be simulated and stored, and \(a(t)\) must be computed along each of them by a backward pass of vector-Jacobian products, so that the endpoint buffer of the previous sections is not available.
Jumps
Consider now a continuous-time Markov chain on a finite set with rate matrix \(R_t(x,y)\) for jumps from \(x\) to \(y\) (the note on discrete diffusions calls it \(Q_t\); the letter \(R\) avoids a clash with \(\mathbb{Q}\)). The Doob note is written for general Markov processes, and for jump processes the reweighting multiplies the rates,
\[ \text{rates of } \mathbb{Q}: \qquad R_t(x,y) \; \frac{h(t,y)}{h(t,x)}. \tag{30}\]
Some details:
With \(f=0\), the conditioned generator of the Doob note reads \(\mathcal{L}^\star \varphi= (\partial_t + \mathcal{L})[h \varphi] / h\). Since \((\partial_t + \mathcal{L}) h = 0\), this equals \(\mathcal{L}[h\varphi]/h - \varphi\, \mathcal{L}h / h\), and with \(\mathcal{L}\varphi(x) = \sum_y R_t(x,y) [\varphi(y) - \varphi(x)]\) one obtains \(\mathcal{L}^\star \varphi(x) = \sum_y R_t(x,y) \frac{h(t,y)}{h(t,x)} [\varphi(y) - \varphi(x)]\).
Gradients become ratios, and the task is to write \(h(t,y)/h(t,x)\) as an average under \(\mathbb{Q}\) for the neighbours \(y\) of \(x\). When the state space is a group such as \(\mathbb{Z}_N^D\) and the reference noise is a translation (uniform rates to neighbouring states), the increment \(X_T - X_t\) is independent of \(X_t\), so that shifting the present by \(y - x\) shifts the endpoint by \(y - x\). Then \(h(t,y) = \mathbb{E}_{\mathbb{P}}[e^{g(X_T + y - x)} \mid X_t = x]\) and
\[ \frac{h(t,y)}{h(t,x)} = \mathbb{E}_{\mathbb{Q}} {\left[ e^{g(X_T + y - x) - g(X_T)} \mid X_t = x \right]} , \tag{31}\]
the discrete analogue of the push representation (Guo et al. 2026). When \(x\) and \(y\) differ in a single coordinate, so do \(X_T + y - x\) and \(X_T\), and only single-coordinate ratios of \(\pi\) are needed. For references without a group structure, such as masked diffusion, one can use instead
\[ \frac{h(t,y)}{h(t,x)} = \mathbb{E}_{\mathbb{P}} {\left[ e^{g(X_T)} \mid X_t = y \right]} \times \mathbb{E}_{\mathbb{Q}} {\left[ e^{-g(X_T)} \mid X_t = x \right]} , \tag{32}\]
since \(\mathbb{E}_{\mathbb{Q}}[e^{-g(X_T)} \mid X_t = x] = \sum_z p(z \mid x) \, e^{g(z)} \, e^{-g(z)} / h(t,x) = 1/h(t,x)\): a product of two independent averages, one obtained by simulating the reference from \(y\) and one along a path of \(\mathbb{Q}\) from \(x\). This is the discrete adjoint of (So et al. 2026), and the backward equation they derive for it is the backward Kolmogorov equation satisfied by \(h\). Rates must remain positive, so the squared loss is replaced by a Bregman divergence; for any Bregman divergence \(D_\varphi\), the minimizer of \(c \mapsto \mathbb{E}[D_\varphi(\xi, c)]\) is \(\mathbb{E}[\xi]\) (Banerjee, Guo, and Wang 2005), and nothing else changes. The usual choice is the generalized Kullback–Leibler divergence \(D(\xi, c) = \xi \log(\xi / c) - \xi + c\). The tilted rates inherit the sparsity of the reference, since \(R_t(x,y) = 0\) implies that the tilted rate vanishes, so that the structure of masked models (one coordinate unmasked at a time) is preserved. Tilting the start also carries over, with the same two devices applied to the ratios of \(\psi\) (Guo et al. 2026). Changing the reference therefore changes a single factor, the way a nudge travels (a number, a matrix, a Jacobian, a shift), and possibly the cost of obtaining bridges.
Where did the stochastic control go?
Many of the papers cited above are written in the language of stochastic optimal control. By the Girsanov theorem, if the controlled process \(\mathbb{P}^u\) is started like \(\mathbb{P}\), then \(\mathbb{E}_{\mathbb{P}^u}\log(d\mathbb{P}^u / d\mathbb{P}) = \tfrac12 \mathbb{E}_{\mathbb{P}^u} \int_0^T \|u(t,X_t)\|^2 \, dt\), and since \(\log(d\mathbb{Q}/d\mathbb{P}) = g(X_T) - \log \mathcal{Z}\),
\[ D_{\text{KL}} {\left( \mathbb{P}^u \,\|\, \mathbb{Q} \right)} = \mathbb{E}_{\mathbb{P}^u} {\left[ \frac 12 \int_0^T \|u(t,X_t)\|^2 \, dt - g(X_T) \right]} + \log \mathcal{Z}. \tag{33}\]
The target \(\mathbb{Q}\) therefore solves the control problem \(\max_u \mathbb{E}[g(X_T) - \tfrac12 \int_0^T \|u\|^2 \, dt]\) of the HJB note, whose solution is the Bellman control \(u^\star = \nabla V\) with \(V = \log h\). If the initial distribution is optimized as well, the optimal one is \(p_0 \, h(0,\cdot)/\mathcal{Z}\): the leak discussed above is this optimal initial distribution.
The direct approach minimizes Equation 33 over a parametrized drift by gradient descent, differentiating through simulated trajectories with the adjoint method; this is what earlier diffusion samplers do (Zhang and Chen 2021; Vargas, Grathwohl, and Doucet 2023). For a Monte Carlo reader, this is variational inference on path-space, the minimization of a reverse Kullback–Leibler divergence, with its mode-seeking behaviour and with a simulation and a backward pass at every optimization step.
The methods described in these notes never minimize Equation 33. Kullback–Leibler control is linearly solvable (Kappen 2005; Todorov 2009): the look-ahead function \(h\) solves the linear backward equation \((\partial_t + \mathcal{L}) h = 0\) of the Doob note, and \(V = \log h\) solves the nonlinear HJB equation only through the Cole–Hopf change of variables discussed in the HJB note. The optimal control is therefore the gradient of the logarithm of an expectation, which Equation 11 writes as a conditional expectation, which Equation 15 turns into a least-squares problem. This is why these methods are called matching methods. Kappen’s path integral control estimates the same conditional expectations state by state by reweighting reference paths (Kappen 2005); the regression amortizes this over all states at once.
The word adjoint has the same origin. The adjoint method note shows that, along a path with its noise held fixed, Pontryagin’s costate \(\lambda(t) = \nabla_{x(t)} g(x(T))\) solves \(\dot\lambda = -\nabla b^\top \lambda\) with \(\lambda(T) = \nabla g(x(T))\). Here \(\lambda(t) = J_{t \to T}^\top \nabla g(X_T)\), and \(\lambda(t) = \alpha_{T-t} \nabla g(X_T)\) for the OU reference. The push representation Equation 11 reads
\[ u^\star(t,x) = \mathbb{E}_{\mathbb{Q}} {\left[ \lambda(t) \mid X_t = x \right]} , \]
which is the stochastic maximum principle (Bismut 1973; Peng 1990), obtained here in one line. Bellman’s \(\nabla V\), Pontryagin’s \(\mathbb{E}[\lambda(t) \mid X_t]\) and Doob’s \(\nabla \log h\) are the same vector field, and adjoint matching regresses the control onto the costate.
Mental models to keep
- To sample from \(\pi\) with a diffusion, reweight a noise process by where it ends. With an OU reference started at \(\pi_{\mathrm{ref}}\), the result is the diffusion model of \(\pi\), which blurs \(\pi\) where annealing would temper it.
- The drift is a delta: a conditional average over where the path is heading, computed pathwise (push the endpoint back) or by likelihood ratio (pull toward the endpoint).
- Samples tell us where \(\pi\) is, the score tells us which way is uphill, and an integration by parts trades one for the other.
- Averages are regressions (least-squares Monte Carlo): train on your own endpoints and redraw the middle from bridges.
- Training on your own samples is adaptive importance sampling without weights: one only learns where one looks, so move slowly.
- A learned sampler is a proposal, and path weights correct it.
- On a finite horizon the start remembers the end, and even the exact drift leaves weights \(h(0, X_0)\): make the reference forget, tilt the start, or pair the ends independently.
- Changing the reference changes how nudges travel: a number, a matrix, a Jacobian, a shift.
- Bellman, Pontryagin and Doob describe the same vector field, and no control problem is ever solved by optimization.
A few things I would try, none of which I have tested. Center (and roughly whiten) the reference at a crude estimate of the mean of \(\pi\), since the leak is driven by \(m_\pi\). Mix the integrands with a time-dependent weight, pulling early and pushing late while always keeping some push. Use the path weights Equation 9 not only to evaluate a learned sampler but to correct it with SMC, which also removes the leak. When a non-trivial prior matters, prefer a moderately long OU horizon to a point start, since the Schrödinger correction then nearly vanishes.