Provably Robust Score-Based Diffusion Posterior Sampling for Plug-and-Play Image Reconstruction

Authors: xingyu  xu
yuejie  chi
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 15:43:20

Abstract: In a great number of tasks in science and engineering, the goal is to infer an unknown image from a small number of noisy measurements collected from a known forward model describing certain sensing or imaging modality. Due to resource constraints, this image reconstruction task is often extremely ill-posed, which necessitates the adoption of expressive prior information to regularize the solution space. Score-based diffusion models, thanks to its impressive empirical success, have emerged as an appealing candidate of an expressive prior in image reconstruction. In order to accommodate diverse tasks at once, it is of great interest to develop efficient, consistent and robust algorithms that incorporate unconditional score functions of an image prior distribution in conjunction with flexible choices of forward models.

This work develops an algorithmic framework for employing score-based diffusion models as an expressive data prior in nonlinear inverse problems with general forward models. Motivated by the plug-and-play framework in the imaging community, we introduce a diffusion plug-and-play method (DPnP) that alternatively calls two samplers, a proximal consistency sampler based solely on the likelihood function of the forward model, and a denoising diffusion sampler based solely on the score functions of the image prior. The key insight is that denoising under white Gaussian noise can be solved rigorously via both stochastic (i.e., DDPM-type) and deterministic (i.e., DDIM-type) samplers using the same set of score functions trained for generation. We establish both asymptotic and non-asymptotic performance guarantees of DPnP, and provide numerical experiments   to illustrate its promise in solving both linear and nonlinear image reconstruction tasks. To the best of our knowledge, DPnP is the first provably-robust posterior sampling method for nonlinear inverse problems using unconditional diffusion priors.

Main paper: 
 introduction: 
 In a great number of sensing and imaging applications, the paramount goal is to infer an unknown image x ⋆ ∈ R d from a collection of measurements y ∈ R m that are possibly noisy, incomplete, and even nonlinear. Examples include restoration tasks such as inpainting, super-resolution, denoising, as well as imaging tasks such as magnetic resonance imaging [LDP07] , optical imaging [SEC + 15], microscopy imaging [HSMC17] , radar and sonar imaging [PEPC10] , and many more.Due to sensing and resource constraints, the problem of image reconstruction is often ill-posed, where the desired resolution of the unknown image overwhelms the set of available observations. Consequently, this necessitates the need of incorporating prior information regarding the unknown image to assist the reconstruction process. Over the years, numerous types of prior information have been considered and adopted, from hand-crafted priors such as subspace or sparsity constraints [Don06, CR12] , to data-driven ones prescribed in the form of neural networks [UVL18, BJPD17] . These priors can be regarded as some sort of generative models for the unknown image, which postulate the high-dimensional image admits certain parsimonious representation in a low-dimensional data manifold. It is desirable that the generative models are sufficiently expressive to capture the Super-resolution (linear) Phase retrieval (nonlinear) Quantized sensing (nonlinear) diversity and structure of the image class of interest, yet nonetheless, still lead to image reconstruction problems that are computationally tractable.Score-based diffusion models as an image prior. Recent years have seen tremendous progress on generative artificial intelligence (AI), where it is possible to generate new data samples -such as images, audio, text -at unprecedented resolution and scale from a target distribution given training data. Diffusion models, originally proposed by [SDWMG15] , are among one of the most successful frameworks, underneath popular content generators such as DALL•E [RDN + 22], Stable Diffusion [RBL + 22], Imagen [SCS + 22], and many others. Roughly speaking, score-based diffusion models convert noise into samples that resemble those from a target data distribution, by forming the reverse Markov diffusion process only using the score functions of the data contaminated at various noise levels [SE19, HJA20, SSDK + 21, SME20]. In particular, [SME20] developed a unified framework to interpret score-based diffusion models as reversing certain Stochastic Differential Equations (SDE) using either SDE or probability flow Ordinary Differential Equations (ODE), leading to stochastic (i.e., DDPM-type) and deterministic (i.e., DDIM-type) samplers, respectively. While the DDIM-type sampler is more amenable to acceleration, the DDPM-type sampler tends to generate images of higher quality and diversity when running for a large number of steps [SME20] .Thanks to the expressive power of score-based diffusion models in generating complex and finegrained images, they have emerged as a plausible candidate of an expressive prior in image reconstruction [SSXE21, CKM + 23, FSR + 23] via the lens of Bayesian posterior sampling. To accommodate diverse applications with various image characteristics and imaging modalities, it is desirable to develop plug-and-play methods that do not require training from scratch or end-to-end training for every new imaging task. Nonetheless, despite a flurry of recent efforts, existing algorithms either are computationally expensive [WTN + 23, CICM23], inconsistent [CKM + 23, KEES22, MSKV24], or confined to linear inverse problems [CICM23, DS24] . Therefore, a natural question arises:Can we develop a practical, consistent and robust algorithm that incorporates score-based diffusion models as an image prior with general (possibly nonlinear) forward models?This paper provides an affirmative answer to this question, by developing an algorithmic framework to sample from the posterior distribution of images, where score-based diffusion models are employed as an expressive image prior in nonlinear inverse problems with general forward models. Specifically, our contributions are as follows.• Diffusion plug-and-play for posterior sampling. Motivated by the plug-and-play [VBW13] framework in the imaging community, we introduce a diffusion plug-and-play method (DPnP) that alternatively calls two samplers, a proximal consistency sampler that aims to generate samples that are more consistent with the measurements, and a denoising diffusion sampler that focuses on sampling from the posterior distribution of an easier problem -image denoising under white Gaussian noise -to enforce the prior constraint. Our method is modular, in the sense that the proximal consistency sampler is solely based on the likelihood function of the forward model, and the denoising diffusion sampler is based solely on the score functions of the image prior. • Posterior sampling for image denoising. While the proximal consistency sampler can be borrowed somewhat straightforwardly from existing literature such as the Metropolis-adjusted Langevin algorithm [RR98] , the denoising diffusion sampler, on the other hand, has not been addressed in the literature to the best of our knowledge. Our key insight is that this can be solved via both stochastic (i.e., DDPM-type) or deterministic (i.e., DDIM-type) samplers by carefully choosing the forward SDEs and discretizing the resulting reversal SDE or ODE using the exponential integrator [ZC22] . Importantly, the denoising diffusion samplers use the same set of unconditional score functions for generation, making it readily implementable without additional training.• Theoretical guarantees. We establish both asymptotic and non-asymptotic performance guarantees of the proposed DPnP method. Asymptotically, we verify the correctness of our method by proving that DPnP converges to the conditional distribution of x ⋆ given measurements y, assuming exact unconditional score estimates of the image prior. We next establish a non-asymptotic convergence theory of DPnP, where its performance degenerates gracefully with respect to the errors of the samplers, due to, e.g., score estimation errors and limited sampling steps. To the best of our knowledge, this provides the first provably-robust method for nonlinear inverse problems using unconditional score-based diffusion priors.We further provide numerical experiments to illustrate its promise in solving both linear and nonlinear image reconstruction tasks, such as super-resolution, phase retrieval, and quantized sensing. Due to its plug-and-play nature, we expect it to be of broad interest to a wide variety of inverse problems.Related works. Given its interdisciplinary nature, our work sits at the intersection of generative modeling, computational imaging, optimization and sampling. Due to space limits, we postpone the discussion of related works to Appendix A.Notation. Let p x denote the probability distribution of x, and p x (•|y) denotes the conditional distribution of x given y. We useX (d)= Y to denote random variables X and Y are equivalent in distribution. The matrix I d denotes an identity matrix of dimension d. For two probability distributions with density p(x) and q(x), the total variation distance between them is TV(p, q) := |p(x) -q(x)|dx. The χ 2 -divergence of p to q is χ 2 (p ∥ q) := (p(x)-q(x)) 2 q(x)dx.2 Score-based generative modelsIn this section, we set up the preliminary on diffusion-based generative models, which we will be relying upon to develop our algorithm. The key components consist of a forward process, which diffuses the data distribution p ⋆ to the standard normal distribution by gradually injecting noise into the samples, and a backward process, which reverses the forward process so that it can transform the standard normal distribution to the data distribution p ⋆ . To facilitate understanding, it will be convenient to formulate these processes in continuous time. For discrete-time formulation and implementation, please refer to Appendix E.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: The abstract and introduction provides explicit description of the contributions and the context.• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
the forward process and score functions: 
 The continuous-time forward diffusion follows the Ornstein-Uhlenbeck (OU) process, defined by the Stochastic Differential Equation (SDE) [SSDK + 21]:EQUATIONwhere (B τ ) τ ≥0 is the standard d-dimensional Brownian motion. It can be shown that [Doo42, Eva12] the marginal distribution of X τ for τ ≥ 0 isEQUATIONIt is then clear that the limiting distribution X ∞ ∼ N (0, I d ) as τ → ∞, i.e., the OU process diffuses X 0 ∼ p ⋆ to the standard normal distribution. The score function of X τ is defined byEQUATIONAn enlightening property [Vin11] of the score function is that it can be interpreted as the minimum mean-squared error (MMSE) estimate of ε t given x t = x, fueled by Tweedie's formula:EQUATIONConsequently, this makes it possible to estimate the score functions via learning to denoise [Hyv05] , by estimating the denoising function ε(τ, •), as typically done in practice [HJA20] .To enable sampling, one needs to "reverse" the forward diffusion process. Fortunately, it is possible to leverage classical theory [And82, AGS05] to reverse the SDE, and apply discretization to the timereversal processes to collect samples. We shall describe two popular approaches below, corresponding to stochastic (i.e., DDPM-type) and deterministic (i.e., DDIM-type) samplers respectively following primarily the framework set forth in [SSDK + 21].Time-reversed SDEs and probability flow ODEs. Let us begin with the more general theory of reversing SDEs, which will be useful in future sections. Consider a SDE given byEQUATIONwhere α ∈ R and β > 0 are constants. For any positive time τ ∞ > 0, define the reversed time parameter τ rev := τ rev (τ ) = τ ∞ -τ.(6) We are now ready to describe the time-reversed processes.1) The time-reversed SDE of (5) on the time interval [0, τ ∞ ] is defined asEQUATION)where B is an independent copy of B, i.e., another Brownian motion. It is a classical result [And82] that the reversed process M rev shares the same path distribution as M , i.e., (M rev τ ) τ ∈[0,τ∞](d) = (M τ ) τ ∈[0,τ∞] . In other words, the joint distribution of (M rev τ1 , M rev τ2 , • • • , M rev τ k ) for any 0 ≤ τ 1 ≤ τ 2 ≤ • • • ≤ τ k ≤ τ ∞ , for any integer k ≥ 1, coincides with that of (M τ1 , M τ2 , • • • , M τ k ).2) In place of the reversed SDE in (7), it is possible to consider the following probability flow ODE [AGS05, SSDK + 21]:dM rev τ rev = -αM rev τ rev + β 2 ∇ log p M τ rev (τ rev , M rev τ rev ) dτ, τ ∈ [0, τ ∞ ], M rev τ∞ ∼ p Mτ ∞ . (8)The reversed ODE satisfies a slightly weaker guarantee than that of the reversed SDE, which nevertheless suffices for most practical purposes [SSDK + 21]:M rev τ (d) = M τ , τ ∈ [0, τ ∞ ].Note that the reversed ODE only guarantees identical marginal distribution for each M rev τ , whereas the reversed SDE guarantees identical joint distribution.Specializing the above to the OU process (1) with proper discretization then leads to popular samplers used for generation, as follows.DDPM-type stochastic samplers. Specializing the time-reversed SDE (7) to the OU process givesdX rev τ rev = X rev τ rev + 2s(τ rev , X rev τ rev ) dτ + √ 2d Bτ , τ ∈ [0, τ ∞ ], X rev τ∞ ∼ p Xτ ∞ .As τ ∞ → ∞, it can be seen from (2) that p Xτ ∞ converges to N (0, I d ). Thus the solution of the above SDE can be approximated by initializing X rev τ∞ ∼ N (0, I d ) instead. The DDPM sampler [HJA20] can be viewed as a discretization of this SDE [SSDK + 21]. DDIM-type deterministic samplers. On the other hand, the probability flow ODE (8) for the OU process reads asEQUATIONAgain, as τ ∞ → ∞, one may approximate the initialization with X rev τ∞ ∼ N (0, I d ). It is known that the popular DDIM sampler [SSDK + 21, SME20] is a discretization of this ODE [ZC22] . The ODE-based deterministic samplers allow more aggressive choice of discretization schedules, as well as fast ODE solvers [LZB + 22], enabling significantly accelerated sampling process compared to the SDE-based stochastic samplers.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: We discuss the scope of our work and provide ablation studies.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
