Paper deep dive

Taming Score-Based Denoisers in ADMM: A Convergent Plug-and-Play Framework

Rajesh Shrestha, Xiao Fu

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 171

Intelligence

Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 97%

Last extracted: 3/13/2026, 1:08:56 AM

Summary

The paper introduces ADMM-PnP with an AC-DC denoiser to address manifold mismatch and convergence issues when integrating score-based generative models into ADMM for solving inverse problems. The AC-DC denoiser uses a three-stage process (auto-correction, directional correction via Langevin dynamics, and score-based denoising) to align iterates with noise-trained manifolds. The authors provide theoretical convergence guarantees for the framework under both constant and adaptive step-size schedules.

Entities (5)

AC-DC denoiser · algorithm · 99%ADMM-PnP · framework · 99%Inverse problems · problem-domain · 98%Tweedie's lemma · mathematical-principle · 97%Score-based generative models · technology · 95%

Relation Signals (3)

ADMM-PnP → utilizes → AC-DC denoiser

confidence 99% · we propose ADMM plug-and-play (ADMM-PnP) with the AC-DC denoiser

AC-DC denoiser → incorporates → Langevin dynamics

confidence 98% · directional correction (DC) using conditional Langevin dynamics

AC-DC denoiser → addresses → manifold mismatch

confidence 95% · To address the manifold mismatch issue, we propose ADMM plug-and-play (ADMM-PnP) with the AC-DC denoiser

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:While score-based generative models have emerged as powerful priors for solving inverse problems, directly integrating them into optimization algorithms such as ADMM remains nontrivial. Two central challenges arise: i) the mismatch between the noisy data manifolds used to train the score functions and the geometry of ADMM iterates, especially due to the influence of dual variables, and ii) the lack of convergence understanding when ADMM is equipped with score-based denoisers. To address the manifold mismatch issue, we propose ADMM plug-and-play (ADMM-PnP) with the AC-DC denoiser, a new framework that embeds a three-stage denoiser into ADMM: (1) auto-correction (AC) via additive Gaussian noise, (2) directional correction (DC) using conditional Langevin dynamics, and (3) score-based denoising. In terms of convergence, we establish two results: first, under proper denoiser parameters, each ADMM iteration is a weakly nonexpansive operator, ensuring high-probability fixed-point $\textit{ball convergence}$ using a constant step size; second, under more relaxed conditions, the AC-DC denoiser is a bounded denoiser, which leads to convergence under an adaptive step size schedule. Experiments on a range of inverse problems demonstrate that our method consistently improves solution quality over a variety of baselines.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

170,274 characters extracted from source content.

Expand or collapse full text

