Principled Probabilistic Imaging using Diffusion Models as Plug-and-Play Priors

Authors: zihui  wu
yu  sun
yifan  chen
bingliang  zhang
yisong  yue
katherine l bouman
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 05:14:01

Abstract: Diffusion models (DMs) have recently shown outstanding capabilities in modeling complex image distributions, making them expressive image priors for solving Bayesian inverse problems. However, most existing DM-based methods rely on approximations in the generative process to be generic to different inverse problems, leading to inaccurate sample distributions that deviate from the target posterior defined within the Bayesian framework. To harness the generative power of DMs while avoiding such approximations, we propose a Markov chain Monte Carlo algorithm that performs posterior sampling for general inverse problems by reducing it to sampling the posterior of a Gaussian denoising problem. Crucially, we leverage a general DM formulation as a unified interface that allows for rigorously solving the denoising problem with a range of state-of-the-art DMs. We demonstrate the effectiveness of the proposed method on six inverse problems (three linear and three nonlinear), including a real-world black hole imaging problem. Experimental results indicate that our proposed method offers more accurate reconstructions and posterior estimation compared to existing DM-based imaging inverse methods.

Main paper: 
 introduction: 
 Inverse problems arise in many computational imaging applications, where the goal is to recover an image x ∈ R n from a set of sparse and noisy measurements y ∈ R m . The relationship between x and y can be described byEQUATIONwhere A(•) : R n → R m is the forward operator (linear or nonlinear) and n is the random measurement noise in R m . Since the sparsity and noisiness of y often lead to significant uncertainty in x, it is preferable to sample the posterior distribution p(x|y) over all possible solutions based on some prior distribution p(x), rather than finding a single deterministic solution. Traditional posterior sampling methods often rely on simple image priors that do not reflect the sophistication of real-world image distributions. On the other hand, diffusion models (DMs) have recently emerged as a powerful tool for modeling highly complex image distributions [33, 62] . Nevertheless, it remains a challenge to turn DMs into reliable imaging inverse solvers, which motivates us to develop a principled Bayesian method that leverages DMs as priors for posterior sampling.Diffusion models generate samples from a distribution by reversing a diffusion process from the target distribution to a simple (usually Gaussian) distribution [33, 62] . In particular, it estimates a clean image x 0 from a noise image x T by successively denoising noisy images, where x t ∼ p t is the intermediate noisy image at time t ∈ [0, T ]. Reversing diffusion requires one to estimate the time-varying gradient log density (score function) ∇ log p t (x t ) along the diffusion process, or ∇ log p t (x t |y) in the case of sampling the posterior p(x|y). Figure 1 : Demonstration of the proposed method, PnP-DM, for posterior sampling using the real data for the M87 black hole from April 6 th , 2017 [21] . The black hole imaging problem is non-convex and highly ill-posed due to severe noise corruption and measurement sparsity. Our method rigorously integrates measurements from a real-world imaging system with an expressive image prior in the form of a diffusion model, which was trained with images from the GRMHD black hole simulation [22] in this case. Besides having high visual quality, our posterior samples accurately capture key features of the M87 black hole such as the bright spot location and ring diameter.To design generic DM-based inverse problem solvers, most existing methods attempt to approximate the time-varying gradient log density ∇ log p t (x t |y) [17, 72, 83, 61, 38, 57, 59, 43, 15, 76, 18, 54, 8] .In particular they first apply Bayes' rule to separate the forward operator from an unconditional prior over the intermediate noisy image x t :∇ log p t (x t |y) = ∇ log p t (y|x t ) + ∇ log p t (x t ).(EQUATIONBy instead aiming to evaluate the right hand side, one can leverage the existing pre-trained DMs for the unconditional term ∇ log p t (x t ). However, the main challenge in this case is that ∇ log p t (y|x t ) is intractable to compute in general, as p t (y|x t ) involves an integral over all possible x 0 's that could give rise to x t [17] . Various methods have been proposed to circumvent the intractability and can mostly be categorized into two groups. One group of methods explicitly approximate ∇ log p t (y|x t ) by making simplifying assumptions [61, 17, 59, 8] . However, even for arguably the finest approximation to date proposed in the recent work [8] , it is exact only when the prior distribution p(x) is Gaussian. For general prior distributions beyond Gaussian, these methods do not sample the true posterior p(x|y). The other group of methods do not make explicit approximations but instead substitute ∇ log p t (y|x t ) with empirically designed updates where y is treated as a guidance signal [72, 83, 38, 57, 43, 15, 76, 18, 54] . Although these methods may have strong empirical performance, they have deviated from the Bayesian formulation and no longer aim to sample the target posterior. In summary, these existing DM-based inverse methods should be best viewed as guidance methods, where the generative process is guided towards the regions where the measurement y is more likely to be observed, not as posterior sampling methods [8] . We also note that some recent work considered combing DMs with Sequential Monte Carlo to ensure asymptotic consistency in posterior sampling [11, 23] , but the investigation has been limited to linear imaging inverse problems.Our contributions In this work, we pursue a different path towards posterior sampling with DM priors by proposing a new Markov chain Monte Carlo (MCMC) algorithm, which we call Plug-and-Play Diffusion Models (PnP-DM). It incorporates DMs in a principled way and circumvents the approximation required when taking the approach in (2) . The proposed algorithm is based on the Split Gibbs Sampler [70] that alternates between two sampling steps that separately involve the likelihood and prior. While the likelihood step can be tackled with traditional sampling techniques, the prior step involves a Bayesian denoising problem that requires careful design. Importantly, we identify a connection between the Bayesian denoising problem and the unconditional image generation problem under a general formulation of DMs presented in [36] (which is referred to as the EDM formulation hereafter). This connection allows us to perform rigorous posterior sampling for denoising using DMs without approximating the generative process and enables the use of a wide range of pretrained DMs through the unified EDM formulation. We present an analysis on the non-asymptotic behavior of PnP-DM by establishing a stationarity guarantee in terms of the average Fisher information. We further demonstrate the strong empirical performance of PnP-DM by investigating three linear and three nonlinear noisy inverse problems, including a black hole interferometric imaging problem involving real data that is both nonlinear and severely ill-posed (see Figure 1 ). Overall, PnP-DM outperforms existing baseline methods, achieving higher accuracy in posterior estimation.Existing works related to SGS Several works have designed algorithms for solving imaging inverse problems based on SGS [52, 19, 6, 27, 77] . The key distinction among these methods lies in their approaches to the prior step. For instance, the works [52, 6, 27] applied Langevin-based updates for sampling π X|Z=z such that the prior information is encoded by either traditional regularizers or off-the-shelf image denoisers. The work [19] tackled the prior step by heuristically customizing a diffusion model (i.e. DDPM [33] ) for sampling π X|Z=z . A concurrent work [77] improved the implementation by devising two diffusion processes that rigorously solve the prior step. Our method differs from [77] by connecting the prior step to the EDM formulation [36] . This connection allows us to seamlessly integrate state-of-the-art DMs as expressive image priors for Bayesian inference through a unified interface, eliminating the need for additional customization for each model and leading to better empirical performance. We also note the recent work [42] that adopted the optimizationbased variable-splitting formulation of HQS and utilized general DMs as image priors. We instead considers the SGS formulation from a Bayesian posterior sampling standpoint. Additionally, while SGS-based methods theoretically accommodate general inverse problems, empirical evidence on real-world nonlinear inverse problems remains scarce in the literature. In this work, we demonstrate our method on three nonlinear inverse problems, including a black hole imaging problem. For a more comprehensive review of related works, see Appendix E. 