posterior sampling via diffusion plug-and-play: 
 We are interested in solving (possibly nonlinear) inverse problems, where the aim is to infer an unknown imagex ⋆ ∈ R d from its measurements y ∈ R m , y = A(x ⋆ ) + ξ,where A : R d → R m is the measurement operator underneath the forward model, and ξ denotes measurement noise. We focus on the Bayesian setting where the prior information of x ⋆ is provided in the form of some prior distribution p ⋆ (•), i.e.,EQUATIONThe posterior distribution given measurements y is defined asEQUATIONHere, L(•; y) is the log-likelihood function of the measurements. Notwithstanding, our framework allows flexible choices of the forward model and the noise distributions. In addition, while this formulation is derived from probabilistic interpretations, it also subsumes the "reward-guided" or "loss-guided" setting [SZY + 23], where L can be viewed as a reward function or a negative loss function, both of which characterize preference over structural properties of x ⋆ .Assumption on the forward model. Throughout the paper, for simplicity, we make the following mild assumption on L, which is applicable to many applications of interest. Assumption 1. We assume L(• ; y) is differentiable almost everywhere, and sup x∈R d L(x; y) < ∞.Goal. Our goal is to sample x from the posterior distribution x ∼ p ⋆ (• | y) given estimates s(τ, x) (resp. ε(τ, x)) of the unconditional score functions s(τ, x) (resp. the noise function ε(τ, x)) in (3), assuming knowledge of the likelihood function L(•; y).We begin with an inspection on one of the most fundamental inverse problems: denoising under white Gaussian noise. As shall be elucidated shortly, the denoising diffusion samplers turn out to be an important building block in our algorithm for general inverse problems.Image denoising under white Gaussian noise. Suppose that we have access to a noisy version of x ⋆ ∼ p ⋆ contaminated by white Gaussian noise, given byEQUATIONwhere η > 0 is the noise intensity assumed to be known. Our goal is to sample from p ⋆ (• | x noisy ) given the score estimates s t (x) (resp. the noise estimates ε t (x)). We will develop our score-based denoising posterior sampler, termed DDS, with two variants, DDS-DDPM and DDS-DDIM, which can be viewed as analogues of the well-known DDPM and DDIM samplers in unconditional scorebased sampling respectively. Before proceeding, it is worth highlighting that the two variants will be derived from different forward diffusion processes, since we observe the resulting variants empirically lead to more competitive performance.A stochastic DDPM-type sampler via heat flow. We begin with a stochastic DDPM-type sampler for denoising, termed DDS-DDPM. We divide our development into the following steps.Step 1: introducing the heat flow. Let us introduce a heat flow with initial distribution p ⋆ , defined by the following SDE:EQUATIONwhere (B τ ) τ ≥0 is the standard d-dimensional Brownian motion. The solution of (13) is simplyEQUATION) Since B τ ∼ N (0, τ I d ), it readily follows that B η 2 (d)= ξ, which together with Y 0 ∼ p ⋆ yield the important observation that x noisy = x ⋆ + ξ can be viewed as an endpoint of the heat flow, in the sense thatx noisy = x ⋆ + ξ (d) = Y η 2 . 2)Step 2: reversing the heat flow. Following similar reasonings in Section 2, the next step boils down to reverse the heat flow (13). The time-reversal of the heat flow SDE ( 13) is (cf. (7)) given byEQUATIONwhere ( Bτ ) τ ≥0 is an independent copy of (B τ ) τ ≥0 . As introduced earlier, the virtue of the time-reversed SDE (15) is that it produces a process Y rev τ with the same path distribution as Y τ , i.e., (Y rev τ ) τ ∈[0,η 2 ] (d) = (Y τ ) τ ∈[0,η 2 ] .In particular, the joint distribution of (Y rev 0 , Y rev η 2 ) is the same as that of (Y 0 , Y η 2 )(d) = (x ⋆ , x noisy ). This implies that the conditional distribution p ⋆ (• | x noisy ) is the same as p Y rev 0 (• | Y rev η 2 = x noisy ). Surprisingly, the latter admits a simple interpretation:p Y rev 0 (• | Y rev η 2 = x noisy )is the distribution of Y rev 0 when we initialize (15) with Y rev η 2 = x noisy ! Therefore, sampling the posterior p ⋆ (• | x noisy ) amounts to solving the following simple SDE:EQUATION3)Step 3: connecting the score functions. It is now immediate to arrive at our proposed stochastic sampler DDS-DDPM by discretization of this SDE (16), which requires knowledge of the score functions ∇ log p Yτ (•). A key observation is that they can in fact be computed from the score function s(τ, x) (cf. (3)), due to the following lemma, whose proof is provided in Appendix D.1.Lemma 1 (Score function of Y τ ). For τ ≥ 0, we have∇ log p Yτ (x) = 1 √ 1 + τ s 1 2 log(1 + τ ), x √ 1 + τ .The resulting sampler, DDS-DDPM, is summarized in Algorithm 2 (deferred in the appendix) using a discretization procedure with an exponential integrator [ZC22] .A deterministic DDIM-type sampler via OU process. We next develop a deterministic DDIM-type sampler for denoising, termed DDS-DDIM.Step 1: introducing a posterior-initialized OU process. To sample from the posterior distribution p ⋆ (•|x noisy ), we first introduce a random variable w which has (unconditional) distributionEQUATIONin the same form of the desired posterior distribution p ⋆ (•|x noisy ). Here, since the noisy observation x noisy is given, we regard it as fixed.We then further introduce z = w -x noisy , which is a "centered" version of w, whose distribution isp z (x) := p w (x + x noisy ) = p ⋆ (x ⋆ = x + x noisy | x ⋆ + ξ = x noisy ).The OU process with initial distribution p z is defined by the SDE:EQUATIONwhere B τ is the standard d-dimensional Brownian motion. As in (2), the marginal distribution of Z τ is given byEQUATION2)Step 2: reversing the OU process. Following similar reasonings in Section 2, reversing the OU process (18) will enable us to generate samples z ∼ p z . Then we can set w = z + x noisy , which, by definition, has distribution p w defined in (17), and is a sample from the desired posterior distribution p ⋆ (•|x noisy ). We are thus led to solve the time-reversed probability flow ODE (cf. (8)) of (18), given bydZ rev τ rev = Z rev τ rev +∇ log p Z τ rev (τ rev , Z rev τ rev ) dτ, τ ∈ [0, τ ∞ ], Z rev τ∞ ∼ N (0, I d ), τ rev = τ ∞ -τ. (20) 3)Step 3: connecting the score functions. We are now one step away from our proposed deterministic sampler DDS-DDIM, which is derived by discretization of the ODE (20). We need to know the score functions ∇ log p Zτ (•), which again can be computed from the score function s(τ, x) (cf. (3)), as documented by the following lemma, whose proof is provided in Appendix D.2.Lemma 2 (Score function of Z τ ). For τ ≥ 0, we haveEQUATIONwhereEQUATION)After plugging this into (20) and solving the ODE for Z rev τ , we see that Z rev 0 + x noisy is the desired sample from the posterior distribution p ⋆ (•|x noisy ), as argued before. Numerically, the ODE (20) is solved by discretization with an exponential integrator [ZC22] , resulting in the sampler DDS-DDIM as summarized in Algorithm 3 (deferred in the appendix).Algorithm 1 Diffusion Plug-and-Play (DPnP) Input: Measurements y ∈ R m , log-likelihood function L(•; y) of the forward model, score estimates s, annealing schedule (η k ) 0≤k≤K . Initialization: Sample x 0 ∼ N (0, η0 4 I d ) Alternating sampling: for k = 0, 1, 2, . . . , K -1 do (1) Proximal consistency sampler:Sample x k+ 1 2 ∝ exp L(• ; y) -1 2η 2 k ∥ • -x k ∥ 2 using sub- routine PCS( x k , y, L, η k ) (Alg. 4).(2) Denoising diffusion sampler:Sample x k+1 ∼ exp log p ⋆ (x) -1 2η 2 k ∥x -x k+ 1 2 ∥ 2 using subroutine DDS-DDPM( x k+ 1 2 , s, η k ) (Alg. 2) or DDS-DDIM( x k+ 1 2 , s, η k ) (Alg. 3). Output: x K .Now we turn to the general setting where the measurement operator A is arbitrary. From the factorization of posterior distribution in (11), one intuitively understands that a posterior sampler must obey two constraints simultaneously: (i) the data prior constraint, corresponding to the first factor p ⋆ (x), which imposes that the posterior sampler should be less likely to sample at those points where p ⋆ (x) is small; (ii) the measurement consistency constraint, corresponding to the second factor e L(x;y) , which imposes that A(x) ≈ y.Diffusion plug-and-play (DPnP). We will apply the idea of alternatively enforcing these two constraints from a sampling perspective in the same spirit of [VDC19, LST21, BB23]. Our algorithm, dubbed diffusion plug-and-play (DPnP), alternates between two samplers, the denoising diffusion sampler (DDS) and the proximal consistency sampler (PCS), which can be viewed as the substitutes for the proximal operator and the gradient step respectively. Given the iterate x k and the annealing parameter η k at the k-th iteration, DPnP proceeds with the following two steps:(i) Proximal consistency sampler to enforce the measurement consistency constraint. DPnP draws a sample x k+ 1 2 from the distribution proportional to exp L(x ; y) -1 2η 2 k ∥x -x k ∥ 2to promote the image to be consistent with the measurements. This step, which we denote as the proximal consistency sampler, can be achieved by small modifications of standard algorithms such as Metropolis-Adjusted Langevin Algorithm (MALA) [RR98] given in Algorithm 4 (deferred in the appendix).(ii) Denoising diffusion sampler to enforce the data prior constraint. DPnP next draws a sample x k+1 from the distribution proportional toEQUATIONto promote the image to be consistent with the prior, where w ∼ N (0, I d ). The last step, which follows from the Bayes' rule, makes it clear that this step can be precisely achieved by the denoising diffusion sampler (developed in Section 3.1) using solely the unconditional score function.Combining both steps lead to the proposed DPnP method described in Algorithm 1. Some comments about the proposed DPnP method are in order.• The proximal consistency sampler PCS can be viewed as a "soft" version of the proximal point method [Dru17] . This can be seen from a first-order approximation: the maximum likelihood of the distribution exp L(•; y) -1 2η 2 k ∥ • -x k ∥ 2 is attained at the point x ′ ∈ R d satisfying ∇ x ′ L(x ′ ; y) - 1 η 2 k (x ′ -x k ) = 0, =⇒ x ′ = x k + η 2 k ∇ x ′ L(x ′ ; y) ≈ x k + η 2 k ∇ x k L( x k ; y),Therefore, the proximal consistency sampler draws random samples "concentrated" around x ′ , which approximates the implicit proximal point update, akin to a gradient step at x k . • On the other end, the denoising posterior sampler DDS can be regarded as a "soft" version of the proximal operator. In particular, when p ⋆ is supported on a low-dimensional manifold M, it forces x to reside in M, like the proximal map. To see this, note that denoising posterior distribution vanishes outside M by (23).• The proximal consistency sampler PCS admits a simple form when the forward model A is linear, i.e. A(x) = Ax for some matrix A ∈ R m×d , and the measurement noise ξ ∼ N (0, Σ) is Gaussian. In this situation, the proximal consistency sampler PCS can be implemented directly byx k+ 1 2 = PCS( x k , y, L, η k ) = x k + Σ 1/2 k w k , w k ∼ N (0, I d ),wherex k = A ⊤ Σ -1 A + 1 η 2 k I d -1 A ⊤ Σ -1 y + 1 η 2 k x k , and Σ k = A ⊤ Σ -1 A + 1 η 2 k I d -1.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer • If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: While our provably robust algorithm may potentially improve the fieldity of diffusion-based methods in solving inverse problems, it awaits future research before anything special can be said about potential societal impacts. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