Taming Score-Based Denoisers in ADMM: A Convergent Plug-and-Play Framework Rajesh Shrestha School of EECS Oregon State University shresthr@oregonstate.edu &Xiao Fu School of EECS Oregon State University xiao.fu@oregonstate.edu Abstract While score-based generative models have emerged as powerful priors for solving inverse problems, directly integrating them into optimization algorithms such as ADMM remains nontrivial. Two central challenges arise: i) the mismatch between the noisy data manifolds used to train the score functions and the geometry of ADMM iterates, especially due to the influence of dual variables, and i) the lack of convergence understanding when ADMM is equipped with score-based denoisers. To address the manifold mismatch issue, we propose ADMM plug-and-play (ADMM-PnP) with the AC-DC denoiser, a new framework that embeds a three-stage denoiser into ADMM: (1) auto-correction (AC) via additive Gaussian noise, (2) directional correction (DC) using conditional Langevin dynamics, and (3) score-based denoising. In terms of convergence, we establish two results: first, under proper denoiser parameters, each ADMM iteration is a weakly nonexpansive operator, ensuring high-probability fixed-point ball convergence using a constant step size; second, under more relaxed conditions, the AC-DC denoiser is a bounded denoiser, which leads to convergence under an adaptive step size schedule. Experiments on a range of inverse problems demonstrate that our method consistently improves solution quality over a variety of baselines. 1 Introduction Inverse problems arise in many fields, including medical imaging (Song et al., 2022b; Jin et al., 2017; Arridge, 1999), remote sensing (Entekhabi et al., 1994; Combal et al., 2003), oceanography (Bennett, 1992), and computational physics (Raissi et al., 2019; Tarantola, 2005). Their solutions typically rely on incorporating prior knowledge or structural assumptions about the target signals, either through explicit regularization or data-driven models. Classical approaches to inverse problems often rely on handcrafted regularizers, such as the ℓ1 _1 norm for sparsity (Yang et al., 2010; Elad and Aharon, 2006; Dabov et al., 2007) and the nuclear norm for low-rank structure (Semerci et al., 2014; Hu et al., 2017). Deep learning introduced a new paradigm of using learned generative models—VAEs, GANs, and normalizing flows—as data-driven regularizers (Ulyanov et al., 2020; Alkhouri et al., 2024; Shah and Hegde, 2018), offering more expressive priors by capturing complex distributions. More recently, pre-trained score functions from diffusion models have gained attention for inverse problems (Song et al., 2022b; Chung et al., 2023), as they effectively approximate data distributions and align solutions with the underlying data geometry (Xiao et al., 2022). The use of pre-trained score functions in inverse problems mainly falls into two categories. The first modifies the MCMC process of diffusion sampling to incorporate observation information, as in DPS (Chung et al., 2023) and DDRM (Kawar et al., 2022), where observations guide unconditional score functions to perform posterior sampling. The second integrates score functions into deterministic or stochastic optimization algorithms; for example, Wang et al. (2024) and Song et al. (2023) use them as “projectors” to keep iterates on the desired data manifold. Furthermore, building on Tweedie’s lemma, which links score functions to signal denoising, works such as (Zhu et al., 2023; Mardani et al., 2024; Renaud et al., 2024b) employ score as denoisers in proximal-gradient-like steps. Challenges. The works in (Zhu et al., 2023; Mardani et al., 2024; Li et al., 2024) present flexible “plug-and-play (PnP)” paradigms that integrate diffusion models with optimization algorithms. However, two challenges remain in this line of work. First, score functions are trained on noisy data manifolds constructed via Gaussian perturbations, whereas optimization iterates need not lie on such manifolds, leading to geometry mismatch and degraded denoising performance. Remedies such as stochastic regularization (Renaud et al., 2024b) or purification (Li et al., 2024) add Gaussian noise to the iterates, but this does not guarantee alignment with the score manifolds. Second, the theoretical understanding of these methods—particularly their convergence properties when combining the score denoisers with various optimization paradigms—remains limited. Contributions. We propose to integrate score-based denoisers with the ADMM framework. Using ADMM iterates with score-based denoising is particularly challenging, as the presence of dual variables further distorts the “noise” geometry—likely explaining why score-based denoising has rarely been combined with primal–dual methods. Nevertheless, ADMM remains attractive for its flexibility in handling diverse inverse problems with multiple regularizers. Our contributions are: ▶ Score-Based AC-DC Denoiser: To mitigate the manifold mismatches, we propose a three-stage denoiser consisting of (1) additive Gaussian noise auto-correction (AC), (2) conditional Langevin dynamics-based directional correction (DC), and (3) score-based denoising. The AC stage pulls ADMM iterates toward neighborhoods of noise-trained manifolds, while DC refines alignment without losing signal information. This combination balances efficiency and accuracy, making score-based denoising effective within ADMM. ▶ Convergence Analysis: We show that, under proper AC-DC parameters, each ADMM iteration is weakly nonexpansive, ensuring convergence to a fixed-point neighborhood under constant step sizes under strongly convex losses. We further relax convexity and prove that an adaptive step-size scheme (Chan et al., 2016) guarantees convergence with high probability. These results extend prior ADMM-PnP convergence theory (Ryu et al., 2019; Chan et al., 2016) to score-based settings. Our method is validated on diverse applications—including inpainting, phase retrieval, Gaussian and motion deblurring, super-resolution, and high dynamic range (HDR). Notation. The detailed notation designation is listed in Appendix A.1. 2 Background Inverse Problems. We consider the typical inverse problem setting where =()+ y=A( x)+ ξ (1) where :ℝd→ℝnA:R^d ^n is the measurement operator, n≤dn≤ d, and ξ is additive noise. In some inverse problems, e.g., signal denoising and image deblurring, we have n=dn=d; while for some other problems, e.g., data compression and recovery, we have n<dn<d. The goal is to recover x from the y, with the knowledge of A. Structural regularization on x is often used to underpin the desired solution: minℓ(||())+h(), _ x~ ( y~||~A( x))+h( x), (2) where ℓ(||()) ( y~||~A( x)) is a divergence term that measures the similarity of y and ()A( x) (e.g., ‖−()‖2\| y-A( x)\|^2), and h()h( x) is a structural regularization term (e.g., ‖1\| x\|_1 for sparse x). Diffusion-Based Inverse Problem Solving. Diffusion models can also be used for solving inverse problems, in ways more subtle than direct regularization. Consider training a diffusion model on 0∼pdata x_0 p_ data via denoising score matching (Song et al., 2021), where the forward process is t|0∼(0,σ(t)),t∈[0,T], x_t| x_0 ( x_0,σ(t) I), t∈[0,T], with variance schedule σ(0)=0σ(0)=0 and σ(t)σ(t) increasing in t. After training, the model provides a score function (,σ(t))≈∇tlog⁡p(t), s_ θ( x,σ(t))≈ _ x_t p( x_t), where p(t)p( x_t) is the marginal density of t x_t. This induces noisy data manifolds ℳσ(t)=supp(t),∀t∈[T],M_σ(t)= supp( x_t), ∀ t∈[T], which are continuous since t x_t is generated by Gaussian perturbations of 0∼pdata x_0 p_ data. These score functions can then be leveraged in different ways to assist inverse problem solving. ▶ Posterior Sampling: Many works formulate inverse problems as posterior sampling from p(|)∝p()p(|)p( x| y) p( x)p( y| x). These methods approximate ∇tlog⁡p(|t) _ x_t p( y| x_t) and combine it with the learned score (t,σ(t))≈∇log⁡p(t) s_ θ( x_t,σ(t))≈∇ p( x_t) to perform stochastic sampling as t decreases (Chung et al., 2023; Song et al., 2022a; Kawar et al., 2021; 2022; Wang et al., 2023). While effective, their performance is often limited by the accuracy of the approximation to ∇log⁡p(|t)∇ p( y| x_t). ▶ Plug-and-Play (PnP) Approaches: Instead of sampling schemes, another line of work employs deterministic or stochastic optimization to solve (2), plugging score functions into the updates as structural regularizers. A representative example is DiffPIR (Zhu et al., 2023), which adopts a variable-splitting reformulation, min,⁡ℓ(∥())+μ2‖−‖22+h(), _ x, z\; ( y\|A( x))+ μ2\| x- z\|_2^2+h( z), where the z-subproblem at iteration k reduces to a standard denoising step: Denoising: (k+1) : z^(k+1) =arg⁡min⁡μ2‖(k+1)−‖22+h(). = _ z μ2 \| x^(k+1)- z \|_2^2+h( z). (3) This step is then tackled using a score-based denoiser: (k+1)←Dσ(k)(~(k+1))=~(k+1)+(σ(k))2(~(k+1),σ(k)) z^(k+1)← D_σ^(k)( x^(k+1))= x^(k+1)+(σ^(k))^2 s_ θ( x^(k+1),σ^(k)) (4) where ~(k+1)=(k+1)+σ(k)(ζ(~(k)−(k+1))/σ(k)+(1−ζ)) x^(k+1)= x^(k+1)+σ^(k)(ζ ( x^(k)- x^(k+1))σ^(k)+(1-ζ) n) with ∼(,) n ( 0, I) and ζ∈[0,1]ζ∈[0,1]; also see (Li et al., 2024) for a similar method. These denoisers are designed following the Tweedie’s lemma (Robbins, 1992) (see Appendix A). The construction of ~(k) x^(k) is meant to make the inputs to the score function closer to a certain ℳσ(k) M_σ^(k). Using score-based denoising, h()h( z) is implicitly reflected in the denoising process and thus does not need to be specified analytically. Another line of approaches explicitly construct regularizers h()h( x) whose gradients correspond to applying the score function. Examples include RED-diff (Mardani et al., 2024) and SNORE (Renaud et al., 2024a). In SNORE, the regularizer is defined as h()=~|log⁡pσ(~)h( x)=E_ x| x p_σ( x), ~←+σϵ x← x+σ ε, ϵ∼(,) ε ( 0, I) and taking its gradient in (2) yields: SNORE update: (k+1) SNORE update: x^(k+1) ←(k)−δ∇ℓ(||((k)))−η(~(k)−Dσi(~(k))), ← x^(k)-δ∇ ( y|| A( x^(k)))-η( x^(k)-D_ _i( x^(k))), (5) where, again, Dσ(~)=~+σ2θ(~,σ)D_σ( x)= x+σ^2 s_θ( x,σ) by the Tweedie’s lemma. Challenges — Manifold Mismatch and Convergence. Score-based PnP methods face two main challenges. First, the score is not trained on algorithm-induced iterates (e.g., (k+1) x^(k+1) in (3)). While both t x_t and (k) x^(k) can be seen as noisy versions of ∼pdata x p_ data, t x_t follows Gaussian noise whereas the distribution of (k) x^(k) is unclear. Many works attempt to bridge this gap by injecting Gaussian noise before applying the score function (cf. ~(k+1) x^(k+1) in (3), ~(k) x^(k) in (5)), or by purification-based schemes (Nie et al., 2022; Alkhouri et al., 2023; Meng et al., 2022). Yet noise injection alone is insufficient, and overfitting to measurement noise remains an issue (Wang et al., 2024). Second, the understanding to convergence of score-based PnP remains limited. Unlike classical denoisers with established theory, the geometry mismatch above makes it unclear whether iterates stabilize or under what conditions convergence can be ensured. Existing analyses mostly cover primal algorithms (see, e.g., (Renaud et al., 2024a)). Primal–dual methods such as ADMM offer greater flexibility for handling multiple regularizers and constraints, but their convergence with score-based denoisers remains unclear, as dual variables further complicate the manifold geometry of the iterates. 3 Proposed Approach In this section, we propose our score-based denoiser and embed it in the ADMM framework. While we focus on ADMM due to its flexibility, the denoiser can be plugged into any other proximal operator based schemes (e.g. proximal gradient or variable-splitting as in DiffPIR). Preliminaries of ADMM-PnP. ADMM-based inverse problem solvers start by rewriting (2) as min,ℓ(||())+γh() s.t. = _ x, z ( y|| A( x))+γ h( z)~ s.t. ~ x= z (6) The augmented Lagrangian of (6) is given by Lρ(,,;)=ℓ(||())+γh()+T(−)+ρ2∥−∥22L_ρ( x, y, λ; y)= ( y|| A( x))+γ h( z)+ λ^T( x- z)+ ρ2 \| x- z \|_2^2, where λ is the dual variable and ρ>0ρ>0 is the penalty parameter. At iteration k, ADMM updates are as follows: (k+1) x^(k+1) =argmin1ρℓ(||())+12∥−(k)+(k)∥22 = _ x 1ρ ( y|| A( x))+ 12 \| x- z^(k)+ u^(k) \|_2^2 (7a) (k+1) z^(k+1) =arg⁡min⁡γρh()+12‖(k+1)−+(k)‖22 = _ z γρh( z)+ 12 \| x^(k+1)- z+ u^(k) \|_2^2 (7b) (k+1) u^(k+1) =(k)+((k+1)−(k+1)) = u^(k)+( x^(k+1)- z^(k+1)) (7c) where (k)=(k)/ρ u^(k)= λ^(k)/ρ is the scaled dual variable. The subproblem (7b) is a denoising problem, and thus can be replaced by (k+1)=Proxγρh((k+1)+(k))=Dσ(k)(~(k)). z^(k+1)=Prox_ γρh( x^(k+1)+ u^(k))=D_σ^(k)( z^(k)). (8) The above is the classical ADMM-PnP method; see (Chan et al., 2016; Ryu et al., 2019). Same as before, Eq. (8) can be replaced by the score-based denoising following the Tweedie’s lemma. However, as ~(k)=(k+1)+(k) z^(k)= x^(k+1)+ u^(k) could be in any of the manifolds ℳσ(t) M_σ(t) on which the score function was trained, such naive replacement does not ensure effective denoising. The existence of the dual variable (k) u^(k) makes the noise distribution in ~(k) z^(k) even harder to understand. Proposed Approach: The AC-DC Denoiser. To address the manifold mismatch issues, we propose a three-stage denoiser. To be specific, in the kkth iteration of the ADMM algorithm, we use the denoising process shown in Algorithm 1. Note that the Tweedie’s lemma step (line 8) can also be substituted by a score ODE based process (Karras et al., 2022) initialized at dc(k) z^(k)_ dc. Our algorithm using these two different denoisers will be referred to as Ours-tweedie and Ours-ode, respectively. Figure 1: Left: direct denoising of ~(k) z^(k) using score functions could lead to unnatural recovered signals with artifacts. Right: AC-DC denoising brings ~(k) z^(k) closer to ℳσ(k)M_σ^(k), and then uses the score function to bring ~(k) z^(k) to the data manifold ℳdataM_ data. Algorithm 1 AC-DC Denoiser at iteration k of ADMM in (7) 1: auto correction (AC): zac(k)←z~(k)+σ(k)n,n∼(,I) z^(k)_ ac← z^(k)+σ^(k) n,~ n N( 0, I) 2: directional correction (DC): 3: (0)←ac(k) w^(0)\;←\; z^(k)_ ac 4: for j=0j=0 to J−1J-1 do 5: (k,j+1)←(k,j)+η(k)(1/σ(k)2(ac(k)−(k,j))+θ((k,j),σ(k)))+2η(k), w^(k,j+1)← w^(k,j)+η^(k) ( 1 _ s^(k)^2( z^(k)_ ac- w^(k,j))+ s_θ( w^(k,j),σ^(k)) )+ 2η^(k) n, 6: end for 7: dc(k)←(J) z^(k)_ dc\;←\; w^(J) 8: Denoising :tw(k)←[0|t=dc(k)]=dc(k)+(σ(k))2θ(dc(k),σ(k)): z^(k)_ tw [ z_0| z_t= z^(k)_ dc]= z^(k)_ dc+(σ^(k))^2 s_θ( z^(k)_ dc,σ^(k)) (Tweedie Denoising) (Alternative: zode(k)←z0 z^(k)_ ode← z_0 by solving: dztdt=λ(t)sθ(zt,t), d z_tdt=λ(t)s_θ( z_t,t), with zσ(k)=zdc(k) z_σ^(k)= z_dc^(k)) (ODE Denoising) The rationale of the AC-DC denoiser is illustrated in Fig.1. Recall that the score function is most effective on the noisy data manifolds ℳσ(t)t=1T\ M_σ(t)\_t=1^T, as it is trained over them. Since ADMM-induced iterates ~(k) z^(k) need not lie on these manifolds, directly applying score-based denoising may be ineffective. The AC step addresses this by adding Gaussian noise, making ac(k) z^(k)_ ac closer to some ℳσ(t) M_σ(t) (see AppendixB). This idea is related to the “purification” step in (Nie et al., 2022; Alkhouri et al., 2023; Li et al., 2024) and noise-added denoising in (Mardani et al., 2024; Renaud et al., 2024a; Zhu et al., 2023) (cf. (5) and (3)). However, AC alone does not guarantee manifold alignment. The proposed DC step, based on Langevin dynamics, further refines dc(k) z^(k)_ dc toward ℳσ(k) M_σ^(k). To see the idea, let us break down the three steps. First, the AC step gives ac(k)=σ(k)+~(k),(k)=~(k)−♮(k),♮(k)∼pdata z^(k)_ ac= z_σ^(k)+ s^(k), s^(k)= z^(k)- z^(k)_ ,\; z^(k)_ p_ data (9) where ♮(k) z^(k)_ is denoised signal of ~(k) z^(k), σ(k)=♮(k)+σ(k)1 z_σ^(k)= z^(k)_ +σ^(k) n_1, ~(k)=2σ(k)2+(k) s^(k)= 2σ^(k) n_2+ s^(k), and 1,2∼i.i.d(,d) n_1, n_2 i.i.d N( 0, I_d). Given a sufficiently large σ(k)σ^(k), ac(k) z_ ac^(k) would have dominated by Gaussian noise—but not necessarily on any of ℳσt M_ _t where the score was trained. Starting from ac(k) z^(k)_ ac, the DC step runs a few iterations of Langevin dynamics targeting the distribution p(σ(k)|ac(k))p( z_σ^(k)| z^(k)_ ac). This is because supp(σ(k)|ac(k))⊆supp(σ(k))=ℳσ(k) supp( z_σ^(k)| z^(k)_ ac) supp( z_σ^(k))= M_σ^(k). In addition, p(σ(k)|ac(k))p( z_σ^(k)| z^(k)_ ac) at the same time retains the information of ac(k) z^(k)_ ac (thereby the information from the measurements). Assume that the forward process used in training the score has sufficiently small time intervals, ℳσ(k) M_σ^(k) is approximately contained in ℳσtt=1T\ M_ _t\_t=1^T. This way, when applying Tweedie’s lemma for denoising, the step is expected to be effective, as the score was trained over ℳσtt=1T\ M_ _t\_t=1^T. Note that the conditional score for the Langevin dynamics step can be expressed as follows ∇log⁡p(σ(k)|ac(k))=(σ(k),σ(k))+∇log⁡p(ac(k)|σ(k)).∇ p( z_σ^(k)| z^(k)_ ac)= s_ θ( z_σ^(k),σ^(k))+∇ p( z^(k)_ ac| z_σ^(k)). (10) Ideally, one would use the exact ∇log⁡p(ac(k)|σ(k))∇ p( z^(k)_ ac| z_σ^(k)) for the DC step—which is unavailable. In practice, we approximate p(ac(k)|σ(k))p( z^(k)_ ac| z_σ^(k)) using a Gaussian distribution. Note that under proper scheduling of σ(k)σ^(k) and mild regularity conditions on (k) s^(k), e.g., when Var((k))1/2≪σ(k)Var( s^(k)) 12 σ^(k), the likelihood can be well-approximated by a locally quadratic form, leading to ∇log⁡p(ac(k)|σ(k))≈−1/σ(k)2(σ(k)−ac(k))∇ p( z^(k)_ ac| z_σ^(k))≈- 1 _ s^(k)^2( z_σ^(k)- z^(k)_ ac) and the DC step in Algorithm 1. 4 Convergence Analysis 4.1 Convergence of Under Weakly Non-Expansive Residuals Following the established convention in ADMM-PnP, e.g., (Buzzard et al., 2018; Sun et al., 2019; Chan, 2019; Teodoro et al., 2017), we aim at understanding the convergence properties when the AC-DC denoiser is used. We will use the following definitions: Definition 1 (Fixed point convergence). Let T:→T:X be the update map of an iterative algorithm, and let (0)∈ x^(0) be arbitrary initialization. The algorithm is said to converge to a fixed point ∗ x^* if for any δ>0δ>0 there exists Kδ>0K_δ>0 such that the sequence generated by the algorithm (k)k∈ℕ+\ x^(k)\_k ^+ satisfies ‖(k)−∗‖2<δ \| x^(k)- x^* \|_2<δ for all k≥Kδk≥ K_δ. Equivalently, limk→∞(k)=∗ _k→∞ x^(k)= x^* with T(∗)=∗T( x )= x . Definition 2 (Sequence convergence to a δ-ball). For a certain δ>0δ>0, a sequence (k)k∈ℕ+\ x^(k)\_k ^+ is said to converge within a δ-ball if there exists K>0K>0 and ∗ x^* such that the following holds for all k≥Kk≥ K . ‖(k)−∗‖2≤δ \| x^(k)- x^* \|_2≤δ (11) Comparing the two definitions, Definition 2 is a weaker statement; that is, even when k→∞k→∞, (k)→∗ x^(k)→ x does not necessarily happen. Nonetheless, convergence to a δ-ball is still meaningful. The notion of δ-ball convergence is often used in numerical analysis for stability characterization; see, e.g., Ren and Argyros (2021); Ren and Wu (2009); Liang (2007). Definition 3 (ADMM convergence to a δ-ball). ADMM is said to converge within a δ ball if the sequences (k)k∈ℕ+\ x^(k)\_k ^+ and (k)k∈ℕ+\ u^(k)\_k ^+ obtained from ADMM converges within a δ-ball. To proceed, consider the following assumption: Assumption 1. For a certain δ>0δ>0 there exists ϵ≤1ε≤ 1 such that for all ~1,~2∈ℝd z_1, z_2 ^d, the following holds: ‖Rσ(~1)−Rσ(~2)‖22≤ϵ2‖~1−~2‖22+δ2 \|R_σ( z_1)-R_σ( z_2) \|_2^2≤ε^2 \| z_1- z_2 \|_2^2+δ^2 (12) where Rσ(~)=(Dσ(k)−I)(~)R_σ( z)=(D_σ^(k)-I)( z) with I being the identity function (i.e., I()=I( z)= z). Here, the notation (Dσ(k)−I)(~)(D_σ^(k)-I)( z) denotes the residual of Dσ(k)D_σ^(k) i.e. Dσ(k)(~)−~D_σ^(k)( z)- z. The next theorem extends the fixed point convergence of ADMM-PnP in Ryu et al. (2019). Unlike Ryu et al. (2019) where RσR_σ needs to be strictly contractive, our result allows RσR_σ to be weakly contractive: Theorem 1. Under Assumption 1, assume that ℓ is μ-strongly convex. Then, there exists ∗ x^*, ∗ u^* and K>0K>0 such that the sequences (k)k∈ℕ+\ x^(k)\_k ^+ and (k)k∈ℕ+\ u^(k)\_k ^+ generated by ADMM-PnP using a fixed step size ρ satisfies ‖(k)−∗‖2≤r \| u^(k)- u^* \|_2≤ r and ‖(k)−∗‖2≤r \| x^(k)- x^* \|_2≤ r with r=(1+ρ+μ)δ¯/1−ϵ¯2r=(1+ ρ+μ) δ/ 1- ε^2 for all k≥Kk≥ K when ϵ/μ(1+ϵ−2ϵ2)<1/ρ εμ(1+ε-2ε^2)< 1ρ where δ¯2=δ2ϵ¯ϵ δ^2= δ^2 ε and ϵ¯=ρ+ρϵ+μϵ+2μϵ2ρ+μ+2μϵ ε= ρ+ρε+με+2με^2ρ+μ+2με. The proof is relegated to Appendix C. Note that when δ=0δ=0, it implies the result in Ryu et al. (2019). 4.2 Convergence Under Weakly Non-Expansiveness with AC-DC In this subsection, we will show that the AC-DC denoiser satisfies Assumption 1 under mild conditions. To this end, consider the following: Assumption 2 (Smoothness of log⁡pdata p_ data). The log data density log⁡pdata p_ data is M-smooth for a constant M>0M>0, i.e., ‖∇log⁡pdata()−∇log⁡pdata()‖2≤M‖−‖2 \|∇ p_ data( x)-∇ p_ data( y) \|_2≤ M \| x- y \|_2 for all ,∈ x, y . Assumption 3 (Coercivity for −log⁡pdata- p_ data). There exists constants c1>0c_1>0 and c2≥0c_2≥ 0 such that ‖∇log⁡pdata()‖22≥−c1log⁡pdata−c2,‖2≤−c1log⁡pdata()+c2,∀∈ \|∇ p_ data( x) \|_2^2≥-c_1 p_ data-c_2,\; \| x \|_2≤-c_1 p_ data( x)+c_2,\;∀ x (13) This coercivity assumption means the negative log-density grows sufficiently fast at infinity, which prevents the Langevin dynamics from "escaping to infinity". This assumption guarantees stability and ensure ergodicity leading to convergence to the stationary distribution (Mattingly et al., 2002). Theorem 2. Suppose that the assumptions in Theorem 1, Assumption 2 and Assumption 3 hold. Further, assume that the DC step reaches the stationary distribution for each k. Let Dσ(k):~(k)↦tw(k)D_σ^(k): z^(k) z^(k)_ tw denote the AC-DC denoiser. Then, we have: (a) With probability at least 1−2e−νk1-2e^- _k, the following holds for iteration k of ADMM-PnP: ‖(Dσ(k)−I)()−(Dσ(k)−I)()‖22≤ϵk2‖−‖22+δk2 \|(D_σ^(k)-I)( x)-(D_σ^(k)-I)( y) \|_2^2≤ε^2_k \| x- y \|_2^2+δ^2_k (14) for any ,∈ x, y and k∈ℕ+k ^+ when σ(k)2+(σ(k))2<1/M _ s^(k)^2+(σ^(k))^2<1/M with ϵk2 ε^2_k =3((2Mσ(k)2/1−σ(k)2M)2+(σ(k))4M2) =3(( 2M _ s^(k)^21- _ s^(k)^2M)^2+(σ^(k))^4M^2) (15) δk2 δ^2_k =3(2(σ(k))2(d+2dνk+2νk)+32dσ(k)2/(1−Mσ(k)2)log⁡2/νk). =3(2(σ^(k))^2(d+2 d _k+2 _k)+ 32d _ s^(k)^2(1-M _ s^(k)^2) 2 _k). (16) In other words, with νk=ln⁡2π/6η+2ln⁡k _k= 2π6η+2 k, the denoiser Dσ(k)D_σ^(k) satisfies part (a) for all k∈ℕ+k ^+ with probability at least 1−η1-η. (b) Assume that σ(k)σ^(k) is scheduled such that limk→∞(σ(k))2νk=0 _k→∞(σ^(k))^2 _k=0 for νk=ln⁡2π/6η+2ln⁡k _k= 2π6η+2 k, ϵ<1ε<1, and ϵ/μ(1+ϵ−2ϵ2)<1/ρ εμ(1+ε-2ε^2)< 1ρ all hold, where ϵ=limk→∞supϵkε= _k→∞ _k with ϵk _k defined in (15). Consequently, δ=limk→∞supδkδ= _k→∞ _k is finite and ADMM-PnP with the AC-DC denoiser converges to an r-ball (see r in Theorem 1) with probability at least 1−η1-η. The proof is relegated to Appendix D. Theorem 2 (a) establishes that the AD-DC denoiser is weakly non-expansive with probability 1−2e−νk1-2e^- _k in iteration k. The (b) part states that when σ(k)σ^(k) is carefully scheduled to approach zero as k grows, then, with high probability, all the iterations satisfy the weakly non-expansiveness together—this leads to the convergence of the ADMM-PnP algorithm. 4.3 Convergence without Convexity of ℓ The weakly non-expansiveness based convergence analysis holds under fixed step size (i.e., ρ) of the ADMM-PnP algorithm, which is consistent with practical implementations in many cases. However, the assumption that the ℓ term is μ-strongly convex is only met by some inverse problems, e.g., signal denoising and deblurring, but not met by others such as signal compression/recovery and data completion. In this subsection, we remove the convexity assumption and analyze the AC-DC denoiser’s properties under the adaptive ρ-scheme following that in (Chan et al., 2016). Theorem 3. Suppose that Assumptions 2-3 hold. Let D:=diam()=sup,∈‖−‖2<∞,S:=inf∈∥∇log⁡pdata()∥2<∞D:= diam(X)= _ x, y \| x- y \|_2<∞,\;S:= _ x \|∇ p_ data( x) \|_2<∞ and define L:=MD+SL:=MD+S. Let Dσ(k):~(k)↦tw(k)D_σ^(k): z^(k) z^(k)_ tw denote the AC-DC denoiser. Also, assume that the DC step reaches the stationary distribution for each k111Note that Theorems 2 and 3 use this stationary distribution assumption for notation conciseness. For their counterparts removing this assumption, see Appendix E.2.. Then, the following hold: (a) (Boundedness) With probability at least 1−2e−νk1-2e^- _k, the denoiser Dσ(k)D_σ^(k) is bounded at each iteration k i.e. 1d‖(Dσ(k)−I)()‖22≤ck2 1d \|(D_σ^(k)-I)( x) \|_2^2≤ c_k^2 whenever σ(k)2+(σ(k))2<1/M, _ s^(k)^2+(σ^(k))^2<1/M, where ck=(σ(k))2(2+4νk+4νk)+16σ(k)2/1−Mσ(k)2log⁡2/νk+2σ(k)4L2+2(σ(k))4L2c_k=(σ^(k))^2(2+4 _k+4 _k)+ 16 _ s^(k)^21-M _ s^(k)^2 2 _k+2 _ s^(k)^4L^2+2(σ^(k))^4L^2, and νk>0 _k>0. Let νk=ln⁡2π26η+2ln⁡k _k= 2π^26η+2 k with η∈(0,1]η∈(0,1]. Consequently, the denoiser Dσ(k)D_σ^(k) is bounded for all k∈ℕ+k _+ with corresponding ckc_k and probability at least 1−η1-η. (b) (Convergence) Assume there exists R<∞R<∞ such that ‖∇ℓ()‖2/d≤R \|∇ ( x) \|_2/ d≤ R. Apply the ρ-increasing rule in (Chan et al., 2016) and schedule (σ(k),σ(k))(σ^(k), _ s^(k)) such that limk→∞(σ(k))2(2+4νk+4νk)=0,limk→∞σ(k)21−Mσ(k)2log⁡2νk=0, _k→∞(σ^(k))^2(2+4 _k+4 _k)=0,\; _k→∞ _ s^(k)^21-M _ s^(k)^2 2 _k=0, limk→∞σ(k)=0,limk→∞σ(k)=0σ(k)2+(σ(k))2<1/M,∀k∈ℕ+ _k→∞σ^(k)=0, _k→∞ _ s^(k)=0\; _ s^(k)^2+(σ^(k))^2<1/M,\;∀ k _+ for νk=ln⁡2π26η+2ln⁡k _k= 2π^26η+2 k with η∈(0,1]η∈(0,1]. Then, the solution sequence converges to a fixed point with probability at least 1−η1-η. The proof is relegated to Appendix E. Theorem 3 (a) shows that, with high probability the denoiser is bounded uniformly across all iterations k. Part (b) further shows that, under the proper scheduling of (σ(k),σ(k))(σ^(k), _ s^(k)), the AD-DC ADMM-PnP algorithm converges to a fixed point with high probability. The condition of D<∞D<∞ implies the data space X has bounded support, which is natural in practice: for images, pixel intensities typically lie within a bounded range such as [0,1][0,1]. Additionally, the condition S<∞S<∞ ensures that there exists at least one point in X where the score norm is finite. This prevents pathological cases where the score diverges everywhere (making the distribution degenerate). Together, these conditions guarantee that the score is “well-behaved”. A remark is that all theoretical results in this section focus on fixed-point convergence, which is not the strongest form of convergence guarantees. Establishing stronger convergence results, e.g., stationary-point convergence, for PnP approaches is considered challenging as the objective function is implicit (more specifically, h(⋅)h(·) in (2) is implicit). Nonetheless, in recent years, some efforts have been made towards establishing stationary-point convergence for PnP methods under certain types of denoisers (see, e.g., (Hurault et al., 2022a; b; Wei et al., 2025; Xu et al., 2025)); more discussions are in Sec. 5. 5 Related Works ADMM-PnP has gained much popularity, due to access to data-driven effective denoisers. It has been used in various applications like image restoration (Chan et al., 2016), data compression (Yuan et al., 2022), hyperspectral imaging (Liu et al., 2022), and medical imaging (Ahmad et al., 2020). Theoretical understanding of PnP algorithms with general “black-box” denoisers remains limited. Unlike classical proximal operators, the implicit regularizer h(⋅)h(·) in (2) handled by data-driven denoisers is typically unknown. Hence, many results are therefore restricted to fixed-point convergence; see, e.g., (Ryu et al., 2019; Chan et al., 2016). Nonetheless, in certain cases where the denoisers have interesting structures, stronger convergence results can be established. For example, (Xu et al., 2025) used classical results from image denoising connecting linear denoisers with quadratic h(⋅)h(·) to show that when linear denoisers are employed, ADMM-PnP converges to KKT points. Hurault et al. (2022a) showed stationary-point convergence of PnP methods gradient-type denoisers, leveraging the fact that this type of denoisers can be written as a proximal operator of a special function (see the nonconvex counterpart in (Hurault et al., 2022b)); Wei et al. (2025) trained denoisers to satisfy a cocoercive conservativity condition, which also ensures convergence of PnP to stationary points associated with an implicit convex h(⋅)h(·). Nonetheless, these results do not cover diffusion-based denoisers. In this work, we generalize the fixed-point convergence proofs in (Chan et al., 2016; Ryu et al., 2019) to accommodate the diffusion score-based AC–DC denoiser. Recent advances in score-based generative modeling have motivated their integration into PnP algorithms. One line of work directly replaces the proximal denoiser with a pre-trained scores (Zhu et al., 2023; Li et al., 2024). Alternatively, others embed the score function as an explicit regularizer with task-specific loss (Mardani et al., 2024; Renaud et al., 2024a). Deterministic version PnP have also been considered. For example, Wang et al. (2024); Song et al. (2023) use unrolled ODE and consistency model-distilled one-step representation to express the target signal, respectively. These methods are similar to (Bora et al., 2017), but with diffusion-driven parameterization. Prior works have emphasized the importance of matching the residual noise to the operating range of the PnP denoiser. D-AMP (Metzler et al., 2016; Eksioglu and Tanc, 2018) achieve this via the Onsager correction, which approximately Gaussianizes the residual under compressive sensing problem structures. Wei et al. (2021) learns a reinforcement learning-based policy to automatically tune all internal parameters, including denoising strength. Unlike AC–DC that provides a generic correction mechanism for a variety inverse problems, these methods either problem specific or require training additional models. Score-based inverse problem solvers have also attempted to “bring” iterates to noisy data manifolds used during training. The work Chung et al. (2022) uses a manifold constraint based on gradient of data-fidelity, while He et al. (2024) uses an off-the-shelf pretrained neural network to impose a manifold constraint. On the other hand, Zirvi et al. (2025) uses the projection of measurement guidance to low-rank subspace, using SVD on the intermediate diffusion state, for similar purposes. The idea of adding noise before evaluating score functions during optimization procedures (similar to our AC step) has been widely considered (Li et al., 2024; Graikos et al., 2022; Renaud et al., 2024b; Mardani et al., 2024). A variant of this called estimation-correction idea proposed in (Karras et al., 2022) is used in (Zhu et al., 2023) for this purpose. 6 Experiments Figure 2: Inpainting under random missings. Figure 3: Inpainting under box missing. Dataset and Evaluation Metrics. For all these tasks, we use two datasets: FFHQ 256×256256× 256 (Karras et al., 2021) and ImageNet 256×256256× 256 (Deng et al., 2009). During testing, we randomly sample 100 images from the validation set of each dataset. All the methods use the pre-trained score model in Chung et al. (2023). We use Peak Signal-to-Noise Ratio (PSNR) as a pixel-wise similarity metric, and Structural Similarity Index (SSIM) and Learned Perceptual Image Patch Similarity (LPIPS) (Zhang et al., 2018) as perceptual similarity metrics. We report these metrics averaged over the 100 test images for each method and inverse problem. Task Description. We consider ∼(,σn2) ξ ( 0, _n^2 I) with σn=0.05 _n=0.05 for all the tasks. (a) For super-resolution, we use cubic interpolation method with kernel size 4 for downsampling the resolution by 44 times. (b) For recovery under Gaussian blurring (Gaussian deblurring), a kernel of size 61 and standard deviation 3 is used. (c) As for recovery under motion blurring (motion deblurring), a kernel of size 61 and standard deviation of 0.5 is used. (d) In inpainting under box mask (box inpainting), an approximately centered mask of size 128×128128× 128 is sampled in image while maintaining the 3232 pixel margin in both spatial dimensions of the input image. (e) For inpainting under random missings (random inpainting), 70%70\% of the pixels are uniformly sampled to be masked, and a scaling of 22 was used in high dynamic range (HDR) before clipping the values. (f) For phase retrieval, similar as in prior works (Wu et al., 2024; Mardani et al., 2024), we use oversampling by factor of 22. g) For deblurring under nonlinear blurring, we use the operator in (Tran et al., 2021) with default settings. Baselines. We use a set of baselines, namely, DPS (Chung et al., 2023), DAPS (Zhang et al., 2024), DDRM (Kawar et al., 2022), DiffPIR (Zhu et al., 2023), RED-diff (Mardani et al., 2024), DPIR (Zhang et al., 2022), DCDP (Li et al., 2025), PMC (Sun et al., 2024). Hyperparameter Settings. We adopt a linear schedule for σ(k)σ^(k) with range [0.1,10][0.1,10] over W decay window i.e. σ(k)=max⁡(0.1,10−(10−0.1)⋅k/W)σ^(k)= (0.1,10-(10-0.1)· k/W). The maximal number of iterations for our proposed method is set to K=W+10K=W+10. At iteration k, we use J=10J=10 DC steps, and the schedules η(k)=5×10−4σ(k)η^(k)=5× 10^-4σ^(k) and σ(k)=0.1/σ(k) _ s^(k)= 0.1 σ^(k). We use gradient descent with Adam optimizer (Kingma and Ba, 2015) for solving each regularized maximum likelihood subproblem (7a). This subproblem is optimized for maximum of 10001000 iterations with convergence detected when the loss value increases more than Δtol=1×10−1 _ tol=1× 10^-1 consecutively for 33 iterations window. We conduct our experiment with two variants based on the third stage: using Tweedie’s lemma (denoted as “Ours-tweedie”) and a 10-step ODE based denoiser (Zhang et al., 2024; Karras et al., 2022) (denoted as “Ours-ode”). We use the preconditioning in Karras et al. (2022) while using the pretrained diffusion models. Qualitative Performance. Figs 2, 3 and 5 show reconstructions under inpainting under random missings, inpainting under box missing, and motion deblurring. It can be seen that our method is able to recover the image that is comparatively natural looking with less noise and artifacts, while being consistent with the measurements. On the other hand, images recovered with DiffPIR appears to suffer from noise and artifacts, whereas DPS leads to measurement-inconsistent reconstructions. Our method outperforms others while other methods appear to either be blurred or contain noisy artifacts in the recovered image. Recovery by DPS is less consistent with the original image; the pattern on the child’s clothing is completely lost. Quantitative Performance. Table 1 summarizes PSNR, SSIM and LPIPS averaged over 100 images on FFHQ and Imagenet datasets. In almost all of the inverse problems, both of our variants (Ours-tweedie and Ours-ode) achieve the best or second-best performance in terms of all metrics. Our method significantly outperforms other PnP baseline methods considered, namely, DDRM, DiffPIR and RED-diff. This demonstrates the effectiveness of our AC-DC denoiser. Figure 4: Recovery under motion blurring. Figure 5: Influence of DC steps in the denoiser. Table 1: Reconstruction metrics (100 images) on FFHQ / ImageNet. Bold: best, blue: 2nd best. FFHQ ImageNet Task Method PSNR↑ SSIM↑ LPIPS↓ PSNR↑ SSIM↑ LPIPS↓ Superresolution (4×) Ours–tweedie 30.439 0.857 0.178 27.318 0.717 0.280 Ours–ode 29.991 0.845 0.156 26.919 0.700 0.276 DAPS 29.529 0.814 0.167 26.653 0.680 0.266 DPS 24.828 0.705 0.257 22.785 0.549 0.411 DDRM 27.145 0.782 0.261 26.105 0.683 0.306 DiffPIR 26.771 0.749 0.208 23.884 0.543 0.336 RED-diff 16.833 0.422 0.547 18.662 0.309 0.519 DPIR 28.849 0.826 0.254 26.524 0.699 0.334 DCDP 27.761 0.639 0.332 24.517 0.525 0.361 PMC 23.774 0.421 0.407 22.534 0.334 0.456 Inpainting (Random) Ours–tweedie 32.844 0.906 0.122 29.564 0.817 0.184 Ours–ode 32.127 0.894 0.095 28.733 0.795 0.148 DAPS 31.652 0.847 0.124 28.137 0.751 0.162 DPS 29.084 0.828 0.181 26.049 0.678 0.318 DDRM 28.969 0.847 0.178 27.883 0.778 0.203 DiffPIR 28.558 0.709 0.230 26.923 0.639 0.222 RED-diff 20.361 0.630 0.275 20.948 0.464 0.315 PMC 23.289 0.755 0.263 25.965 0.636 0.342 Motion Deblur Ours–tweedie 30.003 0.854 0.179 27.149 0.717 0.280 Ours–ode 29.648 0.841 0.154 26.615 0.694 0.275 DAPS 29.051 0.815 0.175 26.571 0.689 0.276 DPS 23.257 0.663 0.265 19.613 0.451 0.451 PMC 19.480 0.590 0.426 21.608 0.480 0.510 FFHQ ImageNet Task Method PSNR↑ SSIM↑ LPIPS↓ PSNR↑ SSIM↑ LPIPS↓ Gaussian Blur Ours-tweedie 30.402 0.853 0.175 27.199 0.705 0.281 Ours-ode 30.019 0.841 0.158 26.899 0.690 0.282 DAPS 29.790 0.813 0.157 26.886 0.678 0.260 DPS 26.106 0.730 0.207 23.995 0.575 0.328 DiffPIR 25.148 0.699 0.230 22.756 0.508 0.374 DPIR 28.875 0.833 0.228 26.702 0.700 0.314 DCDP 16.821 0.171 0.642 15.102 0.136 0.620 PMC 20.172 0.638 0.344 24.103 0.545 0.415 Inpainting (Box) Ours-tweedie 24.025 0.859 0.131 21.626 0.789 0.222 Ours-ode 23.342 0.837 0.136 20.618 0.743 0.227 DAPS 23.643 0.815 0.146 21.303 0.774 0.199 DPS 23.488 0.817 0.164 19.933 0.677 0.309 DiffPIR 20.934 0.561 0.294 19.565 0.562 0.342 RED-diff 18.713 0.523 0.364 18.075 0.499 0.371 DCDP 25.230 0.754 0.163 20.991 0.727 0.195 PMC 14.828 0.697 0.318 15.550 0.666 0.326 Phase Retrieval Ours-tweedie 27.944 0.793 0.209 17.770 0.440 0.471 Ours-ode 27.095 0.757 0.237 16.013 0.339 0.539 DAPS 26.707 0.749 0.230 16.444 0.395 0.512 DPS 11.627 0.366 0.658 9.434 0.216 0.768 RED-diff 15.411 0.490 0.480 12.852 0.204 0.695 DCDP 20.026 0.540 0.424 12.257 0.212 0.665 PMC 10.421 0.287 0.783 8.636 0.129 0.890 Effectiveness of DC. To perform ablation study on the DC stage, we consider the challenging phase retrieval problem. Fig. 5 shows the output of ADMM-PnP with our AC-DC denoiser with different numbers of DC iterations J. With J=0J=0 (disabling DC step), artifacts remain severe. Increasing J progressively results in cleaner images. More Details and Additional Experiments. More details and experiments are in appendices. 7 Conclusion We introduced the AC-DC denoiser, a score-based denoiser designed for integration within the ADMM-PnP framework. The denoiser adopts a three-stage structure aimed at mitigating the mismatch between ADMM iterates and the noisy manifolds on which score functions are trained. We established convergence guarantees for ADMM-PnP with the AC-DC denoiser under both fixed and adaptive step size schedules. Empirical results across a range of inverse problems demonstrate that the proposed method consistently improves solution quality over existing baselines. Limitations. While our analysis provides initial insights, several aspects merit deeper understanding. The second convergence result relaxes convexity by allowing adaptive step sizes, though such schedules are arguably less appealing in practice. Our experiments, however, suggest that constant step sizes also perform well for nonconvex objectives; it is therefore desirable to establish convergence guarantees for constant step sizes in such settings. In addition, our result ensures the stability of the ADMM method, but does not directly explain the reason why the AC–DC denoiser attains high-quality recovery; recoverability and estimation error analyses are also desirable. On the implementation side, the noise schedules used in the AC and DC stages are currently guided by empirical heuristics. Designing problem-adaptive scheduling strategies may further improve both convergence speed and robustness. Additionally, each iteration of AC–DC denoiser needs multiple score evaluations. Reducing the required NFEs could significantly improve its efficiency. Ethics Statement: This work focuses exclusively on the theory and methodology of solving inverse problems. It does not involve human subjects, personal data, or any sensitive procedures. Reproducibility Statement: The source code is provided as a part of the supplementary material. All assumptions, derivations and necessary details regarding the theory and experiments are included in the appendices. Acknowledgment: This work was supported in part by the National Science Foundation (NSF) under Project NSF CCF-2210004. References R. Ahmad, C. A. Bouman, G. T. Buzzard, S. Chan, S. Liu, E. T. Reehorst, and P. Schniter (2020) Plug-and-play methods for magnetic resonance imaging: using denoisers for image recovery. IEEE Signal Processing Magazine 37 (1), p. 105–116. External Links: ISSN 1558-0792, Link, Document Cited by: §5. I. Alkhouri, S. Liang, E. Bell, Q. Qu, R. Wang, and S. Ravishankar (2024) Image reconstruction via autoencoding sequential deep image prior. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, External Links: Link Cited by: §1. I. Alkhouri, S. Liang, R. Wang, Q. Qu, and S. Ravishankar (2023) Robust physics-based deep mri reconstruction via diffusion purification. External Links: 2309.05794, Link Cited by: Appendix B, §2, §3. S. R. Arridge (1999) Optical tomography in medical imaging. Inverse problems 15 (2), p. R41. Cited by: §1. H. H. Bauschke, P. L. Combettes, H. H. Bauschke, and P. L. Combettes (2017) Convex analysis and monotone operator theory in hilbert spaces. Springer. Cited by: §A.3. A. F. Bennett (1992) Inverse methods in physical oceanography. Cambridge university press. Cited by: §1. A. Bora, A. Jalal, E. Price, and A. G. Dimakis (2017) Compressed sensing using generative models. In International conference on machine learning, p. 537–546. Cited by: §5. S. P. Boyd and L. Vandenberghe (2004) Convex optimization. Cambridge university press. Cited by: §A.2. H. J. Brascamp and E. H. Lieb (1976) On extensions of the brunn-minkowski and prékopa-leindler theorems, including inequalities for log concave functions, and with an application to the diffusion equation. Journal of Functional Analysis 22 (4), p. 366–389. External Links: ISSN 0022-1236, Document, Link Cited by: Appendix D. F. E. Browder (1965) Nonexpansive nonlinear operators in a banach space. Proceedings of the National Academy of Sciences of the United States of America 54 (4), p. 1041–1044. External Links: ISSN 00278424, 10916490, Link Cited by: §A.2. G. T. Buzzard, S. H. Chan, S. Sreehari, and C. A. Bouman (2018) Plug-and-play unplugged: optimization-free reconstruction using consensus equilibrium. SIAM Journal on Imaging Sciences 11 (3), p. 2001–2020. Cited by: §4.1. S. H. Chan (2019) Performance analysis of plug-and-play admm: a graph signal processing perspective. (arXiv:1809.00020). Note: arXiv:1809.00020 [eess] External Links: Link, Document Cited by: §4.1. S. Chan, X. Wang, and O. Elgendy (2016) Plug-and-play admm for image restoration: fixed point convergence and applications. IEEE Transactions on Computational Imaging P, p. . External Links: Document Cited by: §E.1, §1, §3, §4.3, §5, §5, Theorem 3. X. Chen, S. S. Du, and X. T. Tong (2020) On stationary-point hitting time and ergodicity of stochastic gradient langevin dynamics. Journal of Machine Learning Research 21 (68), p. 1–41. External Links: Link Cited by: Appendix D. H. Chung, J. Kim, M. T. Mccann, M. L. Klasky, and J. C. Ye (2023) Diffusion posterior sampling for general noisy inverse problems. In International Conference on Learning Representations, External Links: Link Cited by: 2nd item, §H.2, §1, §1, §2, §6, §6. H. Chung, B. Sim, D. Ryu, and J. C. Ye (2022) Improving diffusion models for inverse problems using manifold constraints. Advances in Neural Information Processing Systems 35, p. 25683–25696. Cited by: §5. B. Combal, F. Baret, M. Weiss, A. Trubuil, D. Macé, A. Pragnere, R. Myneni, Y. Knyazikhin, and L. Wang (2003) Retrieval of canopy biophysical variables from bidirectional reflectance: using prior information to solve the ill-posed inverse problem. Remote sensing of environment 84 (1), p. 1–15. Cited by: §1. P. L. Combettes and I. Yamada (2015) Compositions and convex combinations of averaged nonexpansive operators. Journal of Mathematical Analysis and Applications 425 (1), p. 55–70. External Links: ISSN 0022-247X, Document, Link Cited by: §A.2, §A.2, §A.3, §A.3. K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian (2007) Image denoising by sparse 3-d transform-domain collaborative filtering. IEEE Transactions on Image Processing 16 (8), p. 2080–2095. External Links: Document Cited by: §1. A. S. Dalalyan and A. G. Karagulyan (2019) User-friendly guarantees for the langevin monte carlo with inaccurate gradient. Stochastic Processes and their Applications 129 (12), p. 5278–5311. Note: arXiv:1710.00095 [math] External Links: ISSN 03044149, Document Cited by: §E.2. J. Deng, W. Dong, R. Socher, L. Li, K. Li, and L. Fei-Fei (2009) ImageNet: a large-scale hierarchical image database. In 2009 IEEE Conference on Computer Vision and Pattern Recognition, Vol. , p. 248–255. External Links: Document Cited by: §6. A. Dytso, H. V. Poor, and S. Shamai (2021) A general derivative identity for the conditional mean estimator in gaussian noise and some applications. External Links: 2104.01883, Link Cited by: Appendix D. J. Eckstein and D. P. Bertsekas (1992) On the douglas—rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical programming 55, p. 293–318. Cited by: §A.2. E. M. Eksioglu and A. K. Tanc (2018) Denoising amp for mri reconstruction: bm3d-amp-mri. SIAM Journal on Imaging Sciences 11 (3), p. 2090–2109. External Links: Document Cited by: §5. M. Elad and M. Aharon (2006) Image denoising via sparse and redundant representations over learned dictionaries. IEEE Transactions on Image Processing 15 (12), p. 3736–3745. External Links: Document Cited by: §1. D. Entekhabi, H. Nakamura, and E. G. Njoku (1994) Solving the inverse problem for soil moisture and temperature profiles by sequential assimilation of multifrequency remotely sensed observations. IEEE Transactions on Geoscience and Remote Sensing 32 (2), p. 438–448. Cited by: §1. P. Giselsson (2015) Tight global linear convergence rate bounds for douglas–rachford splitting. Journal of Fixed Point Theory and Applications 19, p. 2241 – 2270. External Links: Link Cited by: Lemma 4. A. Graikos, N. Malkin, N. Jojic, and D. Samaras (2022) Diffusion models as plug-and-play priors. In Advances in Neural Information Processing Systems, A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho (Eds.), External Links: Link Cited by: §5. C. Hatsell and L. Nolte (1971) Some geometric properties of the likelihood ratio (corresp.). IEEE Transactions on Information Theory 17 (5), p. 616–618. External Links: Document Cited by: Appendix D. Y. He, N. Murata, C. Lai, Y. Takida, T. Uesaka, D. Kim, W. Liao, Y. Mitsufuji, J. Z. Kolter, R. Salakhutdinov, and S. Ermon (2024) Manifold preserving guided diffusion. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: §5. W. Hu, D. Tao, W. Zhang, Y. Xie, and Y. Yang (2017) The twist tensor nuclear norm for video completion. IEEE Transactions on Neural Networks and Learning Systems 28 (12), p. 2961–2973. External Links: Document Cited by: §1. S. Hurault, A. Leclaire, and N. Papadakis (2022a) Gradient step denoiser for convergent plug-and-play. (arXiv:2110.03220). Note: arXiv:2110.03220 [cs] External Links: Link, Document Cited by: §4.3, §5. S. Hurault, A. Leclaire, and N. Papadakis (2022b) Proximal denoiser for convergent plug-and-play optimization with nonconvex regularization. In Proceedings of the 39th International Conference on Machine Learning, p. 9483–9505 (en). External Links: ISSN 2640-3498, Link Cited by: §4.3, §5. K. H. Jin, M. T. McCann, E. Froustey, and M. Unser (2017) Deep convolutional neural network for inverse problems in imaging. IEEE transactions on image processing 26 (9), p. 4509–4522. Cited by: §1. T. Karras, M. Aittala, S. Laine, and T. Aila (2022) Elucidating the design space of diffusion-based generative models. In Proceedings of the 36th International Conference on Neural Information Processing Systems, NIPS ’22, Red Hook, NY, USA. External Links: ISBN 9781713871088 Cited by: §3, §5, §6, Lemma 20, Lemma 21, Lemma 9, Remark 1. T. Karras, S. Laine, and T. Aila (2021) A style-based generator architecture for generative adversarial networks. IEEE Trans. Pattern Anal. Mach. Intell. 43 (12), p. 4217–4228. External Links: ISSN 0162-8828, Link, Document Cited by: §6. B. Kawar, M. Elad, S. Ermon, and J. Song (2022) Denoising diffusion restoration models. In Advances in Neural Information Processing Systems, A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho (Eds.), External Links: Link Cited by: 1st item, §1, §2, §6. B. Kawar, G. Vaksman, and M. Elad (2021) SNIPS: solving noisy inverse problems stochastically. In Advances in Neural Information Processing Systems, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan (Eds.), External Links: Link Cited by: 1st item, §2. D. P. Kingma and J. Ba (2015) Adam: A method for stochastic optimization. In 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, Y. Bengio and Y. LeCun (Eds.), External Links: Link Cited by: §6. B. Laurent and P. Massart (2000) Adaptive estimation of a quadratic functional by model selection. The Annals of Statistics 28 (5), p. 1302–1338. External Links: ISSN 0090-5364, 2168-8966, Document Cited by: Appendix D, Appendix E. X. Li, S. M. Kwon, I. R. Alkhouri, S. Ravishankar, and Q. Qu (2024) Decoupled data consistency with diffusion purification for image restoration. (arXiv:2403.06054) (en). Note: arXiv:2403.06054 [eess] External Links: Link, Document Cited by: §1, §2, §3, §5, §5. X. Li, S. M. Kwon, S. Liang, I. R. Alkhouri, S. Ravishankar, and Q. Qu (2025) Decoupled data consistency with diffusion purification for image restoration. External Links: 2403.06054, Link Cited by: 7th item, §6. K. Liang (2007) Homocentric convergence ball of the secant method. Applied Mathematics-A Journal of Chinese Universities 22, p. 353–365. Cited by: §4.1. Y. Liu, X. Zhao, Y. Zheng, T. Ma, and H. Zhang (2022) Hyperspectral image restoration by tensor fibered rank constrained optimization and plug-and-play regularization. IEEE Transactions on Geoscience and Remote Sensing 60 (), p. 1–17. External Links: Document Cited by: §5. M. Mardani, J. Song, J. Kautz, and A. Vahdat (2024) A variational perspective on solving inverse problems with diffusion models. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: 4th item, §1, §1, §2, §3, §5, §5, §6, §6. J. C. Mattingly, A. M. Stuart, and D. J. Higham (2002) Ergodicity for sdes and approximations: locally lipschitz vector fields and degenerate noise. Stochastic Processes and their Applications 101 (2), p. 185–232. External Links: ISSN 0304-4149, Document Cited by: Appendix D, §4.2. C. Meng, Y. He, Y. Song, J. Song, J. Wu, J. Zhu, and S. Ermon (2022) SDEdit: guided image synthesis and editing with stochastic differential equations. In International Conference on Learning Representations, External Links: Link Cited by: §2. C. A. Metzler, A. Maleki, and R. G. Baraniuk (2016) From denoising to compressed sensing. (arXiv:1406.4175) (en). Note: arXiv:1406.4175 [cs] External Links: Link, Document Cited by: §5. W. Nie, B. Guo, Y. Huang, C. Xiao, A. Vahdat, and A. Anandkumar (2022) Diffusion models for adversarial purification. In Proceedings of the 39th International Conference on Machine Learning, K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu, and S. Sabato (Eds.), Proceedings of Machine Learning Research, Vol. 162, p. 16805–16827. External Links: Link Cited by: §2, §3. D.P. Palomar and S. Verdu (2006) Gradient of mutual information in linear vector gaussian channels. IEEE Transactions on Information Theory 52 (1), p. 141–154. External Links: Document Cited by: Appendix D. L. Pardo (2018) Statistical inference based on divergence measures. Chapman and Hall/CRC. Cited by: Lemma 5. X. Peng, Z. Zheng, W. Dai, N. Xiao, C. Li, J. Zou, and H. Xiong (2024) Improving diffusion models for inverse problems using optimal posterior covariance. External Links: 2402.02149, Link Cited by: Appendix D. M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics 378, p. 686–707. Cited by: §1. H. Ren and I. K. Argyros (2021) On the complexity of extending the convergence ball of wang’s method for finding a zero of a derivative. Journal of Complexity 64, p. 101526. External Links: ISSN 0885-064X, Document, Link Cited by: §4.1. H. Ren and Q. Wu (2009) Convergence ball and error analysis of a family of iterative methods with cubic convergence. Applied Mathematics and Computation 209 (2), p. 369–378. External Links: ISSN 0096-3003, Document, Link Cited by: §4.1. M. Renaud, J. Prost, A. Leclaire, and N. Papadakis (2024a) Plug-and-play image restoration with stochastic denoising regularization. In Proceedings of the 41st International Conference on Machine Learning, ICML’24. Cited by: Appendix B, §2, §2, §3, §5. M. Renaud, J. Prost, A. Leclaire, and N. Papadakis (2024b) Plug-and-play image restoration with stochastic denoising regularization. CoRR abs/2402.01779. External Links: Link Cited by: §1, §1, §5. H. Robbins (1992) An empirical bayes approach to statistics. Proceedings of the Third Berkeley Symposium on Mathematical and Statistical Probability 1, p. . External Links: ISBN 978-0-387-94037-3, Document Cited by: §2, Lemma 1. E. Ryu, J. Liu, S. Wang, X. Chen, Z. Wang, and W. Yin (2019) Plug-and-play methods provably converge with properly trained denoisers. In Proceedings of the 36th International Conference on Machine Learning, K. Chaudhuri and R. Salakhutdinov (Eds.), Proceedings of Machine Learning Research, Vol. 97, p. 5546–5557. External Links: Link Cited by: §A.3, Appendix C, §C.1, §1, §3, §4.1, §4.1, §5, Lemma 7. O. Semerci, N. Hao, M. E. Kilmer, and E. L. Miller (2014) Tensor-based formulation and nuclear norm regularization for multienergy computed tomography. IEEE Transactions on Image Processing 23 (4), p. 1678–1693. External Links: Document Cited by: §1. V. Shah and C. Hegde (2018) Solving linear inverse problems using gan priors: an algorithm with provable guarantees. In 2018 IEEE international conference on acoustics, speech and signal processing (ICASSP), p. 4609–4613. Cited by: §1. B. Song, S. M. Kwon, Z. Zhang, X. Hu, Q. Qu, and L. Shen (2023) Solving inverse problems with latent diffusion models via hard data consistency. (en). External Links: Link Cited by: §1, §5. J. Song, A. Vahdat, M. Mardani, and J. Kautz (2022a) Pseudoinverse-guided diffusion models for inverse problems. (en). External Links: Link Cited by: §2. Y. Song, L. Shen, L. Xing, and S. Ermon (2022b) Solving inverse problems in medical imaging with score-based generative models. In International Conference on Learning Representations, External Links: Link Cited by: §1, §1. Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole (2021) Score-based generative modeling through stochastic differential equations. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021, External Links: Link Cited by: §2. Y. Sun, B. Wohlberg, and U. S. Kamilov (2019) An online plug-and-play algorithm for regularized image reconstruction. IEEE Transactions on Computational Imaging 5 (3), p. 395–408. Note: arXiv:1809.04693 [cs] External Links: ISSN 2333-9403, 2334-0118, 2573-0436, Document Cited by: §4.1. Y. Sun, Z. Wu, Y. Chen, B. T. Feng, and K. L. Bouman (2024) Provable probabilistic imaging using score-based generative priors. External Links: 2310.10835, Link Cited by: 8th item, §6. A. Tarantola (2005) Inverse problem theory and methods for model parameter estimation. SIAM. Cited by: §1. A. M. Teodoro, J. M. Bioucas-Dias, and M. A. T. Figueiredo (2017) Scene-adapted plug-and-play algorithm with convergence guarantees. (arXiv:1702.02445). Note: arXiv:1702.02445 [cs] External Links: Link, Document Cited by: §4.1. E. C. Titchmarsh and D. R. Heath-Brown (1986) The theory of the riemann zeta-function. Oxford university press. Cited by: §D.1. P. Tran, A. Tran, Q. Phung, and M. Hoai (2021) Explore image deblurring via encoded blur kernel space. 2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), p. 11951–11960. External Links: Link Cited by: §6. G. E. Uhlenbeck and L. S. Ornstein (1930) On the theory of the brownian motion. Phys. Rev. 36, p. 823–841. External Links: Document, Link Cited by: Appendix D. D. Ulyanov, A. Vedaldi, and V. Lempitsky (2020) Deep image prior. International Journal of Computer Vision 128 (7), p. 1867–1888. External Links: ISSN 1573-1405, Link, Document Cited by: §1. C. Villani et al. (2008) Optimal transport: old and new. Vol. 338, Springer. Cited by: §E.2, §E.2. M. J. Wainwright (2019) High-dimensional statistics: a non-asymptotic viewpoint. Vol. 48, Cambridge university press. Cited by: Appendix D. H. Wang, X. Zhang, T. Li, Y. Wan, T. Chen, and J. Sun (2024) DMPlug: a plug-in method for solving inverse problems with diffusion models. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, External Links: Link Cited by: Appendix B, §1, §2, §5. Y. Wang, J. Yu, and J. Zhang (2023) Zero-shot image restoration using denoising diffusion null-space model. In The Eleventh International Conference on Learning Representations, External Links: Link Cited by: §2. D. Wei, P. Chen, H. Xu, J. Yao, F. Li, and T. Zeng (2025) Learning cocoercive conservative denoisers via helmholtz decomposition for poisson inverse problems. (arXiv:2505.08909). Note: arXiv:2505.08909 [cs] External Links: Link, Document Cited by: §4.3, §5. K. Wei, A. Aviles-Rivero, J. Liang, Y. Fu, H. Huang, and C. Schönlieb (2021) TFPnP: tuning-free plug-and-play proximal algorithm with applications to inverse imaging problems. (arXiv:2012.05703). Note: arXiv:2012.05703 [cs] External Links: Link, Document Cited by: §5. Z. Wu, Y. Sun, Y. Chen, B. Zhang, Y. Yue, and K. Bouman (2024) Principled probabilistic imaging using diffusion models as plug-and-play priors. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, External Links: Link Cited by: §6. Z. Xiao, K. Kreis, and A. Vahdat (2022) Tackling the generative learning trilemma with denoising diffusion GANs. In International Conference on Learning Representations, External Links: Link Cited by: §1. L. Xu, L. Cheng, J. Chen, W. Pu, and X. Fu (2025) Radio map estimation via latent domain plug-and-play denoising. External Links: 2501.13472, Link Cited by: §4.3, §5. J. Yang, J. Wright, T. S. Huang, and Y. Ma (2010) Image super-resolution via sparse representation. IEEE Transactions on Image Processing 19 (11), p. 2861–2873. External Links: Document Cited by: §1. Y. Yao, H. Zhou, and Y. Liou (2008) Iterative algorithms with variable anchors for non-expansive mappings. Journal of Applied Mathematics and Computing 28, p. 39–49. Cited by: §A.2, §A.2. X. Yuan, Y. Liu, J. Suo, F. Durand, and Q. Dai (2022) Plug-and-play algorithms for video snapshot compressive imaging. IEEE Transactions on Pattern Analysis and Machine Intelligence 44 (10), p. 7093–7111. External Links: Document Cited by: §5. B. Zhang, W. Chu, J. Berner, C. Meng, A. Anandkumar, and Y. Song (2024) Improving diffusion inverse problem solving with decoupled noise annealing. CoRR abs/2407.01521. External Links: Link Cited by: Appendix F, 5th item, §6, §6. K. Zhang, Y. Li, W. Zuo, L. Zhang, L. Van Gool, and R. Timofte (2022) Plug-and-play image restoration with deep denoiser prior. IEEE Transactions on Pattern Analysis and Machine Intelligence 44 (10), p. 6360–6376. External Links: Document Cited by: 6th item, §6. R. Zhang, P. Isola, A. A. Efros, E. Shechtman, and O. Wang (2018) The Unreasonable Effectiveness of Deep Features as a Perceptual Metric . In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), Vol. , Los Alamitos, CA, USA, p. 586–595. External Links: ISSN , Document, Link Cited by: §I.3, §6. Y. Zhu, K. Zhang, J. Liang, J. Cao, B. Wen, R. Timofte, and L. V. Gool (2023) Denoising diffusion models for plug-and-play image restoration. In IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2023 - Workshops, Vancouver, BC, Canada, June 17-24, 2023, p. 1219–1229. External Links: Link, Document Cited by: Appendix B, 3rd item, §1, §1, §2, §3, §5, §5, §6. R. Zirvi, B. Tolooshams, and A. Anandkumar (2025) Diffusion state-guided projected gradient for inverse problems. In The Thirteenth International Conference on Learning Representations, External Links: Link Cited by: §5. Appendix A Preliminaries A.1 Notation Table 2: Summary of notation Symbol Description ∈ℝd x ^d The unknown signal or image to be recovered ∈ℝn y ^n Measurements, with n≤dn≤ d yiy_i iith element of measurement y :ℝd→ℝnA:R^d ^n Measurement operator X Support of x tX_t Support of t x_t ξ Additive measurement noise t x_t noisy data by using forward diffusion process with noise σ(t)σ(t) σ(k) x_σ^(k) noisy data by using forward diffusion process with noise σ(k)σ^(k) ℓ(||()) ( y||A( x)) data-fidelity loss (e.g. ‖−()‖22 \| y-A( x) \|_2^2) h()h( x) structural regularization prior (enforced via denoiser) ρ>0ρ>0 ADMM penalty parameter Prox(⋅·) Proximal operator (k) u^(k) Scaled dual variable in iteration k of ADMM (k) z^(k) Auxiliary variable in iteration k of ADMM ~(k)=(k+1)+(k) z^(k)= x^(k+1)+ u^(k) Pre-denoising input to the PnP denoiser σ(k)σ^(k) Noise level schedule for the AC-DC denoiser σ(k) _ s^(k) Variance parameter in the AC-DC prior for directional correction. ∼(,) n ( 0, I) Multivariate standard gaussian random variable (,σ)≈∇log⁡p(+σ) s_ θ( x,σ)≈ _ x p( x+σ n) Pretrained score function K Maximum iteration of ADMM J Total iteration of directional correction at for each denoising M Smoothness constant of ∇log⁡pdata∇ p_ data (Assumption 2) MtM_t Smoothness constant of ∇log⁡pt∇ p_t T>0T>0 Maximum time steps used for diffusion ℳσ(t)M_σ(t) Manifold of t x_t ℳσ(k)M_σ^(k) Manifold of t x_t where t∈[0,T]t∈[0,T] such that σ(t)=σ(k)σ(t)=σ^(k) Dσ(k)()D_σ^(k)( z) AC-DC denoiser at kkth iteration Rσ(k)()=Dσ(k)()−R_σ^(k)( z)=D_σ^(k)( z)- z Residual of AC-DC denoiser at kkth iteration I()=I( x)= x Identity mapping function I Identity matrix 0 vector of values 0 T1∘T2()=T1(T2())T_1 T_2( z)=T_1(T_2( z)) Concatenation of two functions T1T_1 and T2T_2 ‖2 \| x \|_2 2-norm of a vector x Cov(⋅)Cov(·) Covariance matrix A.2 Definitions μ-strongly convex function (Boyd and Vandenberghe, 2004). A differentiable function f:ℝm→ℝf:R^m is μ-strongly convex for a certain μ>0μ>0 if f()≥f()+∇f()T(−)+μ2‖−‖22f( y)≥ f( x)+∇ f( x)^T( y- x)+ μ2 \| y- x \|^2_2 (17) for all ,∈ℝm x, y ^m. The notion of of nonexpansive and averaged nonexpansive have been widely used in the convergence analysis of various nonlinear problems (Combettes and Yamada, 2015; Yao et al., 2008; Eckstein and Bertsekas, 1992). We use the generalized form of both nonexpansive and averaged nonexpansive operator for establishing the ball convergence in our method. Nonexpansive function (Browder, 1965). A function T:ℝd→ℝdT:R^d ^d is nonexpansive if T is nonexpansive function if there exists ϵ∈[0,1]ε∈[0,1] such that ‖T()−T()‖22≤ϵ2‖−‖22 \|T( x)-T( y) \|_2^2≤ε^2 \| x- y \|_2^2 (18) for all ,∈ℝd x, y ^d. θ-averaged function (Combettes and Yamada, 2015). A mapping T:ℝd→ℝdT:R^d ^d is defined to be θ-averaged for a constant θ∈(0,1)θ∈(0,1) if there exists a nonexpansive operator R:ℝd→ℝdR:R^d ^d such that T=(1−θ)I+θRT=(1-θ)I+θ R. The notion of relaxed bound ‖Tk()−Tk()‖≤ϵ(k)‖−‖+δ(k),∀,∈ \|T_k( x)-T_k( y) \|≤ε^(k) \| x- y \|+δ^(k),\;∀ x, y was used to study and show the convergence in Yao et al. (2008) when ∑k=1∞|δ(k)|<∞ _k=1^∞|δ^(k)|<∞. We define a similar weaker form of nonexpansive function and θ-averaged functions below. δ-weakly nonexpansive function. A mapping T:ℝd→ℝdT:R^d ^d is said to be δ-weakly nonexpansive for δ≥0δ≥ 0 if there exists ϵ∈[0,1]ε∈[0,1] such that ‖T()−T()‖22≤ϵ2‖−‖22+δ2 \|T( x)-T( y) \|_2^2≤ε^2 \| x- y \|_2^2+δ^2 (19) for all ,∈ℝd x, y ^d. δ-weakly θ-averaged function. A mapping T:ℝd→ℝdT:R^d ^d is defined to be δ-weakly θ-averaged function for a certain δ≥0δ≥ 0 and θ∈(0,1)θ∈(0,1), if there exists a δ-weakly nonexpansive function R:ℝd→ℝdR:R^d ^d such that T=θR+(1−θ)IT=θ R+(1-θ)I. Sub-Gaussian random vector. A random vector ∈ℝd x ^d (with mean []E[ x]) is called sub-Gaussian with parameter σ2σ^2 if its Euclidean norm satisfies a sub-Gaussian tail bound: Pr⁡(‖−[]‖2>ε)≤2exp⁡(−ε22σ2),∀ε>0. \! (\| x-E[ x]\|_2> )≤ 2 \! (- ^22σ^2 ), ∀\, >0. (20) 2-Wasserstein Distance. Let μ and ν be probability measures on ℝdR^d with finite second moments. The 2-Wasserstein distance between μ and ν is defined as W2(μ,ν)=(infγ∈Γ(μ,ν)∫ℝd×ℝd‖x−y‖22γ(x,y))1/2, W_2(μ,ν)= ( _γ∈ (μ,ν) _R^d×R^d\|x-y\|_2^2\,dγ(x,y) )^1/2, (21) where Γ(μ,ν) (μ,ν) denotes the set of all couplings of μ and ν, i.e., Γ(μ,ν)=γ∈(ℝd×ℝd):γ(A×ℝd)=μ(A),γ(ℝd×B)=ν(B),∀A,B⊆ℝd measurable. (μ,ν)= \γ (R^d×R^d):\;γ(A×R^d)=μ(A),\;γ(R^d× B)=ν(B),\;∀ A,B ^d measurable \. (22) A.3 Supporting Lemmas Tweedie’s lemma establishes an important connection between the score of the marginal distribution and expectation of posterior when the likelihood function is gaussian. This allows the score function of the diffusion model to be used as a minimum-mean-square-error (MMSE) denoiser. Lemma 1 (Tweedie’s lemma (Robbins, 1992)). Let p0(0)p_0( x_0) be the prior distribution and then t∼(0,) x_t ( x_0, ) be observed with known. Suppose pt(t)p_t( x_t) be the marginal distribution of t x_t. Then, Tweedie’s lemma computes the posterior expectation of 0 x_0 given t x_t as [0|t]=t+∇log⁡pt(t)E[ x_0| x_t]= x_t+ ∇ p_t( x_t) (23) The lemmas related to θ-averaged from Combettes and Yamada (2015) are used to show fixed point convergence in Ryu et al. (2019). In the following, we extend all these lemmas to a more general δ-weakly θ-averaged cases that will be used later to show our ball convergence. Lemma 2. T:ℝd→ℝdT:R^d ^d be a function. Then, the following statements are equivalent: (a) T is δ-weakly θ-averaged for δ≥0δ≥ 0 and θ∈(0,1)θ∈(0,1). (b) ‖T()−T()‖22+(1−2θ)‖−‖22−2(1−θ)⟨T()−T(),−⟩≤δ2θ2 \|T( x)-T( y) \|_2^2+(1-2θ) \| x- y \|_2^2-2(1-θ) T( x)-T( y), x- y ≤δ^2θ^2, for all ,∈ℝd x, y ^d. (c) (1−1/θ)I+(1/θ)T(1-1/θ)I+(1/θ)T is δ-weakly nonexpansive. (d) ‖T()−T()‖22≤‖−‖22−1−θ‖(I−T)()−(I−T)()‖22+δ2θ \|T( x)-T( y) \|_2^2≤ \| x- y \|_2^2- 1-θ \|(I-T)( x)-(I-T)( y) \|_2^2+δ^2θ, for all ,∈ℝd x, y ^d Proof. Equivalence between (a) and (b): Provided T is δ-weakly θ-averaged, let’s find the LHS - RHS in (b) (1−2θ)‖−‖22+‖T()−T()‖22−2(1−θ)⟨−,T()−T()⟩ (1-2θ) \| x- y \|_2^2+ \|T( x)-T( y) \|_2^2-2(1-θ) x- y,T( x)-T( y) ≤ ≤ (1−2θ)‖−‖22+‖θR()−θR()+(1−θ)(−)‖22 (1-2θ) \| x- y \|_2^2+ \|θ R( x)-θ R( y)+(1-θ)( x- y) \|_2^2 −2(1−θ)⟨−,θ(R()−R())+(1−θ)(−)⟩ 42.67912pt-2(1-θ) x- y,θ(R( x)-R( y))+(1-θ)( x- y) = = (1−2θ)‖−‖22+θ2‖R()−R()‖22+(1−θ)2‖−‖22−2(1−θ)2‖−‖22 (1-2θ) \| x- y \|_2^2+θ^2 \|R( x)-R( y) \|_2^2+(1-θ)^2 \| x- y \|_2^2-2(1-θ)^2 \| x- y \|_2^2 +2θ(1−θ)⟨−,R()−R()⟩−2(1−θ)θ⟨−,R()−R()⟩ 42.67912pt+2θ(1-θ) x- y,R( x)-R( y) -2(1-θ)θ x- y,R( x)-R( y) ≤ ≤ (1−2θ)‖−‖22+θ2‖−‖22+θ2δ2+(1−θ)2‖−‖22−2(1−θ)2‖−‖22 (1-2θ) \| x- y \|_2^2+θ^2 \| x- y \|_2^2+θ^2δ^2+(1-θ)^2 \| x- y \|_2^2-2(1-θ)^2 \| x- y \|_2^2 +2θ(1−θ)⟨−,R()−R()⟩−2(1−θ)θ⟨−,R()−R()⟩ 42.67912pt+2θ(1-θ) x- y,R( x)-R( y) -2(1-θ)θ x- y,R( x)-R( y) = = (1−2θ+θ2+(1−θ)2−2(1−θ)2)‖−‖22+θ2δ2 (1-2θ+θ^2+(1-θ)^2-2(1-θ)^2) \| x- y \|^2_2+θ^2δ^2 = = θ2δ2 θ^2δ^2 For another direction, let us suppose T satisfies (a). Let R=1θ(T−(1−θ)I)R= 1θ (T-(1-θ)I ) so that we have T=θR+(1−θ)IT=θ R+(1-θ)I. Now, we need to show that ‖R()−R()‖22≤‖−‖22+δ2 \|R( x)-R( y) \|_2^2≤ \| x- y \|_2^2+δ^2 for all ,∈ℝd x, y ^d i.e. δ-weakly nonexpansive. ‖R()−R()‖22 \|R( x)-R( y) \|_2^2 = = 1θ2‖T()−T()−(1−θ)(−)‖22 1θ^2 \|T( x)-T( y)-(1-θ)( x- y) \|_2^2 = = 1θ2(‖T()−T()‖22+(1−θ)2‖−‖22−2(1−θ)⟨T()−T(),−⟩) 1θ^2 ( \|T( x)-T( y) \|_2^2+(1-θ)^2 \| x- y \|_2^2-2(1-θ) T( x)-T( y), x- y ) = = 1θ2(‖T()−T()‖22+(1−2θ)‖−‖22−2(1−θ)⟨T()−T(),−⟩+θ2‖−‖22) 1θ^2 ( \|T( x)-T( y) \|_2^2+(1-2θ) \| x- y \|_2^2-2(1-θ) T( x)-T( y), x- y +θ^2 \| x- y \|_2^2 ) ≤ ≤ 1θ2(θ2δ2+θ2‖−‖22) 1θ^2 (θ^2δ^2+θ^2 \| x- y \|_2^2 ) = = ‖−‖22+δ2 \| x- y \|_2^2+δ^2 where, the inequality is due to T satisfying (b). Equivalence between (a) and (c): Note that T is δ-weakly θ-averaged⇔T=θR+(1−θ)IT is δ-weakly θ-averaged T=θ R+(1-θ)I with R being δ-weakly nonexpansive function. Now, we have (1−1/θ)I+(1/θ)T (1-1/θ)I+(1/θ)T =(1−1/θ)I+1/θ⋅(θR+(1−θ)I) =(1-1/θ)I+1/θ· (θ R+(1-θ)I ) =(1−1/θ)I+R−(1−1/θ)I =(1-1/θ)I+R-(1-1/θ)I =R =R Hence, T being δ-weakly θ-average is equivalent to (1−1/θ)I+(1/θ)T(1-1/θ)I+(1/θ)T being δ-weakly nonexpansive. Equivalence between (a) and (d): From equivalence between (a) and (b), we have ‖T()−T()‖22+(1−2θ)‖−‖22−2(1−θ)⟨T()−T(),−⟩≤δ2θ2 \|T( x)-T( y) \|_2^2+(1-2θ) \| x- y \|_2^2-2(1-θ) T( x)-T( y), x- y ≤δ^2θ^2 ⇔ ‖T()−T()‖22+(1−2θ)‖−‖22− \|T( x)-T( y) \|_2^2+(1-2θ) \| x- y \|_2^2- (1−θ)(‖T()−T()‖22+‖−‖22−‖(T−I)()−(T−I)()‖22)≤δ2θ2 42.67912pt(1-θ) ( \|T( x)-T( y) \|_2^2+ \| x- y \|_2^2- \|(T-I)( x)-(T-I)( y) \|_2^2 )≤δ^2θ^2 ⇔ θ‖T()−T()‖22−θ‖−‖22≤δ2θ2−(1−θ)‖(T−I)()−(T−I)()‖22 θ \|T( x)-T( y) \|_2^2-θ \| x- y \|_2^2≤δ^2θ^2-(1-θ) \|(T-I)( x)-(T-I)( y) \|_2^2 ⇔ ‖T()−T()‖22≤‖−‖22−1−θ‖(I−T)()−(I−T)()‖22+δ2θ \|T( x)-T( y) \|_2^2≤ \| x- y \|_2^2- 1-θ \|(I-T)( x)-(I-T)( y) \|_2^2+δ^2θ ∎ Lemma 3 (Concatenation of δ-weakly θ-averaged functions). Assume T1:ℝd→ℝdT_1:R^d ^d and T2:ℝd→ℝdT_2:R^d ^d are δ1 _1-weakly θ1 _1-averaged and δ2 _2-weakly θ2 _2-averaged respectively. Then, T1∘T2T_1 T_2 is δ-weakly θ-averaged, with θ=θ1+θ2−2θ1θ21−θ1θ2θ= _1+ _2-2 _1 _21- _1 _2, and δ2=1θ(δ12θ1+δ22θ2)δ^2= 1θ( _1^2 _1+ _2^2 _2). Proof. Here, we follow the proof structure of Combettes and Yamada (2015). Since θ1(1−θ2)≤(1−θ2) _1(1- _2)≤(1- _2), we have θ1+θ2≤1+θ1θ2 _1+ _2≤ 1+ _1 _2, and therefore, θ=θ1+θ2−2θ1θ21−θ1θ2∈(0,1)θ= _1+ _2-2 _1 _21- _1 _2∈(0,1), and let δ2=δ12θ1+δ22θ2θδ^2= _1^2 _1+ _2^2 _2θ Now, from Lemma (d), for i∈1,2i∈\1,2\, we have, ‖Ti()−Ti()‖22≤‖−‖22−1−θiθi‖(I−Ti)()−(I−Ti)()‖22+δi2θi \|T_i( x)-T_i( y) \|_2^2≤ \| x- y \|_2^2- 1- _i _i \|(I-T_i)( x)-(I-T_i)( y) \|_2^2+ _i^2 _i (24) Then, let us evaluate the composition function using this property. ‖T1∘T2()−T1∘T2()‖22 \|T_1 T_2( x)-T_1 T_2( y) \|_2^2 ≤ ≤ ∥T2()−T2()∥22−1−θ1θ1∥(I−T1)(T2())−(I−T1)(T2()∥22 \|T_2( x)-T_2( y) \|_2^2- 1- _1 _1 \|(I-T_1)(T_2( x))-(I-T_1)(T_2( y) \|_2^2 ≤ ≤ ‖−‖22−1−θ2θ2‖(I−T2)()−(I−T2)()‖22+δ22θ2 \| x- y \|_2^2- 1- _2 _2 \|(I-T_2)( x)-(I-T_2)( y) \|_2^2+ _2^2 _2 −1−θ1θ1∥(I−T1)(T2())−(I−T1)(T2()∥22+δ12θ1 42.67912pt- 1- _1 _1 \|(I-T_1)(T_2( x))-(I-T_1)(T_2( y) \|_2^2+ _1^2 _1 From Bauschke et al. (2017)[Corollary 2.15], we have, for α∈ℝα , ‖α+(1−α)‖22+α(1−α)‖−‖22=α‖22+(1−α)‖22 \|α u+(1-α) v \|_2^2+α(1-α) \| u- v \|_2^2=α \| u \|_2^2+(1-α) \| v \|_2^2 ⟹ α(1−α)‖+‖22≤α‖22+(1−α)‖22 α(1-α) \| u+ v \|_2^2≤α \| u \|_2^2+(1-α) \| v \|_2^2 Now, let =(I−T2)()−(I−T2)() u=(I-T_2)( x)-(I-T_2)( y), =(I−T1)(T2())−(I−T1)(T2()) v=(I-T_1)(T_2( x))-(I-T_1)(T_2( y)), a=1−θ2θ2a= 1- _2 _2, and b=1−θ1θ1b= 1- _1 _1. ‖T1∘T2()−T1∘T2()‖22 \|T_1 T_2( x)-T_1 T_2( y) \|_2^2 ≤ ≤ ‖−‖22−a‖22−b‖22+δ12θ1+δ22θ2 \| x- y \|_2^2-a \| u \|_2^2-b \| v \|_2^2+ _1^2 _1+ _2^2 _2 = = ‖−‖22−(a+b)(a+b‖22+ba+b‖22)+δ12θ1+δ22θ2 \| x- y \|_2^2-(a+b) ( aa+b \| u \|_2^2+ ba+b \| v \|_2^2 )+ _1^2 _1+ _2^2 _2 = = ‖−‖22−(a+b)(a+b‖22+(1−a+b)‖22)+δ12θ1+δ22θ2 \| x- y \|_2^2-(a+b) ( aa+b \| u \|_2^2+ (1- aa+b ) \| v \|_2^2 )+ _1^2 _1+ _2^2 _2 Using the above results, we get, ‖T1∘T2()−T1∘T2()‖22 \|T_1 T_2( x)-T_1 T_2( y) \|_2^2 ≤ ≤ ‖−‖22−(a+b)ab(a+b)2‖+‖22+δ12θ1+δ22θ2 \| x- y \|_2^2-(a+b) ab(a+b)^2 \| u+ v \|_2^2+ _1^2 _1+ _2^2 _2 = = ‖−‖22−ab(a+b)‖(I−T1∘T2)()−(I−T1∘T2)()‖22+δ12θ1+δ22θ2 \| x- y \|_2^2- ab(a+b) \|(I-T_1 T_2)( x)-(I-T_1 T_2)( y) \|_2^2+ _1^2 _1+ _2^2 _2 Let θ=θ1+θ2−2θ1θ21−θ1θ2θ= _1+ _2-2 _1 _21- _1 _2. Then, we can see that aba+b=1−θ aba+b= 1-θ. ‖T1∘T2()−T1∘T2()‖22 \|T_1 T_2( x)-T_1 T_2( y) \|_2^2 ≤‖−‖22−1−θ‖(I−T1∘T2)()−(I−T1∘T2)()‖22 ≤ \| x- y \|_2^2- 1-θ \|(I-T_1 T_2)( x)-(I-T_1 T_2)( y) \|_2^2 +δ2θ 42.67912pt+δ^2θ (25) where, δ2θ=δ12θ1+δ22θ2δ^2θ= _1^2 _1+ _2^2 _2. This implies that T1∘T2T_1 T_2 is δ-weakly θ-averaged with θ=θ1+θ2−2θ1θ21−θ1θ2,δ2 θ= _1+ _2-2 _1 _21- _1 _2,\;δ^2 =1θ(δ12θ1+δ22θ2) = 1θ( _1^2 _1+ _2^2 _2) (26) ∎ Lemma 4 (Proposition 5.4 of Giselsson (2015)). Assume ℓ is μ-strongly convex, closed, and proper. Then, −(2Prox1ρℓ−I)-(2 Prox_ 1ρ -I) is ρ+μ ρ+μ-averaged. Lemma 5 (Pardo (2018)). The KL divergence between two gaussian distributions q1=(1,1)q_1=N( μ_1, _1) and q2=(2,2)q_2=N( μ_2, _2) in ℝdR^d space is given by KL(q1||q2)=12[log|2||1|−d+tr(2−11)+(2−1)T2−1(2−1)]KL(q_1||q_2)= 12 [ | _2|| _1|-d+tr( _2^-1 _1)+( μ_2- μ_1)^T _2^-1( μ_2- μ_1) ] (27) where |⋅||·| denotes the determinant, and tr denotes the trace of the matrix. Appendix B Influence of AC-step on bringing close to ℳσ(t)t=0T\M_σ(t)\_t=0^T Given a noisy image ~(k) z^(k) at each iteration k, the denoising aims to recover the underlying clean image ♮(k)∼p0() z^(k)_ p_0( z) such that ~(k)=♮(k)+(k) z^(k)= z^(k)_ + s^(k), where (k) s^(k) is the noise contained in ~(k) z^(k). The AC-step aims to bring ac(k) z^(k)_ ac closer to the noisy distribution ℳσ(k)M_σ^(k) on which the (⋅,σ(k)) s_ θ(·,σ^(k)) was trained on. Lemma 6 shows that the AC-step tries to match with the distribution induced by the forward diffusion process. Lemma 6. The KL divergence between the target distribution p(σ(k)|♮(k))p( z_σ^(k)| z^(k)_ ) for correction steps and the distribution p(ac(k)|~(k))p( z^(k)_ ac| z^(k)) induced by the approximate correction step in Algorithm 1 is given by KL(p(σ(k)|♮(k))||p(ac(k)|~(k)))=12(σ(k))2∥(k)∥22KL(p( z_σ^(k)| z^(k)_ )||p( z^(k)_ ac| z^(k)))= 12 (σ^(k) )^2 \| s^(k) \|_2^2 (28) Proof. Target distribution: p(σ(k)|♮(k))=(♮(k),(σ(k))2) distribution: p ( z_σ^(k)| z^(k)_ )=N ( z^(k)_ , (σ^(k) )^2 I ) (29) AC induced distribution: p(ac(k)|~(k))=(~(k),(σ(k))2) induced distribution: p ( z^(k)_ ac| z^(k) )=N ( z^(k), (σ^(k) )^2 I ) (30) The KL divergence between these two distribution can be computed in closed form using Lemma 5. KL(q1||q2) (q_1||q_2) =12(log⁡1−d+tr()+(~(k)−♮(k))T(σ(k))−2(~(k)−♮(k))) = 12 ( 1-d+tr( I)+( z^(k)- z^(k)_ )^T (σ^(k) )^-2 I( z^(k)- z^(k)_ ) ) =12(0−d+d+(σ(k))−2‖~(k)−♮(k)‖22) = 12 (0-d+d+ (σ^(k) )^-2 \| z^(k)- z^(k)_ \|_2^2 ) =12(σ(k))2‖(k)‖22 = 12 (σ^(k) )^2 \| s^(k) \|_2^2 where, (k)=~(k)−♮(k) s^(k)= z^(k)- z^(k)_ . ∎ Lemma 6 shows that KL-gap of our approximate AC update; as long as σ(k)σ^(k) is sufficiently large, the two distributions remain close. Alkhouri et al. (2023)[Theorem 1] showed a result with a similar flavor. Larger noise σ(k)σ^(k) makes the posterior nearly indistinguishable, but it also washes out fine structural details originally present (low Signal-to-Noise Ratio with larger σ(k)σ^(k)). Existing works often use annealed scheduling σ(k)↓0σ^(k) 0 (Zhu et al., 2023; Renaud et al., 2024a; Wang et al., 2024) to preserve image details, implicitly assuming ‖(k)‖22 \| s^(k) \|_2^2 decays at least as fast as (σ(k))2(σ^(k))^2. With just the use of annealing σ(k)σ^(k) schedule, it is not sufficient to ensure that ac(k) z^(k)_ ac lands in a desired manifold in each ADMM iteration. To bridge this gap, we propose to use DC-step in addition to the widely used annealing σ(k)σ^(k) schedule that explicitly corrects this gap. Appendix C Proof of Theorem 1 The proof involves showing that the each iteration of ADMM-PnP is weakly non-expansive when the denoiser satisfies Assumption 1. This weakly nonexpansiveness of each step leads to ball convergence of the algorithm. Recall that the subproblems at kkth iteration of ADMM-PnP is given by: (k+1) x^(k+1) =argmin1ρℓ(||())+12∥−(k)+(k)∥22 = _ x 1ρ ( y||A( x))+ 12 \| x- z^(k)+ u^(k) \|_2^2 =Prox1ρℓ((k)−(k)) =Prox_ 1ρ ( z^(k)- u^(k)) (31a) (k+1) z^(k+1) =arg⁡min⁡γρh()+12‖(k+1)−+(k)‖22 = _ z γρh( z)+ 12 \| x^(k+1)- z+ u^(k) \|_2^2 =Proxγρh((k+1)+(k)) =Prox_ γρh( x^(k+1)+ u^(k)) =Dσ(k)((k+1)+(k)) =D_σ^(k)( x^(k+1)+ u^(k)) (31b) (k+1) u^(k+1) =(k)+((k+1)−(k+1)) = u^(k)+( x^(k+1)- z^(k+1)) (31c) Lemma 7 (Ryu et al. (2019)). The steps of ADMM-PnP in (31) can be expressed as (k+1)=T((k)) v^(k+1)=T( v^(k)) with (k)=(k)−(k) v^(k)= z^(k)- u^(k) and T=12I+12(2Dσ(k)−I)(2Prox1ρℓ−I)T= 12I+ 12(2D_σ^(k)-I)(2 Prox_ 1ρ -I) (32) Lemma 8. Dσ:ℝd→ℝdD_σ:R^d ^d satisfies Assumption 1 if and only if 11+2ϵ(2Dσ(k)−I) 11+2ε(2D_σ^(k)-I) is Δ -weakly θ-averaged with θ=2ϵ1+2ϵθ= 2ε1+2ε and Δ2=4δ2(1−θ)2θ2 ^2=4δ^2 (1-θ)^2θ^2. Proof. We follow the similar proof structure as in Ryu et al. (2019). Let θ=2ϵ1+2ϵθ= 2ε1+2ε which implies ϵ=θ2(1−θ)ε= θ2(1-θ). Here, we can clearly see that θ∈[0,1)θ∈[0,1). Let us define G=11+2ϵ(2Dσ(k)−I)G= 11+2ε(2D_σ^(k)-I) which implies Dσ(k)=12(1−θ)G+12ID_σ^(k)= 12(1-θ)G+ 12I. Then, ‖(Dσ(k)−I)()−(Dσ(k)−I)()‖2−ϵ2‖−‖22 \|(D_σ^(k)-I)( x)-(D_σ^(k)-I)( y) \|^2-ε^2 \| x- y \|_2^2 = = ‖(Dσ(k)()−Dσ(k)())−(−)‖22−ϵ2‖−‖22 \| (D_σ^(k)( x)-D_σ^(k)( y) )- ( x- y ) \|_2^2-ε^2 \| x- y \|_2^2 = = ‖(Dσ(k)−12I)()−(Dσ(k)−12I)()‖22+14‖−‖22−θ24(1−θ)2‖−‖22 \| (D_σ^(k)- 12I )( x)- (D_σ^(k)- 12I )( y) \|_2^2+ 14 \| x- y \|_2^2- θ^24(1-θ)^2 \| x- y \|_2^2 −2⟨(Dσ(k)−12I)()−(Dσ(k)−12)(),12(−)⟩ 42.67912pt-2 (D_σ^(k)- 12I )( x)- (D_σ^(k)- 12 )( y), 12( x- y) = = 14(1−θ)2‖G()−G()‖22+14(1−θ2(1−θ)2)‖−‖22 14(1-θ)^2 \|G( x)-G( y) \|_2^2+ 14 (1- θ^2(1-θ)^2 ) \| x- y \|_2^2 −12(1−θ)⟨G()−G(),−⟩ 42.67912pt- 12(1-θ) G( x)-G( y), x- y = = 14(1−θ)2(‖G()−G()‖22−2(1−θ)⟨G()−G(),−⟩+(1−2θ)‖−‖22) 14(1-θ)^2 ( \|G( x)-G( y) \|_2^2-2(1-θ) G( x)-G( y), x- y +(1-2θ) \| x- y \|_2^2 ) Now, 14(1−θ)2(‖G()−G()‖22−2(1−θ)⟨G()−G(),−⟩+(1−2θ)‖−‖22)≤δ2 14(1-θ)^2 ( \|G( x)-G( y) \|_2^2-2(1-θ) G( x)-G( y), x- y +(1-2θ) \| x- y \|_2^2 )≤δ^2 ⇔ ‖G()−G()‖22−2(1−θ)⟨G()−G(),−⟩+(1−2θ)‖−‖22≤4δ2(1−θ)2 \|G( x)-G( y) \|_2^2-2(1-θ) G( x)-G( y), x- y +(1-2θ) \| x- y \|_2^2≤ 4δ^2(1-θ)^2 ⇔ ‖G()−G()‖22−2(1−θ)⟨G()−G(),−⟩+(1−2θ)‖−‖22≤Δ2θ2 \|G( x)-G( y) \|_2^2-2(1-θ) G( x)-G( y), x- y +(1-2θ) \| x- y \|_2^2≤ ^2θ^2 where, Δ2=4δ2(1−θ)2θ2 ^2=4δ^2 (1-θ)^2θ^2. From Lemma (b), this is equivalent to G being Δ -weakly θ-averaged. ∎ C.1 Proof of the theorem We follow the procedures in Ryu et al. (2019) and expand the results in the δ-weakly expansive denoisers. We show that each iteration of PnP ADMM is also weakly nonexpansive when the denoiser satisfies Assumption 1. Proof. From Assumption 1. ‖(Dσ(k)−I)()−(Dσ(k)−I)()‖22≤ϵ2‖−‖22+δ2 \|(D_σ^(k)-I)( x)-(D_σ^(k)-I)( y) \|_2^2≤ε^2 \| x- y \|_2^2+δ^2 (33) From Lemma 4, we have −(2Prox1ρℓ−I)-(2Prox_ 1ρ -I) is ρ+μ ρ+μ-averaged. Then, from Lemma 8, we have 11+2ϵ(2Dσ(k)−I) 11+2ε (2D_σ^(k)-I ) (34) is δ1 _1-weakly θ-averaged with θ=2ϵ1+2ϵθ= 2ε1+2ε and δ12=4δ2(1−θ)2θ2 _1^2=4δ^2 (1-θ)^2θ^2. By Lemma 3, it implies −11+2ϵ(2Dσ(k)−I)(2Prox1ρℓ−I)- 11+2ε(2D_σ^(k)-I)(2Prox_ 1ρ -I) (35) is δ∘ _ -weakly θ∘ _ -averaged with θ∘=ρ+2μϵρ+μ+2μϵ _ = ρ+2μερ+μ+2με and δ∘2=1θ∘⋅4δ22ϵ(1+2ϵ) _ ^2= 1 _ · 4δ^22ε(1+2ε). Now, using the definition of δ∘ _ -weakly θ∘ _ -averagedness, we have (2Dσ(k)−I)(2Prox1ρℓ−I) (2D_σ^(k)-I)(2Prox_ 1ρ -I) =−(1+2ϵ)((1−θ∘)I+θ∘R) =-(1+2ε) ((1- _ )I+ _ R ) =−(1+2ϵ)(μρ+μ+2μϵI+ρ+2μϵρ+μ+2μϵR) =-(1+2ε) ( μρ+μ+2μεI+ ρ+2μερ+μ+2μεR ) where, R is a certain δ∘ _ -weakly nonexpansive function. Plugging this result into ADMM-PnP operator (Lemma 7), we get T T =12I+12(2Dσ(k)−I)(2Prox1ρℓ−I) = 12I+ 12(2D_σ^(k)-I)(2Prox_ 1ρ -I) =12I−12(1+2ϵ)(μρ+μ+2μϵI+ρ+2μϵρ+μ+2μϵR) = 12I- 12(1+2ε) ( μρ+μ+2μεI+ ρ+2μερ+μ+2μεR ) =ρ2(ρ+μ+2μϵ)⏟aI−(1+2ϵ)(ρ+2μϵ)2(ρ+μ+2μϵ)⏟bR = ρ2(ρ+μ+2με)_aI- (1+2ε)(ρ+2με)2(ρ+μ+2με)_bR where, clearly a>0a>0 and b>0b>0. Now, ‖T()−T()‖22=a2‖−‖22+b2‖R()−R()‖22−2⟨a(−),b(R()−R())⟩ \|T( x)-T( y) \|_2^2=a^2 \| x- y \|_2^2+b^2 \|R( x)-R( y) \|_2^2-2 a( x- y),b(R( x)-R( y)) (36) From Young’s inequality, for any γ>0γ>0, we have ⟨a(−),b(R()−R())⟩≤12γa2‖−‖22+γb22‖R()−R()‖22 a( x- y),b(R( x)-R( y)) ≤ 12γa^2 \| x- y \|_2^2+ γ b^22 \|R( x)-R( y) \|_2^2 (37) Plugging this, we get, ‖T()−T()‖22 \|T( x)-T( y) \|_2^2 ≤a2(1+1γ)‖−‖22+b2(1+γ)‖R()−R()‖22 ≤ a^2 (1+ 1γ ) \| x- y \|_2^2+b^2(1+γ) \|R( x)-R( y) \|_2^2 (38) ≤(a2(1+1γ)+b2(1+γ))‖−‖22+b2(1+γ)δ∘2 ≤ (a^2 (1+ 1γ )+b^2(1+γ) ) \| x- y \|_2^2+b^2(1+γ) _ ^2 (39) where, the second inequality is due to δ∘ _ -weak nonexpansiveness of R. Note, that this holds for any γ>0γ>0. When γ=abγ= ab, we have (a2(1+1γ)+b2(1+γ))=(a+b)2 (a^2 (1+ 1γ )+b^2(1+γ) )=(a+b)^2 (40) ‖T()−T()‖22 \|T( x)-T( y) \|_2^2 ≤ ≤ (a+b)2‖−‖22+b2(1+ab)δ∘2 (a+b)^2 \| x- y \|_2^2+b^2 (1+ ab ) _ ^2 = = (ρ+ρϵ+μϵ+2μϵ2ρ+μ+2μϵ)2⏟ϵT2‖−‖22+(ρ+ρϵ+μϵ+2μϵ2)δ2ϵ(ρ+μ+2μϵ)⏟δT2 ( ρ+ρε+με+2με^2ρ+μ+2με )^2_ _T^2 \| x- y \|_2^2+ (ρ+ρε+με+2με^2)δ^2ε(ρ+μ+2με)_ _T^2 (41) Hence, we have ‖T()−T()‖22≤ϵT2‖−‖22+δT2 \|T( x)-T( y) \|_2^2≤ _T^2 \| x- y \|_2^2+ _T^2 (42) This shows that when ϵT≤1 _T≤ 1, then, ∃∗∈ℝd _ v^* ^d, and K>0K>0 such that ∀k≥K∀ k≥ K the following holds: ‖(k)−∗‖22≤ϵT2k‖(0)−∗‖22+δT21−ϵT2 \| v^(k)- v^* \|_2^2≤ _T^2k \| v^(0)- v^* \|_2^2+ _T^21- _T^2 ⟹ limk→∞‖(k)−∗‖≤δT1−ϵT2 _k→∞ \| v^(k)- v^* \|≤ _T 1- _T^2 (43) Hence with this we have the sequence (k)=(k)−(k)k∈ℕ+\ v^(k)= z^(k)- u^(k)\_k ^+ converges within a ball of radius δT1−ϵT2 _T 1- _T^2. Since, −(2Prox1ρℓ−I)-(2Prox_ 1ρ -I) is ρ+μ ρ+μ-averaged, this implies Prox1ρℓ=12ρ+μ(I−R) Prox_ 1ρ = 12 ρ+μ(I-R) (44) for some nonexpansive function R. With this, we have limk→∞‖Prox1ρℓ((k))−Prox1ρℓ(∗)‖22≤(ρ+μ)2‖(k)−∗‖22 _k→∞ \| Prox_ 1ρ ( v^(k))- Prox_ 1ρ ( v^*) \|_2^2≤ ( ρ+μ )^2 \| v^(k)- v^* \|_2^2 ⟹ limk→∞‖(k)−∗‖2≤ρ+μδT1−ϵT2 _k→∞ \| x^(k)- x^* \|_2≤ ρ+μ _T 1- _T^2 (45) where ∗=Prox1ρℓ(∗) x^*= Prox_ 1ρ ( v^*). With these results, we know there exists ∗ u^* such that limk→∞‖k−∗‖2≤(1+ρ+μ)δT1−ϵT2 _k→∞ \| u^k- u^* \|_2≤ (1+ ρ+μ ) _T 1- _T^2 (46) ∎ Appendix D Proof of Theorem 2 Here, we show that our 3-step AC-DC denoiser satisfies Assumption 1 for constants ϵε and δ. In the following, we first show that each step satisfies the weakly nonexpansive assumption. Therefore, the concatenation of these 3 steps meets Assumption 1. Lemma 9. Assume the Variance Exploding (VE) scheduling (Karras et al., 2022) is used in the diffusion model. Given that the log-density log⁡p0 p_0 (i.e. log⁡pdata p_ data) is M-smooth (Assumption 2), the intermediate noisy log-densities log⁡pt\ p_t\ are MtM_t-smooth for t∈[0,T]t∈[0,T] i.e. ‖∇log⁡pt()−∇log⁡pt()‖2≤Mt‖−‖2 \|∇ p_t( x)-∇ p_t( y) \|_2≤ M_t \| x- y \|_2 for all ,∈ x, y . For t such that σ2(t)<1/Mσ^2(t)<1/M, the smoothness constant MtM_t can be upper bounded as Mt≤M1+Mσ2(t)≤M_t≤ M1+Mσ^2(t)≤ M (47) Proof. From Tweedie’s lemma, we have, ⟹ ∇tlog⁡pt(t)=−1σ2(t)(t−[0|t]) _ x_t p_t( x_t)=- 1σ^2(t) ( x_t-E[ x_0| x_t] ) ⟹ ∇t2log⁡pt(t)=−1σ2(t)(−∇t[0|t]) _ x_t^2 p_t( x_t)=- 1σ^2(t)( I- _ x_tE[ x_0| x_t]) (48) Now, let us evaluate the Jacobian ∇t[0|t] _ x_tE[ x_0| x_t], ∇t[0|t] _ x_tE[ x_0| x_t] =∇t∫0∈0p(0|t)0 = _ x_t _ x_0 x_0p( x_0| x_t)d x_0 =∫0∈0(∇tp(t|0)pt(t))p0(0)0 = _ x_0 x_0 ( _ x_t p( x_t| x_0)p_t( x_t) )p_0( x_0)d x_0 =∫0∈0(1pt(t)∇tp(t|0)−p(t|0)1pt2(t)∇tpt(t))p0(0)0 = _ x_0 x_0 ( 1p_t( x_t) _ x_tp( x_t| x_0)-p( x_t| x_0) 1p_t^2( x_t) _ x_tp_t( x_t) )p_0( x_0)d x_0 (49) Given t|0∼(0,σ(t)2) x_t| x_0 ( x_0,σ(t)^2 I) and p(t)=∫0∈p(0)p(t|0)0p( x_t)= _ x_0 p( x_0)p( x_t| x_0)d x_0, we can compute their gradient (similar to Peng et al. (2024)) as: ∇tp(t|0) _ x_tp( x_t| x_0) =−1σ2(t)p(t|0)(t−0) =- 1σ^2(t)p( x_t| x_0)( x_t- x_0) ∇tp(t) _ x_tp( x_t) =∫0∈p0(0)∇tp(t|0)0 = _ x_0 p_0( x_0) _ x_tp( x_t| x_0)d x_0 =∫0∈p0(0)(−1σ2(t)p(t|0)(t−0))0 = _ x_0 p_0( x_0) (- 1σ^2(t)p( x_t| x_0)( x_t- x_0) )d x_0 =−1σ2(t)pt(t)[t−0|t] =- 1σ^2(t)p_t( x_t)E[ x_t- x_0| x_t] (50) Plugging these results in (49) and using integration by parts, ∇t[0|t]= _ x_tE[ x_0| x_t]= 1p(t)∫0∈0∇tp(t|0)p(0)0−∇tp(t)p(t)[0|t] 1p( x_t) _ x_0 x_0 _ x_tp( x_t| x_0)p( x_0)d x_0- _ x_tp( x_t)p( x_t)E[ x_0| x_t] (51) Substituting ∇tp(t|0)=−1/(σ2(t))(t−0)p(t|0) _ x_tp( x_t| x_0)=- 1(σ^2(t))( x_t- x_0)p( x_t| x_0), we have ∇t[0|t]= _ x_tE[ x_0| x_t]= −1σ2(t)p(t)∫0∈0(t−0)p(t|0)p(0)0−∇tp(t)p(t)[0|t] - 1σ^2(t)p( x_t) _ x_0 x_0( x_t- x_0)p( x_t| x_0)p( x_0)d x_0- _ x_tp( x_t)p( x_t)E[ x_0| x_t] =−1σ2(t)[0(t−0)|t]−∇tp(t)p(t)[0|t] = -1σ^2(t)E[ x_0( x_t- x_0)| x_t]- _ x_tp( x_t)p( x_t)E[ x_0| x_t] (52) Substituting ∇tp(t)=∫0∈p(0)∇tp(t|0)0=−pt(t)/σ2(t)[t−0|t] _ x_tp( x_t)= _ x_0 p( x_0) _ x_tp( x_t| x_0)d x_0=- p_t( x_t)σ^2(t)E[ x_t- x_0| x_t], we have ∇t[0|t] _ x_tE[ x_0| x_t] =−1σ2(t)[0(t−0)|t]+1σ2(t)[0|t][t−0t] = -1σ^2(t)E[ x_0( x_t- x_0)| x_t]+ 1σ^2(t)E[ x_0| x_t]E[ x_t- x_0 x_t] =−1σ2(t)Cov(0,t−0) =- 1σ^2(t)Cov( x_0, x_t- x_0) =1σ2(t)Cov(0|t) = 1σ^2(t)Cov( x_0| x_t) (53) Note that this result in (53) is similar to Dytso et al. (2021)[Proposition 1], and has also been derived in Hatsell and Nolte (1971); Palomar and Verdu (2006). From Assumption 2, we have −M⪯∇02log⁡p0(0)⪯M,∀0∈-M I _ x_0^2 p_0( x_0) M I,\;∀ x_0 (54) where, M>0M>0 is a constant. Then, let’s analyze the hessian of the log of posterior distribution p(0|t)p( x_0| x_t), ∇02log⁡p(0|t)=∇02log⁡p0(0)+∇02log⁡p(t|0) _ x_0^2 p( x_0| x_t)= _ x_0^2 p_0( x_0)+ _ x_0^2 p( x_t| x_0) (55) ⟹ −M+1σ2(t)⪯∇02log⁡p(0|t)⪯M+1σ2(t) -M I+ 1σ^2(t) I _ x_0^2 p( x_0| x_t) M I+ 1σ^2(t) I (56) When M<1σ2(t)M< 1σ^2(t), then the distribution log⁡p(0|t) p( x_0| x_t) is strongly log-concave. In this case, the covariance of distribution p(0|t)p( x_0| x_t) can be bounded (Brascamp and Lieb, 1976) as (M+1σ2(t))−1⪯Cov(0|t)⪯(−M+1σ2(t))−1 (M I+ 1σ^2(t) I )^-1 ( x_0| x_t) (-M I+ 1σ^2(t) I )^-1 (57) Combining result from equations (48), (53) and (57), we get ∇t2log⁡pt(t) _ x_t^2 p_t( x_t) =−1σ2(t)(−∇t[0|t]) =- 1σ^2(t) ( I- _ x_tE[ x_0| x_t] ) =−1σ2(t)(−1σ2(t)Cov(0|t)) =- 1σ^2(t) ( I- 1σ^2(t)Cov( x_0| x_t) I ) ‖∇t2log⁡qt(t)‖2 \| _ x_t^2 q_t( x_t) \|_2 ≤1σ2(t)|(1−1σ2(t)⋅σ2(t)Mσ2(t)+1)| ≤ 1σ^2(t) | (1- 1σ^2(t)· σ^2(t)Mσ^2(t)+1 ) | =1σ2(t)⋅Mσ2(t)1+M(σ2(t)) = 1σ^2(t)· Mσ^2(t)1+M(σ^2(t)) =M1+Mσ2(t) = M1+Mσ^2(t) (58) Hence, the smoothness constant of log⁡q(t) q( x_t) is upper bounded as Mt≤M1+Mσ2(t)M_t≤ M1+Mσ^2(t) i.e. Mt≤M_t≤ M. ∎ Lemma 10. Let Hac(k):~(k)↦ac(k)H^(k)_ ac: z^(k) z^(k)_ ac denote the function corresponding to approximate correction in Algorithm 1. Then, with probability at least 1−e−νk1-e^- _k, the following holds for any ,∈ x, y ‖(Hac(k)−I)()−(Hac(k)−I)()‖22≤(δac(k))2+(ϵac(k))2‖−‖22 \|(H^(k)_ ac-I)( x)-(H^(k)_ ac-I)( y) \|_2^2≤(δ^(k)_ ac)^2+(ε^(k)_ ac)^2 \| x- y \|_2^2 (59) where, δac(k)2=2(σ(k))2(d+2dνk+2νk)δ^(k)_ ac^2=2(σ^(k))^2(d+2 d _k+2 _k), and (ϵac(k))2=0(ε^(k)_ ac)^2=0. Proof. For any ,∈ x, y , we have the residuals Rac(k)()=(Hac(k)−I)()=σ(k)1R^(k)_ ac( x)=(H^(k)_ ac-I)( x)=σ^(k) n_1 and Rac(k)()=(Hac(k)−I)()=σ(k)2R^(k)_ ac( y)=(H^(k)_ ac-I)( y)=σ^(k) n_2 where, 1,2∼i.i.d.(,) n_1, n_2 i.i.d. N( 0, I). Then, we can bound the norm of difference of these two residuals as ‖Rac(k)()−Rac(k)()‖22 \|R^(k)_ ac( x)-R^(k)_ ac( y) \|_2^2 =‖σ(k)1−σ(k)2‖22 = \|σ^(k) n_1-σ^(k) n_2 \|_2^2 =(σ(k))2‖12‖22 =(σ^(k))^2 \| n_12 \|_2^2 =2(σ(k))2χd2 =2(σ^(k))^2 _d^2 (60) where, 12=1−2∼(,2) n_12= n_1- n_2 ( 0,2 I) and χd2 _d^2 is standard chi-square distribution with d degree of freedom. From Laurent and Massart (2000)[Lemma 1], the following holds with probability at least 1−e−νk1-e^- _k χd2≤d+2dνk+2νk _d^2≤ d+2 d _k+2 _k (61) Plugging this in proves the lemma. ∎ Lemma 11. Let Hdc(k):ac(k)↦dc(k)H^(k)_ dc: z^(k)_ ac z^(k)_ dc denote the function corresponding to fine correction as defined in Algorithm 1. Then, with probability at least 1−e−νk1-e^- _k, the following holds for any ,∈ x, y if σ(k)2<1Mt _ s^(k)^2< 1M_t: ‖(Hdc(k)−I)()−(Hdc(k)−I)()‖22≤(δdc(k))2+(ϵdc(k))2‖−‖22 \|(H^(k)_ dc-I)( x)-(H^(k)_ dc-I)( y) \|_2^2≤(δ^(k)_ dc)^2+(ε^(k)_ dc)^2 \| x- y \|_2^2 (62) where, (δdc(k))2=32dσ(k)2(1−Mtσ(k)2)log⁡2νk(δ^(k)_ dc)^2= 32d _ s^(k)^2(1-M_t _ s^(k)^2) 2 _k, and (ϵdc(k))2=(2Mtσ(k)21−σ(k)2Mt)2(ε^(k)_ dc)^2= ( 2M_t _ s^(k)^21- _ s^(k)^2M_t )^2. Proof. Recall that the target distribution for this step is given by log⁡p(σ(k)|ac(k))∝log⁡p(ac(k)|σ(k))+log⁡p(σ(k)) p( z_σ^(k)| z^(k)_ ac) p( z^(k)_ ac| z_σ^(k))+ p( z_σ^(k)) (63) where, p(ac(k)|σ(k))=(σ(k),σ(k)2)p( z^(k)_ ac| z_σ^(k))=N( z_σ^(k), _ s^(k)^2 I). Under Assumptions 2 and 3, the target distribution p(σ(k)|ac(k))p( z_σ^(k)| z^(k)_ ac) also inherits smoothness and coercivity properties. These conditions imply the ergodicity of corresponding Langevin diffusion (Mattingly et al., 2002; Chen et al., 2020). In particular, Fokker-Planck equation (Uhlenbeck and Ornstein, 1930) characterizes p(σ(k)|ac(k))p( z_σ^(k)| z^(k)_ ac) as its unique stationary distribution. Consequently, the iterates dc(k) z^(k)_ dc obtained through Langevin dynamics converge to this distribution as the step size η(k)→0η^(k)→ 0 and the number of iterations J→∞J→∞. The gradient and hessian of the log of this desired distribution are given by ∇dc(k)log⁡p(dc(k)|ac(k)) _ z^(k)_ dc p( z^(k)_ dc| z^(k)_ ac) =1σ(k)2(ac(k)−dc(k))+∇dc(k)log⁡pt(dc(k)) = 1 _ s^(k)^2( z^(k)_ ac- z^(k)_ dc)+ _ z^(k)_ dc p_t( z^(k)_ dc) (64) ∇dc(k)2log⁡p(dc(k)|ac(k)) _ z^(k)_ dc^2 p( z^(k)_ dc| z^(k)_ ac) =−1σ(k)2+∇dc(k)2log⁡pt(dc(k)) =- 1 _ s^(k)^2+ _ z^(k)_ dc^2 p_t( z^(k)_ dc) (65) Here, t refers to the noise level such that σ(t)=σ(k)σ(t)=σ^(k). By the MtM_t-smoothness of log⁡pt p_t distribution (Lemma 9), we have −(Mt+1σ(k)2)⪯∇dc(k)2log⁡p(dc(k)|ac(k))⪯(Mt−1σ(k)2)- (M_t+ 1 _ s^(k)^2 ) I _ z^(k)_ dc^2 p( z^(k)_ dc| z^(k)_ ac) (M_t- 1 _ s^(k)^2 ) I (66) When Mt<1σ(k)2M_t< 1 _ s^(k)^2, the hessian is negative semi-definite that implies the distribution being log concave. This also implies that when Mt≪1σ(k)2M_t 1 _ s^(k)^2, the likelihood term dominates in the posterior (63). Using (Wainwright, 2019)[Theorem 3.16], the following holds with probability at least 1−νk1- _k ∥dc(k)−[dc(k)|ac(k)]∥2≤4λtlog⁡2νk \| z^(k)_ dc-E[ z^(k)_ dc| z^(k)_ ac] \|_2≤ 4 _t 2 _k (67) where λt=−Mt+1σ(k)2 _t=-M_t+ 1 _ s^(k)^2. Now, using Tweedie’s lemma, we have [dc(k)|ac(k)]=ac(k)+σ(k)2∇ac(k)log⁡p(ac(k)) [ z^(k)_ dc| z^(k)_ ac]= z^(k)_ ac+ _ s^(k)^2 _ z^(k)_ ac p( z^(k)_ ac) ⟹ [dc(k)|ac(k)]−ac(k)=σ(k)2∇ac(k)log⁡p(ac(k)) [ z^(k)_ dc| z^(k)_ ac]- z^(k)_ ac= _ s^(k)^2 _ z^(k)_ ac p( z^(k)_ ac) (68) Combining the above two results, with probability at least 1−2νk1-2 _k, the difference of residual for x and y can be bounded as ‖Rdc(k)()−Rdc(k)()‖22 \|R^(k)_ dc( x)-R^(k)_ dc( y) \|_2^2 ≤2(24λtlog⁡2νk)22+‖σ(k)2(∇log⁡pac(k)()−∇log⁡pac(k)())‖22 ≤ 2 (2 4 _t 2 _k )_2^2+ \| _ s^(k)^2 ( _ x p_ z^(k)_ ac( x)- _ y p_ z^(k)_ ac( y) ) \|_2^2 ≤32dλtlog⁡2νk+2‖σ(k)2(∇log⁡pac(k)()−∇log⁡pac(k)())‖2 ≤ 32d _t 2 _k+2 \| _ s^(k)^2 ( _ x p_ z^(k)_ ac( x)- _ y p_ z^(k)_ ac( y) ) \|^2 ≤32dλtlog⁡2νk+2σ(k)4Mac(k)2‖−‖2 ≤ 32d _t 2 _k+2 _ s^(k)^4M_ z^(k)_ ac^2 \| x- y \|^2 (69) where, Mac(k)M_ z^(k)_ ac is the smoothness constant of log⁡pac(k) p_ z^(k)_ ac that can be derived using a similar proof procedure as Lemma 9, and λt=−Mt+1σ(k)2 _t=-M_t+ 1 _ s^(k)^2. Following the similar procedure as in proof of Lemma 9, we get, Mac(k)≤Mt1−σ(k)2MtM_ z^(k)_ ac≤ M_t1- _ s^(k)^2M_t (70) Plugging this leads to the lemma. ∎ Lemma 12. Let Htw(k):dc(k)↦tw(k)H^(k)_ tw: z^(k)_ dc z^(k)_ tw denote the projection function using Tweedie’s lemma defined in Algorithm 1. Then, we have the following ‖(Htw(k)−I)()−(Htw(k)−I)()‖22≤(ϵtw(k))2‖−‖22+δtw(k)2 \|(H^(k)_ tw-I)( x)-(H^(k)_ tw-I)( y) \|_2^2≤(ε^(k)_ tw)^2 \| x- y \|_2^2+δ^(k)_ tw^2 (71) for any ,∈ x, y , where (ϵtw(k))2=(σ(k))4Mt2(ε^(k)_ tw)^2=(σ^(k))^4M_t^2, and δtw(k)2=0δ^(k)_ tw^2=0. Proof. From Tweedie’s lemma, we have tw(k) z^(k)_ tw =[0|t=dc(k)] =E[ z_0| z_t= z^(k)_ dc] =dc(k)+(σ(k))2∇log⁡pt(dc(k)) = z^(k)_ dc+(σ^(k))^2∇ p_t( z^(k)_ dc) where, t∈[0,T]t∈[0,T] such that σ(k)=σ(t)σ^(k)=σ(t). Then, the residuals are given by Rtw(k)() R^(k)_ tw( x) =(σ(k))2∇log⁡pt() =(σ^(k))^2∇ p_t( x) (72) Rtw(k)() R^(k)_ tw( y) =(σ(k))2∇log⁡pt() =(σ^(k))^2∇ p_t( y) (73) Now, the norm of the difference of the residuals can be written as ‖Rtw(k)()−Rtw(k)()‖22 \|R^(k)_ tw( x)-R^(k)_ tw( y) \|_2^2 =(σ(k))4‖∇log⁡pt()−∇log⁡pt()‖22 =(σ^(k))^4 \|∇ p_t( x)-∇ p_t( y) \|_2^2 (74) ≤(σ(k))4Mt2‖−‖22 ≤(σ^(k))^4M_t^2 \| x- y \|_2^2 (75) where, MtM_t is the smoothness constant of log⁡pt() p_t( x). ∎ D.1 Main proof: Proof. Part (a) Using Lemma 10, 11, and 12, with probability at least 1−2e−νk1-2e^- _k, we have ‖Rσ(k)()−Rσ(k)()‖22 \|R_σ^(k)( x)-R_σ^(k)( y) \|_2^2 ≤ ≤ 3‖Rac(k)()−Rac(k)()‖22+3‖Rdc(k)()−Rdc(k)()‖22+3‖Rtw(k)()−Rtw(k)()‖22 3 \|R^(k)_ ac( x)-R^(k)_ ac( y) \|_2^2+3 \|R^(k)_ dc( x)-R^(k)_ dc( y) \|_2^2+3 \|R^(k)_ tw( x)-R^(k)_ tw( y) \|_2^2 ≤ ≤ 3((ϵac(k))2+(ϵdc(k))2+(ϵtw(k))2)‖−‖22+3(δac(k)2+δdc(k)2+δtw(k)2) 3((ε^(k)_ ac)^2+(ε^(k)_ dc)^2+(ε^(k)_ tw)^2) \| x- y \|_2^2+3(δ^(k)_ ac^2+δ^(k)_ dc^2+δ^(k)_ tw^2) ≤ ≤ ϵk2‖−‖22+δk2 _k^2 \| x- y \|_2^2+ _k^2 (76) where, ϵk2=3((2Mtσ(k)21−σ(k)2Mt)2+(σ(k))4Mt2) _k^2=3 ( ( 2M_t _ s^(k)^21- _ s^(k)^2M_t )^2+(σ^(k))^4M_t^2 ) δk2=3(2(σ(k))2(d+2dνk+2νk)+32dσ(k)2(1−Mtσ(k)2)log⁡2νk) _k^2=3 (2(σ^(k))^2(d+2 d _k+2 _k)+ 32d _ s^(k)^2(1-M_t _ s^(k)^2) 2 _k ) Using Lemma 9 leads to the final theorem. Part (b). Let us set νk=ln⁡(2π26η)+2ln⁡k _k= ( 2π^26η )+2 k. With this, the above weakly nonexpansiveness holds for all k∈ℕ+k ^+ with probability at least 1−∑k=1∞2e−ln⁡(2π26η)−2ln⁡k 1- _k=1^∞2e^- ( 2π^26η )-2 k = = 1−6ηπ2×∑k=1∞1k2 1- 6ηπ^2× _k=1^∞ 1k^2 (77) Using Riemann zeta function (Titchmarsh and Heath-Brown, 1986) at value 2, we have ζ(2)=∑k=1∞1k2=π26ζ(2)= _k=1^∞ 1k^2= π^26 (78) Plugging this in we get the probability to be at least 1−η1-η. Now, combining the results with Theorem 1 leads to the final proof of this part. ∎ Appendix E Proof of Theorem 3 Here, show that our 3-step AC-DC denoiser is bounded with high probability. We first show that each of 3 steps are bounded, and then combined them to establish the boundedness of our AC-DC denoiser as a whole. And following the boundedness, we show that AC-DC ADMM-PnP converges to a fixed point with proper scheduling of σ(k)σ^(k) and σ(k) _ s^(k). Lemma 13 (Uniform score bound). Suppose Assumption 2 holds. Let D:=diam()=sup,∈‖−‖2<∞ and S:=inf∈∥∇log⁡pdata()∥2<∞.D:= diam(X)= _ x, y \| x- y \|_2<∞ and S:= _ x \|∇ p_ data( x) \|_2<∞. Then, with L=MD+SL=MD+S, we have sup∈‖∇log⁡pdata()‖∞≤L. _ x \|∇ p_ data( x) \|_∞≤ L. Proof. From Assumption 2, we have ‖∇log⁡pdata()−∇log⁡pdata()‖2≤M‖−‖2,∀,∈ \|∇ p_ data( x)-∇ p_ data( y) \|_2≤ M \| x- y \|_2,\;∀ x, y (79) Fix any 0∈ x_0 . By the triangle inequality, for all ∈ x , ‖∇log⁡pdata()‖2 \|∇ p_ data( x) \|_2 ≤‖∇log⁡pdata()−∇log⁡pdata(0)‖2+‖∇log⁡pdata(0)‖2 ≤ \|∇ p_ data( x)-∇ p_ data( x_0) \|_2+ \|∇ p_ data( x_0) \|_2 ≤M‖−0‖2+‖∇log⁡pdata(0)‖2 ≤ M \| x- x_0 \|_2+ \|∇ p_ data( x_0) \|_2 (80) Taking the supremum over ∈ x and then the infimum over 0∈ x_0 yields sup∈‖∇log⁡pdata()‖2 _ x \|∇ p_ data( x) \|_2 ≤sup∈M‖−0‖2+inf0∈‖∇log⁡pdata(0)‖2 ≤ _ x M \| x- x_0 \|_2+ _ x_0 \|∇ p_ data( x_0) \|_2 (81) ≤MD+S ≤ MD+S (82) Because ‖∞≤‖2 \| u \|_∞≤ \| u \|_2 for any vector, sup∈‖∇log⁡pdata()‖∞≤MD+S _ x \|∇ p_ data( x) \|_∞≤ MD+S (83) which proves the above lemma. ∎ Lemma 14. Assuming ∥∇logpdata()∥∞≤,∀∈ \|∇ p_ data( x) \|_∞≤,∀ x , the score of intermediate noisy distributions ptt∈[0,T]\p_t\_t∈[0,T] are bounded as. ‖∇log⁡pt()‖2≤dL \|∇ p_t( x) \|_2≤ dL (84) for all ∈ x . Proof. We have t=0+σ(t) x_t= x_0+σ(t) n with ∼(,) n ( 0, I). Let us denote 1=σ(t) n_1=σ(t) n. Then, the marginal distribution is given by the convolution of two distributions. pt()=∫1∈p1(1)p0(−1)1 p_t( x)= _ x_1 p_ n_1( x_1)p_0( x- x_1)d x_1 (85) Then, the score is given by ∇log⁡pt() ∇ p_t( x) =∇pt()pt() = ∇ p_t( x)p_t( x) =1pt()∫1∈p1(1)∇p0(−1)1 = 1p_t( x) _ x_1 p_ n_1( x_1) _ xp_0( x- x_1)d x_1 =1pt()∫1∈p1(1)p0(−1)∇p0(−1)p0(−1)1 = 1p_t( x) _ x_1 p_ n_1( x_1)p_0( x- x_1) _ xp_0( x- x_1)p_0( x- x_1)d x_1 =1pt()∫1∈p1(1)p0(−1)∇log⁡p0(−1)1 = 1p_t( x) _ x_1 p_ n_1( x_1)p_0( x- x_1) _ x p_0( x- x_1)d x_1 (86) Now the norm can be bounded as ‖∇log⁡pt()‖2 \|∇ p_t( x) \|_2 ≤1pt()∫1∈p1(1)p0(−1)‖∇log⁡p0(−1)‖21 ≤ 1p_t( x) _ x_1 p_ n_1( x_1)p_0( x- x_1) \| _ x p_0( x- x_1) \|_2d x_1 ≤sup2∈‖∇log⁡p0(2)‖21pt()∫1∈p1(1)p0(−1)1 ≤ _ x_2 \| _ x p_0( x_2) \|_2 1p_t( x) _ x_1 p_ n_1( x_1)p_0( x- x_1)d x_1 =dL = dL (87) The final equality is due to the fact ‖2≤d‖∞ \| x \|_2≤ d \| x \|_∞. ∎ Lemma 15. Let Hac(k):~(k)↦ac(k)H^(k)_ ac: z^(k) z^(k)_ ac denote the function corresponding to approximate correction to noise level σ(k)σ^(k) defined in Algorithm 1. Then, with probability at least 1−e−ν1-e^-ν, the following holds for any ,∈ x, y 1d‖(Hac(k)−I)()‖22≤(σ(k))2(1+2ν+2ν) 1d \|(H^(k)_ ac-I)( x) \|_2^2≤(σ^(k))^2(1+2 ν+2ν) (88) Proof. For any ∈ x , we have the residual Rac(k)()=(Hac(k)−I)()=σ(k),∼(,)R^(k)_ ac( x)=(H^(k)_ ac-I)( x)=σ^(k) n,\; n ( 0, I). Then, ‖Rac(k)()‖22=(σ(k))2‖22=(σ(k))2χd2 \|R^(k)_ ac( x) \|_2^2=(σ^(k))^2 \| n \|_2^2=(σ^(k))^2 _d^2 (89) where χd2 _d^2 is standard chi-square distribution with d degrees of freedom. From Laurent and Massart (2000)[Lemma 1], the following holds with probability at least 1−e−ν1-e^-ν χd2≤d+2dν+2ν _d^2≤ d+2 dν+2ν (90) This implies 1d∥Hac(k)−I)()∥22≤(σ(k))2(1+21ν+2ν) 1d \|H^(k)_ ac-I)( x) \|_2^2≤(σ^(k))^2(1+2 1ν+2ν) with probability at least 1−e−ν1-e^-ν due to d≥1d≥ 1. ∎ Lemma 16. Let Hdc(k):ac(k)↦dc(k)H^(k)_ dc: z^(k)_ ac z^(k)_ dc denote the function corresponding to fine correction defined in Algorithm 1. Assume ‖∇log⁡pdata()‖∞≤L \|∇ p_ data( x) \|_∞≤ L. Then, with probability at least 1−e−ν1-e^-ν, the following holds for any ,∈ x, y : 1d‖(Hdc(k)−I)()‖22≤8λtlog⁡2ν+σ(k)4L2 1d \|(H^(k)_ dc-I)( x) \|_2^2≤ 8 _t 2ν+ _ s^(k)^4L^2 (91) where, λt=−Mt+1σ(k)2 _t=-M_t+ 1 _ s^(k)^2. Proof. From (67), with probability at least 1−2ν1-2ν, the norm can be bounded as ‖Rdc(k)()‖22 \|R^(k)_ dc( x) \|_2^2 ≤8dλtlog⁡2ν+σ(k)4‖∇log⁡pac(k)()‖22 ≤ 8d _t 2ν+ _ s^(k)^4 \| _ x p_ z^(k)_ ac( x) \|_2^2 (92) Then, using Lemma 14, we have 1d‖Rdc(k)()‖22≤8λtlog⁡2ν+σ(k)4L2 1d \|R^(k)_ dc( x) \|_2^2≤ 8 _t 2ν+ _ s^(k)^4L^2 (93) ∎ Lemma 17. Let Htw(k):dc(k)↦tw(k)H^(k)_ tw: z^(k)_ dc z^(k)_ tw denote the projection function using Tweedie’s lemma defined in Algorithm 1. Assume ‖∇log⁡pdata()‖∞≤L \|∇ p_ data( x) \|_∞≤ L. Then, we have the following 1d‖(Htw(k)−I)()‖22≤(σ(k))4L2 1d \|(H^(k)_ tw-I)( x) \|_2^2≤(σ^(k))^4L^2 (94) for any ∈ x . Proof. From Tweedie’s lemma, we have tw(k) z^(k)_ tw =[0|t=dc(k)] =E[ z_0| z_t= z^(k)_ dc] =dc(k)+(σ(k))2∇log⁡pt(dc(k)) = z^(k)_ dc+(σ^(k))^2∇ p_t( z^(k)_ dc) (95) where, t∈[0,T]t∈[0,T] such that σ(k)=σ(t)σ^(k)=σ(t). Now, the norm of residual can be written as ‖Rtw(k)()‖22 \|R^(k)_ tw( x) \|_2^2 =‖(σ(k))2∇log⁡pt()‖22 = \|(σ^(k))^2∇ p_t( x) \|_2^2 ≤(σ(k))4L2⋅d ≤(σ^(k))^4L^2· d (96) where, L bound of gradient from Lemma 14. ∎ E.1 Main proof Combining Lemmas 15, 16 and 17 leads to the proof of part (a) of Theorem 3. With probability at least 1−2e−νk1-2e^- _k, the denoiser satisfies the bounded residual condition. /1d‖(Dσ(k)−I)()‖22≤ck2. \|(D_σ^(k)-I)( x) \|_2^2≤ c_k^2. (97) Let’s define the relative residue as: βk:=1d(‖(k)−(k−1)‖2+‖(k)−(k−1)‖2+‖(k)−(k−1)‖2) _k:= 1 d ( \| x^(k)- x^(k-1) \|_2+ \| z^(k)- z^(k-1) \|_2+ \| u^(k)- u^(k-1) \|_2 ) (98) For any η∈[0,1)η∈[0,1) and a constant γ>1γ>1, the penalty parameter ρk _k is adjusted at each iteration k according to following rule (Chan et al., 2016): ρk+1=γρk if βk+1≥ηβk (Case 1)ρk else (Case 2) _k+1= casesγ _k& if _k+1≥η _k (Case 1)\\ _k& else (Case 2) cases (99) The PnP-ADMM with adaptive penalty involves two cases as shown above. At iteration k, if Case 1 holds, then by Lemma 18 we have βk+1≤6ck+2ck−1+2Rρk _k+1≤ 6c_k+2c_k-1+ 2R _k (100) On the other hand if Case 2 holds, then, βk+1≤ηβk _k+1≤η _k (101) Define ak=6ck+2ck−1+2Rρka_k=6c_k+2c_k-1+ 2R _k. Combining two cases, we get βk+1≤δβk+ak,δ=η,if Case 2 holds at iteration k0,if Case 1 holds at iteration k _k+1≤δ _k+a_k, δ= casesη, &if Case 2 holds at iteration $k$\\ 0, &if Case 1 holds at iteration $k$ cases (102) Note that for Case 1 ρk+1=γρk _k+1=γ _k and with γ>1γ>1, we get limk→∞cρk=0 _k→∞ c _k=0. In addition, with νk=ln⁡2π26η+2ln⁡k _k= 2π^26η+2 k, and the scheduling of σ(k)σ^(k), σ(k) _ s^(k) that satisfies limk→∞(σ(k))2(2+4νk+4νk)=0,limk→∞σ(k)21−Mσ(k)2log⁡2νk=0,limk→∞σ(k)=0,limk→∞σ(k)=0 _k→∞(σ^(k))^2(2+4 _k+4 _k)=0,\; _k→∞ _ s^(k)^21-M _ s^(k)^2 2 _k=0,\; _k→∞σ^(k)=0, _k→∞ _ s^(k)=0 results in limk→∞ak=0 _k→∞a_k=0. As k→∞k→∞, 3 different scenarios could occur. Let’s analyze each of the scenarios one by one. • Scene 1: Case 1 occurs infinitely many times and Case 2 occurs finitely many times. When Scene 1 occurs, then, there exists a constant K1>0K_1>0 such that βk+1≤ak _k+1≤ a_k for all k≥K1k≥ K_1. Since limk→∞ak=0 _k→∞a_k=0, this leads to limk→∞βk=0 _k→∞ _k=0. • Scene 2: Case 2 occurs infinitely many times and Case 1 occurs finitely many times. Similarly, there exists a constant K2>0K_2>0 such that βk+1≤ηβk _k+1≤η _k for all k≥K2k≥ K_2. Then, we have βk≤ηk−K2βK2 _k≤η^k-K_2 _K_2 (103) And with η∈[0,1)η∈[0,1), we have limk→∞βk=0 _k→∞ _k=0. • Scene 3: Both Case 1 and Case 2 occurs infinitely many times. With the Scene 1 and Scene 2 converging, the sequence limk→∞βk=0 _k→∞ _k=0 under this Scene as well. This proves the part (b) of Theorem 3. Lemma 18. For any iteration k that falls into Case 1, the following holds βk+1≤6ck+2ck−1+2Rρk _k+1≤ 6c_k+2c_k-1+ 2R _k (104) Proof. Consider the subproblem (7a) with adaptive penalty parameter ρk _k (as defined in (99)), (k+1)=argmin1ρkℓ(||())+12∥−(k)+(k)∥22 x^(k+1)= _ x 1 _k ( y|| A( x))+ 12 \| x- z^(k)+ u^(k) \|_2^2 (105) With the first order optimality condition, the solution (k+1) x^(k+1) satisfies 1ρk∇ℓ(||())|=(k+1)+((k+1)−(k)+(k))= 1 _k _ x ( y|| A( x)) |_ x= x^(k+1)+( x^(k+1)- z^(k)+ u^(k))= 0 (106) ⟹1d∥(k+1)−(k)+(k)∥2=1ρkd∥∇ℓ(||())|=(k+1)∥2 1 d \| x^(k+1)- z^(k)+ u^(k) \|_2= 1 _k d \| _ x ( y|| A( x)) |_ x= x^(k+1) \|_2 (107) Using the assumption of existence of R<∞R<∞ such that ∥∇ℓ(||())∥2/d≤R,∀∈ \| _ x ( y|| A( x)) \|_2/ d≤ R,\;∀ x , we get 1d‖(k+1)−(k)+(k)‖2≤Rρk 1 d \| x^(k+1)- z^(k)+ u^(k) \|_2≤ R _k (108) Since the denoiser Dσ(k)D_σ^(k) is bounded with probability at least 1−2e−νk1-2e^- _k, we have 1d‖(Dσ(k)−I)()‖2≤ck 1 d \|(D_σ^(k)-I)( x) \|_2≤ c_k (109) Now, 1d‖(k+1)−(k+1)+(k)‖2 1 d \| x^(k+1)- z^(k+1)+ u^(k) \|_2 =1d‖(k+1)+(k)−Dσ(k)((k+1)+(k))‖2 = 1 d \| x^(k+1)+ u^(k)-D_σ^(k)( x^(k+1)+ u^(k)) \|_2 (110) =1d‖(Dσ(k)−I)((k+1)+(k))‖2 = 1 d \|(D_σ^(k)-I)( x^(k+1)+ u^(k)) \|_2 (111) ≤ck ≤ c_k (112) Now, using triangle inequality, we can bound ‖(k+1)−(k)‖2 \| z^(k+1)- z^(k) \|_2 as 1d‖(k+1)−(k)‖2 1 d \| z^(k+1)- z^(k) \|_2 =1d‖(k+1)−(k+1)−(k)+(k+1)+(k)−(k)‖2 = 1 d \| z^(k+1)- x^(k+1)- u^(k)+ x^(k+1)+ u^(k)- z^(k) \|_2 (113) ≤Rρk+ck ≤ R _k+c_k (114) Similarly, it can be shown that 1d‖(k+1)‖2 1 d \| u^(k+1) \|_2 =1d‖(k)+(k+1)−(k+1)‖2 = 1 d \| u^(k)+ x^(k+1)- z^(k+1) \|_2 (115) =1d‖(k)+(k+1)−Dσ(k)((k+1)+(k))‖2 = 1 d \| u^(k)+ x^(k+1)-D_σ^(k)( x^(k+1)+ u^(k)) \|_2 (116) =1d‖(Dσ(k)−I)((k+1)+(k))‖2 = 1 d \|(D_σ^(k)-I)( x^(k+1)+ u^(k)) \|_2 (117) =ck =c_k (118) This implies 1d‖(k+1)−(k)‖2≤2ck 1 d \| u^(k+1)- u^(k) \|_2≤ 2c_k. Finally, we use (k+1)=(k+1)−(k)+(k+1) x^(k+1)= u^(k+1)- u^(k)+ z^(k+1) to obtain 1d‖(k+1)−(k)‖2 1 d \| x^(k+1)- x^(k) \|_2 (119) = = 1d‖(k+1)−(k)+(k+1)−(k)+(k−1)−(k)‖2 1 d \| u^(k+1)- u^(k)+ z^(k+1)- u^(k)+ u^(k-1)- z^(k) \|_2 (120) ≤ ≤ 1d‖(k+1)−(k)‖2+1d‖(k)−(k−1)‖2+1d‖(k+1)−(k)‖2 1 d \| u^(k+1)- u^(k) \|_2+ 1 d \| u^(k)- u^(k-1) \|_2+ 1 d \| z^(k+1)- z^(k) \|_2 (121) ≤ ≤ 2ck+2ck−1+Rρk+ck 2c_k+2c_k-1+ R _k+c_k (122) = = 3ck+2ck−1+Rρk 3c_k+2c_k-1+ R _k (123) Combining all the bounds results using triangle inequality results in βk+1≤6ck+2ck−1+2Rρk _k+1≤ 6c_k+2c_k-1+ 2R _k (124) where, ck=(σ(k))2(2+4νk+4νk)+16σ(k)2/1−Mσ(k)2log⁡2/νk+2σ(k)4L2+2(σ(k))4L2c_k=(σ^(k))^2(2+4 _k+4 _k)+ 16 _ s^(k)^21-M _ s^(k)^2 2 _k+2 _ s^(k)^4L^2+2(σ^(k))^4L^2 ∎ Remark 1. While the proposed method and its theoretical results are based on Variance Exploding (VE) scheduling, they can be easily extended to Variance Preserving (VP) scheduling case (Karras et al., 2022). E.2 Theoretical Results with finite DC steps J Lemma 19. Let Hdc(k):ac(k)↦dc(k)H^(k)_ dc: z^(k)_ ac z^(k)_ dc denote the function corresponding to fine correction as defined in Algorithm 1 with finite J and η(k)≤2σ(k)2η^(k)≤ 2 _ s^(k)^2. Also, let π(k)=p(σ(k)|ac(k))π^(k)=p( z_σ^(k)| z^(k)_ ac) be the stationary target distribution and π~0(k) π^(k)_0 be initial distribution used for the DC at iteration k. Then, with probability at least 1−e−νk1-e^- _k, the following holds for any ,∈ x, y if 1/σ(k)2<Mt 1 _ s^(k)^2<M_t: ‖(Hdc(k)−I)()−(Hdc(k)−I)()‖22≤(δdc(k))2+(ϵdc(k))2‖−‖22 \|(H^(k)_ dc-I)( x)-(H^(k)_ dc-I)( y) \|_2^2≤(δ^(k)_ dc)^2+(ε^(k)_ dc)^2 \| x- y \|_2^2 (125) where,0<κ<10<κ<1, C>0C>0, (δdc(k))2=64dσ(k)2(1−Mtσ(k)2)log⁡2νk+C(1−κ)2J22(π~0(k),π(k))+O((η(k))2)(δ^(k)_ dc)^2= 64d _ s^(k)^2(1-M_t _ s^(k)^2) 2 _k+C(1-κ)^2JW_2^2( π^(k)_0,π^(k))+O ((η^(k))^2 ), and (ϵdc(k))2=(22Mtσ(k)21−Mtσ(k)2)2(ε^(k)_ dc)^2= ( 2 2M_t _ s^(k)^21-M_t _ s^(k)^2 )^2. Proof. With finite J, the langevin dynamics doesn’t necessarily converge to the stationary distribution π(k)=p(σ(k)|ac(k))π^(k)=p( z_σ^(k)| z^(k)_ ac). Let π~0(k) π^(k)_0 be the initial distribution used to initialize the finite step langevin dynamics and π~(k) π^(k) be the distribution of the iterate after running finite J steps of langevin dynamics. Using (Dalalyan and Karagulyan, 2019)[Theorem 1], the following holds for η(k)≤2σ(k)2η^(k)≤ 2 _ s^(k)^2 2(π~(k),π(k))≤(1−κ)J2(π~0(k),π(k))+O(η(k)) _2( π^(k),π^(k))≤(1-κ)^JW_2( π^(k)_0,π^(k))+O (η^(k) ) (126) where, 0<κ<10<κ<1 and 2W_2 is the 2-Wasserstein distance. Using Kantorovich and Rubinstein dual representation (Villani and others, 2008), ‖π~(k)[]−π(k)[]‖2 \|E_ π^(k)[ z]-E_π^(k)[ z] \|_2 =‖∫dc(k)π~(k)()−∫π(k)()‖2 = \| z^(k)_ dcd π^(k)( z)- zdπ^(k)( z) \|_2 (127) ≤1(π~(k),π(k)) _1( π^(k),π^(k)) (128) ≤2(π~(k),π(k)) _2( π^(k),π^(k)) (129) The last inequality is due to the Holder’s inequality Villani and others (2008). Using triangle inequality leads to following: ‖dc(k)−π(k)[]‖2 \| z^(k)_ dc-E_π^(k)[ z] \|_2 ≤‖dc(k)−π~(k)[]‖2+‖π~(k)[]−π(k)[]‖2 ≤ \| z^(k)_ dc-E_ π^(k)[ z] \|_2+ \|E_ π^(k)[ z]-E_π^(k)[ z] \|_2 (130) ≤‖dc(k)−π~(k)[]‖2+(1−κ)J2(π~0(k),π(k))+O(η(k)) ≤ \| z^(k)_ dc-E_ π^(k)[ z] \|_2+(1-κ)^JW_2( π^(k)_0,π^(k))+O (η^(k) ) (131) Using this result and following the same procedure as in Lemma 11, we get the theorem. ∎ Theorem 4. Suppose that the assumptions in Theorem 1, Assumption 2, and Assumption 3 hold. Further, assume that the DC steps finite steps J and η(k)≤2σ(k)2η^(k)≤ 2 _ s^(k)^2. Also, let π(k)=p(σ(k)|ac(k))π^(k)=p( z_σ^(k)| z^(k)_ ac) be the stationary target distribution and π~0(k) π^(k)_0 be initial distribution used for the DC at iteration k. Let Dσ(k):~(k)↦tw(k)D_σ^(k): z^(k) z^(k)_ tw denote the AC-DC denoiser. Then, we have: With probability at least 1−2e−νk1-2e^- _k, the following holds for iteration k of ADMM-PnP: ‖(Dσ(k)−I)()−(Dσ(k)−I)()‖22≤ϵk2‖−‖22+δk2 \|(D_σ^(k)-I)( x)-(D_σ^(k)-I)( y) \|_2^2≤ε^2_k \| x- y \|_2^2+δ^2_k (132) for any ,∈ x, y , k∈ℕ+k ^+, a constant 0<κ<10<κ<1 and C>0C>0, when σ(k)2+(σ(k))2<1/M _ s^(k)^2+(σ^(k))^2<1/M with ϵk2 ε^2_k =3((22Mσ(k)2/(1−Mσ(k)2))2+(σ(k))4M2) =3(( 2 2M _ s^(k)^2(1-M _ s^(k)^2))^2+(σ^(k))^4M^2) (133) δk2 δ^2_k =3(2(σ(k))2(d+2dνk+2νk)+64dσ(k)2/(1−Mσ(k)2)log⁡2/νk)+ =3(2(σ^(k))^2(d+2 d _k+2 _k)+ 64d _ s^(k)^2(1-M _ s^(k)^2) 2 _k)+ C(1−κ)2J22(π~0(k),π(k))+O((η(k))2). C(1-κ)^2JW_2^2( π^(k)_0,π^(k))+O ((η^(k))^2 ). (134) In other words, with νk=ln⁡2π/6η+2ln⁡k _k= 2π6η+2 k, the denoiser Dσ(k)D_σ^(k) satisfies part (a) for all k∈ℕ+k ^+ with probability at least 1−η1-η. Proof. By substituting Lemma 11 with Lemma 19 leads to the theorem. ∎ Theorem 5. Suppose that Assumptions 2-3 hold. Let D:=diam()=sup,∈‖−‖2<∞,S:=inf∈∥∇log⁡pdata()∥2<∞D:= diam(X)= _ x, y \| x- y \|_2<∞,\;S:= _ x \|∇ p_ data( x) \|_2<∞ and define L:=MD+SL:=MD+S. Let Dσ(k):~(k)↦tw(k)D_σ^(k): z^(k) z^(k)_ tw denote the AC-DC denoiser. Further, assume that the DC steps finite steps J and η(k)≤2σ(k)2η^(k)≤ 2 _ s^(k)^2. Also, let π(k)=p(σ(k)|ac(k))π^(k)=p( z_σ^(k)| z^(k)_ ac) be the stationary target distribution and π~0(k) π^(k)_0 be initial distribution used for the DC at iteration k. Then, the following hold: (Boundedness) With probability at least 1−2e−νk1-2e^- _k, the denoiser Dσ(k)D_σ^(k) is bounded at each iteration k i.e. 1d‖(Dσ(k)−I)()‖22≤ck2 1d \|(D_σ^(k)-I)( x) \|_2^2≤ c_k^2 whenever σ(k)2+(σ(k))2<1/M, _ s^(k)^2+(σ^(k))^2<1/M, where ck=(σ(k))2(2+4νk+4νk)+32σ(k)2/(1−Mσ(k)2)log⁡2/νk+C(1−κ)2J22(π~0(k),π(k))+O((η(k))2)+4L2σ(k)4+2(σ(k))4L2c_k=(σ^(k))^2(2+4 _k+4 _k)+ 32 _ s^(k)^2(1-M _ s^(k)^2) 2 _k+C(1-κ)^2JW_2^2( π^(k)_0,π^(k))+O ((η^(k))^2 )+4L^2 _ s^(k)^4+2(σ^(k))^4L^2, 0<κ<10<κ<1, C>0C>0 and νk>0 _k>0. Let νk=ln⁡2π26η+2ln⁡k _k= 2π^26η+2 k with η∈(0,1]η∈(0,1]. Consequently, the denoiser Dσ(k)D_σ^(k) is bounded for all k∈ℕ+k _+ with corresponding ckc_k and probability at least 1−η1-η. Proof. The proof follows similar as in Theorem 3 by incorporating the effect of finite J in Lemma 16 as done in Lemma 19. ∎ Appendix F Theoretical results for ODE based denoiser Refer to the Zhang et al. (2024) for details on ODE based denoiser. F.1 Theoretical results equivalent to Theorem 2 Lemma 20. Let Hode(k):dc(k)↦ode(k)H^(k)_ ode: z^(k)_ dc z^(k)_ ode denote the projection function using ode based denoiser (Karras et al., 2022) in Algorithm 1. Then, we have the following ‖(Hode(k)−I)()−(Hode(k)−I)()‖22≤(ϵode(k))2‖−‖22+δode(k)2 \|(H^(k)_ ode-I)( x)-(H^(k)_ ode-I)( y) \|_2^2≤(ε^(k)_ ode)^2 \| x- y \|_2^2+δ^(k)_ ode^2 (135) for any ,∈ x, y with (ϵode(k))2=2(∫t=tσ(k)0(σ(t)σ′(t)Mt)2t)(ε^(k)_ ode)^2=2 ( _t=t_σ^(k)^0 (σ(t)σ (t)M_t )^2dt ), and δode(k)2=0δ^(k)_ ode^2=0. Proof. Then, the difference of residual of ode projection i.e. Rode(k)=Hode(k)−IR^(k)_ ode=H^(k)_ ode-I can be bounded as ‖Rode(k)()−Rode(k)()‖22= \|R^(k)_ ode( x)-R^(k)_ ode( y) \|_2^2= ‖∫t=tσ(k)0−σ(t)σ′(t)(∇log⁡pt()−∇log⁡pt())dt‖22 \| _t=t_σ^(k)^0-σ(t)σ (t)(∇ p_t( x)-∇ p_t( y))dt \|_2^2 ≤ ≤ 2‖∫t=tσ(k)0−σ(t)σ′(t)(∇log⁡pt()−∇log⁡pt())dt‖22 2 \| _t=t_σ^(k)^0-σ(t)σ (t)(∇ p_t( x)-∇ p_t( y))dt \|_2^2 ≤ ≤ 2∫t=tσ(k)0(σ(t)σ′(t))2‖(∇log⁡pt()−∇log⁡pt())‖22t 2 _t=t_σ^(k)^0 (σ(t)σ (t) )^2 \|(∇ p_t( x)-∇ p_t( y)) \|_2^2dt ≤ ≤ 2∫t=tσ(k)0(σ(t)σ′(t))2Mt2‖−‖22t 2 _t=t_σ^(k)^0 (σ(t)σ (t) )^2M_t^2 \| x- y \|_2^2dt ≤ ≤ 2(∫t=tσ(k)0(σ(t)σ′(t))2Mt2t)‖−‖22 2 ( _t=t_σ^(k)^0 (σ(t)σ (t) )^2M_t^2dt ) \| x- y \|_2^2 (136) ∎ Theorem 6. Suppose that the assumptions in Theorem 1, Assumption 2 and Assumption 3 hold. Further, assume that the step size satisfies η(k)→0η^(k)→ 0 and the number of iterations J→∞J→∞. Let Dσ(k):~(k)↦tw(k)D_σ^(k): z^(k) z^(k)_ tw denote the AC-DC denoiser. Then, we have: (a) With probability at least 1−2e−νk1-2e^- _k, the following holds for iteration k of ADMM-PnP: ‖(Dσ(k)−I)()−(Dσ(k)−I)()‖22≤ϵk2‖−‖22+δk2 \|(D_σ^(k)-I)( x)-(D_σ^(k)-I)( y) \|_2^2≤ε^2_k \| x- y \|_2^2+δ^2_k (137) for any ,∈ x, y and k∈ℕ+k ^+ when σ(k)2+(σ(k))2<1/M _ s^(k)^2+(σ^(k))^2<1/M with ϵk2 ε^2_k =3(2Mσ(k)2/1−σ(k)2M)2+6∫t=tσ(k)0(σ(t)σ′(t))2Mt2dt)) =3( 2M _ s^(k)^21- _ s^(k)^2M)^2+6 _t=t_σ^(k)^0 (σ(t)σ (t) )^2M_t^2dt)) (138) δk2 δ^2_k =3(2(σ(k))2(d+2dνk+2νk)+32dσ(k)2/(1−Mσ(k)2)log⁡2/νk). =3(2(σ^(k))^2(d+2 d _k+2 _k)+ 32d _ s^(k)^2(1-M _ s^(k)^2) 2 _k). (139) In other words, if νk=ln⁡2π/6η+2nk _k= 2π6η+2nk, the denoiser Dσ(k)D_σ^(k) satisfies part (a) for all k∈ℕ+k ^+ with probability at least 1−η1-η. (b) Assume that σ(k)σ^(k) is scheduled such that limk→∞(σ(k))2νk=0 _k→∞(σ^(k))^2 _k=0 for νk=ln⁡2π/6η+2nk _k= 2π6η+2nk, ϵ<1ε<1, and ϵ/μ(1+ϵ−2ϵ2)<1/ρ εμ(1+ε-2ε^2)< 1ρ all hold, where ϵ=limk→∞supϵkε= _k→∞ _k with ϵk _k defined in (138). Consequently, δ=limk→∞supδkδ= _k→∞ _k is finite and ADMM-PnP with the AC-DC denoiser with ode based denoiser converges to an r-ball (see r in Theorem 1) with probability at least 1−η1-η. Proof. The proof follow similar to the proof of Theorem 2 in Appendix D.1 with the residual bound of Tweedie’s lemma replaced by Lemma 20. ∎ F.2 Theoretical Results equivalent to Theorem 3 Lemma 21. Let Hode(k):dc(k)↦ode(k)H^(k)_ ode: z^(k)_ dc z^(k)_ ode denote the projection function using ode based denoiser (Karras et al., 2022) in Algorithm 1. Assume ‖∇log⁡pdata()‖∞≤L,∀∈ \|∇ p_ data( x) \|_∞≤ L,\;∀ x . Then, we have the following 1d∥(Hode(k)−I)()∥22≤L2∫t=tσ(k)0(σ(t)σ′(t)2dt 1d \|(H^(k)_ ode-I)( x) \|_2^2≤ L^2 _t=t_σ^(k)^0(σ(t)σ (t)^2dt (140) for any ∈ x . Proof. Then, the residual of ode projection i.e. Rode(k)=Hode(k)−IR^(k)_ ode=H^(k)_ ode-I can be bounded as ‖Rode(k)()‖22= \|R^(k)_ ode( x) \|_2^2= ‖∫t=tσ(k)0−σ(t)σ′(t)∇log⁡pt()dt‖22 \| _t=t_σ^(k)^0-σ(t)σ (t)∇ p_t( x)dt \|_2^2 ≤ ≤ d⋅L2‖∫t=tσ(k)0−σ(t)σ′(t)dt‖22 d· L^2 \| _t=t_σ^(k)^0-σ(t)σ (t)dt \|_2^2 ≤ ≤ d⋅L2‖∫t=tσ(k)0−σ(t)σ′(t)dt‖22 d· L^2 \| _t=t_σ^(k)^0-σ(t)σ (t)dt \|_2^2 ≤ ≤ d⋅L2∫t=tσ(k)0(σ(t)σ′(t)2dt d· L^2 _t=t_σ^(k)^0(σ(t)σ (t)^2dt (141) ∎ A theorem analogous to Theorem 3 can also be obtained for ODE-based denoiser. The only difference lies in the expression of the constant ckc_k, which in this case becomes ck=(σ(k))2(2+4νk+4νk)+16σ(k)2/1−Mσ(k)2log2/νk+2σ(k)4L2+2L2∫t=tσ(k)0(σ(t)σ′(t)2dt.c_k=(σ^(k))^2(2+4 _k+4 _k)+ 16 _ s^(k)^21-M _ s^(k)^2 2 _k+2 _ s^(k)^4L^2+2L^2 _t=t_σ^(k)^0(σ(t)σ (t)^2dt. Appendix G Usage of Large Language Models (LLM) An LLM was used solely to assist with polishing the writing. LLM played no part in the experiments, results and conclusion. Appendix H Experimental Details Table 3: Hyperparameter settings for each task Task ρ W lr of Adam in (7a) Superresolution (4×4×) 100100 100100 3×10−23× 10^-2 Gaussian Deblur 100100 100100 5×10−25× 10^-2 HDR 500500 100100 3×10−23× 10^-2 Inpainting (Random) 500500 100100 1×10−11× 10^-1 Inpainting (Box) 500500 100100 1×10−11× 10^-1 Motion Deblur 100100 100100 1×10−11× 10^-1 Nonlinear Deblur 300300 400400 3×10−13× 10^-1 Phase Retrieval 100100 400400 1×10−11× 10^-1 H.1 Details on task specific data-fidelity loss ℓ We use mean square error (MSE) as the data-fidelity loss for every task i.e. ℓ(||)=−logp(|)=12σn2∥−()∥22 ( y|| x)=- p( y| x)= 12 _n^2 \| y-A( x) \|_2^2 (142) H.2 Details on pretrained diffusion models The pretrained models provided in Chung et al. (2023) are used in our experiment. Refer to Chung et al. (2023) for more details on these pretrained models. H.3 Baseline Details Unless mentioned otherwise, we conduct the experiments in the default settings of their original implementation except for maintaining consistency within the measurement operators. • DDRM (Kawar et al., 2021): We use 2020 steps DDIM with η=0.85η=0.85 and ηb=1 _b=1 as specified in Kawar et al. (2022). • DPS (Chung et al., 2023) : The original implementation is ran in their default settings. • DiffPIR (Zhu et al., 2023): The default settings are adopted in the experiments. • RED-diff (Mardani et al., 2024): We use λ=0.25λ=0.25 and lr=0.5lr=0.5 as specified in the paper. • DAPS (Zhang et al., 2024): We use the best performing DAPS-4K version as proposed in the paper. • DPIR (Zhang et al., 2022): We employ "drunet_color" as PnP denoiser, while keeping all the other settings at their default values. • DCDP (Li et al., 2025): All the settings are set to their default values. • PMC (Sun et al., 2024): PMC was proposed using different score models for two different tasks with relatively high measurement SNR. For a fair comparison, we used our own implementation with the same score model checkpoints as our methods, and further tuned this method accordingly. H.4 Evaluation Metrics For all the methods, we use the implementation of PSNR, SSIM, and LPIPS provided in piq python package. The default settings for these metrics are used except the average pooling enabled for LPIPS. H.5 Computation Resource Details All the experiments were run on a instance equipped with one Nvidia H100 GPU, 20 cores of 2.0 Ghz Intel Xeon Platinum 8480CL CPU, and 64 GB of RAM. Appendix I Additional Experimental Results I.1 Illustration of proposed Denoiser Figure 6 illustrates the effect of the proposed correction-denoising procedure. The noisy input image ~(k) z^(k) typically lies far away from the Gaussian noise manifold, leading to poor denoising performance if directly used. To address this mismatch, our method first performs correction to effectively gaussianize the noise which is then denoised using Tweedie’s lemma or ode-style score integration, producing a high-quality clean reconstruction. Figure 6: Illustration of correction and denoising step in proposed method. I.2 Empirical validation of Assumption 2 and Assumption 3 To assess the practicality of the smoothness and coercivity assumptions used in Theorem 2, we conduct two diagnostic experiments using a pretrained score model on the validation split of the FFHQ dataset. These experiments are designed to evaluate (i) the empirical Lipschitz behavior of the score function ∇log⁡pdata()∇ p_data( x) (Assumption 2), and (i) the coercivity of the energy landscape −log⁡pdata()- p_data( x) (Assumption 3). Empirical smoothness of the score. We randomly 1000 samples of 1,2 x_1, x_2 and compute score differences ‖∇log⁡pdata(1)−∇log⁡pdata(2)‖2\|∇ p_data( x_1)-∇ p_data( x_2)\|_2 and image differences ‖1−2‖2\| x_1- x_2\|_2. Figure 7 plots the histogram of their ratio. The distribution concentrates around a finite value (mostly between 5050 and 160160), indicating that the score behaves approximately M-Lipschitz with a moderate empirical constant. This supports the smoothness requirement in Assumption 2. Figure 7: Histogram of score difference norm ratio: ‖∇log⁡pdata(2)−∇log⁡pdata(1)‖2/‖1−2‖2\|∇ p_data( x_2)-∇ p_data( x_1)\|_2/\| x_1- x_2\|_2, illustrating empirical smoothness (Assumption 2). Empirical coercivity. To evaluate coercivity, we scale images by factors c∈1,1.5,2,3c∈\1,1.5,2,3\ and measure the quantity ⟨,−∇log⁡pdata()⟩ x,-∇ p_data( x) as a function of the squared image norm ‖22\| x\|_2^2. As shown in Figure 8, the inner product grows approximately linearly with ‖22\| x\|_2^2, indicating that the learned score consistently pulls large-norm images back toward the data manifold. This behavior is consistent with the coercivity structure assumed in Assumption 3. Figure 8: Empirical coercivity test: relationship between ⟨,−∇log⁡pdata()⟩ x,-∇ p_data( x) and ‖22\| x\|_2^2. The strong positive correlation indicates coercive energy behavior (Assumption 3). Together, these empirical diagnostics demonstrate that the theoretical assumptions employed in our analysis hold approximately in practice and therefore justify the use of DC correction in our AC–DC algorithm. I.3 Illustration of Usage of Additional Regularization To further demonstrate the flexibility of integrating diffusion-based PnP denoisers within the ADMM framework, we present an example where we employ an additional perceptual regularization term which will be handled in the maximum-likelihood (ML) step. In particular, the x-update step of ADMM with an LPIPS perceptual regularization (Zhang et al., 2018) becomes: (k+1)=arg⁡min⁡1ρℓ(y∥())+12‖−(k)+(k)‖22+λlpipsLPIPSVGG(,ref), x^(k+1)= _ x 1ρ (y\|A( x))+ 12\| x- z^(k)+ u^(k)\|_2^2+ _lpips\,LPIPS_VGG( x, x_ref), (143) where ref x_ref is the reference image and λlpips _lpips controls the perceptual strength. This example highlights the flexibility of the proposed method: unlike traditional diffusion-based PnP approaches that struggle in the presence of dual variables, our design enables seamless incorporation of additional regularization terms. In Fig. 9 we illustrate box inpainting reconstruction task with the perceptual LPIPS-VGG regularization which enhances semantic content consistency while allowing visual style transfer from the reference images. Figure 9: Demonstration of incorporating additional perceptual regularization. I.4 Results on additional tasks The results on additional two tasks: HDR and nonlinear deblurring are presented in the Table 4. Table 4: Reconstruction metrics (100 images) on FFHQ / ImageNet for additional tasks. Bold: best, blue: 2nd best. FFHQ ImageNet Task Method PSNR↑ SSIM↑ LPIPS↓ PSNR↑ SSIM↑ LPIPS↓ HDR Ours–tweedie 27.425 0.853 0.164 26.515 0.817 0.182 DAPS 26.94 0.852 0.154 26.848 0.816 0.172 RED-diff 26.815 0.836 0.241 20.794 0.771 0.232 PMC 21.582 0.707 0.291 22.745 0.707 0.290 Nonlinear Deblur Ours–tweedie 29.326 0.823 0.185 27.837 0.725 0.212 DAPS 28.598 0.782 0.172 27.745 0.739 0.201 DPS 23.746 0.668 0.276 22.724 0.543 0.394 RED-diff 26.9 0.72 0.234 25.488 0.72 0.207 PMC 21.102 0.0623 0.354 22.347 0.533 0.430 I.5 Ablation study We perform the ablation study on the significance of our proposed correction steps. The results are presented in the Table 5. Table 5: Comparison of our method with and without correction on FFHQ. Best results are highlighted in bold. Ours-tweedie without correction Ours-tweedie with correction Tasks PSNR↑ SSIM↑ LPIPS↓ PSNR↑ SSIM↑ LPIPS↓ Superresolution (4x) 26.915 0.730 0.314 30.439 0.857 0.178 Gaussian Blur 28.896 0.788 0.275 30.402 0.853 0.175 Inpainting (Box) 15.604 0.617 0.361 24.025 0.859 0.131 Motion Deblur 25.123 0.538 0.370 30.003 0.854 0.179 Nonlinear Deblur 21.731 0.561 0.375 29.326 0.823 0.185 Phase Retrieval 11.978 0.181 0.726 27.944 0.793 0.209 I.6 Influence of decay schedule and NFE efficiency In our ADMM-PnP scheme, the size of the decay window for σ(k)σ^(k) determines the total number of iterations – and thus the speed of convergence. A shorter window (small W) drives σ(k)σ^(k) down more quickly, often reaching convergence in fewer steps but at the risk of settling in a suboptimal local minimum. To study this trade-off, we sweep W∈5,10,50,100,200,300,400,500W∈\5,10,50,100,200,300,400,500\ for each task. Since each iteration of Ours-tweedie uses 1111 score evaluations (10 for the DC update and 11 for the Tweedie’s lemma based denoiser), these W value translate to Number of Function Evaluations (NFE)=55,110,550,1100,2200,3300,4400,5500Number of Function Evaluations (NFE)=\55,110,550,1100,2200,3300,4400,5500\ By contrast, each Ours-ode iteration costs 2020 NFEs, giving Number of Function Evaluations (NFE)=100,200,1000,2000,4000,6000,8000,10000Number of Function Evaluations (NFE)=\100,200,1000,2000,4000,6000,8000,10000\ Figures 10-15 plot mean± . (standard deviation) performance of our methods and all baselines against NFE over 100100 images of FFHQ dataset. For most tasks, quality saturates after just 1010 iterations (110 NFE for Ours-tweedie, 200 NFE for Ours-ode), showing a rapid decay schedule suffices to achieve near-peak results. However, on the hardest inverse problems (phase retrieval and nonlinear blur), gradually decaying noise (larger W) and more NFEs yield significantly better reconstructions–far outpacing every baseline. Thus, while aggressive schedules excel on simple tasks, challenging problems benefit from extended iteration and gentler annealing; given enough NFEs, our approach establishes state-of-the-art performance across most of the tasks. Figure 10: Performance with respect to NFE for Superresolution task (FFHQ) Figure 11: Performance with respect to NFE for Gaussian deblurring task (FFHQ) Figure 12: Performance with respect to NFE for Inpainting with random missings (FFHQ) Figure 13: Performance with respect to NFE for Inpainting with box missing (FFHQ) Figure 14: Performance with respect to NFE for Motion blur (FFHQ) Figure 15: Performance with respect to NFE for Phase retrieval (FFHQ) I.7 More qualitative results Figure 16: Recovery under 4×4× superresolution task on FFHQ Figure 17: Recovery under inpainting with random missings on FFHQ Figure 18: Recovery under motion blur task on FFHQ

Instant research discovery