method: 
 A schematic diagram for the proposed method is shown in Figure 2 . Our method, dubbed PnP-DM, builds upon the SGS framework with rigorous implementations of the two sampling steps and an annealing schedule for the coupling parameter ρ. We start with our implementations of the first step for solving both linear and nonlinear inverse problems. Prior step: denoising with noise level < l a t e x i t s h a 1 _ b a s e 6 4 = " c F Y B E x c r t j + I i J N g M p X I J t 9 a e M A = " > A A A C W n i c b Z H L S g M x F I b T 8 T 7 e W n X n Z r A I r s q M S H U l R T c u F W w r 2 F L O p G k b m 8 u Y n F G H M u / g V t 9 M 8 G F M a x e 2 9 k D g z / e f A z l / 4 k R w i 2 H 4 V f C W l l d W 1 9 Y 3 / M 2 t 7 Z 3 d Y m m v Y X V q K K t T L b R 5 i M E y w R W r I 0 f B H h L D Q M a C N e P h 9 d h v v j B j u V b 3 m C W s L a G v e I 9 T Q I c a L T P Q n W G n W A 4 r 4 a S C / y K a i j K Z 1 m 2 n V L h s d T V N J V N I B V j 7 G I U J t k d g k F P B c r + V W p Y A H U K f P T q p Q D L b H k 2 e m w f H j n S D n j b u K A w m 9 O / E C K S 1 m Y x d p w Q c 2 H l v D B d 6 7 F m B M Z A t N G O Z z 9 7 n g U w F c q N f Z y l C n A o w b 7 l / P L M U 9 i 7 a I 6 6 S F J G r / u 9 S v V Q E q I N x z k G X G 0 Z R Z E 4 A N d z l E t A B G K D o f s N 3 i U f z + f 4 X j d N K V K 1 U 7 8 7 K t a t p 9 u v k k B y R E x K R c 1 I j N + S W 1 A k l T + S d f J D P w r f n e R v e 5 m + r V 5 j O 7 J O Z 8 g 5 + A D o 2 t 2 I = < / l a t e x i t > ⇢k start from at noise level < l a t e x i t s h a 1 _ b a s e 6 4 = " s B j + T S l L s N g E n D 3 7 E q w o K 1 l V s n g = " > A A A C W n i c b Z H L T g I x F I b L e A W 8 g L p z M 5 G Q u C I z x q A r Y 3 T j E h O 5 J E D I m V K g 0 s v Y n l E J 4 R 3 c 6 p u Z + D C W y 0 K Q k z T 5 + / 3 n J D 1 / o 1 h w i 0 H w n f I 2 N r e 2 d 3 b T m e z e / s F h L n 9 U s z o x l F W p F t o 0 I r B M c M W q y F G w R m w Y y E i w e j S 8 n / r 1 V 2 Y s 1 + o J R z F r S + g r 3 u M U 0 K F a y w x 0 J + j k C k E p m J X / X 4 Q L U S C L q n T y q Z t W V 9 N E M o V U g L X N M I i x P Q a D n A o 2 y b Q S y 2 K g Q + i z p p M K J L P t 8 e y 5 E 7 / o S N f v a e O O Q n 9 G / 0 6 M Q V o 7 k p H r l I A D u + p N 4 V q P v S g w B k Z r z U h O l u + r Q C Y C u d F v y x Q h S g S Y 9 0 m m u L Q U 9 q 7 b Y 6 7 i B J G r / n y p X i J 8 1 P 4 0 Z 7 / L D a M o R k 4 A N d z l 4 t M B G K D o f i P j E g 9 X 8 / 0 v a h e l s F w q P 1 4 W b u 8 W 2 e + S U 3 J G z k l I r s g t e S A V U i W U P J M P 8 k m + U j + e 5 6 W 9 7 L z V S y 1 m j s l S e S e / x 5 y 3 J w = = < / l a t e x i t > Our method alternates between a likelihood step that enforces data consistency and a prior step that solves a denoising posterior sampling problem by leveraging the Split Gibbs Sampler [70] . An annealing schedule controls the strength of the two steps at each iteration to facilitate efficient and accurate sampling. A crucial part of our design is the prior step, where we identify a key connection to a general diffusion model framework called the EDM [36] . This connection allows us to easily incorporate a family of state-of-the-art diffusion models as priors to conduct posterior sampling in a principled way without additional training. Our method demonstrates strong performance on a variety of linear and nonlinear inverse problems.For the likelihood step at iteration k, we sampleEQUATION)Linear forward model and Gaussian noise We first consider a simple yet common case where the forward model A is linear and the noise distribution is zero-mean Gaussian, i.e. A := A ∈ R m×n and n ∼ N (0, Σ). In this case, the potential function of the likelihood term is f (x; y) = 1 2 ∥y -Ax∥ 2(up to an additive constant that does not depend on x and y) where∥ • ∥ 2 Σ := ⟨•, Σ -1 •⟩. It is then straightforward to show that π Z|X=x = N (m(x), Λ -1 )whereΛ := A T Σ -1 A + 1 ρ 2 I and m(x) := Λ -1 (A T Σ -1 y + 1 ρ 2 x).The problem of sampling from Gaussian distributions has been systematically studied [71] . We refer readers to Appendix C.1 for a more detailed discussion.General case For general nonlinear inverse problems, the likelihood step is not sampling from a Gaussian distribution anymore. Nevertheless, since we have access to π Z|X=x in closed form up to a multiplicative factor, we can use Monte Carlo methods based on Langevin dynamics to draw samples from it as long as the likelihood potential is differentiable. Specifically, we first set up the following Langevin SDE that admits π Z|X=x as the stationary distributiondz t = ∇ log π Z|X=x (z t )dt + √ 2dw t = -∇f (z; y) - 1 ρ 2 (z -x) dt + √ 2dw t .We then initialize the SDE at z 0 = x and run it with Euler discretization. The pseudocode is provided in Appendix C.1.For the prior step at iteration k, we sampleEQUATIONA closer examination of (6) reveals that this prior step is essentially to draw posterior samples for a Gaussian denoising problem, where the "measurement" is z (k) , the noise level is ρ, and the prior distribution is p(x) ∝ exp(-g(x)).We tackle this denoising posterior sampling problem within SGS using DMs as image priors. In particular, we leverage the EDM framework [36] , which was originally proposed to unify various formulations of DMs for unconditional image generation. To see the connection of the EDM framework to (6) , consider a family of mollified distributions p(x; σ) given by adding i.i.d Gaussian noise of standard deviation σ to the prior distribution p(x), i.e. x + σϵ ∼ p(x; σ). The core idea of the EDM framework is that a variety of state-of-the-art DMs can be unified into the following reverse SDE:EQUATIONwhere wt is an n-dimensional Wiener process running backward in time, σ(t) > 0 is a pre-defined noise level schedule with σ(0) = 0, s(t) is a pre-defined scaling schedule, and σ(t), ṡ(t) are their time derivatives. As shown in [36] , the defining property of ( 7) is that x t /s(t) ∼ p(x; σ(t)) for any time t. Therefore, solving this SDE backward in time allows us to travel from any noise level σ(t) to the clean image distribution at t = 0. This means that we can use (7) to solve (6) with arbitrary noise level ρ as long as ρ is within the range of σ(t). Indeed, the distribution of x 0 conditioned on x t isp(x 0 |x t ) ∝ p(x t |x 0 )p(x 0 ) ∝ N (s(t)x 0 , s(t) 2 σ(t) 2 I) exp(-g(x 0 )) ∝ exp -g(x 0 ) - 1 2σ(t) 2 ∥x 0 -x t /s(t)∥ 2 2 .We highlight that the last line exactly matches (6) when x t = s(t)z (k) and σ(t) = ρ. Therefore, we can naturally design a practical algorithm that samples (6) by following these three steps: (1) find t * such that σ(t * ) = ρ, (2) initialize at x t * = s(t * )z (k) , and (3) solve ( 7) backward from t * to 0 by choosing the discretization time steps and integration scheme. Through this unified interface, any DMs, once converted to the EDM formulation, can be directly turned into a rigorous solver for (6) .Leveraging the connection with EDM, our prior step implementation comes with a large design space that encompasses a variety of existing DMs, such as DDPM (or VP-SDE) [33] , VE-SDE [62] , and iDDPM [51] . In our experiments, we conduct posterior sampling with all these different models within our framework and all of them provide high-quality samples. The pseudocode of our implementation and more details on the EDM formulation for the prior step is given in Appendix C.2.The pseudocode of PnP-DM in complete form is presented in Algorithm 1. PnP-DM alternates between the two sampling steps with an annealing schedule {ρ k } for the coupling parameter. We find that the annealing schedule on ρ accelerates the mixing time of the Markov chain and prevents the algorithm from getting stuck in bad local minima for solving highly ill-posed inverse problems. This is a common practice in both Langevin-based [39, 34, 64] and SGS-based [6, 77] MCMC algorithms to improve the empirical performance in solving inverse problems.Our work shares some similarities with PnP-SGS [19] but contains three main key differences. First, as demonstrated in our experiments, we investigate three nonlinear inverse problems, while nonlinear inverse problems are beyond the scope of [19] . Our experiments show that PnP-SGS struggles with challenging nonlinear inverse problems such as Fourier phase retrieval. Second, we adopt the EDM formulation to ensure that the prior step of PnP-DM is a rigorous mapping from the image manifold with the desired noise level to the clean image manifold, aligning with the theory of SGS. In contrast, the prior step of PnP-SGS [19] is heuristic (which is also pointed out by [77] ) and not rigorously designed to sample (6) . Third, unlike PnP-SGS [19] that uses a constant ρ, we consider an annealing schedule {ρ k } for the coupling parameter, which is important for highly ill-posed inverse problems.Algorithm 1 Plug-and-Play Diffusion Models (PnP-DM)Input: initialization x 0 ∈ R n , total number of iterations K > 0, coupling strength schedule {ρ k > 0} K-1 k=0 , likelihood potential f ( • ; y) with measurements y ∈ R m , pretrained model D θ ( • ; • ) that approximates ∇ log p (x; σ) with (D θ (x; σ) -x)/σ 2 . 1: for k = 0, ..., K -1 do 2: z (k) ← LikelihoodStep(x (k) , ρ k , f ( • ; y)) ▷ Section 3.1 3: x (k+1) ← PriorStep(z (k) , ρ k , D θ ( • ; • )) ▷ Section 3.2 4: end for 5: return x (k+1) 3.4 Theoretical insights < l a t e x i t s h a 1 _ b a s e 6 4 = " o 7 a Q L 1 G D 1 P K N b B c + x 4 w 4 N U H 5 v D Q = " > A A A C W 3 i c b Z H L T g I x F I b L 4 A V H V N C 4 c j O R k L g i M 8 a g K 0 N 0 4 x I T u S S A 5 E w p 2 N B 2 x v a M S g g P 4 V a f z I X v Y r k s 5 H K S J n + / / 5 y k 5 2 8 Y C 2 7 Q 9 3 9 S T n p r e 2 c 3 s + f u Z w 8 O j 3 L 5 4 7 q J E k 1 Z j U Y i 0 s 0 Q D B N c s R p y F K w Z a w Y y F K w R D u + n f u O N a c M j 9 Y S j m H U k D B T v c w p o U a O t k q 7 / 3 O z m C n 7 J n 5 W 3 L o K F K J B F V b v 5 1 G 2 7 F 9 F E M o V U g D G t w I + x M w a N n A o 2 c d u J Y T H Q I Q x Y y 0 o F k p n O e P b e i V e 0 p O f 1 I 2 2 P Q m 9 G / 0 + M Q R o z k q H t l I A v Z t W b w o 0 e e 1 W g N Y w 2 m q G c L N 9 X g U w E c h 2 9 L 1 O E M B G g P y Z u c W k p 7 N 9 0 x l z F C S J X g / l S / U R 4 G H n T o L 0 e 1 4 y i G F k B V H O b i 0 d f Q A N F + x 2 u T T x Y z X d d 1 C 9 L Q b l U f r w q V O 4 W 2 W f I G T k n F y Q g 1 6 R C H k i V 1 A g l Q / J J v s h 3 6 t d J O 6 6 T n b c 6 q c X M C V k q 5 / Q P i k q 3 g Q = = < / l a t e x i t > ⌫ X 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " c w x S e n a / B S 0 D 6 0 4 l 8 M 2 h Q 6 w M y + o = " > A A A C W 3 i c b Z H L S g M x F I Y z U y 9 1 v C u u 3 A y W g q s y I 6 K u R H T j U s F e s K 3 l T J r W 0 C Q z J m f U o c x D u N U n c + G 7 m N Y u 7 O V A 4 M / 3 n w M 5 f 6 J E c I N B 8 O 2 4 h a X l l d X i m r e + s b m 1 v b O 7 V z N x q i m r 0 l j E u h G B Y Y I r V k W O g j U S z U B G g t W j w c 3 I r 7 8 y b X i s H j B L W F t C X / E e p 4 A W 1 V s q 7 Q R P j 5 2 d U l A J x u X P i 3 A i S m R S d 5 1 d 5 7 L V j W k q m U I q w J h m G C T Y H o J G T g X L v V Z q W A J 0 A H 3 W t F K B Z K Y 9 H L 8 3 9 8 u W d P 1 e r O 1 R 6 I / p / 4 k h S G M y G d l O C f h s Z r 0 R X O i x F w V a Q 7 b Q j G Q + f Z 8 F M h X I d f w 2 T R G i V I B + z 7 3 y 1 F L Y u 2 g P u U p S R K 7 6 f 0 v 1 U u F j 7 I + C 9 r t c M 4 o i s w K o 5 j Y X n z 6 D B o r 2 O z y b e D i b 7 7 y o n V T C s 8 r Z / W n p 6 n q S f Z E c k i N y T E J y T q 7 I L b k j V U L J g H y Q T / L l / L g F 1 3 M 3 / l p d Z z K z T 6 b K P f g F j i y 3 g w = = < / l a t e x i t > ⌫ Z 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 H / u M 3 8 N M 1 u 1 j g i G D q P m + V R b y m Y = " > A A A C W 3 i c b Z H L T g I x F I b L 4 A V H V N C 4 c j O R k L g i M 8 a g K 0 N 0 4 x I T u S S A 5 E w p 2 N B 2 x v a M S g g P 4 V a f z I X v Y r k s 5 H K S J n + / / 5 y k 5 2 8 Y C 2 7 Q 9 3 9 S T n p r e 2 c 3 s + f u Z w 8 O j 3 L 5 4 7 q J E k 1 Z j U Y i 0 s 0 Q D B N c s R p y F K w Z a w Y y F K w R D u + n f u O N a c M j 9 Y S j m H U k D B T v c w p o U a O t k m 7 w 3 O z m C n 7 J n 5 W 3 L o K F K J B F V b v 5 1 G 2 7 F 9 F E M o V U g D G t w I + x M w a N n A o 2 c d u J Y T H Q I Q x Y y 0 o F k p n O e P b e i V e 0 p O f 1 I 2 2 P Q m 9 G / 0 + M Q R o z k q H t l I A v Z t W b w o 0 e e 1 W g N Y w 2 m q G c L N 9 X g U w E c h 2 9 L 1 O E M B G g P y Z u c W k p 7 N 9 0 x l z F C S J X g / l S / U R 4 G H n T o L 0 e 1 4 y i G F k B V H O b i 0 d f Q A N F + x 2 u T T x Y z X d d 1 C 9 L Q b l U f r w q V O 4 W 2 W f I G T k n F y Q g 1 6 R C H k i V 1 A g l Q / J J v s h 3 6 t d J O 6 6 T n b c 6 q c X M C V k q 5 / Q P j D 2 3 g g = = < / l a t e x i t > ⌫ X 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " F Y 3 E c d L c s o I Z z Y o E c o e E T 0 n o W U c = " > A A A C X 3 i c b Z H L S g M x F I b T 8 V b r r e p K 3 A w W w Y 1 l R k R d S d G N 4 E b B t m J b y 5 k 0 b Y O 5 j M k Z t Q z z G G 7 1 u V z 6 J q a 1 C 9 t 6 I P D n + 8 + B n D 9 R L L j F I P j K e X P z C 4 t L + e X C y u r a + k Z x c 6 t m d W I o q 1 I t t L m P w D L B F a s i R 8 H u Y 8 N A R o L V o 6 f L o V 9 / Y c Z y r e 5 w E L O W h J 7 i X U 4 B H W o 0 V d J O r w / D 7 P G h X S w F 5 W B U / q w I x 6 J E x n X T 3 s y d N z u a J p I p p A K s b Y R B j K 0 U D H I q W F Z o J p b F Q J + g x x p O K p D M t t L R m z N / 3 5 G O 3 9 X G H Y X + i P 6 d S E F a O 5 C R 6 5 S A f T v t D e G / H n t W Y A w M / j U j m U 3 e p 4 F M B H K j X y c p Q p Q I M G 9 Z Y X 9 i K e y e t V K u 4 g S R q 9 7 v U t 1 E + K j 9 Y d h + h x t G U Q y c A G q 4 y 8 W n f T B A 0 X 1 J w S U e T u c 7 K 2 p H 5 f C k f H J 7 X K p c j L P P k 1 2 y R w 5 I S E 5 J h V y R G 1 I l l G j y T j 7 I Z + 7 b W / L W v e J v q 5 c b z 2 y T i f J 2 f g D 8 7 L k c < / l a t e x i t > ⌫ Z K 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 b a 6 X 5 0 2 w R B v C r G y Q A D 3 o + T 2 N Z I = " > A A A C W 3 i c b Z H N T g I x F I X L + I c j K m h c u Z l I S F y R G W P Q l T G 6 M X G D i Q g J I L l T C j S 0 n b G 9 o 5 I J D + F W n 8 y F 7 2 J B F v J z k y a n 3 7 k 3 6 T 0 N Y 8 E N + v 5 3 x l l b 3 9 j c y m 6 7 O 7 n d v f 1 8 4 e D J R I m m r E Y j E e l G C I Y J r l g N O Q r W i D U D G Q p W D 4 e 3 E 7 / + y r T h k X r E U c z a E v q K 9 z g F t K j e U k n n / r n R y R f 9 s j 8 t b 1 k E M 1 E k s 6 p 2 C p m r V j e i i W Q K q Q B j m o E f Y z s F j Z w K N n Z b i W E x 0 C H 0 W d N K B Z K Z d j p 9 7 9 g r W d L 1 e p G 2 R 6 E 3 p f 8 n U p D G j G R o O y X g w C x 6 E 7 j S Y y 8 K t I b R S j O U 4 / n 7 I p C J Q K 6 j t 3 m K E C Y C 9 P v Y L c 0 t h b 3 L d s p V n C B y 1 f 9 b q p c I D y N v E r T X 5 Z p R F C M r g G p u c / H o A D R Q t N / h 2 s S D x X y X x d N Z O a i U K w / n x e u b W f Z Z c k x O y C k J y A W 5J n e k S m q E k i H 5 I J / k K / P j r D m u k / t r d T K z m U M y V 8 7 R L 7 7 r t 5 w = < / l a t e x i t > ⌫ X K Non-stationary processPrior step < l a t e x i t s h a 1 _ b a s e 6 4 = " m 0D J h V h p g M O U V m 4 9 X m p k U y c U 1 + I = " > A A A C W X i c b Z H L T g I x F I Y 7 4 w 3 H G 8 j S z U R C 4 o r M G I O u D N G N S 0 z k k g C S M 6 V A Q 9 s Z 2 z M q I T y D W 3 0 0 4 8 t Y L g t B T t L k 7 / e f k / T 8 j R L B D Q b B t + N u b e / s 7 m X 2 v Y P D o + O T b O 6 0 b u J U U 1 a j s Y h 1 M w L D B F e s h h w F a y a a g Y w E a 0 S j + 5 n f e G X a 8 F g 9 4 T h h H Q k D x f u c A l p U a y f 8 u d n N F o J S M C / / v w i X o k C W V e 3 m n N t 2 L 6 a p Z A q p A G N a Y Z B g Z w I a O R V s 6 r V T w x K g I x i w l p U K J D O d y f y 1 U 7 9 o S c / v x 9 o e h f 6 c / p 2 Y g D R m L C P b K Q G H Z t 2 b w Y 0 e e 1 G g N Y w 3 m p G c r t 7 X g U w F c h 2 / r V K E K B W g 3 6 d e c W U p 7 N 9 0 J l w l K S J X g 8 V S / V T 4 G P u z m P 0 e 1 4 y i G F s B V H O b i 0 + H o I G i / Q z P J h 6 u 5 / t f 1 C 9 L Y b l U f r w q V O 6 W 2 W f I G T k n F y Q k 1 6 R C H k i V 1 A g l n H y Q T / L l / L i O m 3 G 9 R a v r L G f y Z K X c / C 8 S c b b U < / l a t e x i t > ⇡ X< l a t e x i t s h a 1 _ b a s e 6 4 = " z 8 I A A Z 3 y H r M s c Q 8 + F r u p o d s 7 G 9 o = " > A A A C W X i c b Z H N S g M x F I U z 4 1 8 d f 6 t L N 4 O l 4 K r M i F R X I r p x W c G 2 Y j u W O 2 n a B p P M m N x R S 5 l n c K u P J r 6 M a Z 2 F r V 4 I n H z n X s g 9 i V P B D Q b B p + M u L a + s r p X W v Y 3 N r e 2 d 3 f J e y y S Z p q x J E 5 H o u x g M E 1 y x J n I U 7 C 7 V D G Q s W D t + v J r 6 7 W e m D U / U L Y 5 T F k k Y K j 7 g F N C i Z j f l D / e 9 3 U p Q C 2 b l / x V h I S q k q E a v 7 J x 3 + w n N J F N I B R j T C Y M U o w l o 5 F S w 3 O t m h q V A H 2 H I O l Y q k M x E k 9 l r c 7 9 q S d 8 f J N o e h f 6 M / p 6 Y g D R m L G P b K Q F H Z ti G F s B V H O b i 0 9 H o I G i / Q z P J h 4 u 5 v t X t I 5 r Y b 1 W v z m p X F w W 2 Z f I A T k k R y Q k p + S C X J M G a R J K O H k j 7 + T D + X I d t +D J h V h p g M O U V m 4 9 X m p k U y c U 1 + I = " > A A A C W X i c b Z H L T g I x F I Y 7 4 w 3 H G 8 j S z U R C 4 o r M G I O u D N G N S 0 z k k g C S M 6 V A Q 9 s Z 2 z M q I T y D W 3 0 0 4 8 t Y L g t B T t L k 7 / e f k / T 8 j R L B D Q b B t + N u b e / s 7 m X 2 v Y P D o + O T b O 6 0 b u J U U 1 a j s Y h 1 M w L D B F e s h h w F a y a a g Y w E a 0 S j + 5 n f e G X a 8 F g 9 4 T h h H Q k D x f u c A l p U a y f 8 u d n N F o J S M C / / v w i X o k C W V e 3 m n N t 2 L 6 a p Z A q p A G N a Y Z B g Z w I a O R V s 6 r V T w x K g I x i w l p U K J D O d y f y 1 U 7 9 o S c / v x 9 o e h f 6 c / p 2 Y g D R m L C P b K Q G H Z t 2 b w Y 0 e e 1 G g N Y w 3 m p G c r t 7 X g U w F c h 2 / r V K E K B W g 3 6 d e c W U p 7 N 9 0 J l w l K S J X g 8 V S / V T 4 G P u z m P 0 e 1 4 y i G F s B V H O b i 0 + H o I G i / Q z P J h 6 u 5 / t f 1 C 9 L Y b l U f r w q V O 6 W 2 W f I G T k n F y Q k 1 6 R C H k i V 1 A g l n H y Q T / L l / L i O m 3 G 9 R a v r L G f y Z K X c / C 8 S c b b U < / l a t e x i t > ⇡ X< l a t e x i t s h a 1 _ b a s e 6 4 = " z 8 I A A Z 3 y H r M s c Q 8 + F r u p o d s 7 G 9 o = " > A A A C W X i c b Z H N S g M x F I U z 4 1 8 d f 6 t L N 4 O l 4 K r M i F R X I r p x W c G 2 Y j u W O 2 n a B p P M m N x R S 5 l n c K u P J r 6 M a Z 2 F r V 4 I n H z n X s g 9 i V P B D Q b B p + M u L a + s r p X W v Y 3 N r e 2 d 3 f J e y y S Z p q x J E 5 H o u x g M E 1 y x J n I U 7 C 7 V D G Q s W D t + v J r 6 7 W e m D U / U L Y 5 T F k k Y K j 7 g F N C i Z j f l D / e 9 3 U p Q C 2 b l / x V h I S q k q E a v 7 J x 3 + w n N J F N I B R j T C Y M U o w l o 5 F S w 3 O t m h q V A H 2 H I O l Y q k M x E k 9 l r c 7 9 q S d 8 f J N o e h f 6 M / p 6 Y g D R m L G P b K Q F H Z ti G F s B V H O b i 0 9 H o I G i / Q z P J h 4 u 5 v t X t I 5 r Y b 1 W v z m p X F w W 2 Z f I A T k k R y Q k p + S C X J M G a R J K O H k j 7 + T D + X I d t +D J h V h p g M O U V m 4 9 X m p k U y c U 1 + I = " > A A A C W X i c b Z H L T g I x F I Y 7 4 w 3 H G 8 j S z U R C 4 o r M G I O u D N G N S 0 z k k g C S M 6 V A Q 9 s Z 2 z M q I T y D W 3 0 0 4 8 t Y L g t B T t L k 7 / e f k / T 8 j R L B D Q b B t + N u b e / s 7 m X 2 v Y P D o + O T b O 6 0 b u J U U 1 a j s Y h 1 M w L D B F e s h h w F a y a a g Y w E a 0 S j + 5 n f e G X a 8 F g 9 4 T h h H Q k D x f u c A l p U a y f 8 u d n N F o J S M C / / v w i X o k C W V e 3 m n N t 2 L 6 a p Z A q p A G N a Y Z B g Z w I a O R V s 6 r V T w x K g I x i w l p U K J D O d y f y 1 U 7 9 o S c / v x 9 o e h f 6 c / p 2 Y g D R m L C P b K Q G H Z t 2 b w Y 0 e e 1 G g N Y w 3 m p G c r t 7 X g U w F c h 2 / r V K E K B W g 3 6 d e c W U p 7 N 9 0 J l w l K S J X g 8 V S / V T 4 G P u z m P 0 e 1 4 y i G F s B V H O b i 0 + H o I G i / Q z P J h 6 u 5 / t f 1 C 9 L Y b l U f r w q V O 6 W 2 W f I G T k n F y Q k 1 6 R C H k i V 1 A g l n H y Q T / L l / L i O m 3 G 9 R a v r L G f y Z K X c / C 8 S c b b U < / l a tE i 2 C c 0 0 = " > A A A C X 3 i c b Z H L S g M x F I b T 8 V b r r e p K 3 A y W g h v L j I i 6 k q I b w Y 2 C v U B b y 5 k 0 r c F c x u S M W o Z 5 D L f 6 X C 5 9 E 9 P a h W 0 9 E P j z / e d A z p 8 o F t x i E H z l v I X F p e W V / G p h b X 1 j c 6 u 4 v V O 3 O j G U 1 a g W 2 j Q j s E x w x W r I U b B m b B j I S L B G 9 H Q 1 8 h s v z F i u 1 T 0 O Y 9 a R M F C 8 z y m g Q 6 2 2 S r r p z V G Y P T S 7 x V J Q C c b l z 4 t w I k p k U r f d 7 d x F u 6 d p I p l C K s D a V h j E 2 E n B I K e C Z Y V 2 Y l k M 9 A k G r O W k A s l s J x 2 / O f P L j v T 8 v j b u K P T H 9 O 9 E C t L a o Y x c p w R 8 t L P e C P 7 r s W c F x s D w X z O S 2 f R 9 F s h E I D f 6 d Z o i R I k A 8 5 Y V y l N L Y f + 8 k 3 I V J 4 h c D X 6 X 6 i f C R + 2 P w v Z 7 3 D C K Y u g E U M N d L j 5 9 B A M U 3 Z c U X O LC 6 Z K m / v B / k K u R o = < / l a t e x i t > ⌫ X K 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " m 0 D J h V h p g M O U V m 4 9 X m p k U y c U 1 + I = " > A A A C W X i c b Z H L T g I x F I Y 7 4 w 3 H G 8 j S z U R C 4 o r M G I O u D N G N S 0 z k k g C S M 6 V A Q 9 s Z 2 z M q I T y D W 3 0 0 4 8 t Y L g t B T t L k 7 / e f k / T 8 j R L B D Q b B t + N u b e / s 7 m X 2 v Y P D o + O T b O 6 0 b u J U U 1 a j s Y h 1 M w L D B F e s h h w F a y a a g Y w E a 0 S j + 5 n f e G X a 8 F g 9 4 T h h H Q k D x f u c A l p U a y f 8 u d n N F o J S M C / / v w i X o k C W V e 3 m n N t 2 L 6 a p Z A q p A G N a Y Z B g Z w I a O R V s 6 r V T w x K g I x i w l p U K J D O d y f y 1 U 7 9 o S c / v x 9 o e h f 6 c / p 2 Y g D R m L C P b K Q G H Z t 2 b w Y 0 e e 1 G g N Y w 3 m p G c r t 7 X g U w F c h 2 / r V K E K B W g 3 6 d e c W U p 7 N 9 0 J l w l K S J X g 8 V S / V T 4 G P u z m P 0 e 1 4 y i G F s B V H O b i 0 + H o I G i / Q z P J h 6 u 5 / t f 1 C 9 L Y b l U f r w q V O 6 W 2 W f I G T k n F y Q k 1 6 R C H k i V 1 A g l n H y Q T / L l / L i O m 3 G 9 R a v r L G f y Z K X c / C 8 S c b b U < / l a t e x i t > ⇡ X< l a t e x i t s h a 1 _ b a s e 6 4 = " Q j 5 M 6 K g 7 + F T R c 8 b R u 2 + 9 j 1 + P 8 p c = " We provide some theoretical insights on the non-asymptotic behavior of PnP-DM. We start with the following definitions.> A A A C W n i c b Z H L S g M x F I b T 8 d a L l 1 b d u R k s B V d l R q S 6 U Y p u X F a w F 2 h L O Z O m b W y S G Z M z a i l 9 B 7 f 6 Z o I P Y 3 p Z 2 M u B w J / v P w d y / g S R 4 A Y 9 7 y f h b G 3 v 7 O 4 l U + n M / s H h U T Z 3 X D N h r C m r 0 l C E u h G A Y Y I r V k W O g j U i z U A G g t W D 4 c P U r 7 8 x b X i o n n E U s b a E v u I 9 T g E t q r U Q 4 l u v k 8 1 7 R W 9 W 7 r r w F y J P F l X p 5 B J 3 r W 5 I Y 8 k U U g H G N H 0 v w v Y Y N H I q 2 C T d i g 2 L g A 6 h z 5 p W K p D M t M e z 5 0 7 c g i V d t x d q e x S 6 M / p / Y g z S m J E M b K c E H J h V b w o 3 e u x V g d Y w 2 m g G c r J 8 X w U y F s h 1 + L 5 M E Y J Y g P 6 Y p A t L S 2 H v p j 3 m K o o R u e r P l + r F w s X Q n e b s d r l m F M X I C q C a 2 1 x c O g A N F O 1 v p G 3 i / m q + 6 6 J 2 W f R L x d L T V b 5 8 v 8 g + S c 7 I O b k g P r k m Z f J I K q R K K H k h n + S L f C d + H c d J O Z l 5 q 5 N Y z J y Q p X J O / w C H a L c G < / l a t e x i t > ⌧ = 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " D G + b a F K Q Q 3 B N H 5 k b a m k j + E u c t Q o = " > A A A C W n i c b Z H L T g I x F I b L e O P i B d S d m 4 m E x B W Z M Q b d a I h u X G L C L Q F C z p Q C l b Y z t m d U Q n g H t / p m J j 6 M 5 b I Q 5 C R N / n 7 / O U n P 3 y A S 3 K D n f S e c r e 2 d 3 b 1 k K p 3 Z P z g 8 y u a O 6 y a M N W U 1 G o p Q N w M w T H D F a s h R s G a k G c h A s E Y w e p j 5 j V e m D Q 9 V F c c R 6 0 g Y K N 7 n F N C i e h s h v q 1 2 s 3 m v 6 M 3 L / S / 8 p c i T Z V W 6 u c R d u x f S W D K F V I A x L d + L s D M B j Z w K N k 2 3 Y 8 M i o C M Y s J a V C i Q z n c n 8 u V O 3 Y E n P 7 Y f a H o X u n P 6 d m I A 0 Z i w D 2 y k B h 2 b d m 8 G N H n t R o D W M N 5 q B n K 7 e 1 4 G M B X I d v q 1 S h C A W o N + n 6 c L K U t i / 6 U y 4 i m J E r g a L p f q x c D F 0 Z z m 7 P a 4 Z R T G 2 A q j m N h e X D k E D R f s b a Z u 4 v 5 7 v f 1 G / L P q l Y u n p K l + + X 2 a f J G f k n F w Q n 1 y T M n k k F V I j l D y T D / J J v h I / j u O k n M y i 1 U k s Z 0 7 I S j m n v 8 1 M t y o = < / l a t e x i t > ⌧ = T Time from initialization ( ) < l a t e x i t s h a 1 _ b a s e 6 4 = " R m Y I R w j A C x 5 v 1 B s I Q 0 g 8 s o V M B l Y = " > A A A E M X i c r V P L b t N A F J 3 Y P I p 5 t I U l m x F R p E Q R l d N F Y R N o Y c O y S K S t F I f o e j x x R h 2 P 3 f F 1 2 s j 1 1 / A J f A S f g M o K s e U n G D t Z 5 N F I I H E l a + 6 c c 9 + + 4 y d S p O i 6 N z X L v n P 3 3 v 2 t B 8 7 D R 4 + f b O / s P j 1 J 4 0 w z 3 m O x j P W Z D y m X Q v E e C p T 8 L N E c I l / y U / / 8 f c m f T r h O R a w + 4 T T h g w h C J U a C A R p o u F v 7 0 f D 8 K J 8 W X S 8 C H D O Q + V H R z L 0 q c q 5 5 U J T 0 V V G 0 2 q W i C q d R G e o o D w r q + b E M 0 m l k D m M z x K 4 n + Q j 7 W R N b a 9 z L i U E / 7 1 N P g S + B e j I O a T L E y q O 5 a u x p E Y 6 x N T s G d C E l t i d V 9 M U i Q O e L A S 4 L E 8 H U + T e J r m f d z / P R 7 q b y S q P r D V X S 9 j 8 1 t b G y / 5 n D Q 8 i G O 3 V 3 z 6 2 E r i u d u V J / + 8 3 p J l + + O 8 d m G 9 5 4 Q c y y i C t k E t K 0 3 3 E T H O S g U T D J C 8 f L U p 4 A O 4 e Q 9 4 2 q I O L p I K 9 2 p a A N g w R 0 F G v z K a Q V u u i R Q 5 S W l R r L 8 u e l q 1 w J 3 s r x C w V a w / R W 0 o + K 5 f s q E G U S h Y 4 v l 1 E E P 5 O g r 8 p 9 X m w K R 6 8 H u V B J h i h U O G t q l E m K M S 2 f D w 2 E 5 g z l 1 C j A t D B z o W w M G h i a R + a Y i X d W 5 7 u u n O z v d Q 7 2 D j 6 6 9 c N 3 Z C Z b 5 D l 5 Q Z q k Q 1 6 R Q / K B H J M e Y d a R F V q J d W F / t W /For two probability measures µ and µ such that µ ≪ µ, the Kullback-Leibler (KL) divergence and Fisher information (or Fisher divergence) of µ with respect to µ are defined, respectively, asKL(µ|| µ) := µ log µ µ and FI(µ|| µ) := µ ∇ log µ µ 2 2.Both divergences are equal to zero if and only if µ = µ. KL divergence is a common metric for quantifying the difference of one distribution with respect to another. Fisher information has been used for analyzing the stationarity of sampling algorithms [3, 65] .We analyze PnP-DM via a continuous-time perspective, leveraging the interpolation techniques introduced for Langevin Monte Carlo [67, 3, 65] . We assume that the likelihood step ( 5) can be implemented exactly and the prior step ( 6) involves running the reverse diffusion process (7) with an approximated score function s t ≈ ∇ log p t := ∇ log p( • ; σ(t)). Let ν X 0 be the distribution of the initialization x (0) . Let ν Z k and ν X k+1 be the distributions of z (k) and x (k+1) at the k th iteration. Recall that the stationary distributions are π X and π Z . Our analysis is concerned with two continuous-time processes: (1) the non-stationary process from ν X 0 , a non-stationary initialization, to ν X K where ( 7) is run with the approximated score function s t and (2) the stationary process that alternates between stationary distributions π X and π Z . These two processes are the interpolation PnP-DM in non-stationary and stationary states and define continuous transitions over discrete iterations. A conceptual illustration of the two processes is provided in Figure 3 with the exact formulations in Appendix A. Now we present our main result: Theorem 3.1. Consider running K iterations of PnP-DM with ρ k ≡ ρ > 0 and a score estimate s t ≈ ∇ log p t := ∇ log p( • ; σ(t)). Let t * > 0 be such that σ(t * ) = ρ and δ := inf t∈[0,t * ] v(t) where v(t) := s(t) 2 σ(t)σ(t). Define ν τ and π τ as the distributions at time τ of the non-stationary and stationary process, respectively. Then, for over K iterations of PnP-DM, or equivalently over τ ∈ [0, T ] with T := K(t * + 1), we haveEQUATION)where we assume that the score estimation error ϵ score :=t * +1 1 v(τ ) 2 E πτ ∥s τ -∇ log p τ ∥ 2 2 dτ < ∞.The proof is provided in Appendix A. This theorem states that the average distance (measured by Fisher information) of the non-stationary process with respect to the stationary process over K iterations of PnP-DM goes to zero at a rate of O(1/K) under certain conditions up to the score approximation error. Note that our theory only requires L 2 -accurate score estimate under the measure π τ , which is a relatively weaker condition than the common L ∞ -accurate score estimateStd. deviation ( ) [17] .< l a t e x i t s h a 1 _ b a s e 6 4 = " f K D p b 1 f x 6 i 7 j C 7 X a R Y F k 3 Y S v l h 4 = " > A A A M Y H i c j Z Z b b 9 s 2 F M f V 7 p Y l 2 d p s b 9 s L k S J A B 2 S Z X e z y N K y J 3 S Y B g s D L p S 1 a d Q E l U x I R S V R I K r Y j 6 G P s d f t c e 9 0 n 2 a E k H k m x N 0 y A 7 a P f O Y e H / P N Q l p f F X O n B 4 K 8 H D z / 4 8 K O P P 1 n 7 d H 1 j 8 7 P P H z 3 e + u K V E r n 0 2 a U v Y i H f e F S x m K f s U n M d s z e Z Z D T x Y v b a u x 4 Z / + t b J h U X 6 Y V e Z O x 9 Q s O U B 9 y n G t A 7 N 6 E 6 U k F x P i 6 v H j 8 Z 7 A 2 q i y w b w 8 Z 4 4 j T X 5 G p r 4 x d 3 K v w 8 Y a n 2 Y 6 r U u + E g 0 + 8 L K j X 3 Y 1 a u u 7 l i G f W v a c g K T 4 h r T T 1 V k v r a I Z k U A V N m Y j T + 9 i a n M d c L A i E x U 7 1 M m q h A p L q T 6 c X g 8 g S V U 2 K m T 9 Q i 8 U S 8 l G W w T Y K s W y p T o S O e h r 3 A l P s s k N R v h / d F A k 5 t h y W B k G T 4 3 b N d w r S / 1 8 t N u C + F U b W 0 u Z a I U N I s W v Q l S D q z 4 R o W 7 B M o M O 0 F + V y 3 o 4 V S 5 B k B V G 1 W f 4 F V A e 7 P y z o 0 4 G E u 7 y l n x N E C l m B j Y C H b p m V S d n O z 3 Z c r D o X k O o L 5 l y Y U 7 9 X q s P K + 2 I b Z a e u I C c n u p b I 0 T 2 B t S T M X L U h E b x m J m d Z M E u N l c s U y k z z W X I p Z u z 2 Q G b L U B L N V 7 Q I o j 6 m c l / 8 z P o r M 0 b m 3 H r m Y q i h O m 7 l O B c x x S k y c I j w l t B 6 n l z K v j m K 5 p L 9 U g S r X 1 3 d 2 d s h E C h F A c S u P X s S s y G L K o Q 4 c B 2 O Q n 5 v G 2 C V T F v C U G 0 E A S p H Q d N 1 N 2 a x J L u D 0 Z E J V / r K Y d G 5 6 U e 0 g Z T F u 7 V 5 M 8 1 s W F 4 3 R 8 8 Y s S W h Z n F Q / P Y 8 v p I h B 6 U V Z j N D s R c A j I U + y u v p + a / d H g W V D 7 V H 1 U + s k M r N Z Q o J U M H / z S Q t 3 1 7 W N 6 + 6 W t Y P K M O H g M j K L r L C P M s D u L j j K s h N H 5 / 8 S R + c Y F / N k K Q g Y + G F 7 w O C 6 j l t V t K 5 H m o A V 1 e p C T Y D K s 6 U A Y D g T n g Z L f m D o h 8 2 f o 8 P c 4 B J E 2 E 5 d h B b f M Y k Y b I s D 3 o 4 C t s V T F S B 2 o R E t H 6 e j 0 1 O q d F G 4 m s 1 1 U d 2 X v 7 m A y m b r f D i + m s L D u s 4 4 M U o U J 1 f t c J q n C 1 J r 0 a T 0 H o T V d n k J V I A v i u v 1 L P I Q + R b 5 i J h F D F F k U Y R I W a Q Q z S 1 q e y G 0 K E S 0 s G i B 6 M 6 i O 0 T S I o n o x q I b R K l F b Z f m F u W I T i w 6 Q X R r 0 S 2 i m U U z R C J h I c p Y 3 6 E P D p 5 u f d V d 2 e k T Y S h 8 B p b t 2 9 h 9 D B t b N E Z 0 b N H x q o 0 l M / j P I N A n T P I 5 a l 6 3 U i W 7 b S I r N L o W f d d d x 3 X X b 7 2 l c h H V W A p s S J r x K T O W 3 e o m U 9 G 0 6 d c L a F H s 1 Q t c 3 H k X n 7 c y d H E r x c s u f o n 4 q I u P E J 9 2 8 S n i w y 4 + b G X u 4 m P E k y 6 e I D 7 r 4 j P E l 1 1 8 + R 8 T r O X p v G 1 V 4 j 5 V T C s i A g J / 2 R 6 8 W X 5 T p 4 / q X M 8 r R l g f 0 Z l F L x C 9 Q A U Q n V r 0 F t F b f L L A 7 j a v P E 0 9 I M 1 8 j T V q F e 0 7 O p r 2 H a 2 q B 3 3 H A X Z G + 1 A a e 1 Y c r 7 v X Rassumption in prior analysis of sampling methods involving score estimates [5, 65] . This result resembles the first-order stationarity for Langevin Monte Carlo [3] . Unlike the non-asymptotic analysis in [77] , we utilize the average Fisher information instead of the total variation distance, enabling us to obtain an explicit convergence rate. Here δ is the infimum of the diffusion coefficient along the reverse diffusion in (7) ; see further discussions on the role of δ in Appendix A.3. Our theory shows that the accurate implementations of the two sampling steps lead to a sampler that provably converges to the stationary process that alternates between the two target stationary distributions. A test image is randomly generated from the prior (see top left of Figure 4 ), and the measurement is calculated according to (1) with n ∼ N (0, 0.01 2 I). We compare our method with the popular DM-based method DPS [17] . We draw 1,000 samples and visualize the empirical mean and per-pixel standard deviation for both algorithms. Compared with the true posterior (second column), we find that the both methods accurately estimate the mean. However, the standard deviation image estimated by DPS significantly deviates from the ground truth. In contrast, our standard deviation image matches the ground truth in terms of both absolute magnitude and spatial distribution. These results highlight the accuracy of our method over DPS by taking a more principled Bayesian approach. 