theoretical analysis: 
 In this section, we establish both asymptotic and non-asymptotic performance guarantees of DPnP.Asymptotic consistency. We begin with the asymptotic consistency of DPnP in the theorem below.Theorem 1 (Asymptotic consistency of DPnP). Assume the score function estimate s(τ, •) is accurate, i.e., s(τ, x) = s(τ, x), and assume the ODE/SDEs in DDS and PCS are solved exactly. Let (ε l ) l≥0 be a decreasing sequence of positive numbers satisfying lim l→∞ ε l = 0, and (k l ) l≥0 be an increasing sequence of integers with k 0 = 0. Set the annealing scheduleη k = ε l , for k l-1 ≤ k < k l , l = 1, 2, • • • Let min l ′ =1,2,••• |k l ′ -k l ′ -1 | → ∞, the output x k l of DPnP converges in distribution to the posterior distribution p ⋆ (•|y) for l → ∞.In words, Theorem 1 establishes the asymptotic consistency of DPnP under fairly mild assumptions on the forward model (cf. Assumption 1): as long as the sampled distributions of DDS and PCS are exact, then running DPnP with a slowly diminishing annealing schedule of {η k } will output samples approaching the desired posterior distribution p ⋆ (•|y) when the number of iterations l goes to infinity.Non-asymptotic error analysis. We now step away from the idealized setting when the sampled distributions of DDS and PCS are exact. In practice, there are many sources of errors that can influence the sampled distributions of DDS and PCS, e.g., the discretization error arising from numerically solving ODE/SDE, and the score estimation error. In effect, these non-idealities will make PCS and DDS inexact. That is, the distribution they generate will slightly deviate from the distribution they ought to sample from. In this paper, we model such deviations by the total variation distance from the distribution generated by PCS (resp. DDS) to the ideal distribution proportional to exp(L(x; y) -1 2η 2 k ∥x -x k ∥ 2 ) (resp. p ⋆ (x ⋆ = x|x ⋆ + η k ε = x k+ 1 2 ))uniformly over all iterations. Analyzing these errors is out of the scope of this paper, and we point the interested readers to parallel lines of works, e.g., [LWCC23, MV19, CLA + 21], among many others. In our analysis, we will assume a black-box bound for the total variation errors of PCS and DDS, which can be combined with existing analyses of the respective samplers to bound the iteration complexity of DPnP. Theorem 2 (Non-asymptotic robustness of DPnP). With the notation in DPnP (Algorithm 1), set η k ≡ η > 0. Under Assumption 1, there exists λ := λ(p ⋆ , L, η) ∈ (0, 1), such that the following holds. Define a stationary distribution π η by π η (x) ∝ p ⋆ (x)q η (x), where q η is defined byEQUATIONwhere * denotes convolution. If PCS has error at most ε PCS in total variation and DDS has error at most ε DDS in total variation per iteration, then for any accuracy goalEQUATION)Before interpreting Theorem 2, we observe that q 0 (x) = e L(x;y) , thus π 0 (x) ∝ p ⋆ (x)e L(x;y) coincides with the desired posterior distribution p ⋆ (•|y). Thus Theorem 2 tells us that, assuming a constant annealing schedule η k = η, the output of DPnP converges in total variation to the distribution π η , which is a distorted version of the desired posterior distribution up to level η, with sufficiently many iterations. A few remarks are in order.Non-diminishing η. It can be seen from Theorem 2 that even with a nonzero η, DPnP already enforces the data prior strictly. On the other hand, the measurement consistency is distorted by an order of η. This is usually tolerable, since the measurements are themselves contaminated by Ground truth noise, thus when η is smaller than the noise level, the distortion would be tolerable. In practice, it is beneficial to choose an annealing schedule of {η k }, which will be elaborated in Section 5.Provable robustness. Theorem 2 indicates the performance of DPnP degenerates gracefully in the presence of sampling errors. To the best of our knowledge, this is the first provably consistent and robust posterior sampling method for nonlinear inverse problems using score-based diffusion priors. 
numerical experiments: 
 We provide preliminary numerical evidence to corroborate the promise of DPnP in solving both linear and nonlinear image reconstruction tasks. We denote DPnP with the subroutines DDS-DDPM and DDS-DDIM as DPnP-DDPM and DPnP-DDIM respectively.  
