Diffusion-Based Posterior Sampling: A Feynman-Kac Analysis of Bias and Stability
Matias G. Delgadino∗
Sebastien Motsch⋆
arXiv:2605.06538v1 [cs.LG] 7 May 2026
William Porteous∗
Advait Parulekar†
Sanjay Shakkottai†
Abstract Diffusion-based posterior samplers use pretrained diffusion priors to sample from measurement- or reward-conditioned posteriors, and are widely used for inverse problems. Yet their theoretical behavior remains poorly understood: even with exact prior scores, their outputs are biased, and in low-temperature regimes their discretizations can become unstable. We characterize this bias by introducing a tractable surrogate path connecting the true posterior to a standard Gaussian and comparing it to the sampler’s path. Their density ratio satisfies a parabolic PDE whose reaction term measures the accumulated bias. A Feynman-Kac representation then expresses the RadonNikodym correction as an explicit path expectation, identifying which posterior regions are over- or under-sampled. We apply this framework to DPS and STSL, a related sampler. For DPS, the correction is an Ornstein-Uhlenbeck path expectation coupling the data conditional covariance with the reward curvature, revealing where DPS overor under-samples. Next, we reinterpret STSL as an auxiliary drift that steers trajectories toward low-uncertainty regions, flattening the spatially varying part of the DPS reaction term. Finally, we characterize early guidancestopping, a common mitigation for low-temperature instabilities caused by forward-Euler integration of the vector field. Together, these results clarify sampler bias, explain existing correctives, and guide stable variant designs.
1
Introduction
Diffusion and score-based generative models [Sohl-Dickstein et al., 2015, Ho et al., 2020, Song and Ermon, 2020, Song et al., 2021] have become the workhorse of modern generative modeling, powering text-to-image systems [Rombach et al., 2022, Ramesh et al., 2021, Saharia et al., 2022, Dhariwal and Nichol, 2021] and an expanding range of scientific and medical inverse problems [Song et al., 2022]. Their flexibility hinges on a single learned object (the score ∇ log ρt of the noised data distribution) which can be repurposed across downstream tasks without retraining. A canonical such task is sampling from a posterior of the form µy (x) ∝ eRy (x) ρ∗ (x), encompassing both classical inverse problems y = A(x)+ϵ and reward-tilted generation [Daras et al., 2024]. Even granting access to a perfect score oracle, posterior sampling is computationally ∗ Department of Mathematics, UT Austin. † Department of Electrical and Computer Engineering, UT Austin. ⋆ School of Mathematical and Statistical Sciences, Arizona State University.
Emails: {matias.delgadino, wporteous}@utexas.edu; {advaitp, sanjay.shakkottai}@utexas.edu; [email protected]. Preprint.
intractable in the worst case [Gupta et al., 2024], so practical algorithms rely on heuristic guidance that approximates the time-dependent posterior score. An early and influential such heuristic is Diffusion Posterior Sampling (DPS) [Chung et al., 2023]: it replaces the intractable conditional score ∇xt log p(y | xt ) by the gradient of the reward evaluated at the Tweedie posterior mean x̂0 (xt ) = E[X0 | Xt = xt ] [Robbins, 1956], yielding a plug-and-play guidance compatible with any pretrained score network. Its simplicity has made DPS the de-facto baseline for inverse problems and inspired a line of research targeting its known weaknesses: manifold-constrained gradients [Chung et al., 2022], denoising restoration [Kawar et al., 2022], pseudoinverse-guided diffusion [Song et al., 2023], latent-space extensions [Rout et al., 2023b], second-order Tweedie corrections [Rout et al., 2023a, Boys et al., 2024], proximal approaches to decrease the gradient computation burden [Rout et al., 2025], filtering and SMC-based reweightings [Dou and Song, 2024, Wu et al., 2024, MOUFAD et al., 2025], and recent drift-control schemes [Ren et al., 2026, Guo et al., 2026, Anil et al., 2026]. Yet despite this flurry of activity, two basic questions remain open. First, the DPS approximation is biased even for Gaussian-mixture priors with quadratic rewards; but which samples does this bias over- or under-represent, and why do correctives like STSL improve performance? Existing analyses establish convergence under restrictive assumptions on the prior or measurement operator [Xu and Chi, 2024, Parulekar et al., 2025, Moitra et al., 2026] or treat the algorithm as a black box, leaving its preferred classes unexplained. Second, in the low-temperature regime needed for hard measurement constraints in image inverse problems, standard DPS is numerically unstable. Practitioners routinely fall back on early guidance-stopping and trajectory-dependent step sizes, but the effect of these heuristics on the sampled distribution has never been quantified. Contributions. We close both gaps with a unified analysis based on the classical FeynmanKac formula [Karatzas and Shreve, 1991], complementing recent stochastic-analytic perspectives on guidance [Bruna and Han, 2024, Ren et al., 2026, Guo et al., 2026]. (i) An exact bias formula for DPS. In Section 3, we derive a pointwise Radon–Nikodym weight ω(x) relating the DPS-induced distribution to the true posterior. Using trajectory reversal, this weight can be written as an expectation over Ornstein–Uhlenbeck paths. The spatially varying part of the reaction term cDPS captures the alignment between conditional covariance and reward curvature, identifying where DPS over- or under-samples. (ii) STSL-type bias reduction. We identify the spectral structure of the DPS bias: it is amplified where the data manifold has high conditional uncertainty along rewardsensitive directions. This motivates an auxiliary potential drift ∇U that steers trajectories toward lower-uncertainty regions and flattens the spatially varying part of the DPS reaction term. The trace-of-covariance choice U (t, x) = tr(Σt (x)) recovers the empirically successful STSL correction [Rout et al., 2025] and connects naturally to recent neural drift-control approaches [Ren et al., 2026, Guo et al., 2026]. (iii) Quantifying low-temperature instability and early stopping. Finally, in Section 5 we show that the standard implementation of DPS systematically violates the stability condition of the forward-Euler of the bias vector field, leading to oscillations. Practitioners have implemented early-guidance-stopping2 as a way to mitigate them. We are the first to characterize the early-guidance-stopping heuristic as a weighted version of the prior.
2
Background and Related Work
Score Based Generative Models. We consider the problem of sampling from a distribution whose density is given by ρ∗ (x). Score based generative models use a trained score network sθ∗ (x, t) ≈ σt2 ∇ log pt (x) + x to approximate the denoising process dXt = (Xt + 2∇ log ρt (Xt )) dt + 2With early-guidance-stopping,
√
2 dBt
(1)
the diffusion evolves with a drift consisting of the (score + reward guidance); after a pre-specified time, the guidance term is removed and the drift is purely due to the score.
2
The key to implementing this is that the score network sθ∗ (x, t) can be trained from samples 2 p . X ∼ p as: θ∗ = argminθ EX,η X − sθ (e−t X + 1 − e−2t η, t) Throughout, we use a subscript t to denote a noised distribution, so pt := edt p(et x) ∗ N (0, 1 − e−2t ) is the marginal of the standard Ornstein-Uhlenbeck (OU) noising process √ dXt = −Xt dt + 2 dBt , X0 ∼ ρ∗ , (2) which interpolates between ρ∗ at t = 0 and the standard Gaussian γ = N (0, I) as t → ∞. Equation (1) is the Anderson reversal of Equation (2) [Anderson, 1982a], and sampling via this reverse process runs in polynomial time [Rombach et al., 2022, Song and Ermon, 2020, Ramesh et al., 2021, Saharia et al., 2022], provided the score network has been trained in advance on samples from the target distribution. This compares favorably to classical approaches such as Langevin dynamics [Vempala and Wibisono, 2022], whose convergence rate is instance-dependent and can be arbitrarily slow, see Section B.1. Posterior Sampling. A natural application of score-based models is to inverse problems and posterior sampling. The score network characterizes a prior ρ∗ , and at test time one tilts the samples by a log-likelihood Ry (x) to target the posterior µy := eRy ρ∗ /Z. The main challenge is that the posterior score ∇ log(µy )t cannot be easily computed from ∇ log pt .3 A range of approximate algorithms have been proposed to circumvent this. A central theme is the use of the prior scores ∇ log pt through Tweedie’s formula to obtain realistic-looking samples even when posterior sampling. Specifically, E[X0 | Xt = xt ] is used4 as a computationally tractable proxy for the initial condition X0 (i.e., the value that would result if the reverse diffusion were run to completion starting from Xt = xt ) and is fed into the reward model Ry (·) when modifying the drift at test time, see Section B.2. Although the resulting samples are not formally drawn from µy , this heuristic performs well in practice. Feynman-Kac formulas. The sampling literature often focuses on error bounds for approximate sampling algorithms [Lee et al., 2023, Chen et al., 2023, Vempala and Wibisono, 2022], see Section B.1. These are instantiated as upper bounds on the KL/TV/χ2 distance between the distribution of the sampling algorithm and the ground truth, and under favorable circumstances can be shown to be polynomially or exponentially small in the parameters of the instance. In posterior sampling, such an error is known to be large As such, a KL bound is often vacuous, unless it is accompanied by strong assumptions about the instance. Rather than focusing on bounding this error, we apply machinery that allows us to explicitly track the Radon-Nikodym derivative of approximate posterior sampling algorithm with respect to the true posterior. In particular, we will exploit the Feynman-Kac representation. Consider two time-dependent densities evolving under possibly different transport and reaction fields: ∂t πt = −∇ · (vt πt ) + ∆πt + ft πt (3) ∂t πt′ = −∇ · (vt′ πt′ ) + ∆πt′ + ft′ πt′ |{z} |{z} | {z } transport
diffusion
reaction
PDEs containing only the transport and diffusion terms are Fokker-Planck equations, and their solutions can be represented as the marginal densities of an SDE with corresponding drift and diffusion. The reaction term ft πt introduces a path-dependent weighting. In the special case ft = −κt with κt ≥ 0, this corresponds to killing, or early termination, of the SDE at rate κt . When ft is positive, the reaction term may instead be interpreted as spawning, birth, or branching at rate ft . Thus the resulting solution is generally not a probability density: it is an unnormalized measure, whose total mass evolves according to the cumulative effect of killing and spawning. After normalization, it gives the density of the corresponding weighted process at the terminal time. We can write the PDE for the ratio of the marginals gt := πt′ /πt as ∂t gt = ∆gt + bt · ∇gt − ct gt , 3 The noising process and the tilt do not commute: (eRy ρ
(4)
R ∗ )t ̸= e y ρt . 4 From Tweedie’s formula, E[X | X = x ] = c t t 0 0,t + c1,t ∇ log pt is affine in the prior score, and thus can
be efficiently estimated at test-time.
3
for an appropriate choice of bt , ct . Letting (Zs )s∈[0,t] be the diffusion process associated to the stochastic characteristics √ dZs = bs (Zs ) ds + 2 dWs , (5) the Feynman-Kac representation of (4) reads Z t ct−s (Zs ) ds Z0 = x , gt (x) = E g0 (Zt ) exp −
(6)
0
Please see Appendix A for some elaboration of these techniques.
3
Surrogate path and the Bias of DPS
This section develops a general surrogate-path framework for analyzing diffusion-based posterior samplers. Given a reward Ry : Rd → R our goal is to sample from the posterior Ry that arises as an exponential tilt of the prior: µy = e Zρ∗ , where Z ∈ R is a normalization − constant. Our starting point is to create a surrogate path → µ t : [0, ∞) → P(Rd ), that interpolates µy with the standard Gaussian. This path is designed such that we can track the Radon-Nikodym derivative between the marginals of this path, and the marginals of the sampler using the Feynman-Kac machinery. As we will see, the algorithm can often be fruitfully instantiated as the SDE that results from dropping the reaction term from the PDE describing the evolution of the surrogate path. − A natural (and almost exhaustive) family of paths is given by t 7→ → µ := ht ρt , where the t
Zt
function h· : [0, ∞) × Rd → R only needs to satisfy h0 = eRy and h∞ ≡ C to match the end −: points of the interpolation. We first describe the evolution for the surrogate path t 7→ ← µ t − = Lemma 1 (Informal, see Lemma 3). For any time horizon T , the reverse trajectory ← µ t → − µ T −t satisfies ← − = −∇ · (x← − ) − 2∇ · (∇ log → − ← − ) + ∆← − − c[h − − ∂ µ µ µ µ µ ,→ ρ ]← µ t
t
T −t
t
t
t
T −t
T −t
t
← − =→ − µ µT 0 (Surrogate Path) − ρ t and the specific where c[ht , ρt ] is an appropriate scalar field that depends on noised prior → choice of ht . Because we have access to a score network sθ (x, t) = ∇ log ρt (x), fixing a specific surrogate − − trajectory → µt = µ t , or equivalently a function ht , we have direct access to the score ∇ log → → − ∇ log ht + ∇ log ρ t . The algorithm path we consider does not contain a reaction term and solves directly ← − ∂t ν−t = −∇ · (x← ν−t ) − 2∇ · (∇ log → µ T −t ← ν−t ) + ∆νt (Algorithm Path) ← ν−0 = N (0, I). The solution ← ν− = Law(Y ) is obtained as the law of the associated SDE t
t
√ − dYt = Yt + 2∇ log → µ T −t (Yt ) + 2dBt Y0 ∼ N (0, I).
(Algorithm SDE)
The difference between (Surrogate Path) and (Algorithm Path) is merely the presence/absence of the reaction term. Using the Feymann-Kac formula (6) we can express their ratio as weighted expectation over the paths (Algorithm SDE): " ! # Z T µT (x) − ≈ EY ∼(Algorithm SDE) exp − c[hT −t , → ρ T −t ](Yt ) dt YT = x , νT (x) 0 and characterizes how the reaction term creates a mismatch between the output of the algorithm Law(YT ) and the true posterior µy . Note that due to this modification in (Algorithm Path), some amount of bias is unavoidable. Indeed in the worst case, any path beginning at the posterior µy and ending in a tractable distribution like N (0, I), is generated 4
Algorithmic path bias
Surrogate path
tilt pretrained
T
time
Figure 1: The blue dotted line illustrates the path taken by the standard forward OU process − −OU , whose ρ t . The violet line illustrates the OU process ← µ from ρ∗ , ρ ⃗t and its reversal ← t ← − reversal µ t we cannot track at inference time. The red line illustrates the surrogate path R (x̂ ) → − S µ DP = e y t ρ∗/Z we construct, with the same beginning and end points as µ ⃗ t . The orange t line denotes the algorithm path νtDPS which disregards the reaction term results in a sample from νyDPS with an unavoidable bias. by a evolution that we cannot compute in polynomial time, see Remark 1 and [Gupta et al., 2024]. Nevertheless, some paths inspire useful approximations, that yield good empirical results; see for instance [Bruna and Han, 2024, Parulekar et al., 2025, Ren et al., 2026]. OU Interpolation. A canonical (Surrogate Path) is given by the solution to the OU OU − dynamics: → µ OU = ht ρt/Zt from µy to N (0, I). This is in fact, the only case where the t − reaction coefficient c[ht , → ρ t ] = 0. In this case, we can use the Feymann-Kac formula to − → − OU express the quotient ht = → µ OU t / ρ t through the representation h i Ry (X0 ) hOU (x) = E e X = x . (7) t t − −OU = → → − µ OU To solve for (Algorithm Path) ← µ t T −t , we need access to the score ∇ log µ t = → − OU ∇ log ht + ∇ log ρ t . Solving (7) at test-time is not tractable, and therefore we cannot efficiently get an acceptable approximation to ∇ log hOU t . − This discussion motivates the design problem: create a surrogate path t 7→ → µ t with both a tractable score and a small reaction term. Diffusion Posterior Sampling. We show below that the DPS algorithm [Chung et al., 2023] can be interpreted as stemming from the following (Surrogate Path). 1 Ry (x̂t (x)) → − S µ DP (x) = e ρt (x), t Zt − − which retains the correct endpoints → µ =µ ,→ µ 0
y
∞ = γ.
(DPS Surrogate path) − S Crucially, the score ∇ log → µ DP t
can be written in terms of x̂s (x) := E[X0 | Xs = x],
Σs (x) := Cov(X0 | Xs = x),
both of which are computable from the diffusion score oracle sθ and its Jacobian ∇sθ via Tweedie’s formula [Robbins, 1956], see Appendix D.1. Heuristically, the (DPS Surrogate path) is obtained from the OU interpolation by swapping the conditional expectation inside the exponential (7) for S R(X0 ) hDP (x) = eR(E[X0 |Xt =x]) ̸= hOU |Xt = x]. t t (x) = E[e
We can identify the true reversal of the (DPS Surrogate path), −DP S = ∆← −DP S − 2 ∇· ∇ log → − −DP S − ∇· (x ← −DP S ) − c ← −DP S , S← ∂t ← µ µ µ DP µ DP S (T − t, x) µ t t t T −t µ t t (DPS Surrogate PDE) 5
with an explicit reaction coefficient 1 2 2 d tr Σt (x)(D Ry )(x̂t (x))Σt (x) + Σt (x)∇Ry (x̂t (x)) − dt cDP S (t, x) = − log Zt , (et − e−t )2 (8) The DPS algorithmic path. The difficulty with implementing Equation (DPS Surrogate PDE) as an SDE is the reaction term. We can construct an alternate PDE that has the same transport and diffusion term but no reaction term − ∂← ν− = ∆← ν− − 2 ∇· ∇ log → µ DP S ← ν− − ∇· (x ← ν− ), (DPS path) t
t
t
T −t
t
t
− S Using the identities of Appendix D.1, for the score ∇ log → µ DP , we get the (Algorithm SDE) t approximated by the DPS algorithm as ( √ 2 dYt = Yt + 2∇ log ρT −t (Yt ) + et −e 2 dBt , −t ΣT −t (Yt ) ∇Ry (x̂T −t (Yt )) dt + Y0 ∼ γ, ← − ν DP S
(DPS SDE)
T Applying the Feynman-Kac formula (6) to the quotient ← −DP S , we obtain the following µ T characterization of the bias of the DPS algorithm. S Theorem 1. The terminal law νyDP S := ← ν−DP of the DPS-SDE (DPS SDE) differs from T the true posterior µy by a pointwise multiplicative weight:
µy (x) = ω(x) νyDP S (x).
(9)
The weight ω admits two equivalent Feynman–Kac representations in terms of the reaction term cDP S defined in (8): (i) Backward path (condition on the DPS denoising process arriving at YT = x): # " Z T → − µ T (Y0 ) YT = x . (10) ω(x) = EY ∼(DPS SDE) exp − cDP S (T − s, Ys ) ds γ(Y0 ) 0 (ii) Forward path (condition on the OU process (2) starting at X0 = x): " # Z T 1 γ(XT ) = EX∼OU → exp cDP S (s, Xs ) ds X0 = x . − ω(x) µ T (XT ) 0
(11)
Both path functionals are expressible in terms of quantities obtainable from the score oracle and its Jacobian via Tweedie’s formula. Importance-weighting DPS samples by ω recovers µy exactly. Discussion of Theorem 1. Equations (10) and (11) give us an explicit handle on the 1 distribution of the DPS sampler. Writing Equation (9) as ω(x) µy (x) = νyDP S (x) shows that, relative to ground truth, DPS under-samples points x where ω(x) > 1 and oversamples where ω(x) < 1. We illustrate this with a simple mixture-of-gaussians prior in Fig. 2. We can also simplify the expression for ω to get an approximate expression with a geometric − interpretation. Using the fact that Zt does not depend on x, and considering that → µT ≈ γ for large T , we get the following approximation: " # Z T 1 Z0 ≈ EX∼OU exp c̃DP S (s, Ys ) ds X0 = x ω(x) ZT 0 tr Σs (x)(D 2 Ry )(x̂s (x))Σs (x) 2 where c̃DP S (s, x) = + Σs (x)∇Ry (x̂s (x)) . c̃DP S formalizes an (es −e−s )2 interplay between the prior and the reward model. Concretely, diagonalizing the conditional covariance, Σt (x) =
d X
λi (t, x) ui (t, x) ui (t, x)⊤ ,
λi (t, x) ≥ 0, {ui (t, x)}di=1 ⊂ Rd orthonormal,
i=1
6
Figure 2: True Posterior versus DPS Samples: Dashed line is measurement constraint Ax = y, A(x1 , x2 ) = (0, x2 ), y = (0, −2.5). (a) Prior (analytic): ρ0 , 4-component, equal-weight, 2 ρ0 where R(x1 , x2 ) = −2 ∥Ax − y∥ ; Gaussian mixture; (b) Posterior (analytic): µy = exp(R) Z 1 (c) Weight (log-scale): ω(x) , 20 trajectory estimate, darkest is undersampling, lightest is oversampling, gray background not computed; (d) DPS Samples: 5 × 105 samples show x1 extremal modes are nearly absent while x2 < y2 is over-sampled and x2 > y2 undersampled. we can rewrite the reaction coefficient (8) as c̃DP S (t, x) =
d X 1 i λ2 (t, x) γR (t, x), (et − e−t )2 i=1 i
(12)
where the coefficients i γR (t, x) := ui (t, x)⊤ (D2 R)(x̂t (x)) ui (t, x) + ui (t, x) · ∇R(x̂t (x))
2
quantify how sharply the reward R varies along the eigendirection ui . λi (t, x) is large along directions of high posterior uncertainty about X0 given Xt = x (the local tangent directions i of the data manifold at x̂t (x)), while γR measures the reward sensitivity along those same directions. The term c̃DP S is hence amplified precisely where the data manifold is broad and the reward landscape is active along the same axes.
4
Bias Reduction
We see in Theorem 1 that the ratio between the density of thehDPS R sampler and the true i
posterior can be expressed (approximately) as: EX∼DPS SDE e− cDPS (s,Xs ) ds |X0 = x . This gives a clear design goal: it is beneficial to design paths that result in small variations in cDPS over the trajectories. This would correspond to a smaller reaction term in Equation (Surrogate Path), and a smaller bias when we implement the corresponding algorithm. For instance, we can add an extra potential vector field ∇U to the SDE (Algorithm SDE), that drives trajectories to regions where cDP S has small oscillations. In terms of the algorithm (Algorithm Path), this amounts to solving − − ← − S← ∂t ← ν−t = ∆← ν−t − ∇ · (x ← ν−t ) − 2∇ · (∇ log → µ DP T −t ν t ) + r ∇ · (∇U ν t ), {z } | {z } | Equation (DPS path)
extra guidance
where the drift intensity r ≥ 0 is a hyperparameter. As we see below, such a change in drift can readily be matched with a corresponding change in reaction term that reinterprets the Surrogate Path with the updated drift and a modified reaction term. Interpreting the drift as a reaction term. Just as diffusion can be recast as a drift − − − involving the score, ∆→ ρ t = ∇ · (∇ log → ρt→ ρ t ), the additional drift ∇U can be recast as a reaction term through the tautological identity: ∇ · (∇U ← ν−t ) cU = = ∆U + ∇U · ∇ log ← ν−t . ← − νt 7
In other words, we can rewrite the DPS Surrogate PDE as: −DP S = ∆← −DP S −2 ∇· (∇ log → − −DP S −(c ∂← µ µ µ DP S +x+r∇U )← µ t
t
t
T −t
t
← −DP S DP S (T −t, x)+r cU ) µ t
In the language of Theorem 1, this modifies the reaction term to ceff = cDP S + r cU . As a consequence, excessively large r is counterproductive: the reaction term becomes dominated by r cU and the original bias structure is lost. Remark 1. There exists in principle a potential U ∗ satisfying cDP S + cU ∗ = 0, which would eliminate the bias exactly. Computing U ∗ directly is exponentially slow; recent work instead approximates it via a variational characterizations, training a non-linear [Guo et al., 2026] or linear [Ren et al., 2026] neural network for each specific reward. STSL as a special case. STSL [Rout et al., 2025] chooses a potential U that drives the trajectory toward low-uncertainty regions of the initial condition X0 . Up to constants, the choice is d X U (t, x) = tr(Σt (x)) = λi (t, x) ≥ 0. i=1
Since the λi are non-negative, smaller values of U correspond to smaller spread in the dominant eigendirections of Σt (x), which in turn flattens the spatially varying part of cDP S in (12). In practice this yields a better algorithm with reduced output uncertainty [Rout et al., 2023a].
5
Numerical Instabilities of the DPS Algorithm
As discussed in earlier sections, the algorithm evolution introduces a bias when compared to the surrogate evolution. In this section, we study a different issue with the actual implementation of the DPS algorithm, namely instability of the evolution close to the constraint manifold. In the context of the DPS algorithm proposed in Chung et al. [2024], we show that the instability unavoidably occur in the space parallel to the data manifold due to the systematic violation of the forward Euler stability condition. This phenomenon has indeed been observed in practice. To mitigate this, a common practice is to “turn off” reward guidance close to the data manifold, which we refer to as early guidance stopping. In other words, the stochastic evolution starts with both the score (corresponding to the untilted prior) and reward guidance drift terms until some intermediate time tstop ∈ (0, T ), after which the diffusion proceeds with only the untilted score. We show that early guidance stopping can be explicitly characterized as an appropriately weighted tilt of the prior. Instability of DPS. We first examine the algorithmic implementation of the DPS algorithm, which is conceptualized as an approximation to the solution of the SDE (DPS SDE). The exact algorithm proposed by the authors of [Chung et al., 2024] is given in Appendix C. The DPS algorithm progressively denoises over discrete time-steps, with the reward guidance weighted at each time-step through a guidance schedule {ζi }N i=1 that is taken to be trajectorydependent, α ζi = , (13) ∥y − A(x)∥2 where A : Rd → RL is a general observation operator, y ∈ RL is the observation, and α ∈ [0.2, 1] is a hyperparameter chosen depending on the inverse problem to be solved. In practice, the choice of bias schedule significantly affects the performance of the algorithm. A first observation is that (13) does not account for the time discretization ∆ti of the SDE (DPS SDE); effectively, this corresponds to multiplying the biasing vector field by a time-dependent factor. In terms of the surrogate (Surrogate Path), this corresponds to the curve e−αηt ∥A(x̂t (x))−y∥2 ρt − t 7→ → µt = , (14) Zt with annealing schedule (for the linear noising schedule {βi }1000 i=1 used in the classical DDPM, see Appendix F for details) given by ηt ≈
105 √ . 1 + 300 t 8
(15)
Oscillations of ( f(xt) − 1k) ⋅ 1
observed on MNIST close for t → 0
Last 4 iterates, recentered at their mean
Figure 3: (Top Left) A pictorial depiction of instability - as the trajectory approaches the data manifold, the large effective guidance schedule triggers oscillations in the trajectory. (Top Right) An exhibition of these oscillation on a posterior sampling task with an MNIST prior. (Bottom) A plot of the last four iterates of DPS, re-centered about their mean. The guidance tilted the distribution towards the digit 3. We observe periodic oscillations in pixel space (the deviations from the mean at alternate time steps are similar to each other). Please see Figure 4 and Appendix H for details. The key takeaway is that that schedule weight is large for a reasonable choice of hyperparameters, with the qualitatative implication of strong enforcement of measurement constraints as we approach the data manifold (small t / low-temperature regime). Indeed notice that this path-dependent bias schedule yields the target density of: 5
µTarget = y
e−α 10 ∥A(x̂t (x))−y∥2 ρt , Z
(16)
in which the reward is unsquared and the constraint {A(x) = y} has a large weight. Inevitable Oscillations. The unsquared residual in target (16) distorts the dynamics in a way that no choice of step size can repair. To see this, consider the one-dimensional example where gradient flow on |x| under forward Euler is: xn+1 = xn − ∆t sign(xn ). The gradient sign(x) has unbounded Lipschitz constant at the origin, so any ∆t > 0 produces a limit cycle of amplitude ∼ ∆t around the minimum; this is bounded but non-convergent. The DPS bias integration is the multidimensional analogue. As Yt approaches the constraint {A(Y ) = y}, the gradient ∇∥A(Y )−y∥2 = ∇A(Y )⊤ (A(Y )− y)/∥A(Y ) − y∥2 does not vanish, while the annealing schedule (15) multiplying the drift cancels the Euler step size exactly (see Appendix G). The forward Euler stability criterion is therefore inevitably violated near the constraint, and the iteration enters a limit cycle of amplitude ∼ σmax (∇A PT M )2 tangent to M. The advantage of ∥ · ∥2 over ∥ · ∥22 as the reward is that the iterates remain semi-stable in the sense of Lyapunov: they settle into a limit cycle at distance utmost α ∥∇A∥op from the constraint manifold. By contrast, when the forward Euler stability criterion is violated for ∥ · ∥22 , the oscillations diverge. Early Guidance Stopping. To avoid these numerical instabilities, practitioners apply early guidance stopping Algorithm 2, terminating the guidance at some intermediate time tstop ∈ [0, T ]. Combining this with the bias result of Theorem 1, we recover the output of the standard DPS algorithm with early guidance stopping. Theorem 2. [Early Guidance Stopping] If guidance is stopped at time tstop = T − t∗ , the output of the DPS algorithm is given by EXt ∼OU wt∗ (Xt∗ ) eηt∗ Ry (x̂t∗ (Xt∗ )) X0 = x νyDP S,t∗ (x) = ρ∗ (x), Z∗ 9
where " wt∗ (x) := EOU
γ(XT −t∗ ) exp → − µ T (XT −t∗ )
Z T −t∗
! S cDP (t∗ + s, Xs ) ds ∗
# X0 = x ,
(17)
0
with
dηt , dt where ηt is annealing schedule (15) and α > 0 is a hyper-parameter. c∗DP S (t, x) = cDP S (t, x) ηt + α∥A(x) − y∥2
See for instance [Huang et al., 2026, Proposition 2.6] for the effect of early guidance stopping in the simpler linear-quadratic case.
6
Acknowledgments
The research of AP and SS has been partially supported by NSF Grants 2019844, 2505865 and 2112471, and the UT Austin Machine Learning Lab. The research of MGD was partially supported by NSF-DMS-2205937.
10
A
Feynman-Kac representations of the Radon-Nikodym derivative
We consider πt′ , πt ∈ C 2 ((0, T ) × Rd ) satisfying a problem of the form (4): ← ν−t (Surrogate − (Algorithm Path) are such an example. Path) and ← µ t Lemma 2 (Feynman-Kac for Density Ratio). Let T > 0 and α ∈ (0, 1). Consider two initial measures, with Lebesgue densities π0 (x)(x), π0′ (x)dx with π0 , π0′ ∈ Cl2+α oc(Rd ). Suppose π0 (x) > 0 for all x ∈ Rd , and the Radon-Nikodym derivative dπt′ /dπt = g0 (x) ∈ C 2 (Rd ) with 2 sub-Gaussian growth |g0 (x)| + |∇g0 (x)| + |D2 g0 (x)| ≤ C eλ|x| for some C, λ. Suppose also ′ ′ that v, v , f, f satisfy, α/2, 1+α
v, v ′ ∈ Cloc
α/2, α
f, f ′ ∈ Cloc
(0, T ) × Rd )
(0, T ) × Rd )
and growth conditions, for some fixed K1 , K2 (independent of t), |vt (x)| + |vt′ (x)| ≤ K1 (1 + |x|),
(t, x) ∈ (0, T ) × Rd
|ft (x)| + |ft′ (x)| ≤ K2 (1 + |x|2 ),
(t, x) ∈ (0, T ) × Rd
. Then we have the following: (i) There exist unique classical solutions π, π ′ ∈ C 1,2 (0, T ) × Rd to ( ∂t πt = −∇ · (vt πt ) + ∆πt + ft πt , π|t=0 = π0 , ∂t πt′ = −∇ · (vt′ πt′ ) + ∆πt′ + ft′ πt′ ,
(18)
π ′ |t=0 = π0′ ,
πt , πt′ > 0 for all t ∈ (0, T ). (ii) The ratio gt (x) := πt′ (x)/πt (x) belongs to C 1,2 (0, T )×Rd and is the unique classical solution to ∂t gt = ∆gt + bt · ∇gt − ct gt , g|t=0 = g0 , (19) with bt (x) := 2 ∇ log πt (x) − vt′ (x) + vt (x), ct (x) := ft′ (x) − ft (x) + ∇ · vt (x) − vt′ (x) . (iii) The ratio admits the Feynman–Kac representation Z t x gt (x) = E g0 (Xt ) exp − ct−s (Xs ) ds ,
(20)
0
where (Xs )s∈[0,t] is the unique strong solution of the SDE √ dXs = bs (Xs ) ds + 2 dWs , X0 = x, and W is a standard d-dimensional Brownian motion. Proof. That the assumptions on the coefficients and measures π0′ , π0 imply (i) is a classical result [Ladyženskaja et al., 1968, Ch. IV, Thm. 5.1]. Consequently, gt is well defined, positive, and C 2 ((0, T ) × Rd ): we seek to prove (ii) and (iii). First, write gt = exp(log gt ): in the broader context of Hamilton-Jacobi-Bellman equations, this is sometimes called the Cole-Hopf transformation. For never-vanishing φ ∈ C 2 (Rd ), this transformation yields the 2 identity ∆φ φ = ∆ log φ + |∇φ| . Together with (4), the Laplacian-identity gives the equations ′ for ∂t log πt and ∂t log πt , taking the difference to obtain ∂t log gt = ∆ log gt + (−vt′ + vt )∇ log gt + |∇ log πt′ |2 − |∇ log πt |2 + (−∇ · vt′ + ∇ · vt + ft′ − ft ) Introduce ∇ log πt to write |∇ log πt′ |2 = |∇ log gt |2 + 2∇ log πt ∇ log gt + |∇ log πt |2 and thus ∂t log gt = ∆ log gt + |∇ log gt |2 + (2∇ log πt − vt′ + vt ) ·∇ log gt + (−∇ · vt′ + ∇ · vt + ft′ − ft ) | {z } | {z } bt (x)
11
−ct (x)
where bt (x) and ct (x) are spatially dependent coefficients. Multiply by gt and apply again t the identity for ∆g gt to conclude ∂t gt = ∆gt + bt (x) · ∇gt − c(t, x)gt (x) π′ g0 (t, x) = 0 (t, x) π0 Now consider consider the SDE dXs = b(s, Xs ) ds +
√
2 dWs ,
X0 = x,
which has generator Ls = ∆ + b(s, ·) · ∇. Fix now t ∈ (0, T ) and define the process, which depends on the whole trajectory, Z s Ms := gt−s (Xs ) exp − c(t − r, Xr ) dr . 0
Applying the Ito formula and the equation for gt (x), the drift term vanishes: √ − R s c(t−r,Xr ) dr dMs = 2 e 0 ∇gt−s (Xs ) · dWs , so M is a martingale (see [Karatzas and Shreve, 1991, Theorem 5.7.6] for standard presentation). At endpoints, the martingale property gives Ex [Mt ] = M0 = gt (x) (20).
B
Bounds and Identities on the OU process
B.1
The effective sample backward path
The starting point for sampling from a prior in score-based generative models is the forward − path t 7→ → ρ t that interpolates between the prior ρ0 = ρ∗ and the Gaussian ρ∞ = N (0, I). This is obtained by solving the OU process ρt = Law(Xt ), where {Xt }t≥0 satisfies the SDE √ dXt = −Xt dt + 2 dBt , X0 ∼ ρ∗ , or, equivalently, the Fokker–Planck equation → − − ∂t − ρ t = ∆→ ρ t + ∇ · (x → ρ t ), → − ρ =ρ . 0
(21)
∗
By the log-Sobolev inequality for ρ∞ and the relative-entropy decay along the OU flow, we have the exponentially decaying bound for t > 1, H(ρt | ρ∞ ) ≤ C e−t ,
(22)
where C is a universal constant independent of dimension. To obtain approximate samples from ρ∗ , we fix a time horizon T and reverse the path: set ← − − ρt =→ ρ T −t . The reverse path satisfies ← − − − − ∂t − ρ t = ∆← ρ t − 2∇ · (∇ log → ρ T −t ← ρ t ) − ∇ · (x ← ρ t ), (23) ← − → − ρ = ρ , 0
T
where we have used the diffusion-to-drift identity − − − ∆→ ρ = −∆← ρ + 2∇ · (∇ log → ρ t
t
← −
T −t ρ t ).
− In practice, the initial condition is replaced by a standard Gaussian ← ρ eff 0 = N (0, I). The computationally intensive part of this strategy is obtaining a good approximation sθ (t, x) ≈ − ∇ log → ρ t (x) from samples; see Section 2. The effective samples are then obtained by approximating the solution of ← ← −eff ← −eff ← −eff ∂t − ρ eff t = ∆ ρ t − 2∇ · (sθ (T − t, ·) ρ t ) − ∇ · (x ρ t ), (24) ← − ρ eff = N (0, I), 0
12
which arises as the law of √ dX̃t = X̃t dt + 2 sθ (T − t, X̃t ) dt + 2 dBt , X̃0 ∼ N (0, I). − ← − d H(← Differentiating the relative entropy dt ρ eff t | ρ t ) along the flow yields the relative-entropy ← − bound between the effective samples X̃T ∼ ρ eff T and the true distribution X0 ∼ ρ∗ , Z TZ − → − ∥∇ log ρt (x) − sθ (t, x)∥2 ρt (x) dx dt . + H(← ρ eff | ρ ) ≤ H( ρ | ρ ) (25) ∗ T ∞ T {z } | d 0 R | {z } ≤ Ce−T Score-approximation error
In what follows, we assume the score-approximation error is negligible, so the samples obtained by denoising are for practical purposes indistinguishable from ρ∗ . B.2
Tweedie’s identities
A widely used heuristic is to estimate X0 from a noisy observation Xt via the conditional expectation x̂t (x) := E[X0 | Xt = x], (26) which, under the negligible-score-approximation assumption of Section B.1, also equals E[X̃T | X̃T −t = x] for the effective backward process. Tweedie’s formula expresses (26) in closed form via the score: x̂t (x) = et x + (et − e−t ) ∇ log ρt (x). Derivation.
(27)
The OU semigroup admits the Gaussian kernel 1 ∥x − e−t x0 ∥2 ρt|0 (x | x0 ) = exp − , 2(1 − e−2t ) (2π(1 − e−2t ))d/2
so ∇x log ρt|0 (x | x0 ) = −(x − e−t x0 )/(1 − e−2t ). x0 ) ρ∗ (x0 ) dx0 in x and dividing by ρt (x), ∇ log ρt (x) = −
Differentiating ρt (x) =
(28) R
ρt|0 (x |
x − e−t x̂t (x) . 1 − e−2t
Solving for x̂t (x) and using et (1 − e−2t ) = et − e−t yields (27).
□
Second-order identity. Differentiating (27) in x, ∇x̂t (x) = et I + (et − e−t ) ∇2 log ρt (x).
(29)
The right-hand side admits a probabilistic interpretation as a rescaled conditional covariance: Σt (x) := Cov(X0 | Xt = x) = (et − e−t ) ∇x̂t (x).
(30)
Proof of (30). Differentiating ρt twice via (28) and subtracting (∇ log ρt )(∇ log ρt )⊤ to convert from ∇2 ρt /ρt to ∇2 log ρt , ∇2 log ρt (x) = −
I e−2t + Σt (x). 1 − e−2t (1 − e−2t )2
(31)
Substituting into (29), the identity (et − e−t )/(1 − e−2t ) = et cancels the et I contribution, and the identity (et − e−t ) e−2t /(1 − e−2t )2 = e−t /(1 − e−2t ) = 1/(et − e−t ) collapses the covariance term, giving Σt (x) ∇x̂t (x) = t , e − e−t which rearranges to (30). □ 13
B.2.1
The zero-noise limit
Throughout this subsection we assume ρ∗ is supported on a smooth, compact, k-dimensional submanifold M ⊂ Rd of positive reach τM > 0, with a smooth positive density with respect to the volume measure on M. The orthogonal projection PM : x 7→ arg minx0 ∈M ∥x−x0 ∥ is then well-defined and smooth on the tubular neighborhood NτM := {x ∈ Rd : dist(x, M) < τM }, and the case x ∈ M corresponds to dist(x, M) = 0. We compute the small-t behavior of x̂t and ∇x̂t via Laplace’s method on M; standard references include Varadhan [1966], Dembo and Zeitouni [2010]. Setup. As t → 0+ , 1 − e−2t = 2t + O(t2 ) and e−t = 1 + O(t), so the OU kernel (28) concentrates as 1 + O(t) ∥x − x0 ∥2 ρt|0 (x | x0 ) = exp − + O(1) as t → 0+ . 4t (4πt)d/2 The integrals defining x̂t (x) are of Laplace form on M with phase Φ(x0 ) = ∥x − x0 ∥2 at temperature 4t. Local geometry. Fix x ∈ NτM , set p := PM (x), and let T := T Mp with PT the orthogonal projection onto T . Parametrize M near p by tangent vectors, φ(v) = p + v + O(∥v∥2 ),
v ∈ T,
where the O(∥v∥2 ) correction lies in the normal space N := T ⊥ and is bounded by the second fundamental form. Since x − p ∈ N is orthogonal to v ∈ T , ∥x − φ(v)∥2 = ∥x − p∥2 + ∥v∥2 + O(∥v∥3 ),
(32)
2
so the phase Φ(x0 ) = ∥x − x0 ∥ is locally quadratic in tangent coordinates with Hessian 2 IT at the minimizer p. Limit of x̂t . On NτM , the phase Φ|M has unique global minimizer p with Φ(p) = ∥x − p∥2 . Laplace’s method on M, applied with the local expansion (32), gives the asymptotic Z f (x0 ) e−Φ(x0 )/(4t) ρ∗ (x0 ) dvol(x0 ) = (4πt)k/2 e−Φ(p)/(4t) ρ∗ (p) f (p) + O(t) M
for any C 1 function f on M. Recall that the Tweedie estimate is the ratio R x e−Φ(x0 )/(4t) ρ∗ (x0 ) dvol(x0 ) R 0 . x̂t (x) = M e−Φ(x0 )/(4t) ρ∗ (x0 ) dvol(x0 ) M Applying the asymptotic to numerator and denominator, the common prefactor (4πt)k/2 e−Φ(p)/(4t) ρ∗ (p) cancels, leaving x̂t (x) √ = p + O(t), where the O(t) remainder collects the next-order Laplace corrections. The O( t) contributions from tangential fluctuations vanish by Gaussian symmetry on T , since the linear function v 7→ v has zero mean under a centered Gaussian. Hence x̂t (x) = PM (x) + O(t), Limit of ∇x̂t .
lim x̂t (x) = PM (x).
t→0+
(33)
By (30) and et − e−t = 2t + O(t3 ), ∇x̂t (x) =
Σt (x) t e − e−t
=
Σt (x) 1 + O(t2 ) , 2t
(34)
so it suffices to compute Σt (x) to leading order. By (32), the conditional law ρ0|t is asymptotically Gaussian on T with covariance 2t IT , and the normal component of X0 − p is of order ∥v∥2 = O(t) and contributes only at order t2 . Hence Σt (x) = 2t PT + O(t2 ),
(35)
lim ∇x̂t (x) = PT MPM (x) .
(36)
and substituting into (34) yields t→0+
14
Consequence. In the small-noise regime relevant near the constraint set, x̂t acts as the orthogonal projection onto the data manifold and ∇x̂t acts as the projection onto the corresponding tangent space. This is the geometric structure that drives the manifold-tangent oscillations of the DPS guidance analyzed in Section 5.
C
The DPS Algorithm [Chung et al., 2024]
Algorithm 1 DPS N Require: N , y, {ζi }N i=1 , {σ̃i }i=1 1: xN ∼ N (0, I) 2: for i = N − 1 to 0 do 3: ŝ ← sθ (xi , i) 4: x̂0 ← √1ᾱ (xi + (1 − ᾱi )ŝ) i 5: z ∼ N (0, I) √ √ α (1−ᾱ
)
ᾱi−1 βi
′ i i−1 xi−1 ← xi + 1−ᾱi x̂0 + σ̃i z 1−ᾱi ′ 7: xi−1 ← xi−1 − ζi ∇xi ∥y − A(x)∥22 8: end for 9: return x̂0
6:
Here A : Rd → RL is a general linear or non-linear observation operator and y ∈ RL is the observation. The bias schedule {ζi }N i=1 is taken to be trajectory-dependent, α ζi = , (37) ∥y − A(x)∥2 where α ∈ [0.2, 1] is a hyperparameter chosen depending on the inverse problem to be solved. In practice, the choice of bias schedule significantly affects the performance of the algorithm.
D
Proof of Theorem 1
In this section, we will prove the following theorem S Theorem 1. The terminal law νyDP S := ← ν−DP of the DPS-SDE (DPS SDE) differs from T the true posterior µy by a pointwise multiplicative weight: µy (x) = ω(x) νyDP S (x).
(9)
The weight ω admits two equivalent Feynman–Kac representations in terms of the reaction term cDP S defined in (8): (i) Backward path (condition on the DPS denoising process arriving at YT = x): " # Z T → − µ T (Y0 ) ω(x) = EY ∼(DPS SDE) exp − cDP S (T − s, Ys ) ds YT = x . (10) γ(Y0 ) 0 (ii) Forward path (condition on the OU process (2) starting at X0 = x): " # Z T γ(XT ) 1 = EX∼OU → exp cDP S (s, Xs ) ds X0 = x . − ω(x) µ T (XT ) 0
(11)
Both path functionals are expressible in terms of quantities obtainable from the score oracle and its Jacobian via Tweedie’s formula. Importance-weighting DPS samples by ω recovers µy exactly. Sketch. The proof has three steps (elaborated below). Step 1. We record an evolution equation satisfied by any tilted prior path t 7→ µt := ht ρt /Zt . ∂t µt = ∆(µt ) + ∇ · (xµt ) + c[h](t, x) − ∂t log Zt µt , 15
for an appropriate choice of c[h] (Lemma 3). Step 2. Specializing to the DPS tilt ht = eRy ◦x̂t and using the Kolmogorov backward equation for the conditional mean x̂t (x) = E[X0 | Xt = x], we identify the reaction term as exactly cDP S from (8) (Lemma 4). Step 3. In Lemma 5, we use Lemmas 3 and 4 and apply the Feynman-Kac formula Lemma 2 µ to the resulting evolution equation to obtain an expression for the weights µDPy S . Reversing the path integral, in Lemma 6 we recast the path expectation along the DPS reverse SDE (DPS SDE) with the forward OU path measure. D.1
Two Tweedie identities and the Kolmogorov backward equation
We collect three facts used repeatedly. Throughout, ρt denotes the marginal of the OU forward process (2) starting at X0 ∼ ρ∗ , and x̂t , Σt are the conditional mean and covariance of X0 given Xt . (F1) First-order Tweedie. A direct integration of the OU semigroup yields x̂t (x) = et x + (et − e−t )∇ log ρt (x).
(38)
(F2) Second-order Tweedie. Differentiating (38) in x and using the standard identity Σt (x) = (et − e−t )2 D2 log ρt (x) + (et − e−t )et I (which follows from a second-order expansion of the OU posterior, or equivalently from differentiating Tweedie under Bayes’ rule): ∇x̂t (x) =
1 et − e−t
Σt (x).
(39)
In particular, ∇x̂t is symmetric and positive semidefinite. (F3) Kolmogorov backward equation for x̂t . By Anderson’s reversal [Anderson, 1982b], the time-reversed OU process X̃s := XT −s satisfies the SDE √ dX̃s = (X̃s + 2∇ log ρT −s (X̃s )) ds + 2 dB̃s , (40) with generator L̃t f := ∆f + (x + 2∇ log ρt (x)) · ∇f . Since x̂t (x) = E[X0 | Xt = x] = E[X̃T | X̃T −t = x], x̂t is a Kolmogorov backward solution along X̃, and hence ∂t x̂t (x) = ∆x̂t (x) + ∇x̂t (x) · x + 2∇ log ρt (x) . D.2
(41)
Step 1: Evolution of tilted prior paths
Lemma 3 (Tilted-prior Fokker-Planck). Let h ∈ C 1,2 ([0, T ] × Rd ) be positive with ht ∈ L1 (ρt dx), and define the tilted prior path Z ht (x)ρt (x) πt (x) := , Zt := ht (x)ρt (x) dx. Zt Then πt obeys
∂t πt = L† πt + c[h](t, x) − ∂t log Zt πt ,
(42)
†
where L π := ∆π + ∇ · (xπ) is the OU Fokker-Planck operator and the reaction term is c[h](t, x) := ∂t log ht − 2∇ log ρt + x · ∇ log ht − |∇ log ht |2 − ∆ log ht . (43) Proof. Using the product rule and writing Zt terms as log-derivatives: ∂t πt =
(∂t ht )ρt ht (∂t ρt ) Żt + − πt . Zt Zt Zt
Writing ∂t ht = ht ∂t log ht and substituting the OU Fokker–Planck equation ∂t ρt = L† ρt : ∂t πt = (∂t log ht ) πt +
ht † L ρt − (∂t log Zt ) πt . Zt
16
The OU Fokker–Planck operator is L† ρ = ∆ρ+∇·(xρ). Using ∆ρt = ρt (|∇ log ρt |2 +∆ log ρt ) and ∇ · (xρt ) = ρt (d + x · ∇ log ρt ): h i ht L† ρt = ht ρt |∇ log ρt |2 + ∆ log ρt + d + x · ∇ log ρt . Set φ = ht /Zt so that πt = φρt . We compute L† πt = L† (φρt ) via the product rule applied to each term of L† = ∆ + ∇ · (x · ): ∆(φρt ) = φ∆ρt + 2∇φ · ∇ρt + ρt ∆φ, ∇ · (xφρt ) = φ ∇ · (xρt ) + ρt x · ∇φ. Summing, we have L† (φρt ) = φ L† ρt + 2∇φ · ∇ρt + ρt ∆φ + ρt x · ∇φ. Rearranging to isolate φL† ρt : φ L† ρt = L† (φρt ) − 2∇φ · ∇ρt − ρt (∆φ + x · ∇φ) . Since Zt does not depend on x, we have ∇φ = ∇ht /Zt and ∆φ = ∆ht /Zt . Substituting φρt = πt , ∇ρt = ρt ∇ log ρt , and ρt /Zt = πt /ht , and using the identities ∇ht /ht = ∇ log ht and ∆ht /ht = |∇ log ht |2 + ∆ log ht : h i ht † L ρt = L† πt − πt 2∇ log ρt · ∇ log ht + |∇ log ht |2 + ∆ log ht + x · ∇ log ht . Zt Collecting all πt terms: ∂t πt = L† πt + (c[h](t, x) − ∂t log Zt ) πt , where the reaction coefficient is c[h](t, x) = ∂t log ht − (2∇ log ρt + x) · ∇ log ht − |∇ log ht |2 − ∆ log ht
D.3
Step 2: Specialization to the DPS tilt
− Lemma 4 (Reaction term for the DPS tilt). With ht (x) = exp(Ry (x̂t (x))), so that πt = → µt in (DPS Surrogate path), the reaction c[h] from (43) reduces to h i 1 2 2 c[h](t, x) = − t tr Σ (x) D R (x̂ (x)) Σ (x) + |Σ (x)∇R (x̂ (x))| . (44) t y t t t y t (e − e−t )2 Equivalently, c[h](t, x) − ∂t log Zt = cDP S (t, x) with 1 2 2 d cDP S (t, x) = − tr Σt (x)(D Ry )(x̂t (x))Σt (x) + Σt (x)∇Ry (x̂t (x)) − dt log Zt (et − e−t )2 . Proof. Set log ht = Ry ◦ x̂t and apply the chain rule component-wise: ∂t log ht = ∇Ry (x̂t ) · ∂t x̂t , ∇ log ht = (∇x̂t ) ∇Ry (x̂t ), ∆ log ht = ∇Ry (x̂t ) · ∆x̂t + tr ∇x̂t D2 Ry (x̂t ) ∇x̂t , using symmetry of ∇x̂t from (39). Substituting into (43) and grouping terms by their dependence on ∇Ry and D2 Ry , we get the following expression for c[h]: c[h](t, x) = ∇Ry (x̂t ) · ∂t x̂t − (2∇ log ρt + x) · (∇x̂t ) ∇Ry (x̂t ) − |(∇x̂t ) ∇Ry (x̂t )|2 − ∇Ry (x̂t ) · ∆x̂t + tr ∇x̂t D2 Ry (x̂t ) ∇x̂t Grouping the terms of ∇R together, we have c[h](t, x) = ∇Ry (x̂t ) · [∂t x̂t − ∆x̂t − (2∇ log ρt + x) · (∇x̂t )] − |(∇x̂t ) ∇Ry (x̂t )|2 + tr ∇x̂t D2 Ry (x̂t ) ∇x̂t By the Kolmogorov backward equation (41), we have ∂t x̂t − ∇x̂t (x + 2∇ log ρt ) − ∆x̂t = 0. This is the key cancellation underlying the bias formula. The remaining contributions are −|∇ log ht |2 and the trace piece of −∆ log ht : −|(∇x̂t )∇Ry (x̂t )|2 − tr ∇x̂t D2 Ry (x̂t ) ∇x̂t . Using (39) to substitute ∇x̂t = Σt /(et − e−t ) and pulling out the common scalar factor yields (44). 17
D.4
Step 3: Feynman-Kac two ways
We now combine Lemmas 3 and 4 to prove Theorem 1(ii). S Lemma 5. The terminal law νyDP S := ← ν−DP of the DPS-SDE (DPS SDE) differs from the T true posterior µy by a pointwise multiplicative weight: µy (x) = ω(x) νyDP S (x). The weight ω admits a Feynman–Kac representations in terms of the reaction term cDP S defined in (8): " # Z T → − µ T (Y0 ) YT = x . ω(x) = EY ∼(DPS SDE) cDP S (T − s, Ys ) ds exp − γ(Y0 ) 0 − := → − Proof. Setting ← µ µ T −t and applying (42) together with the Anderson identity t − = ∆← − − ∇· (x + 2∇ log ← − )← − −L† ← µ µ µ µ t
t
t
t
gives exactly ← − = ∆← − − 2 ∇· ∇ log → − − − ∇· (x ← − )−c ← − ∂t µ µ µ T −t ← µ µ t t t t DP S (t, x) µ t , ← − → − µ 0 (x) = µ T (x),
(45)
− = with reaction −cDP S by Lemma 4, which is (DPS Surrogate PDE). By construction, ← µ t → − µ 0 = µy . Equation (45) is the Fokker-Planck equation associated with the DPS reverse SDE ( √ 2 2 dBt , dYt = Yt + 2∇ log ρT −t (Yt ) + et −e −t ΣT −t (Yt ) ∇Ry (x̂T −t (Yt )) dt + (46) Y0 ∼ γ, augmented by a multiplicative reaction −cDP S . The DPS algorithm itself directly simulates S (46) from Y0 ∼ γ, that is, without the source. This produces a marginal ← ν−t with ← ν−T = µDP . y Effectively, the DPS algorithm solves ( − ∂t ← ν−t = ∆← ν−t − 2∇· (∇ log → µ T −t ← ν−t ) − ∇· x ← ν−t , (47) S µDP (x) = γ(x), 0 Two operators (45) and (47) differ only by a multiplicative reaction and an initial condition. They can be related by a Feynman-Kac formula of Lemma 2. The ground truth satisfies − (x), for any test function φ(x) we have µy (x) = ← µ t Z Z φ(x)µy (x) dx = φ(x)νTGT (x) dx Rd Rd R → − µ T (Y0 ) − T cDP S (T −s,Ys ) ds e 0 = EY ∼(46) φ(YT ) γ(Y0 ) → R − µ T (Y0 ) − T cDP S (T −s,Ys ) ds = EY ∼(46) φ(YT ) E e 0 |YT γ(Y0 ) | {z } w(YT )
Z = Rd
S φ(x)w(x)µDP (x) dx, y
where we conditioned on the value of YT , used the law of total expectation, and the observation YT ∼ ← ν−T . Concretely, → R ← − (x) − µy (x) µ µ T (Y0 ) − T cDP S (T −s,Ys ) ds t 0 = = E e Y = x . (48) T Y ∼(46) ← −DP S (x) S (x) µDP γ(Y0 ) µ y T − − −DP S = The boundary factor → µ T /γ accounts for the mismatch between ← ν−GT =→ µ T and ← µ 0 0 γ. 18
The formula for the weight using OU, can now be derived using Anderson the time reversal of SDEs. For clarity and brevity, we instead provide a derivation using PDE satisfied by the ratio of the algorithmic path to the PDE for the ratio of the algorithmic path to the surrogate path. As in Surrogate Path and Algorithm Path we here define ← − := ← −DP S = 1 exp(R ◦ x̂ µ µ t y T −t (x))ρT −t (x) t ZT −t ← ν−t := ν DP S = Law(Yt ). and for 0 ≤ t < T , as in A define the density ratio ψt (x), − =← ψt ← µ ν−t t 1 noting limt→T ψt = ψT = ω(x) as in (11).
(49)
Lemma 6. The ratio ψt (x) solves the parabolic initial value problem, ( S ∂t ψt = ∆ψt − x∇ψt + cDP T −t (x)ψt γ(x) , ψ0 (x) = ← − ν GT (x) 0
and thus, by Feynman-Kac formula, with t < T , Z t γ(Xt ) ψt (x) := EOU → cT −t+s (Xs )ds X0 = x exp − µ T (Xt ) 0
(50)
1 := ψT (x) as the ratio. Taking the limit t → T yields the statement of Theorem 1 with ω(x)
− (x) and Proof. The following calculations are justified classically, since for t < T , both ← µ t ← − ← − ← − DP S νt (x) are smooth, positive densities. Note log ψt = log ν t − log µ t , using (45) and (47) we have − ∂t log ← ν−t = ∆ log ← ν−t + |∇ log ← ν−t |2 − d − 2∆ log ← µ t ← − ← − ← − − x∇ log ν − 2∇ log µ ∇ log ν t
t
t
− = ∆ log ← − + |∇ log ← − |2 − d − 2∆ log ← − ∂t log ← µ µ µ µ t t t t ← − ← − 2 DP S − x∇ log µ − 2|∇ log µ | − c t
t
T −t
Taking the difference and completing the square shows − ∂t log ψt = ∂t log ← ν−t − ∂t log ← µ t S = ∆ log ψt + |∇ log ψt |2 − x∇ log ψt + cDP T −t exp(·)
(51)
Applying the Cole-Hopf transformation (log ψt 7→ ψt ) obtains (51). Finally, applying Feynman-Kac to an initial-value problem (rather than terminal-value) induces time-reversal S of the multiplier cDP T −t+s (Xs ) in the path-functional.
E
Early Guidance Stopping, Proof of Theorem 2
Algorithm 2 DPS with Early Guidance Stopping N Require: istop , N , y, {ζi }N i=1 , {σ̃i }i=1 1: xN ∼ N (0, I) 2: for i = N − 1 to 0 do 3: ŝ ← sθ (xi , i) 4: x̂0 ← √1ᾱ xi + (1 − ᾱi )ŝ i 5: z ∼ N (0, I) √ √ α (1−ᾱ
)
ᾱi−1 βi
′ i i−1 6: xi−1 ← xi + 1−ᾱi x̂0 + σ̃i z 1−ᾱi 7: if i > istop then ′ 8: xi−1 ← xi−1 − ζi ∇xi ∥y − A(xi )∥22 9: else ′ 10: xi−1 ← xi−1 11: end if 12: end for 13: return x̂0
19
The output Algorithm 2, up to discretization error, is characterized in the following result. Theorem 2. [Early Guidance Stopping] If guidance is stopped at time tstop = T − t∗ , the output of the DPS algorithm is given by EXt ∼OU wt∗ (Xt∗ ) eηt∗ Ry (x̂t∗ (Xt∗ )) X0 = x νyDP S,t∗ (x) = ρ∗ (x), Z∗ where " wt∗ (x) := EOU
γ(XT −t∗ ) exp → − µ T (XT −t∗ )
Z T −t∗
#
! S cDP (t∗ + s, Xs ) ds ∗
X0 = x ,
(17)
0
with
dηt , dt where ηt is annealing schedule (15) and α > 0 is a hyper-parameter. c∗DP S (t, x) = cDP S (t, x) ηt + α∥A(x) − y∥2
Proof. We apply (Surrogate Path) to the annealed path (14): e − t 7→ → µt =
−αηt ∥A(x̂t (x))−y∥2
Zt
By Lemma 3, the reverse equation reads − = ∆← − − ∇ · (x ← − ) − 2∇ · (∇ log → − ∂← µ µ µ ρ t
t
ρt
.
← − − d µ t, ηT −t R ◦ x̂T −t ← − αηT −t cDP S − α dt {z } | t
t
← −
T −t µ t ) + 2αηT −t ∇ · (∇(R ◦ x̂T −t ) µ t )
S (T −t) cDP ∗
where R(x) = ∥A(x) − y∥2 for the (possibly non-linear) observation operator A. The corresponding algorithmic SDE (Algorithm Path) with early stopping at time tstop := T − t∗ is ← − ∆ ν−t − ∇ · (x ← ν−t ) − 2∇ · (∇ log → ρ T −t ← ν−t ) + 2αηT −t ∇ · (∇(R ◦ x̂T −t ) ← ν−t ), t ∈ (0, tstop ), ∂t ← ν−t = ← − ← − → − ← − ∆ ν − ∇ · (x ν ) − 2∇ · (∇ log ρ ν ), t ∈ [t , T ). t
T −t
t
t
stop
(52) − and ← On (0, tstop ), both ← µ ν−t satisfy the same equation, so Theorem 1 applies directly and t yields # ! " Z T −t∗ ← ν−tstop (x) γ(XT −t∗ ) DP S X0 = x . = ωt∗ (x) = EOU → exp c∗ (t∗ + s, Xs ) ds ← − − µ µ T (XT −t∗ ) tstop (x) 0 − Substituting the explicit form of ← µ tstop , ← − − αηt∗ R◦x̂t∗ (x) ← ν−T −t∗ (x) = ωt∗ (x) ← µ ρ T −t∗ (x). tstop (x) = ωt∗ (x) e
(53)
On [tstop , T ), the SDE for ← ν−t is the unbiased reverse OU equation, started from (53). Setting s = T − t for the corresponding forward time, the Radon–Nikodym derivative of the initial − condition with respect to the OU forward marginal → ρ t∗ is d← ν−T −t∗ gt∗ (x) := ← (x) = ωt∗ (x) eαηt∗ R(x̂t∗ (x)) . d− ρ T −t∗
Applying the Feynman–Kac identity for ratios (Lemma 2), h i d← ν−T d← ν−T (x0 ) = ← (x0 ) = E gt∗ (Xt∗ ) X0 = x0 , − dρ∗ dρT
(54)
where {Xs }s≥0 is the OU process started from X0 ∼ ρ∗ . Substituting the expression for gt∗ yields the claim. 20
F
Time discretization of DDPM
To establish the correspondence between the discrete variance schedule βi used in Denoising Diffusion Probabilistic Models (DDPM) Ho et al. [2020] and the continuous time steps ∆ti of the underlying Ornstein-Uhlenbeck (OU) process, we compare their respective transition kernels. The forward Markov jump process in DDPM defines the transition from step i to i + 1 as: p q(xi+1 |xi ) = N (xi+1 ; 1 − βi xi , βi I) (55) The continuous-time reverse SDE under consideration is given by: √ dXt = −Xt dt + 2dBt (56) For a finite time increment ∆ti , the exact solution to this SDE yields the transition: (57) p(xt+∆ti |xt ) = N (xt+∆ti ; e−∆ti xt , (1 − e−2∆ti )I) For the discrete Markov chain to exactly discretize the continuous SDE, the coefficients of the mean and variance must be consistent across regimes: p (58) e−∆ti = 1 − βi 1 − e−2∆ti = βi . Solving for ∆ti we obtain 1 ∆ti = − ln(1 − βi ) (59) 2 The linear noise schedule of DDPM is given by βmax − βmin βi = βmin + i ≈ 10−4 + 2i 10−5 for i = 1, . . . , 1000, N with the choices βmin = 10−4 , βmax = 0.02, and N = 1000 steps. Applying the first-order Taylor expansion ln(1 − ϵ) ≈ −ϵ, we obtain the approximately linear relationship between ∆ti and βi : 1 ∆ti ≈ βi (60) 4 Next, we derive a relationship in time between the discrete steps ti and the continuous time t by summing over the increments: i i i X 1 1X 1X i(i + 1) βj = 10−4 + 2j 10−5 = ti = ∆tj ≈ 10−4 i + 2 · 10−5 . (61) 4 j=1 4 j=1 4 2 j=1 Next, we solve for a functioni(t) that maps continuous time to discrete steps by inverting
= t: the quadratic relationship 41 10−4 i + 2 · 10−5 i(i+1) 2 √ 121 + 16t 105 − 11 i(t) = (62) 2 This function i(t) provides a mapping from continuous time t to the corresponding discrete step index i in the DDPM framework, allowing us to understand the time step behavior in the continuum limit. Substituting i(t) into (60), we obtain 1 1 ∆t(t) ≈ βi(t) = 10−4 + 2 · 10−5 i(t) 4 4 p 1 −4 = 10 + 10−5 121 + 16 · 105 t − 11 4 10−5 p = 121 + 16 · 105 t − 1 . (63) 4 √ A naive large-t approximation ∆t(t) ∼ 10−5 t would incorrectly vanish at t = 0. To preserve the nonzero constant floor at the origin, we drop the small −1 term (negligible compared to √ −5 121 = 11) but keep the constant 121 inside the square root. Pulling the prefactor 104 inside the radical yields the compact form q √ (64) ∆t(t) ≈ 10−5 t + ∆t20 ≈ 3 · 10−5 + 3 · 10−3 t. 21
G
Forward Euler instability
In terms of implementation, DPS Algorithm 1 integrates the different terms of (DPS SDE) in √ different ways. The denoising step, corresponding to the terms Yt + 2∇ log ρT −t (Yt ) + 2 dBt , is integrated implicitly via DDPM in Step 6, avoiding numerical instabilities. The bias term 2 αηt et −e −t ΣT −t (Yt ) ∇Ry (x̂T −t (Yt )), however, is integrated explicitly via forward Euler in Step 7. The annealing schedule ηt = 1/∆t(t) in (15) is an auxiliary quantity we introduce to compensate for the missing time step ∆t in Step 7. In place of ∆t, the algorithm multiplies by the path-dependent factor α , ζi = ∥A(x) − y∥2 when the reward is Ry (x) = ∥A(x) − y∥22 . In our analysis, this is equivalent to using the modified reward Ryeff (x) = 2∥A(x) − y∥2 together with the annealing schedule (15). To derive the oscillations, we assume that the prior ρ∗ is a smooth distribution supported on a smooth lower-dimensional manifold M ⊂ Rd embedded in the ambient space. In the limit t → 0+ , x̂t (x) → PM (x), ∇x̂t (x) → PT MPM (x) , where PM denotes the orthogonal projection onto M and PT MPM (x) denotes the orthogonal projection onto the tangent space at PM (x), see Section B.2.1. Consequently, as t → 0+ , the bias guidance is well-approximated by a flow on M, dYt = −
1 PT MYt ∇∥AYt − y∥2 dt, ∆t0
where PT MYt denotes the orthogonal projection onto the tangent space at Yt . Taking a forward Euler step of size ∆t0 , the prefactor 1/∆t0 cancels the step size exactly, so one Euler step corresponds to one full projected-gradient step on M, independently of ∆t0 . Local Lipschitz constant. The local Lipschitz constant of the projected drift on M scales as 2 1 σmax (A PT MY ) Llip ∼ , (65) ∆t0 ∥AY − y∥2 where the residual in the denominator originates from the gradient of the unsquared norm, which renormalizes to unit magnitude as Y approaches the constraint. Stability criterion and inevitable oscillations. Forward Euler stability requires ∆t0 · Llip ≤ 2, i.e., 2 σmax (A PT MY ) ≤ 2 ∥AY − y∥2 . (66) As Y approaches the constraint set {Y : AY = y}, the right-hand side tends to zero while the left-hand side depends only on A and the local geometry of M. The criterion is therefore inevitably violated near the constraint. This is the standard pathology of forward Euler applied to the unsquared norm: the gradient does not vanish as the residual shrinks, but merely renormalizes to unit magnitude along A⊤ (AY − y)/∥AY − y∥2 . The iteration overshoots and oscillates around the constraint, and no choice of step size can restore stability/convergence. We note that these oscillations occur only parallel to the data manifold. Implicit integration of the bias drift, as in Rout et al. [2025], avoids these numerical instabilities.
H
Empirical Evidence of Instability
Conditional Guidance for MNIST digits We consider the setting of posterior sampling with the MNIST prior. This dataset consists of paired images of handwritten digits and their corresponding labels, denoted (x, y). We train a simple MLP classifier softmax(f (x)) ≈ 1y over the MNIST dataset, and set the reward to be R(x) = ∥f (x) − 1k ∥ for some fixed target k. We run DPS with guidance schedule constant at 0.1. The evolution of ∥f (xt ) − 1k ∥ as well as (f (xt ) − 1k ) · 1 is plotted below. 22
Figure 4: We plot a projected discrepancy (f (xt ) − 1k ) · 1 ((Columns 1 and 3)) and the reward ∥f (xt ) − 1k ∥ ((Columns 2 and 4)) across t where denoising proceeds from left (most noise) to right (least noise). The second row depicts a close-up plot of just the last 10 steps to highlight the oscillations. (Columns 1 and 2) are run with a constant guidance schedule Algorithm 1, while (Column 3 and 4) are run with early guidance stopping Algorithm 2 with parameter istop = 100.
Figure 5: Top Are plots associated to the standard DPS algorithm Algorithm 1, Top Left We plot αt = (Pt:k δt ) · (Pt:k δt−1 ) along the DPS trajectories for a constant guidance schedule ζ = 0.1. Top Middle A close-up of steps 525 → 500. Top Right A close-up of steps 25 → 0. Note that αt is close to 0 at the intermediate noise levels, but drops to ≈ −1 towards the low noise levels. Bottom Are the same plots associated to Algorithm 2 with istop = 100. We observe that αt remains close to 0 both at intermediate noise levels and low noise levels. We see a distinct oscillatory pattern is sustained throughout the trajectory when the guidance schedule is constant Algorithm 1. Turning off the guidance schedule Algorithm 2 at time-step istop = 100 eliminates the oscillations in that period, though now the reward is not pulled toward 0. This indicates that the instability is associated with the reward guidance. We see in either case that a softmax applied to the logits of the classifier results in a very high confidence prediction of the correct class despite these oscillations. We also plot the alignment between consecutive steps of the algorithm δt = xt − xt−1 . Because these vectors lie in 784 dimensions, to emphasize the step over step alignment we maintain a subspace described by the most recent ℓ = 50 such steps. In particular, let Pt = [δt−ℓ+1 , δt−ℓ+2 , · · · δt ] ∈ R784×ℓ , and let Pt:k denote just the projection onto the top k principle axis. We plot αt = (Ps:k δs−1 ) · (Ps:k δs ) over the trajectory in Fig. 5. For a purely oscillating trajectory, we expect δt ≈ −δt−1 , resulting in αt ≈ −1. When the δt is “unrelated” to δt , we expect αt ≈ 0. All experiments were run in a few minutes on a single NVIDIA H100 GPU.
23
References Brian D. O. Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Applications, 12:313–326, 1982a. URL https://api.semanticscholar.org/CorpusID: 3897405. Brian D.O. Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Applications, 12(3):313–326, 1982b. ISSN 0304-4149. doi: https://doi.org/10.1016/ 0304-4149(82)90051-5. URL https://www.sciencedirect.com/science/article/pii/ 0304414982900515. Gautham Govind Anil, Shaan Ul Haque, Nithish Kannen, Dheeraj Nagaraj, Sanjay Shakkottai, and Karthikeyan Shanmugam. Fine-tuning diffusion models via intermediate distribution shaping, 2026. URL https://arxiv.org/abs/2510.02692. Benjamin Boys, Mark Girolami, Jakiw Pidstrigach, Sebastian Reich, Alan Mosca, and O. Deniz Akyildiz. Tweedie moment projected diffusions for inverse problems, 2024. URL https://arxiv.org/abs/2310.06721. Joan Bruna and Jiequn Han. Posterior sampling with denoising oracles via tilted transport, 2024. URL https://arxiv.org/abs/2407.00745. Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R. Zhang. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions, 2023. URL https://arxiv.org/abs/2209.11215. Hyungjin Chung, Byeongsu Sim, Dohoon Ryu, and Jong Chul Ye. Improving diffusion models for inverse problems using manifold constraints. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL https://openreview.net/forum?id=nJJjv0JDJju. Hyungjin Chung, Jeongsol Kim, Michael T Mccann, Marc L Klasky, and Jong Chul Ye. Diffusion posterior sampling for general inverse problems. In International Conference on Learning Representations, 2023. Hyungjin Chung, Jeongsol Kim, Michael T. Mccann, Marc L. Klasky, and Jong Chul Ye. Diffusion posterior sampling for general noisy inverse problems, 2024. URL https: //arxiv.org/abs/2209.14687. Giannis Daras, Hyungjin Chung, Chieh-Hsin Lai, Yuki Mitsufuji, Jong Chul Ye, Peyman Milanfar, Alexandros G Dimakis, and Mauricio Delbracio. A survey on diffusion models for inverse problems. arXiv preprint arXiv:2410.00083, 2024. Amir Dembo and Ofer Zeitouni. Large Deviations Techniques and Applications, volume 38 of Stochastic Modelling and Applied Probability. Springer, Berlin, Heidelberg, 2nd edition, 2010. doi: 10.1007/978-3-642-03311-7. Corrected reprint of the second (1998) edition. Prafulla Dhariwal and Alex Nichol. Diffusion models beat gans on image synthesis. 2021. URL https://arxiv.org/abs/2105.05233. Zehao Dou and Yang Song. Diffusion posterior sampling for linear inverse problem solving: A filtering perspective. In The Twelfth International Conference on Learning Representations, 2024. URL https://openreview.net/forum?id=tplXNcHZs1. Zhengyi Guo, Wenpin Tang, and Renyuan Xu. Conditional diffusion guidance under hard constraint: A stochastic analysis approach. arXiv preprint arXiv:2602.05533, 2026. Shivam Gupta, Ajil Jalal, Aditya Parulekar, Eric Price, and Zhiyang Xun. Diffusion posterior sampling is computationally intractable. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp, editors, Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 17020–17059. PMLR, 21–27 Jul 2024. URL https://proceedings.mlr.press/v235/gupta24a.html. 24
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models, 2020. URL https://arxiv.org/abs/2006.11239. Jerry Y Huang, Justin Lin, Sheel Shah, Kartik Nair, and Nicholas M Boffi. How to guide your flow: Few-step alignment via flow map reward guidance. arXiv preprint arXiv:2604.27147, 2026. Ioannis Karatzas and Steven E. Shreve. Brownian motion and stochastic calculus, volume 113 of Graduate Texts in Mathematics. Springer-Verlag, New York, second edition, 1991. ISBN 0-387-97655-8. doi: 10.1007/978-1-4612-0949-2. URL https://doi.org/10.1007/ 978-1-4612-0949-2. Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song. Denoising diffusion restoration models. In Advances in Neural Information Processing Systems, volume 35, pages 23593–23606, 2022. O. A. Ladyženskaja, V. A. Solonnikov, and N. N. Ural’ceva. Linear and Quasi-linear Equations of Parabolic Type, volume 23 of Translations of Mathematical Monographs. American Mathematical Society, Providence, RI, 1968. Holden Lee, Jianfeng Lu, and Yixin Tan. Convergence for score-based generative modeling with polynomial complexity, 2023. URL https://arxiv.org/abs/2206.06227. Ankur Moitra, Andrej Risteski, and Dhruv Rohatgi. Steering diffusion models with quadratic rewards: a fine-grained analysis, 2026. URL https://arxiv.org/abs/2602.16570. Badr MOUFAD, Yazid Janati, Lisa Bedin, Alain Oliviero Durmus, randal douc, Eric Moulines, and Jimmy Olsson. Variational diffusion posterior sampling with midpoint guidance. In The Thirteenth International Conference on Learning Representations, 2025. URL https://openreview.net/forum?id=6EUtjXAvmj. Advait Parulekar, Litu Rout, Karthikeyan Shanmugam, and Sanjay Shakkottai. Efficient approximate posterior sampling with annealed langevin monte carlo, 2025. URL https: //arxiv.org/abs/2508.07631. Aditya Ramesh, Mikhail Pavlov, Gabriel Goh, Scott Gray, Chelsea Voss, Alec Radford, Mark Chen, and Ilya Sutskever. Zero-shot text-to-image generation. In International Conference on Machine Learning, pages 8821–8831. PMLR, 2021. Yinuo Ren, Wenhao Gao, Lexing Ying, Grant M. Rotskoff, and Jiequn Han. Driftlite: Lightweight drift control for inference-time scaling of diffusion models. In The Fourteenth International Conference on Learning Representations, 2026. URL https://openreview. net/forum?id=l01eG3Qikl. H. Robbins. An empirical bayes approach to statistics. Proc. 3rd Berkeley Symp. Math. Statist. Probab., 1956, 1:157–163, 1956. URL https://cir.nii.ac.jp/crid/ 1572824500694511232. Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer. High-resolution image synthesis with latent diffusion models, 2022. URL https://arxiv. org/abs/2112.10752. Litu Rout, Yujia Chen, Abhishek Kumar, Constantine Caramanis, Sanjay Shakkottai, and Wen-Sheng Chu. Beyond first-order tweedie: Solving inverse problems using latent diffusion, 2023a. URL https://arxiv.org/abs/2312.00852. Litu Rout, Negin Raoof, Giannis Daras, Constantine Caramanis, Alexandros G Dimakis, and Sanjay Shakkottai. Solving inverse problems provably via posterior sampling with latent diffusion models. In Thirty-seventh Conference on Neural Information Processing Systems, 2023b. Litu Rout, Yujia Chen, Nataniel Ruiz, Abhishek Kumar, Constantine Caramanis, Sanjay Shakkottai, and Wen-Sheng Chu. RB-modulation: Training-free personalization using stochastic optimal control. In The Thirteenth International Conference on Learning Representations, 2025. URL https://openreview.net/forum?id=bnINPG5A32. 25
Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily Denton, Seyed Kamyar Seyed Ghasemipour, Burcu Karagol Ayan, S Sara Mahdavi, Rapha Gontijo Lopes, et al. Photorealistic text-to-image diffusion models with deep language understanding. arXiv preprint arXiv:2205.11487, 2022. Jascha Sohl-Dickstein, Eric A. Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics, 2015. URL https://arxiv. org/abs/1503.03585. Jiaming Song, Arash Vahdat, Morteza Mardani, and Jan Kautz. Pseudoinverse-guided diffusion models for inverse problems. In International Conference on Learning Representations (ICLR), May 2023. Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution, 2020. URL https://arxiv.org/abs/1907.05600. Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 2021. Yang Song, Liyue Shen, Lei Xing, and Stefano Ermon. Solving inverse problems in medical imaging with score-based generative models, 2022. URL https://arxiv.org/abs/2111. 08005. S. R. S. Varadhan. Asymptotic probabilities and differential equations. Communications on Pure and Applied Mathematics, 19(3):261–286, 1966. doi: 10.1002/cpa.3160190303. Santosh S. Vempala and Andre Wibisono. Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices, 2022. URL https://arxiv.org/abs/1903.08568. Luhuan Wu, Brian L. Trippe, Christian A. Naesseth, David M. Blei, and John P. Cunningham. Practical and asymptotically exact conditional sampling in diffusion models, 2024. URL https://arxiv.org/abs/2306.17775. Xingyu Xu and Yuejie Chi. Provably robust score-based diffusion posterior sampling for plug-and-play image reconstruction, 2024. URL https://arxiv.org/abs/2403.17042.
26