benchmark experiments: 
 Dataset and inverse problems We test our proposed algorithm and several baseline methods on 100 images from the validation set of the FFHQ dataset [37] for five inverse problems: (1) Gaussian deblur with kernel size 61×61 and standard deviation 3.0, (2) Motion deblur with kernel size 61×61 and intensity of 0.5, (3) Super-resolution with 4× downsampling ratio, (4) the coded diffraction patterns (CDP) reconstruction problem (nonlinear) in [10, 50] (phase retrieval with a phase mask), and (5) the Fourier phase retrieval (nonlinear) with 4× oversampling. We add i.i.d. Gaussian noise to all the simulated measurements y. In particular, i.e. n ∼ N (0, σ 2 y I). For all problems except for Fourier phase retrieval, the noise standard deviation is set as σ y = 0.05. Due to the severe ill-posedness of Fourier phase retrieval, we consider a smaller noise standard deviation σ y = 0.01.We consider four variants of DMs as plug-in priors for our method, namely VP-SDE (VP) [33] , VE-SDE (VE) [62] , iDDPM [51] , and EDM [36] . We compare our method with various baselines, including (1) optimization-based methods: PnP-ADMM [13] , DPIR [79] ; (2) conditional DMs: DDRM [38] , DPS [18] ; and (3) SGS-based method: PnP-SGS [19] , DPnP [77] . For fair comparison, we use the same pre-trained score function checkpoint for all DM-based methods. Since the pre-trained score function was trained with the DDPM formu-lation (VP-SDE) [33] , we convert it to the EDM formulation by applying the VP preconditioning [36] . We use the Peak Signal-to-Noise Ratio (PSNR), the Structural Similarity Index Measure (SSIM), and the Learned Perceptual Image Patch Similarity (LPIPS) distance for quantitative comparison. For each sampling method, we draw 20 randoms samples, calculate their mean, and report the metrics on the mean image. More experimental details are provided in Appendices B, C, D. L M H H v a 0 G v e 3 U h k S 3 A 0 e U = " > A A A M a H i c j Z b f b 9 s 2 E M f V 7 l e X Z F u y P Q z D X o g E A T o g z e x i 7 Z 6 G N X H a J E A Q e G n S F p 2 6 g J I p i Y g k K i Q V 2 1 H 1 p + x 1 + 5 v 2 L + y v 2 F E S T 1 L s D T N g 6 / S 5 O x 7 5 5 V G W l 8 V c 6 c H g r 3 v 3 P / j w o 4 8 / e f D p y u r a Z 5 9 / s b 7 x 5 S s l c u m z C 1 / E Q r 7 x q G I x T 9 m F 5 j p m b z L J a O L F 7 L V 3 N T L + 1 z d M K i 7 S c z 3 P 2 L u E h i k P u E 8 1 o M v 1 j f e u R 2 X h z r y k f G R + 3 1 + u b w 1 2 B 9 W H L B r D x t h y m s / 4 c m P 1 Z 3 c i / D x h q f Z j q t S v w 0 G m 3 x V U a u 7 H r F x x c 8 U y 6 l / R k B W e E F e a e q o k 9 W e b Z F I E T J k J 0 v j R d U 5 j r u c E Q m K m e p k 0 U Y F I d S f T i 8 H l C S o n J K E 6 I m q e e C J e y D L Y J k H W D Z W p 0 B F P w 1 5 g y n 0 W S O q 3 w / s i A a e 2 w 5 J A S D L 8 / v E O Y d r f 7 e U m 3 J f C q F v a X E t E K G k W z f s S J J 3 Z c A 0 L 9 g k U m P S C f K 7 b 0 U I p 8 o w A q j a t v 8 C q A P d n Z R 0 a 8 D C X d 5 Q z 4 m g B S 7 A x s J B N 0 z o p u 7 7 e 7 M s V h 0 J y H c H 8 S x O K 9 2 p 5 W H l X b M P s t H X E h G R 3 U l m a J 7 C 2 p J m L F i S i N 4 z E T G s m i f E y u W S Z S R 5 r L s W 0 3 R 7 I D F l q g t m y d g G U x 1 T O y v 8 Z H 0 X m C N 1 Z j 5 x P V B S n z V w n A u Y 4 I S Z O E Z 4 S W o / T S 5 l V R 7 J c 0 F + q Q J U r K 9 v b 2 2 Q s h Q i g u J V H z 2 N W Z D H l U A e O g z H I T 0 1 j 7 J A J C 3 j K j S A A p U h o u u K m b N o k F 3 B 6 M q E q f 1 m M O z e 9 q H a Q s j h o 7 V 5 M c y 2 L 8 8 b o e W O W J L Q s T q p L z + M L K W J Q e l 4 W I z R 7 E f B I y J O s r r 7 X 2 v 1 R Y N l Q e 1 R d a p 1 E Z j Z L S J A K 5 m + + a e H u u L Z x 3 Z 2 y d l A Z J h x c R m a R 1 V c V w O a F 7 g 4 4 y r I T R 2 f / E k d n G B f z Z C E I G P h h e 8 D g u o 5 b V r S u R 5 q A J d X q Q k 2 A y r O F A G A 4 E 5 4 G C 3 5 g 6 I f N n 6 H D 3 O A S R N h O X Y Q W 3 z K J G G y L A 9 6 O A r b F E x U g d q E R L T 9 I R 6 e n V O m i c D W b 6 a K 6 L 3 9 z A Z X N 1 v l w f D W F h 3 W d c W K U K E 4 u 2 + E 0 T + e k 1 q J J 6 T 0 I q + 3 y E q g A P x T X 6 1 n k I f I t 8 h E x i x i i y K I I k b J I I Z p Z 1 P Z C a F G I a G 7 R H N G t R b e I p E U S 0 b V F 1 4 h S i 9 o u z S 3 K E Z 1 Y d I L o x q I b R F O L p o h E w k K U s b 5 D H x w 8 3 f q q u 7 L T J 8 J Q + A 4 s 2 7 O x e x h 2 Y N E B o m O L j p d t L J n C f w a B P m G S z 1 D z u p U q 2 W 0 T W a H R N e + 7 b j u u 2 3 7 r L Z S L q M Z S Y E P S l E + Y s e x W N 5 m K p k 2 / n k O L Y q + e 4 + J e d v H L V o Y u b q V 4 0 c U v E B 9 1 8 R H i 0 y 4 + R X z Y x Y e t z F 1 8 j H j c x W P E Z 1 1 8 h v i i i y / + Y 4 K 1 P J 2 3 r U r c h 4 p p R U R A 4 C / b g z f M 7 + r 0 U Z 3f C q F v a X E t E K G k W z f s S J J 3 Z c A 0 L 9 g k U m P S C f K 7 b 0 U I p 8 o w A q j a t v 8 C q A P d n Z R 0 a 8 D C X d 5 Q z 4 m g B S 7 A x s J B N 0 z o p u 7 7 e 7 M s V h 0 J y H c H 8 S x O K 9 2 p 5 W H l X b M P s t H X E h G R 3 U l m a J 7 C 2 p J m L F i S i N 4 z E T G s m i f E y u W S Z S R 5 r L s W 0 3 R 7 I D F l q g t m y d g G U x 1 T O y v 8 Z H 0 X m C N 1 Z j 5 x P V B S n z V w n A u Y 4 I S Z O E Z 4 S W o / T S 5 l V R 7 J c 0 F + q Q J U r K 9 v b 2 2 Q s h Q i g u J V H z 2 N W Z D H l U A e O g z H I T 0 1 j 7 J A J C 3 j K j S A A p U h o u u K m b N o k F 3 B 6 M q E q f 1 m M O z e 9 q H a Q s j h o 7 V 5 M c y 2 L 8 8 b o e W O W J L Q s T q p L z + M L K W J Q e l 4 W I z R 7 E f B I y J O s r r 7 X 2 v 1 R Y N l Q e 1 R d a p 1 E Z j Z L S J A K 5 m + + a e H u u L Z x 3 Z 2 y d l A Z J h x c R m a R 1 V c V w O a F 7 g 4 4 y r I T R 2 f / E k d n G B f z Z C E I G P h h e 8 D g u o 5 b V r S u R 5 q A J d X q Q k 2 A y r O F A G A 4 E 5 4 G C 3 5 g 6 I f N n 6 H D 3 O A S R N h O X Y Q W 3 z K J G G y L A 9 6 O A r b F E x U g d q E R L T 9 I R 6 e n V O m i c D W b 6 a K 6 L 3 9 z A Z X N 1 v l w f D W F h 3 W d c W K U K E 4 u 2 + E 0 T + e k 1 q J J 6 T 0 I q + 3 y E q g A P x T X 6 1 n k I f I t 8 h E x i x i i y K I I k b J I I Z p Z 1 P Z C a F G I a G 7 R H N G t R b e I p E U S 0 b V F 1 4 h S i 9 o u z S 3 K E Z 1 Y d I L o x q I b R F O L p o h E w k K U s b 5 D H x w 8 3 f q q u 7 L T J 8 J Q + AM P s t H X E h G R 3 U l m a J 7 C 2 p J m L F i S i N 4 z E T G s m i f E y u W S Z S R 5 r L s W 0 3 R 7 I D F l q g t m y d g G U x 1 T O y v 8 Z H 0 X m C N 1 Z j 5 x P V B S n z V w n A u Y 4 I S Z O E Z 4 S W o / T S 5 l V R 7 J c 0 F + q Q J U r K 9 v b 2 2 Q s h Q i g u J V H z 2 N W Z D H l U A e O g z H I T 0 1 j 7 J A J C 3 j K j S A A p U h o u u K m b N o k F 3 B 6 M q E q f 1 m M O z e 9 q H a Q s j h o 7 V 5 M c y 2 L 8 8 b o e W O W J L Q s T q p L z + M L K W J Q e l 4 W I z R 7 E f B I y J O s r r 7 X 2 v 1 R Y N l Q e 1 R d a p 1 E Z j Z L S J A K 5 m + + a e H u u L Z x 3 Z 2 y d l A Z J h x c R m a R 1 V c V w O a F 7 g 4 4 y r I T R 2 f / E k d n G B f z Z C E I G P h h e 8 D g u o 5 b V r S u R 5 q A J d X q Q k 2 A y r O F A G A 4 E 5 4 G C 3 5 g 6 I f N n 6 H D 3 O A S R N h O X Y Q W 3 z K J G G y L A 9 6 O A r b F E x U g d q E R L T 9 I R 6 e n V O m i c D W b 6 a K 6 L 3 9 z A Z X N 1 v l w f D W F h 3 W d c W K U K E 4 u 2 + E 0 T + e k 1 q J J 6 T 0 I q + 3 y E q g A P x T X 6 1 n k I f I t 8 h E x i x i i y K I I k b J I I Z p Z 1 P Z C a F G I a G 7 R H N G t R b e I p E U S 0 b V F 1 4 h S i 9 o u z S 3 K E Z 1 Y d I L o x q I b R F O L p o h E w k K U s b 5 D H x w 8 3 f q q u 7 L T J 8 J Q + AL M H H v a 0 G v e 3 U h k S 3 A 0 e U = " > A A A M a H i c j Z b f b 9 s 2 E M f V 7 l e X Z F u y P Q z D X o g E A T o g z e x i 7 Z 6 G N X H a J E A Q e G n S F p 2 6 g J I p i Y g k K i Q V 2 1 H 1 p + x 1 + 5 v 2 L + y v 2 F E S Tm z C 1 / E Q r 7 x q G I x T 9 m F 5 j p m b z L J a O L F 7 L V 3 N T L + 1 z d M K i 7 S c z 3 P 2 L u E h i k P u E 8 1 o M v 1 j f e u R 2 X h z r y k f G R + 3 1 + u b w 1 2 B 9 W H L B r D x t h y m s / 4 c m P 1 Z 3 c i / D x h q f Z j q t S v w 0 G m 3 x V U a u 7 H r F x x c 8 U y 6 l / R k B W e E F e a e q o k 9 W e b Z F I E T J k J 0 v j R d U 5 j r u c E Q m K m e p k 0 U Y F I d S f T i 8 H l C S o n J K E 6 I m q e e C J e y D L Y J k H W D Z W p 0 B F P w 1 5 g y n 0 W S O q 3 w / s i A a e 2 w 5 J A S D L 8 / v E O Y d r f 7 e U m 3 J f C q F v a X E t E K G k W z f s S J J 3 Z c A 0 L 9 g k U m P S C f K 7 b 0 U I p 8 o w A q j a t v 8 C q A P d n Z R 0 a 8 D C X d 5 Q z 4 m g B S 7 A x s J B N 0 z o p u 7 7 e 7 M s V h 0 J y H c H 8 S x O K 9 2 p 5 W H l X b M P s t H X E h G R 3 U l m a J 7 C 2 p J m L F i S i N 4 z E T G s m i f E y u W S Z S R 5 r L s W 0 3 R 7 I D F l q g t m y d g G U x 1 T O y v 8 Z H 0 X m C N 1 Z j 5 x P V B S n z V w n A u Y 4 I S Z O E Z 4 S W o / T S 5 l V R 7 J c 0 F + q Q J U r K 9 v b 2 2 Q s h Q i g u J V H z 2 N W Z D H l U A e O g z H I T 0 1 j 7 J A J C 3 j K j S A A p U h o u u K m b N o k F 3 B 6 M q E q f 1 m M O z e 9 q H a Q s j h o 7 V 5 M c y 2 L 8 8 b o e W O W J L Q s T q p L z + M L K W J Q e l 4 W I z R 7 E f B I y J O s r r 7 X 2 v 1 R Y N l Q e 1 R d a p 1 E Z j Z L S J A K 5 m + + a e H u u L Z x 3 Z 2 y d l A Z J h x c R m a R 1 V c V w O a F 7 g 4 4 y r I T R 2 f / E k d n G B f z Z C E I G P h h e 8 D g u o 5 b V r S u R 5 q A J d X q Q k 2 A y r O F A G A 4 E 5 4 G C 3 5 g 6 I f N n 6 H D 3 O A S R N h O X Y Q W 3 z K J G G y L A 9 6 O A r b F E x U g d q E R L T 9 I R 6 e n V O m i c D W b 6 a K 6 L 3 9 z A Z X N 1 v l w f D W F h 3 W d c W K U K E 4 u 2 + E 0 T + e k 1 q J J 6 T 0 I q + 3 y E q g A P x T X 6 1 n k I f I t 8 h E x i x i i y K I I k b J I I Z p Z 1 P Z C a F G I a G 7 R H N G t R b e I p E U S 0 b V F 1 4 h S i 9 o u z S 3 K E Z 1 Y d I L o x q I b R F O L p o h E w k K U s b 5 D H x w 8 3 f q q u 7 L T J 8 J Q + A 4 s 2 7 O x e x h 2 Y N E B o m O L j p d t L J n C f w a B P m G S z 1 D z u p U q 2 W 0 T W a H R N e + 7 b j u u 2 3 7 r L Z S L q M Z S Y E P S l E + Y s e x W N 5 m K p k 2 / n k O L Y q + e 4 + J e d v H L V o Y u b q V 4 0 c U v E B 9 1 8 R H i 0 y 4 + R X z Y x Y e t z F 1 8 j H j c x W P E Z 1 1 8 h v i i i y / + Y 4 K 1 P J 2 3 r U r c h 4 p p R U R A 4 C / b g z f M 7 + r 0 U Z 3F + q Q J U r K 9 v b 2 2 Q s h Q i g u J V H z 2 N W Z D H l U A e O g z H I T 0 1 j 7 J A J C 3 j K j S A A p U h o u u K m b N o k F 3 B 6 M q E q f 1 m M O z e 9 q H a Q s j h o 7 V 5 M c y 2 L 8 8 b o e W O W J L Q s T q p L z + M L K W J Q e l 4 W I z R 7 E f B I y J O s r r 7 X 2 v 1 R Y N l Q e 1 R d a p 1 E Z j Z L S J A K 5 m + + a e H u u L Z x 3 Z 2 y d l A Z J h x c R m a R 1 V c V w O a F 7 g 4 4 y r I T R 2 f / E k d n G B f z Z C E I G P h h e 8 D g u o 5 b V r S u R 5 q A J d X q Q k 2 A y r O F A G A 4 E 5 4 G C 3 5 g 6 I f N n 6 H D 3 O A S R N h O X Y Q W 3 z K J G G y L A 9 6 O A r b F E x U g d q E R L T 9 I R 6 e n V O m i c D W b 6 a K 6 L 3 9 z A Z X N 1 v l w f D W F h 3 W d c W K U K E 4 u 2 + E0 T + e k 1 q J J 6 T 0 I q + 3 y E q g A P x T X 6 Q R D H M N d v P s X D o t / x Q = " > A A A C Y H i c b V F N S N A E N E r x o / u p N L F S F Q S k e p J R C e F a w t m C K T a Y u m i k Q t I X / D q / t r / S t U H / p g c M z D z N s E N x g E X t L y y u l Z Z z Y t a r t f r O n U l z T V m H p i L V v R g M E y x D n I U r J d p B j I W r B s / X Y g v T h q f q F k c Z s Y K p w C m i l K J K A j y Y p D A K h o j a A U T + I s k n J I G m e L o e c R Y O U p I p p A K M u Q + D D P s F a O R U s N K L c s M y o E w Z P e W K p D M I v J q X f t M r A T J t n J / o v d K E A a M K x Z w s O e + N x X q x A a x j a a y n K n B Z k L D p n V U R l y A f i u s x R m J z C y H J G r e R S S T P x v A a Z R j C w B q r n N x a e P o I G i / R P P J h O t I o a Y b v V v j l u n F M s + Q f X J A D k l I T s g u S L X p E M o y c g + S C f z r d b c a t u / b f V d a Y z u Q G t P b a K U Q = = < / l a t e x i t > std < l a t e x i t s h a 1 _ b a s e 6 4 = " V H O V N 8 p L e o + Q 2 t S j A W Z l P Q L H p W Y = " > A A A C e H i c b V F N T x s x E H W 2 h U J K I S l H L h Z R B B x I d x G i P V W I X n o E q S G R s l E 0 d r z B w v Y u 9 m y b a N l 7 f 0 2 v 7 V / p X + k J 5 + P Q f I x k + f m 9 G c n v D c u U d B i G f y v B q 9 db 2 2 9 2 d q t v 9 9 7 t H 9 T q 7 + 9 d m l s u 2 j x V q e 0 y c E J J I 9 o o U Y l u Z g V o p k S H P X 6 Z 6 p 3 v w j q Z m m 8 4 y U R f w 8 j I R H J A T w 1 q x 8 8 x A 1 v E T B f j s j y f 3 8 Results: linear problems A quantitative comparison is provided in Table 1 . PnP-DM generally outperforms the baseline methods and that the VE and EDM variants consistently outperform the other two variants on these linear problems. Figure 5 contains visual examples for the motion deblur problem (see Appendix F.2 for the other two linear problems). PnP-DM provides high-quality reconstructions that are both sharp and consistent with the ground truth image. We also provide an uncertainty quantification analysis based on pixel-wise statistics in Figure 6 PnP-SGS, our approach captures a higher percentage (97.46%) of ground truth pixels than the baselines (96.20% and 88.77%). If the true posterior were truly Gaussian, 99% of the ground-truth pixels should lie within the 3-sigma CI; however, as the posterior is not Gaussian with a DM-based prior, we do not necessarily expect to reach 99% coverage.8 f Y g 3 4 4 J L C 4 b A c 1 B p h K 5 w V X Q f R A j T I o m 4 H 9 c r n e J j y X A u D X I F z v S j M s F + A R c m V K K t x 7 k Q G / B F G o u e h A S 1 c v 5 i Z K W n T M 0 O a p N Y f g 3 T G / j 9 R g H Z u o p n v n H 1 y V Z u S G z X x Z M B a m G w U m S 6 X 3 6 u E z h V K m / 5 Y Z h F Y r s C O y 2 p z y R Q m n / q F N F m O K M 1 o b i r J F c W U T r d A h 9 I K j m r i A X A r f S 6 U P 4 A F j n 5 X V Z 9 4 t J r v O r i / a E V XResults: nonlinear problems We provide a quantitative comparison in Table 2 . For the CDP reconstruction problem, PnP-DM performs on par with DPS but outperforms other SGS-based methods. We then consider the Fourier phase retrieval (FPR) problem, which is known to be a challenging nonlinear inverse problem. One challenge lies in its invariance to 180 • rotation, so the posterior distribution have two modes, one with upright images and another with 180 • -rotated images, that equally fit the measurement. To increase the chance of getting properly-oriented reconstructions, we run each algorithm with four different random initializations and report the metrics for the best run, following the practice in [18] . We find that PnP-DM significantly outperforms the baselines on this highly ill-posed inverse problem. As shown in Figure 7 (a), our method can provide high-quality reconstructions for both orientations, while the baseline methods fail to capture at least one of the two modes. We further run our method for a test image with 100 different random initialization and collect reconstructions in both orientations that are above 28dB in PSNR (90 out of 100 runs). The percentage of upright and rotated reconstructions are visualized by the pie chart in Figure 7 (b) . With a prior on upright face images, our method generate mostly samples with the upright orientation. Nevertheless, it can also find the other mode that has an equal likelihood, demonstrating its ability to capture multi-modal posterior distributions.Problem setup We finally validate PnP-DM on a real-world nonlinear imaging inverse problem: black hole imaging (BHI) (see Appendix B for more details). A visual illustration of BHI is provided in Figure 8 (a). This BHI inverse problem is severely ill-posed. Even with an Earth-sized telescope, only a small fraction of the Fourier frequencies of the target black hole can be measured (region within the red box); in reality, this region is further subsampled with a highly sparse pattern (black lines). Additionally, the atmospheric noise causes nonlinearity of this BHI problem that sometimes results in a multi-modal posterior distribution of the reconstructed image [63] . Here we demonstrate the effectiveness of PnP-DM in capturing a multi-modal posterior distribution. For brevity, we restrict our choice of diffusion models in PnP-DM to EDM and use DPS as the baseline.Results on simulated data We use the simulated data from [63] where the measurements are generated assuming that the ground-truth black hole image were at the location of the Sagittarius A * black hole. Figure 8 (b) visually compares the results obtained by PnP-DM and DPS. We use the t-SNE method [66] to cluster the generated samples (100 for each method) and identify two modes in the samples generated by PnP-DM and three modes in those generated by DPS. We visualize the mean and three samples for each image mode. A metric for quantifying the degree of data mismatch is labeled on the top right corner of each image. As illustrated by both the mean and sample images, PnP-DM successfully captures the two modes previously identified for this dataset [63] . Note that PnP-DM generates high-fidelity samples from both modes with sharp details of the flux ring, and its samples from "Mode 1" align well with the ground truth image. In contrast, two out of the three modes sampled by DPS fail to exhibit a meaningful black hole structure and do not correspond with the observed measurements, as indicated by the significantly larger data mismatch values.Results on real data Finally, we apply PnP-DM to the real M87 black hole data from April 6 th , 2017 [21] , with the results shown in Figure 1 . By leveraging an expressive DM-based image prior, PnP-DM generates high-quality samples that are both visually plausible and consistent with the ring diameters observed in the official EHT reconstruction. These results highlight the robustness and effectiveness of our method in tackling a highly ill-posed real-world inverse problem. 