discussion: 
 This paper sets forth a rigorous and versatile algorithmic framework called DPnP for solving nonlinear inverse problems via posterior sampling, using image priors prescribed by score-based diffusion models with general forward models. DPnP alternates between two sampling steps implemented by DDS and PCS, to promote consistency with the data prior constraint and the measurement constraint respectively. We provide both asymptotic and non-asymptotic convergence guarantees, establishing DPnP as the first provably consistent and robust score-based diffusion posterior sampling method for general nonlinear inverse problems.Algorithmic unrolling and plug-and-play image reconstruction. Composite optimization algorithms, which aim to minimize the sum of a measurement fidelity term and a regularization term promoting desirable solution structures, have been the backbone of inverse problem solvers. To unleash the power of deep learning, [GL10] advocates the perspective of algorithmic unrolling, which turns an iterative algorithm into concatenations of linear and nonlinear layers like in a neural network.[VBW13] recognized that the proximal mapping step in many composite optimization algorithms can be regarded as a denoiser or denoising operator with respect to the given prior, and proposed to "plug in" alternative denoisers, in particular state-of-the-art deep learning denoisers, leading to a class of popular algorithms known as plug-and-play methods [BCSB18] ; see [MLE21] for a review.Regularization by denoising and score matching. [Vin11] pointed out a connection between score matching and image denoising, which is a consequence of the Tweedie's formula [Efr11] . The regularization by denoising (RED) framework [REM17] follows the plug-and-play framework to minimize a regularized objective function, where the regularizer is defined based on the plug-in image denoiser; [RS18] later clarified that the RED framework can be interpreted as score matching by denoising using the Tweedie's formula. [KVE21] developed a stochastic image denoiser for posterior sampling of image denoising using annealed Langevin dynamics. [FBS24] provided a framework to learn exact proximal operators for inverse problems.Plug-and-play posterior sampling. Motivated by the need to characterize the uncertainty, tackling image reconstruction as posterior sampling from a Bayesian perspective is another important approach. Our method is inspired by the plug-and-play framework but takes on a sampling perspective, exploiting the connection between optimization and sampling [Wib18] . Along similar lines, [LBA + 22, BB23] proposed Bayesian counterparts of plug-and-play for posterior sampling, where they leveraged the connection to score matching for sampling from the image prior, but did not consider score-based diffusion models for the image prior, which is a key aspect of ours; see also [SWC + 23]. [CDC23] extended the split Gibbs sampler [VDC19] in the plug-and-play framework, and advocated the use of score-based diffusion models such as DDPM [HJA20] for image denoising based on heuristic observations. In contrast, we rigorously derive the denoising diffusion samplers from first principles, unraveling critical gaps from naïve applications of the generative samplers to denoising, and offer theoretical guarantees on the correctness of our approach.Score-based diffusion models as image priors. Several representative methods for solving inverse problems using score-based diffusion priors alternates between taking steps along the diffusion process and projecting onto the measurement constraint, e.g. The discrete-time forward and backward diffusion processes can be understood as two processes constructed in the following manner:1) a forward processx 0 → x 1 → • • • →x T that starts with samples from the target image distribution and diffuses into a noise distribution (e.g., standard Gaussians) by gradually injecting noise into the samples;2) a reverse processx rev T → x rev T -1 → • • • → x rev 0that starts from pure noise (e.g., standard Gaussians) and converts it into samples whose distribution is close to the target image distribution.Consider the forward Markov process in R d that starts with a sample from the data distribution p X , and adds noise over the trajectory according toEQUATIONEQUATIONwhere {w t } 1≤t≤T 's are independent standard Gaussian vectors, i.e., w t i.i.d.∼ N (0, I d ), and {β t ∈ (0, 1)} describes the noise-injection rates used in each step. Therefore, we can write x t equivalently asEQUATIONHere, (ᾱ t ) t=0,1,••• ,T is the schedule of diffusion given byEQUATIONClearly, it verifies that1 ≥ ᾱ0 > ᾱ1 > • • • > ᾱT > 0.As long as ᾱT is vanishing, it is easy to observe that the distribution of x T approaches N (0, I d ).Score functions. As will be seen, in order to sample from p ⋆ , it turns out to be sufficient to learn the score functions of p xt at each step of the forward process, defined asEQUATIONAs in the continuous time, the score function can be viewed as a MMSE estimate:EQUATIONComparing ( 27) and (2), it can be checked that the discrete-time diffusion process can be embedded into the continuous-time one via the time changet → 1 2 log 1 ᾱt , in the sense that x ⋆ t (d) = X 1 2 log 1 ᾱt .This establishes a correspondence between the continuous-time formulation and the discrete-time formulation.Within the discrete-time formulation, we assume the score function estimates are given as s t (•) : R d → R d , t = 1, . . . , T such that s t ≈ s ⋆ t in analogy with the continuous-time counterpart.With the discrete-time perspective established in Appendix B, we are now ready to present the detailed description and the implementation of our algorithms DDS-DDPM (Algorithm 2), DDS-DDIM (Algorithm 3), and PCS (Algorithm 4).Proof. The marginal distribution (14) of the heat flow can be written asEQUATIONComparing ( 2) and (31), it is not hard to check thatY τ (d) = √ 1 + τ X 1 2 log(1+τ ) .Algorithm 2 Denoising Diffusion Sampler (stochastic) DDS-DDPM(x , s, η)Input: noisy data x noisy ∈ R d , score estimates s := { s t (•) : R d → R d , t = 1, . . . , T } or noise estimates ε = { ε t (•) : R d → R d , t = 1, .. . , T }, and noise level η > 0.Scheduling: Compute the diffusion schedule (τ t ) 0≤t≤T ′ byEQUATIONT ′ := max t : 0 ≤ t ≤ T, ᾱt > 1 η 2 + 1.Initialization: Set x T ′ = x noisy . Diffusion: for t = T ′ , T ′ -1, . . . , 1 do x t-1 = x t -2( √ τ t - √ τ t-1 ) ε t + τ t -τ t-1 w t , w t ∼ N (0, I d ).whereε t := ε t ( √ ᾱt x t ) = - 1 √ 1 -ᾱt s t √ ᾱt x t .Output: x 0 .Algorithm 3 Denoising Diffusion Sampler (deterministic) DDS-DDIM(x noisy , s, η)Input: noisy data x noisy ∈ R d , score estimates s := { s t (•) : R d → R d , t = 1, . . . , T } or noise estimates ε = { ε t (•) : R d → R d , t = 1, .. . , T }, and noise level η > 0.Scheduling: Compute the diffusion schedule(ū t ) 0≤t≤T ′ by ūt = (η 2 + 1) ᾱt -1 η 2 + ᾱt -1 , 0 ≤ t ≤ T ′ ,whereT ′ := max t : 0 ≤ t ≤ T, ᾱt > 1 η 2 + 1.Initialization: Draw z T ′ ∼ N (0, I d ).Diffusion: for t = T ′ , T ′ -1, . . . , 1 doz t-1 = (η 2 -1)ū t-1 + 1 (η 2 -1)ū t + 1 z t + (η 2 -1)ū t-1 + 1 • h(η, ūt-1 ) -h(η, ūt ) ε t , where h(η, u) := -arctan η √ u -1 -1 , ε t := ε t √ ᾱt x noisy + η 2 √ ūt ᾱt z t (η 2 -1)ū t + 1 = - 1 √ 1 -ᾱt s t √ ᾱt x noisy + η 2 √ ūt ᾱt z t (η 2 -1)ū t + 1.Output: x noisy + z 0 .Denote θ = 1 2 log(1 + τ ) as a short-hand. We havep Yτ (x) = p √ 1+τ X θ (x) ∝ p X θ 1 √ 1 + τ x .Therefore it follows that∇ log p Yτ (x) = ∇ x log p X θ 1 √ 1 + τ x = 1 √ 1 + τ s θ, 1 √ 1 + τ x ,where we used the definition s(θ, •) = ∇ log p X θ (•). Plugging the definition θ = 1 2 log(1 + τ ) into the above equation yields the desired result. (1) One step of discretized Langevin: Set r = e -γ/η 2 , andz n+ 1 2 = rz n + (1 -r)x + η 2 (1 -r)∇ zn L(z n ; y) + η 1 -r 2 w n , w n ∼ N (0, I d ). This is equivalent to drawing z n+ 1 2 from a distribution with density Q(•; z n ), whereQ(z ′ ; z) = 1 (2π(1 -r 2 )) d/2 exp - z ′ -rz + (1 -r)x + η 2 (1 -r)∇ z L(z; y) 2 2(1 -r 2 ).(2) Metropolis adjustment: Computeq = exp L(z n+ 1 2 ; y) -1 2η 2 ∥z n+ 1 2 -x∥ 2 exp L(z n ; y) -1 2η 2 ∥z n -x∥ 2 • Q(z n ; z n+ 1 2 ) Q(z n+ 1 2 ; z n ) ,and setz n+1 = z n+ 1 2 , with probability min(1, q), z nwith probability 1 -min(1, q).Output: z N .Proof. We first compute the probability density function of z. Recall that z = w -x noisy , thus applying Bayes rule yieldsp z (x) = p w (x + x noisy ) = p ⋆ (x ⋆ = x + x noisy |x ⋆ + ξ = x noisy ) = p ⋆ (x + x noisy )p ξ (-x) p x ⋆ +ξ (x noisy ) ∝ p ⋆ (x + x noisy )p ξ (-x),where ξ ∼ N (0, η 2 I d ). It is straightforward to computeEQUATION)We proceed to compute the probability density function of Z τ . According to (19), it follows thatp Zτ (x) = p e -τ z * p √ 1-e -2τ ε (x) = p e -τ z (x ′ ) p √ 1-e -2τ ε (x -x ′ )dx ′ ∝ p z (e τ x ′ ) exp - 1 2(1 -e -2τ ) ∥x -x ′ ∥ 2 dx ′ ∝ p ⋆ (x noisy + e τ x ′ ) exp - 1 2η 2 ∥e τ x ′ ∥ 2 exp - 1 2(1 -e -2τ ) ∥x -x ′ ∥ 2 dx ′ , ∝ p ⋆ (x ′ ) exp - 1 2η 2 ∥x ′ -x noisy ∥ 2 exp - 1 2(1 -e -2τ )∥x -e -τ (x ′ -x noisy )∥ 2 dx ′ , (33) where * denotes convolution, the penultimate line follows from (32) and the last line follow from the change of variable x ′ → e -τ (x ′ -x noisy ). One may exercise some brute force to verify thatexp - 1 2η 2 ∥x ′ -x noisy ∥ 2 exp - 1 2(1 -e -2τ ) ∥x -e -τ (x ′ -x noisy )∥ 2 = exp - e 2τ ∥x∥ 2 2(η 2 + e 2τ -1) exp - 1 2(1 -e -2τ ) e -τ x noisy + e τ -τ η 2 x η 2 + e -1 -e -τ x ′ 2 ∝ exp - e 2τ ∥x∥ 2 2(η 2 + e 2τ -1) p √ 1-e -2τ ε e -τ x noisy + e τ -τ η 2 x η 2 + e 2τ -1 -e -τ x ′ ,where τ is as defined in (22). Plug this back into (33), we seep Zτ (x) ∝ exp - e 2τ ∥x∥ 2 2(η 2 + e 2τ -1) p ⋆ (x ′ )p √ 1-e -2τ ε e -τ x noisy + e τ -τ η 2 x η 2 + e 2τ -1 -e -τ x ′ dx ′ ∝ exp - e 2τ ∥x∥ 2 2(η 2 + e 2τ -1) p ⋆ (e τ x ′ )p √ 1-e -2τ ε e -τ x noisy + e τ -τ η 2 x η 2 + e 2τ -1 -x ′ dx ′ ∝ exp - e 2τ ∥x∥ 2 2(η 2 + e 2τ -1) p e -τ x0 * p √ 1-e -2τ ε e -τ x noisy + e τ -τ η 2 x η 2 + e 2τ -1 ∝ exp - e 2τ ∥x∥ 2 2(η 2 + e 2τ -1)p Xτ e -τ x noisy + e τ -τ η 2 x η 2 + e 2τ -1 , where the second line applies the change of variable x ′ → e τ x ′ in the integral, the penultimate line follows from p e -τ x0 (x ′ ) ∝ p ⋆ (e τ x ′ ) (since x 0 ∼ p ⋆ ), and the last line follows from X τ (d) = e -τ x 0 + √ 1 -e -2τ ε.Finally, from the above formula, we obtain∇ log p Zτ (x) = ∇ x - e 2τ ∥x∥ 2 2(η 2 + e 2τ -1) + ∇ x log p Xτ e -τ x noisy + e τ -τ η 2 x η 2 + e 2τ -1 = - e 2τ x η 2 + e 2τ -1 + e τ -τ η 2 η 2 + e 2τ -1 s τ , e -τ x noisy + e τ -τ η 2 x η 2 + e 2τ -1 ,where we used the definition s(τ , •) = ∇ log p Xτ (•).Consider a SDE of the form:dM τ = v(τ )M τ + f (τ, M τ ) dτ + βdB τ , τ ∈ [0, τ ∞ ], M 0 ∼ p M0 , where v : [0, τ ∞ ] → R, f : [0, τ ∞ ] × R d → R d are deterministic functions, and β > 0 is a constant. Given discretization time points 0 = τ 0 ≤ τ 1 ≤ • • • ≤ τ k ≤ τ ∞ , a naïve way to discretize the SDE is M τi+1 -M τi ≈ v(τ i )M τi + f (τ i , M τi ) (τ i+1 -τ i ) + β τ i+1 -τ i ε i , i = 0, 1, • • • , k -1,where ε i ∼ N (0, I d ) is a standard d-dimensional Gaussian random vector which is independent of M τi . Although this approach is straightforward, it has the drawback that the linear term v(τ )M τ is discretized rather crude. For example, for the OU process where v ≡ -1, f ≡ 0, β = 2, the SDE can be solved analytically as in (2), while the above approach still has a discretization error.A more accurate discretization, known to significantly improve the quality of score-based generative models, is given by the exponential integrator [ZC22] , which preserves the linear term and discretizes the SDE tod M τ = v(τ ) M τ + f (τ i , M τi ) dτ + βdB τ , τ ∈ [τ i , τ i+1 ], i = 0, 1, • • • , k, with initialization M 0 ∼ p M0 . On each time interval [τ i , τ i+1 ],this is simply a linear SDE, which can be explicitly solved byM τ (d) = e V (τ )-V (τi) M τi + τ τi e V (τ )-V (τ ) dτ f (τ i , M τi ) + β τ τi e 2(V (τ )-V (τ )) dτ 1/2 ε i ,where V is the antiderivative of v:V (τ ) = τ 0 v(τ )dτ .Taking τ = τ i+1 , we obtainEQUATIONwhich provides an iterative formula to compute M τi+1Plug the expression of ∇ log p Yτ in Lemma 1 into (16), and use the notation τ rev = η 2 -τ , we obtain, for τ ∈ [0, η 2 ], thatdY rev τ rev = 1 √ 1 + τ rev s 1 2 log(1 + τ rev ), Y rev τ rev √ 1 + τ rev dτ + d Bτ = - 1 √ τ rev ε cont 1 2 log(1 + τ rev ), Y rev τ rev √ 1 + τ rev dτ + d Bτ .Choosing discretization time points. To discretize this SDE, we first choose the discretization time points. Recalling (??), it is most reasonable to discretize at those time points 0 ≤ τ rev0 ≤ • • • ≤ τ rev T ′ ≤ η 2 which satisfy 1 2 log(1 + τ rev t ) = 1 2 log 1 ᾱt , 0 ≤ t ≤ T ′ .This solves toEQUATIONThe requirement that τ rev t ≤ η 2 translates to ᾱt ≥ 1 1+η 2 , which yields the following choice of T ′ :EQUATION)Applying the exponential integrator. Now we apply the exponential integrator to discretize the SDE on each time interval τ rev ∈ [τ t-1 , τ t ], t = 1, • • • , T ′ as follows:d Y rev τ rev = - 1 √ τ rev ε cont 1 2 log(1 + τ rev t ), Y rev τ rev t √ 1 + τ rev t dτ + d Bτ , = - 1 √ τ rev ε ⋆ t Y rev τ rev t √ 1 + τ rev t dτ + d Bτ = - 1 √ τ rev ε ⋆ t √ ᾱt Y rev τ rev t dτ + d Bτ .The SDE can be integrated directly on τ rev ∈ [τ t-1 , τ t ] (see also (34), with v ≡ 0), yieldingY rev τ rev t-1 = Y rev τ rev t -2( τ rev t -τ rev t-1 ) • ε ⋆ t √ ᾱt Y rev τ rev t + η 2 -τt-1 η 2 -τt d Bτ dτ (d) = Y rev τ rev t -2( τ rev t -τ rev t-1 ) • ε ⋆ t √ ᾱt Y rev τ rev t + τ rev t -τ rev t-1 w t ,whereEQUATIONwhich is exactly the update equation in Algorithm 2, except that ε ⋆ t is replaced by the noise estimate ε t .Plug in the expression of s Z in Lemma 2 into the probability flow ODE (20), we obtaindZ rev τ rev = η 2 -1 η 2 + e 2τ rev -1 Z rev τ rev dτ + e τ rev -τ (τ rev ) η 2 η 2 + e 2τ rev -1 s τ (τ rev ), e -τ (τ rev ) noisy + e τ rev -τ (τ rev ) η 2 x η 2 + e 2τ rev -1 dτ = η 2 -1 η 2 + e 2τ rev -1 Z rev τ rev dτ - e 2τ rev e 2τ rev -1 ε cont τ (τ rev), e -τ (τ rev ) x noisy + e τ rev -τ (τ rev ) η 2 x η 2 + e 2τ rev -1 dτ,where the second line used the definition (22).Choosing discretization time points. Similar to the derivation in Appendix E.2, we discretize at time pointsEQUATION)which solves toEQUATION)To make this well-defined, we requireη 2 + ᾱt -1 (η 2 + 1) ᾱt -1 > 0,which is equivalent toᾱt > 1 1 + η 2 .This leads to the same choice of T ′ as in (36). We also setτ ∞ = τ rev T ′ .It is convenient to introduce a notation for the corresponding discrete schedule of τ rev t , denoted byūt = e -2τ rev t = (η 2 + 1) ᾱt -1 η 2 + ᾱt -1 , t = 0, 1, • • • , T ′ .Applying the exponential integrator. Now we apply the exponential integrator, which discretizes the ODE on each time intervalτ rev ∈ [τ t-1 , τ t ], t = 1, • • • , T ′ , as d Z rev τ rev = η 2 -1 η 2 + e 2τ rev -1 Z rev τ rev dτ - e 2τ rev e 2τ rev -1 ε cont τ (τ rev t ), e -τ (τ rev t ) x noisy + e τ rev t -τ (τ rev t ) η 2 Z rev τ rev η 2 + e 2τ rev t -1 dτ = η 2 -1 η 2 + e 2τ rev -1 Z rev τ rev dτ - e 2τ rev e 2τ rev -1 ε ⋆ t √ ᾱt x noisy + e τ rev t √ ᾱt η 2 Z rev τ rev η 2 + e 2τ rev t -1 dτ, = η 2 -1 η 2 + e 2τ rev -1 Z rev τ rev dτ - e 2τ rev e 2τ rev -1 ε ⋆ t √ ᾱt x noisy + √ ūt √ ᾱt η 2 Z rev τ rev (η 2 -1)ū t + 1 dτ,where the second line follows from (38), and the last line follows from dividing both the denominator and the numerator in the fraction inside ε t by e 2τ rev t . This is a first-order linear ODE on τ rev ∈ [τ t-1 , τ t ], which can be solved explicitly (cf. (34)) byZ rev τ rev = (η 2 -1)e -2τ rev + 1 (η 2 -1)ū t + 1 Z rev τ rev t + (η 2 -1)e -2τ rev + 1 • h(η, e -2τ rev ) -h(η, ūt ) • ε ⋆ t √ ᾱt x noisy + √ ūt √ ᾱt η 2 Z rev τ rev (η 2 -1)ū t + 1 , for τ rev ∈ [τ t-1 , τ t ], where h(η, u) := -arctan η √ u -1 -1 .Plug in τ rev = τ t-1 in the above solution, and set z t = Z rev τ rev t , we obtainz t-1 = (η 2 -1)ū t-1 + 1 (η 2 -1)ū t + 1 z t (40) + (η 2 -1)ū t-1 + 1 • h(η, ūt-1 ) -h(η, ūt ) • ε ⋆ t √ ᾱt x noisy + √ ūt √ ᾱt 2 z t (η 2 -1)ū t + 1 .The initialization, which should ideally be z T ′ = Z rev τ∞ ∼ p Zτ ∞ , is approximated by z T ′ ∼ N (0, I d ). This is exactly the update equation and the initialization in Algorithm 3, except that ε ⋆ t is replaced by the noise estimate ε t .We first note that the Metropolis-adjustment step in PCS (cf. Algorithm 4) is standard following the classical form of MALA [RR98] . Therefore, we focus on explaining the Langevin step. Recall the continuous-time Langevin dynamics for sampling from the distribution exp(L(•; y) -1 2η 2 ∥ • -x∥ 2 ):EQUATION)The classical form of MALA, as in [RR98] , performs one step of a straightforward discretization of (41) as the Langevin step, as follows:z n+ 1 2 ≈ z n -γ∇ zn L(z n ; y) + γ η 2 (z n -x) + 2γw n , w n ∼ N (0, I d ).In our setting, due to the presence of the linear drift term 1 η 2 (Z τ -x), which can be quite large when η is small, we apply the exponential integrator instead. Set the discretization time points τ n = nγ, the exponential integrator reads asdZ τ = -∇ Znγ L(Z nγ ; y)dτ + 1 η 2 (Z τ -x)dτ + √ 2dB τ , nγ ≤ τ ≤ (n + 1)γ.Solve this linear SDE on nγ ≤ τ ≤ (n + 1)γ directly (see also (34)) to obtainZ (n+1)γ (d) = rZ nγ + (1 -r)x + η 2 (1 -r)∇ Znγ L(Z nγ ; y) + η 1 -r 2 w n , w n ∼ N (0, I d ),where r := e -γ/η 2 . This is the same as the update equation for the Langevin step in PCS (cf. Algorithm 4).F.1 Proof of Theorem 1 Proof. We first collect the asymptotic correctness of our subroutines PCS and DDS in the following two lemmas. The correctness of PCS is actually well-known, see e.g., [Tie94, Corollary 2].Lemma 3 (Correctness of PCS). Under Assumption 1, with notation in Algorithm 4, in the continuous-time limit: γ → 0, N → ∞, the algorithm PCS outputs samples with distribution ∝ exp(L(•; y)+ 1 2η ∥ • -x∥ 2 ).The next lemma guarantees the correctness of DDS with exact unconditional score functions.Lemma 4 (Correctness of DDS). Assume the score function estimation s t is accurate, i.e. s t = s ⋆ t . In the continuous-time limit: T → ∞, ᾱT → 0, ᾱt-1 ᾱt → 1, unifomly in t, both DDS-DDIM and DDS-DDPM output x obeying the posterior distribution p⋆ (x ⋆ = x | x ⋆ + ηε = x noisy ), ε ∼ N (0, I d ).The proof of Theorem 1 is based on two lemmas on the one-step transition kernel of DPnP and the asymptotic behavior of the transition kernel, which we will present soon. First, we set up some notations. Denotep η (x) := p x ⋆ ∼p ⋆ ,ε∼N (0,I d ) (x ⋆ + ηε = x) = 1 (2π) d/2 η d p ⋆ (z)e -1 2η 2 ∥x-z∥ 2 dz.From the first equality, it is clear that p η → p ⋆ when η → 0 + . We will also use the notation q η defined in (24), which we recall here:q η (x) := 1 (2π) d/2 η d e L(z;y)-1 2η 2 ∥x-z∥ 2 dz.In virtue of the Assumption 1, we know that q η is finite for all x ∈ R d .For convenience, we introduce a notation for application of transition kernels. For a probability distribution p(x) and a probability transition kernel K(x, x ′ ), denote by p • K the probability distribution given byp • K(x ′ ) = p(x)K(x, x ′ )dx.The first lemma characterizes the one-step behavior of DPnP in terms of Markov transition kernels.Lemma 5. Under the settings of Lemma 4 and Lemma 3, the one-step transition kernel of DPnP with η k = η is given by:K DPnP,η (x, x ′ ) = q 0 (z) p η (z) e -1 2η 2 ∥z-x∥ 2 -1 2η 2 ∥z-x ′ ∥ 2 dz p ⋆ (x ′ ) q η (x) .In other words, if x k has distribution p x k , then the distribution of x k+1 isp x k+1 (x ′ ) = p x k • K DPnP,η (x) = p x k (x)K DPnP,η (x, x ′ )dx.The proof is postponed to Appendix F.3. The next lemma analyzes the ergodic properties of the Markov chain with transition kernel K DPnP,η . These properties are known [BB23] but scattered in different literatures, so we will provide a brief proof to be self-contained.Lemma 6. The Markov transition kernel K DPnP,η has the following properties:(i) (Stationary distribution.) Let π η be the probability distribution defined byπ η (x) = c η p ⋆ (x)q η (x),where c η > 0 is the normalization constant such that π η (x)dx = 1. Then K DPnP,η is reversible with stationary distribution π η .(ii) (Convergence.) For any initial distribution p, the distribution of the Markov chain with kernel K DPnP,η converges to π η :EQUATIONwhereK (n) DPnP,η is the n-step transition kernel of K DPnP,η .The proof is postponed to Appendix F.4. We now show how to prove Theorem 1 with the above two lemmas. With the annealing schedule in Theorem 1, between steps k l-1 ≤ k < k l , which consist of consecutive (k l -k l-1 ) steps, the transition kernel of one-step of DPnP is K DPnP,ε l . As (k l -k l-1 ) → ∞, Lemma 6 implies thatTV(p x k l , π ε l ) = TV(p x k l-1 • K (k l -k l-1 ) DPnP,ε l , π ε l ) → 0.Under the assumption in Theorem 1 that ε l → 0, we let l → ∞ to seelim l→∞ π ε l = c 0 p ⋆ (•)e L(•;y) = p ⋆ (•|y), thus p x k l → p ⋆ (•|y), as claimed.Proof. For DDS-DDPM, we note that under the continuous-time limit in Lemma 4, the discretization time points given by (35) verifyτ rev 0 = 0, sup 0≤t≤T ′ -1 |τ rev t -τ rev t+1 | → 0, τ rev T ′ → 1 1 + η 2 -1 -1 = η 2 , T ′ → ∞.Therefore, these discretization time points 0 = τ rev 0 ≤ • • • ≤ τ rev T ′ ≤ η 2 form a partition of [0, η 2 ], which becomes infinitely fine in the continuous-time limit. Thus the discretized integrator (37) converges to the solution of the SDE (16), which, as we have already argued in Appendix D, produces samples obeying the denoising posterior distribution p ⋆ (•|x noisy ), as claimed.The proof for DDS-DDPM follows similarly, by observing that the discretization time points in (35) form an infinitely fine partition of [0, ∞) in the continuous-time limit.Proof. The proof is based on computing the transition kernel of the two subroutines. We claim that (i) Sampling with probability density proportional to exp(L(•; y) -1 2η 2 ∥ • -x∥ 2 ) is equivalent to applying the following Markov transition kernelK PCS,η (x, x ′ ) = 1 q η (x) e L(x ′ ;y)-1 2η 2 ∥x ′ -x∥ 2 . (ii) Sampling with probability p ⋆ (x ⋆ | x ⋆ + ηε = x), where ε ∼ N (0, I d ), is equivalent to applying the following Markov transition kernel:K DDS,η (x, x ′ ) = 1 p η (x) p ⋆ (x ′ )e -1 2η 2 ∥x ′ -x∥ 2 .It is then clear thatK DPnP,η (x, x ′ ) = K PCS,η (x, z)K DDS,η (z, x ′ )dz = q 0 (z) p η (z) e -1 2η 2 ∥z-x∥ 2 -1 2η 2 ∥z-x ′ ∥ 2 dz p ⋆ (x ′ ) q η (x) ,as desired. We now prove the above two claims. For (i), note that by (23), we know K DDS,η (x, •) ∝p ⋆ (•)e -1 2η 2 ∥•-x∥ 2. Thus it suffices to compute the normalization constant, which isp ⋆ (x ′ )e -1 2η 2 ∥x ′ -x∥ 2 dx ′ = p η (x),by the definition of p η . ThereforeK DDS,η (x, x ′ ) = 1 p η (x) p ⋆ (x ′ )e -1 2η 2 ∥x ′ -x∥ 2, as claimed. The proof of (ii) follows similarly.Proof. We first introduce a fundamental lemma [PW24] , which provides a simple method to bound the total variation between two distributions.Lemma 7 (Data-processing inequality). Let p, q be two probability distributions, and K be a probability transition kernel. ThenTV(p • K, q • K) ≤ TV(p, q).We now prove the two items in Lemma 6 separately.Proof of (i). We first show that π η is well-defined, i.e., p ⋆ (x)q n η(x)dx < ∞. This can be seen from Assumption 1, which implies q η (x) ≲ eTo show that K DPnP,η is reversible with stationary distribution π η , it suffices to verifyπ η (x)K DPnP,η (x, x ′ ) = π η (x ′ )K DPnP,η (x ′ , x), ∀x, x ′ ∈ R d .However, it is easily checked that both sides are equal toc η q 0 (z) p η (z) e -1 2η 2 ∥z-x∥ 2 -1 2η 2 ∥z-x ′ ∥ 2 dz p ⋆ (x ′ )p ⋆ (x).Proof of (ii). We define an auxiliary Markov transition kernelK aux,η = K DDS,η • K PCS,η . More explicitly, K aux,η (x, x ′ ) = K DDS,η (x, z)K PCS,η (z, x ′ )dz = p ⋆ (z) q η (z) e -1 2η 2 ∥z-x∥ 2 -1 2η 2 ∥z-x ′ ∥ 2 dz e L(x ′ ;y) p η (x) . (43) It is easy to see that p • K (n) DPnP,η = p • K PCS,η • K (n-1) aux,η • K DDS,η .(44) Thus we are led to investigate the ergodic properties of K aux,η . Similar to the proof of item (i) above, it is not hard to show that K aux,η is reversible with respect to the stationary distribution µ η (x) := c η p η (x)q 0 (x) = c η p η (x)e L(x;y) .π η = µ η • K DDS,η . (45) It is apparent that µ(x ′ ) > 0 and K aux,η (x, x ′ )/µ η (x ′ ) > 0 for all x, x ′ ∈ R d . By [Tie94, Corollary 1], such a Markov transition kernel obeys, for any probability distribution q, thatTV(q • K (n) aux,η , µ η ) → 0, n → ∞.In view of ( 44) and (45), we set q = p • K PCS,η and invoke the data-processing inequality to obtainTV(p • K (n) DPnP,η , π η ) = TV(q • K (n-1) aux,η • K DDS,η , µ η • K DDS,η ) ≤ TV(q • K (n-1) aux,η , µ η ) → 0,as n → ∞. This completes the proof.Proof. Denote by KPCS,η and KDDS,η and the transition kernels for PCS and for DDS, respectively. Note that these may deviate from the transition kernels K PCS,η and K DDS,η defined for the idealized asymptotic setting in Appendix F. We haveTV(p x N , π η ) = TV(p x N -1 2 • KDDS,η , π η ) ≤ TV(p x N -1 2 • K DDS,η , π η ) + TV(p x N -1 2 • K DDS,η , p x N -1 2 • KDDS,η ) ≤ TV(p x N -1 2 • K DDS,η , π η ) + ε DDS ,where the second line is triangle inequality, and the third line follows from the assumption in Theorem 2 that DDS has error at most ε DDS in total variation, by taking the input of DDS to be x N -1 2 . Similarly, from p x N -1 2 = p x N -1 • KPCS,η and the assumption that PCS has error at most ε PCS in total variation, we can showTV(p x N -1 2 •K DDS,η , π η ) ≤ TV(p x N -1 •K PCS,η •K DDS,η , π η )+ε PCS = TV(p x N -1 •K DPnP,η , π η )+ε PCS .The above two inequalities together implyTV(p x N , π η ) ≤ TV(p x N -1 • K DPnP,η , π η ) + ε DDS + ε PCS .Iterating this process, we obtainEQUATIONIt remains to boundTV(p x1 • K (N -1)DPnP,η , π η ). For this, we need the following two lemmas.Lemma 8 (Comparing TV and χ 2 -divergence, [PW24] ). For any two distributions p, q, we haveTV(p, q) ≤ χ 2 (p ∥ q).Lemma 9 (χ 2 -contractivity of K DPnP,η ). There exists some λ := λ(p ⋆ , L, η) ∈ (0, 1), such that for any probability distribution p(x), we haveχ 2 (p • K (N ) DPnP,η ∥ π η ) ≤ λ 2N χ 2 (p ∥ π η ).A form of Lemma 9 is well-known for Markov chains with countable state spaces, but relatively few sources provide a complete proof for the abstract setting we consider here with continuous state space. For sake of completeness, we prove Lemma 9 in Appendix F.6.Combining the above two lemmas, we obtainTV(p x1 • K (N -1) DPnP,η , π η ) ≤ χ 2 (p x1 • K (N -1) DPnP,η ∥ π η ) ≤ λ N -1 χ 2 (p x1 ∥ π η ).Plug this into (46), we obtainTV(p x N , π η ) ≤ λ N -1 χ 2 (p x1 ∥ π η ) + (N -1)(ε DDS + ε PCS ).With N ≍ log(1/εacc) 1-λ such that λ N -1 ≤ exp -(N -1)(1 -λ) ≤ ε acc , the desired result readily follows.Proof. We need a few fundamental properties of reversible Markov chains, which are collected below.First we set up some notations. Define the Hilbert space L 2 (π) to be the space of square-integrable functions with respect to measure π, i.e., those functions f :R d → C such that ∥f ∥ L 2 (π) := |f (x)| 2 π(x)dx 1/2 < ∞.The first well-known property [SC97] offers a way to represent a reversible transition kernel as a self-adjoint operator (infinite-dimensional symmetric matrix).Lemma 10 (Self-adjoint representation of reversible Markov operator). Assume K(x, x ′ ) is a Markov transition kernel that is reversible with respect to the stationary distribution π(x). Then the integral operator K : L 2 (π) → L 2 (π) defined byKf (x) = K(x, x ′ )f (x ′ )dx ′is self-adjoint and compact. For any probability distribution p(x) such that p 2 (x) π(x) dx < ∞, we havep(x) • Kf (x)dx = p • K(x ′ )f (x ′ )dx ′ .Moreover, the eigenvalues of K are the same as those of K.The following theorem is a generalization of the classical Perron-Frobenius theory for finitedimensional transition matrix to strictly positive operators. The form we present here can be found in [Sch12, Theorem V.6.6]; see also [Bou23b, Theorem III.6 .7] for a more elementary treatment which can also be adapted to the form we need.Theorem 3 (Jentzsch). Let K(x, x ′ ) be a Markov transition kernel. If K(x, x ′ ) > 0 for any x, x ′ ∈ R d , then K has a unique stationary distribution π. Moreover, 1 is a simple eigenvalue of K, with π being the only left eigenfunction, and the constant function 1 being the only right eigenfunction. In addition, there exists λ ∈ (0, 1) such that any other eigenvalue of K has modulus no larger than λ.We are now ready to prove Lemma 9. We divide the proof into the following steps.Step 1: controlling the eigenvalues of K DPnP,η . Recall the auxiliary kernel K aux,η defined in (43). It is a standard result in linear algebra or function analysis [Bou23a] that K aux,η = K PCS,η • K DDS,η has same eigenvalues as K DPnP,η = K DDS,η • K PCS,η . From (43), it is easy to check K aux,η (x, x ′ ) > 0, thus Theorem 3 implies 1 is a simple eigenvalue of K DPnP,η . Moreover, there exists λ := λ(p ⋆ , L, η) ∈ (0, 1), such that any other eigenvalue of K aux,η has modulus no larger than λ.Since K DPnP,η has the same eigenvalues as K aux,η , and, by Lemma 10, the operator K DPnP,η also has the same eigenvalues as these two, we conclude that K DPnP,η is a self-adjoint compact operator on L 2 (π η ), of whom 1 is a simple eigenvalue. Moreover, any other eigenvalue of K DPnP,η has modulus no larger than λ.Step 2: establishing the contractivity of K DPnP,η in L 2 (π η ). It is easy to verify that the constant function 1, which takes value 1 for any x ∈ R d , is a eigenfunction of K DPnP,η associated to the simple eigenvalue 1, thus is the only (up to scaling) eigenfunction associated to that eigenvalue. It is also a unit-length eigenfunction, since ∥1∥ L 2 (πη) = ( 1 • π η (x)dx) 1/2 = 1. Therefore, the operator K DPnP,η -11 ⊤ is a self-adjoint operator whose eigenvalues have modulus no larger than λ, where 11 ⊤ is the orthogonal projection onto 1 in L 2 (π η ), defined by11 ⊤ f (x) ≡ f (x ′ )π η (x ′ )dx ′ , ∀x ∈ R d .Using the fact that K DPnP,η 11 ⊤ = 11 ⊤ K DPnP,η = 11 ⊤ , one may show(K DPnP,η -11 ⊤ ) N = K (N )DPnP,η -11 ⊤ by expanding the product, see e.g. [SC97] . Consequently, K N DPnP,η -11 ⊤ is a self-adjoint operator whose eigenvalues have modulus no larger than λ N , i.e.,EQUATIONwhere ∥ • ∥ L 2 (πη)→L 2 (πη) denotes the operator norm on L 2 (π η ).Step 3: bounding the inner product of p • K (N )DPnP,η -π η with any square-integrable function. Note that when χ 2 (p ∥ π η ) = ∞, the conclusion is trivially true. For the rest part of the proof, we assume χ 2 (p ∥ π η ) < ∞. Now, for any f ∈ L 2 (π η ), by applying Lemma 10 iteratively, we obtainp • K (N ) DPnP,η (x)f (x)dx = p(x ′ )K N DPnP,η f (x ′ )dx ′ = p(x ′ ) • (K N DPnP,η -11 ⊤ )f (x ′ )dx ′ + p(x ′ )11 ⊤ f (x ′ )dx ′ = p(x ′ ) • (K N DPnP,η -11 ⊤ )f (x ′ )dx ′ + f (x ′ )π η (x ′ )dx ′ ,where the last line follows from the definition of 11 ⊤ and p(x ′ )dx ′ = 1. Rearrange the terms to seeEQUATIONIn particular, taking p = π η yieldsEQUATIONSubstract (49) from (48), and then take absolute value, we obtain p • K (N ) DPnP,η (x) -π η (x) f (x)dx = p(x ′ ) -π η (x ′ ) • (K N DPnP,η -11 ⊤ )f (x ′ )dx ′ ≤ p(x ′ ) -π η (x ′ ) 2 π η (x) dx 1/2 • (K N DPnP,η -11 ⊤ )f (x ′ ) L 2 (πη)In the Algorithm 1, the initial guess x 0 is set to be a properly scaled Gaussian random vector. However, from Theorem 2 it can be infered that using a heuristic posterior sampler as the initializer could decrease χ 2 (p x1 ∥ π η ), hence potentially improve the convergence speed of DPnP. By using existing algorithms like DPS or LGD-MC as initializer, DPnP can improve upon the results of existing algorithms towards the correct posterior distribution efficiently and provably. In our experiments, we find it helpful to initialize DPnP with LGD-MC, which accelerates the algorithm significantly.We discuss the choice of the annealing schedule η k in DPnP (Algorithm 1). As seen in the theoretical analysis (Theorem 2), if we set all the η k ≡ η for some constant η > 0, then DPnP converges to a distribution π η , which can be regarded as a version of the posterior distribution p ⋆ (•|y) distorted by an order of O(η). The smaller η is, the more accurate the final distribution will be. On the other hand, it was also seen that in many cases, the spectral gap is Ω(η), hence the convergence time is O( 1 η ). Therefore, smaller η would make it take longer to converge. 4To strike a balance between the accuracy and the convergence rate, we find it empirically successful to adapt an gradually decreasing schedule for η k , similar to [BB23] . In the first few iterations, we set η k to be a large constant. After this initial phase, we decrease η k slowly, eventually to η N which is chosen to be a small constant. An example of such an annealing schedule isη 0 = η 1 = • • • = η K0 , η 0 > 0 is a large constant, η k = (η K /η 0 ) k-K 0 K-K 0 η 0 , K 0 < k ≤ K, η K > 0 a small constant,where K 0 < K is the length of the initial phase, which can be chosen as, e.g. K 0 = K/5. For all the numerical experiments, we set η 0 = 0.4, η N = 0.15, K 0 = 4, K = 20. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: No data or models is released by this work.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: All assets are credited properly.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.2η 2 ∥x-z∥ 2 dz ≲ 1, hence p ⋆ (x)q η (x)dx ≲ p ⋆ (x)dx = 1.https://github.com/DPS2022/diffusion-posterior-samplingStrictly speaking, while the number of iterations required to converge increases as η gets smaller, the computational complexity per iteration will decrease. However, in experiments, we observed that the latter effect is not strong enough to compensate for the increase in overall complexity caused by the former. 
 
 Appendix:  


