Paper deep dive
Scale-Consistent Posterior Dynamics for Diffusion Inverse Problems
Zhaoqiang Liu, Tongyao Pang, Ruibing Wang, Yang Zheng
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 8/18/2026, 5:45:34 AM
Summary
This paper proposes a scale-consistent posterior dynamics framework for solving diffusion inverse problems. It introduces an ideal one-parameter Stochastic Differential Equation (SDE) family where a stochasticity parameter controls exploration without altering posterior marginals. To make this tractable, the authors derive a scale-consistent surrogate in clean-image coordinates, organize it via log-SNR, and project diffusion uncertainty to form a noise-conditioned covariance path. A frozen-target Langevin corrector is interleaved with the transport to ensure the surrogate SDE tracks the moving targets. The method is discretized using a Lie-Trotter splitting and a variance-matched split-step IMEX predictor. Theoretical proofs cover marginal invariance, posterior convergence, and weak error bounds. Experiments on FFHQ and ImageNet demonstrate competitive reconstruction fidelity for super-resolution and deblurring.
Entities (10)
Relation Signals (9)
Zhaoqiang Liu → authored → Scale-Consistent Posterior Dynamics
confidence 99% · Title block: Zhaoqiang Liu ... Tongyao Pang ... Ruibing Wang ... Yang Zheng
Tongyao Pang → authored → Scale-Consistent Posterior Dynamics
confidence 99% · Title block: Zhaoqiang Liu ... Tongyao Pang
Scale-Consistent Posterior Dynamics → evaluatedon → ImageNet
confidence 99% · Abstract: Experiments on FFHQ and ImageNet
Scale-Consistent Posterior Dynamics → evaluatedon → FFHQ
confidence 99% · Abstract: Experiments on FFHQ and ImageNet
Scale-Consistent Posterior Dynamics → appliedto → Super-Resolution
confidence 95% · Abstract: competitive reconstruction fidelity for super-resolution
Scale-Consistent Posterior Dynamics → appliedto → Deblurring
confidence 95% · Abstract: competitive reconstruction fidelity for ... deblurring
Scale-Consistent Posterior Dynamics → solves → Diffusion Inverse Problems
confidence 95% · Abstract: Scale-Consistent Posterior Dynamics for Diffusion Inverse Problems
Scale-Consistent Posterior Dynamics → uses → IMEX Predictor
confidence 92% · Abstract: We discretize this model with ... a variance-matched split-step IMEX predictor
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Posterior sampling with a pretrained diffusion prior is governed by a conditional score whose intermediate likelihood component is generally intractable. We begin from an ideal one-parameter posterior SDE family in which a stochasticity parameter controls probability-flow transport and stochastic exploration without changing the posterior marginals. To obtain a tractable model, we express the likelihood in a rescaled clean-image coordinate and use log-SNR to organize the resulting posterior proxies. Projecting the diffusion uncertainty through the forward operator then yields a noise-conditioned covariance path whose targets approach the clean posterior. Because endpoint consistency of these targets does not ensure that a surrogate transport follows them, we interleave the transport with a frozen-target Langevin corrector, producing a continuous surrogate SDE. We discretize this model with an outer Lie--Trotter splitting and a variance-matched split-step IMEX predictor that treats the learned prior explicitly, the linear likelihood implicitly, and the stochastic innovation after the implicit solve. We prove marginal invariance of the ideal family, posterior convergence of the continuous surrogate under mixing and transport-defect conditions, and a first-order weak error bound for the discrete algorithm. Experiments on FFHQ and ImageNet with 100 score evaluations demonstrate competitive reconstruction fidelity for super-resolution and deblurring. A controlled 100-image ablation separates scale consistency from the finite-step effects of stochastic-increment placement, continuation, and corrector allocation. A separate noiseless box-inpainting study shows that large exploration reaches a performance plateau only when the matched innovation is injected after the stiff likelihood solve.
Tags
Links
- Source: https://arxiv.org/abs/2608.15144v1
- Canonical: https://arxiv.org/abs/2608.15144v1
Trouble viewing inline? Open PDF directly →
Full Text
93,967 characters extracted from source content.
Expand or collapse full text
Scale-Consistent Posterior Dynamics for Diffusion Inverse Problems Zhaoqiang Liu Thanks: University of Electronic Science and Technology of China, Chengdu, China. Tongyao Pang Email: typang@tsinghua.edu.cn Thanks: Corresponding author. Tsinghua University, Beijing, China (). Ruibing Wang11footnotemark: 1 Yang Zheng11footnotemark: 1 Abstract Posterior sampling with a pretrained diffusion prior is governed by a conditional score whose intermediate likelihood component is generally intractable. We begin from an ideal one-parameter posterior SDE family in which a stochasticity parameter controls probability-flow transport and stochastic exploration without changing the posterior marginals. To obtain a tractable model, we express the likelihood in a rescaled clean-image coordinate and use log-SNR to organize the resulting posterior proxies. Projecting the diffusion uncertainty through the forward operator then yields a noise-conditioned covariance path whose targets approach the clean posterior. Because endpoint consistency of these targets does not ensure that a surrogate transport follows them, we interleave the transport with a frozen-target Langevin corrector, producing a continuous surrogate SDE. We discretize this model with an outer Lie–Trotter splitting and a variance-matched split-step IMEX predictor that treats the learned prior explicitly, the linear likelihood implicitly, and the stochastic innovation after the implicit solve. We prove marginal invariance of the ideal family, posterior convergence of the continuous surrogate under mixing and transport-defect conditions, and a first-order weak error bound for the discrete algorithm. Experiments on FFHQ and ImageNet with 100 score evaluations demonstrate competitive reconstruction fidelity for super-resolution and deblurring. A controlled 100-image ablation separates scale consistency from the finite-step effects of stochastic-increment placement, continuation, and corrector allocation. A separate noiseless box-inpainting study shows that large exploration reaches a performance plateau only when the matched innovation is injected after the stiff likelihood solve. keywords diffusion inverse problems, posterior sampling, scale consistency, stochastic differential equations, Langevin dynamics, implicit-explicit methods †runningheads: Scale-Consistent Posterior Dynamics / Z. Liu, T. Pang, R. Wang, and Y. Zheng MSC 65J22, 65C30, 62C10, 68T07, 94A08 1 Introduction Many imaging tasks, including deblurring, super-resolution, tomography, and accelerated magnetic resonance imaging, can be modeled as linear inverse problems. We consider measurements of the form =+, y= A x+ n, (1) where ∈ℝm×d A ^m× d is known, ∈ℝd x ^d is the unknown clean image, and ∼(,σd2) n ( 0, _d^2 I). Because A is commonly ill-conditioned or rank deficient, reconstruction requires prior information. Classical choices include total variation and nonlocal self-similarity [31, 39, 7, 5, 10, 18, 15]. Learned restoration networks can instead infer task-specific regularities from paired or corrupted data [48, 26, 46, 4, 23, 24]. A pretrained generative model offers a different route: it can serve as a reusable prior without retraining for each sensing operator. In a Bayesian formulation, the target is the posterior p(∣)∝p(∣)p0(),p( x y) p( y x)p_0( x), where the Gaussian likelihood follows from (1) and p0p_0 is represented implicitly by the generative prior. Diffusion models are particularly suitable for this role because they learn the scores of a family of progressively smoothed data distributions [37, 17, 38]. In continuous time, the forward process is dt=(t,t)dt+g(t)dt,d x_t= f( x_t,t)\,dt+g(t)\,d W_t, (2) with 0∼p0 x_0 p_0 and a terminal law close to a tractable Gaussian. Its reverse-time SDE is dt=[(t,t)−g(t)2∇tlogpt(t)]dt+g(t)d¯t,t:T→0,d x_t= [ f( x_t,t)-g(t)^2 _ x_t p_t( x_t) ]dt+g(t)\,d W_t, t:T→ 0, (3) where the score ∇logpt∇ p_t is approximated by a neural network trained with denoising score matching [3, 38]. Posterior sampling also admits a stochastic optimal-control interpretation. Let ℙprP pr denote the reverse path law induced by (3). Conditioning on the measurement is the relative-entropy regularized control problem ℚ⋆=argminℚ≪ℙprKL(ℚ∥ℙpr)−ℚ[logp(∣0)],dℚ⋆dℙpr=p(∣0)p().Q_ y = *arg\,min_Q pr \KL(Q\,\|\,P pr)-E_Q\! [ p( y x_0) ] \, dQ_ y dP pr= p( y x_0)p( y). (4) Thus, the posterior path is the likelihood-tilted reference path. For drift-controlled laws with the same diffusion coefficient and initial law, Girsanov’s theorem identifies the path-space KL with the quadratic control cost. The associated Doob potential is ht()=ℙpr[p(∣0)∣t=]=pt(∣).h_t( x)=E_P pr\! [p( y x_0) x_t= x ]=p_t( y x). In the reverse-time convention of (3), the optimal h-transform adds the drift −g(t)2∇tloght(t).-g(t)^2 _ x_t h_t( x_t). Equivalently, it replaces the prior score with ∇tlogpt(t∣)=∇tlogpt(t)+∇tlogpt(∣t). _ x_t p_t( x_t y)= _ x_t p_t( x_t)+ _ x_t p_t( y x_t). The first term is supplied by the pretrained model, but the optimal control, or equivalently the intermediate likelihood score, is generally unavailable because pt(∣t)=∫p(∣0)p(0∣t)d0.p_t( y x_t)= p( y x_0)p( x_0 x_t)\,d x_0. (5) Approximating this integral at every noise level is the central difficulty in diffusion inverse problems. Existing methods use denoised point estimates, operator-specific decompositions, variational approximations, or auxiliary Monte Carlo procedures. These choices trade generality and computational cost against the accuracy and numerical stability of the likelihood guidance. An exact reference model clarifies what is lost when this likelihood is approximated. If the conditional score were available, the posterior reverse dynamics would belong to a one-parameter SDE family with identical time marginals. Its parameter ζ≥0ζ≥ 0 controls stochasticity: ζ=0ζ=0 is probability-flow transport, ζ=1ζ=1 is the usual reverse SDE, and larger values increase exploration without changing the ideal posterior path. This family provides the starting point for our method, but it cannot be implemented using an unconditional diffusion model alone. We obtain a tractable model in three steps. First, we express the measurement likelihood in the rescaled variable =t/αt z= x_t/ _t. This scale-consistent construction places the learned prior and data-fidelity terms in a common clean-image coordinate and prevents a diffusion-time-dependent mismatch between their scales. In the log-SNR variable, the resulting posterior proxies form a continuation toward the clean posterior. Second, because the rescaled state remains noisy away from the data endpoint, we augment the measurement covariance by the diffusion uncertainty projected through the forward operator. This covariance path softens unreliable measurement directions at high noise and recovers the clean likelihood at the endpoint. The associated targets converge in total variation to the Bayesian posterior; under additional regularity, their scores converge locally uniformly as well. Target consistency, however, is not sampler convergence: a transport driven by an approximate score need not have the prescribed continuation targets as its marginals. Our third step therefore interleaves the surrogate transport with a frozen-target Langevin corrector along the entire log-SNR path. The corrector equilibrates toward the current target while the transport advances that target, yielding a single continuous surrogate SDE. This distinction between a convergent target path and a dynamics that tracks it is central to both the construction and the analysis. For computation, we apply a Lie–Trotter predictor-corrector splitting to the continuous model. The transport predictor uses a variance-matched split-step IMEX update: a backward-Euler likelihood solve damps stiff deterministic transients, after which a calibrated stochastic innovation is added without being filtered by the likelihood resolvent. The frozen-target corrector uses a centered Crank–Nicolson update to reduce invariant-law bias. Both stages avoid an implicit neural-network solve. We prove posterior convergence of the continuous surrogate under log-Sobolev and transport-defect conditions, and then establish a first-order weak error bound for the complete discrete scheme. Together, these results separate the modeling error caused by the surrogate score from the numerical error caused by time discretization. Experiments on FFHQ and ImageNet use a fixed budget of 100 score evaluations for super-resolution and deblurring. The results show competitive reconstruction fidelity and expose the finite-budget trade-off between distortion and perceptual quality. Controlled paired ablations isolate scale consistency, stochastic-increment placement, noise-conditioned continuation, and the allocation of score evaluations between predictor and corrector steps. We treat noiseless box inpainting, where discretization and exploration interact more strongly, in a separate controlled study. Contributions The main contributions are as follows: • We formulate an exact one-parameter posterior SDE family and show that its stochasticity parameter changes the balance between exploration and exploitation without changing the posterior marginals. • We derive a scale-consistent surrogate in the clean-image coordinate and organize it as a noise-conditioned covariance path whose targets converge to the exact posterior. • We construct a continuously corrected surrogate SDE and a practical variance-matched split-step IMEX predictor whose post-solve innovation retains exploration in stiff measured directions. • We establish continuous posterior convergence and first-order weak accuracy, evaluate the framework on FFHQ and ImageNet, and separately quantify the contributions of its main modeling and numerical components, including exploration and splitting in noiseless box inpainting. 2 Related Work Diffusion models provide expressive image priors through the time-dependent score of the noised data distribution. This structure makes it possible to incorporate measurement information directly into the reverse diffusion dynamics, extending the plug-and-play principle from learned denoisers to stochastic generative processes [41, 30, 19, 36]. For an inverse problem with observation y, the posterior score at diffusion time t decomposes into the prior score ∇tlogpt(t) _ x_t p_t( x_t) and the likelihood score ∇tlogpt(∣t) _ x_t p_t( y x_t). The former is supplied by a pretrained score network, whereas the latter is generally intractable because it requires marginalizing the clean-image likelihood over p(0∣t)p( x_0 x_t). Consequently, most diffusion inverse solvers differ primarily in how they approximate or avoid this intermediate likelihood score [11]. A widely used strategy replaces the unknown clean image by a denoised estimate obtained from Tweedie’s formula and evaluates measurement consistency at that estimate. Diffusion Posterior Sampling (DPS) applies this construction directly, using the gradient of the measurement error evaluated at the denoised estimate to guide each reverse diffusion step [8]. Related methods refine the approximation by accounting for the local data-manifold geometry or the conditional covariance of the clean estimate, as in Manifold-Constrained Gradient (MCG) and Pseudoinverse-Guided Diffusion Models (Π ) [9, 34]. These approaches are flexible and require only differentiation through the forward operator, but the resulting guidance depends on the accuracy of the denoised estimate. At high noise levels, errors in this estimate can be amplified by the forward operator, while an improperly scaled likelihood gradient can dominate or vanish relative to the learned prior score. For linear inverse problems, additional structure in the sensing operator can be used to enforce measurement consistency more directly. SNIPS and Denoising Diffusion Restoration Models (DDRM) decompose the dynamics through the singular vectors of the sensing operator and treat its observed and unobserved components separately [22, 21]. Denoising Diffusion Null-Space Models (DDNM) similarly preserve the measurement-determined range-space component while using the diffusion prior to generate the complementary null-space component [42]. These approaches are effective when the required pseudoinverse or singular-value decomposition is available, but their implementation is tied to the structure of the forward operator. Denoising Diffusion Models for Plug-and-Play Image Restoration (DiffPIR) instead embeds the diffusion denoiser in a half-quadratic splitting scheme and alternates prior denoising with data-consistency updates [50]. Such splitting schemes are easy to interpret, although their stability and accuracy can remain sensitive to the data-consistency weight and to the number of inner corrections. Other methods avoid a direct pointwise approximation of the likelihood score. RED-Diff adopts a variational formulation in which a tractable distribution is optimized against the diffusion posterior using a score-matching surrogate [27]. Midpoint-Guided Posterior Sampling (MGPS) decomposes the reverse transition and evaluates the intractable guidance at an intermediate state, providing a trade-off between the complexity of the guidance approximation and that of the prior transition [28]. Monte Carlo approaches instead introduce auxiliary variables, filtering recursions, or particles to obtain more faithful posterior samples [44, 14, 6, 43, 40]. Decoupled Annealing Posterior Sampling (DAPS) separates noise annealing from the reverse diffusion trajectory, which permits larger nonlocal updates at the cost of additional inner sampling steps [47]. Consistency-based inverse solvers provide another route to reducing the number of score evaluations [45]. The continuation viewpoint predates diffusion models. Numerical continuation embeds a difficult problem in a parameterized family and tracks its solutions from an easier starting point, often using predictor and corrector steps [2]. Bayesian annealing applies the same idea to probability laws by bridging a tractable reference and the posterior with intermediate distributions [29, 12]. Score-based models already provide a natural path of progressively less smoothed distributions through the Gaussian noise scale [35], and DAPS uses a related annealing principle for diffusion inverse problems [47]. In our construction, this viewpoint arises after scale alignment: log-SNR indexes clean-coordinate posterior proxies, and projected diffusion noise defines an auxiliary covariance path for the likelihood. The continuation regularizes intermediate targets, while scale alignment and the implicit discretization address the separate problems of coordinate mismatch and numerical stiffness. Against this background, we use an exact one-parameter posterior SDE family as the reference model and explicitly separate three approximation layers. The scale-consistent likelihood first replaces the intractable conditional score; the covariance path then moderates this likelihood away from the data endpoint; and a frozen-target Langevin corrector addresses the fact that endpoint consistency of the targets does not imply convergence of the surrogate sampler law. Predictor-corrector sampling is also standard in unconditional score-based generation [38], but here the two components have distinct roles: the predictor transports samples along a changing posterior path, whereas the corrector equilibrates toward a fixed target at each level. The resulting discretization evaluates the learned prior explicitly, treats the linear likelihood implicitly, and separates the stochastic update from the implicit resolvent, retaining stable linear solves without requiring an implicit network evaluation. 3 Method We develop the sampler by moving from an ideal posterior model to a tractable numerical method. The ideal model first identifies the posterior path and separates the choice of stochastic exploration from the target distribution. Because its conditional score is unavailable, we then construct a scale-aligned continuation surrogate and add a Langevin correction that tracks the moving targets. A split-step IMEX discretization turns the resulting continuous model into the practical algorithm, whose continuous convergence and finite-step error are analyzed at the end of the section. 3.1 Ideal Posterior Dynamics We begin with the posterior SDE that would be available if the exact conditional score were known. Besides providing a reference model for the approximations that follow, this formulation introduces a stochasticity parameter ζ that controls the balance between deterministic transport and stochastic exploration without changing the exact posterior marginals. We identify the unknown clean signal with the diffusion endpoint 0∈ℝd x_0 ^d and observe noisy linear measurements =0+, y= A x_0+ n, (6) where ∈ℝm×d A ^m× d is known and ∼(,σd2) n ( 0, _d^2 I). For the convergence analysis we assume σd>0 _d>0; singular noiseless constraints are not covered by the total-variation theorem below. The prior p0p_0 is represented by the variance-preserving SDE obtained from (2) with (,t)=−12β(t) f( x,t)=- 12β(t) x and g(t)=β(t)g(t)= β(t): dt=−12β(t)tdt+β(t)dt,0∼p0,d x_t=- 12β(t) x_t\,dt+ β(t)\,d W_t, x_0 p_0, (7) where β(t)>0β(t)>0 is a prescribed noise schedule. Its marginal admits the closed form t=αt0+σtϵ,ϵ∼(,), x_t= _t x_0+ _t ε, ε ( 0, I), (8) with αt:=exp(−∫0tβ(τ)dτ),σt2:=1−αt. _t:= \! (-\! _0^tβ(τ)\,dτ ), _t^2:=1- _t. (9) For sufficiently large T, pTp_T is close to the standard Gaussian density. The corresponding reverse-time prior dynamics are dt=[−12β(t)t−β(t)∇tlogpt(t)]dt+β(t)d~t,t:T→0,d x_t= [- 12β(t) x_t-β(t) _ x_t p_t( x_t) ]dt+ β(t)\,d W_t, t:T→ 0, (10) where the score is approximated in practice by a pretrained network. Conditioning the initial law on the observation changes the target to p(0∣)=p0(0)p(∣0)p().p( x_0 y)= p_0( x_0)p( y x_0)p( y). Let pt(⋅∣)p_t(· y) be the marginal obtained by initializing (7) with 0∼p0(⋅∣) x_0 p_0(· y). If its score were available, the stochastic posterior reverse SDE would be dt=[−12β(t)t−β(t)∇tlogpt(t∣)]dt+β(t)d¯t,t:T→0.d x_t= [- 12β(t) x_t-β(t) _ x_t p_t( x_t y) ]dt+ β(t)\,d W_t, t:T→ 0. (11) This stochastic process is one member of a larger family with the same posterior marginals. Proposition 1 (Marginal invariance of the posterior SDE family). Let p0(⋅∣)p_0(· y) be a Borel probability measure. Assume that for every ε∈(0,T) ∈(0,T) the reverse-time SDE below and its Fokker–Planck equation are well posed on [ε,T][ ,T]. For any ζ≥0ζ≥ 0, consider dt=(−12β(t)t−1+ζ2β(t)∇tlogpt(t∣))dt+ζβ(t)d¯t,t:T→0,d x_t= (- 12β(t)\, x_t- 1+ζ2\,β(t)\, _ x_t p_t( x_t y) )\,dt+ ζ\,β(t)\,d W_t, t:T→ 0, (12) initialized with T∼pT(⋅∣) x_T p_T(· y). Then the marginal density of t x_t is pt(⋅∣)p_t(· y) for every t∈(0,T]t∈(0,T]. Moreover, pt(⋅∣)p_t(· y) converges weakly to p0(⋅∣)p_0(· y) as t↓0t 0. The choice ζ=1ζ=1 gives the stochastic posterior reverse SDE (11), whereas ζ=0ζ=0 gives the posterior probability-flow ODE. Proof. See Appendix A. Within this exact family, decreasing ζ reduces the injected noise and approaches deterministic probability-flow transport, thereby emphasizing exploitation of the score field. Increasing ζ strengthens the Brownian component and promotes exploration. Proposition 1 shows that this pathwise choice does not alter the ideal marginal law. For sufficiently large T, pT(⋅∣)p_T(· y) is close to the standard Gaussian density, motivating the practical initialization T∼(,) x_T ( 0, I). Near the data endpoint, we stop at a sufficiently small ε>0 >0 and use ε x_ as an approximate posterior sample. The family (12) is an ideal reference because its conditional score is not available from the unconditional diffusion model. The next subsection replaces this score by a tractable continuation and then corrects the resulting mismatch at the level of the sampler law. 3.2 Continuous Surrogate The conditional score in (12) decomposes into the unconditional prior score and an intermediate likelihood score. The latter is intractable because pt(∣t)p_t( y x_t) requires averaging the clean likelihood over p(0∣t)p( x_0 x_t), as shown in (5). We first replace this integral by a scale-aligned clean-coordinate proxy. In the rescaled log-SNR variable, the proxy naturally becomes a continuation toward the clean posterior, which in turn suggests a noise-dependent continuation of the likelihood covariance. The resulting score defines a tractable transport SDE, but its law need not follow the moving continuation targets. A frozen-target Langevin component repairs this remaining distributional mismatch and completes the continuous surrogate model. Scale alignment A direct use of the clean likelihood at the noisy state would compare t A x_t with y, even though the VP marginal (8) scales the signal component by αt _t. To remove this coordinate mismatch, we first retain the deterministic scaling in (8) and use the clean-coordinate proxy 0≈t/αt x_0≈ x_t/ _t. This gives the scale-consistent surrogate likelihood Ltsc(t):=exp(−12σd2‖tαt−‖22),L_t sc( x_t):= \! (- 12 _d^2 \| A x_t _t- y \|_2^2 ), (13) whose score is ∇tlogLtsc(t)=−1αtσd2⊤(tαt−). _ x_t L_t sc( x_t)=- 1 _t _d^2 A ( A x_t _t- y ). (14) Equivalently, define the positive residual term tsc(t,):=1αtσd2⊤(tαt−)=−∇tlogLtsc(t). d_t sc( x_t; y):= 1 _t _d^2 A ( A x_t _t- y )=- _ x_t L_t sc( x_t). (15) Substitution into the one-parameter reverse family yields the baseline scale-consistent surrogate dt=[−12β(t)t−1+ζ2β(t)(∇tlogpt(t)−tsc(t,))]dt+ζβ(t)d¯t,t:T→0.d x_t= [- 12β(t) x_t- 1+ζ2β(t) ( _ x_t p_t( x_t)- d_t sc( x_t; y) ) ]dt+ ζβ(t)\,d W_t, t:T→ 0. (16) The factor 1/αt1/ _t in (14) is the chain-rule factor required when the clean likelihood is evaluated through the rescaled state. Thus, the learned prior and the data-fidelity term are both interpreted in the same clean-image coordinate. This scale alignment is the mechanism evaluated by the scale-consistency ablations in Section 4. To expose this structure more clearly, define the rescaled log-SNR coordinate λ:=tαt=0+τλϵ,τλ:=e−λ=1−αtαt, z_λ:= x_t _t= x_0+ _λ ε, _λ:=e^-λ= 1- _t _t, (17) where λ=12logαt1−αt.λ= 12 _t1- _t. Let πλ _λ denote the density of λ z_λ. In this coordinate, the scale-consistent target and its score take the simple forms π~λsc(∣) π_λ sc( z y) ∝πλ()exp(−12σd2‖−‖22), _λ( z) \! (- 12 _d^2\| A z- y\|_2^2 ), (18) ~λsc(,) s_λ sc( z; y) =∇logπλ()−1σd2⊤(−). = _ z _λ( z)- 1 _d^2 A ( A z- y). Thus, the scale-consistent approximation becomes a family of clean-coordinate posterior proxies indexed by λ. Since πλ=p0∗(,e−2λ) _λ=p_0*N( 0,e^-2λ I), increasing λ progressively removes the Gaussian smoothing and exposes the clean posterior. In this sense, the change of variables turns the scale-consistent posterior proxy into a continuation path. It still treats the perturbed variable λ=0+τλϵ z_λ= x_0+ _λ ε as clean when evaluating the likelihood, which motivates a second continuation within the likelihood itself. Continuation path Away from the data endpoint, the residual −λ y- A z_λ contains both measurement noise and the projected diffusion perturbation. We account for both contributions through the covariance λ:=σd2+sλ⊤,sλ≥0,sλ⟶0. _λ:= _d^2 I+s_λ A A , s_λ≥ 0, s_λ 0. (19) The log-SNR λ is the continuation parameter. The nominal choice sλ=τλ2s_λ= _λ^2 is the projected diffusion variance; a smooth schedule with the same zero-noise limit can be used to control how rapidly the data term is exposed. The additional covariance softens uncertain measurement directions at high noise and vanishes as λ increases, thereby recovering the target likelihood. This defines a covariance continuation of probability densities, rather than a deterministic solution path or a scalar power-posterior path. The corresponding surrogate likelihood is Lλ():=exp(−12‖−‖λ−12).L_λ( z):= \! (- 12\| A z- y\|_ _λ^-1^2 ). (20) This likelihood is a continuation surrogate rather than the exact conditional density of y given λ z_λ. It preserves the clean-coordinate scaling in (18), but replaces σd2 _d^2 I by λ _λ to weaken unreliable measurement directions at high diffusion noise. Since λ→σd2 _λ→ _d^2 I as λ→∞λ→∞, it recovers the scale-consistent likelihood and score at the data endpoint. The nominal schedule is also aligned with the prior transport. Indeed, if PqP_q denotes convolution with (,q)N( 0,q I) and q=τλ2q= _λ^2, then πλ=Pqp0 _λ=P_qp_0 and, up to a scalar independent of z, Lλ=PqL∞L_λ=P_qL_∞. Both factors therefore satisfy ∂λf=−e−2λΔf _λf=-e^-2λ f. Their normalized product is not generally an exact transport marginal because the heat semigroup does not preserve products; the remaining interaction is quantified by the transport defect in (37). We define the noise-conditioned continuation target π¯λ(∣):=1Zλπλ()Lλ(),Zλ:=∫πλ()Lλ(). π_λ( z y):= 1Z_λ _λ( z)L_λ( z), Z_λ:= _λ( z)L_λ( z)\,d z. (21) Its score is ¯λ(,)=∇logπλ()−λ+λ,λ:=⊤λ−1,λ:=⊤λ−1. s_λ( z; y)= _ z _λ( z)- H_λ z+ b_λ, H_λ:= A _λ^-1 A, b_λ:= A _λ^-1 y. (22) Because λ→σd2 _λ→ _d^2 I and Gaussian smoothing vanishes as λ→∞λ→∞, this target converges to p(⋅∣)p(· y) in total variation. If, in addition, p0p_0 has a strictly positive C1C^1 density and its Gaussian mollifications satisfy ∇logπλ→∇logp0∇ _λ→∇ p_0 locally uniformly, then the target score converges locally uniformly to the clean-posterior score; see Appendix B. Thus, target consistency is already built into the continuation path before a sampling dynamics is chosen. To turn this target path into a computable transport SDE, we evaluate its prior component in the noise-prediction parameterization: ∇logπλ()≈θ,λpr():=−eλϵ^θ(αλ,λ). _ z _λ( z)≈ s_θ,λ pr( z):=-e^λ ε_θ( _λ z,λ). (23) Substituting (22) into the one-parameter construction gives the continuation-driven surrogate dynamics dλ=(1+ζ)e−2λ¯λ(λ,)dλ+2ζe−λdλ,ζ≥0.d z_λ=(1+ζ)e^-2λ s_λ( z_λ; y)\,dλ+ 2ζ\,e^-λ\,d W_λ, ζ≥ 0. (24) Unlike the exact family in Proposition 1, the marginal law of (24) need not track π¯λ π_λ. Endpoint consistency of the targets alone therefore does not guarantee convergence of the sampler law. This distinction leads to the final component of the continuous model. Corrected transport We supplement the transport SDE (24) with a frozen-target corrector. For fixed λ, consider the isotropic Langevin diffusion in an auxiliary time r, dr=¯λ(r,)dr+2dr.d u_r= s_λ( u_r; y)\,dr+ 2\,d B_r. (25) For an exact prior score, π¯λ π_λ is invariant for (25). The distinction between (24) and (25) is structural. Equation (24) evolves in the physical continuation clock λ and transports samples while the target π¯λ π_λ itself changes; consequently, it does not generally preserve that target. By contrast, (25) freezes λ and evolves in a separate MCMC clock r. Its matched drift and noise satisfy the fluctuation–dissipation relation, so π¯λ π_λ is invariant and the corrector only equilibrates toward the current target. Thus, the corrector is not a second copy of the predictor even though both drifts contain the same continuation score. To distribute this correction along the entire continuation path, we combine the two processes at the generator level. Let γλ≥0 _λ≥ 0 denote the corrector intensity, meaning that an auxiliary corrector time γλdλ _λdλ is run during an infinitesimal continuation interval dλdλ. Combining the generators of (24) and (25) gives the continuously interleaved SDE dλ= d z_λ= [(1+ζ)e−2λ+γλ]¯λ(λ,)dλ [(1+ζ)e^-2λ+ _λ ] s_λ( z_λ; y)\,dλ (26) +2ζe−λdλ+2γλdλ, + 2ζ\,e^-λ\,d W_λ+ 2 _λ\,d B_λ, where W and B are independent Brownian motions. This is the continuous surrogate SDE model underlying our predictor-corrector sampler. It retains the exploration parameter inherited from the ideal family, uses the continuation score as a tractable transport direction, and allocates an auxiliary Langevin time γλdλ _λdλ to control target tracking. 3.3 Split-Step Discretization We retain the continuous surrogate SDE. At the outer level, a Lie–Trotter splitting advances the continuation predictor and then runs a corrector at the new noise level. Within the predictor, the learned prior is explicit and the linear likelihood is implicit. A direct backward-Euler discretization would place the Brownian increment inside the same likelihood resolvent as the drift. Although this is stable, stability alone need not preserve the finite-step stochastic law. Implicit SDE solvers can bias nondegenerate invariant distributions when stiff scales are unresolved [25]; even for an Ornstein–Uhlenbeck process, different implicit Euler choices have different law-preservation properties [33]. Split-step backward Euler separates the implicit drift solve from the stochastic update [16], and the ordering of stochastic subflows can materially affect distributional accuracy [1]. These observations lead to the variance-matched split-step IMEX predictor below. We use an increasing log-SNR grid λ0<λ1<⋯<λK,τk:=e−λk, _0< _1<·s< _K, _k:=e^- _k, which follows the reverse generative direction from high noise to low noise. For k=0,…,K−1k=0,…,K-1, define hk h_k :=λk+1−λk, := _k+1- _k, Δ1,k _1,k :=τk−τk+1, := _k- _k+1, (27) Δ2,k _2,k :=τk2−τk+12, := _k^2- _k+1^2, ωk _k :=1+ζ2Δ2,k. := 1+ζ2 _2,k. Let ℛkpredR_k pred denote the exact transition kernel of (24) from λk _k to λk+1 _k+1, and let λrC_λ^r denote the corrector semigroup of (25) run for auxiliary time r. Define rk+1:=∫λkλk+1γss.r_k+1:= _ _k _k+1 _s\,ds. (28) Freezing the corrector target at the right endpoint gives μk+1≈(μkℛkpred)λk+1rk+1. _k+1≈ ( _kR_k pred )C_ _k+1^r_k+1. (29) The first substep transports the law to the next continuation level, and the second substep equilibrates it toward the frozen target at that level. Under standard smoothness conditions, this ordering is a first-order splitting of (26). Variance-matched split-step predictor Write ϵ^k:=ϵ^θ(αλkk,λk),^0,k:=k−τkϵ^k,k+1:=λk+1,k+1:=λk+1, ε_k:= ε_θ( _ _k z_k, _k), x_0,k:= z_k- _k ε_k, H_k+1:= H_ _k+1, b_k+1:= b_ _k+1, where ^0,k=k+τk2kpr(k) x_0,k= z_k+ _k^2 s_k pr( z_k) is the denoised clean estimate. The predictor SDE has the integral form k+1= z_k+1= k−(1+ζ)∫λkλk+1e−λϵ^θ(αλλ,λ)λ z_k-(1+ζ) _ _k _k+1e^-λ ε_θ( _λ z_λ,λ)\,dλ −(1+ζ)∫λkλk+1e−2λ(λ−λ)dλ+2ζ∫λkλk+1e−λdλ. -(1+ζ) _ _k _k+1e^-2λ( H_λ z_λ- b_λ)\,dλ+ 2ζ _ _k _k+1e^-λ\,d W_λ. For the deterministic stage, we freeze the denoised estimate ^0,k x_0,k over the step and the quadratic likelihood at the right endpoint. Define ak:=e−(1+ζ)hk,qkvm:=τk+12(1−e−2ζhk),k:=(+ωkk+1)−1.a_k:=e^-(1+ζ)h_k, q_k vm:= _k+1^2 (1-e^-2ζ h_k ), R_k:= ( I+ _k H_k+1 )^-1. (30) The predictor consists of one deterministic IMEX stage followed by one variance-matched stochastic stage: (+ωkk+1)k+1det ( I+ _k H_k+1 ) z_k+1 det =^0,k+ak(k−^0,k)+ωkk+1, = x_0,k+a_k ( z_k- x_0,k )+ _k b_k+1, k+1pred z_k+1 pred =k+1det+qkvmk,k∼(,). = z_k+1 det+ q_k vm\, ξ_k, ξ_k ( 0, I). (31) The right-hand side of the first line is the exact mean of the frozen-denoiser prior flow, followed by a backward-Euler likelihood solve. The second line is applied after that solve, so the likelihood resolvent does not attenuate the new stochastic directions. Its variance is the exact conditional variance of the same frozen-denoiser prior flow. The term “variance matched” refers to this prior conditional variance; it does not assert exact integration of the full nonlinear posterior SDE over a finite step. The next proposition isolates the finite-step effect and the asymptotic consistency of this choice. Proposition 2 (Covariance retention and first-order consistency). Suppose k+1 H_k+1 is symmetric positive semidefinite and let ~k+1:=^0,k+ak(k−^0,k)+ωkk+1. z_k+1:= x_0,k+a_k( z_k- x_0,k)+ _k b_k+1. Conditional on k z_k, the split-step update in (31) has covariance qkvmq_k vm I. By contrast, inserting the same innovation inside the backward-Euler solve, k+1dir=k(~k+1+qkvmk), z_k+1 dir= R_k ( z_k+1+ q_k vm ξ_k ), gives covariance qkvmk2q_k vm R_k^2. Hence a likelihood eigenmode of curvature ν retains only the fraction (1+ωkν)−2(1+ _kν)^-2 of the injected variance under direct backward Euler. Moreover, for every fixed ζ≥0ζ≥ 0, qkvm=ζΔ2,k+O(hk2)as hk→0.q_k vm=ζ _2,k+O(h_k^2) h_k→ 0. (32) Thus the variance matching changes the finite-step covariance but not the first-order weak limit of the original predictor SDE. Proof. The two covariance formulas follow directly from the affine Gaussian updates and the symmetry of k R_k. Since τk+1=τke−hk _k+1= _ke^-h_k, Taylor expansion gives qkvm=2ζτk2hk+O(hk2)q_k vm=2ζ _k^2h_k+O(h_k^2) and ζΔ2,k=2ζτk2hk+O(hk2)ζ _2,k=2ζ _k^2h_k+O(h_k^2). This distinction between deterministic stability and distributional accuracy is already visible in the normalized Ornstein–Uhlenbeck test of [25]. If c denotes the time step divided by the fast relaxation scale, the stationary variance produced by their θ-method is [1−1−2θ2c]−1. [1- 1-2θ2c ]^-1. Crank–Nicolson (θ=1/2θ=1/2) therefore retains the unit target variance, whereas backward Euler (θ=1θ=1) gives 2/(2+c)2/(2+c), which tends to zero in the unresolved stiff limit. Proposition 2 is a one-step, mode-wise counterpart of the same damping mechanism: the implicit resolvent is useful for the deterministic likelihood drift, but placing the new innovation inside that resolvent can suppress the nondegenerate random directions. The analogy concerns covariance damping; our predictor remains a nonstationary, learned-score discretization rather than a fixed OU sampler. The same reference shows that the trapezoidal rule does not generically retain the correct invariant distribution for nonlinear stiff SDEs. Accordingly, we use the OU calculation only to motivate covariance preservation for the linear–Gaussian likelihood subflow, not to claim exact invariant-law preservation for the full learned dynamics. Frozen-target corrector Let k+1pr:=θ,λk+1pr s_k+1 pr:= s_θ, _k+1 pr, (0)=k+1pred u^(0)= z_k+1 pred, η=ηk+1η= _k+1, and let (j)∼(,) ξ^(j) ( 0, I) be independent. One corrector step is (+η2k+1)(j+1) ( I+ η2 H_k+1 ) u^(j+1) =(−η2k+1)(j)+η[k+1pr((j))+k+1]+2η(j). = ( I- η2 H_k+1 ) u^(j)+η [ s_k+1 pr( u^(j))+ b_k+1 ]+ 2η\, ξ^(j). (33) The predictor and corrector share the same IMEX decomposition but serve different purposes. The backward-Euler likelihood stage in the predictor is L-stable and damps stiff deterministic transients, while the post-solve innovation prevents that damping from also collapsing exploration. The Crank–Nicolson corrector is A-stable rather than L-stable, but its centered drift and matched noise preserve the covariance of the likelihood-only Gaussian dynamics at a fixed level [33]. This invariant-law property is more relevant to correction than stronger deterministic damping. The remaining corrector bias comes from the explicit learned score, the finite corrector step size, and score error. When ζ=0ζ=0, the predictor noise vanishes and the transport update is deterministic. If the corrector at λk+1 _k+1 uses step sizes ηk+1,jj=0Mk+1−1\ _k+1,j\_j=0^M_k+1-1, its total numerical auxiliary time is rk+1num:=∑j=0Mk+1−1ηk+1,j.r_k+1 num:= _j=0^M_k+1-1 _k+1,j. (34) For the discrete scheme to approximate a fixed finite-intensity process (26) as the log-SNR grid is refined, the budgets should satisfy rk+1num=∫λkλk+1γss+o(hk),maxjηk+1,j⟶0.r_k+1 num= _ _k _k+1 _s\,ds+o(h_k), _j _k+1,j 0. (35) For the constant step size used in Algorithm 1, this reduces to Mk+1ηk+1=γλk+1hk+o(hk).M_k+1 _k+1= _ _k+1h_k+o(h_k). (36) The practical algorithm may deliberately use a larger corrector budget. Its effective intensity is γk+1eff:=Mk+1ηk+1/hk _k+1 eff:=M_k+1 _k+1/h_k. If Mk+1ηk+1M_k+1 _k+1 is held fixed while hk→0h_k→ 0, then γk+1eff→∞ _k+1 eff→∞. This regime is an increasingly strong, or over-corrected, splitting that approaches equilibration at each level; it is not a consistent discretization of a fixed finite γλ _λ. Algorithm 1 summarizes the split-step IMEX realization of (29). Each predictor requires one score evaluation and each corrector step requires one additional score evaluation. At a finite trained endpoint, the implementation uses one final score evaluation to form the denoised estimate ^0=K−τKϵ^K x_0= z_K- _K ε_K; this call is included in the reported NFE budget. Algorithm 1 Noise-Conditioned Split-Step IMEX Predictor-Corrector Sampler 0: Measurement y, operator A, noise level σd _d, score network ϵ^θ ε_θ, increasing log-SNR grid λkk=0K\ _k\_k=0^K, parameter ζ≥0ζ≥ 0, corrector steps Mk\M_k\, step sizes ηk\ _k\ 0: Approximate posterior sample ^0 x_0 1: Compute αλk _ _k, hkh_k, and Δ2,k _2,k from (27) 2: Sample λ0∼(,) x_ _0 ( 0, I) and set 0←λ0/αλ0 z_0← x_ _0/ _ _0 3: for k=0k=0 to K−1K-1 do 4: Form λk+1 _ _k+1, λk+1 H_ _k+1, and λk+1 b_ _k+1 from (19) and (22) 5: Compute ^0,k x_0,k, aka_k, qkvmq_k vm, and the predictor k+1pred z_k+1 pred using (31) 6: (0)←k+1pred u^(0)← z_k+1 pred 7: for j=0j=0 to Mk+1−1M_k+1-1 do 8: Compute (j+1) u^(j+1) from (33) with λ=λk+1λ= _k+1 and η=ηk+1η= _k+1 9: end for 10: k+1←(Mk+1) z_k+1← u^(M_k+1) 11: end for 12: Kalg←K z_K alg← z_K 13: Compute ϵ^K ε_K at the lowest finite trained level and set ^0←Kalg−τKϵ^K x_0← z_K alg- _K ε_K 14: return ^0 x_0 3.4 Convergence and Weak Error Subsection 3.2 establishes that the continuation targets approach the desired posterior, while also showing why this endpoint property does not by itself control the sampler law. We now close that gap at the continuous level and then quantify the additional error introduced by the split-step IMEX discretization. The first result proves posterior convergence of the continuous surrogate SDE; the second gives the finite-interval weak error of its discrete realization. To quantify the fact that the moving continuation targets are not exact marginals of the predictor, define their transport defect by rλ() r_λ( z) :=∂λlogπ¯λ()+e−2λΔπ¯λ()π¯λ() := _λ π_λ( z)+e^-2λ π_λ( z) π_λ( z) (37) =∂λlogπ¯λ()+e−2λ[Δlogπ¯λ()+‖¯λ(,)‖22]. = _λ π_λ( z)+e^-2λ [ π_λ( z)+\| s_λ( z; y)\|_2^2 ]. The defect vanishes precisely when the target path satisfies ∂λπ¯λ=−e−2λΔπ¯λ _λ π_λ=-e^-2λ π_λ, the marginal equation associated with the predictor. For the nominal schedule sλ=τλ2s_λ= _λ^2, the product rule gives the explicit identity rλ=2e−2λ(λpr⋅λlik−π¯λ[λpr⋅λlik]),λpr:=∇logπλ,λlik:=−λ+λ.r_λ=2e^-2λ ( s_λ pr\!· s_λ lik-E_ π_λ[ s_λ pr\!· s_λ lik] ), s_λ pr:=∇ _λ, s_λ lik:=- H_λ z+ b_λ. (38) Thus the nominal prior and likelihood continuations are separately compatible with the predictor, while their score interaction is the only residual transport mismatch. In particular, this mismatch vanishes near the endpoint if the centered interaction remains uniformly integrable. We assume that, for every law ν in the evolution with finite relative entropy, −∫rλdν≤cλKL(ν∥π¯λ)+ελ,cλ,ελ≥0.- r_λ\,dν≤ c_λKL(ν\,\|\, π_λ)+ _λ, c_λ, _λ≥ 0. (39) This is a structural error of the continuous surrogate target, not a time-discretization error. For example, the entropy variational inequality shows that (39) holds, for any θ>0θ>0 for which the exponential moment is finite, with cλ=θ−1,ελ=θ−1log∫e−θrλdπ¯λ.c_λ=θ^-1, _λ=θ^-1 e^-θ r_λ\,d π_λ. (40) Hence, if this exponential moment tends to one for some fixed θ, then ελ→0 _λ→ 0; the corrector intensity can subsequently be chosen so that 2γλρλ2 _λ _λ dominates θ−1θ^-1. Theorem 3 (Posterior convergence of the continuously corrected dynamics). Assume that the score in (26) is exact and that its coefficients and the positive smooth target family are sufficiently regular for the SDE and the entropy identities below to be well posed. Suppose that (39) holds and that, for every λ, π¯λ π_λ satisfies the log-Sobolev inequality Entπ¯λ(f2)≤2ρλ∫‖∇f‖22dπ¯λ,ρλ>0.Ent_ π_λ(f^2)≤ 2 _λ \|∇ f\|_2^2\,d π_λ, _λ>0. (41) Let μλ _λ be the law of (26), assume H(λ0):=KL(μλ0∥π¯λ0)<∞H( _0):=KL( _ _0\,\|\, π_ _0)<∞, and define kλ:=2γλρλ−cλ.k_λ:=2 _λ _λ-c_λ. Then, for every λ≥λ0λ≥ _0, H(λ)≤ H(λ)≤ e−∫λ0λkuduH(λ0) e^- _ _0^λk_u\,duH( _0) (42) +∫λ0λe−∫sλkuduεsds. + _ _0^λe^- _s^λk_u\,du _s\,ds. Consequently, if ∫λ0∞kudu=∞,∫λ0λe−∫sλkuduεsds⟶0, _ _0^∞k_u\,du=∞, _ _0^λe^- _s^λk_u\,du _s\,ds 0, (43) then ∥μλ−p(⋅∣)∥TV⟶0as λ→∞.\| _λ-p(· y)\|_ TV 0 λ→∞. (44) A simple sufficient condition is that, eventually, kλ≥k∗>0k_λ≥ k_*>0 and ελ→0 _λ→ 0. Proof. See Appendix C. Theorem 3 makes the role of the corrector budget explicit: γλ _λ should be increased when the mixing constant ρλ _λ deteriorates or when the transport defect grows. This mechanism distributes mixing along the continuation path, although a local Langevin corrector cannot in general remove a genuinely multimodal mixing bottleneck. The theorem concerns the exact continuous process (26); score and time-discretization errors require a separate analysis. The continuation path enters this result through its endpoint, transport, and mixing properties. The consistency statement in Subsection 3.2 ensures that the targets approach the posterior, rλr_λ measures their departure from the predictor’s natural marginal equation, and ρλ _λ controls the rate at which the corrector reduces this departure. A useful continuation should therefore keep the transport defect small without causing an unnecessary deterioration of the mixing constant. If the target path satisfies the predictor marginal equation exactly, then rλ=0r_λ=0 and the defect terms cλc_λ and ελ _λ can be removed. The continuation changes the constants in the numerical analysis below but not its first-order rate, provided that the coefficients remain sufficiently smooth. We now analyze the numerical error of the complete Algorithm 1 on a fixed finite interval. This is a different limit from Theorem 3: the terminal log-SNR is fixed while the splitting grid and the corrector steps are refined. The theorem concerns the transported state Kalg z_K alg before the finite-level terminal denoiser. The latter is the reported point estimate in our experiments. It approaches the identity map as τK→0 _K→ 0 under the same score regularity, so it does not change the asymptotic posterior limit. Theorem 4 (First-order weak error of the complete split-step scheme). Fix Λ>λ0 > _0 and ζ≥0ζ≥ 0, and let μΛ _ be the law at Λ of the continuously interleaved SDE (26). Let Kalg z_K alg be generated by Algorithm 1 on a grid with λK=Λ _K= , using the exact continuation score in both substeps. Define h h :=max0≤k<Khk, := _0≤ k<Kh_k, η∗ _* :=max1≤k≤Kηk, := _1≤ k≤ K _k, (45) δbud _ bud :=∑k=0K−1|Mk+1ηk+1−rk+1|, := _k=0^K-1 |M_k+1 _k+1-r_k+1 |, rk+1 r_k+1 :=∫λkλk+1γsds. := _ _k _k+1 _s\,ds. Assume that on [λ0,Λ][ _0, ]: (i) The continuation score and the coefficients of the predictor and corrector are globally Lipschitz, have sufficiently many bounded derivatives in the state and continuation variables, and have at most linear growth. The likelihood coefficients, including λ H_λ, are uniformly bounded on this fixed interval. (i) The corrector intensity γλ _λ is bounded and sufficiently smooth. (i) The exact and numerical processes have uniform moments of sufficiently high order, and the numerical linear systems are uniformly stable for the grids under consideration. Moreover, the total numerical corrector time is uniformly bounded: for some R∗<∞R_*<∞, ∑k=0K−1Mk+1ηk+1≤R∗ _k=0^K-1M_k+1 _k+1≤ R_*. Then, for every φ∈Cb6(ℝd) ∈ C_b^6(R^d), there is a constant Cφ,ΛC_ , independent of h, η∗ _*, and δbud _ bud such that |φ(Kalg)−∫φdμΛ|≤Cφ,Λ(h+η∗+δbud). |E ( z_K alg)- \,d _ |≤ C_ , (h+ _*+ _ bud ). (46) In particular, Algorithm 1 is first-order weakly accurate when η∗=O(h) _*=O(h) and δbud=O(h) _ bud=O(h). Proof. See Appendix D. The three terms in (46) have distinct origins. The term h contains the outer Lie–Trotter and split-step IMEX predictor errors, η∗ _* is the accumulated weak error of the numerical corrector, and δbud _ bud measures the discrepancy between the numerical and continuous corrector times. Corollary 5 (Weak posterior convergence of the discrete algorithm). Assume the conditions of Theorems 3 and 4. Then, for every φ∈Cb6(ℝd) ∈ C_b^6(R^d), |φ(Kalg)−∫φdp(⋅∣)|≤ |E ( z_K alg)- \,dp(· y) |≤ Cφ,λK(h+η∗+δbud) C_ , _K (h+ _*+ _ bud ) (47) +2∥φ∥∞∥μλK−p(⋅∣)∥TV. +2\| \|_∞\| _ _K-p(· y)\|_ TV. Consequently, consider any sequence for which λK→∞ _K→∞, the continuous convergence conditions (43) hold, and, for every φ∈Cb6(ℝd) ∈ C_b^6(R^d), Cφ,λK(h+η∗+δbud)⟶0,C_ , _K (h+ _*+ _ bud ) 0, (48) then the expectations converge for every Cb6C_b^6 test function. Since this class contains the convergence-determining class Cc∞(ℝd)C_c^∞(R^d), the law of the output of Algorithm 1 converges weakly to p(⋅∣)p(· y). Proof. Apply the triangle inequality, Theorem 4, and |∫φd(μ−ν)|≤2‖φ‖∞‖μ−ν‖TV| \,d(μ-ν)|≤ 2\| \|_∞\|μ-ν\|_ TV. The second term in (47) vanishes by Theorem 3, while the first vanishes by (48). Theorem 4 and Corollary 5 assume an exact score. A learned-score error is a separate modeling error and does not generally vanish when the numerical grid is refined. 4 Experiments The experiments follow the construction in Section 3. We first evaluate the complete sampler under a shared protocol. We then use controlled comparisons to isolate scale consistency, noise-conditioned continuation, pathwise correction, stochastic exploration, and innovation placement. This separates end-to-end reconstruction quality from the role of each component. 4.1 Experimental protocol We evaluate on 100 images from FFHQ 256×256256× 256 [20] and 100 images from ImageNet 256×256256× 256 [32]. For FFHQ, we use the pretrained prior from [8]; for ImageNet, we use the unconditional prior from [13]. We consider bicubic ×4× 4 super-resolution (SR), Gaussian deblurring with a 61×6161× 61 kernel of standard deviation 3.03.0, and motion deblurring with a 61×6161× 61 kernel of intensity 0.50.5. Independent Gaussian measurement noise with σd=0.05 _d=0.05 is added in every main task. We report PSNR, SSIM, and LPIPS [49]; higher PSNR and SSIM and lower LPIPS are better. Our main sampler uses 100 score-network evaluations (NFEs), including predictor, corrector, and terminal-denoiser calls. We write P+CP+C for the predictor and corrector allocation. The baselines use 100 outer diffusion or annealing levels, except for the released 1000-step DPS protocol. Since MGPS and DAPS call the score network more than once within an outer level, we report both levels and actual NFEs in Table 1. Images are ordered numerically, and all configurations use indices from 0 to 99 with fixed per-image seeds. The main sampler uses the VM-IMEX predictor, a uniform log-SNR grid, P+C=80+20P+C=80+20, target-change corrector allocation, η=×10−4η=3\!×\!10^-4, and no spectral preconditioner. Its noise-conditioned likelihood covariance is (τ)=σd2+s(τ)⊤,s(τ)=τ4τ2+τc2,τc=6. (τ)= _d^2 I+s(τ) A A , s(τ)= τ^4τ^2+ _c^2, _c=6. (49) We use ζ=0.1ζ=0.1 for FFHQ and ζ=0.01ζ=0.01 for ImageNet; all remaining settings are shared across tasks and datasets. 4.2 Reconstruction results We compare against DPS [8], DiffPIR [50], DDRM [21], RED-Diff [27], MGPS [28], and DAPS [47]. All methods use the same pretrained priors and measurement settings. The comparison is matched in outer levels rather than wall-clock cost. Corrector calls are included in our 100-NFE budget, while repeated score calls internal to MGPS and DAPS are reported explicitly. Method Levels/ NFEs↓ Gaussian deblur Motion deblur SR (×4× 4) PSNR↑ SSIM↑ LPIPS↓ PSNR↑ SSIM↑ LPIPS↓ PSNR↑ SSIM↑ LPIPS↓ FFHQ DPS [8] 1000/1000 25.41 .708 .144 23.69 .658 .178 24.21 .673 .194 DiffPIR [50] 100/100 28.24 .774 .184 27.09 .704 .194 27.67 .772 .118 DDRM [21] 100/100 29.12 .837 .132 29.23 .846 .089 29.30 .843 .141 RED-Diff [27] 100/100 28.46 .705 .242 26.62 .768 .235 27.72 .715 .361 MGPS [28] 100/814 27.88 .792 .149 26.82 .769 .147 27.55 .783 .115 DAPS [47] 100/500 28.71 .792 .155 28.39 .795 .142 28.66 .764 .142 Ours (VM-IMEX) 100/100 29.55 .848 .207 30.57 .869 .137 29.52 .848 .200 ImageNet DPS [8] 1000/1000 21.81 .526 .314 21.15 .503 .352 21.58 .507 .404 DiffPIR [50] 100/100 24.18 .597 .493 23.16 .514 .422 23.71 .577 .338 DDRM [21] 100/100 25.59 .714 .336 26.61 .761 .189 25.64 .720 .321 RED-Diff [27] 100/100 24.52 .661 .491 23.76 .635 .429 24.46 .609 .537 MGPS [28] 100/674 24.68 .657 .397 24.14 .642 .352 24.77 .667 .305 DAPS [47] 100/500 24.98 .661 .349 25.33 .685 .303 25.08 .669 .319 Ours (VM-IMEX) 100/100 25.62 .716 .395 27.41 .764 .180 25.35 .712 .391 Table 1: Reconstruction results with σd=0.05 _d=0.05, averaged over 100 images per dataset. Best and second-best values within each dataset and task are bold and underlined. Computational cost is reported as outer levels/actual score-network evaluations. MGPS uses dataset-dependent inner variational updates, and DAPS uses five reverse-diffusion calls per annealing level. All corrector calls are included in our 100-NFE budget. Our sampler gives the best PSNR and SSIM on all three FFHQ tasks. On ImageNet, it gives the best Gaussian-deblurring PSNR and SSIM and the best motion-deblurring result under all three metrics. DDRM remains strongest on ImageNet SR. The perceptual results are more mixed: DDRM and MGPS retain the best LPIPS on several FFHQ tasks. The shared sampler therefore favors distortion without uniformly dominating the fidelity and perception trade-off. Figure 1 compares every method in Table 1 on the three operators. The images are arranged without intermediate gutters so that structural and textural differences can be compared directly. Measurement DPS DiffPIR DDRM RED-Diff MGPS DAPS Ours Ground truth Gaussian Motion SR Figure 1: Representative FFHQ reconstructions for all methods in Table 1. Columns share the same measurement within each row. 4.3 Ablation study We use FFHQ Gaussian deblurring and ImageNet motion deblurring as representative noisy tasks, followed by noiseless box inpainting as a stiff discretization stress test. The ablations use the following abbreviations. SC Scale-consistent likelihood in the clean-image coordinate. NC Noise-conditioned continuation in (49). PC Frozen-level pathwise Langevin correction. VM Variance-matched innovation placed after the implicit likelihood solve. EX Stochastic exploration controlled by ζ. FFHQ Gaussian deblur ImageNet motion deblur ID Configuration P+CP+C PSNR↑ SSIM↑ LPIPS↓ PSNR↑ SSIM↑ LPIPS↓ SC0 Unscaled likelihood, t-BE 100+0100+0 19.409 .6623 .3179 11.192 .1314 .9364 SC Scale-consistent likelihood, t-BE 100+0100+0 29.409 .8361 .1858 24.969 .5825 .2076 VM0 Direct BE, no NC 100+0100+0 29.488 .8472 .2058 24.799 .5735 .2039 VM VM-IMEX, no NC 100+0100+0 29.502 .8473 .2060 24.729 .5694 .2060 NC VM-IMEX + NC 100+0100+0 29.523 .8472 .2075 27.585 .7786 .1968 PC NC + corrector 80+2080+20 29.547 .8476 .2067 27.406 .7640 .1802 Table 2: Controlled 100-image component ablation. Best and second-best values within each task are bold and underlined. SC0 to SC changes only scale consistency; VM0 to VM changes only innovation placement; VM to NC changes only the likelihood path; NC to PC reallocates 20 of the fixed 100 NFEs to correction. The SC and VM blocks use different coordinate discretizations. SC gives the dominant structural improvement, increasing PSNR by 10.00110.001 dB on FFHQ and 13.77613.776 dB on ImageNet. NC is deliberately separated from SC: it changes FFHQ PSNR by only 0.0210.021 dB but improves ImageNet motion deblurring by 2.8562.856 dB. This confirms that continuation is a conditioning device whose finite-budget benefit depends on the task rather than a universal source of improvement. At the same 100-NFE budget, PC slightly improves FFHQ PSNR and reduces ImageNet LPIPS from .1968.1968 to .1802.1802, while lowering ImageNet PSNR by .179.179 dB. The corrector therefore provides the law-level mechanism used by the convergence theorem, but replacing predictor levels can introduce a finite-budget fidelity and perception trade-off. On these two noisy tasks, VM has only a small and nonuniform metric effect, consistent with Proposition 2: its purpose is covariance retention, not monotone improvement for every operator. Exploration. Figure 2 varies only EX, while Figure 5 shows matched visual samples. The two noisy tasks favor small, task-dependent values of ζ. In box inpainting, the interval 2≤ζ≤122≤ζ≤ 12 produces a partial-refresh regime, followed by recovery and a plateau beyond approximately ζ=64ζ=64. Thus ζ controls a finite-step exploration and exploitation trade-off, but increasing it does not monotonically improve either distortion or perceptual quality. Figure 2: Effect of EX. Each panel reports PSNR on the left axis and LPIPS on the right axis over 100 images. Only ζ changes within a panel. Box inpainting uses missing-square metrics and the two deblurring panels use the complete 80+2080+20 sampler. Split-step stress test. Noiseless centered-box inpainting imposes a singular constraint with a large null space, making the intermediate noise geometry visible. We remove a 128×128128× 128 square from the same 100 FFHQ images and use 100 predictor NFEs, τc=24 _c=24, identical Gaussian draws, and no terminal projection. At ζ=256ζ=256, VM-IMEX obtains 19.41019.410 dB missing-square PSNR and .2195.2195 LPIPS, whereas direct BE obtains 6.5276.527 dB and .8435.8435. The two schemes are otherwise identical. Figure 3 now compares the same two discretizations in both panels. At ζ=0ζ=0 their innovations vanish and the implementations agree. At 100 NFEs the PSNR gap increases with ζ. Refining the grid from 25 to 400 NFEs does not close the gap at ζ=256ζ=256 because the likelihood resolvent remains stiff on the tested grids. Direct BE filters the new observed-coordinate innovation through this resolvent; VM-IMEX injects the same matched innovation after the solve. This explains why the distinction is small in the low-ζ noisy rows of Table 2 but large in the singular box test. Figure 3: Aligned comparison of VM-IMEX and direct BE for noiseless box inpainting. (a) The two schemes use the same 100-NFE grid and differ only in innovation placement. (b) The same pair is compared under grid refinement at ζ=256ζ=256. Each point averages 30 images with paired random draws. Method Levels/NFEs↓ M-PSNR↑ M-SSIM↑ M-LPIPS↓ BGE↓ DPS [8] 1000/1000 18.358 .5280 .1825 .02262 DiffPIR [50] 100/100 19.457 .5889 .2191 .02318 DDRM [21] 100/100 18.199 .5424 .2138 .02775 RED-Diff [27] 100/100 19.038 .5804 .2371 .02757 MGPS [28] 100/814 19.403 .5595 .2257 .01747 DAPS [47] 100/500 18.710 .5462 .1901 .02662 Ours (VM-IMEX) 100/100 19.455 .5888 .2191 .02360 Table 3: Fresh noiseless box-inpainting comparison on the same 100 FFHQ images. M-PSNR, M-SSIM, and M-LPIPS are computed only on the missing 128×128128× 128 square. BGE is the mean absolute error of the reconstructed gradient across the mask boundary. Best and second-best values are bold and underlined. Figure 4: Matched reconstructions for every method in Table 3. From top to bottom, the rows show images 38, 15, and 80. All runs use the same FFHQ prior, centered mask, and noiseless observed pixels. Figure 5: Matched-seed visual effect of EX on image 38. Each row fixes ζ∈0,8,24,256ζ∈\0,8,24,256\, and the five central columns show consecutive samples. The measurement and ground truth are repeated for alignment. All runs use the same 100-NFE VM-IMEX continuation sampler without terminal projection; only ζ changes. Table 3 and Figure 4 extend the stress test to every baseline in Table 1. Ours and DiffPIR have nearly identical missing-region distortion, differing by less than .002.002 dB in PSNR and 10−410^-4 in SSIM. DPS and DAPS obtain lower missing-region LPIPS, while MGPS gives the smallest boundary-gradient error. The comparison therefore exposes a fidelity, perception, and boundary trade-off rather than a uniform ranking. It also separates two claims: the paired VM-IMEX comparison above identifies the effect of innovation placement within our discretization, whereas the cross-method table includes additional differences in posterior approximation, optimization, and computational cost. Ours uses no terminal projection in this experiment; the missing-region metrics are unaffected by the observed-pixel replacement used by some native baselines. Figure 5 complements the aggregate sweep in Figure 2. At ζ=0ζ=0, the outputs are similar and smooth; ζ=8ζ=8 produces unstable high-frequency patterns across all five samples. Setting ζ=24ζ=24 restores coherent faces with visible variation, and ζ=256ζ=256 yields a similar stable regime. Thus, the matched visualization exposes a nonmonotone finite-step transition rather than estimating diversity. 5 Conclusion We introduced a scale-consistent predictor-corrector framework for diffusion inverse problems. The construction places the learned prior and measurement likelihood in a common rescaled coordinate. Within this framework, noise-conditioned continuation accounts for projected diffusion uncertainty in the likelihood covariance. A variance-matched split-step predictor stabilizes the likelihood drift without filtering the stochastic innovation, while Langevin corrections track the evolving target. This viewpoint separates three questions that are often conflated: whether the target approaches the posterior, whether the continuous sampler tracks that target, and whether a finite-step discretization approximates the continuous dynamics. Our analysis establishes corresponding consistency, convergence, and weak-error results. Under a fixed 100-NFE budget, the method achieves strong reconstruction fidelity on FFHQ and ImageNet, particularly for motion deblurring. The controlled ablations show that scale consistency is the dominant structural requirement, while continuation is strongly task dependent. Allocating part of the fixed budget to corrector steps can improve perceptual quality at the expense of predictor resolution and distortion metrics. On the noisy deblurring tasks studied here, moving the matched innovation outside the likelihood resolvent has only a small aggregate effect. The separate noiseless box study explains this contrast: increasing ζ improves reconstruction until residual refresh saturates, whereas filtering that same large innovation through the stiff mask resolvent destroys the noise geometry seen by the next denoiser. Cross-task ζ sweeps further show that its finite-budget distortion–perception optimum is operator and prior dependent. The aligned stochasticity and grid-refinement comparisons connect the VM-IMEX advantage to attenuation of the matched innovation by an unresolved implicit likelihood solve. Appendix A Proof of the Posterior SDE Family Proof of Proposition 1. Write the forward VP drift as (,t):=−12β(t). b( x,t):=- 12β(t) x. Conditioning the initial distribution on the fixed observation y does not change the forward transition kernel. In particular, for t>0t>0, pt():=pt(∣)=∫(,αt0,σt2)p0(d0∣).p_t y( x):=p_t( x y)= \! ( x; _t x_0, _t^2 I )p_0(d x_0 y). Since σt2>0 _t^2>0 for every t>0t>0, Gaussian smoothing makes ptp_t y strictly positive and smooth, even when p0(⋅∣)p_0(· y) does not admit a density. The density ptp_t y satisfies the forward Fokker–Planck equation ∂tpt=−∇⋅((⋅,t)pt)+β(t)2Δpt. _tp_t y=-∇\!·\! ( b(·,t)p_t y )+ β(t)2 p_t y. (50) Let qtq_t denote the marginal density of the reverse-time process in (12), and define its drift by ζ(,t):=(,t)−1+ζ2β(t)∇logpt(). b_ζ( x,t):= b( x,t)- 1+ζ2β(t) _ x p_t y( x). Because the process is evolved from T to 00, its Fokker–Planck equation, written in the original time variable t, has a negative diffusion term: ∂tqt=−∇⋅(ζ(⋅,t)qt)−ζβ(t)2Δqt. _tq_t=-∇\!·\! ( b_ζ(·,t)q_t )- ζβ(t)2 q_t. (51) Indeed, this equation follows directly by introducing the increasing reverse clock s=T−ts=T-t and applying the standard forward Fokker–Planck equation to T−s x_T-s. Fix ε∈(0,T) ∈(0,T). We verify that qt=ptq_t=p_t y solves (51) on [ε,T][ ,T]. Using pt∇logpt=∇ptp_t y∇ p_t y=∇ p_t y gives −∇⋅(ζ(⋅,t)pt)−ζβ(t)2Δpt -∇\!·\! ( b_ζ(·,t)p_t y )- ζβ(t)2 p_t y =−∇⋅((⋅,t)pt)+1+ζ2β(t)Δpt−ζ2β(t)Δpt =-∇\!·\! ( b(·,t)p_t y )+ 1+ζ2β(t) p_t y- ζ2β(t) p_t y =−∇⋅((⋅,t)pt)+β(t)2Δpt=∂tpt, =-∇\!·\! ( b(·,t)p_t y )+ β(t)2 p_t y= _tp_t y, (52) where the last equality follows from (50). Since qT=pTq_T=p_T y by initialization, uniqueness of the Fokker–Planck solution implies qt=ptq_t=p_t y for all t∈[ε,T]t∈[ ,T]. Since ε>0 >0 is arbitrary, the equality holds for every t∈(0,T]t∈(0,T]. Finally, continuity of the forward VP semigroup at t=0t=0 gives pt(⋅∣)⇒p0(⋅∣)p_t(· y) p_0(· y) as t↓0t 0. Appendix B Proof of Terminal Target and Score Consistency Proof of the endpoint consistency statement in Subsection 3.2. The rescaled prior density is the Gaussian convolution πλ=p0∗(,e−2λ). _λ=p_0*N( 0,e^-2λ I). Since the Gaussian kernels form an approximate identity, ‖πλ−p0‖1→0\| _λ-p_0\|_1→ 0. Define L∞():=exp(−12σd2‖−‖22).L_∞( z):= \! (- 12 _d^2\| A z- y\|_2^2 ). Because λ−1→σd−2 _λ^-1→ _d^-2 I, we have Lλ()→L∞()L_λ( z)→ L_∞( z) pointwise, and both likelihood factors take values in (0,1](0,1]. Consequently, ‖πλLλ−p0L∞‖1 \| _λL_λ-p_0L_∞\|_1 ≤‖πλ−p0‖1+∫πλ|Lλ−L∞| ≤\| _λ-p_0\|_1+ _λ|L_λ-L_∞|\,d z ≤2‖πλ−p0‖1+∫p0|Lλ−L∞|⟶0, ≤ 2\| _λ-p_0\|_1+ p_0|L_λ-L_∞|\,d z 0, (53) where the last term converges by dominated convergence. Hence Zλ→Z∞:=∫p0L∞>0Z_λ→ Z_∞:= p_0L_∞\,d z>0. Normalizing the two unnormalized densities gives ∥π¯λ(⋅∣)−p(⋅∣)∥1⟶0,\| π_λ(· y)-p(· y)\|_1 0, which proves the stated total-variation convergence. For the score statement, the Gaussian measurement model gives ∇logp(∣)=−1σd2⊤(−). _ z p( y z)=- 1 _d^2 A ( A z- y). Bayes’ rule therefore yields ∇logp(∣)=∇logp0()−1σd2⊤(−). _ z p( z y)= _ z p_0( z)- 1 _d^2 A ( A z- y). Subtracting this identity from (22) gives ¯λ(,)−∇logp(∣)=∇logπλ()−∇logp0()−⊤(λ−1−σd−2)(−). s_λ( z; y)- _ z p( z y)= _ z _λ( z)- _ z p_0( z)- A ( _λ^-1- _d^-2 I)( A z- y). The first difference converges locally uniformly by assumption. The second difference does so because λ−1→σd−2 _λ^-1→ _d^-2 I and − A z- y is bounded on compact subsets. This proves the stated score convergence. Appendix C Proof of Continuous Posterior Convergence Proof of Theorem 3. Write qλ=e−2λq_λ=e^-2λ, π¯=π¯λ π= π_λ, μ=μλμ= _λ, and g=μ/π¯g=μ/ π. The Fokker–Planck equation of (26) can be written as ∂λμ= _λμ= −qλ∇⋅(μ∇logπ¯)+ζqλ∇⋅(μ∇logg) -q_λ∇\!·(μ∇ π)+ζ q_λ∇\!·(μ∇ g) (54) +γλ∇⋅(μ∇logg). + _λ∇\!·(μ∇ g). Indeed, the first two terms on the right equal the predictor Fokker–Planck operator because (1+ζ)qλ=qλ+ζqλ(1+ζ)q_λ=q_λ+ζ q_λ, while the last term is the frozen-target corrector operator. By (37), ∂λπ¯=−qλΔπ¯+rλπ¯. _λ π=-q_λ π+r_λ π. (55) Differentiating H(λ)=∫μloggH(λ)= μ g and integrating by parts in (54) gives H′(λ)= H (λ)= −ζqλ∫∥∇logg∥22dμ−γλ∫∥∇logg∥22dμ -ζ q_λ \|∇ g\|_2^2\,dμ- _λ \|∇ g\|_2^2\,dμ (56) −∫rλdμ. - r_λ\,dμ. For completeness, the contribution of the first term in (54) is qλ∫∇logπ¯⋅∇loggdμ=−qλ∫Δπ¯π¯dμ.q_λ ∇ π·∇ g\,dμ=-q_λ π π\,dμ. It cancels the corresponding term from −∫∂λlogπ¯dμ- _λ π\,dμ by (55), leaving the final defect term in (56). This cancellation is also why rλr_λ, rather than a numerical local error, is the relevant continuous tracking quantity. Let ℐ(μ∥π¯):=∫∥∇logg∥22dμ.I(μ\,\|\, π):= \|∇ g\|_2^2\,dμ. Applying (41) to f=gf= g gives ℐ(μ∥π¯)≥2ρλH(λ).I(μ\,\|\, π)≥ 2 _λH(λ). Dropping the nonpositive first term in (56) and using (39), we obtain H′(λ)≤−(2γλρλ−cλ)H(λ)+ελ=−kλH(λ)+ελ.H (λ)≤- (2 _λ _λ-c_λ )H(λ)+ _λ=-k_λH(λ)+ _λ. The integrating-factor form of Grönwall’s inequality is precisely (42). Conditions (43) therefore imply H(λ)→0H(λ)→ 0. The simpler sufficient condition stated in the theorem follows by splitting the convolution at a fixed large value of s and using kλ≥k∗>0k_λ≥ k_*>0 on the remaining tail. Finally, Pinsker’s inequality, the triangle inequality, and the endpoint consistency stated in Subsection 3.2 yield ∥μλ−p(⋅∣)∥TV \| _λ-p(· y)\|_ TV ≤∥μλ−π¯λ∥TV+∥π¯λ−p(⋅∣)∥TV ≤\| _λ- π_λ\|_ TV+\| π_λ-p(· y)\|_ TV ≤12H(λ)+∥π¯λ−p(⋅∣)∥TV⟶0, ≤ 12H(λ)+\| π_λ-p(· y)\|_ TV 0, which proves (44). Appendix D Proof of the Weak Error Bound for the Complete Scheme Proof of Theorem 4. Let kQ_k denote the exact transition kernel of the continuously interleaved SDE (26) from λk _k to λk+1 _k+1. Let ℛkpredR_k pred and λk+1rC_ _k+1^r be the exact predictor and frozen-corrector kernels introduced before (29), and define the exact split kernel k:=ℛkpredλk+1rk+1.S_k:=R_k predC_ _k+1^r_k+1. We write ℛ^k R_k for the semi-implicit predictor kernel in (31) and ^λη C_λ^η for one semi-implicit corrector step in (33). The one-level transition of Algorithm 1 is therefore ^k:=ℛ^k(^λk+1ηk+1)Mk+1. S_k:= R_k ( C_ _k+1 _k+1 )^M_k+1. We first collect the local weak estimates. Under the smoothness and moment assumptions of the theorem, the Duhamel expansion for the two generators gives the local Lie–Trotter estimate |f(Xk)−f(Xk)|≤Cf(1+‖X‖q)hk2. |Ef(XQ_k)-Ef(XS_k) |≤ C_f (1+E\|X\|^q )h_k^2. (57) Here and below, XXK denotes a random variable whose law is obtained by applying the kernel K to the law of X. A stochastic Taylor expansion of the predictor gives |f(Xℛkpred)−f(Xℛ^k)|≤Cf(1+‖X‖q)hk2. |Ef(XR_k pred)-Ef(X R_k) |≤ C_f (1+E\|X\|^q )h_k^2. (58) To verify this order for the variance-matched split step, use τk+1=τke−hk _k+1= _ke^-h_k in (30) to obtain ak a_k =1−(1+ζ)hk+O(hk2), =1-(1+ζ)h_k+O(h_k^2), (59) ωk _k =(1+ζ)τk2hk+O(hk2), =(1+ζ) _k^2h_k+O(h_k^2), qkvm q_k vm =2ζτk2hk+O(hk2). =2ζ _k^2h_k+O(h_k^2). These are precisely the drift and diffusion coefficients of the predictor generator frozen at λk _k, up to O(hk2)O(h_k^2). Moreover, k=+O(hk) R_k= I+O(h_k) under the bounded-coefficient assumption. Hence moving an O(hk)O(h_k) innovation covariance from inside the likelihood resolvent to after the solve changes its covariance only by O(hk2)O(h_k^2). Taylor expansion of the propagated test function therefore yields (58). This is a fixed-coefficient consistency statement; it is not uniform in a regime where ωk‖k+1‖ _k\| H_k+1\| remains order one as the grid changes. Similarly, the centered implicit treatment of the corrector likelihood and the explicit treatment of its smooth prior score give |f(Xλη)−f(X^λη)|≤Cf(1+‖X‖q)η2. |Ef(XC_λ^η)-Ef(X C_λ^η) |≤ C_f (1+E\|X\|^q )η^2. (60) The constants are uniform on the fixed interval [λ0,Λ][ _0, ] for the propagated test functions arising below. For the corrector at level k+1k+1, telescoping (60) over Mk+1M_k+1 steps shows that replacing the exact corrector of duration Mk+1ηk+1M_k+1 _k+1 by the numerical corrector costs at most CMk+1ηk+12CM_k+1 _k+1^2. Smooth dependence of the corrector semigroup on its auxiliary time also gives |f(Xλk+1Mk+1ηk+1)−f(Xλk+1rk+1)| |Ef(XC_ _k+1^M_k+1 _k+1)-Ef(XC_ _k+1^r_k+1) | (61) ≤Cf(1+‖X‖q)|Mk+1ηk+1−rk+1|. ≤ C_f (1+E\|X\|^q ) |M_k+1 _k+1-r_k+1 |. We now telescope the exact kernels k\Q_k\ against the numerical kernels ^k\ S_k\ over the full interval. The backward Kolmogorov regularity and the assumed moment stability allow (57)–(61) to be applied to every propagated test function. Consequently, |φ(Kalg)−∫φdμΛ| |E ( z_K alg)- \,d _ | ≤Cφ,Λ[∑k=0K−1hk2+∑k=0K−1Mk+1ηk+12+∑k=0K−1|Mk+1ηk+1−rk+1|]. ≤ C_ , [ _k=0^K-1h_k^2+ _k=0^K-1M_k+1 _k+1^2+ _k=0^K-1 |M_k+1 _k+1-r_k+1 | ]. (62) Since ∑khk=Λ−λ0 _kh_k= - _0, ∑khk2≤h(Λ−λ0). _kh_k^2≤ h( - _0). Moreover, ∑kMk+1ηk+12≤η∗∑kMk+1ηk+1≤R∗η∗. _kM_k+1 _k+1^2≤ _* _kM_k+1 _k+1≤ R_* _*. Substitution into (62) proves (46). References [1] A. Abdulle, G. Vilmart, and K. C. Zygalakis (2015) Long time accuracy of Lie–Trotter splitting methods for Langevin dynamics. SIAM Journal on Numerical Analysis 53 (1), p. 1–16. External Links: Document Cited by: §3.3. [2] E. L. Allgower and K. Georg (2003) Introduction to numerical continuation methods. Classics in Applied Mathematics, Vol. 45, Society for Industrial and Applied Mathematics, Philadelphia, PA. External Links: Document Cited by: §2. [3] B. D. Anderson (1982) Reverse-time diffusion equation models. Stochastic Processes and their Applications 12 (3), p. 313–326. Cited by: §1. [4] J. Batson and L. Royer (2019) Noise2Self: blind denoising by self-supervision. In Proceedings of the International Conference on Machine Learning, p. 524–533. Cited by: §1. [5] A. Buades, B. Coll, and J. Morel (2005) A non-local algorithm for image denoising. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, Vol. 2, p. 60–65. Cited by: §1. [6] G. Cardoso, Y. J. E. Idrissi, S. L. Corff, and E. Moulines (2023) Monte carlo guided diffusion for bayesian linear inverse problems. arXiv preprint arXiv:2308.07983. Cited by: §2. [7] T. Chan, A. Marquina, and P. Mulet (2000) High-order total variation-based image restoration. SIAM Journal on Scientific Computing 22 (2), p. 503–516. Cited by: §1. [8] H. Chung, J. Kim, M. T. Mccann, M. L. Klasky, and J. C. Ye (2022) Diffusion posterior sampling for general noisy inverse problems. arXiv preprint arXiv:2209.14687. Cited by: §2, §4.1, §4.2, Table 1, Table 1, Table 3. [9] H. Chung, B. Sim, D. Ryu, and J. C. Ye (2022) Improving diffusion models for inverse problems using manifold constraints. In Advances in Neural Information Processing Systems, Vol. 35, p. 25683–25696. Cited by: §2. [10] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian (2007) Image denoising by sparse 3-d transform-domain collaborative filtering. IEEE Transactions on Image Processing 16 (8), p. 2080–2095. Cited by: §1. [11] G. Daras, H. Chung, C. Lai, Y. Mitsufuji, J. C. Ye, P. Milanfar, A. G. Dimakis, and M. Delbracio (2024) A survey on diffusion models for inverse problems. External Links: 2410.00083, Document Cited by: §2. [12] P. Del Moral, A. Doucet, and A. Jasra (2006) Sequential Monte Carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68 (3), p. 411–436. External Links: Document Cited by: §2. [13] P. Dhariwal and A. Nichol (2021) Diffusion models beat gans on image synthesis. Advances in neural information processing systems 34, p. 8780–8794. Cited by: §4.1. [14] Z. Dou and Y. Song (2024) Diffusion posterior sampling for linear inverse problem solving: a filtering perspective. In The Twelfth International Conference on Learning Representations, Cited by: §2. [15] S. Gu, L. Zhang, W. Zuo, and X. Feng (2014) Weighted nuclear norm minimization with application to image denoising. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 2862–2869. Cited by: §1. [16] D. J. Higham, X. Mao, and A. M. Stuart (2002) Strong convergence of euler-type methods for nonlinear stochastic differential equations. SIAM Journal on Numerical Analysis 40 (3), p. 1041–1063. External Links: Document Cited by: §3.3. [17] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems 33, p. 6840–6851. Cited by: §1. [18] H. Ji, C. Liu, Z. Shen, and Y. Xu (2010) Robust video denoising using low rank matrix completion. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 1791–1798. Cited by: §1. [19] U. S. Kamilov, C. A. Bouman, G. T. Buzzard, and B. Wohlberg (2023) Plug-and-play methods for integrating physical and learned models in computational imaging: theory, algorithms, and applications. IEEE Signal Processing Magazine 40 (1), p. 85–97. Cited by: §2. [20] T. Karras, S. Laine, and T. Aila (2019) A style-based generator architecture for generative adversarial networks. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, p. 4401–4410. Cited by: §4.1. [21] B. Kawar, M. Elad, S. Ermon, and J. Song (2022) Denoising diffusion restoration models. External Links: 2201.11793 Cited by: §2, §4.2, Table 1, Table 1, Table 3. [22] B. Kawar, G. Vaksman, and M. Elad (2021) Snips: solving noisy inverse problems stochastically. Advances in Neural Information Processing Systems 34, p. 21757–21769. Cited by: §2. [23] S. Laine, T. Karras, J. Lehtinen, and T. Aila (2019) High-quality self-supervised deep image denoising. In Advances in Neural Information Processing Systems 32 (NeurIPS 2019), p. . Cited by: §1. [24] J. Lehtinen, J. Munkberg, J. Hasselgren, S. Laine, T. Karras, M. Aittala, and T. Aila (2018) Noise2Noise: learning image restoration without clean data. In Proceedings of the International Conference on Machine Learning, p. 2971–2980. Cited by: §1. [25] T. Li, A. Abdulle, and W. E (2008) Effectiveness of implicit methods for stiff stochastic differential equations. Communications in Computational Physics 3 (2), p. 295–307. External Links: Document Cited by: §3.3, §3.3. [26] J. Liang, J. Cao, G. Sun, K. Zhang, L. Van Gool, and R. Timofte (2021) SwinIR: image restoration using swin transformer. In Proceedings of the IEEE/CVF International Conference on Computer Vision, p. 1833–1844. Cited by: §1. [27] M. Mardani, J. Song, J. Kautz, and A. Vahdat (2023) A variational perspective on solving inverse problems with diffusion models. arXiv preprint arXiv:2305.04391. Cited by: §2, §4.2, Table 1, Table 1, Table 3. [28] B. Moufad, Y. Janati, L. Bedin, A. Durmus, R. Douc, E. Moulines, and J. Olsson (2024) Variational diffusion posterior sampling with midpoint guidance. arXiv preprint arXiv:2410.09945. Cited by: §2, §4.2, Table 1, Table 1, Table 3. [29] R. M. Neal (2001) Annealed importance sampling. Statistics and Computing 11 (2), p. 125–139. External Links: Document Cited by: §2. [30] Y. Romano, M. Elad, and P. Milanfar (2017) The little engine that could: regularization by denoising (red). SIAM Journal on Imaging Sciences 10 (4), p. 1804–1844. Cited by: §2. [31] L. I. Rudin, S. Osher, and E. Fatemi (1992) Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena 60 (1–4), p. 259–268. Cited by: §1. [32] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, et al. (2015) Imagenet large scale visual recognition challenge. International journal of computer vision 115, p. 211–252. Cited by: §4.1. [33] H. Schurz (1999) Preservation of probabilistic laws through euler methods for Ornstein–Uhlenbeck process. Stochastic Analysis and Applications 17 (3), p. 463–486. External Links: Document Cited by: §3.3, §3.3. [34] J. Song, A. Vahdat, M. Mardani, and J. Kautz (2023) Pseudoinverse-guided diffusion models for inverse problems. In International Conference on Learning Representations, Cited by: §2. [35] Y. Song and S. Ermon (2019) Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems 32. Cited by: §2. [36] Y. Song, L. Shen, L. Xing, and S. Ermon (2022) Solving inverse problems in medical imaging with score-based generative models. In International Conference on Learning Representations, Cited by: §2. [37] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole (2020) Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456. Cited by: §1. [38] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole (2020) Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456. Cited by: §1, §1, §2. [39] D. Strong and T. Chan (2003) Edge-preserving and scale-dependent properties of total variation regularization. Inverse Problems 19 (6), p. S165. Cited by: §1. [40] Y. Sun, Z. Wu, Y. Chen, B. T. Feng, and K. L. Bouman (2024) Provable probabilistic imaging using score-based generative priors. IEEE Transactions on Computational Imaging. Cited by: §2. [41] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg (2013) Plug-and-play priors for model based reconstruction. In 2013 IEEE global conference on signal and information processing, p. 945–948. Cited by: §2. [42] Y. Wang, J. Yu, and J. Zhang (2022) Zero-shot image restoration using denoising diffusion null-space model. External Links: 2212.00490 Cited by: §2. [43] L. Wu, B. Trippe, C. Naesseth, D. Blei, and J. P. Cunningham (2023) Practical and asymptotically exact conditional sampling in diffusion models. Advances in Neural Information Processing Systems 36, p. 31372–31403. Cited by: §2. [44] Z. Wu, Y. Sun, Y. Chen, B. Zhang, Y. Yue, and K. Bouman (2024) Principled probabilistic imaging using diffusion models as plug-and-play priors. Advances in Neural Information Processing Systems 37, p. 118389–118427. Cited by: §2. [45] T. Xu, Z. Zhu, J. Li, D. He, Y. Wang, M. Sun, L. Li, H. Qin, Y. Wang, J. Liu, et al. (2024) Consistency model is an effective posterior sample approximation for diffusion inverse solvers. arXiv preprint arXiv:2403.12063. Cited by: §2. [46] S. W. Zamir, A. Arora, S. Khan, M. Hayat, F. S. Khan, and M. Yang (2022) Restormer: efficient transformer for high-resolution image restoration. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, p. 5728–5739. Cited by: §1. [47] B. Zhang, W. Chu, J. Berner, C. Meng, A. Anandkumar, and Y. Song (2025) Improving diffusion inverse problem solving with decoupled noise annealing. In Proceedings of the Computer Vision and Pattern Recognition Conference, p. 20895–20905. Cited by: §2, §2, §4.2, Table 1, Table 1, Table 3. [48] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang (2017) Beyond a gaussian denoiser: residual learning of deep cnn for image denoising. IEEE Transactions on Image Processing 26 (7), p. 3142–3155. Cited by: §1. [49] R. Zhang, P. Isola, A. A. Efros, E. Shechtman, and O. Wang (2018) The unreasonable effectiveness of deep features as a perceptual metric. In Proceedings of the IEEE conference on computer vision and pattern recognition, p. 586–595. Cited by: §4.1. [50] Y. Zhu, K. Zhang, J. Liang, J. Cao, B. Wen, R. Timofte, and L. V. Gool (2023) Denoising diffusion models for plug-and-play image restoration. External Links: 2305.08995 Cited by: §2, §4.2, Table 1, Table 1, Table 3.