conclusion: 
 We have introduced PnP-DM, a posterior sampling method for solving imaging inverse problems. The backbone of our method is a split Gibbs sampler that iteratively alternates between two steps that separately involve the likelihood and prior. Crucially, we establish a link between the prior step and a general DM framework known as the EDM formulation. By leveraging this connection, we seamlessly integrate a diverse range of state-of-the-art DMs as priors through a unified interface.Experimental results demonstrate that our method outperforms existing DM-based methods across n n 3 2 + / P i L 1 0 q k 0 m V n r g i F f O t Q x U I e s z P N d c j e J p L R y A n Z G + d i V N j f X D G p u I h P 9 T x h 7 y P q x 9 z j L t W A z p e / t d N 4 w q Q j q c s y + z K l k z t / 8 v P l 1 U F / s L W 9 / X S b D P q b m 8 O d j Q E M d n a 2 N w c 7 Z N g f l J 9 V q / 6 M z x 8 / + s W e C D e N W K z d k C r 1 2 3 C Q 6 P c Z l Z q 7 I c u X 7 F S x h L o X 1 G e Z I 8 S F p o 7 K S f V Z I 4 k U H l P F 0 m n 4 H a w g 5 H p O w C V k q h N J I + W J W L c i n R B M j q B y Q i K q A 6 L m k S P C h a g C m y C I u q I y F j r g s d 9 both linear and nonlinear inverse problems, including a nonlinear and severely ill-posed black hole interferometric imaging problem.x j L n L P N C n m d 4 V E R i 1 m Z Z 4 Q p L h 9 x v r h G m 3 3 4 m N u C t F o X t u Y g 0 R v q R J M O 9 K E L V W w z V s 2 C W Q Y N J x c r l u Z v O l S B M C q C x n d 4 N l A u 7 O 8 s r V 4 3 4 q b y h X i K M F b M H 4 w E Z W i q a KZ F N M + O y 0 v H 4 g o p Q l B 6 n m c j H H Y 8 4 J G Q R k m V f b c Z d 2 e B b U P u U X m p d B J J U S w h Q S p Y f / G N M 3 v d N o 1 r r + e V g U o / 4 m A q Z B Z J d V U e F M + 3 1 8 G Q 5 y 0 / O r v D j 8 7 Q L + T R g h M w s E N 5 Y M B 1 5 X d b 0 i o f q R 1 u y V Y l q h 1 U m i w 4 A M O V 8 N h b s A N D O x R / h o b i B r c g / G b p w j f 4 m k n E M D b Y 4 8 0 s M D Z 4 o j z E N j S i 4 f v x 6 O S E K p 1 l t m Y z n Z X 3 + e 8 2 o L w u n Q v H V 1 N 4 W F c R x 4 U S 2 f F 5 M 5 3 m 8 Z x U W t Q h n Q d h W S 4 n g g z w Q 3 G / j k E O I t c g F x E z i C E K D A o Q K Y M U o p l B T S / 4 B v m I 5 g b N E V 0 b d I 1 I G i Q R X R p 0 i S g 2 q O n S 1 K A U 0 b F B x 4 i u D L p C N D V o i k h E z E c Z q z u 0 w c H T j a 2 8 y 1 t 9 I g o K 3 4 F h u 8 Z 3 F 9 3 2 D d p H d G T Q 0 W 2 F J V P 4 z y D Q J 0 z y G W p e t V I p u 2 k i I z S a 5 l 3 T d c t 0 3 W 2 9 h X Q B 1 Z g K x h A 0 5 R N W j E y p 6 0 h F 4 7 p f T 6 F F s V d P c X O v 2 v h V I 0 M b N 1 K 8 a O M X i A / b + B D x S RK + j Y a r f F 1 R q 7 k W s H E w z x V L q X d G A F a 4 Q V 5 q 6 q i T 1 t U 1 S K X y m z H R o 9 O 1 1 R i O u c w I h E V O 9 T B oO + u R + U y F U d L M d S Z g j j N i 4 h T h C a H 1 O L 2 U R X U A y y X 9 p f J V O R h s b 2 + T i R T C h + J W H p 1 H r E g j y q E O H A d j k J + b x t g h M + b z h B t B A E o R 0 2 Q w T d i 8 S S 7 g 9 K R C V f 6 y m H Q e e l H t I G V x 0 N q 9 m O Z e F u e N 0 f N G L I 5 p W Z x U t 5 7 H E 1 J E o H R e F m M 0 e x H w S s j i t K 6 + 1 9 r 9 U W D Z U H t c 3 W q d R G o 2 S 0 i Q C u Z v / p J i u j O 1 j T v d K W s H l U H M w W V k F m l 9 V z 5 s X j D d A U d Z d u L o 4 l / i 6 A L j I h 4 v B Q E D P 2 w P G F z X c a u K 1 v V I E 7 C i W l 2 o C V B Z u h Q A D G f C E 3 / J D w z 9 s P k L d J g H X I I I 2 q m L w O J b J h G D b b H P 2 1 H A t n i m f M R T a E T L D 5 L x 6 S l V u i i m m i 1 0 U T 2 X H 6 a A y m b r P D i + m s L L u s 4 4 M U o U J 5 f t c J o n O a m 1 a F J 6 L 8 J q u 9 w Y K s A / i u t 1 L X I R e R Z 5 i J h F D F F o U Y h I W a Q Q L S x q e y G w K E C U W 5 Q j u r X o F p G 0 S C K 6 t u g a U W J R 2 6 W Z R R m i E 4 t O E N 1 Y d I N o b t E c k Y h Z g D L W T + i D g 6 d b X / V U d v p E G A p / Q 8 v 2 b O w e h h 1 Y d I D o 2 K L j V R t L 5 v C b Q a B P m O Q L 1 L x u p U p 2 2 0 R W a H T l f d d t xV v f + O z z L z Y f f / l a i U x 6 7 N I T k Z B v X a p Y x B N 2 q b m O 2 N t U M h q 7 E X v j X o + M / 8 2 M S c V F c q H z l L 2 P a Z B w n 3 t U A 3 o 5 u 9 r c H u w N q o s s G 8 P G 2 H a a a 3 z 1 e P 2 X y V R 4 W c w S 7 U V U q V + H g 1 S / L 6 j U 3 I t Y u T b J F E u p d 0 0 D V r h C X G v q q p L U 1 w 5 J p f C Z M t O h 0 b c 3 G Y 2 4 z g m E R E z 1 M m m s f J H o T q Y b g c s V V E 5 J T H V I V B 6 7 I l r K M t g m Q d a M y k T o k C d B L z D h H v M l 9 d r h P R G D U 9 t h i S 8 k G X 7 3 d J c w 7 e 3 1 c m P u S W G 0 L G 2 u J S K Q N A 3 z v g R x Z z Z c w 4 I 9 A g W m v S C P 6 3 a 0 Q I o s J Y C q L e o v s C r A v U V Z h / o 8 y O Q d 5 Y w 4 W s A S b A w s Z M s 0 S s J u b r b 6 c k W B k F y H M P / S h O K z W h 1 W 3 h X b M D t t H T I h 2 Z 1 U l m Q x r C 1 u 5 q I F C e m M k Y h p z S Q x X i Z X L D P O I s 2 l m L f b A 5 k B S 0 w w W 9 U u g L K I y k X 5 P + P D 0 B y Y O + u R + V S F U d L M d S p g j l N i 4 h T h C a H 1 O L 2 U R X U A y y X 9 p f J V u b a 2 s 7 N D x l I I H 4 p b e X Q e s S K N K I c 6 c B y M Q X 5 u G m O X T J n P E 2 4 E A S h F T J O 1 S c L m T X I B p y c V q v K X x b j z 0 I t q B y m L w 9 b u x T T 3 s r h o j J 4 3 Y n F M y + K 0 u v U 8 n p A i A q X z s h i h 2 Y u A V 0 I W p 3 X 1 / d b u j w L L h t q j 6 l b r J F K z W U K C V D B / 8 5 c U k 9 2 J b d z J b l k 7 q A x i D i 4 j s 0 j r u / J h 8 4 L J L j j K s h N H F / 8 S R x c Y F / F 4 K Q g Y + G F 7 w O C 6 j l t V t K 5 H m o A V 1 e p C T Y D K 0 q U A Y D g T n v h L f m D o h 8 1 f o M M 8 4 B J E 0 E 5 d B B b f M o k Y b I t 9 3 o 4 C t s V T 5 S O e Q C N a f p i M z s 6 o 0 k U x 0 W y h i + q 5 / D A B V D Z b 5 8 H x 1 R R e 1 n X G q V G i O L 1 q h 9 M 8 y U m t R Z P S e x F W 2 + X G U A H + U V y v a 5 G L y L P I Q 8 Q s Y o h C i 0 J E y i K F a G F R 2 w u B R Q G i 3 K I c 0 a 1 F t 4 i k R R L R j U U 3 i B K L 2 i 7 N L M o Q n V p 0 i m h m 0 Q z R 3 K I 5 I h G z A G W s n 9 A H B 0 + 3 v u q p 7 P S J M B T + B p b t 2 9 h 9 D D u 0 6 B D R i U U n q z a W z O E 3 g 0 C f M M k X q H n d S p X s t o m s 0 O j K + 6 7 b j u u 2 3 3 p L 5 U K q s R T Y k D T n U 2 Y s u 9 V N p q J J 0 6 8 X 0 K L Y q x e 4 u F d d / K q V o Y t b K V 5 0 8 Q v E x 1 1 8 j P i s i 8 8 Q H 3 X x U S tQ = " > A A A M g n i c j Z b h b t s 2 E M f V r t u 6 1 N v S 7 e O + E A k C d F v m 2 W m T J s C G N X H a J E A Q e G n S F p 2 6 g J I o i Y g k K i Q V 2 x H 0 E H u a f d 0 e Y 2 + z o y R S U u w M M 2 C L + t 0 d j / z f U Z a T R l T I w e C f e / c / e v D x J 5 8 + / G z p U e / z L 7 5 c f v z V G 8 E y 7 p J z l 0 W M v 3 O w I B F N y L m k M i L v U k 5 w 7 E T k r X M 5 U v a 3 1 4 Q L y p I z O U v J h x g H C f W p i y W g i + X v 7 S z x C H c 4 d k l u X 2 X Y u / O n u F h e H f R 3 t p 9 v D 7 b Q o L + x s z n Y 2 V a D Z 0 8 3 N w d o 2 B + U n 1 W r / o w v H j / 6 x f a Y m 8 U k k W 6 E h f h t O E j l h x x z S d 2 I F E t 2 J k i K 3 U s c k N x h 7 F J i R x S o + q y h l D O f C L V 0 H P 0 A K 4 i o n C F w i Y j o R O J Y + C y R r U g n A p P D M P d Q j G W I x C x 2 W D Q X p b A O g q h r z B M m Q 5 o E H c e E u s Q H f Z r p X R a D U e p p k c 8 4 G v 6 4 s Y 6 I d P u d 2 J i 6 n C n d C x 2 r C Q s 4 T s N Z V 4 K 4 t R o q Y c M u g g R e x 8 m l s p k t 4 C x L E a C y n N 0 N l g m o O y 0 q V 5 8 G G b + l n B J H M t i C 9 o G N r K i m S s j V 1 U p X r i h g n M o Q 1 l 8 o V 3 M v F r s V t 8 V W T C 9 b h o R x c i u U J F k M e 4 v r t U i G Q n x N U E S k J B w p K + E L t h l n k a S c T Z r y Q G R A E u V M F r U L o C z C f F r 8 T / 8 w V I f r 1 n 7 4 z B N h l N R r 9 R i s 0 U P K T y C a I F z N 0 w m Z l o e 1 m N O f C 1 8 U S 0 t r a 2 t o z B n z I b m W R 8 4 i k q c R p p A H j o M a o J / r x l h H H v F p Q p U g A D m L c b J k J 2 R S B + d w e l I m S n u R j 1 s 3 H a 9 m k i L f b 8 Y d n / p a 5 G f 1 o G O N S B z j I j 8 u L x 2 L y z i L Q O l Z k Y / M s O M B j 4 Q s T q v s u 8 2 4 O w t s G 3 K P y k u l E 0 t V s R g H q W D 9 6 p v k 9 r q t G 9 d e L y o D 5 k F M w a R k Z m l 1 F T 4 U L 7 D X w V A U L T 8 8 v c M P T 4 1 f R O M 5 J 2 B g h / L A g M r K b 1 H S K h + q H R Z k q x L V D i J L 5 x y A m Z X Q x J + z A z N 2 K P 7 U G N S N 2 Q I L m q W z Q O M b w g 2 G s c Y + b W a B s c a e 8 A 2 2 o R E 1 3 0 9 G J y d Y y D y 3 J Z n K v L w v f r c B F X X p X D i + E s P D u o o 4 V k r k x x f N d J I m M 1 R p U Y d 0 H o R l u Z w Y M s A P N v t 1 N H I M c j V y D S I a E Y N C j U K D h E b C o K l G T S 8 E G g U G z T S a G X S j 0 Y 1 B X C N u 0 J V G V w Y l G j V d m m m U G X S s 0 b F B 1 x p d G z T R a G I Q i 0 l g Z K z u j A 0 O n m x s 5 V 3 R 6 h O m K H w H m u 1 q 3 1 3 j t q / R v k F H G h 0 t K iLimitations PnP-DM can be further improved in the following two aspects. First, PnP-DM currently requires evaluating the likelihood and prior steps for the entire image at a time. This potentially poses computational challenges in solving large-scale inverse problems (e.g. 3D imaging) or those with expensive likelihood evaluation (e.g. PDE inverse problems). Second, the current theoretical analysis does not consider the approximation error introduced in the likelihood step for general nonlinear inverse problems when running Langevin MCMC for finite iterations. Explicit incorporation of this error would offer further insights into the empirical performance of PnP-DM.Broader impacts We expect this work to make a positive impact in computational imaging and related application domains. For many imaging problems, there is a need to facilitate image reconstruction with expressive image priors and quantify uncertainty, which could lead to better imaging systems that enables further understanding of the imaging target. Nonetheless, as we are introducing DMs as priors into the imaging process, it is inevitable to inherent the potential bias of these models. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our work does not involve models that have a high risk for misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: We have properly credited the owners of all the assets used in the paper. Their licenses are listed in the footnotes and Appendix G.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [Yes]Justification: We have documented the code we provide in the supplement.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model 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. 
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 as details about compensation (if any)?Answer: [NA] Justification: Our research does not involve crowdsourcing.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
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: Our research does not involve human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor 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 NeurIPS 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. 
 
 Appendix: In this section, we formally introduce the interpolation of PnP-DM. We consider the case where the coupling strength ρ is constant, i.e. ρ k ≡ ρ and make the following assumption. Assumption A.1. There exists a unique t * such that σ(t * ) = ρ.This assumption is satisfied for common diffusion models. Popular choices of the noise level schedule include σ(t) = t or σ(t) = √ t, which are monotonically increasing functions of t. We first present two propositions showing that the two steps in SGS can be implemented by running two SDEs. Proposition A.2 (Brownian bridge for the likelihood step). For iteration k with iterate x (k) , the likelihood step of SGS is equivalent to solving the following SDE from t = 0 to t = 1:EQUATION)where x 0 = x (k) and ϕ t (x) := exp[-f (z; y) -1 2ρ 2 (1-t) ∥x -z∥ 2 2 ]dz.Proof. This proposition is due to the Brownian bridge construction presented in Lemma 4 of [78] . This SDE satisfies that p(x1 |x 0 ) ∝ exp -f (x 1 ; y) -1 2ρ 2 ∥x 0 -x 1 ∥ 2 2 .Therefore, solving (9) from t = 0 to t = 1 is equivalent to taking a likelihood step.Proposition A.3 (EDM reverse diffusion for the prior step). For iteration k with iterate z (k) , the prior step of SGS is equivalent to solving the following SDE from t = t * to t = 0:EQUATIONwherex t * = s(t * )z (k) , u(t) := ṡ(t) s(t) , v(t) := s(t) 2 σ(t)σ(t), and p t is the distribution of s(t)x + s(t)σ(t)ϵ with x following the prior distribution p(x) ∝ exp(-g(x)) and ϵ ∼ N (0, I).Proof. First note that ( 10) is exactly (7) written in terms of u(t) and v(t). We know that the (10) is the reverse SDE of the following SDE dx t = u(t)x t dt + v(t)dw t .(11) where x 0 ∼ p(x) and p t is the marginal distribution of x t . As we showed in the main text, it holds for (11) thatp(x 0 |x t ) ∝ exp -g(x 0 ) - 1 2σ(t) 2 ∥x 0 -x t /s(t)∥ 2 2 .As (10) is the time-reversed process of (11) , they share the same path distribution and thus the same conditional distribution p(x 0 |x t ). So, if we set x t * = s(t * )z (k) , we have thatp(x 0 |x t * ) ∝ exp -g(x 0 ) - 1 2σ(t * ) 2 ∥x 0 -z (k) ∥ 2 2 ∝ exp -g(x 0 ) - 1 2ρ 2 ∥x 0 -z (k) ∥ 2 2 ,which is the desired conditional distribution of the prior step. Therefore, solving (10) from t = t * to t = 0 is equivalent to taking a prior step.Due to Proposition A.2 and Proposition A.3, the SDEs ( 9) and ( 10) implement the two desired conditional distributions in SGS. In PnP-DM, the prior step involves a network that approximates the score function of the prior distribution, i.e. s t ≈ ∇ log p t , so the continuous-time process for the actual update isEQUATION)We can then interpolate PnP-DM by considering a dynamic that alternates between running ( 9) and (12) .Since each likelihood step takes 1 unit of time and each prior step takes t * unit of time, the total time of the interpolating process for K iterations of PnP-DM is T := K(t * + 1). We use τ to denote the time that has elapsed from initializing PnP-DM with x (0) . We define {ν τ } and {π τ } as the distributions at time τ of the non-stationary process initialized at x (0) ∼ ν X 0 (Figure 3 top) and the stationary process initialized at x (0) ∼ π X (Figure 3 bottom), respectively. Therefore, we have• ν τ = ν X k , π τ = π X for τ = k(t * + 1) with k = 0, • • • , K, and • ν τ = ν Z k , π τ = π Z for τ = k(t * + 1) + 1 with k = 0, • • • , K -1.A.2 Proof of Theorem 3.1Before proving our main result, we present a key lemma for our analysis, which quantifies the time-derivative of the KL divergence in terms of the Fisher information along a pair of general diffusion processes. Lemma A.4. Given the following pair of diffusion processes dx t = b(x t , t)dt + c(t)dw t (13)EQUATION)where b(•,•) : R n × R → R n , b(•, •) : R n × R → R n, and c(•) : R → R. Let µ t be the distribution of x t initialized with x 0 ∼ µ 0 for (13), and let µ t be the distribution of x t initialized with x 0 ∼ µ 0 for (14) . Then we haveEQUATION)Proof of Lemma A.4. Writing b(•, t) as b t and b(•, t) as b t , by the Fokker-Planck equations of ( 13) and ( 14), we have that∂ t µ t = div c(t) 2 2 ∇ log µ t -b t µ t and ∂ t µ t = div c(t) 2 2 ∇ log µ t -b t µ t .Now we are ready to prove Theorem 3.1.Proof of Theorem 3.1. We first consider the likelihood steps over K iterations of PnP-DM. Applying Lemma 2 of [78] to the likelihood steps ( 9) of the non-stationary and stationary processes, we have that∂ τ KL(π τ ||ν τ ) = - ρ 2 2 FI(π τ ||ν τ ) ≤ - ρ 2 4 FI(π τ ||ν τ ),for τ ∈ [k(t * + 1), k(t * + 1) + 1] with k = 0, ..., K -1. Integrating both sides over τ ∈ [k(t * + 1), k(t * + 1) + 1], we getEQUATIONfor k = 0, ..., K -1.Then, applying Lemma A.4 to the prior steps (12) withb(x t , t) := u(t)x t -v(t) 2 ∇ log p t (x t ) b(x t , t) := u(t)x t -v(t) 2 s t (x t ) c(t) := v(t) δ := inf t∈[0,t * ] v(t),we have that∂ τ KL(π τ ||ν τ ) ≤ - v(τ ) 2 4 FI (π τ ||ν τ ) + 1 v(τ ) 2 v(τ ) 2 (s τ -∇ log p τ ) 2 2 π τ ≤ - v(τ ) 2 4 FI (π τ ||ν τ ) + v(τ ) 2 ∥s τ -∇ log p τ ∥ 2 2 π τ ≤ - δ 2 4 FI (π τ ||ν τ ) + v(τ ) 2 E πτ ∥s τ -∇ log p τ ∥ 2 2 ,for τ ∈ [k(t * + 1) + 1, (k + 1)(t * + 1)] with k = 0, ..., K -1. Integrating both sides over τ ∈ [k(t * + 1) + 1, (k + 1)(t * + 1)], we getEQUATIONwhereϵ score := (k+1)(t * +1) k(t * +1)+1 v(τ ) 2 E πτ ∥s τ -∇ log p τ ∥ 2 2 dτ = t * +1 1 v(τ ) 2 E πτ ∥s τ -∇ log p τ ∥ 2 2 dτ.Finally, combining ( 16) and ( 17) for k = 0, ..., K -1, we obtainT 0 FI (π τ ||ν τ ) dτ ≤ 4[KL(π X ||ν X 0 ) -KL(π X ||ν X K )] min(ρ, δ) 2 + 4Kϵ score δ 2 ≤ 4KL(π X ||ν X 0 ) min(ρ, δ) 2 + 4Kϵ score δ 2 .The proof is concluded by dividing T = K(t * + 1) on both sides.To facilitate the discussion, we first present the following proposition. Proposition A.5. Define a weighting function λ(τ ) over τ ∈ [0, T ] such that for k = 0, ..., K -1,λ(τ ) = ρ 2 if τ ∈ [k(t * + 1), k(t * + 1) + 1], v(τ ) 2 if τ ∈ [k(t * + 1) + 1, (k + 1)(t * + 1)].Then, under the same settings of Theorem 3.1, we haveEQUATION)where ϵ score :=t * +1 1 v(τ ) 2 E πτ ∥s τ -∇ log p τ ∥ 2 2 dτ .Proof. With the definition of λ(τ ), we can apply Lemma 2 of [78] to the likelihood steps and obtainEQUATIONfor k = 0, ..., K -1. Similarly, we can apply Lemma A.4 to the prior steps and obtainEQUATION)whereϵ score := (k+1)(t * +1) k(t * +1)+1 v(τ ) 2 E πτ ∥s τ -∇ log p τ ∥ 2 2 dτ = t * +1 1 v(τ ) 2 E πτ ∥s τ -∇ log p τ ∥ 2 2 dτ.Together, for τ ∈ [0, T ]. We can then get (18) by combining ( 19) and ( 20) for k = 0, ..., K -1 and dividing by T := K(t * + 1).Unlike Theorem 3.1, this proposition calculates the weighted average of the Fisher information along the two processes with the weighting function λ(τ ). The bound in Theorem 3.1 on the unweighted average of Fisher information can be obtained by further lower-bounding the left hand side of ( 18) using the infimum of λ(τ ) over τ ∈ [0, T ]. Given this observation, we can see the role of δ in Theorem 3.1. With a strictly positive δ, the weighting function λ(τ ) is always strictly positive, so the (unweighted) average Fisher information must converge to 0. This is precisely the case for the VP-and VE-SDE [62] . On the other hand, if δ = 0, the Fisher information FI (π τ ||ν τ ) may be increasingly large as λ(τ ) gets closer to 0. For iDDPM and EDM, this could happen near t = 0 in the reverse diffusion at v(0) = 0 for these diffusion processes. Nevertheless, we can instead consider a slightly adjusted diffusion coefficient ṽ(t) := v(t) + ϵ with ϵ > 0. Using the relation between scores and diffusions div(p∇ log p) = ∆p, we get the following reverse SDE which has the same law as (7) at each t:dx t = ṡ(t) s(t) x t + ϵ 2 2 -2s(t) 2 σ(t)σ(t) ∇ log p x t s(t) ; σ(t) dt+ s(t) 2 σ(t)σ(t) + ϵ d wt .In this case, ṽ(t) = s(t) 2 σ(t)σ(t) + ϵ is strictly positive, so the convergence on the unweighted average Fisher information is also guaranteed.Data usage We list the data we have used for our experiments:• For the synthetic prior experiment, we took images from the CelebA dataset [44] , turned them into grayscale, rescaled them to [-1, 1], and resized them to 32 × 32 pixels for efficient computation. We then found the empirical mean and covariance of the images to construct the Gaussian image prior. The test image was randomly drawn from this Gaussian prior. • For the benchmark experiments, we used the first 100 images (index 00000 to 00099) in the FFHQ dataset [37] . For all linear inverse problems, the test images were in RGB and normalized to range [-1, 1]. For all nonlinear problems, the test images were in grayscale and normalized to range [0, 1]. • For the black hole experiments, we used the simulated data used in [63] and the publicly available EHT 2017 data1 that was used to produce the first image of the M87 black hole.The forward model is defined asy ∼ N (Bx, σ 2 y I)where B ∈ R n×n is a circulant matrix that effectively implements a convolution with kernel k under the circular boundary condition. For the Gaussian deblurring problem, we fixed the kernel k as a Gaussian kernel with standard deviation 3.0 and size 61 × 61. For the motion deblurring problem, we randomly generated the kernel k for each test image using the code2 with intensity of 0.5 and size 61 × 61. For fair comparison, the blur kernel for each test image was set the same for all compared methods.The forward model is defined asy ∼ N (P f x, σ 2 y I)where P f ∈ R n f ×n is a matrix that implements a block averaging filter to downscale the images by a factor of f . Specifically, we set f = 4 and used the SVD implementation from the code3 of [38] .Coded diffraction patterns (CDP) CDP is a measurement model originally proposed in [10] . The target x is illuminated by a coherent source and modulated by a phase mask D. The light field then undergoes the far-field Fraunhofer diffraction and is measured by a standard camera. Mathematically, the forward model of CDP is defined as y ∼ N (|F Dx|, σ 2 y I) where F denotes the 2D Fourier transform. We followed [74] to set D as a diagonal matrix with entries drawn randomly from the complex unit circle.Fourier phase retrieval We adopted a similar setting as [18] . In particular, the forward model is defined as y ∼ N (|F P x|, σ 2 y I), where P denotes the oversampling matrix that effectively pads x in 2D matrix form with zeros. We considered a 4× oversampling ratio for grayscale images of size 256 × 256, so P x has a size of 512 × 512.Black hole imaging We adopted the same BHI setup as in [63, 65] . The relationship between the black hole image and each interferometric measurement, or so-called visibility, is given byEQUATIONwhere a and b denote a pair of telescopes, t represents the time of measurement acquisition, i is the imaginary unit, and F t a,b (x) is the Fourier component of the image x corresponding to the baseline between telescopes a and b at time t. In practice, there are three main sources of noise in (21) : gain error g a and g b at the telescopes, phase error ϕ t a and ϕ t b , and baseline-based additive white Gaussian noise η a,b . The gain and phase errors stem from atmospheric turbulence and instrument miscalibration and often cannot be ignored. To correct for these two errors, multiple noisy visibilities can be combined into data products that are invariant to these errors, which are called closure phase and log closure amplitude measurements [12] y cph t,(a,b,c) = ∠(V a,b V b,c V a,c ) := A cph t,(a,b,c) (x), y logcamp t,(a,b,c,d) = log |V t a,b ||V t c,d | |V a,c ||V t b,d | := A logcamp t,(a,b,c,d) (x),where ∠ computes the angle of a complex number. Given a total of M telescopes, there are in total (M -1)(M -2)2closure phase and M (M -3) 2 log closure amplitude measurements at time t, after eliminating repetitive measurements. In our experiments, we used a 9-telescope array (M = 9) from the Event Horizon Telescope (EHT) and constructed the data likelihood term based on these nonlinear closure quantities. Additionally, because the closure quantities do not constrain the total flux (i.e. summation of the pixel values) of the underlying black hole image, we added a constraint on the total flux in the likelihood term. The overall potential function of the likelihood is given byEQUATIONIn this equation, y flux is the total flux of the underlying black hole, which can be accurately measured. We use y := (y cph , y logcamp , y flux ) to denote all the measurements and c, d as the indices for the closure phase and log closure amplitude measurements. Parameters σ cph , σ logcamp were given by the telescope system and σ flux was set to √ 2 in our experiments to constrain the total flux. The data mismatch metric reported in Figure 8 is defined as the sum of the reduced χ 2 values for the closure phase and log closure amplitude measurements, which are calculated using the ehtim.obsdata.Obsdata.chisq function of the ehtim package 4 . Both χ 2 values should ideally be around 1 for data with high signal-to-noise ratio (SNR). Therefore, a data mismatch value around 2 to 3 is considered as fitting the measurements well.Linear forward model and Gaussian noise As we showed in the main paper, in case of linear forward models and Gaussian noise, the likelihood step isπ Z|X=x = N (m(x), Λ -1 )whereΛ := A T Σ -1 A+ 1 ρ 2 I and m(x) := Λ -1 (A T Σ -1 y + 1 ρ 2 x).The bottleneck here is that both the mean and the covariance involve the matrix inverse Λ -1 , which can be prohibitive to compute directly for high-dimensional problems. Nevertheless, the computational cost can be significantly alleviated when the noise is i.i.d. Gaussian, i.e. Σ = σ 2 y I, and A can be efficiently decomposed. For example, if one can efficiently calculate the SVD of the forward model A, i.e. A = U SV T , one can find the Cholesky decomposition of Λ -1 asΛ -1 = LL T where L := V 1 σ 2 S 2 + 1 ρ 2 I -1/2 .Since S is a diagonal matrix, the second term can be calculated with only O(n) complexity. Then, leveraging the property of multivariate Gaussian distribution, we can sample η ∼ N (0, I) and calculate z = m(x)+Lη as a sample that exactly follows the target Gaussian distribution N (m(x), Λ -1 ).An analogous derivation with Fourier transform can be done when A is a circulant convolution matrix.Nonlinear forward model We provide the pseudocode of the LMC algorithm for sampling the likelihood step with general differentiable forward models in Algorithm 2.Algorithm 2 Langevin Monte Carlo for the likelihood step under general A Input: state x, coupling strength ρ > 0, likelihood potential f ( • ; y) with measurements y Hyperparameter: step size γ > 0, number of iterations J > 01: u 0 ← x 2: for j = 0, • • • , J -1 do 3:u j+1 ← u j -γ∇f (u j ; y) -γ ρ 2 (u jx) + √ 2γϵ j where ϵ j ∼ N (0, I) 4: end for 5: return u J Table 3 summarizes the hyperparameters we used for solving the nonlinear inverse problems considered in this work. The EDM framework We formally introduce the EDM formulation [36] using our notations. The forward diffusion process is defined as the following linear Itô SDEEQUATIONwhere u(t) : R → R, v(t) : R → R are the drift and diffusion coefficients. The generative process is the time-reversed version of (23) . According to [2] , it is another Itô SDE of the formEQUATIONwhere p t (x t ) is the marginal distribution of x t . There also exists a reverse probability flow ODEEQUATIONwhich shares the same marginal distributions as (24) . Based on (23), we havep(x t |x 0 ) = N (s(t)x 0 , s(t) 2 σ(t) 2 I)where s(t) := exp t 0 u(ξ)dξ and σ(t) := t 0 v(ξ) 2s(ξ) 2 dξ. We also have x t /s(t) ∼ p(x; σ(t)) where p(x; σ(t)) is the distribution obtained by adding i.i.d. Gaussian noise of standard deviation σ(t)) to the prior data. The idea of the EDM formulation is to write the reverse diffuison directly in terms of the scaling and noise level of x t with respect to x 0 , which are more important than the drift and diffusion coefficients. With the relations between u(t), v(t), p t and s(t), σ(t), p(•; σ(t)), we can rewrite ( 24) and ( 25) asEQUATIONandEQUATIONNote that ( 26) is precisely the SDE (7) we considered for the prior step. Finally, due to the Tweedie's formula [24] , we can approximate∇ xt log p ( • ; σ(t)) by a denoiser [D θ ( • ; σ(t)) -• ]/σ(t) 2 trainedto minimize the ℓ 2 error of a denoising objective. Substituting the score function by the approximation with the denoiser and using the chain rule, we can further rewrite ( 26) and ( 27) asEQUATIONandEQUATION)Pseudocode We provide the pseudocode for our prior step in Algorithm 3. Note that the update rule is precisely the Euler discretization of ( 28) and (29) . The discretization time steps {t i } N i=0 , scaling schedule s(•), and noise schedule σ(•), are kept the same as in Table 1 of [36] . For all experiments, we set the total number of time steps to 100, i.e. N = 100. We note that this does not imply that each prior step has a number of function evaluations (NFE) equal to 100. Since ρ is to a small value as the algorithm runs, the number of steps in later iterations of the algorithm are much fewer than 100. The prior step is similar to the image synthesis process in SDEdit [49] that starts from the middle of the reverse diffusion process. We used the pf-ODE solver for the CDP problem and the SDE solver for all other problems. Part of our code implementation is based on the repository 5 .Algorithm 3 EDM for the prior step (Bayesian denoising with noise level ρ) Input: noisy image z ∈ R n , assumed noise level ρ > 0, pretrained model D θ ( • ; • ) that approximates ∇ log p (x; σ) with (D θ (x; σ)x)/σ2 Hyperparameter: discretization time steps {t i } N i=0 (monotonically decreasing to t N = 0), scaling schedule s(•), noise schedule σ(•), solver (SDE or pf-ODE)1: i * ← smallest i such that σ(t i ) ≤ ρ ▷ Find the starting point of the reverse diffusion 2: v i * ← s(t i * )z ▷ Initialize at time t i * 3: for i = i * , • • • , N -1 do 4: λ ← 2 if solver is SDE else 1 5: d i ← λ σ(ti) σ(ti) + ṡ(ti) s(ti) v i -λ σ(ti)s(ti) σ(ti) D θ vi s(ti) ; σ (t i ) 6: v i+1 ← v i + (t i+1 -t i )d i ▷ Drift 7:if i ̸ = N -1 and solver is SDE then 8:v i+1 ← v i+1 + s(t) 2 σ(t)σ(t)(t i -t i+1 )ϵ i where ϵ i ∼ N (0, I) ▷ Diffusion 9:end if 10: end for 11: return v N Model checkpoint For experiments with FFHQ color images, we used the pretrained checkpoint from [16] available at the repository 6 . For experiments with synthetic data, FFHQ grayscale images, and black hole images, we trained our own models using the same repository. The model network is based on the U-Net architecture in [51] with BigGAN [9] residual blocks, multi-resolution attention, and multi-head attention with fixed channels per head. See the appendix of [16] for architecture details. Specifically, we changed the input and output channels to 1 and 2, respectively, to accommodate grayscale inputs, and reduced the number of down-pooling and up-pooling levels in the U-Net for smaller images. We trained all models until convergence using an exponential moving average (EMA) rate of 0.9999, 32-bit precision, and the AdamW optimizer [45] . Here is a list of training data we used for each model:• For the Gaussian prior model, we randomly generated images from the constructed Gaussian prior distribution. • For the FFHQ grayscale model, we used the images with index 01000 to 69999 in the FFHQ dataset. • For the black hole model, we used 3068 simulated black hole images from the GRMHD simulation, which stands for general relativistic magnetohydrodynamic simulation [22] . See Figure 9 for some example training images. We applied data augmentation with random flipping and resizing, so that the flux spin rotation and the ring diameter vary from sample to sample. Preconditioning Since the checkpoints we used are all trained based on the DDPM (or VP-SDE) formulation [33] , we converted them to the denoiser D θ under the EDM formulation via the VP preconditioning [36] . Specifically, if we denote the pretrained model as F θ ( • ; • ), the model we used for Algorithm 3 isEQUATIONwherec skip (σ) = 1, c out (σ) = -σ, c in (σ) = 1/ √ σ 2 + 1, and c noise (σ) = 999σ -1 VP (σ).Here σ -1 VP (•) is the inverse of the VP-SDE noise schedule defined as σ VP (t) := e Connection to DDS-DDPM in [77] A concurrent work [77] introduced a rigorous implementation of the prior step, called DDS-DDPM, by converting the DDPM [33] (or VP-SDE [62] ) sampler into a reverse diffusion based on the VE-SDE [62] . The diffusion process after the conversion can be used to solve (6) rigorously by properly choosing the starting point. In fact, our formulation admits DDS-DDPM as a special case with the VP preconditioning and reverse diffusion based on the VE-SDE.Here we explicitly show this connection. For the VE-SDE, we have s VE (t) = 1, σ VE (t) = √ t, u VE (t) = 0, and v VE (t) = 1. So (28) becomesEQUATIONApplying the VP preconditioning (30) to (31) , we obtainEQUATION)We can then rescale the time range from [0, 1] to [0, 1000], discretize (32) backward in time over the time steps {τ t } from [77] , and apply the exponential integrator [82] to the drift term, resulting in the following update rule:x t-1 = x t -2( √ τ t - √ τ t-1 )F θ x t √ τ t + 1 ; σ -1 VP ( √ τ t ) + τ t -τ t-1 ϵ where ϵ ∼ N (0, I)Based on the definition τ t := ᾱ-1 t -1 = σ VP (t) 2 in DDS-DDPM, we getEQUATIONThis is exactly the update rule of DDS-DDPM with F θ ( • ; t) denoting the noise estimate ϵ t (•) of DDPM. One can also verify that the initialization in DDS-DDPM is equivalent to ours by checking that ᾱt ≥ 1 η 2 +1 is equivalent to τ t = σ VP (t) 2 ≤ η 2 where η ≡ ρ is the assumed noise level in (6) . As one can see, DDS-DDPM is equivalent to our prior step by choosing the VP-preconditioning, VE reverse diffusion, and a particular integration scheme. In fact, our prior step allows for more general definitions of diffusion processes and includes both the ODE and SDE solvers.Annealing schedule for ρ In this work, we considered an exponential annealing schedule for the coupling strength ρ. We note that the schedule can be more general than exponential decay and we leave the investigation of other decay schedules in future work. Specifically, we specified a starting level ρ 0 , decay rate α, and a minimum value ρ min . Then we set 4 summarizes the annealing hyperparameters that we used for all the inverse problems considered in this work. Initialization For the linear inverse problems, we used the zero initialization, i.e. x (0) = 0 ∈ R n . For the CDP and Fourier phase retrieval problems, we used the Gaussian initialization, i.e.ρ k = max(α k ρ 0 , ρ min ) for k = 0, • • • , K -1. Tablex (0) ∼ N (0, I). For black hole imaging experiments, we used the uniform random initialization between 0 and 1 for each pixel. We found that PnP-DM, as an MCMC algorithm, is insensitive to the initialization. Except for the black hole experiments where we found the negative values would cause problems, any reasonable initialization would lead to comparable results. This observation corroborates our convergence result, which holds for any initialization ν X 0 .We ran 500 iterations for the synthetic prior experiments, 200 iterations for the black hole experiments, and 100 for all other experiments. The numbers were chosen so that the algorithm was fully converged.To collect multiple samples using our method, there are two main approaches:(1) Run a single Markov chain and collect samples after a certain number of iterations, known as the burn-in period, to ensure the chain has converged. (2) Run several independent Markov chains and collect one sample from each chain after convergence. The first approach is more efficient, but the collected samples are not entirely independent and thus may have a small effective sample size. The second approach ensures all samples are fully independent but takes longer to run. In our experiments, we used the first approach for all tests involving 256 × 256 images to enhance efficiency. Specifically, we set the burn-in period to 40 iterations and collected 20 random samples from the remaining 60 iterations (one every 3 iterations). For other experiments, due to the smaller image sizes, we employed the second approach to obtain fully independent samples.Compute All experiments were performed on NVIDIA RTX A6000 and A100 GPUs. The runtime per image depends on several factors, such as the choice of GPU, the total number of iterations and the coupling strength schedule {ρ k } (as it takes more network evaluations for larger ρ for our EDM-based denoiser). In our actual experiments, we ran each image for at least 100 iterations to ensure convergence, which took around 1 minute for a single Markov chain. Here we present a comparison of computational efficiency with the major baselines on a linear super-resolution and a nonlinear coded diffraction patterns problem in Table 5 . The clock time in seconds and number of function evaluations (NFE) are calculated for each method to measure its computational efficiency. All hyperparameters are kept the same for each method as those used for Table 1 and Table 2 in the manuscript. As expected, DM-based approaches (DDRM & DPS) generally yield shorter runtimes due to their lower NFEs. Nevertheless, our PnP-DM method significantly outperforms these methods while achieving comparable runtimes with DPS (≈ 1.5×), despite its larger NFEs (≈ 3×). This is primarily due to two factors: 1) PnP-DM avoids running the full diffusion process by adapting the starting noise level to ρ k at each iteration, and 2) the runtime is further reduced by using an annealing schedule of ρ k . We also note that the runtime reported for DDRM and DPS below is the time it takes to generate one sample. For the linear inverse problem experiments, where we generated 20 samples for each sampling method, PnP-DM was faster than DPS because we took 20 samples that PnP-DM generated along one Markov chain of batch size 1 (hence same runtime as below, around 50 seconds) but DPS requires running a diffusion process with batch size 20, which was significantly slower (around 330 seconds). PnP-ADMM We set the ADMM penalty parameter as 2 and ran for 500 iterations to ensure convergence. We used the pretrained DnCNN denoiser [80] available at the deepinv library7 .DPIR We followed the annealing schedule in [79] and ran for 40 iterations. We used the pretrained DRUNet denoiser [81] available at the deepinv library.DDRM We ran all the experiments with the default parameters: η B = 1.0, η = 0.85, and 20 steps for the DDIM sampler [58] . For the Gaussian deblur problem, we used the SVD-based forward model implementation based on separable 1D convolution. We ran it with an additional batch dimension to collect multiple samples.DPS We followed the original paper to use a 1000-step DDPM sampler backbone. For the linear inverse problems, we used the step size given in [17] , i.e. ζ ′ = 1. For the nonlinear inverse problems, we optimized the step size ζ ′ by performing a grid search, which led to ζ ′ = 3 for CDP and Fourier phase retrieval and ζ ′ = 0.001 for black hole imaging. For the synthetic prior experiments, we also optimized the step size and used ζ ′ = 0.1 for compressed sensing and ζ ′ = 1 for Gaussian deblur. We ran it with an additional batch dimension to collect multiple samples.PnP-SGS We performed a grid search for the coupling parameter ρ and found that ρ = 0.1 worked the best for all problems. We followed the practice in [19] to have a burn-in period of 20 iterations during which the reverse diffusion is early-stopped. We ran the algorithm for 100 iterations in total and collect 20 samples in the 80 iterations after the burn-in period.DPnP We implemented the DDS-DDPM sampler for the prior step. For fair comparison, we used the same annealing schedule for the coupling strength (denoted as η k in [77] ) as PnP-DM. We ran it for the same number of iterations for each inverse problem with the same way of collecting samples as our method.HIO We set β = 0.7 and applied both the non-negative constraint and the finite support constraint.To mitigate the instability of reconstruction depending on initialization, we first repeatedly ran the algorithm with 100 different random initializations and chose the reconstruction that has the best measurement fit. Then we ran another 10,000 iterations with the chosen reconstruction to ensure convergence and report the metrics on the final reconstruction.Image reconstruction with plug-and-play priors Plug-and-Play priors (PnP) [68] is an algorithmic framework that leverages off-the-shelf denoisers for solving imaging inverse problems.Recognizing the equivalence between the proximal operator and finding the maximum a posteriori (MAP) solution to a denoising problem, PnP substitutes the proximal update in many optimization algorithms, such as ADMM [13, 55] and HQS [81, 79] , with generic denoising algorithms, particularly those based on deep learning [48, 81, 79] . The PnP framework enjoys both convergence guarantees [64, 55] and strong empirical performance [75, 1] due to its compatibility with stateof-the-art learning-based denoising priors. Recent works have also proposed learning-based PnP frameworks that have direction interpretations from an optimization perspective [20, 26] . See [35] for a comprehensive review on the theory and practice of PnP.Posterior sampling with MCMC and learning priors Learning-based priors have also been considered in the Bayesian context, where one seeks to sample the posterior distribution defined under a learned prior. An important technique is denoising score matching (DSM) [69] , which connects image denoising with learning the score function of an image distribution. Based on DSM, prior works have incorporated deep denoising priors into MCMC formulations, particularly focusing the Langevin Monte Carlo and its variants as they involve the score function of the target distribution [39, 34, 40, 65, 53] . Recently, methods based on SGS have also gained increasing popularity [52, 19, 6, 27, 77] . Unlike PnP methods based on optimization, these sampling methods possess the ability to generate diverse solutions and quantify the uncertainty of solution space.The remarkable performance of diffusion models [33, 62] on modelling image distributions makes them desirable choices as images priors for solving inverse problems. One popular approach is to leverage a pretrained unconditional model and modify the reverse diffusion process during inference to enforce data consistency [72, 18, 17, 59, 38, 60, 61, 8, 43, 83, 57] . Despite of the promising performance of these methods, they usually involve approximations and empirically driven designs that are hard to justify theoretically and may lead to inconsistent sample distributions. Another line of work learns task-specific models, which achieves higher accuracy at the cost of re-training models for new problems [4, 56, 43] . Methods based on Particle Filtering and Sequential Monte Carlo are also considered to ensure asymptotic consistency [11, 23, 73] . Diffusion models have also been considered as a prior for variational inference [28, 46] and pluy-and-play image reconstruction [32, 47] . In addition to the compressed sensing experiment presented in the main paper, we show another comparison on a Gaussian deblurring problem in Figure 10 . Here, the linear forward model A ∈ R m×n is a 2D convolution matrix with a Gaussian blur kernel of size 7×7 and standard deviation 3.0. Similar to the compressed sensing experiment, both methods yield accurate reconstructions of the mean. However, in terms of the posterior standard deviation, DPS exhibits a notable difference from the ground truth, whereas our method achieves a significantly more accurate result.G L K o Q 4 c B 2 O Q n 5 v G 2 C V T F v C U G 0 E A S p H Q d N 1 N 2 a x J L u D 0 Z E J V / r K Y d G 5 6 U e 0 g Z T F u 7 V 5 M 8 1 s W F 4 3 R 8 8 Y s S W h Z n F Q / P Y 8 v p I h B 6 U V Z j N D s R c A j I U + y u v p + a / d H g W V D 7 V H 1 U + s k M r N Z Q o J U M H / z S Q t 3 1 7 W N 6 + 6 W t Y P K M O H g M j K L r L C P M s D u L j j K s h N H 5 / 8 S R + c Y F / N k K Q g Y + G F 7 w O C 6 j l t V t K 5 H m o A V 1 e p C T Y D K s 6 U A Y D g T n g Z L f m D o h 8 2 f o 8 P c 4 B J E 2 E 5 d h B b f M Y k Y b I s D 3 o 4 C t s V T F S B 2 o R E t H 6 e j 0 1 O q d F G 4 m s 1 1 U d 2 X v 7 m A y m b r f D i + m s L D u s 4 4 M U o U J 1 f t c J q n C 1 J r 0 a T 0 H o T V d n k J V I A v i u v 1 L P I Q + R b 5 i J h F D F F k U Y R I W a Q Q z S 1 q e y G 0 K E S 0 s G i B 6 M 6 i O 0 T S I o n o x q I b R K l F b Z f m F u W I T i w 6 Q X RWe provide visual comparisons for the Gaussian deblur and super-resolution problems in Figure 11 . We provide visual comparisons for the CDP reconstruction problem in Figure 15 , where we visualize one sample for each method. As shown by the red zoom-in boxes, PnP-DM can recover fine-grained features such as the hair threads that are missing in the reconstructions by the baselines. Additional reconstruction examples are given in Figure 16 for the CDP reconstruction problem.We then show some additional reconstruction examples with comparison to DPS in Figure 17 . For each method, we visualize the best reconstruction out of four runs for each test image according to the PSNR value. While DPS failed on around half of the test images, our proposed method provided high-fidelity reconstructions on almost all test images. This comparison highlights the better robustness of our method over DPS. well with data mismatch values around 2. These samples exhibit a ring diameter consistent with the official EHT reconstruction in Figure 1 and share a common bright spot location at the lower half of the ring.Sensitivity analysis on the annealing schedule {ρ k } In Figure 20 , we present a sensitivity analysis on the annealing schedule {ρ k }. In particular, we show the PSNR curves of x k with different exponential decay rates α (left) and minimum coupling levels ρ min (right) for one linear (super-resolution) and one nonlinear (coded diffraction patterns) problem. We have the following conclusions based on the results. First, different decay rates lead to different rates of convergence, which corroborates with our theoretical insights that ρ plays the same role as the step size. The final level of PSNR is not sensitive to different decay rates, as all curves converge to the same level. Second, as ρ min decreases, the final PSNR becomes higher. This is as expected because the stationary distribution of the x k , π X should converge to the true target posterior, p(x|y), as ρ decreases. coupling strength is achieved. Despite being run for 100 iterations in total, our method generates good images in around 40 iterations, which is around 30 seconds and 1600 NFEs.We list the licenses of all the assets we used in this paper: The authors should answer "Yes" if the results are accompanied by error bars, confidence 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 are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run 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). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: We have added a "Broarder Impacts" section in the main manuscript. Guidelines:• Data -CelebA [44]: Unknown Mode 2 Mode 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 M Z p K b T t j W Y 7 + N g J z h 0 9 s x d 8 O l E = " > A A A M g n i c j Z b / b t x E E M f d U q C k B 6 T w J / + s E k U q E I 6 7 l I R E A t H k 0 i a R o u h o k 7 Y q L t H a t 7 Z XP i L 1 0 q k 0 m V n r g i F f O t Q x U I e s z P N d c j e J p L R y A n Z G + d i V N j f X D G p u I h P 9 T x h 7 y P q x 9 z j L t W A z p e / t d N 4 w q Q j q c s y + z K l k z t / 8 v P l 1 U F / s L W 9 / X S b D P q b m 8 O d j Q E M d n a 2 N w c 7 Z N g f l J 9 V q / 6 M z x 8 / + s W e C D e N W K z d k C r 1 2 3 C Q 6 P c Z l Z q 7 I c u X 7 F S x h L o X 1 G e Z I 8 S F p o 7 K S f V Z I 4 k U H l P F 0 m n 4 H a w g 5 H p O w C V k q h N J I + W J W L c i n R B M j q B y Q i K q A 6 L m k S P C h a g C m y C I u q I y F j r g s d 9 x j L n L P N C n m d 4 V E R i 1 m Z Z 4 Q p L h 9 x v r h G m 3 3 4 m N u C t F o X t u Y g 0 R v q R J M O 9 K E L V W w z V s 2 C W Q Y N J x c r l u Z v O l S B M C q C x n d 4 N l A u 7 O 8 s r V 4 3 4 q b y h X i K M F b M H 4 w E Z W i q a K 2 e X lZ F N M + O y 0 v H 4 g o p Q l B 6 n m c j H H Y 8 4 J G Q R k m V f b c Z d 2 e B b U P u U X m p d B J J U S w h Q S p Y f / G N M 3 v d N o 1 r r + e V g U o / 4 m A q Z B Z J d V U e F M + 3 1 8 G Q 5 y 0 / O r v D j 8 7 Q L + T R g h M w s E N 5 Y M B 1 5 X d b 0 i o f q R 1 u y V Y l q h 1 U m i w 4 A M O V 8 N h b s A N D O x R / h o b i B r c g / G b p w j f 4 m k n E M D b Y 4 8 0 s M D Z 4 o j z E N j S i 4 f v x 6 O S E K p 1 l t m Y z n Z X 3 + e 8 2 o L w u n Q v H V 1 N 4 W F c R x 4 U S 2 f F 5 M 5 3 m 8 Z x U W t Q h n Q d h W S 4 n g g z w Q 3 G / j k E O I t c g F x E z i C E K D A o Q K Y M U o p l B T S / 4 B v m I 5 g b N E V 0 b d I 1 I G i Q R X R p 0 i S g 2 q O n S 1 K A U 0 b F B x 4 i u D L p C N D V o i k h E z E c Z q z u 0 w c H T j a 2 8 y 1 t 9 I g o K 3 4 F h u 8 Z 3 F 9 3 2 D d p H d G T Q 0 W 2 F J V P 4 z y D Q J 0 z y G W p e t V I p u 2 k i I z S a 5 l 3 T d c t 0 3 W 2 9 h X Q B 1 Z g K x h A 0 5 R N W j E y p 6 0 h F 4 7 p f T 6 F F s V d P c X O v 2 v h V I 0 M b N 1 K 8 a O M X i A / b + B D x S R u f I D 5 o 4 4 N G 5 j Y + Q j x u 4 z H i l 2 3 8 E v F Z G 5 / 9 x w I r e V p v W 6 W 4 T x T T i g i P w F + 2 A + + e X 1 f h o y r W c b I R 5 k f 0 0 q D n i J 6 j A o h O D H q H 6 B 0 + W a C 6 9 S t P n Q 9 I v d 5 i N G o U 7 RQ = " > A A A M g n i c j Z b h b t s 2 E M f V r t u 6 1 N v S 7 e O + E A k C d F v m 2 W m T J s C G N X H a J E A Q e G n S F p 2 6 g J I o i Y g k K i Q V 2 x H 0 E H u a f d 0 e Y 2 + z o y R S U u w M M 2 C L + t 0 d j / z f U Z a T R l T I w e C f e / c / e v D x J 5 8 + / G z p U e / z L 7 5 c f v z V G 8 E y 7 p J z l 0 W M v 3 O w I B F N y L m k M i L v U k 5 w 7 E T k r X M 5 U v a 3 1 4 Q L y p I z O U v J h x g H C f W p i y W g i + X v 7 S z x C H c 4 d k l u X 2 X Y u / O n u F h e H f R 3 t p 9 v D 7 b Q o L + x s z n Y 2 V a D Z 0 8 3 N w d o 2 B + U n 1 W r / o w v H j / 6 x f a Y m 8 U k k W 6 E h f h t O E j l h x x z S d 2 I F E t 2 J k i K 3 U s c k N x h 7 F J i R x S o + q y h l D O f C L V 0 H P 0 A K 4 i o n C F w i Y j o R O J Y + C y R r U g n A p P D M P d Q j G W I x C x 2 W D Q X p b A O g q h r z B M m Q 5 o E H c e E u s Q H f Z r p X R• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).https://eventhorizontelescope.org/blog/public-data-release-event-horizon-telescope-2017-observationshttps://github.com/LeviBorodenko/motionblur (license unknown)https://github.com/bahjat-kawar/ddrm (MIT license)https://github.com/achael/eht-imaging (GPL-3.0 license)https://github.com/NVlabs/edm/tree/main (Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International Public license)βdt 2 +βmint -1 with β d = 19.9 and β min = 0.1. This adaption allows us to make a fair comparison with other DM-based methods using the same pretrained models. One can also incorporate DMs trained with other formulations into PnP-DM by properly setting the preconditioning parameters.6 https://github.com/jychoi118/P2-weighting (MIT license)https://github.com/deepinv/deepinv (BSD-3-Clause license) 