Decision: In this paper, the authors present a plug-and-play algorithm based on score-based diffusion models as a signal prior for a variety of linear and nonlinear inverse problems.  In contrast to existing diffusion based approaches for inverse problems, the authors prove a robustness guarantee for their proposed method.  The proposed method also demonstrates superior performance in a variety of metrics on certain inverse problems.  The strengths of this paper are the novel theoretical guarantee of robustness in a diffusion based algorithm, and the flexibility of the proposed algorithm to handle a variety of diffusion models.  A weakness of the paper is that the experimental investigations could be more thorough.  Overall, the strengths outweigh the weaknesses, and the paper is recommended for acceptance.

Reviews: review: paper_summary: paper_summary: This paper focuses on developing a plug-and-play algorithm for using score-based diffusion models as an expressive prior for solving nonlinear inverse problems with general forward models. While going from current state $x_{k+1}$ to $x_{k}$, this paper makes two gradient updates: (1) from $x_{k+1}$ to $x_{k+\frac{1}{2}}$ using gradients from the measurement error and (2) from $x_{k+\frac{1}{2}}$ to $x_k$ using denoising diffusion score function. Theoretically, the authors prove robustness of the proposed algorithm in sampling from the posterior, and empirically, they show that the proposed algorithm outperforms one commonly used baseline DPS.
strength_weakness: strength_weakness: ### **Weaknesses and comments**

