Provably Robust Score-Based Diffusion Posterior Sampling for Plug-and-Play Image Reconstruction
Authors: Xingyu Xu
Yuejie Chi
Venue: NeurIPS 2024 poster (2024)
Keywords: score-based generative models
nonlinear inverse problems
posterior sampling
plug-and-play
Primary area: diffusion_based_models

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.

Provably Robust Score-Based Diffusion Posterior Sampling for Plug-and-Play Image Reconstruction
Xingyu Xu' Carnegie Mellon University
Yuejie Chit Carnegie Mellon University
Abstract
In great number of applications, the goal is to infer an unknown image from small number of noisy measurements collected from a known and possibly nonlinear forward model describing certain sensing or imaging modality, which is often ill-posed. Score-based diffusion models, thanks to their impressive empirical success, have emerged as appealing candidate of an expressive prior in image reconstruction. In order to accommodate diverse tasks at once, 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 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 aad---addllyyy play framework in the imaging community, we introduce diffusion plug-and and-play method (DPnP) that alternatively calls two samplers, 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 various tasks. To the best of our knowledge, DPnP is the first provably-robus posterior sampling method for nonlinear inverse problems using unconditional diffusion priors.
Introduction
In a great number of sensing and imaging applications, the paramount goal is to infer an unknown image x= Rd from collection of measurements y E IR that are possibly noisy, incomplete, and even nonlinear. Examples include restoration tasks such as inpainting, super-resolution, denoising, well as imaging tasks such magnetic resonance imaging [DPP07] optical imaging (SEC 15). microscopy imaging [SMCI7), radar and sonar imaging [EPCCOO. 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 [Dn066, [R12], to data-driven ones prescribed in the form of neural networks (UVLI8) [JJPD17]. These priors can be regarded as some sort of generative models for the unknown image, which postu- late 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
x xiinyuxu@andrew.cmm cmu. edu "eeejiee@aale leceandrew. cmu. edu
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Super-resolution (linear)
Phase retrieval (nonlinear)
Quantized sensing (nonlinear)
Figure 1: Solving linear and nonlinear inverse problems with Diffusion Plug- Plug-and-Pa and- (DPnP).
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 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 [SDWMG13), 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 HHA20 SSDK 21 SME20). In particular, (SME220 developed 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-ttpe) 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 large number of steps [SME20).
among one of the most successful
Thanks to the expressive power of score based diffusion models in generating complex and fine- grained images, they have emerged as a plausible candidate of an expressive prior image reconstruc- tion [SSXE22] 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+ 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?
1.1 Our contribution
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, contributions are follows.
Diffusion plug-and-ple for posterior sampling. Motivated by the plug-and-play [BW133 frame- work 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 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 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- adjuustd 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-asympiotic 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 non-asymptotic convergence theory of DPnP, where its performance degenerates gracefully with respect to the errors of the samplers, due to, c.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 wide variety of inverse problems.
Related works. Given its interdisciplinary nature, work sits at the intersection of generative modeling, computational imaging, optimization and sampling. Due to space limits, we postpone the discussion of related works Appendix A
Notation. Let Pe denote the probability distribution of x, and p=(-1)) denotes the conditional
(d) Y to denote random variables X and Y are equivalent in distribution. The matrix La denotes an identity matrix of dimension d. For two probability distributions with density p(z) and q(x), the total variation distance between them is TV(p, q) = J [p(x) - (dr. The x² diverregee of p to q is x2 Il q) = (p(a)>g(a)) dæ. ((m)
distribution of æ given We use X
2 Score-based generative models
In this section, we set up the preliminary on diffusion culfusion-bbasd generative models, which we will be relying upon 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 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 B
2.1 The forward process and score functions
The continuous-time forward diffusion follows the Omsttin-Uhlenbeee (OU) process, defined by the Stochastic Differential Equation (SDE) [SSDK
dx,- = dx + v2dB,, T T Xo ~
where (B+), r is the standard d-dimensional Brownian motion. can be shown that (Do422 Eva 12 the marginal distribution of X, for 2 0 is
Xo + v1 e Xo Xo~p" 5~ N(
It is then clear that the limiting distribution Xxx ~ (() as T oo, i.c. the OU process diffuses Xo ~ p² to the standard normal distribution. The score function of X, is defined by
s(T, Plg.xx (r))
An enlightening property [Vin of the score function is that it can be interpreted as the minimum mean-squared error (MMSE) estimate of e, given 22 = x, fueled by Tweedie's 's formula:
S(T, x) = VI Exo=p"..~N00 (E) (Ele Xu + v1 - e -- x
(4)
=:0(7,a)
Consequently, this makes it possible to estimate the score functions via learning denoise (Hyv05). by estimating the denoising function E(T, :). as typically done in practice (HA20)]
2.2 The reverse process and sampling
To enable sampling, one needs to "reverse" the forward diffusion process. Fortunately, it is possible to leverage classical theory [nnd8 AGSO5) to reverse the SDE, and apply discretization to the time- reversal 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 future sections. Consider SDE given by
dM, = aM,dr + VddB,, T 02 Mg ~ PM
(5)
where a E R and 3 > 0 are constants. For any positive time Too > 0, define the reversed time parameter
= - (T) Too
We are now ready describe the time-reversed processes.
1) The time-reversed SDE of (5) on the time interval [0, Tool is defined
dM = (-aM BV log LMg (M.
dr + v BdB,, T E [0, Tool, M ~ PM (7)
where B is independent of B, i.e., another Brownian motion. It is classical result And82
that the reversed process r shares the same path distribution as M, i.c., (M²") [e(0
(M).)_] In other words, the joint distribution of (M .M .... - for any 0 S S
T2 ..... T I Too for any integer is 2 1, coincides with that of (Mr, M M,
2) In place of the reversed SDE in it is possible to consider the following probability flow ODE AGSO5 SSDK+211
pro ( (-aM
dM
M.e Mrev
log pM (T M dr, T [ool, MTe PM
(8)
The reversed ODE satisfies slightly weaker guarantee than that of the reversed SDE, which
nevertheless suffices for most practical purposes [SSDK*21] Mrev (d) M,, E |0, Too). Note that the reversed ODE only guarantees identical marginal distribution for each M whereas the reversed SDE guarantees identical joint distribution.
Specializing the above to the OU process with proper discretization then leads to popular samplers used for generation, as follows.
DDPM-type stochastic samplers. Specializing the time-reversed SDE
to the OU process gives
d. play (X + 2s(r X )) dr + V2dB,, [0, Tool:
X ~ Px+
As Too - it can be seen from that px converges to (0 Ia). Thus the solution of the above SDE can be approximated by initializing X NNO Ia) instead. The DDPM sampler (HA20))
can be viewed as discretization of this SDE [SSDK+2] +21
DDIM-type deterministic samplers. On the other hand, the probability flow ODE 8 for the OU process reads as
dx rey
(X + 8 X )dd,, T E [0,
X .x, Too
Again, as Too + DC, one may approximate the initialization with X ~ NO Ia). is known that the popular DDIM sampler [SSDK 21, SME20) is a discretization of this ODE [C222. The ODE-based deterministic samplers allow more aggressive choice of discretization schedules, well as fast ODE solvers [LZB +2), enabling significantly accelerated sampling process compared to the SDE-based stochastic samplers.
Posterior sampling via diffusion plug-and-play
We are interested in solving (possibly nonlinear) inverse problems, where the aim is to infer an unknown image x² E B from its measurements y R
y A(x")) +
where A : d + R 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.,
x² p" (x),
(10)

---

The posterior distribution given measurements is defined as
p" (x)l) v) x (x) plylæ* = ==) =x)=
(11)
Here, L(:;8) y) is the log- log-liiel 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, also subsumes the "reward-guided" or "loss-guided" setting [3]Y=23), where L can be viewed as reward function 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 L(x; < co.
Goal. Our goal is to sample T\ from the posterior distribution T p" y) given estimates ((T, (resp. z)r of the unconditional score functions s(T, x:) (resp. the noise function c(T, x)) in assuming knowledge of the likelihood function L( y).
3.1 Key ingredient: score-based denoising posterior sampling
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 by
Tnoisy x² + +., {~ {~N(0, 17"
(12)
where 7) > is the noise intensity assumed to known. Our goal is to sample from p" (-1 Tnoisy) given the score estimates 5<(x) (resp. the noise estimates & (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 score- based 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 stochastic DDPM- type sampler for denoising, termed DDS-DDPM. We divide our development into the following steps.
1) Step I: introducing the heat flow. Let us introduce a heat flow with initial distribution p", defined by the following SDE:
dy, dB,, T 2 0, Yo ~~mp"
(13)
where (Bz)r20 is the standard d-dimensional Brownian motion. The solution of (13 is simply
Y, - Yo + B,,
T T20
(14)
(d)
Since Br ~ NO rTa). it readily follows that B.p E, which together with Yn ~ p" yield the important observation that Fniisy xx + can be viewed as endpoint of the heat flow, in the
(d)
sense that Tnoisy = xx +
Y,p.
2) Step 2: reversing the heat flow. Following similar reasonings in Section the next step boils down to reverse the heat flow
(13). The reverseversal of the heat flow SDE (13) is (cf. (7) given by
dy VlogPy,, (Y,5",, )dr + dB,, T TE [0.
Y PY,2
(15)
where (B=) is an independent copy of (Br)r2o. As introduced earlier, the virtue of the time-reversed SDE (15) is that produces process yrev with the same path distribution as Y,,
i.e., (yrev) ree(0.")) (Y=)rejo.v?)) In particular, the joint distribution of (Yrev, Y!S") is the same
as that of (Yo, X,,)) (d) (x" x. This implies that the conditional distribution p" p- (noisy) is the same Py (- y = xn Surprisingly, the latter admits simple interpretation: Pyge (- y @noisy) is the distribution of Y when we initialize (5) with y I
Therefore, sampling the posterior p² |/noisy)) amounts to solving the following simple SDE:
dy = V log loggp,,,,, (X,"",)d + dB,, T 10, Y Znoisy.
(16)
5

---

3) Step 3: connecting the score functions. It is now immediate to arrive at our proposed stochastic sampler DDS- by discretization of this SDE (T6), which requires knowledge of the score functions V log PY, (+). A key observation is that they can in fact be computed from the score
function s(T, x) (ef. (3) due to the following lemma, whose proof provided in Appendix D.1
Lemma (Score function of Y,). For 2 O, we have
V logpy, (x)
log +T),
VI+
VI+7
The resulting sampler, DDS-DDPM, is summarized in Algorithm 2 (deferred in the appendix)
using a discretization procedure with an exponential integrator (ZC22
deterministic DDIM-type sampler via OU process. We next develop a deterministic DDIM-type sampler for denoising. termed DDS-DDIM.
1) Step I: introducing posterior-initializer OU process. To sample from the posterior distribution p' (-]Znossy), we first introduce a random variable w which has (unconditional) distribution
Pw (x) = (= x I + Inoisy):
(17)
in the same form of the desired posterior distribution (-pattioon"((((ooooop(ooa (-)yyo)s)) Here, since the noisy observation T noisy is given, we regard it as fixed. We then further introduce = w - Inoisy. which a "centered"
version of w, whose distribution is
(2(x) (u + Tnoisy) p- (x* + Fooisy x + Inoisy). )
The OU process with initial distribution p, is defined by the SDE:
dz. = -2dr dB,, 7 2 0, Zo ~ Pz:
(18)
where B, is the standard d-dimensional Brownian motion. As in ( the marginal distribution of Z, is given by
Z, (d) - Zp + e1- Zo~ ~ Pa, ~~ NO), T 2 0.
(19)
2) Step 2: reversing the OU process. Following similar reasonings in Section reversing the OU process (88) will enable us to generate samples 2 Pz. Then
+ Tnoisy which, by definition, has distribution Pu defined (17). and is a sample from the desired posterior distribution p" (-|Znnoisy). We are thus led to solve the time-reversed probability flow ODE (cf. of (18) given by
1
7
V Pz, er (.pt) dr, E [0, Too]:
Zrev
N(0, Ia),
(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 V log pz, (.), which again can be computed from the score function s(T, x)
(cf. ) documented by the following lemma, whose proof is provided in Appendix D.2)
Lemma 2 (Score function of 77), For T 2 O, have
2Tx
Enoisy + 0--,1p20
V log (x)
172 ²TT -
112 + 2 -
8
7,
172 + c -
(21)
where
12((22 11) 172 + 2 +1
11)
7 = 7(T) = log
(22)
After plugging this into 20 and solving the ODE for Ztey, we see that Z Tnoisy is desired sample from the posterior distribution p* (+)xnoisy)) 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)
Measurements E Rm, log-likelihood function L(5v) of the forward model, score estimates 3, annealing schedule (1k JoccK Initialization: Sample To ~ (00 4 Ia) Alternating sampling: for 0, 1.,,..., K do
Input:
(1) Proximal consistency sampler: Sample Tx+ oc exp (Z(::)) y) 2:1 -xull using sub-
routine PCS(0 y, C, m)) (Alg-)
(2) Denoising disfusion sampler: Sample Te+1 exp (logr* (2) शमत || Fx+ 12
using
subroutine DDS-DDPM
(Alg. - 2 or DDS-DDIM( 8, ma) (Alg. 3
Output: Tg.
3.2 Our algorithm: diffusion plug-and-play
Now we turn to the general setting where the mcasurement operator A is arbitrary. From the factorization of posterior distribution in (1I), one intuitively understands that posterior sampler must obey two constraints simultaneously: (i) the data prior constraint, corresponding to the first factor p" (Ir), which imposes that the posterior sampler should be less likely to sample at those points where p² (ar) is small; (ii) the measurement consistency constraint, corresponding to the second factor ef(wyy) which imposes that A(x) R y-
Diffusion plug-and-play (DPnP). We will apply the idea of alternatively enforcing these two constraints from sampling perspective in the same spirit of [VDCI9 LST21 [ST2]|BB3]]] dubbed diffusion auu-and--lay (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 Tg and the annealing parameter 77 at the k-th iteration, DPnP proceeds with the following two steps:
[B23]. Our algorithm.
(i) Proximal consistency sampler to enforce the measurement consistency constraint. DPnP draws sample ē from the distribution proportional to exp [[x; y) - - to promote the image be consistent with the measurements. This step, which we denote the proximal consistency sampler, can be achieved by small modifications of standard algorithms such as Metropolis-Adjusted Langevin Algorithm (MALA) [R998] given in Algorithm (deferred in the appendix).
(ii) Denoising diffusion sampler enforce the data prior constraint. DPnP next draws sample Tx from the distribution proportional to
exp (-(- log (logp((x): (-(-logp"(x)
+ 2n | T Tx+ (1)) x px = x + Th w Fx++z)
(23)
to promote the image to be consistent with the prior, where NN0, I2). The last step. which follows from the Bayes' rule, makes clear that this step can be precisely achieved by the denoising diffusion sampler (developed in Section 3.I using solely the unconditional score function.
Combining both steps lead to the proposed DPnP method described in Algorithm ! about the proposed DPnP method are in order.
Some comments
The proximal consistency sampler PCS can be viewed "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
->xel2)) is attained at the point x' TR satisfying
2n
V Ve (z' xx) 0,
x' = E T V C(x's 8) ~ T4 + m V a2 L(
Therefore, the proximal consistency sampler draws random samples "concentrated" around r', which approximates the implicit proximal point update, akin gradient step at Zy.
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 to reside in M, like the proximal map. To see this, note that denoising posterior distribution vanishes outside M by 23
7

---

The proximal consisteney sampler PCS admits a simple form when the forward model A is linear,
i.e. A(x)) Az for some matrix A E R mmx and the measurement noise & ~ (N Σ) is Gaussian. In this situation, the proximal consistency sampler PCS can be implemented directly by
xR+$ = PCS (TK, U: L, n1) = T. + 2 Wkk
Wk ~ ~N Ia),
where 3x+ (AT- #4)) (A + +xx and Σn (AT2 A +
4 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 A(T, is accurate, i..., s(T,,=) s(T, x), and assume the ODE/SDES in DDS and PCS are solved exactly. Let (E:)r2o be a decreasing sequence of positive numbers satisfying lim; Et O, and (k))220 be an increasing sequence of integers with kg 0. Set the annealing schedule Ti = Et. for
- S. k < Ki, = 1,2, ... Let minp:=1,2... ik ky 1| + co, the output Th, of DPnP converges in distribution to the posterior distribution p² (-ly) for I 00.
In words, Theorem establishes the asymptotic consistency of DPnP under fairly mild assumptions on the forward model (cf. Assumption I: as long as the sampled distributions of DDS and PCS are exact, then running DPnP with slowly diminishing annealing schedule of (n²) will output samples approaching the desired posterior distribution p" (-1y) when the number of iterations I goes to infinity.
Non-asymptott error analysis. We now step away from the idcalized 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. 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(C(x; y) = ||x
Fx 12) (resp. p² (a* 2/x* + e = * )) 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 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 aoompssymto robustness of DPnP). With the notation in DPnP (Algorithm 7 set Th = 1) O. Under Assumption [ there exists X = A(p", f. E (0, 1), such that the following holds. Define a stationary distribution Tn by (n((x) x p² (x) (I), where On defined by
(y(x)) == e( * PPne(a) = (2x)d/2mpd
O(d1:v)- -te 11== da',
6 6~ *~N(0, Ta):
(24)
where * denotes convolution. If PCS has error at most EPCS in total variation and DDS has error at most EDDS in total variation per iteration, then for any accuracy goal Eacc > O, with K x (og(1/
we have
1
TV(pax mn) S Eace v (ps, [1)) + 1-X (EDOS + EPCS) log
(25)
Eacc
Before interpreting Theorem 2. we observe that (z(x) = ((aa)) thus To(x) x p (z)ee((r))) coin- cides with the desired posterior distribution p² (.1y). Thus Theorem tells us that, assuming constant annealing schedule 7/&- = y), the output of DPnP converges in total variation to the distribution Ty, which is distorted version of the desired posterior distribution up level n, with sufficiently many iterations. A few remarks are in order.
-in-dinin 7/. It can be seen from Theorem 2 that even with nonzero 7, DPnP already enforces the data prior strictly. On the other hand, the measurement consistency distorted by an order of 17. This is usually tolerable, since the measurements are themselves contaminated by
8

---

Table Samples of different algorithms for phase retrieval, where DPnP generate images of higher quality and recover fine details of the image more faithfully than the state-of-thee DPS [CKM 23 and LGD-MC [SZY 23] algorithms. Pixel-based ReSample [SKZ* did not generate meaningful
images, possibly due to high nonlinearity of phase retrieval.
DPnP-
DPnP- DDIM
Ground
Input
DPS
LGD-MC
ReSample
DDPM
truth
Table 2: Samples of different algorithms for quantized sensing.
DPnP- DDPM
DPnP-
Ground truth
Input
DPS
LGD-MC
ReSample
DDIM
ANME
USE
USE
USE
USE
USE
noise, thus when is smaller than the noise level, the distortion would be tolerable. In practice, is beneficial to choose annealing schedule of {m) J. which will be elaborated in Section
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.
5 Numerical experiments
We provide preliminary numerical evidence 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 DPnP-DDPM and DPnP-DDIM respectively.
Table 3:
Evaluation of solving inverse problems FFHQ 256 256 validation dataset (lk sam- ples). Despite considerable efforts to optimize parameters, pixel-based ReSample did not generate meaningful results for phase retrieval.
Super-resolution (4x, lincar)
Phase retrieval (nonlincar)
Quantized sensing (nonlincar)
Time
per sample
Algorithm
PSNR
LPIPS
PSNR
LPIPS
PSNR
DPnP-DDIM (ours)
0.301
24.2
0.376
22.4
0.293
24.2
90s
DPS [CKM*29
0.331
23.1
0.490
17.4
0.367
21.7
~ 60s
LGD-MC (n
23
0.318
23.9
0.522
16.4
0.317
23.9
60s
ReSample spixcl-bused SKZ 23
0.313
23.9
0.318
22.6
70s

---

Table 4: Evaluation of solving inverse problems on ImageNet 256 256 validation dataset (1k samples), Despite considerable efforts to optimize parameters, pixel pixel-based ReSample did not generate meaningful results for phase retrieval.
Super-resolution (4x, linear)
Phase retrieval
Quantized sensing (nonlinear)
Time
(nonlinear)
per sample
Algorithm
LPIPS
PSNR f
LPIPS 4
PSNR t
LPIPS :
PSNR
DPnP-DDIM (ours)
0.416
21.6
0.562
13.4
0.363
23.0
240s
DPS [C
0.473
20.2
0.677
13.4
0.542
18.7
150s
LGD-MC (n 5)
23
0.416
20.9
0.592
12.8
0.384
22.3
150s
ReSample (pixel- based) ISKZ 23
0.464
20.1
0.414
19.8
180s
Experimental setups. We compare DPnP with the state-of-the art algorithms including DPS (3)M*23), LGD-MC 231. and pixel-based ReSample [SKZ*23] for super-resolution (lin- car), phase retrieval (nonlinear), and quantized sensing (nonlinear). Definitions of these forward measurement models are in Appendix G.2 The annealing schedule {Th, J of DPnP is fixed across all tasks (Appendix H.2), while DPS, LGD-MC, and ReSample are fine-tuned with reasonable effort for best performance, All experiments are run on a single Nvidia L40 GPU. More details and experiments are in Appendix G
Sample images. We present the sample images of different algorithms for the most complicated task of phase retrieval. For phase retrieval, Fourier transform is performed to the image with coded mask, and only the magnitude of the Fourier transform is taken as the measurement [SEC* 15). Due to the nonlinearity of the operation of taking magnitude. the forward model is nonlinear. The samples generated by different algorithms are shown in Table [
We also present the sample images for the nonlinear problem of quantized sensing. In quantized sensing, each pixel of the image randomly dithered and then quantized to one bit per channel. Nonlinearity of quantizing renders this forward model nonlinear. The samples generated by different algorithms are shown in Table 2
Evaluation. We evaluate the performance of DPnP on the FFHQ validation dataset [KLA19] and the ImageNet validation dataset [RDS+15) Since DPnP-DDIM has similar performance with DPnP- DDPM but admits much faster implementation, only DPnP-DDIM is evaluated. The LPIPS and PSNR are shown in Table 3 and Table These two metrics are arguably the more relevant ones for solving inverse problems. For comparison under other metrics such as FID, SSIM, cf. Appendix G.4
It can be that, DPnP is capable of solving both linear and nonlinear problems, and, in comparison with prior tfat-oof-the-ar performs better in recovering fine and crisper details.
6 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 hetween 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.
Acknowledgments and Disclosure of Funding
This work is supported in part by Office of Naval Research under N00014-19-1-240 and by National Science Foundation under DMS-213400 344800 and ECCS-21 X. Xu is also gratefully supported by the Axel Berny Presidential Graduate Fellowship Carnegie Mellon University.
References
[AGSO5]
L. Ambrosio, N. Gigli, and G. Savaré. Gradient flows: in metric spaces and in the space of probability measures. Springer Science & Business Media, 2005.
[And82] B. D. Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Applications, 12(3):313-326, 1982.
10

---

[BB23] C. A. Bouman and G. T. Buzzard. Generative plug and play: Posterior sampling for inverse problems. arXiv preprint arXiv:2306.072:3 072333 2023.
[BCSB18] G. T. Buzzard, S. H. Chan, S. Sreehari, and C. A. Bouman. Plug-and-play unplugged: Optimization-Fres reconstruction using consensus equilibrium. SIAM Journal on Imaging Sciences, 113)::001-20200 2018.
[BDBDD24] J. Benton, V. De Bortoli, A. Doucet, and G. Deligiannidis. Nearly d-linear convergence bounds for diffusion models via stochastic localization. In The Twelfth International Conference Learning Representations, 2024.
[BJPD17] A. Bora, A. Jalal, E. Price, and A. G. Dimakis. Compressed sensing using generative models. In International conference machine learning, pages 537-546. PMLR, 2017.
[Bou23a] N. Bourbaki. Théories spectrales Chapitres 1 el 2. Springer, 2023.
[Bou23b] N. Bourbaki. Théories spectrales: Chapitres 3 5. Springer Nature, 2023.
[CCL+ S. Chen, S. Chewi, J. Li, Li, A. Salim, and A. R. Zhang. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. arXiv preprint arXiv:2209. /1215, 2022.
[CCL+233 S. Chen, S. Chewi, H. Lee, Y. Li, J. Lu, and A. Salim. The probability flow ODE is provably fast. Neural Information Processing Systems, 2023.
[CDC23] F. Coeurdoux, N. Dobigeon, and P. Chainais. Plug-and-play split Gibbs sam- pler: embedding deep generative priors in Bayesian inference. arXiv preprint arXiv:2304. 11.304.111333 11134, 2023.
[CDD23] S. Chen, G. Daras, and A. Dimakis. Restoration-degraaattoo beyond linear diffusions: A non-asymptotic analysis for DDIM-type samplers. In International Conference on Machine Learning, pages 4462-4484, 2023.
[CICM23] G. Cardoso, Y. E. Idrissi, S. L. Corff, and E. Moulines. Monte carlo guided diffusion for Bayesian linear inverse problems. arXiv preprint XXXv:23088.09833 2023.
[CKM+233 H. Chung, J. Kim, M. T. Mccann, M. L. Klasky, and C. Ye. Diffusion posterior sampling for general noisy inverse problems. In International Conference on Learning Representations, 2023.
[CLA+21] S. Chewi, C. Lu, K. Ahn, X. Cheng, T. Le Gouic, and P. Rigollet. Optimal dimension dependence of the metropolis- langevin algorithm. In Conference on Learning Theory, pages 1260-1300. PMLR, 2021.
[CLL23] H. Chen, H. Lee, and Lu. Improved analysis of score-based generative model- ing: User-friendly bounds under minimal smoothness assumptions. In International Conference on Machine Learning, pages 4735 4763, 2023.
[CLS15] E. J. Candes, X. Li, and M. Soltanolkotabi. Phase retrieval via Wirtinger flow: Theory and algorithms. IEEE Transactions on Information Theory, 11842:7985 2015.
[CLY23] H. Chung, S. Lee, and J. C. Ye. Fast diffiusion sampler for inverse problems by geometric decomposition. arXiv preprint arXiv:2303.00555 2023.
[CR12] E. Candes and B. Recht. Exact matrix completion via convex optimization. Communi- cations of the ACM, 55(6):111-119) 2012.
[Don06] D. L. Donoho. Compressed sensing. IEEE Transactions on Information Theory.
52(4):1289-1308 2006.
[Doo42] I.. Doob. The Brownian movement and stochastic equations. Annals of Mathematics. 43(2):351-369, 1942.

---

[Drul7] D. Drusvyatskiy. The proximal point method revisited.
arXiv preprint
arXiv:1712.06035, 2017.
[DS24] Z. Dou and Y. Song. Diffusion posterior sampling for linear inverse problem solving: filtering perspective. In International Conference on Learning Representations, 2024.
[Efr1 B. Efron. Tweedie's formula and selection bias. Journal of the American Statistical Association, 106(496):1602-1114, 2011.
[val12] L. Evans. An introduction to stochastic differential equations, volume 82. American Mathematical Soc., 2012.
[FBS24] Z. Fang. Buchanan, and J. Sulam. What's in a prior? Learned proximal net- works for inverse problems. In The Twelfth International Conference on Learning Representations, 2024.
[3SR+23) B. T. Feng, J. Smith, M. Rubinstein, H. Chang, K. L.. Bouman, and W. T. Freeman. Score-based diffusion models as principled priors for inverse imaging. In 2023 IEEE/CVF International Conference on Computer Vision (ICCV), pages 10486-10497, 2023.
[4] S. Gupta, A. Jalal, A. Parulekar, E. Price, and Z. Xun. Diffusion posterior sampling is computationally intractable. arXiv preprint arXiv:2402.12727, 2024.
[GL10] K. Gregor and Y. LeCun. Learning fast approximations of sparse coding. In Pro- ceedings of the 27th international conference on international conference on machine learning, pages 399 406, 2010.
[GMJS22] A. Graikos, N. Malkin, N. Jojic, and D. Samaras. Diffusion models as plug-and- and- priors. Advances Neural Information Processing Systems, 35: 15:14715-14728, 2022.
[HJA20] J. Ho, A. Jain, and P. Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 33:6840-6851, 2020.
[HSMC17) J. Huang, M. Sun, J. Ma, and Y. Chi. Super-resolution image reconstruction for high-density three-dimensional single-molecule microscopy. IEEE Transactions on Computational Imaging, 3(4):763-7733 2017.
[Hyv05] A. Hyvärinen. Estimation of non-normalized statistical models by score matching.
Journal of Machine Learning Research, 6(4), 2005.
[KEES2]] B. Kawar, M. Elad, S. Ermon, and Song. Denoising diffusion restoration models. Advances in Neural Information Processing Systems, 35:23593-23600 2022.
[Key81] R. Keys. Cubic convolution interpolation for digital image processing. IEEE transac- tions on acoustics, speech, and signal processing, 29(6):1133-111 1981.
[KHR23] F. Koehler, A. Heckett, and A. Risteski. Statistical efficiency of score matching: The view from isoperimetry. International Conference on Learning Representations, 2023.
[KLA19] T. Karras, S. Laine, and T. Aila. A style-based generator architecture for generative adversarial networks, In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 4401-440 2019.
[КЕЕ21] B. Kawar, G. Vaksman, and M. Elad. Stochastic image denoising by sampling from the posterior distribution. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 1866-1875. 2021.
[LBA+22] R. Laumont, V. D. Bortoli, A. Almansa, J. Delon, A. Durmus, and M. Pereyra. Bayesian imaging using plug & play priors: when Langevin meets Tweedie. SIAM Journal on Imaging Sciences, 15(2):701-737, 2022.
12

---

[LDP07] M. Lustig, D. Donoho, and J. M. Pauly. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magnetic Resonance in Medicine: An Officcall Journal of the International Society for Magnetic Resonance in Medicine, 58(6): 188-1195, 2007.
[LHE*24] G. Li, Y. Huang, Efimov, Y. Wei, Y. Chi, and Y. Chen. Accelerating convergence of score-basee diffusion models, provably. arXiv preprint arXiv:2403.0:88 2024.
[LHW24] G. Li, Z. Huang, and Y. Wei. Towards a mathematical theory for consistency training in diffusion models. arXiv preprint 2xiv:2402.07802. 2024.
[LST21] Y. T. Lee, R. Shen, and K. Tian. Structured logconcave sampling with restricted Gaussian oracle. In Conference on Learning Theory, pages 2993-3050. PMLR, 2021.
[LWCC23] G. Li, Y. Wei, Y. Chen, and Y. Chi. Towards faster non-asympiotic convergence for diffusion-based generative models. arXiv preprint arXiiv2330600925/ 2023.
[2Z C. Lu, Y. Zhou, F. Bao, J. Chen, C. Li, and J. Zhu. DPM-Solver: A fast ODE solver for diffusion probabilistic model sampling in around 10 steps. Advances in Neural Information Processing Systems, 35:5775-5787, 2022.
[MLE21] V. Monga, Y. Li, and Y. C. Eldar. Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing, IEEE Signal Processing Magazine, 38(2):18- 44, 2021.
[MSKV24] M. Mardani, J. Song, J. Kautz, and A. Vahdat. A variational perspective on solving inverse problems with diffusion models. In The Twelfth International Conference on Learning Representations, 2024.
[MV19] O. Mangoubi and N. K. Vishnoi. Nonconvex sampling with the metropolii-adjussed langevin algorithm. In Conference on learning theory, pages 2259-2293, PMLR, 2019.
[PEPC10] L. C. Potter, E. Ertin, J. T. Parker, and M. Cetin. Sparsity and compressed sensing in radar imaging. Proceedings of the IEEE, 98(6): 1006-10000 2010.
[PRS+24] C. Pabbaraju, D. Rohatgi, A. Sevekari. Lee, A. Moitra, and A. Risteski. Provable benefits of score matching. Advances in Neural Information Processing Systems, 36, 2024.
[PW24] Y. Polyanskiy and Y. Wu. Information theory: From coding to learning. Cambridge university press, 2024.
[2BL+22] R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer. High-resolution image synthesis with latent diffusion models. In JEEE/CVV Conference on Computer Vision and Pattern Recognition, pages 10684-10695, 2022.
[RDN+22] A. Ramesh, P. Dhariwal, A. Nichol, C. Chu, and M. Chen. Hierarchical text- conditional image generation with CLIP latents. arXiv preprint 22411:2204.06125, 2022.
[RDS 15] O. Russakovsky, J. Deng, H. Su, J. Krause, Satheesh, Ma, Z. Huang. A. Karpathy, A. Khosla, M. Bernstein, A. C. Berg, and L. Fei-Fei. ImageNet Large Scale Visual Recognition Challenge. International Journal of Computer Vision (IJCV), 115(3):211- 252, 2015.
[REM17] Y. Romano, M. Elad, and P. Milanfar. The little engine that could: Regularization by denoising (RED). SIAM Journal on Imaging Sciences, 10(4):1804-184) 2017.
[RR98] G. O. Roberts and J. S. Rosenthal. Optimal scaling of discrete approximations to Langevin diffusions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255-268, 1998.
[RS18] E. Reehorst and P. Schniter. Regularization by denoising: Clarifications and new interpretations. IEEE transactions computational imaging, 5(1):52-67, 2018.
13

---

[SC97] L. Saloff-Coste. Lectures on finite Markov chains, pages 301-413. Springer Berlin Heidelberg, Berlin, Heidelberg, 1997.
[Sch12] H. Schaefer. Banach Lattices and Positive Operators, volume 215. Springer Science & Business Media, 2012.
[SCS+22] C. Saharia, W. Chan, S. Saxena, L. Li, Whang, E. L. Denton, K. Ghasemipour, R. Gontijo Lopes, B. Karagol Ayan, T. Salimans, et al. Photorealistic text-to-image diffusion models with deep language understanding. Advances in Neural Information Processing Systems, 35:36479 0944449999944,,, 2022.
[SDWMG15] Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International Conference Machine Learning, pages 2256-2265, 2015.
[SE19] Y. Song and S. Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 32, 2019.
[SEC*15] 15] Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev. Phase retrieval with application to optical imaging: contemporary overview. IEEE signal processing magazine, 32(3):87-109, 2015.
[3KZ*23] B. Song, S. M. Kwon, Z. Zhang, X. Hu, Q. Qu, and L. Shen. Solving inverse problems with latent diffusion models via hard data consistency. arXiv preprint M.Xi:22307.08123, 2023.
[SME20] Song, C. Meng, and S. Ermon. Denoising diffusion implicit models, In International Conference on Learning Representations, 2020.
[SSDK+21] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, Ermon, and B. Poole. Score- based generative modeling through stochastic differential equations. International Conference on Learning Representations, 2021.
[SSXE21] Y. Song, L. Shen, L. Xing, and S. Ermon. Solving inverse problems medical imaging with score-basee based generative models. In International Conference Learning Representations, 2021.
[SVMK22] J. Song, A. Vahdat, M. Mardani, and J. Kautz. Pseudoinverss-guidee diffusion models for inverse problems. In International Conference on Learning Representations, 2022.
[3] Y. Sun, Z. Wu, Y. Chen, B. T. Feng, and K. L. Bouman. Provable probabilistic imaging using score-based generative priors. arXiv preprint arXiv:2310.108: 2023.
[3Z Song, Q Zhang, H. Yin, M. Mardani, M.-Y. Liu, J. Kautz, Y. Chen, and A. Vah- dat. goss-guided diffusion models for plug-and-play controllable generation. In International Conference on Machine Learning, pages 32483-32498. PMLR, 2023.
[Tie94] L. Tiemnyy Markov Chains for Exploring Posterior Distributions. The Annals of Statistics, 22(4):1701 1728, 1994.
[TYT+222 B. I.. Trippe, J. Yim, D. Tischer, D. Baker, T. Broderick, R. Barzilay, and T. Jaakkola. Diffusion probabilistic modeling of protein backbones in 3d for the motif-scaffolding problem. arXiv preprint arXiv:2206.04//9 19, 2022.
[TZ24] W. Tang and H. Zhao. Contractive diffusion probabilistic models. arXiv preprint arXiv:2401.13115, 2024.
[UVL18] D. Ulyanov, A. Vedaldi, and V. Lempitsky. Deep image prior. In Proceedings the IEEE conference on computer vision and pattern recognition, pages
2018.
[VBW13] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg. Plug-and-play priors for model based reconstruction. IEEE Global Conference on Signal and Information Processing, pages 945-949 IEEE, 2013.
14

---

[VDCI9] M. Vono, N. Dobigeon, and P. Chainais. Split-and-augmente Gibbs sampler- application large-scale inference problems. IEEE Transactions on Signal Processing. 1766)1164-1661, 2019.
[Vin1 P. Vincent. A connection between score matching and denoising autoencoders. Neural computation, 23(7):661111774, 2011.
[Wib18] A. Wibisono. Sampling as optimization in the space of measures: The Langevin dynamics composite optimization problem. In Conference Learning Theory, pages 2093-3027. PMLR, 2018.
[3T L. Wu, B. Trippe, C. Naesseth, D. Blei, and J. P. Cunningham. Practical and asymptot- ically exact conditional sampling in diffusion models, Advances in Neural Information Processing Systems, 36, 2023.
[7C22] Q Zhang and Y. Chen. Fast sampling of diffusion models with exponential integrator. In The Eleventh International Conference on Learning Representations, 2022.
15

---

A Related works
Algorithmic unrolling and plug-and-play image reconstruction. Composite optimization algo- rithms, 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. [G.IO] advocates the perspective of algorithmic unrolling, which turns an iterative algorithm into concatenations of linear and nonlinear layers like in neural network. [VWW13 recognized that the proximal mapping step in many composite optimization algorithms can be regarded as a denoiser 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-ani and-play methods [BCSB18]; see [MLE2]] for review.
Regularization by denoising and score matching. [ pointed out a connection between score matching and image denoising, which is consequence of the Tweedie's formula [Eiil] 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; [SS18 later clarified that the RED framework can be interpreted as score matching by denoising using the Tweedie's 's formula. [VE22] developed 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, tack- ling 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 perspec- tive, 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]. [CDC233 extended the split Gibbs sampler [VDCI9] in the plug-and-play framework, and advocated the use of score-based diffusion models such as DDPM [HA)20 for image denoising based on heuristic observations. In contrast, rigorously derive the denoising diffusion samplers from first principles, unraveling critical gaps from naive applications of the generative samplers to denoising, and offer theoretical guarantees 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 pro- cess and projecting onto the measurement constraint, e.g., [CKM+23 KEES22 S3 CLY23 GMJS22 SVMK22). However, these approaches do not possess asymptotic consistency guaran- tees. +SZY+23) proposed to use multiple Monte Carlo samples to reduce bias. On the other hand, [CICM23] developed Monte Carlo guided diffusion methods for Bayesian linear inverse problems which tend to be computationally expensive, and [DS24] recently introduced a filtering perspective and applied particle filtering. Although asymptotically consistent, these approaches are limited to linear inverse problems. 2T WTN 23) introduced sequential Monte Carlo (SMC) algo- rithms for conditional sampling using unconditional diffusion models that are asymptotically exact. [SKV244 developed a variational perspective that connects to the regularization by denoising framework. (GJP+24) showed that the worst-case complexity of diffusion posterior sampling can take super-polynomial time regardless of the algorithm in use.
Theory of diffusion models and score matching. A number of recent papers have studied the non- asymptotic convergence rates of popular diffusion samplers, including but not limited to stochastic DDPM-type typee samplers [CCL (CL+22 CLL23 LWCC23 BDBDD24| TZ24), deterministic DDIM-type
samplers LWCC23 [CCL CDD23], and accelerated samplers [LHE" 22 LHW24. In addition,
the statistical efficiency of score matching has also been investigated [KHR23 PRS+24
B Discrete-time formulation of diffusion processes
The discrete-time forward and backward diffusion processes can be understood as two processes constructed in the following manner:
1) a forward process
ZT
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;
16

---

2) reverse process
2/ 2--11
that starts from pure noise (e.g., standard Gaussians) and converts it into samples whose distribution close to the target image distribution.
Consider the forward Markov process in R that starts with sample from the data distribution px. and adds noise over the trajectory according
Fo ~ p- It /1 Be x1-1 VB wi,
(26a)
1 1
(26b)
where {wy (whis'''s are independent standard Gaussian vectors, i.e., we (0, 1)} describes the noise-injection rates used in each step. Therefore, as
iid NO Ia), and (3e
write zı equivalently
Vā zo + VI- EtEt, Et Lo), 0,0, Here, (0,1 is the schedule diffusion given by
(27)
0t == 1- Be,
II ok
1STST.
(28)
Clearly, it verifies that 1 2 0g > á1 ... oT > 0. As long as är is vanishing, it is easy to observe that the distribution æT approaches N(O,
Score functions. As will be seen, in order to sample from p², it turns out to be sufficient to learn the score functions of pm, each step of the forward process, defined
s (x)
V Po, Px,())) (z),
t= 0,1, ,T.
(29)
As in the continuous time, the score function can be viewed as MMSE estimate:
s (z)
at Enowp' Exo-p' sc- Vā - VI âūE; = x)
(30)
(x)
Comparing 27 and 2
be checked that the discrete-time diffusion process be embedded into the continuous-time one via the time change
f +
log
2
in the sense that
log
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 (+))) Rd Rd, = 1...., T such that 8 = st in analogy with the continuous time counterpart.
C Details of algorithm subroutines
With the discrete-time perspective established in Appendix B description and the implementation of our algorithms DDS-DDPM (Algorithm 2b. DDS-DDIM
ready present the detailed
(Algorithm 3) and PCS (Algorithm 4
D Score functions of diffusion denoising samplers
D.1 Proof of Lemma i
Proof. The marginal distribution (14) of the heat flow can be written as
(d) Yo + Vre, Yo pr 6 N((0
Y,
(31)
Comparing and 31 it is not hard to check that
Y, (d)
VI+TX log(1+r)
17

---

Algorithm 2 Denoising Diffusion Sampler (stochastic) DDS-DDPMI
Tnnisy 8,11) 17)
Input: noisy data Tooisy E R score estimates 3 == (-((-) R + R 1...,, T or noise
estimates غ = {+e(-) IR + R t = I...., T), and noise level r >
Scheduling: Compute the diffusion schedule (Ta)oSAST) by
Tt = ū 1, 0S1ST,
where
1
T = max t:0 T, ai 172 +
Initialization:
Set DTV = I nnoisy.
Diffusion: for 1. T", T ,..., 1 do
Tt-1 Tt 22( Te + VT. Te-1 wa,
Wt ~ N0
where
& = (v 84)
Se (Vār
VI- at
Output:
Z0-
Algorithm 3 Denoising Diffusion Sampler (deterministic) DDS-DDIM DZnMMM 17)
Input: noisy data Tnoisy E Rdd score estimates 3 == (-x(-) Rd + Rd t 1, I,..., T3} or noise
estimates ē = (-)(-)) IRd + R4 1= I,.... T), and noise level ry > 0.
Scheduling: Compute the diffusion schedule (ū.)oSesT) by
(1²2 1)ā
" 1/2 + Cg
OSIST,
where
1 T' = max 1:00t5 8
++1
Initialization:
Draw 2T N (0, la).
Diffusion: for t = T T" I,..., do
2-11 Vln² 1)44 2x + (nn² 1) 14-11 1 + 1. (h(n, u._1) h(n. ū.)) E..
v1p² 1) 11
where
h(1, u) arctan
Vu-
2/2 21
)
1 St (ooitnoiay +
112 /
E, & Vortinoizy (1) 1) 1) +
v1 - ai
(172 1)8s +
Output: Tnoisy + Z0-
Denote log(1 + r) as a short-hand. We have
(x) (x)
(z) xxx
(
Therefore it follows that
(x+=²)
VI+7* (0. (1)
V llggyy (x) V. logpx.
where we used the definition s00,-) = V loggxx.)) (.). Plugging the definition = 2 log(1 7) into the above equation yields the desired result.
0
18

---

Algorithm 4 Proximal Consistency Sampler PCS(x, C, (adapted from Metropolis-Adjuste Langevin Algorithm [R998]
starting point R. measurements IR" log-likelihood function of the forward model v): proximal parameter v > O.
Input:
Hyperparameter: Langevin stepsize 7, and the number of iterations N. Initialization: 20 I.
Update: for 0.1,... .N - 1 do One step of discretized Langevin: Set e~7/"" and
2 = + (1 ()=-))=- 1 r)v ((=n: + 7V r2
Wn ~ NO, Ta).
This is equivalent to drawing 2 from a distribution with density Q(: kn), where
| - (r= + (1 r)x + 1 (
Q(2112)= (2 (r) exp
2(1 r2)
2) Metropolis adjustment: Compute
exp ([(= (i2)
212² 11 xx1²)
3.
Q inn 2n
(Z
1 1l - 2:11?
1++ 55 Fm.
cxp
20j
and set
Entz"
with probability min(1, g):
Zn+1
Zn
with probability 1 - min moin(1,00 g).
Output: EN-
D.2 Proof of Lemma 2
Proof. We first compute the probability density function of z. Recall that 2 Enoisy. thus applying Bayes rule yields
p=(x)= = (w + Fnoisy)
= =x" = x + Tnaisyl 1x + E = Innoisy)
p" (x + Knoisy )PE (-x)
p" p²(x0 + a nnnisy )PE) (p((-r):
Per + (Fnoisy)
where & TT)) straightforward to compute
2n ||11
(=(-z) (2
therefore
(x) x p=(x=-=ppp (x + T e 2 1|=11
(32)
We proceed to compute the probability density function of 2r. According to (19 it follows that Pz. (a)
Pe- PV
Pe-rz(x') P. PVTeo
zic(1 2i
x x')da'
2(1 2e--e-27) Pz p2(et exp
||x - x²(12) do'
(- 2np l1e a'll exp
pt (nnoisy + +e-x')
|x a'l12 dz',
2(1 - c-2r)
x f p" (x') (xx))exp exp
(- 2n2
I| Znoisy| 12)
exp
2(1 - e
T xnoisy )1? 113
dx',
(33)
where denotes convolution, the penultimate line follows from 332 and the last line follow from the change of variable x' + e-r(x= (x' Fnnoisy.. ). One may exercise some brute force to verify that
1 I|x' - xnoisy|l 112 exp (- (-
exp 2(1 - 22) 1 e
e (x' Eno
19

---

22r ||=||]
c- i/-11p0
= exp
exp
noisy
+
-
217 + e2t - -1)
2(1 e ---2)
172 + 62t
c²r |
( (e Inoisy
c np x
exp
2(12 + e²2
PVI-e-zz
+
11² e
-1
F x
where F is as defined in (22). Plug this back into
(33 we see
e²r |
e $
pz. (2) x exp
p (x' (1'')-1 (1'i e
Tnois
+
e : de
2(122 + eaT 1))
112 e2
1
e |1:11
e--31²
x exp
p¹ p'(e'="(p/11-- P/1------
Tnoisy
+
x/ da
2(17 + 622
172 + e
27 ||:||
e² 1
x exp
2(72 + e2T
Po-tzo
PV
Tnoisy
+
n72 + 2 -
e² |
X exp
px+ e Tnoisy +
11 + e2T -1)
2(7 + e
where the second line applies the change of variable x' + d the integral, the penultimate
(d)
line follows from Po-tao (x') x xe"x') (since to ~ p²), and the last line follows from X; e- To vI - e g
Finally, from the above formula, we obtain
2T |
e 57-
V logpz, v
+ Veloxxx e 3 Inoisy
+
2(122 + e - 1)
112 + e27
2 x
6" 91=
8² 02-21120
+
F, 5
&nnisy +
e2r -1 112 + e -1
112 + e2T -
where we used the definition s(F,-) = V Vlog ()).
E
Discretization via the exponential integrator
E.1
General form of the exponential integrator
Consider a SDE of the form:
dM, (u(r)M, f(T, Mr) VBBB, [0,To) Mo PMo:
where v [0, Tx] + R, f 10, T x Rd + Rdd are deterministic functions, and a constant. Given discretization time points To S ... S Th S Toor naïve way to discretize the SDE is
1, where e; (0, Ia) a standard d-dimensional Gaussian random vector which is independent of Mr. Although this approach is straightforward, it has the drawback that the linear term u(T M, is discretized rather crude. For example, for the OU process where v = 1, f = can be solved analytically as in (2), while the above approach still has a discretization error.
Mr Mr. (u(T))M,, f(T,, M,.) (T++1 - T3) /B TEE.
0,
2, the SDE
A more accurate discretization, known significantly improve the quality of cosre-based generative models, is given by the exponential integrator [7C22] which preserves the linear term and discretizes the SDE
d,, (e(T) + f(T;, Mr.))dr + /dd, 7 [Ta, T41), i= 0, 1......... 1. k.
with initialization M, ~ PM-- On each time interval [T,, T44a), this is simply a linear SDE, which can be explicitly solved by
M (d) (()) + ([CCto--VEaa))) f(T,, Mr) M(T,
VB ( (2((6)-
1/2
Es.
where V is the antiderivative u:
V(r) = [ (F))d..
20

---

Taking r = TH+1 we obtain
(fr (l (((nt))-V((dd
Man
e (V
S(Ti, MTr)
VBc (T441)
Ei.
(34)
which provides an iterative formula to compute M
E.2 Discretization DDS-DDPM
Plug the expression of V log py, Lemma into (66), and use the notation -rev = .72 - we obtain, for 7 [0, that
yrev
dy
log + rre)),
dr + dB,
I reev
2
VI+ rrev Yrev
1 rev
log(1 -lee),
dr + ddB.
v1 + rev
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 70% S S rey 172 which satisfy
1
log(1 + ) = log OstST.
2
2 at
This solves to
7t ág 1
(35)
The requirement that r! < translates to a. 2 It
which yields the following choice of T':
T' max t:0SITT, 172 +1
(36)
Applying the exponential integrator. Now we apply the exponential integrator to discretize the SDE on each time interval rev E [T-1, Tel, 1.... T" follows:
Prev
dy
rev
cont
log(1 T(
dr + dB,,
Vyrev
VT+ +.ey
Prev
dr + dß,
Virev"
v1 + T.ev
Prey T!
dr + dBr
Vrtev"
The SDE can be integrated directly on rev
[Te-11 Ti] (see also 34
with v = 0), yielding
rev
rey
2(V Tt
T(=1
rey a.
+
dB,df
rev
2(V
Trey
Cit
Frey
T Wt:
where we N/ Ia) is independent of yrev Set x4
ÿ rey rlo,
we obtain
a क 2(VT.
VT-1) -
(Vars,) + 11- W, A/ 12),
(37)
which is exactly the update equation in Algorithm 2 E.
except that ot is replaced by the noise estimate
21

---

E.3 Discretization of DDS-DDIM
Plug the expression of sz in Lemma into the probability flow ODE 20 we obtaín
1j2 -1
- mp
-F(r 11
d d
112 + 22 - Z dr +
1 F(r e 7(( noisy +
dr
n1 + 2
rp + 22 -
112 172 + e2
227
to f(r" 1pp200
7 dr
econt
F((tev), e -- noiss
dr,
22rm
112 + 2
where the second line used the definition
Choosing discretization time points. Similar to the derivation in Appendix time points
E.2
we discretize at
T
7/ I S T/p < 112 which obey
1 t= 1=0......,
F(((-) log
(38)
which solves to
11 + 0t
rev
log
(39)
2
(1) + 1)ā+ 1
To make this well-defined, we require
112 + αe 1 0,
(17² + 1)a;
which is equivalent
αt >
1 1+ 1/2
This leads to the same choice of T as in 36). We also set
T
It is convenient to introduce notation for the corresponding discrete schedule of ,(" denoted by
(1)2 + 1)ā. 172 + ā
Úe
= 2T me
1 1= 1,0,1....,
Applying the exponential integrator. Now we apply the exponential integrator, which discretizes the ODE on each time interval rrev
[,],,,7],1 1,... T",
1
2
e m0 - -((r("), 2
rrev
2rev Zq
con't
F(T("v)
F(rt
Tnoisy
+
dr
172 + 2r
172 + 2 -
112 2-1
m
/atl) 7rev 1/2
Frev dr
v Tt
2: -1 dr,
112 + g2r"m
-1
22rm
112
2
rev
+ 2
Zrex dr
22
E
(17 1)ux,
+
dr,
where the second line follows from and the numerator in the fraction inside E, by e²r?". This is a first-order linear ODE on rev E [T) which be solved explicitly (cf. (34) by
38
and the last line follows from dividing both the denominator
(np² - 1)e-2 2 +1
Grev tytorr
2rev
Vn1² - 1)44 +1
+ (T² -iee-2t 27 +1. (h(n, -h(n, -h(,, -et (Varenotey + Vuivan,1 Zrev
(172 1)ū, +
for rov E [T--1. Te], where
h(1, u) == arctan
22

---

Plug in 7Te TT-1 in the above solution, and set 2 Zaze, we obtain
Vnn (1-1) 1)4411 2t
2--1
(40)
V(1² 1) +1
V11 11)u-1 1
(h(n)ū=-1) 1(1, ū+)) ot Vo (Voar)nni
Vuz ax
(172 1) 77 +
The initialization, which should ideally be 2r" Zree pz, approximated by 24" N No, Ia). This is exactly the update equation and the initialization in Algorithm except that E, is replaced by the noise estimate ,.. E.4
Discretization of PCS
We first note that the Metropolis- Metropolis-adjustmenr 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-timg time Langevin dynamics for sampling from the distribution exp (xp(C(-i))
- 3):
27
1 dz, v, ((z:))Z) + ,,2 (z, 122 z)dr VdB,, 720, 0, Zo N/ Ia).
(41)
The classical form of MALA, in [R998] performs step of straightforward discretization of as the Langevin step, as follows:
2 ~ Zn > (( + 172
(zn x) +1 22 Wn ~ NO
In our setting, due the presence of the linear drift term + ( I), which can be quite large when 7 is small, apply the exponential integrator instead. Set the discretization time points Tn nay the exponential integrator reads as
dz, - 2n. ((Zn y 1 (Z, 2)dr + V2dB,, ny ST (n+ 1)7.
1/2
Solve this linear SDE on my < T < (n 1)y directly (see also to obtain
(n+1)y rny + (1 r)x + n/2 r) Z.., L( ((Z v) + r
Wn, N(0, Ia),
where = = 7/1² This is the same the update equation for the Langevin step in PCS (cf. Algo- rithm
F Proof of main theorems
F. Proof Theorem
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 (Correctness of PCS). Under Assumption 1 with notation in Algorithm f in the continuous-time limit: ? N - oo, the algorithm PCS outputs samples with distribution x exp(C(:iv) + 2 12)
The next lemma guarantees the correctness of DDS with exact unconditional score functions.
Lemma 4 (Correctness of DDS). Assume the score function estimation 3, is accurate, i.e. 3, s In & - unifomly in t, both DDS-DDIM and DDS- the continuous limit: T - co, oer 0, à.
DDPM output x obeying the posterior distribution p" (x* x x² + TE = xnoisy). ) ~ N La).
The proof of Theorem 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. Denote
(x) = Purmp" soM/0,,0 + mE = x) (2m))/227d / (x)) - 1|= dz.
23

---

From the first equality, it is clear that Py + p" when 17 + 0+. We will also use the notation In defined in (24). which we recall here:
(x(()) = (2x)a/211p) / ((xiv)) 20 =1 de.
In virtue of the Assumption we know that 9, is finite for all E Rd.
For convenience, we introduce a notation for application transition kernels. For probability distribution p(z) and probability transition kernel K (z, x'), denote by K the probability distribution given by
po K(xx)=)==)) pl=)
The first lemma characterizes the one-step behavior of DPnP in terms of Markov transition kernlss
Lemma 5. Under the settings of Lemma 4 and Lemma the one-step transition kernel of DPnP
with Ti = 1) is given by:
KDPnP (n (2,,²)) (ʃ 90(2) (z)
27 1|
11 2n 1==111111 dz
(x')
4n (1(x)
In other words, Ty has distribution pay, then the distribution of Tx+1 14+10. is
Pa. (x') = Pua KDPpP,, = f Pag (x) K (2, x' dr.
The proof is postponed Appendix E.3 Markov chain with transition kernel DPPP... These properties are known [B233] but scattered in different literatures, will provide brief proof to be self-contained
The next lemma analyzes the ergodic properties of the
Lemma 6. The Markov transition kernel KDPnP,,, has the the following properties:
(i) (Stationary distribution.) Let Ty be the probability distribution defined by
(T((x) C,,p" (z))
where On > the normalization constant such that J (r(x)dd 1. Then KDPnP.,, is reversible with stationary distribution Te.
(ii) (Convergence.) For any initial distribution p, the distribution of the Markov chain with kernel KDPnP., converges to q.
TV(p o k<(n) p Tr) + 0, n + xc,
(42)
{(n)) where DPnP.n is the n-step transition kernel of KDPaP..
The proof is postponed Appendix F.4
We show how to prove Theorem 1 with the above two lemmas. With the annealing schedule in Theorem [ between steps k < k < ki which consist of consecutive (k: ki-1) steps, the transition kernel of one-step of DPnP is KDPoP.s. As (ki - 1)-1) + co, Lemma implies that
TV Te,) = TV(Pa,,
K- (
DPnP.51
ne,) - 0.
Under the assumption in Theorem thate: e + 0, we let oo to see lime ne, = cop" (-)e p' (-ly), thus pax, p" (:ly), as claimed.
24

---

F.2 Proof of Lemma
Proof. For DDS-DDPM, we note that under continuous-time limit in Lemma the discretization time points given by 35 verify
0
0,
sup 0S
Tev 2++1
+ 0,
rev TT"
1 11 T - co.
Therefore, these discretization time points 0 == TD .... TTp 17² form partition of [0, n7
which becomes infinitely fine in the continuous-time limit. Thus the discretized integrator (37 converges to the solution of the SDE (I6). which, samples obeying the denoising posterior distribution p" (-|Enoisy). as claimed.
have already argued in Appendix D produces
The proof for DDS-DDPM follows similarly, by observing that the discretization time points in (35 form an infinitely fine partition of [0, oo) in the continuous-time limit.
F.3 Proof of Lemma
Proof. The proof is based on computing the transition kernel of the two subroutines. We claim that
(i) Sampling with probability density proportional exp(A(;v) 27 ll -=]l?² is equivalent applying the following Markov transition kernel
1
C(a w) 2/. || =x1l
(I,
9n (x)
Sampling with probability p" (a*1 xx* + TE = x), where La), is equivalent to applying the following Markov transition kernel:
(x) p' (x')e ス|/ 121
1
(x,s,,(xx2'')
It is then clear that
KDPpP (2r, x') = f KDDDs,
= ( g0(2)
e 20 110- 1- 2., =1==0'11
p (x²)
dz
(z)
(n(x))
as desired. We now prove the above two claims. For (i), note that by 23 we know K DDDS, (2, p=(-)e -111 Thus it suffices compute the normalization constant, which is
f (=())o zhe |s =I1 dr' =P,(
by the definition of Pty Therefore
(o (x')e zńं lo aal
KDDS x') (r(())
as claimed. The proof of (ii) follows similarly.
F.4 Proof of Lemma 6
Proof. We first introduce fundamental lemma PW224] which provides simple method to bound the total variation between two distributions.
Lemma 7 processro inequality). Let q be two probability distributions, and K be
probability transition kernel. Then
TV(p - K, K) & TV(p.q).
We now prove the two items in Lemma separately.
Proof of (i). We first show that n, well-defined, i.e., S (2))))))))) < co. This can be seen from Assumption which implies (x) $ fe + 1|1-=11 dz $ 1, hence
f (1())), ('(x)dx= 1
25

---

To show that KDPpP,, is reversible with stationary distribution Kyp it suffices to verify (=) (n x') n (x') KDPnP, (2', x), Vz, x' e R..
However, it is easily checked that both sides are equal to
Gn ( u0(2) a =|=-==|+ a 1=='x dz p¹ (x')p"(x) (2)
Proof of (ii). We define an auxiliary Markov transition kernel Kaux,n = KDDS. a KPcs. More explicitly,
(x, (DDS..( (t KPCs.,
(frite (2(2)
2 11=-==11= 212 || x/
dz
Pr (2))
(43)
It is easy to see that
po ,PnP,, = o KPCS, K(m=1) aux.. KDDS..
K(n)
(44)
Thus are led investigate the ergodic properties of Kaus,n- Similar the proof of item (i) above, is not hard to show that K Kauxn is reversible with respect to the stationary distribution
M = (,),(2))00e CyPr (x)
Moreover, one may check that
Th = un KDDS.
(45)
is apparent that u(x') > 0 and 1]. such a Markov transition kernel obeys, for any probability distribution q, that
for all 2, e Rd, By [Tie94] Corollary
TV(g K(n) aismm HG))
0,
In view of 44 and 45
we set q p KPCS,,, and invoke the data-processing inequality to obtain
TV(p [(n)) Dp Tn) TV(g K(m=)) KODS My KDDS.")
TV(g K(a-1), Mn)
0,
00. This completes the proof.
F5
Proof of Theorem 2
Proof. Denote by KeCS., and KDDS,, and the transition kernels for PCS and for DDS, respectively, Note that these may deviate from the transition kernels Kpos,,, and KODS,,, defined for the idcalized asymptotic setting in Appendix [
We have
TV(pa TT,) = TV(Ps x- KDDS..: TT,)
< TV(psx=> KDDS. π.) + TV (1pax_3" KDDS Pa a KDDS.,)
TV(Pa KDDS. Tr)) EDDS.
where the second line is triangle inequality, and the third line follows from the assumption in Theorem 2 that DDS has error at most EDDS total variation, by taking the input DDS to be -
Similarly, from P2,- P o KpCS,,, and the assumption that PCS has error at most EPCS in total variation, we can show
TV(pa, 1 .Kppss. mT) S T KoDs.n x,) +EPCS -KPpe., tECCS
The above two inequalities together imply TV T.y) S TV(pav-1 KDPnP.n: nT) + EDDS + EPCS-
Iterating this process, obtain
TV(pa n S TV(pa, o DPnP.n mn) + (N - () + EPCS).
(N-1)
(46)
(N1 DPnP,n Tn). For this, need the following two lemmas.
It remains to bound TV(pa,
26

---

Lemma 8 (Comparing TV and -iveergeeee,, PW24 ) For any two distributions we have
TV(p, S V (p
Lemma 9 ( cycttttrryyy of KDPmP..). There exists some À = X(p", ,,,)) E (0, 1), such that for any probability distribution p(a:), we have
x (2 o K DPmP,, in
|| m.) ANN x2(11 m.))
A form of Lemma 9 is well-known for Markov chains with countable state spaces, but relatively few sources provide complete proof for the abstract setting we consider here with continuous state space. For sake of completeness, we prove Lemma in Appendix F.6
Combining the above two lemmas, obtain
TV(pa,
K (N-1)
DPnP, TT) V x2 DPPn. mT) S AN-1 Xx (q)
Plug this into 46 obtain
TV(p2x. TT) S N-1 V X2 (T)) + (N ()DDDS + EPOS).
With N x (ogea/) 1-X such that N S exp ( (N - -1) - ))) - Eace, the desired result readily
follows.
F.6 Proof of Lemma
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 (=x to be the space of square-integrable functions with respect to measure π, i.e., those functions f Rdd - C such that
|l = ( ((1)))2))))
1/2
80.
The first well-known property (SC97 offers way to represent reversible transition kernel self-adjoint operator (infinite-dine dimensional symmetric matrix).
Lemma 10 (Self-adjoint representation of reversible Markov operator). Assume K K(x, transition kernel that is reversible with respect to the stationary distribution π T(x). Then the integral operator K. (T) L22 (T) defined by
(x,xx') x²) is Markov
Kf(z)= ((,,)=
self-adjoint and compact. For any probability distribution p(x) such that f
p(() x(x)
dx < oo, we have
f p(x) (f(x)dd== = f K(z))(z')dd'
Moreover, the eigenvalues of K. are the same as those of K.
The following theorem is generalization of the classical Perron-Frobenius theory for finite- dimensional transition matrix to strictly positive operators. The form we present here can found in [2ch12 Theorem V.6.6]; see also [Bou23b] Theorem I1I.6.7] for more elementary treatment which can also be adapted to the form we need.
Theorem 3 (Jentzsch). Let (x,xx) be Markov transition kernel. If K(xxx') for any x. E Rdd, then K has a unique stationary distribution w. Moreover, 1 is simple eigenvalue of K, with being the only left eigenfunction, and the constant fanction being the only right eigenfunction. In addition, there exists X (0, 1) such that any other eigenvalue of K. has modulus no larger than X.
27

---

We are now ready to prove Lemma 9 We divide the proof into the following steps.
Step I: controlling the eigenvalues of KDPPPP,, Recall the auxiliary kernel Kaux.n defined in (43) It a standard result in linear algebra or function analysis [Bou23a| that Kaux, = Kpcs. DPPnP.,
KoDs., KPcs. From (43) it is easy check K Ka (T, x') > 0, thus Theorem implies 1 is simple eigenvalue of KDPnP,.. Moreover, there exists Mpp,C.n) (0, 1), such that any other eigenvalue of Kaux,n has modulus no larger than X.
pDDS,) has same eigenvalues
Since KDPnP.,, has the same eigenvalues Kaux.y and, by Lemma 10 the operator KDDpnP,, also has the same eigenvalues as these two, we conclude that KDPnP.,, is a self-adjoint compact operator on T,2 (TTy). 1 of whom is a simple eigenvalue. Moreover, any other eigenvalue of KDpnn,, has modulus no larger than X.
Step 2: establishing the contractivity of KDPnP.,, in (T,,,, is easy to verify that the constant function 1, which takes value for any x E RRd is eigenfunction of KDPnP., 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 ((]]) = (f T., (z)dx))// = 1. Therefore, the operator KDPaP,,, 11 is self-adjoint operator whose eigenvalues have modulus no larger than X, where 11 T is the orthogonal projection onto 1 in 1,2 (T,)) defined by
11 11 f(x) = f(x'())))))))) (zz)
Vo EEER
Using the fact that C, = 11 KDPnP. 11", one may show (KDPnP,n - 11 117 x 11 by expanding the product, e.g. [SC97). Consequently, KD 11 is self-adjoint operator whose eigenvalues have modulus no larger than XN, ie..
|KDPn.., - 11 (1,(.)))1((T.) S SA,,
(47)
where || (lLP((,))+LP(T denotes the operator norm on L² (T(TT.).
(NN) Step 3: bounding the inner product of o D - my with any square-interrable function. Note that when xx |I n.)) = oc, the conclusion is trivially true. For the rest part of the proof, we assume x² ).,) < co. Now, for any f E L (TT), by applying Lemma 10 iteratively, we obtain
(N) f p Kppn,, (x) (x)((x)r
/ p(x' (an'))oPPP, on (z')dx'
/ p(x'). . (KDPnP.n
11 11 ff + p(x))11 f(x')dd"
/ p(x'). (KDPpp,, 11 T of("')d'' + (a))))(("))dt)
where the last line follows from the definition of 11 T and p(pa)dd = Rearrange the terms to see
/
(N) (poKkPAn.,, (x) - (x) (=))) = / p(z') (KDPnP.n - 11 T ) ((a')d''.
(48)
In partieular, taking p= n. yields
0 / (x()")) (K 11 )f(''dd''
(49)
Substract 49
from 48 and then take absolute value, obtain
/f (p K-OMP, (x) - xx,())) (f (z))
// (p(=') - n. (x')) (KGPpp..n
11 ))(a')da'
x (/ (p(x') - n,((x') dx
1/2
||(xCDDmP.n - 11 (TTT)((())) )f(2')llR())))
(x)
28

---

(50)
Step choosing appropriate square-integrable function. Now, set
f(x) = po po KDnnP., (x) (x)
(x))
is easily checked that
I (p oPnp,,, (x) mT(x)) = x² KDPp.,, || Tn):
[1) v ( K(N)
DPnP,n mT))
Plug these equations into 50) we obtain
2p K DPnP,,, || m,) S x2 m,))
as claimed.
G Additional numerical results
G.1 Implementation details
Score functions used. We use the same pre-trained score functions as in CKM+23 3
Normalization. All images normalized in the usual way to fit into the range [-1, 1].
Parameters of our algorithm. We choose the same annealing schedule across all tasks. Please see Appendix for detailed discussion.
Parameters of comparison methods. We made our best effort to fine-tune the other algorithms within a reasonable amount of time for each task. We list the paramters, following the notations in the original paper [CKM+23 SZY+233, as follows.
Super-resolution: For DPS, the learning rate is set to 0.6. For LGD-MC, the MC sampling variance Л; 0.05, the loss coefficient X 10-3 and the learning rate is 60.0. For ReSample, the stochastic resampling variance parameter = 10 for FFHQ and y = 4.0 for ImageNet.
Phase retrieval: For DPS, the stepsize is set to 0.8. For LGD-MC, the MC sampling variance rt = 0.05, the loss coefficient X = 10-3 and the learning rate is set to 400.0.
Quantized sensing: For DPS, the stepsize is set to 100.0. For LGD-MC, the MC sampling 10-5 and the learning rate set 500.0. For ReSample, the stochastic resampling variance parameter y = 4 for FFHQ and 3.5 for ImageNet.
variance 7; 0.05, the loss coefficient À
G.2 Forward measurement operators
Super-resolution. The forward model for super-resolution is the usual bicubic downsampling operator [Key8]], which is a linear operator (in fact, a block Hankel matrix). We use a downsampling ratio of 4 in all our experiments. The measurement noise is set to be white Gaussian, with variance 0.2. Note that the noise variance moderately larger than that in [CM 23 to better reflect the scenario in practical inverse problems.
Phase retrieval. We consider phase retrieval with coded mask, which a classical inverse problem (CLS15). For 256 x 256 image (for each color channel) in our experiments, we first generate random mask E 1226x256 (which is shared across color channels), then apply Fourier transform F to M @ where @ denotes the Hadarmard (entrywise) product, and finally preserve only the magnitudes of the Fourier transform. Formally, the forward measurement operator is A(x) = mag(F(M M @x)), where mag(-) computes the entrywise magnitude of matrix with complex entries. The measurement noise is again set to be white Gaussian, with variance 0.2.
29

---

Table 5: Samples of different algorithms for super-resolution (4x).
DPnP. DDPM
DPnP- DDIM
Ground
Input
DPS
LGD-MC
ReSample
truth
Quantized sensing. this work, quantized sensing refers to the task of reconstructing an image from its low-bit quantized version. The forward measurement operator is one-bit per channel, dithered quantization operator. More precisely, the forward measurement operator an entrywise application of the following stochastic function Q with dithering level 0 >> 0:
(pixel) = -1,, with probability
14
1-1, with probability 14
where pixel E [- is the value of each pixel in each channel. The measurements quantized sensing therefore mue-bitpper-hanne images. The dithering level is set to 0.4 in our experiments. G.3 Sample images for other inverse problems Super-resolution. The samples generated by different algorithms are shown in Table
Across different tasks, linear and nonlinear, can be seen that DPnP has stronger capability of reconstructing the image with higher fidelity the fine details.
G.4 Additional performance metrics
Computation time in terms of Neural Function Estimations (NFEs). additional to the clock time statistics in the main text, also measure the computational cost per sample of different algorithms in terms of the number of Neural Function Estimations, i.c., the number of calls to score functions. The results are in Table 6 Note that the NFEs for DPnP depends the initialization, the annealing schedule (and the number of timesteps for (PnM-D We provide typical numbers of NFEs with the choice of parameters given in Appendix H and with suitable number of timesteps.
Table 6: Number of NEEs for different algorithms.
Algorithm
DPnP-DDIM DPnP-DDPM
DPS
LGD-MC
ReSample
NFEs
1500
3000
1000
1000
1000
Frechet Inception Distance (FID) and Structural Similarity Index Measure (SSIM). We also compare the FID and SSIM of different algorithms across different tasks. The results are shown in Table 7 and Table should be pointed out that FID arguably not very relevant notion to measure quality of solving inverse problems, as accurately solving inverse problems means that the generated distribution is close to conditional, i.e., posterior distribution of the image, while FID only measure the closeness to the unconditional, i.c., prior distribution of the image.
30

---

Table 7: FID and SSIM of solving inverse problems on FFHQ 256 x 256 validation dataset (1k samples).
Super-resolution (4x, linear)
Phase retrieval
Quantized sensing
(nonlinear)
(nonlinear)
Algorithm
FID
SSIM
FID L
SSIM
FID
SSIM ?
DPnP-DDIM (ours)
36.3
0.668
46.5
0.631
37.3
0.712
DPS [CKM 23
38.6
0.636
52.0
0.494
42.1
0.601
LGD-MC (n 5) SSZY=23 23
36.8
0.651
82.3
0.414
40.3
0.639
ReSample (pixel-based) |SK 23
40.2
0.641
40.0
0.657
Table 8: FID and SSIM of solving inverse problems on ImageNet 256 x 256 validation dataset (Ik samples).
Super-resolution (4x, linear)
Phase retrieval (nonlinear)
Quantized sensing
(nonlinear)
Algorithm
FID 4
SSIM
FID +
SSIM
FID 4
SSIM f
DPnP-DDIM (ours) DPS C3
47.5
0.510
73.5
0.289
43.2
0.623
61.4
0.496
92.7
0.318
82.4
0.459
LGD-MC (n 5) SZY+23
46.2
0.503
89.8
0.234
46.8
0.563
ReSample (pixel-based) 23
68.3
0.427
53.7
0.491
H
Ablation studies
H.1 Initialization
In the Algorithm T1 the initial guess Ty is set be properly scaled Gaussian random vector. However, from Theorem it can be infered that using a heuristic posterior sampler as the initializer could decrease x² (pa. T, hence potentially improve the convergence speed 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. our experiments, we find it helpful to initialize DPnP with LGD-MC, which accelerates the algorithm significantly.
H.2 Annealing schedule
We diseuss the choice of the annealing schedule mes in DPnP (Algorithm 1. As seen in the theoretical analysis (Theorem 2 if we set all the 7/2 = r for some constant >> O. then DPnP converges to distribution Ty. which can be regarded as version the posterior distribution p" (-ly) distorted by
an order of (n)) The smaller 1j is, the more accurate the final distribution will be. On the other hand, it was also seen that in many cases, the spectral gap (()). hence the convergence time is O(2).
Therefore, smaller 7) would make it take longer to converge
To strike balance between the accuracy and the convergence rate, find it empirically successful to adapt an gradually decreasing schedule for 7ą. similar to [BB23. In the first few iterations, we set mn be large constant. After this initial phase, we decrease Mi slowly, eventually 7ι which is chosen to be a small constant. An example of such an annealing schedule is
TJu = 171 = = 7Ka
170 > 0 is a large constant,
b- Ky
7] (n no)
K K- 1)0.
Ko < k S K,
ЧパK 0 a small constant,
where Ko < K is the length of the initial phase, which can be chosen as, e.g. Ko K/5. For all the numerical experiments, we set fjo 0.4, y/n 0.15, Ko 4, 20.
"Strictly speaking, while the number of iterations required to converge increases as y gets smaller, the computational complexity 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.
31

---

NeurIPS Paper Checklist
1. Claims
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.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made 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 is clear that these goals attained by the paper.
2. Limitations
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 limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors encouraged to create separate "Limitations" section in their paper.
The paper should point any strong assumptions and how robust the results are 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 few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails 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 faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, 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 impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
32

---

Justification:
We provide clear statement of assumptions for our theory.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We provide algorithm tables and experimental setups with a reasonable level of detail.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
33

---

Answer:
[No]
Justification: The code is still preparation to meet publication standard. We promise to release it after this work is accepted.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips... cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to to reproduce the results. See the NeurlPS code and data submission guidelines (https
ssicnpolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Experimental settings are given in detail. Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
Justification: The computational cost of reporting error bars is too high. Most of the closely related works do not include them.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
34

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: All experiments done a single Nvidia L40 GPU, The time of execution is reported.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: We fully honor the NeurIPS Code of Ethics this work.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration 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-basse methods in solving inverse problems, it awaits future research before anything special can 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.
35

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is 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 being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there 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 monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
No data or models is released by this work.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 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 best faith effort.
12. Lieenses for existing assets
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?
Question:
Answer: [Yes]
Justification: All assets credited properly.
Guidelines:
The answer NA means that the paper does 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 URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.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, paperswi thcode. com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of 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.
36

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets,
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
37


Reviews:
[{"review_id": "6c4zpw6Eu8", "reviewer": "Reviewer_43oG", "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.", "questions": "Please see the weakness section above.", "limitations": "Yes, the authors have addressed the limitations.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. This paper provides a robust algorithm for solving nonlinear inverse problems using unconditional score-based diffusion priors. \n2. The paper is nicely written and theoretical analysis in a simple setting clearly demonstrates the major contributions of the paper.", "weaknesses": "### **Weaknesses and comments**\n\n1. Line 171: \"Assumption on forward model is applicable to *many* applications of interest\" what are the applications?\n2. Assumption 1: what are some examples of $\\mathcal{L(\\cdot, y)}$ that are differentiable almost everywhere and used in practice?\n3. Step 1 in **A stochastic  DDPM-type sampler via heat flow**: How do you know the noise level $\\eta$? This is typically unknown.\n4. 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. \n5. 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?\n6. 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).\n7. 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?\n8. 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.\n9. 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.\n10. The authors are encouraged to cite the published version of the papers where applicable (see for instance CCL+22).\n\n\nReferences\n\n1. Solving Linear Inverse Problems Provably via Posterior Sampling with Latent Diffusion Models.\n2. Prompt-tuning latent diffusion models for inverse problems\n3. Beyond First-Order Tweedie: Solving Inverse Problems using Latent Diffusion\n4. Tweedie Moment Projected Diffusions for Inverse Problems", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review and valuable comments. Below we address your comments in a point-to-point manner. \n\n**Examples for Assumption 1 on the forward model and the likelihood (point 1,2).**\n- A typical scenario is where the measurement noise $\\xi$ in the measurement model $y=\\mathcal{A}(x^\\star)+\\xi$ is Gaussian, e.g., $\\xi\\sim\\mathcal N(0,\\sigma^2I_m)$. In such case, we have$$\\mathcal L(x;y)=\\log p(y|x^\\star=x)=-\\frac1{2\\sigma^2}\\Vert y-\\mathcal A(x)\\Vert^2 - \\frac{m}2\\log(2\\pi\\sigma^2).$$\nTherefore Assumption 1 holds for *all* applications where the forward model $\\mathcal A$ is almost everywhere differentiable, in particular, if $\\mathcal A$ is linear. This includes inpainting, deblurring, super-resolution, phase retrieval, etc. It is noteworthy that most of the existing works (e.g. [CKM+23, SKZ+23, SZY+23, DS24, MSKV24], and [1-4] you listed), have the same differentiability assumption or requires the stronger condition that $\\mathcal A$ is linear. \n\n**The noise level $\\eta$ (point 3).**\n- The noise level $\\eta$ in the DPnP framework is a hyperparameter used for annealing. It is not related to, and should not be confused with, the *measurement noise level* (the variance of $\\xi$ in the measurement model $y=\\mathcal A(x^\\star)+\\xi$). The annealing noise level at the $k$-th iteration, denoted by $\\eta_k$ in Algorithm 1, is set manually, hence is always known.\n\n**How to know the likelihood term (point 4).**\n- As the example above (please refer to our response to point 1 and 2) demonstrates, the likelihood can be computed if the forward model $\\mathcal{A}$ and the noise distribution is known. This setting is adopted by most of the existing works (please refer to our response to point 1,2), including [1-4] you listed (in fact, most of them assumed $\\mathcal{L}(x;y)=-\\frac1{2\\sigma^2}\\Vert y-\\mathcal A(x)\\Vert^2 + \\rm const$, which is a special case of ours with Gaussian measurement noise). It might also be helpful to note that our log-likelihood function $\\mathcal L$ is fully, simply determined by the measurement model, and should not be confused with $\\log p(y|x_t)$, which was also called likelihood by some of the previous works.\n\n**Computation cost of inverting $\\Sigma$ (point 5).**\n- While most of the existing works (including [1-4]) simply assumed (implicitly or explicitly) $\\Sigma=\\frac1{\\sigma^2}I$ is scalar (please refer to our response to point 4), it is possible to handle general $\\Sigma$ with standard numerical tricks in our work. For example, one can pre-compute, *once and for all*, the vector $A^\\top\\Sigma^{-1}y$ (which amounts to solving a linear equation) and a prefactorizization of the PSD matrix $A^\\top\\Sigma^{-1}A$. Using such precomputed information, each iteration only needs to compute several matrix-vector products, which has a negligible computation cost compared to Neural Function Evaluations. The cost of computing the prefactorization of $A^\\top\\Sigma^{-1} A$ has also been observed to be acceptable [SVMK22, KEES22]. Another option is to use MALA for the PCS step as we did for general non-linear inverse problems, so that there is no need to invert any matrix.\n\n**Discretization error and TV bound for subsamplers (point 6,7).**\n- As Theorem 2 indicates, discretization error affects our algorithm only through total variation error of subsamplers DDS and PCS, so we respond to point 6 and 7 together. The reason we did not include an analysis of TV error of subsamplers is, as explained in the paper, that both subsamplers can be analyzed with existing techniques in the references we provided. For example, applying the results in [LWCC23], one may show that DDS-DDPM has total variation error $$\\varepsilon_{\\sf DDS-DDPM}\\le \\tilde C\\frac{d^2}{\\sqrt{T'}}+\\tilde C\\sqrt d\\varepsilon_{\\sf score},$$ where $\\tilde C$ hides logarithmic dependence on parameters, $T'$ is the number of steps (cf. Algorithm 2), $\\varepsilon_{\\sf score}$ is score estimation error defined in [LWCC23]. Similar bound for DDS-DDIM can also be proved using the results therein. On the other hand, utilizing results on the convergence of MALA [CLA+21], one may prove that PCS has total variation error $$\\varepsilon_{\\sf PCS}\\le C\\exp(-C\\sqrt{N}\\eta^{-1/6}d^{-1/4}).$$\n\n**Comparison with works other than DPS (point 8).**\n- Thank you for the references. As we mentioned in General Response to All Reviewers, we compare our algorithm not only with DPS, but **also with more recent, highly competitive** algorithm like LGD-MC [SZY+23], and plan to add comparison with ReSample (pre version available in the rebuttal pdf). Note that comparison with such stronger baseline as we have done is actually not too common in the literature, given the highly active status of the field and the well-established position of DPS, e.g. in all the above references and [1-4] you listed, the strongest non-latent baseline is DPS (latent models are discussed in the next paragraph). In addition, none of these previous works come with a provable robustness guarantee, thus our algorithm is also of independent theoretical interest.\n- Concerning [1-4], all of them except [2] are confined to linear inverse problems, while our focus is on general, non-linear inverse problems (e.g. phase retrieval and quantized sensing). [1-3] rely on the use of latent diffusion models ([2] further requires prompting the model) which is again a narrower setting than ours as discussed in General Response to All Reviewers. Unavailability of code also hinders comparison with [2-4].\n \n**Whether the dataset is the same with DPS (point 9).**\n- We have strictly followed the original DPS paper on the usage of dataset. We used the same dataset without picking smaller subset than DPS. Our experiments are done for different inverse problems, though, as our focus is on more non-linear problems. In doing so, we made our best effort to ensure fair comparison by fine-tuning all competing algorithms and listing explicitly our problem parameters in the paper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "t415mmmJab", "reviewer": "Reviewer_tiqf", "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.", "questions": "- Can the findings be extended to latent domain samplers? This would greatly improve the efficiency of the technique.\n- How does the method compare to other samplers such as ReSample?\n- How does performance scale with NFEs?\n- I would recommend changing DDS to some other abbreviation to avoid confusion with the Decomposed Diffusion Sampling method [2].\n- What does G denote in line 136?\n\n[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 are not clearly addressed in the paper.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The theoretical analysis is a valuable contribution. A lack of robustness to sampling errors and the resulting error accumulation has been a key challenge of diffusion-based solvers, especially in highly nonlinear tasks such as phase retrieval.\n\n- The paper is well-written overall and the structure is logical. \n\n- The experimental results are promising. In particular the proposed plug-and-play sampler achieves significant improvement over DPS, a well-established technique in the literature.", "weaknesses": "- 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.\n- 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. \n\n[1] Song, Bowen, et al. \"Solving Inverse Problems with Latent Diffusion Models via Hard Data Consistency.\" The Twelfth International Conference on Learning Representations.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review and insightful suggestions. Below we address your concerns in a point-to-point manner.\n\n**Comparison with ReSample [SKZ+23].**\n- Thank you for your suggestion. There are a few recent algorithms that included a stochastic correction step for better robustness. The LGD-MC algorithm [SZY+23], which we compared against in our paper, is also one of them. We chose it due to its simplicity for implementation, but we are also more than happy to include comparison with the non-latent version of ReSample in the final version; comparison of latent versions, on the other hand, requires developing a latent version of our algorithm, which we feel more appropriate to leave to future work as discussed in General Response to All Reviewers. There are some preliminary comparisons in General Response to All Reviewers and in the rebuttal pdf, and we plan to include a more detailed comparison in the final version.\n\n**Computational cost.**\n- Thank you for your suggestion. In our paper, \"~3000 NFEs\" is the computational cost of the DDPM version of DPnP; we also proposed a DDIM version in the paper, which costs ~1500 NFEs, only 1.5x the cost of DPS (1000 NFEs). Moreover, the DDIM version allows to incorporate acceleration techniques for unconditional DDIM sampling, e.g. DPMSolver++, which can reduce the computation significantly. We would be more than happy to add more discussions in the final version.\n\n**Performance scaling with NFEs.**\n- Thank you for raising this good point. We have included a few samples with different NFEs in the rebuttal, Figure 3, and are willing to investigate it more systematically in the final version if possible.\n\n**The name DDS for Denoising Diffusion Sampler.**\n- Thank you for your suggestion. We will change it to an unambiguous name in the final version.\n\n**Notation $G$ in line 136.**\n- Thank you for your very careful reading and sharp observation. This is a typo for $M$, the forward diffusion process. We will proofread the paper carefully and make corrections in the final version.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "XCoVHsc39m", "reviewer": "Reviewer_ppwi", "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.", "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. \n\nCould 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": "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.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper establishes both asymptotic and non-asymptotic performance guarantees for DPnP and provides numerical experiments to demonstrate its effectiveness across various tasks. The theoretical analysis presented is a valuable contribution to the field.", "weaknesses": "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. \n\nMinor 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.\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive feedback on our theoretical contribution and for your insightful comments. It seems that the reference items [1-2] in the review are unfortunately missing, so this rebuttal will be based on our understanding of general related literatures. We would really appreciate it and would be happy to discuss in more details if you could provide these references during discussion period.\n\n**Distinction from previous work and new insights.**\n- All of the previous works that used both score-based model and plug-and play, e.g. [RS18, LBA+22, BB23, CDC23]  are more or less heuristic and do not have provable consistency or robustness. While some of them realized the potential of split Gibbs sampling as a plug-and-play posterior sampler, none of them found a way to solve the denoising step rigorously. One of our new insights is to show that the denoising step can be solved in full rigor with a reverse diffusion process using only unconditional score function, and it is even possible to make this diffusion process a DDIM one which enables application of many acceleration techniques. Moreover, the gain of our denoising diffusion step is not only mathematical rigor, but also improved robustness, since by running a full, multiple-step reverse diffusion process instead of a one-step heuristic denoiser as in previous works, we are able to better compensate for the drifts from the data manifold induced by the consistency step. It is with this insight that we are able to establish the first provably consistent and robust diffusion posterior sampling algorithm. \n\n**Demonstration of convergence.**\n- To further demonstrate the convergence of our algorithm as shown by our theory, we run a proof-of-concept experiment with data generated from a Gaussian mixture model, so that the convergence of distribution can be examined directly. The detailed settings can be found in General Response to All Reviewers, and the results can be found in the rebuttal pdf. It can be clearly verified from these results that our algorithm does converge to the true posterior distribution. \n\n**Additional noise in DDS-DDPM.**\n- 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, DDS-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].\n\n**Typos.**\n- Thank you for your careful reading. We will proofread our paper more carefully and correct this typo among a few others in the final version.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "EBN5e2MHK7", "reviewer": "Reviewer_juip", "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.", "questions": "See Weaknesses.", "limitations": "The authors adequately discussed the limitations of their method.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The paper is well-written and easy to follow/understand (barring some notation issues, see Weaknesses).\n2. Each step in the process of developing the DPnP algorithm is well-motivated and explained thoroughly.\n3. The provided proofs for the main theorems in Appendix F are well-written and seem correct.\n4. DPnP outperforms competitors in non-linear inverse problems in nearly all tested metrics in Section 5 and Appendix G.\n5. The overall contribution is impactful with respect to non-linear inverse problems.", "weaknesses": "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).\n2. 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.\n3. 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.\n4. 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.\n5. 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive feedback and insightful comments. We will make sure to take your invaluable suggestions into account in revising the paper and in our future work. Below we address the two questions you raise.\n\n**Connections with projection-type approaches.**\n- Thank you for your thoughtful question. We are willing to clarify our connection and difference with projection-type approaches as follows. \n> While our algorithm draws inspiration from projection-type approaches, the essential difference between our algorithm and these approaches persists even for linear problems. Roughly speaking, our proximal consistency step can be reduced (cf. Line 274-277 in our paper) to a slightly more complicated, yet theoretically rigorous version of the (proximal) projection step in these approaches for linear problems, but the same _cannot_ be said about our denoising diffusion step. \n> \n> A manifestation of such essential difference is that our denoising diffusion step requires multiple steps of reverse diffusion in a single iteration of outer loop, and has an accurate, simple expression for the (stationary) distribution it aims to sample from. In contrast, for projection-type approaches, it is typical that only one step of reverse diffusion is run in a single iteration of outer loop, and the distribution of its output, even in the ideal continuous setting, cannot be expressed in a simple way, which significantly hinders theoretical analysis. From a big picture, by running more steps of reverse diffusion every time after we do proximal consistency sampling, we stabilize the process since we better compensate for the drifts from the data manifold induced by (proximal) projection onto the null space of consistency equations. This is the heuristic reason for the improved robustness of our algorithm, which is backed by rigorous mathematical proofs.\n\n\n**Order of proximal consistency step and denoising diffusion step.**\n- Thank you for your sharp observation. The first reason for such an order is that the algorithm would have visually better outputs if ending with the denoising diffusion step instead of the proximal consistency step, as denoising diffusion by design outputs visually satisfactory images while proximal consistency sampler only cares about consistency. Such a difference was not present in previous projection-type approaches, as the denoising stage there runs only one step of reverse diffusion, thus it is insignificant whether it comes as the penultimate step or the last step. This again constitutes an example of the difference of our algorithm. The second reason is pedagogical, since this order better parallels that of projected gradient descent, which may facilitate the understanding of the logic of our algorithm.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
# General Response to All Reviewers

We would like to express our cordial thanks to all the reviewers for their careful review and constructive feedback. Below we address some common concerns raised by the reviewers. Our point-to-point response can be found in the separate rebuttal to each reviewer.

**Comparison with latent diffusion.** Our work is complementary (rather than parallel) to previous literatures on latent diffusion. It is mostly straightforward to incorporate latent diffusion models into our framework (it suffices to modify the Denoising Diffusion Sampler in our framework to a latent version), which will hopefully further allow our algorithm to take advantage of latent diffusion models. On the other hand, such pre-trained latent model is not always available for imaging tasks (e.g. the experiments of our work and many previous works [CKM+23, CLY23, MSKV24]), and our framework still works when only non-latent models are available. We opt to present only the non-latent version in the paper to better highlight the key novel ideas of our work (together with its theory) within limited pages and without compromising generality of settings, which is already quite lengthy. However, we agree that incorporating latent diffusion is a promising direction for future research.

**Experimental evaluation.** We compare our algorithm with the well-established DPS algorithm **as well as more recent, highly competitive algorithm** such as LGD-MC [SZY+23] which already demonstrates significant advantage over DPS. We abide to the experimental settings in these previous works with best efforts, and the size of testing dataset is kept strictly the same as that in DPS. Following the suggestions of reviewers, **we have also included a brief comparison with ReSample [SKZ+23] below and in the rebuttal pdf**, which we will try to make more complete in the final version. In addition, we would like to emphasize that none of the previous works have a provable robustness guarantee, thus our work is also of independent theoretical interest.

------------------------------
# Rebuttal pdf
Please find in attached our one-page pdf containing the additional experimental results cited in our responses. The settings of these results are explained below. 

**Comparison with ReSample.** For reasons explained above, we compare our algorithm with the non-latent (also termed "pixel-based" in the [SKZ+23]) version of ReSample, and leave the comparison of latent versions for future work after we developed a latent version of our algorithm. For ReSample, we use DDIM with $T=1000$ steps, which is the same as all the other algorithms. We let ReSample run without stochastic resampling for the first $250$ steps as specified in the original paper. After that, we run stochastic resampling every $10$ steps as in the original paper. All other parameters, including the resampling noise factor $\gamma$, the (conjugate) gradient stepsize, early stopping criterion, etc., has been fine-tuned with reasonable effort for best performance. We evaluate its performance on the same set of inverse problems on FFHQ-1k dataset as in our paper, i.e. super-resolution, phase retrieval, and quantized sensing (problem parameters was given in our paper). **The results are shown in Table 1 in the rebuttal pdf.** We make a few comments on the results here.
- Non-latent ReSample does not work on phase retrieval, despite our considerable effort in trying different parameters. This may be accounted by the complicated non-linear, non-convex landscape of the phase retrieval objective, as in such case the hard consistency enforcing step in ReSample may easily stuck at local minimum. It is possible that using the latent version of ReSample may alleviate this issue, but as we said, this awaits future research. 
- Overall, on the tasks where ReSample works, i.e., super-resolution and quantized sensing, the performance of ReSample is better than DPS and close to that of LGD-MC which we already compared with. All of DPS, LGD-MC and ReSample are surpassed by our algorithm in quantitative metrics.

**Demonstration of convergence.** We verify our theoretical claim of convergence of our algorithm on a Gaussian mixture model as a toy example. We consider the setting where the unconditional distribution $p^\star(x)$ is a 2D GMM given by
$$x = [x_1, x_2] \sim 0.6 \cdot\mathcal N([-3, -1]^\top, I_2) + 0.4 \cdot\mathcal N([1, 1]^\top, I_2).$$
We consider a linear inverse problem where $\mathcal{A}(x) = x_1-x_2$, and suppose the measurement $y=-0.5$. For simplicity of demonstration, we further consider the noiseless setting, so that the posterior distribution $p(x|y)$ is a degenerate distribution supported on the 1D line $\\{x\in\mathbb{R}^2: y=\mathcal{A}(x)\\}=\\{x\in\mathbb{R}^2: x_1-x_2=-0.5\\}$. **A depiction of the setting is given in Figure 2.(a) in the rebuttal pdf.** We run DPnP on this linear inverse problem assuming the unconditional score function is exactly known (which can be computed from the GMM expression). We use an annealing schedule with $K=120$, $\eta_0=0.5$, $\eta_K=10^{-4}$. **The distribution of $\hat x_k$, the $k$-th iterate of DPnP, is shown in Figure 1 in the rebuttal pdf.** We further compute the total variation distance between the distribution of $\hat x_k$ and the true posterior distribution. **The result of total variation error is shown in Figure 2.(b) in the rebuttal pdf.** It can be inferred clearly that the distribution of the iterates of DPnP eventually converges (in total variation) to the true posterior distribution.

Meta-review: 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.
Decision: Accept (poster)