Decision: This paper received 5 reviews with balanced views but overall positive. The discussion period has enabled to reach a consensus for acceptance.

All reviewers recognized the quality of the work :
- "the paper is well-written, provides clear and concise explanations of the methods and contributions"
- "The strength of the paper is mainly its experimental section, which covers a range of both linear and non-linear inverse problems."
- "Unlike some previous works, the proposed method is effective for both linear and non-linear inverse problems."
- "The paper provides sufficient theoretical analysis to support the proposed method."
- "The theoretical insights provided, including the stationarity guarantee in terms of average Fisher information, add depth to the method."
- "Experiments show the efficacy of the proposed method compared to existing methods"
- "the experiments are very good; the authors first consider a toy example in which the posterior has a closed form. This serves as a nice introduction to the experimental section ..."
- "this is an original way of solving the diffusion posterior sampling problem."

The authors/reviewers discussion period was useful for resolving the reviewers concerns. Yet, it is required for the final version of the paper that the authors :
- provide a more detailed comparison with Ceurdoux et. al. and clearly state the similarities/differences,
- provide an extended theoretical analysis with approximate score (at least in supplementary material).

Reviews: review: paper_summary: paper_summary: In this paper, the authors treat the problem of sampling from posterior of image inverse problems. Their formulation is based on the Split Gibbs Sampler (SGS) algorithm, which alternates between sampling from Moreau regularized versions of the prior and of the likelihood. The main contribution of the paper is to use the EDM formulation of diffusion models for sampling the prior term. The EDM formulation was 
proposed to unify various diffusion models in a common general formulation.
strength_weakness: strength_weakness: 1. **Significance of Contribution**: My primary concern is the overall significance of the contribution. When compared to the works of Coeurdoux et al. (2023) and Xu & Chi (2024), the advancements appear to be relatively minor. Specifically, the contribution seems to involve using a more general formulation of diffusion models (EDM) instead of the more specific DDPM or DDIM. Are there additional conceptual differences? Furthermore, given this similarity, what accounts for the significant performance disparity with these methods?