1. Line 171: "Assumption on forward model is applicable to *many* applications of interest" what are the applications?
2. Assumption 1: what are some examples of $\mathcal{L(\cdot, y)}$ that are differentiable almost everywhere and used in practice?
3. Step 1 in **A stochastic  DDPM-type sampler via heat flow**: How do you know the noise level $\eta$? This is typically unknown.
4. The paper is overloaded with lots of notations. To sample from the posterior, you need the measurement conditional score. By Bayes' theorem you can write this term as the combination of prior and likelihood. How do you get the unknown likelihood term exactly? Please explain the key idea without overloading with notations. 
5. Inverting $\Sigma$ is difficult in high-dimension posterior sampling problems. How do you get around this issue unless you are making further approximations like diagonal covariance or scalar? In that case, how do you sample from the posterior exactly?
6. One of the interesting parts of the theoretical results in this line of research is the characterization of the discretization error. But this seems to be out of scope of this paper according to the authors (line 299).
7. Theoretical results are valid provided you get a tight TV bound for both the consistency and diffusion samplers. Any ideas how to obtain these bounds?
8. Experimental results are compared with DPS which the authors claim to be state-of-the-art. Many of the recently developed methods [1,2,3,4] outperform DPS and there is no comparison with any stronger baselines. See citations below and references therein.
9. For evaluation which dataset is used. Is it the same subset used in the original paper and other follow-up papers? The results differ quite a lot if you pick some smaller subset.
10. The authors are encouraged to cite the published version of the papers where applicable (see for instance CCL+22).