2. **Algorithm Efficiency**: The algorithms seem to require considerable computational time, necessitating the diffusion process to be run at each iteration. While this is mathematically sensible, it appears excessive. Given the annealing process chosen, how many calls to the denoiser are typically made in practice? Additionally, how does this number compare to methods that perform a single diffusion using an approximation of p(y|x_k), such as DPS? Including a comparative table in the paper would be valuable.

3. **Theoretical Results**: The theoretical results presented, while interesting, are somewhat limited as they assume a fixed rho parameter, which is not the case in practical applications. Additionally, the algorithm is presumed to run with the true score, which is also not realistic. In addition, although it may be beyond the scope of the current paper, providing theoretical insights into the convergence of the proposed Gibbs Sampler (with approximate score) to the true posterior would be beneficial.
questions: questions: - In Figure 4, the diffusion process is represented in a closed-form. Did you train a denoiser, or did you use the closed-form score instead?

- How did you tune the hyperparameters for your method and for each method in the comparison? Did you use their default hyperparameters? For example, DPIR provides hyperparameters fine-tuned for various inverse problems. However, the SR x4 task is not one of the inverse problems for which DPIR's authors fine-tuned their hyperparameters, likely making the default hyperparameters sub-optimal for this task.

- Did you employ any regularization parameters to balance the data-fidelity and regularization terms?