References

1. Solving Linear Inverse Problems Provably via Posterior Sampling with Latent Diffusion Models.
2. Prompt-tuning latent diffusion models for inverse problems
3. Beyond First-Order Tweedie: Solving Inverse Problems using Latent Diffusion
4. Tweedie Moment Projected Diffusions for Inverse Problems
questions: questions: Please see the weakness section above.
limitations: limitations: Yes, the authors have addressed the limitations.
paper_summary: paper_summary: This paper introduces a diffusion-based sampling framework closely related to plug-and-play methods for solving general inverse problems. The technique alternates between two steps: calling a proximal consistency sampler that enforces data-fidelity, and regularization via a denoising diffusion sampler leveraging strong diffusion-based image priors.  Theoretical results demonstrate asymptotic consistency and robustness to sampling errors. Numerical experiments show promising reconstruction quality.
strength_weakness: strength_weakness: - The experimental evaluation is somewhat lacking. It would be interesting to see comparison with more contemporary solvers such as ReSample [1], which has improved robustness against sampling errors due to a posterior mean correction scheme. Moreover, in-depth ablation studies on the multiple hyperparameters of the algorithm are missing. Thus, it is unclear how much hyperparameter tuning is necessary.
- The proposed technique appears to have a very high compute cost (3000 NFEs). More discussion on the compute requirements and possible ways to accelerate the algorithm would be valuable. 

[1] Song, Bowen, et al. "Solving Inverse Problems with Latent Diffusion Models via Hard Data Consistency." The Twelfth International Conference on Learning Representations.
questions: questions: - Can the findings be extended to latent domain samplers? This would greatly improve the efficiency of the technique.
- How does the method compare to other samplers such as ReSample?
- How does performance scale with NFEs?
- I would recommend changing DDS to some other abbreviation to avoid confusion with the Decomposed Diffusion Sampling method [2].
- What does G denote in line 136?