- For the backward diffusion process, did you use the same number of iterations for all values of  $\rho_k$ ? Specifically, is the step size of the diffusion process consistent, or is it adjusted according to the value of $\rho_k$ ?
limitations: limitations: No limitations.
paper_summary: paper_summary: This paper proposes a Markov Chain Monte Carlo algorithm for posterior sampling in both linear and non-linear inverse problems. The core of the proposed method is based on a Split Gibbs Sampler that alternates between two steps: one involving the likelihood and the other the prior. Additionally, the paper connects the Bayesian denoising problem with unconditional generation using Diffusion Models. The proposed method is validated on a range of linear and non-linear inverse problems, with an additional real-world application in black hole imaging.
strength_weakness: strength_weakness: There are no apparent weaknesses. However, I am curious about the reconstruction speed comparison (seconds/image) between DPS and the proposed method. It would be practically very attractive for the community to use this algorithm if a rigorous comparison of reconstruction speed is made.
questions: questions: No questions.
limitations: limitations: No limitations other than those mentioned in the paper.
paper_summary: paper_summary: The paper introduces a method to address Bayesian inverse problems in computational imaging. It leverages the generative capabilities of diffusion models (DMs) to sample the posterior distribution over all possible solutions from noisy and sparse measurements. The method combines a Markov chain Monte Carlo (MCMC) algorithm with a general DM formulation to perform rigorous posterior sampling, effectively integrating state-of-the-art DMs as expressive image priors. The method deviates from current methods that rely on approximating the intractable posterior via separating the forward operator from an unconditional prior over the intermediate noisy image. The approach is validated on six inverse problems, demonstrating superior accuracy and posterior estimation compared to existing DM-based methods.
strength_weakness: strength_weakness: - While the method demonstrates superior performance, the computational cost and efficiency are not thoroughly discussed. A comparison of computational resources required compared to other methods would be beneficial.

The impact of various parameters, such as the annealing schedule for the coupling parameter ρ, on the method's performance is not deeply investigated. A sensitivity analysis could provide insights into the method's robustness.
questions: questions: 1. Can the authors compare the computational resources and time required for PnP-DM versus existing DM-based methods?
      
2. How sensitive is the performance of PnP-DM to the choice of the annealing schedule for the coupling parameter ρ? Is there an optimal range or strategy for selecting these parameters?
limitations: limitations: Yes.
paper_summary: paper_summary: This paper proposes a novel posterior sampling algorithm for solving inverse problems with Diffusion models. The proposed algorithm is based on the Split Gibbs sampling scheme and consists in, first introducing an extended distribution $\pi_\rho(x, z)$ that admits the original posterior $\pi$ as $x$ marginal after letting $\rho$ tend to $0$. Then, for a fixed $\rho$ this joint distribution can be sampled via Gibbs sampling; the form of the joint distribution $\pi_\rho(x, z)$ allows sampling sampling both $X$ and $Z$ easily. For $Z$, this can be done exactly in the case of linear inverse problems, or with Langevin Monte Carlo otherwise. For $X$, the particular structure of the joint distribution is exploited and as it turns out, it can be sampled by simply running the learned backward Diffusion. Notably, this is the only step that requires using the learned denoiser. The theoretical properties of the ideal sampler are inverstigated and extensive experiments are considered.
strength_weakness: strength_weakness: - **Computational cost**:  I believe that the proposed method is much slower than existing alternatives. This is not that much of an issue for me (but of course this depends on the reader and practitioner). In my opinion the compute time is not properly addressed; although the authors discuss it in the bottom of page 21, I would like to see a proper runtime comparison with existing methods. I would like to emphasize that  this is only a minor weakness.

- **Originality**: I believe that the proposed method is related to SDEdit [1]; it can be thought of as a variant of SDEdit. Indeed, to see why this is the case, when the present algorithm is applied to inpainting, the first iteration, the likelihood step, will result in a noisy image with the observation already present In it. This image is then denoised using the backward process during the prior step. This is in fact very similar to sdedit which starts with some initial image, say the pseudo-inverse of the observation, noises it to some given noise level and denoises it using the backward Diffusion process. See Algorithm 3, [1]. I believe that the main difference with the algorithm presented in your paper is the use of the decreasing noise schedule, which makes sense theoretically. For a fair treatment I suggest you mention this in your paper. 

[1] Meng, C., He, Y., Song, Y., Song, J., Wu, J., Zhu, J.Y. and Ermon, S., 2021. Sdedit: Guided image synthesis and editing with stochastic differential equations.
questions: questions: see above
limitations: limitations: see above
paper_summary: paper_summary: It is proposed to solve image inverse problems with diffusion priors by sliced Gibbs sampling (SGS). To sample from the product of prior and likelihood distributions, SGS assumes a pair of variables $x,z$ coupled by a Gaussian and alternates conditional sampling steps: on $x$ using the prior density multiplied with the conditional Gaussian and on $z$ using the likelihood multiplied with the conditional Gaussian. When the prior is a diffusion model, the two steps can be implemented by modifying the given pretrained score model and by Langevin dynamics (or exactly in the case of a linear problem), respectively. This procedure is evaluated on a few standard linear and nonlinear IPs and on a black home imaging problems, and the reconstructions found by the proposed method are typically visually and quantitatively better than the ones found by algorithms from prior work.
strength_weakness: strength_weakness: - On fairness and strength of experimental comparisons:
  - The results do not come with confidence intervals, which makes it hard to assess significance.
  - I may have missed/misunderstood this, but are all methods using the same pretrained prior? Appendix D seems to suggest otherwise, which would be a problem (it isn't fair to compare posterior sampling methods with different priors). 
    - Similarly, different methods seem to be assuming different numbers of sampling steps and otherwise incompatible choices. How were hyperparameters, such as the $\rho_k$ annealing schedule, chosen?
  - Computation cost is not sufficiently discussed.
    - The method is run for at least a hundred Gibbs iterations in each experiment, and each one involves diffusion sampling initialized from some intermediate noise level. This makes the number of model evaluations quite large. It would be necessary to compare the number of function calls, as well as wall time, for all the methods.
- Convergence rate:
  - Theorem 3.1 tells us the rate of convergence to the stationary process for a fixed $\rho$, but not how close the stationary marginal is to the true posterior. They are not equal, in general, for positive $\rho$, and the rate of convergence shown has a $\frac1\rho$ factor, which means the convergence guarantees with small $\rho$ are weaker.
  - To understand the method, it would be important to show the dependence on the number of iterations, both by showing examples at different iterations in an illustrative case and by tracking the convergence of some metric at a function of the number of steps.
- Please also see questions below.
questions: questions: - After equation (1), writing $y=A(x)+n$ where $n\in\mathbb{R}^n$ means that $y$ depends deterministically on $x$. In fact $n$ should be a random variable taking values in $\mathbb{R}^n$, which makes $y$ also a random variable.
- I am surprised not to see any reference to the paper whose title is a suffix of this one's, ["Diffusion models as plug-and-play priors", NeurIPS'22], also abbreviated "PnP" in some subsequent work. There a stochastic optimisation is used to find modes of the posterior.
  - Related to this, "PnP-DM" is not a very informative name for the proposed algorithm when so many of the baseline methods already have "PnP" in the name in some other configuration of letters.
limitations: limitations: Yes, but it would be good to compare computation costs, see above.

Comments: title: Response to authors
Response to rebuttal
comment: We appreciate your constructive feedback on our work. We will include a more detailed comparison with Coeurdoux et al. and the theoretical analysis with an approximate score function in the final version. We will also clarify that the current theory is only with fixed $\rho$.
I thank the authors for their careful answer. 
I think that the manuscript would benefit from a more detailed comparison with Ceurdoux et al. and from an extended theoretical analysis with approximate score. Also if should be made clear that the theory is only with fixed rho parameter.   With these clarifications I am ready to recommend acceptance.
Thank you for your response. We will provide additional clarifications on the similarities and differences between our work and prior works in the final version. Please also see [responses W1.1--W1.3](https://openreview.net/forum?id=Xq9HQf7VNV&noteId=O9TzVFCQmB) to Reviewer u5tS for more information on other differences.
Thanks for the clarification and additional experiments. This makes sense to me.
Thank you for sharing the sampling speed results. The shared results make sense to me. I strongly feel that this paper provides novel enough contributions to be considered for acceptance. I have raised my confidence score to 4. Hope this helps AC decide about the paper.
Thank you for your response. 

Reviewer u5tS mentions that this work bears significant similarities with the work of Coeurdoux et al. and I agree that this is the case with the main difference being that they use a constant $\rho$, which is also more or less the main difference between your work and that of [1] that I am mentioning in my initial review. While reading the paper I did not notice these similarities with the work of Coeurdoux et al., this means that paper in its current form does not adequately address the related work and might mislead the reader into thinking that the paper is much more novel than what it actually is. 

Nonetheless, the empirical results are promising and show that the decreasing schedule is indeed a good idea. I maintain my score but i strongly recommend the authors to modify their paper and clearly state the similarities with previous works.
I'm happy to recommend acceptance given the clarifications. Thank you.
Thank you for your response and comments. Below we answer them point by point.

> **W1.1**

Indeed, we have tried to apply PnP-SGS [1] to nonlinear inverse problems by combining our likelihood step with its prior step design, which was how we obtained the PnP-SGS results for Table 2. While it could still handle the coded diffraction pattern problem, it failed on the more challenging Fourier phase retrieval problem. We also tried to include the annealing strategy for PnP-SGS but found that the method diverged with large $\rho$, probably due to its heuristic design of the prior step, and thus did not benefit from annealing. Overall, our method significantly outperforms PnP-SGS by at least 1dB in PSNR for all linear problems and coded diffraction patterns, and 15dB for Fourier phase retrieval. Our experimental results indicate that PnP-SGS struggles with challenging nonlinear inverse problems, such as Fourier phase retrieval. 

> **W1.2**

There are mainly two aspects. Here we explain using the notations in the original DDPM paper. First, to rigorously implement the prior step of SGS, one need to assume the input as an observation of the *unscaled* image $\boldsymbol{x}_0$ with additive white Gaussian noise. However, in DDPM, the mean of state $\boldsymbol{x}\_t$ is not $\boldsymbol{x}\_0$, but $\sqrt{\bar{\alpha}\_t}\boldsymbol{x}\_0$, which is a *down-scaled* version of $\boldsymbol{x}\_0$. So, it is inaccurate to directly use DDPM as a denoiser. This mismatch is particularly significant when starting from a large $t^\ast$ as $\sqrt{\bar{\alpha}\_{t^\ast}}$ is close to 0. Second, according to the SGS formulation, the noise estimation module is unnecessary, and one should always denoise $\boldsymbol{z}^{(k)}$ (the $\boldsymbol{z}$ iterate of SGS at iteration $k$) assuming a noise level $\rho$. However, PnP-SGS does not take into account the hyperparameter $\rho$ in the denoising problem. It is unclear how the generated $\boldsymbol{x}^{(k+1)}$ relates to the target conditional distribution $\pi^{X|Z=\boldsymbol{z}^{(k)}}(\boldsymbol{x})$ for each prior step of SGS. Therefore, the posterior distribution sampled from the denoiser in [1] is not the desired posterior distribution even with a perfect score function. 

> **W1.3**

We never claimed that the idea of annealing is our contribution. Nevertheless, we believe that it is our contribution to propose a general framework that accommodates annealing in an easy yet rigorous way for SGS and leads to significant performance improvement. Our focus is more on the algorithm design and experimental validation. To the best of our knowledge, no existing theory on SGS provides a non-asymptotic convergence bound with annealing, so this is a a challenging open question. Although our current theory assumes fixed $\rho$, it still provides some theoretical insights on the algorithm behavior, such as the interpretation of $\rho$ as step size, and potentially opens up new theoretical directions for SGS based on the Fisher information. We hope to extend the analysis to the annealing case in the future.

> **W2**

We respectfully disagree with your last statement. DPS requires backpropagating through the entire denoiser network for each diffusion step after the forward pass of the denoiser, introducing a significant computational overhead. On the other hand, our method does not require doing so. This is the main reason why the clock time does not scale linearly with NFE for these two methods. Moreover, unlike DPS that applies a likelihood update for every function evaluation. Our method does so only every multiple function evaluations. 

> **W3.1**

As we showed in Appendix C.3 "Annealing schedule for $\rho$" paragraph, our annealing schedule decreases to a fixed minimum level at $\rho_{\min}$, so we indeed fixed $\rho$ at small values after certain numbers of iterations. Furthermore, we did run our algorithms for the synthetic prior experiments with a fixed schedule of $\rho$ throughout the process, as shown in Table 4. Experimental results show that our method can accurately sample the posterior, which corroborates our theory. For experiments on linear inverse problems and coded diffraction patterns with FFHQ images, the results with fixed $\rho$ are on par with those with annealing, probably because the problems are not highly non-convex. For more challenging nonlinear problems, such as the Fourier phase retrieval, we empirically find that an annealing schedule is essential to overcome the non-convexity of the problem and provide accurate reconstruction.
 
> **W3.2**

By combining the above argument with those in the manuscript, we will get the following bound with score function error: $$\frac{1}{T} \int_0^{T} \mathsf{FI}\left(\pi_\tau || \nu_\tau\right) \mathrm{d} \tau \leq \frac{4\mathsf{KL}(\pi^X||\nu_0^X)}{K(1+t^\ast) \min(\rho, \delta)^2} + \frac{1}{K(1+t^\ast)\delta^2}\sum_{k=1}^K\epsilon_{\text{score}}(k).$$
I thank the authors for their detailed answers. Here are additional remarks / questions. 

**W1:**

1. Although not demonstrated in the experiments, I believe the approach in [1] could be equally applied to nonlinear inverse problems.
   
2. Could you clarify why you consider the denoiser design based on diffusion models in [1] to be heuristic?

3. The annealing strategy, despite not being adopted in [1], is commonly used in Plug-and-Play literature. However, incorporating this strategy introduces a discrepancy between your theoretical framework and experimental results, which is concerning.

4. I acknowledge that using EDM as the diffusion model offers a more flexible and optimized framework.

**W2:**

I do not understand why the clock time does not scale linearly with NFE. It appears that the operations outside the denoiser evaluation in both DPS and your method are relatively insignificant in time.

**W3:**

1. I strongly disagree with your statement, as you did not run the algorithms with fixed $\rho$ even for smaller values.

2. Can you explicitly state the final convergence bound analogous to Theorem 3.1, including the time integral from 0 to T, that would result from this additional theoretical consideration?