[2] Chung, Hyungjin, Suhyeon Lee, and Jong Chul Ye. "Fast diffusion sampler for inverse problems by geometric decomposition." arXiv preprint arXiv:2303.05754 3.4 (2023).
limitations: limitations: Limitations are not clearly addressed in the paper.
paper_summary: paper_summary: This paper introduces a diffusion plug-and-play method (DPnP) that uses score-based diffusion models as expressive data priors for nonlinear inverse problems with general forward models. By combining a proximal consistency sampler and a denoising diffusion sampler, the method offers provably robust posterior sampling, with performance guarantees and demonstrated effectiveness across various tasks.
strength_weakness: strength_weakness: Although I appreciate the theoretical aspect of this paper, as mentioned in the abstract, “this paper develops an algorithmic framework for employing score-based diffusion models as an expressive data prior in nonlinear inverse problems.” However, there are already existing works on embedding denoising diffusion models into plug-and-play frameworks as data priors, such as [1-2]. It would be helpful if the author could clearly highlight any new insights within this paper to distinguish it from previous work; otherwise, the novelty of this paper may appear somewhat incremental. 

Minor one: There are several typos, e.g., in the abstract, the sentence "Score-based diffusion models, thanks to its impressive empirical success, have emerged as an appealing candidate of an expressive prior in image reconstruction." The correct pronoun should be "their" instead of "its" to match the plural subject "Score-based diffusion models."
questions: questions: I wonder if the asymptotic consistency and non-asymptotic error analysis of DPnP, as established in this paper, demonstrate convergence. If so, it is necessary to verify this through numerical experiments. 

Could you please explain why the result of DPnPDDPM shown in Table 1 is smooth, while the one shown in Table 5 is noisy with a lot of noticeable noise? Thank you.
limitations: limitations: As noted in the weaknesses section, the primary limitation is that diffusion-based PnP methods have already been proposed in [1-2] and applied to various inverse imaging problems.
paper_summary: paper_summary: The authors propose 'diffusion plug-and-play' (DPnP), a plug-and-play diffusion framework that alternatively calls what amounts to a consistency sampler based on the likelihood of the forward model, followed by (essentially) and unconditional diffusion step using the score function. The authors provide many subsequent theoretical results for this framework.
strength_weakness: strength_weakness: 1. This paper would benefit from improvements to the mathematical notation. Namely, it would be easier to read the math if scalar quantities were better differentiated from vector quantities (e.g., bold-faced vectors).
2. The formatting structure of the main paper is, at times, very granular with many lists. Personally, I appreciate this, but I feel that others may think that things could be ordered/structured better. This is a very minor issue, but worth keeping in mind.
3. I acknowledge that the primary contribution of this paper is the theoretical results associated with the DPnP framework, but I think that more robust experimental evaluation would have benefited this work. In particular, it would have been nice to see how the DPnP framework performs on more 'typical' linear inverse problems (e.g., inpainting, deblurring). This would have given a better sense of the overall performance of the approach, even if the focus is performance on non-linear inverse problems. I do not expect you to perform these experiments in the revision/for the rebuttal, I just wanted to note that the results would be more convincing if there were more of them.
4. It would be nice to discuss the connections between the proximal consistency sampler and projection-type approaches (e.g., that leveraged by Chung et al.'s 'Come-Closer-Diffuse-Faster) when the inverse problem is linear. In fact, I suspect that, for linear inverse problems, DPnP can be reduced to a simple projection step, followed by the reverse diffusion step. It would be worthwhile to discuss these connections, even if the discussion is relegated to an appendix.
5. In Algorithm 1, proximal consistency sampling is done before the denoising diffusion sampling. Why is this? The previously mentioned projection methods have the steps flipped, so I wonder why you have decided to structure the DPnP algorithm this way.
questions: questions: See Weaknesses.
limitations: limitations: The authors adequately discussed the limitations of their method.

Comments: title: Thank you for your response, and further clarification
Discussion with Authors
Thank you for your review. Have we addressed your concerns?
Thank you for your prompt reply, and clarification on our rebuttal
Further comments that may be helpful to improve the paper
Response to the author rebuttal
comment: Thank you for engaging with us! We are happy to hear that many of your concerns have been addressed successfully, and appreciate your detailed comments that provide us an opportunity to clarify further the remaining points.

Regarding your comments on Q4, we would like to clarify two points:
- Our approach does not involve approximating $\log p(y|x_t)$ as in the previous algorithms. Our $\mathcal L(x; y)$ is not an approximation to $\log p(y|x_t)$; it is simply the likelihood $\log p(y|x_0)$ (where $x_0$ is the ground truth signal), using the notation in DPS paper, which is assumed known in most of the previous works. Our approach deviates significantly from DPS; the reason we were able to bypass this approximation is that we directly tackle the whole posterior distribution $\propto p^\star(x_0)p(y|x_0) = p^\star(x_0) \exp(\mathcal{L}(x_0; y))$ using a split Gibbs sampler, which is made practical by developing proximal samplers for both $p^\star(x_0)$ and $\exp(\mathcal{L}(x_0; y))$ with diffusion and MALA.
- The measurement noise level $\sigma$ is assumed known or is a tunable parameter, which is consistent with most of the previous works (e.g. DDRM, DPS, LGD-MC, ReSample). As can be seen from the experimental results (including our rebuttal pdf), our algorithm demonstrates significant improvement for highly non-linear problems like phase retrieval, over recent works like LGD-MC and ReSample.

Regarding your comments on Q5, we would like to note that (i) it is a common choice as in many previous works (e.g. the references above) that $\Sigma$ is chosen as a scalar identity, in which case it is not necessary to perform a matrix inversion; (ii) The storage cost of prefactorization has been found to be managable (of $O(n)$ order) with memory efficient SVD in many practical inverse problems, e.g. denoising, inpainting, super resolution, deblurring, and colorization, cf. DDRM [KEES22]; (iii) If memory is really of concern, there is also the option to simply use MALA for the proximal consistency step that avoids direct inversion, which is still theoretically sound as our theory tolerates errors for both subsamplers.

Regarding your comments on Q6 and Q7, we would be happy to include more discussion in the final paper when space permits, which we suppressed in the submission and left a few references. Note that the dependence on $d$ in $\varepsilon_{\sf DDS-DDPM}$ can be further improved to $\tilde{O}(\sqrt{d/T}) + O(\varepsilon_{\sf score})$ using sharper results in [BDBDD24]. However, such dependency with $d$ is known to be tight and generally non-avoidable when plain diffusion models are used.

Regarding your comments on Q8, since there is only one day left before the discussion deadline, it is challenging to provide more experiments result in time before the discussion period ends. Nonetheless, we are committed to include more algorithm evaluation in the final version. We also want to provide a bit more discussion regarding the additional experimental results regarding the full evaluation of ReSample on FFHQ dataset, which is already in the rebuttal pdf. Therein, it can be seen that LGD-MC, one of the baseline in our original submission, is a competitive baseline with performance close to that of ReSample. Our algorithm demonstrates significant advantages over both, especially on the phase retrieval task. We expect similar conclusions will hold when we compare ReSample on other datasets/tasks.

Thank you again for your careful review. We appreciate your constructive feedback and are happy to discuss more.
The reviewer thanks the authors for the detailed response. The reviewer is satisfied with the clarifications in Q1, Q2, Q3 and Q9. However, the major concerns still remain. Regarding Q4, the term $\log p(y|x_t)$ is typically not computed and existing methods (e.g. DPS) approximate this using $\log p(y|E[x_0|x_t])$, which is similar to the term $L(x;y)$ used in this paper except the additive score function and a cross-term from the Tweedie's formula. Since the noise level $\sigma$ is anyway not known and needs to be tuned in stepsize, the reviewer suspects that the proposed method doesn't offer any major advantages over existing methods at the cost of more compute. This observation is also supported by experiments in Tables 2 and 3.

For Q5, how would you precompute $A^T\Sigma^{-1}y$ for inverse problems typically considered in practice, such as Gaussian deblur, motion deblur or super-resolution? What would be the storage space complexity in high-dimensional applications where images could be of size 1024x1024? The reviewer thinks that these issues have not been properly addressed in the paper. 

For Q6 and Q7, there is no discussion regarding this important piece of information in the main paper, which would essentially help the reader when to choose this algorithm over others. Especially, the dependence on $d$ in $\varepsilon_{DDS-DDPM}$ seems problematic for large-scale applications with $d$ of the order $10^6$ as in recent state-of-the-art inverse solvers. 

For Q8, the compared baselines are weak and do not adequately justify the claims of the paper. The reviewer was referring to more recently developed pixel space diffusion based inverse solvers such as TMPD or ReSample. The reviewer thanks the authors for providing some preliminary experiments on ReSample. 

The reviewer will follow the guidelines in revising the score if needed after all the questions have been addressed properly.
Dear Reviewer 43oG,

We've taken your initial feedback into careful consideration in our response. Could you kindly confirm whether our responses have appropriately addressed your concerns?

If you find that we have properly addressed your concerns, could you please kindly consider increasing your initial score accordingly? Please let us know if you have further comments.

Thank you for your time and effort in reviewing our work!

Many thanks, Authors
Thank you for addressing my concerns. I think that this work has value and I will keep my original score of 7.
Thank you so much for responding to our rebuttal in a timely manner! We really appreciate it. We also appreciate your voicing of concerns regarding the experiments, and would like to clarify further on the performance DPnP-DDPM. In fact, we have addressed this under the headline of "Additional noise in DDS-DDPM" (we apologize the our algorithm name DPnP-DDPM was misspelled by DDS here due to auto-correction). We have to condense the review due to the character limit, and it is likely that you have missed it and therefore we want to repeat it below:
> - Thank you for your sharp observations. As far as we can see, only the second row in Table 5 has noticeable noise for DPnP-DDPM. However, we also note that in this row, DPnP-DDPM also recovers visibly finer details and textures of the board and the text in the original image. This tradeoff between the capability of reconstructing finer details and the risk of introducing additional noise is indeed a general phenomenon that has been observed in previous works, e.g., in [SKZ+23, page 21].

In addition, we want to highlight that Table 1 is results for phase retrieval, while Table 5 is for quantized sensing, which are **different measurement forward models**. Therefore, their results are not directly comparable, and visually they may appear quite different.

Again, we are happy to include new experiments to further substantiate our paper, if you are willing to provide further feedback. Thank you again for engaging with us!
The experiments presented are inadequate, and certain experimental outcomes are peculiar and fail to substantiate the claims made in the paper. For instance, the results for DPnPDDPM depicted in Table 1 are smooth, whereas those in Table 5 exhibit significant noise. The author did not address these issues during the rebuttal phase.
The reviewer has read the authors' rebuttal as well as the comments from other reviewers. Based on these, the reviewer prefers to maintain the initial score.