Convergent Bregman Plug-and-Play Image Restoration for Poisson Inverse Problems
Authors: Samuel Hurault
Ulugbek Kamilov
Arthur Leclaire
Nicolas Papadakis
Venue: NeurIPS 2023 poster (2023)
Keywords: Plug-and-Play
Poisson Inverse Problems
Bregman distance
Proximal Gradient Descent
nonconvex and nonsmooth optimization
Poisson inverse problems
TL;DR: We propose a Bregman generalization of the plug-and-play framework to solve image inverse problems involving non-Lipschitz gradient data-fidelity terms.

Abstract: Plug-and-Play (PnP) methods are efficient iterative algorithms for solving ill-posed image inverse problems. PnP methods are obtained by using deep Gaussian denoisers instead of the proximal operator or the gradient-descent step within proximal algorithms. Current PnP schemes rely on data-fidelity terms that have either Lipschitz gradients or closed-form proximal operators, which is not applicable to Poisson inverse problems. Based on the observation that the Gaussian noise is not the adequate noise model in this setting, we propose to generalize PnP using the Bregman Proximal Gradient (BPG) method. BPG replaces the Euclidean distance with a Bregman divergence that can better capture the smoothness properties of the problem. We introduce the Bregman Score Denoiser specifically parametrized and trained for the new Bregman geometry and prove that it corresponds to the proximal operator of a nonconvex potential. We propose two PnP algorithms based on the Bregman Score Denoiser for solving Poisson inverse problems. Extending the convergence results of BPG in the nonconvex settings, we show that the proposed methods converge, targeting stationary points of an explicit global functional. Experimental evaluations conducted on various Poisson inverse problems validate the convergence results and showcase effective restoration performance.

Convergent Bregman Plug-and-Play Image Restoration for Poisson Inverse Problems
Samuel Hurault
Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251 samuel huraul tmath u-boordeaxx...
Ulugbek Kamiloy
Washington University in St. Louis
kamilov@uustl.em
Arthur Leclaire
Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251 LTCI, Télécom Paris, IP Paris
arthur htt.:uu.....iireteeeeeeopprrrrr
Nicolas Papadakis
Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251
nicolas, papadakis u u-bordeaux.fF
Abstract
Plug-and-Play (PnP) methods are efficient iterative algorithms for solving ill-posed image inverse problems. PnP methods are obtained by using deep Gaussian denois- ers instead of the proximal operator the gradient-descent step within proximal algorithms. Current PnP schemes rely on data- fidelity terms that have either Lip- schitz gradients closed-form proximal operators, which is not applicable Poisson inverse problems. Based on the observation that the Gaussian noise is not the adequate noise model in this setting, propose generalize PnP using the Bregman Proximal Gradient (BPG) method. BPG replaces the Euclidean distance with Bregman divergence that can better capture the smoothness properties of the problem. We introduce the Bregman Score Denoiser specifically parametrized and trained for the new Bregman geometry and prove that it corresponds the proximal operator of nonconvex potential. We propose two PnP algorithms based on the Bregman Score Denoiser for solving Poisson inverse problems. Extending the convergence results of BPG in the nonconvex settings, we show that the pro- posed methods converge, targeting stationary points of an explicit global functional. Experimental evaluations conducted on various Poisson inverse problems validate the convergence results and showease effective restoration performance.
Introduction
Ill-posed image inverse problems classically formulated with minimization problem of the form arg min (f(x) + g(x)
(I)
=ER*
where ff is data-fidh lata-fidelity term, 9 regularization term, and À regularization parameter. The data-fidelity term generally written as the negative log- likelihoohd fxx = - (t)l)))) of the
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

probabilistic observation model chosen to describe the physics of an acquisition y R" from linear measurements Az of an image z R". In applications such as Positron Emission Tomography (PET) or astronomical CCD cameras [Bertero et al., 2009], where images are obtained by counting particles (photons or electrons), it is common to use the Poisson noise model y ~ P(a with parameter > 0. The corresponding negative log-likelihood corresponds to the Kullback-Leible divergence
m
f(x) - - y log
y/
+ a(Ax); - ys.
((Ax))
The minimization of (1) can be addressed with proximal splitting algorithms [Combettes and Pesquet, 2011]. Depending on the properties of the functions f and g. they consist in alternatively evaluating the proximal operator and/or performing gradient-descent step f and g.
Plug-and-Play (PnP) |Venkatakrishnan et al., 2013] and Regularization-by beiiiiiiiinnbb---Diiiiiinngg (RED) [Romano et al., 2017] methods build on proximal splitting algorithms by replacing the proximal or gradient descent updates with off-the-shelf Gaussian denoising operators, typically deep neural networks trained to remove Gaussian noise. The intuition behind PnP/RED is that the proximal (resp. gradient- descent) mapping of g writes as the maximum-a-posterioar (MAP) (resp. posterior mean) estimation from a Gaussian noise observation, under prior exp(-g). A remarkable property of these priors is that they are decoupled from the degradation model represented by f in the sense that one learned prior p can serve as regularizer for many inverse problems. When using deep denoisers corresponding to exact gradient-step [Hurault et al., 2021, Cohen et al., 2021] or proximal [Hurault et al., 2022] maps, RED and PnP methods become real optimization algorithms with known convergence guarantees.
However, for Poisson noise the data- fide-fideiity tenm (2) is neither smooth with a Lipschitz gradient, nor proximable for A + Id (meaning that its proximal operator cannot be computed in a closed form), which limits the use of standard splitting procedures. Bauschke et al. [2017] addressed this issue by introducing Proximal Gradient Descent (PGD) algorithm in the Bregman divergence paradigm, called Bregman Proximal Gradient (BPG). The benefit of BPG is that the smoothness condition on f for sufficient decrease of PGD is replaced by the "NoLip" condition "Lh - f is convex" for a convex potential h. For instance, Bauschke et al. [2017] show that the data-fidelity term (2) satisfies the NoLip condition for the Burg's 's entropy h(())
> log(z:).
Our primary goal is to exploit the BPG algorithm for minimizing (1) using PnP and RED priors. This requires an interpretation of the plug-in in denoiser as a Bregman proximal operator.
Prox (v) = arg min gzz) + (r,
where (1r y/ is the Bregman divergence associated with h. In Section 3, we show that by selecting a suitable noise model, other than the traditional Gaussian one, the MAP denoiser can be expressed as a Bregman proximal operator. Remarkably, the corresponding noise distribution belongs to an exponential family, which allows for a closed-forn posterior mean (MMSE) denoiser generalizing the Tweedie's 's formula. Using this interpretation, we derive RED prior tailored to the Bregman geometry. By presenting a prior compatible with the noise model, we highlight the limitation of the decoupling between prior and data-fidelity suggested in the existing PnP literature.
In order to safely use our MAP and MMSE denoisers in the BPG method, we introduce the Bregman Score denoiser, which generalizes the Gradient Step denoiser from [Hurault et al., 2021, Cohen et al., 2021]. Our denoiser provides an approximation of the log prior of the noisy distribution of images. Moreover, based on the characterization of the Bregman Proximal operator from Gribonval and Nikolova [2020], we show under mild conditions that our denoiser can be expressed as the Bregman proximal operator of an explicit nonconvex potential.
In Section 4, we use the Bregman Score Denoiser within RED and PnP methods and propose the B-RED and B-PnP algorithms. Elaborating on the results from Bolte et al. [2018] for the BPG algorithm in the nonconvex setting, we demonstrate that RED-BPG and PnP-BPG are guaranteed to converge towards stationary points of an explicit functional. We finally show in Section 5 the relevance of the proposed framework in the context of Poisson inverse problems.
2 Related Works
Poisson Inverse Problems variety of methods have been proposed to solve Poisson image restoration from the Bayesian variational formulation (I) with Poisson data-fidelit term (2), Although

---

f' is convex, this is a challenging optimization problem as f is non mon-smooth and does not have closed- form proximal operator for + Id, thus precluding the direct application of splitting algorithms such as Proximal Gradient Descent or ADMM. Moreover, we wish to regularize (1) with explicit denoising prior, which generally nonconvex. PIDAL [Figueiredo and Bioucas-Dias, 2009, 2010] and related methods INg et al., 2010, Setzer et al., 2010] solve (2) using modified versions of the alternating direction method of multipliers (ADMM), As f is proximable when = Id, the idea is to add supplementary constraint in the minimization problem and to adopt augmented Lagrangian framework. Figueiredo and Bioucas- Dias [2010] prove the convergence of their algorithm with convex regularization. However, no convergence is established for nonconvex regularization. Boulanger et al. [2018] adopt the primal-dual PDHG algorithm which also splits A from the of update. Once again, there is not convergence guarantee for the primal-dus algorithm with nonconvex regularization.
Plug-and-Play (PnP) PnP methods were successfully used for solving variety of IR tasks by including deep denoisers into different optimization algorithms, including Gradient Descent [Romano et at., 2017], Half- Quadratic-Splitinng [Zhang et al., 2017, 2021], ADMM [Ryu et al., 2019, Sun et al., 2021], and PGD [Kamilov et al., 2017, Terris et al., 2020]. Variations of PnP have been proposed to solve Poisson inverse problems. Rond et al. [2016] use PnP-ADMM and approximates, at each iteration, the non-tractable proximal operator of the data-fidelity term (2) with an inner minimization procedure. Sanghvi et al. [2022] propose a PnP version of the PIDAL algorithm which is also unrolled for image deconvolution. Theoretical convergence of PnP algorithms with deep denoisers has recently been addressed by a variety of studies [Ryu et al., 2019, Sun et al., 2021, Terris et al., 2020] (see also review in Kamilov et al. [2023]). Most of these works require non-rcalistic or sub-optimal constraints on the deep denoiser, such as nonexpansiveness. More recently, convergence was addressed by making PnP genuine optimization algorithms again. This is done by building deep denoisers as exact gradient descent operators (Gradient-Stee denoiser) [Cohen et al., 2021, Hurault al., 2021] or exact proximal operators [Hurault et al., 2022, 2023]. These PnP algorithms thus minimize an explicit functional (2) with an explicit (nonconvex) deep regularization.
Bregman optimization Bauschke et al. [2017] replace the smoothness condition of PGD by the NoLip assumption (4) with a Bregman generalization of PGD called Bregman Proximal Gradient (BPG). In the nonconvex setting, Bolte et al. [2018] prove global convergence of the algorithm to critical point of (1). The analysis however requires assumptions that are not verified by the Poisson data-fidelity term and Burg's 's entropy Bregman potential. Al-Shabili et al. [2022] considered unrolled Bregman PnP and RED algorithms, but without any theoretical convergence analysis. Moreover, the
interaction between the data fidelity, the Bregman potential, and the denoiser not explored.
Bregman denoising prior
The overall objective of this work is to efficiently solve ill-posed image restoration (IR) problems involving a data-fidelity term f verifying the NoLip assumption for some convex potential h
NoLip
There is L > 0 such that Lh ff is convex on int dom h.
(4)
PnP provides an elegant framework for solving ill-posed inverse problems with a denoising prior. However, the intuition and efficiency of PnP methods inherit from the fact that Gaussian noise is well suited for the Euclidean T2 distance, the latter naturally arising in the MAP formulation of the Gaussian denoising problem. When the Euclidean distance is replaced by a more general Bregman divergence, the noise model needs to be adapted accordingly for the prior.
In Section 3.1, first discuss the choice of the noise model associated leading to Bregman formulations of the MAP and MMSE estimators. Then we introduce in Section 3.2 the Bregman Score Denoiser that will be used to regularize the inverse problem (1).
Bregman divergence,
Notations For convenience, we assume throughout our analysis that the convex potential Ch S
R - IR U (+oc) is C² and of Legendre type (definition in Appendix A. Its convex conjugate h* then also C of Legendre type. Dr (Ir, y) denotes its associated Bregman divergence
[(x) Mly) - (Vh(y), - v) if E dom(h) (5) +00
Dp RR x int dom h + (0, +oo]: : (z, y) +
otherwise.
3

---

3.1 Bregman noise model
We consider the following observation noise model, referred to as Bregman noise',
for I, dom(h) int do p(ylz) == exp (-> y) + p(z)).
(6)
We assume that there is > 0 and a normalizing function p dom h)) + R such that the expres- sion (6) defines probability measure. For instance, for h(x) ز retrieve the Gaussian noise model with variance o². As shown Section 5, for h given by Burg's entropy. p(y(z) corresponds to multivariate Inverse Gamma (IG) distribution.
|||||||||| 12 - σz and p 0, we
Given noisy observation y int dom(h), i.e. a realization of random variable Y with conditional probability p)y) we now consider two optimal estimators of x, the MAP and the posterior mean.
Maximum-A A-Poteriio (MAP) estimator The MAP denoiser selects the mode of the a-posteriori probability distribution p(aly)) Given the prior px, it writes
AMAP arg min log p( (xly) arg min log px (x) log logp(ylz)) = Prox" (pt log )
(v). (7)
Under the Bregman noise model (6), the MAP denoiser writes as the Bregman proximal operator (see relation (3)) of - (log Px + p). This acknowledges for the fact that the introduced Bregman noise is the adequate noise model for generalizing PnP methods within the Bregman framework.
Posterior mean (MMSE) estimator The MMSE denoiser is the expected value of the posterior probability distribution and the optimal Bayes estimator for the L2 score. Note that our Bregman noise conditional probability (6) belongs to the regular exponential family of distributions
(8) >VV())).<x)> -())= >h(x)) p(x) and po(y) = exp (yh(1) - >Vhhyy), v)). is shown in [Efron, 2011] (for T = ld and generalized in [Kim and Ye, 2021] for T of Id) that the corre- sponding posterior mean estimator verifies a generalized Tweedie formula VT (y)..MMMM((V)) = -V lolo (v + V logg (y), which translates (see Appendix B for details)
p(yla) (v) exp (()) +xx))
with (τ)) = -Vhh(y)
?yMSE(y) = B[xlv] = y (v*h())) -1 V log
(9)
Note that for the Gaussian noise model, we have h(x) = 11 1/0 and (9) falls back to the more classical Tweedie formula of the Gaussian posterior mean denoiser ã = y o (i--o'V(-logaaa log (v).
Therefore, given an off-the-shelf "Bregman demoiser" By specially devised to remove Bregman noise (6) of level y, if the denoiser approximates the posterior mean B, (v) = (MME then provides an approximation of the score - V py Rs >>>M(yy) (y B,(y))
3.2 Bregman Score Denoiser
Based on previous observations, we propose to define a denoiser following the form of the MMSE (9)
B, (v) = y (V2h(y)) (v)
(10)
with B, Rd IR a nonconvex potential parametrized by neural network. Such denoiser corresponds to the Bregman generalization of the Gaussian Noise Gradient-Step denoiser proposed in [Hurault al., 2021, Cohen et al., 2021].
When B.y is trained as a denoiser for the associated Bregman noise (6) with L2 loss, approximates the optimal estimator for the Lg score, precisely the MMSE (9). Comparing (10) and (9), get Vg, ~ 2 V logpy- í.e. the score is properly approximated with an explicit conservative vector field. We refer to this denoiser as the Bregman Score denoiser.
Is the Bregman Score Denoiser Bregman proximal operator? We showed in relation (7) that the optimal Bregman MAP denoiser M AP is a Bregman proximal operator. We want to generalize this
'The Bregman divergence being non- -ommemtme,e,,, the order of the variables (x, y) in Dr. is important. Distributions of the form (6) with reverse order in Dú have been characterized in [Banerjee et al., 2005] but this analysis does not apply here.

---

property to our Bregman denoiser (10), When trained with L2 loss, the denoiser should approximate the MMSE rather than the MAP. For Gaussian noise, Gribonval [2011] re-conciliates the two views by showing that the Gaussian MMSE denoiser actually writes as an Fuclidean proximal operator.
Similarly, extending the characterization from [Gribonval and Nikolova, 2020] of Bregman proximal operators, we now prove that, under some convexity conditions, the proposed Bregman Score Denoiser (10) explicitly writes as the Bregman proximal operator of nonconvex potential.
Proposition 1 (Proof in Appendix C). Let h be C² and of Legendre type. Let 9r R U (+00) proper and differentiable and (y) dom(h) + IR" defined in (10). Assume Im(B,) C int do).hh) With V/y I IR U {+00} defined by
-h(u) + (Vh(y). v) (v) iy E int dom(() otherwise.
(v)(() {
(11)
+00
suppose that wy Vh is convex on int dom(ht). Then for dy :B - -RU(+00) defined by
0,(x) (r(1) Dn(x, for E (x) if x E <m((B,)
(12)
+00
otherwise
we have that for each y int dom(h)
B-(y) E arg min
(x, ( ox(3)
(13)
*ERR"
Remark 1. Note that the order in the Bregman divergence is important in order to fit the definition of the Bregman proximal operator (3). In this order, [Gribonval and Nikolova, 2020, Theorem 3] does not directly apply. We propose instead in Appendix a new version of their main theorem.
Remark 2. As shown in Appendix C (equation (64)), (y) can be writlen as
(B) Vib, o Vh²) Vh(y) so that the Jacobian JB, (y) = v2h(y). V o (Vh(())
The hypothesis of convexity of thy, Vh" on int dom(h") is thus equivalent to the fact that the Jacobian of B, is positive semi-definite on int dommhh.
This proposition generalizes the result of [Hurault et al., 2022, Prop. 1] to any Bregman geometry, which proves that the Gradient-Step Gaussian denoiser writes as an Euclidean proximal operator when bby V (*)( 1 || 9- (z) is convex. More generally, as exhibited for Poisson inverse problems Section 5, such convexity condition translates constraint on the deep potential Ør-
To conclude, the Bregman Score Denoiser provides, via exact gradient or proximal mapping, two distinct explicit nonconvex priors 9, and o, that can be used for subsequent PnP image restoration.
4 Plug-and-Play (PnP) image restoration with Bregman Score Denoiser
We now regularize inverse problems with the explicit prior provided by the Bregman Score De- noiser (10). Properties of the Bregman Proximal Gradient (BPG) algorithm recalled in Section 4.1. We show the convergence of our two B-RED and B-PnP algorithms Sections 4.2 and 4.3.
4.1 Bregman Proximal Gradient (BPG) algorithm
Let F and R be two proper and lower semi cemi-continuous functions with F of class C on int dom(h). Bauschke et al. [2017] propose to minimize V = + R using the following BPG algorithm:
k+1 E arg min(R(x) + (xr - h VF(x*)) + _ D ())- =ER"
(14)
Recalling the general expression of proximal operators defined in relation (3), when Vh(ax) - (()(())) E domm(""), the previous iteration can be written as (see Appendix D)
x k+1 E Prox OVh"(V/i -VF))=V)))
(15)
With formulation (15), the BPG algorithm generalizes the Proximal Gradient Descent (PGD) algo- rithm in a different geometry defined by h.
Convergence of BPG If F verifies the NoLip condition (4) for some LF > 0 and if < b one prove that the objective function & decreases along the iterates (15). Global convergence of the
5

---

iterates for nonconvex F and R. is also shown in [Bolte et al., 2018]. However, Bolte et al. [2018] take assumptions on F and that are not satisfied in the context of Poisson inverse problems. For instance, h is assumed strongly convex on the full domain I which is not satisfied by Burg's entropy. Additionally; F is assumed to have Lipschitz gradiennt on bounded subset of RR. which not true for F the Poisson data-fidelity term (2). Following the same structure of their proof, we extend in Appendix D.1 (Proposition 3 and Theorem 3) the convergence theory from [Bolte et al., 2018] with the more general Assumption 1 below, which is verified for Poisson inverse problems (Appendix E.4).
Application to the IR problem In what follows, we consider two variants of BPG for minimizing (1) with gradient updates on the data-fidelity term f. These algorithms respectively correspond to Bregman generalizations the RED Desseent (RED-GD) (ED-GD) [Romano al., 2017] and the Plug-and-Hy und-Play PGD algorithms. For the rest of this section, we consider the following assumptions Assumption 1.
(i) h: Cn + UU (+o0} is of class C2 and of Legendre-type.
(ii) / R. + IR ( (+x0} is proper, lower-bounded, coercive, of class C2 on int dom(l), with dom(h) dom and is subanalytic.
(iii) NoLip: Lyhh f is convex on int dom(h).
(iv) h is assumed strongly convex on any bounded convex subset of its domain and for all a >>> 0,
h and f are Lipschitz continuous on {a E dom(h), Vo (x) S a).
(v) 9-y given by the Bregman Score Denoiser (10) and its associated Φy obtained from Proposi- tion I are lower-bounder bounded and subanalytic.
Even though the Poisson data-fidelity term (2) convex, our convergence results also hold for more general nonconvex data-fidelity terms. Assumption (iv) generalizes [Bolte et al., 2018, Assumption DJ and allows prove global convergence of the iterates. The subanalytic assumption, defined in Appendix A.3, is sufficient condition for the Kurdyka--Lojssewwicz (KL) property to be veri- fied [Bolte et al., 2007]. The latter, also defined Appendix A.3, can be interpreted as the fact that, up to a reparameterization, the function is sharp. The KL property is widely used in nonconvex optimization [Attouch et al., 2013, Ochs et al., 2014, Bolte et al., 2018]. As detailed in Appendix A.3, the subanalytic assumption of f. 9- and Φ, allows for the Kurdyka-Lojasiews (KL) property to be verified by the objective functions f + 9- and f + Φ7.
4.2 Bregman Regularixatioon-y--eeee
(B-RED)
We first generalize the RED Gradient-Descent (RED-GD) algorithm [Romano et al., 2017] in the Bregman framework. Classically, RED-GD is simple gradient-descent algorithm applied to the functional >f + 9,, where the gradient Vg,, is assumed to be implicitly given by an image denoiser B, (parametrized by 7) via Vg.- ld B,, Instead, our Bregman Score Denoiser (10) provides an explicit regularizing potential 9, whose gradient approximates the score via the Tweedie formula (9). We propose minimize Fx,. >f + 9r on dom(() using the Bregman Gradient Descent algorithm
T++ 1 = Vh (Vh- - (VF)))((R))
(16)
which writes in a more general version as the BPG algorithm (14) with R=
1k+1 arg min Ty. f + (zx)) + (n(a,
(17)
=E3
As detailed in Appendix E.4, in the context of Poisson inverse problems. for h being Burg's entropy, FA,- = Af + 9-, verifies the NoLip condition only on bounded convex subsets of dom(h). Thus we select C' non-empty closed bounded convex subset of dom(h). For the algorithm (17) to be well-posed and to verify a sufficient decrease of (x),) the iterates need to verify I E C. We propose modify (17) as the Bregman version of Projected Gradient Descent, which corresponds to the BPG algorithm (14) with R ic, the characteristic function of the set C::
(B-RED)
A+1 E T, = arg min (ic(x) + (x - k FFn + Dn(a, (18) I
xERR
For general convergence of B-RED, need the following assumptions

---

Assumption 2.
(i) R R= with C a non-empty closed, bounded, convex and semi-algebraic subset of dom(h) such that C n int dom(h) * O.
(ii) gry has Lipschitz continuous gradient and there is Ly > O such that L 9, is convex int dom(h).
In [Bolte et al.. 2018, Bauschke et al.. 2017], the NoLip constant L needs to be known to set the stepsize the BPG algorithm as =L In practice, the NoLip constant which depends on f. 9y and C' is either unknown or over-estimated. In order to avoid small stepsize, we adapt the backtracking strategy of [Beck, 2017, Chapter 10] automatically adjust the stepsize while keeping convergence guarantees. Given y E (0. 1), 1) E [0, 1) and an initial stepsize To > O, the following backtracking update rule on is applied each iteration k:
while (11) Fx, (T, (xx)) (T,
(TA,, T - TT. T
(19)
Using the general nonconvex convergence analysis of BPG realized in Appendix D.I, we can show sufficient decrease of the objective and convergence of the iterates of B-RED.
Theorem 1 (Proof in Appendix D.2). Under Assumption and Assumption 2, the iterates (zx) given by the B- B-RED algorithm (18) with the backtracking procedure (19) decrease Fx, and converge to critical point of Vp = ic + Fx,, with rate minogkgK Dr (n(x ,x*) = O(1KK).
4.3 Bregman Plug-and-Play (B-PnP)
We now consider the equivalent of PnP Proximal Gradient Descent algorithm in the Breg- man framework. Given a denoiser By with Im(B,) C dom(h) and A > 0 such that Im 1 - xV.) S dom(Vh"), it writes
(B-PnP)
1++ (Vh (Vh-)V(((a))
(20)
We use again as B. the Bregman Score Denoiser (10). With Q/ defined from 9, as in (11) and assuming that s/y Vh is convex on int dom(h*), Proposition I states that the Bregman Score denoiser By is the Bregman proximal operator of some nonconvex potential Φ, verifying (12). The algorithm B-PnP (20) then becomes x++1 E Prox ohh" (Vh (xx)) which writes as a Bregman Proximal Gradient algorithm, with stepsize T = 1,
x A+1 E arg min ( &b (x) + (x - (V + Da(x,=")).
(21)
2E3"
With Proposition 1, we have (v) E Prox d. (v) iee a proximal step on Φ, with stepsize 1. We are thus forced to keep a fixed stepsize T = in the BPG algorithm (21) and no backtracking is possible. Using Appendix D.
show that B-PnP converges towards stationary point of Xf + ÓT.
Theorem 2 (Proof in Appendix D.3). Assume Assumption / and bby o Vh strictly convex on domm(h") Im(B.) S dom(h) and XL5 1 (with L, specified in Assumption 1), the iterates In given by the B-PnP algorithm (20) decrease Af + Φ, and converge to a critical point of >f + Φ. with rate mingeks (x 1 O(1/K).
int dom(h"). Then for Im(Vh. XV
Remark 3. The condition Im Im(B.)) C dom(h) and the required convexity of th,, Vh come from Proposition while the condition Im Vh xV S. dom( Vh") allows the algorithum B-PnP (20) to be well well-posed. These assumptions will be discussed with more details in the context of Poisson image restoration in Section 5.
Application to Poisson inverse problems
We consider ill-posed inverse problems involving the Poisson data-fidelity term f introduced in (2).
The Euclidean geometry (i.e. h(x) = |x||| does not suit for such f, as it does not have a Lipschitz
gradient. In [Bauschke et al., 2017, Lemma 7J. it is shown that an adequate Bregman potential h (in
the sense that there exists Ly such that Lyh - is convex) for (2) is the Burg's 's entropy
1
h(x)) = Σ log(x;),
(22)
7

---

for which dom(h) = + and Lyh-f is convex on int dom (h) = R for L, 2 1lll For further computation, note that the Burg's entropy (22) satisfies Vh(x:) = Vh*(z) - and <² = # The Bregman score denoiser associated to the Burg's s entropy is presented in Section 5.1. The corre- sponding Bregman RED and PnP algorithms are applied to Poisson Image deblurring in Section 5.2.
5.1 Bregman Score Denoiser with Burg's entropy
We now specify the study of Section 3 to the case of the Burg's entropy (22). In this case, the Bregman noise model (6) writes (see Appendix E.1 for detailed calculus)
p(yl)) exp(p(x) m))
II i=1
Ti
exp
(-7
(23)
yo
Yt
For 1, this is product of Inverse Gamma (IG(a, B)) distributions with parameter 3, = 224 and a; = - 1. This noise model has mean (for y > 2) - x and variance (for > 3) (x-2) x². In particular, for large y, the noise becomes centered x with signal-dependent variance a²² //- Furthermore, using Burg's entropy (22), the optimal posterior mean (9) and the Bregman Score
x2
Denoiser (10) respectively write, for y E ++
1 -1² (VV(-oggPY) log
@MMSE(V) (v) = y
(24)
?
(v) = v
(25)
Denoising in practice As Hurault al. [2021]. we chose to parametrize the deep potential 9r as
1
(r(x)
2 a - (z)1
(26)
where N₂ is the deep convolutional neural network architecture DRUNet [Zhang et al., 2021] that contains Softplus activations and takes the noise level y as input. Vg. is computed with automatic differentiation. We train By to denoise images corrupted with random Inverse Gamma noise of level 7, sampled from clean images via p(yla) = -II= IG ZG B:) (1zg). To sample y; ~ IG( 3r),
we sample 2 ~ G(a,, 3,) and take = 11/2. 24 Denoting as p the distribution of database of clean images, training is performed with the 2 loss
((7)) = Ex-n-y-T,, (3) (1) ==1
(27)
Denoising performance We evaluate the performance of the proposed Bregman Score DRUNet
(B-DRUNet) denoiser (26). It is trained with the loss (27), with 1/7 uniformly sampled in (0, 0.1). More details on the architecture and the training can be found in Appendix E.2. We compare the performance of B-DRUNet (26) with the same network DRUNet directly trained to denoise inverse Gamma noise with L2 loss. Qualitative and quantitative results presented in Figure and Table show that the Bregman Score Denoiser (B-DRUNet), although constrained to be written as (25) with a conservative vector field Vgy, performs on par with the unconstrained denoiser (DRUNet).
Y
10
25
50
100
200
DRUNet
28.42
30.91
32.80
34.76
36.79
B-DRUNet
28.38
30.88
32.76
34.74
36.71
Table 1: Average denoising PSNR performance of Inverse Gamma noise denoisers B-DRUNet and DRUNet on 256 x 256 center-cropped images from the CBSD68 dataset, for various noise levels y-
Bregman Proximal Operator Considering the Burg's 's entropy (22) in Proposition I we get that, if "h : 2 1/ o Vh* (a) = (-1) is convex on dom(h") R -, the Bregman Score Denoiser (25) satisfies (v) = Prax for the nonconvex potential ,/ (y) = Σ log(us) - 94 (v) 1. The
convexity of sh. o Vh can be verified using the following characterization, which translates to a condition on the deep potential 9,, (see Appendix E.2 for details)
Vx E dom((h)) Vd R", (V2n 0
(28)
AA R" YdER", R", ( (V'''g,(o)d < Σ (1r Vo
(29)
8

---

(a) Clean
(b) Noisy (17.8 85 dB)
(c) DRUNet (32.56 dB)
(d) B-DRUNet (32.56 dB)
Figure I: Denoising of a 256 x 256 image corrupted with Inverse Gamma noise of level 25.
is difficult to constrain gy to satisfy the above condition. Appendix E.2, we proceed to a detailed experimental validation of this assumption. We empirically observe images that the trained deep potential By satisfies the above inequality. This suggests that the convexity condition for Proposition true least locally, condition holds for the theoretical MMSE denoiser (24). As explained below, our denoiser is trained by minimizing the T.2 cost for which the optimal Bayes estimator is the MMSE, The MMSE is then a theoretical denoiser that our denoiser tries to approximate. This explains why our denoiser, after training, naturally satisfies this condition without necessitating supplementary constraints.
the image manifold. Indeed, prove in Appendix E.3 that the convexity
5.2 Bregman Plug-and-Piag for Poisson Image Deblurring
We now derive the explicit B-RED and B-PnP algorithms in the context Poisson image restoration. Choosing C [0, Rt" for some R >> O, the B-RED and B-PnP algorithms (18) and (20) write
(B-RED)
k+1 2
arg min
v Fx
+
T log
I I
[0, R}} Isisn n (30)
2
if 1tr2
S R
V (ott)e
isisn (31)
R else
x
(B-PnP) 2 k+1 B.
(32)
I 1++T*\S(0
Verification of the assumptions of Theorems 1 and 2 for the convergence of both algorithms is discussed in Appendix E.4.
Poisson deblurring Equipped with the Bregman Score Denoiser, we now investigate the practical performance of the boug-and-play B-RED and B-PnP algorithms for image deblurring with Poisson noise. In this context, the degradation operator A is convolution with blur kernel. We verify the efficiency both algorithms over variety of blur kernels (real-world camera shake, uniform and Gaussian). The hyper-parameters 7, À are optimized for each algorithm and for each noise level by grid search and are given in Appendix F.5. practice, we first initialize the algorithm with 100 steps with large 7 and y so as to quickly initialize the algorithm closer to the right stationary point.
Note that the constraint AL,- 1 for convergence of the B-PnP algorithm (Theorem 1) may not be respected. The global NoLip constant L, can indeed be locally very lose. As explained in the Appendix (32), can adopt backtracking-like strategy on the regularization parameter À to ensure convergence. Nevertheless. with the proposed default value of X, this backtracking algorithm was never activated over the variety of blur kernels and noise levels experimented.
We show in Figures 2 and that both B-PnP and B-RED algorithms provide good visual reconstruction. Moreover, observe that, in practice, both algorithms satisfy the sufficient decrease property of the objective function as well as the convergence the iterates. Additional quantitative performance and comparisons with other Poisson deblurring methods are given in Appendix E.4.
6 Conclusion
In this paper, we derive complete extension of the plug-and-play framework in the general Bregman paradigm for non-smooth image inverse problems. Given convex potential h adapted to the geometry

---

(a) Clean
(b) Observed (14.91dB)
(c) B-RED (23.01dB)
(d) B-PnP (22.96dB)
6-E0
5 RED
20000
10t4
PM
-22
2400
2000
10-
-9.9
(e) (Af + (TL)
(I) (Af + (ZL)
100 150 200
B-RED
B-PnP
(g) |2541 2:1
Figure 2: Deblurring from the indicated motion kernel and Poisson noise with a = 40.
(a) Clean
(b) Observed (16.21dB)
(c) B-RED (24.50BB)
(d) B-PnP (24.51dB)
D-ED
B RED
8 PxP
10
10-
135
(e) (Af + (12) B-RED
(1) (Af (xL) B-PnP
(g) 20:4 2211
Figure 3: Deblurring from the indicated Gaussian blur kernel and Poisson noise with a = 60.
of the problem, we propose a new deep denoiser, parametrized by h, which provably writes as the Bregman proximal operator of nonconvex potential. We argue that this denoiser should be trained on a particular noise model, called Bregman noise, that also depends on h. By plugging this denoiser in the BPG algorithm, we propose two new plug-and-play algorithms, called B-PnP and B-RED, and show that both algorithms converge to stationary points of explicit nonconvex functionals. We apply this framework to Poisson image inverse problem. Experiments on image deblurring illustrate numerically the convergence and the efficiency of the approach.
The central significance of work stems from its theoretical study but we recognize certain limits within our experimental results. First, when applied to deblurring with Poisson noise, our proposed algorithms do not outperform existing methods in terms of PSNR. Second, while we prove that B-RED is convergent without restriction, the convergence of B-PnP depends on a specific convexity condition. Despite being confirmed with experiments and having robust theoretical foundations, this assumption could potentially be not verified when applied to non-natural images that significantly differ from those in the training dataset. Finally, due to the nonconvex nature of our proposed prior, the practical performance of the algorithms can be sensitive to their initialization.
Acknowledgements This work was funded by the French ministry of research through a CDSN grant of ENS Paris-Saclay. It has also been carried out with financial support from the French Research Agency through the PostProdI.FAA and Mistic projects (ANR-19-CF2330007--0) and ANR- 19- (9-CE40-005). It has also been supported by the NSF CAREER award under grant CCF-2043134.
10

---

References
Abdullah H Al-Shabili, Xiaojian Xu, Ivan Selesnick, and Ulugbek Kamilov. Bregman plug-and-plaa priors. In 2022 IEEE International Conference Image Processing (ICIP), pages 241-245. IEEE, 2022.
Hédy Attouch, Jéröme Bolte, Patrick Redont, and Antoine Soubeyran. Proximal alternating mini- mization and projection methods for nonconvex problems: An approach based on the kurdyka- lojasiewicz inequality. Mathematics of operations research, 35(2):438 4474457. 2010.
Hedy Attouch, Jérôme Bolte, and Benar Fux Svaiter. Convergence of descent methods for semi- algebraic and tame problems: proximal algorithms, forward-backe backward splitting, and regularized gauss- seidel methods. Mathematical Programming. 137(1-22:91-122 2013.
Arindam Banerjee, Srujana Merugu, Inderjit S Dhillon, Joydeep Ghosh, and John Lafferty. Clustering with bregman divergences. Journal of machine learning research, 6(10), 2005.
Heinz H Bauschke, Jérôme Bolte, and Marc Teboulle. A descent lemma beyond lipschitz gradient continuity: first-order methods revisited and applications. Mathematics of Operations Research, 42(2):330-348, 2017.
Amir Beck. First-order methods in optimization. SIAM, 2017.
Mario Bertero, Patrizia Boccacci, Gabriele Desiderà, and Giuseppe Vicidomini. Image deblurring with poisson data: from cells to galaxies. Inverse Problems, 25(12): 123006, 2009.
Jérôme Bolte, Aris Danillidis, and Adrian Lewis. The tojasiewicz inequality for nonsmooth subana- lytic functions with applications to subgradient dynamical systems. SIAM Journal on Optimization, 1754):1005-1223, 2007.
Jérôme Bolte, Aris Daniilidis, Olivier Ley, and Laurent Mazet. Characterizations of lojasiewicz inequalities: subgradient flows, talweg, convexity. Transactions of the American Mathematical Society, 362(6):3319 2010.
Jérôme Bolte, Shoham Sabach, Marc Teboulle, and Yakov Vaisbourd. First order methods beyond convexity and lipschitz gradient continuity with applications quadratic inverse problems. SIAM Journal on Optimization, 28(3):1131211 -1551,, 2018.
Jérôme Boulanger, Nelly Pustelnik, Laurent Condat, Lucie Sengmanivong, and Tristan Piolot. Nonsmooth convex optimization for structured illumination microscopy image reconstruction. Inverse problems, 34(9):095004, 2018
Regev Cohen, Yochai Blau, Daniel Freedman, and Ehud Rivlin. It has potential: Gradient-driven denoisers for convergent solutions to inverse problems. Advances in Neural Information Processing Systems, 34, 2021.
Patrick Combettes and Jean San-Christoophe Pesquet. Proximal splitting methods in signal processing. In Print Algorithms for Inverse Problems in Science and Engineering, pages 185-212. Springer, 2011.
Michel Coste. An introduction to semialgebraic geometry, 2000.
Bradley Efron. Tweedie's s formula and selection bias. Journal of the American Statistical Association, 106(496):1602-1611 2011.
Mario AT Figuciredo and Jose Bioucas-Dias. Deconvolution of poissonian images using variable splitting and augmented lagrangian optimization. In 2009 IEEE/SP 15th Workshop on Statistical Signal Processing, pages 733-736. IEEE, 2009.
Mário AT Figueiredo and José M Bioucas- Dias. Restoration of poissonian images using alternating direction optimization. IEEE transactions on Image Processing, 19(12)331333115, 2010.
Rémi Gribonval. Should penalized least squares regression be interpreted as maximum posteriori estimation? IEEE Transactions Signal Processing, 2905)22405- 2011.

---

Rémi Gribonval and Mila Nikolova. A characterization of proximity operators. Journal of Mathemat- ical Imaging and Vision, 62(6):773-789, 2020.
Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis. Gradient step denoiser for convergent plug-and-pla arXiv preprint arXiv:2110.03220, 2021.
Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis. Proximal denoiser for convergent plug- and-play optimization with nonconvex regularization. In International Conference on Machine Learning, pages 9483-9505. PMLR, 2022.
Samuel Hurault, Antomin Chambolle, Arthur Leclaire, and Nicolas Papadakis. A relaxed proximal gradient descent algorithm for convergent plug-and-play with proximal denoiser. Scale Space and Variational Methods in Computer Vision: 9th International Conference, SSVM 2023, Santa Margherita di Pula, Itly, May 21-25, 2023, Proceedings, pages 379-392. Springer, 2023.
U. S. Kamilov, C. A. Bouman, T. Buzzard, and Wohlberg, Plug-and-play methods for integrating physical and learned models in computational imaging. IEEE Signal Process. Mag., 40855988 January 2023.
Ulugbek S. Kamilov, Hassan Mansour, and Brendt Wohlberg. A plug-and-play priors approach for solving nonlinear imaging inverse problems. IEEE Signal Process. Lett., 24(12):182-1187, December 2017.
Kwanyoung Kim and Jong Chul Ye. Noise2score: tweedie's 's approach to self-supervised image denoising without clean images. Advances in Neural Information Processing Systems, 34:864-874, 2021.
Steven Krantz and Harold R Parks. A primer of real analytic functions. Springer Science Business Media, 2002.
Anat Levin, Yair Weiss, Fredo Durand, and William T Freeman. Understanding and evaluating blind deconvolution algorithms. In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 1964-1971. IEEE, 2009.
Michael K Ng, Pierre Weiss, and Xiaoming Yuan. Solving constrained total-variation image restora- tion and reconstruction problems via alternating direction methods. SIAM journal on Scientific Computing, 32(5):2710-2777 2010.
Peter Ochs, Yunjin Chen, Thomas Brox, and Thomas Pock. ipiano: Inertial proximal algorithm for nonconvex optimization. SIAM Journal on Imaging Sciences, 7(2):1388-141) 2014.
Tyrrell Rockafellar. Convex analysis, volume 11, Princcton university press, 1997.
Yaniv Romano, Michael Elad, and Peyman Milanfar. The little engine that could: Regularization by denoising (red). SIAM Journal on Imaging Sciences,
2017.
Arie Rond, Raja Giryes, and Michael Elad. Poisson inverse problems by the plug-and-plan scheme. Journal of Visual Communication and Image Representation, 41:96-108, 2016.
Ernest Ryu, Jialin Liu, Sicheng Wang, Xiaohan Chen, Zhangyang Wang, and Wotao Yin. Plug-and- play methods provably converge with properly trained denoisers. In International Conference on Machine Learning, pages 5546-5557. PMLR, 2019.
Yash Sanghvi, Abhiram Gnanasambandam, and Stanley Chan. Photon limited non-blind deblurring using algorithm unrolling. IEFE Transactions on Computational Imaging, 8:851-864, 2022.
Simon Setzer, Gabriele Steidl, and Tanja Teuber. Deblurring poissonian images by split bregman techniques. Journal of Visual Communication and Image Representation, 21(3):193-19,, 2010.
Masahiro Shiota. Geometry of subanalytic and semialgebraic sets, volume 150. Springer Science Business Media, 2012.
Yu Sun, Zihui Wu, Xiaojian Xu, Brendt Wohlberg, and Ulugbek S Kamilov. Scalable plug-and-play admm with convergence guarantees. IEFE Transactions on Computational Imaging, 7:849-86. 2021.
12

---

Matthieu Terris, Audrey Repetti, -annCCristtophe Pesquet, and Yves Wiaux. Building firmly nonex- pansive convolutional neural networks. In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 8658 8662. IEFE, 2020.
Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg. Plug-and-play priors for model based reconstruction. In 2013 IEEE Global Conference on Signal and Information Processing, pages 945-948. IEEE, 2013.
Kai Zhang, Wangmeng Zuo, Shuhang Gu, and Lei Zhang. Learning deep emn denoiser prior for image restoration. In Proceedings of the IEEE conference on computer vision and paltern recognition, pages 3929-3938, 2017.
Kai Zhang, Yawei Li, Wangmeng Zuo, Lei Zhang, Luc Van Gool, and Radu Timofte. Plug-and play image restoration with deep denoiser prior. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021.
13

---

A Definitions
A.1 Legendre functions, Bregman divergence
Definition 1 (Legendre function, Rockafellar [1997]). Let h : C R" + R +RU(+O) be a proper lower coniiconniuuuss convex convex function. is called:
(i) essentially smooth, if h is differentiable on int dom(h), with moreover [[h(<>>) - Do for every sequence {a" Juets of int dom(h) converging towards a boundary point of dom(h).
(ii) of Legendre type if h is essentially smooth and strictly convex on int (l)
We also recall the following properties of Legendre functions that are used without justification in our analysis (see [Rockafellar, 1997, Section 26] for more details):
h is of Legendre type and only if its convex conjugate h* is of Legendre type.
For h of Legendre type, (oom(Vh) int dom(h), Vh is bijection from ind dom(h) to int dom(m") (h²) and (Vh)-1= Vh.
Examples
We are interested in the two following Legendre functions
Le potential: : h(x)- = ||1=11 C = dom(h) (h) = dom(h") = R2" Dx 3) = 2110z= y/1
. Burg's entropy h(x)
Σ - log(2;), C dom(h) -R dom(h²) = -2___ and
D,(x, Ei= yi
log
22
1.
A.2 Nonconvex subdifferential
Following Attouch et al. [2013], we use as notion of subdifferential of a proper, nonconvex function & the limiting subdifferential
(z) {w RR" Eez - x, ((xx)) ((x), w. E )0()))
(33)
with d& the Fréchet subdifferential of Φ. We have
{w 3x + x, ((x)) - D(x), We + www E (0(x)) e @2(x).
(34)
A.3 Kurdyka-L.ojsievier (KL) property and subanalytic functions.
Definition 2 (Kurdyka-Lojasiewiic) (KL) [Attouch et al., 2013]), A funccion J Rn - R +00
said to have the hardylk--oojsiiwiic. property at x² dom(f) if there exists n E (0, +00),
neighborhood U of x" and continuous concave function 1/ [0. 7)) - R. such that 1/( = 0, v)
C. on (0, 11), v/ > on (0,1) () and Vx E U [f(x)(()))))))) < f(x" + n). the Kundyka-Lojasiene
inequality holds:
('((x)) - f(x='( ddistt(,,,, 0f(x)) 21
(35)
Proper lower semicontinuous functions which satisfy the point of dom((f) are called KL. functiiss
inequality at each
This condition can be interpreted the fact that, up to reparameterization, the function is locally sharp around its critical points i.e. we can bound its subgradients away from O. For more details and interpretations, we refer to [Attouch et al., 2010] and [Bolte et al., 2010]. As shown in Bolte et al. [2007]. a large class of functions that have the KL-property is given by continuous subanalytic functions.
Definition 3 (Subanalytic [Bolte al., 2007]).
A subset S of IR is semianalytic set if each point of Rd admits a neighborhood V for which there exists finite number of real analytic functions figj Rig : B - R. such that
Sn L snv=UP= no EI fu, 9i. < 0}
(36)
14

---

A subset S of IR" is a subanalytic set if each point of R" admit a neighborhood V for which
snv= R", (I. E EU}
(37)
where U is a bounded semianalytic subset of R R with p 2
A function f: _n + U {+oo} (resp. f R RT is called subanalytic if its graph
{(x,y) R x y [()) (resp. {(x, v) E R Rmm y f()))) is subanalytic
subset of RD x (resp. R" x J.
Subanalytic functionn encompass both analytic functions and semialgebraic [Coste, 2000] functions [Shiota, 2012]. Stability results of subanalytic functions by sum, product, and composition are given in [Shiota, 2012, Section 1.2].
B More details the Bregman Denoising Prior
We detail here the calculations realized in Section 3.
We first verify that the Bregman noise conditional probability belongs to the regular ex- ponential family of distributions. po(y) = exp (yh(y) - ((Vh(y) v))
Indeed for T(y) >Vh(v), v(xx) = (h(x) - p(z) and
x)
Po(y) exp ((r, T(y) +xx))
exp (7h(1) (Vh(y)), exp (7(2, Vh(y)) - -h(x) + p(x))
(38)
exp (->(1(z) - h(w) - -Vh((), u\) + p(x))
exp (=yDD> y) + p(x))
which corresponds to the Bregman noise conditional probability introduced Equation (6).
The Tweedie formula of the posterior mean is then [Kim and Ye, 2021]
VT(y).&MMME(( - log (ol)o(()> + V Py
(39)
Using
V Vlogpolu) = -Vh(v) - - (oA))) - (40)
we get
(MMS)))) = (Vkk) Volgpr)))
(41)
As is strictly convex, hyy is invertible and
{MMSE()) = y + 1 (V2h())) V py (v).
(42)
C Proof of Proposition
characterization of Bregman proximal operators We first derive in Proposition 2 an extension of [Gribonval and Nikolova, 2020, Corollary 5 a)] when h is strictly convex and thus V/ is invertible. We will then see that Proposition is a direct application of this result.
Proposition 2. Let h of Legendre type on IR.. Let Ç int dom(h) - R... Suppose Im(() dom(h). The following properties are equivalent:
(i) There is ob RR" + RR {+00} such that Im(() C do and for each y E int dom
((y) E arg ming Dx(x, + d(x)). #ER"
(43)
(ii) There is Ls.c s Ra U (+0) proper convex on int dom(hi) such that <(Vh' (2)) for each int dom(h").
15

---

When (i) holds, 2b can be chosen given ø with
(c(Vh(())) -
if a dom((")))
((2)
(44)
+00
otherwise.
When (ii) holds, & can be chosen given vp with
(C(W). h(((#)) (Vh(yy)) for (x) if E Im(()
((x)
(45)
+00
otherwise.
Remark 4. |Gribonval and Nikolova, 2020, Corollary 5 a)] with y = int dom(h) states that (i) is equivalent to
(ii) There is Ls.c g R" + R U {+00} convex such that (x(s- E 0g(x) for each x E Lm(
However (ii) and (iii) are not equivalent. We show here that (ūi implies (iii) but the converse not true. Let dp convex defined from (ii). Thanks to Segendre-Fencheee identity, we have that
Vz int dom (h*), <(Vh" (=)) E 0v((2)
(46)
22 E ind dom (h"), E 0ui (s(Vh"())))
(47)
Vy E int dom (h), Vh(y) E Oip (C(y))
(48)
Vz E Im((), (x())((z)) E abp (x).
(49)
Therefore (ii) implies (iti) with g = ... However, the last lign is just an implication.
Proof. We follow the same order of arguments from [Gribonval and Nikolova, 2020]. We first prove a general result reminiscent to [Gribonval and Nikolova, 2020, Theorem 3] for a general form of divergence function and then apply this result to Bregman divergences.
R - RU A y + 2 bijection from y Lemma 1. Let a : - Ru R.
to Z (with C J. Consider Ç:) R". Let D(z, v) == a(y) (x. A(y)) b(z). Suppose Im({) dom(b). The following properties are equivalent:
(i) There is R + RU (+00) such that Im(() C dom( and each E y
((y) arg min o(x)).
(50)
#ER"
(ii) There is Lx. :b Rp + RU (+00} proper convex such that ((A-1(2)) E Orb(2) each z Z.
When (i) holds, v can be chosen given with
(12) (s),
((((A)))))
if Z otherwise.
(51)
100
When (ii) holds, & can be chosen given vp with
(xx) = { (C()),A()) - b((())) - (A(()) for E (x) if E Im({)
(52)
t00
otherwise.
Before proving the Lemma, we can directly see that Proposition 2 the specialization of Lemma with Bregman divergences. Given h of Legendre-type, the divergence D(x,y) = a(y) (Ir, A(y)) + b(x)) becomes the Bregman divergence Da y) defined in (5) for y int dom(h)) 2 = ind domih")) A(y) = Vh(u), a(y) = (Vh(y), h(y) and b(x) = M(x) if dom(h) and + oo otherwise. For h of Legendre-type, Vh is bijection from int dom(h) to
int dom(h") and (Vh)-1= V
0
16

---

Proof. We now prove Lemma 1. We follow the same arguments as the proof of [Gribonval and Nikolova, 2020, Theorem 3.c)]. (i) (ii)
Define
(s(A- (((((())))))(((AA((()))))) ((C((-
if E 2 else
((z) = d=z)=
(53)
+00
As we assume Im(() C dom(b) and lm(() C dom(t) (c), have p R" + IR U [+oo} and dom(p) = Z. Let E 2 and y = A (z). From (.. ((y) is the global minimizer of x + D(I,y) + +(x) or of + - (x, A(v)) b(z) + 6(x) and V2' E 2 A (z').
p(z') p(z) ((
(C()), +b((())) + +(((v)) (54)
-((()),A0') 103-A11 A(v)) ((t'))(()))-- (S(o') - M(((')) *((()))
- -((()) +b((6)))+((()))))
2 (C(o), -oo0)-A(00) ((v)) (=)A- a' -)--
By definition of the subdifferential, this shows that
(z) E 0p(z).
(55)
Let f the lower convex envelope of p (pointwise supremum of all the convex Iss.c functions below p). is proper convex I.s.c. and V2 2, 0p(z) + W. By [Gribonval and Nikolova, 2020, Proposition 3], V2 E Z, p(z) = p(2) and 0p(z)= = ())... Thus for ai = م. we get (ii).
(ii) > (i)
Define 1) y + R by
<vy) = (C(v), A(v)) - (A())).
(56)
The previous definition is valid because by (ii), Im(A) Z c dom(@v) C. dom( and therefore (A(()) < +00. By (ii), V2, E Z.
(1=) - <1x') 2 (C(A (=))- -=== =²)
(57)
which gives Vy, s EV.
(((())) *(A()) 2 (C(1)) A(V'))
(58)
This yields
7lo') - 73)) = (C()). - -((()))))))))))) - (C(v), + (AA())
(59)
2 (c(r/)(((w)) -
A(y)).
We define RR" + RU { { obeying dom = lm(() with
0(x) = { n(y) for y E (x) i Im(() otherwise.
(60)
400
For y, y E (z). as ((y) ((v), we have by (59) nlv') - n(v) 2 0 and n(y) - n(v') 2 0 and
thus 7/5/) = 7(y). The definition of thus independent of the choice of y E (x).
For x' E Tm((), x'= = ((v'). Using the previous inequality with 7, we get
0(x') - e(((v)) = B(((v')) o(((v))
nTt//) - n(y)
(61)
2 (C(1/) - C(v), A(r))
(a' ((1))
that is to say, Va E Im(()
0(x') A(w)) > O(((v)) A(v))
(62)
Given the definition of 0, this is also true for x' & Im(().
17

---

We set Ø = 0 - b. As Im(() C dom(b). b(((v)) < too and Φ: RR - R U [+oo). Adding a(y) on both sides, we get
Va' E R", b(x') + o(x') - (a', A(y)) + a 2 + + ((()))))(((()))))))))))))))))))) aly) Vo' E RR, d(x') + D 3) 2 + + D(C(v),
(63)
((w) E arg min ((x) + (
x
As this is true for all E V we get the desired result.
Eventually, Proposition is direct application Proposition with B, : R² RR" defined on int dh) by
(v) ( Vh*) Vh(y). The function By verifies Va E int (h").
(64)
(Vh" = V( o V+")=)
(65)
and ii = v Vh is assumed convex on int domm(h"). From Proposition 2, we get that there is by R >> TR U (+00) such that for each y E int dom(h)
(v) arg mig Dr (x, y) + o (x))
(66)
Moreover, for y int (h),
(y) V Vh") Vh(y)
(67)
(V Vh Vh(y)
(68)
= (Vh(y)))))
(69)
As h is assumed strictly convex on int dom(h), for y E int dhm the Hessian of h, denoted as hyy) is invertible. By differentiating
VI((Vh()) -
(70)
get
<2h'(hh(m) = (V2)()))-
(71)
so that
By = (v (v))
(72)
With the definition (11), we directly get
(2) (V2h(y))- = -V
(73)
Finally, Proposition also indicates that Φy be chosen given 8/y with
(x(x) = {
(B, (v), Vh(y)) - h(B, bh VA (Vh(y)) for y E B7 i E Im(B.)
+00
otherwise.
(74)
This gives Vy E int dom
(x), (B), Vh(y)) -4 VV
(B,(() - y. Vh(y)) h(B, + My) + (y, Vh(y)) h(y) e
(75)
-B,, + (y, Vh(n)) h(() -())-11)))))) (v)
(v) (1)
D
The Bregman Proximal Gradient (BPG) algorithm
D.I Convergence analysis of the nonconvex BPG algorithm
We study in this section the convergence of the BPG algorithm
k+1 E TTxx) arg (in( + (I k VF(x*)) + I Dnle, .
(76)
for minimizing V = F + R. with nonconvex functions F and/or R.
For the rest of the section, take the following general assumptions.
18

---

Assumption 3.
(i) h IR" + R is of Legendre-type.
(ii) : R + is proper, C on int dom(h), with dom(h) C dom dom(F).
(iii) R. : Rd - is proper, lower semi-continuous with dom R. n int dom(h) + 0.
(iv) v = F + R is lower-bounded, coercive and verifies the Rundyka-Lojcssewiiz (KL) property (defined in Appendix A.3).
(v) For x E int int dom(h), T (x) nonempty and included in int dom(h).
Note that, since R. is nonconvex, the mapping T, is not in general single-valued.
Assumption (v) required for the algorithm to be well-posed. As shown in [Bauschke et al., 2017, Bolte et al., 2018], one sufficient condition for (x) & 0 is the supercoercivity of function h + XR for all À 0, that is | true when dom(h) is open (which is the case for Burg's entropy for example).
+()+))(()
+00. As T, (xx) dom()), T, (+) int dnt
|||||
The convergence of the BPG algorithm in the nonconvex setting is studied by the authors of [Bolte al., 2018]. Under the main assumption that Lh is convex on inf dom(h), they show first the sufficient decrease property (and thus convergence) of the function values, and second, global convergence the iterates. However, as we will develop in Appendix D.3, Lh F is not convex on the full domain of int dom(h) but only on the compact subset (R).
One can verify that all the iterates (76) belong to the convex
T (a)) Conv(dom R) O int dom (h),
(77)
where Conv E)) stands for the convex envelope of E. We argue that it is enough assume Lh F convex on this convex subset with the following assumption.
Assumption 4. There is L > O such that, Lh F is convex on Conv(dom R) n inl dom(h).
We can now prove result similar than [Bolte et al., 2018, Proposition 4.1].
Proposition 3. Under Assumptions 3 and 4, let (x* kKE be a sequence generated by (18) with 0 < TL < 1. Then the following properties hold
(i) (V(x*)) is non-increusing and converges.
(ii) Σg (a(kkt1 x²) < co and minosks (n(x+) ,xa) = O(1/K).
Proof. We adapt here the proof from [Bolte al., 2018] to the case where Lh F is not globally convex but only convex on the convex subset Conv(dom R) n int dom(h).
Sufficient decrease property We first show that the sufficient decrease property of v/ (zx) holds. This is true because the following characterisation of C convex functions holds on convex subset. We recall this classical proof for the sake of completeness.
Lemma Let f X T be of class C2,, then f locally convex on convex subset of X df if and only if C, Dy(x, v) f(z) fly) (Vi(y), v) 0.
Proof. Let E C, for all E (0, fly x) C and by convexity of f
fly + 1(x v)) (v) ((((x)) {()))
(78)
i.e.
f(y + t(x w)) f(w) < f(x) -11)
(79)
1
and
fly + t(x - v)) f(w)
(V)()), v) lim
< f(x) - f(v).
(80)
6-30+
19

---

+ Let L, C, (0, and 2 = y E C. We have
f())+ (Vf(2), y 2) [(x) 2 f(=)+ (=) 2)-
(81)
(82)
Combining both equations gives
s (00-
z))
= =[()) + (1
(83)
[(tx ()-)
Therefore, we have Lh F convex on C, if and only if, Vx,y E C. DLN-F(Z.V) 2 0, i.e. Dp(r,)) y) S (z, y).
The rest of the proof identical to the of [Bolte al., 2018] and recall it here. Given the optimality conditions in (76), all the iterates 2%. E C satisfy
R(x*+1> + (x 2t, VF(xk)) + Du( s R(x*)
(84)
using DF(x,)) v) S LDj (2, v), get
(R(x + (() - (R(a*) + F(a*)) 1 - -1 Di -n ( )++)+ LDD
(85)
= + ) (2 xx) 00
which, together with the fact that V is lower bounded, proves (i). Summing the previous inequality from 0 to K gives
K 0 E
(a k+1 xx)
T I-TL"
T (Zz)) S 1- TL
(T (w inf v(x)) < +o0. (86)
ZEC
K=
Thus (Di(x ))x is summable and converges to 0 when lt + to0. Finally
1
min Da Da(a k+ xty) 0skSK
S K
Datatt1
,ts K K+11 7L
(L) (Voz) - inf (x)). (87)
ZEC
k=
To prove global convergence of the iterates upon the Kurdyka--Lojaiewiiez (KL) property, [Bolte et al., 2018, Theorem 4.1] is based on the hypotheses (a) dom(h) = R" and h is strongly convex on I2n and (b) Vh and F are Lipschitz continuous on any bounded subset of R". These assumptions are clearly not verified for h being the Burg's entropy (22) or F the Poisson data-fidelity term (2). Indeed, in that case, dom(h) = I and is strongly convex only on bounded sets. Moreover F and h are not Lipschitz continuous near 0.
However, thanks the proven decrease of the iterates and as is assumed coercive, the iterates remain bounded. We can adopt the following weaker assumptions to ensure that the iterates do not tend to +oo or
Assumption 5.
(i) h is strongly convex on any bounded convex subset of its domain.
(ii) For all 0, Vh and F are Lipschitz continuous on (V(x) S a).
Under these assumptions, we prove the equivalent of [Bolte et al., 2018, Theorem 4.1].
Theorem Under Assumption 3, and 5, the sequence (a* )kEN generated by (18) with < T <1 converges to a critical point of V.
20

---

We follow the proof given in [Bolte et al., 2018] with some updates to adapt to our weaker Assump- tion In particular, we show Lemma 3 that the sequence (œx)a21 is still "gradient-like" i.e. verifies the assumptions H1, H2 and H3 from [Attouch et al., 2013]. Once these conditions are verified, the result directly follows from [Bolte al., 2018, Theorem 6.2].
Let us first note that by coereivity of y and decrease of the iterates V (I+) (see Proposition 3), the iterates remain in the set
Vk 1 1, as C(ra) {ar dom(h), (x) (2o))-
(88)
Lemma 3. The sequence (11))22 satisfies the following three conditions:
H1 (Sufficient decrease condition)
((xx) (Tx++)) 2 0l - xx 1x
(89)
H2 (Relative error condition) Yk 2 there exists Wk+11 (x*+1) such that
|w++1||| b|| xxl1,
(90)
H3 (Continuity condition) Any subsequence (x"t) converging towards x" verifies
(a + (x=").
(91)
Proof. H1: Sufficient decrease condition. From (85), we get
(<xx) - (I 2 - -1)
Du (n x*t)
(92)
Besides, h is assumed to be strongly convex on any bounded convex subset of its domain. Furthermore, notice that Conv(C(za)) n dom(h) convex subset of dom(h) as intersection of convex sets. Therefore, there is 0 such that
Vx, y E Conv(C(za)) dom(h)) Da(2, v) 2 ox1|> - y/1?.
(93)
With the convention (x, v) = +oo if x & dom(h) or y & int dom
Vx, y E Conv(C(ra)). Da(a,y) 2 ox1|= v|1?.
(94)
As vk 2 1, E C(xa), we get
(xx) (Z Th (- -1) || mull?
(95)
which proves (HI).
H2: Relative error condition. Given (76), the optimality condition for the update of Ø641 is
0 E (R(**+1) VF(**) _ (V/(at+1)- - vh(**)).
1
(96)
For
Wk+1 VF(a*tl)- - VF(a*)+ + [(V((=t)- - vh(2)+11)
(97)
have
Wh+1 E (1((1++1))+ - (F(att1)
(98)
and
| $ [|VF(***l) - VF + -1V(((>)) vh(a++1111
(99)
By assumption, VF and Vh are Lipschitz continuous C(zo) = (V(x)) < (zzo)). As seen before, Ak 2 1, Fy E C(zo). Thus, there is such that
|| |1 < [+)++ VF(**)I I| + vh(ak+1)]] <
(100)
21

---

H3: Continuity condition. Let subsequence (a%) converging towards x". Using the optimality in the update of z we have
R(x*) + (xt x VF 1)) + Du (ot, 21-11)
(101)
S R(x))) + (x - x-1 VF -- ))> + Dn (x" :
(102)
+ R(x²)= S R(xx) + (x' xh (F((*=1)) + DD ( - Dn (2t 22 (103)
From (95) and the fact that (V(x))) converges, we have that || 14 + O. Thus (x also converges towards x²_ In addition, since h is continuously differentiable, (x 1b-1) h(x*) - h(**-1) - (Vh(* ,,==" 1 + O. Passing to the limit in (103), we get
lim R(x*t) S R(x").
(104)
+0+
By lower semicontinuity of R and continuity of F, we get the desired result:
R(x + F(x*t) - R(x") F(x").
(105)
Backtracking The convergence actually requires to control the NoLip constant. In order to avoid small stepsizes, we adapt the backtracking strategy of [Beck, 2017, Chapter 10] the BPG algorithm.
Given y E (0, 1), 7 E [0, 1) and an initial stepsize To > 0, the following backtracking update rule on is applied each iteration A:
while ((xx) - ((T, <
- (D)(C(( T TJT.
(106)
Proposition 4. At each iteration of the algorithm, the backtracking procedure (106) is finite and with backtracking, the convergence results of Proposition 3 and Theorem 3 still hold.
Proof. For a given stepsize т, we showed in equation (85) that
$ (2x) 1 (T+()))) 2 ( - (H(T)
(107)
Taking T < 1ny, we get 1 L > 2 so that
(xx)) - + (T+(T)) > 7 Dr (Tx),,
(108)
T
Hence, when < 1-y, the sufficient decrease condition is satisfied and the backtracking procedure (T 4- 177) must end.
Replacing the former sufficient decrease condition (107) with (108), the rest of the proofs from Proposition 3 and Theorem 3 are identical.
D.2 Proof of Theorem
Proof. We recall the B-RED algorithm
1 arg min(ic(x) + (x /h, (x),(x")) + Dx (109)
(B-RED)
2 A++
E T, (xx))
wER=
It corresponds to the BPG algorithm (14) with F = Fx,, = Xf + 9, and R = ic.
Theorem 1 is direct application of Proposition 4 that is to say of the convergence results of Proposition and Theorem with backtracking. Given Assumptions and 2, we verify that Assumptions 3, 4 and 5 verified:
Assumption 3. R iG verifies dom(R) n int dom(h) C int dom(() + 0. Moreover, R. is semi-algebraic and thus subanalytic as the indicator function of a closed semi-algebraic set. g, and
22

---

being both assumed subanalytic and bounded by [Shiota, 2012. 1.2.1.9], their sum Xf + Je is then subanalytic (up to adding a constant make and 9n non-negative). v Fx,, i is then also subanalytic and thus KI.. V/ is also lower-bounded and coercive f is lower-bounded and coercive and 9y is lower-bounded. Finally, for a E int do)( T,(x) is non-empty as + Xic supercoercive.
Assumption By summing convex functions, using Assumption ( and Assumption 2(ii). L ALy + Ly verifies Lh - (Af + 9y) convex on Conv(dom R) n int dom(h) C' n int dom(h).
Assumption 5. As 9, is assumed to have globally Lipschitz continuous gradient, this follows directly from Assumption l(iv).
D.3 Proof of Theorem
(B-PnP)
k By Vh' XV
(110)
corresponds to the BPG algorithm (14) with F Af and R o.
Theorem is direct application of the convergence results of Proposition and Theorem 3. We now denote V Af + .0y- Given Assumptions 1, we verify that Assumptions 3, and are verified.
Assumption For R p.,, we have Im (B.) C do and as for y int dom(h), Im(B,) C int dom(h) (by equation 13), we get dom (R) n int (h) * 0. d, and f being both assumed subanalytic and bownee by [Shiota, 2012, 1.2.1.91, their sum Aʃ + &er is then subanalytic (up to adding constant to make f' and &e non-negative) and thus KI., Vr is coercive because f is coercive and 0, is lower-bounded. Finally, (1) well-posed is ensured via Proposition
Assumption 4. This is the L-smad property of given by Assumption 1. Assumption 5. This is directly given by Assumption 1(iv).
E
Application to Poisson Inverse Problems
E.1 Burg's entropy Bregman noise model
For h being the Burg's entropy (22), the Bregman noise model (6) writes for x, RT as
p(yla) exp (->Di (x, 2) + p(x))
exp(p(x)) exp
(h(x)) - h(y) (Vh(y), x -1))
cxp(p(x)) exp
log(<;) ) + llg(ys) 3ig 9ii)
(111)
exp(p(x) + +
II
7
exp
-
Zi
3/
y/b
E.2 Inverse Gamma Bregman denoiser
Derivation of (29) We first derive here the condition for convexity n7: 1yy Vh introduced in Section 5.1.
We have
(n)()) = VV
(112)
and
(V
(113)
which gives
The (Vh(y)) =
(v) (VV)))
(114)
For Burg's entropy and y Vh*(x) - with < 0, this writes
(z) = ( (5) 2y Diag Vel, (v)).
(115)
23

---

Using
(y(()) -h(y) - (y( -1,
(116)
Vey(() (4yy Vgg(1) 1
(g
(117)
and
(v) = (2h(0) v2g((0) = y²
v2 gv)),
(118)
we get
2y (x)) = (y
+ = =1²² - -
(119)
= (1- -yyV (1)
For 17y to be convex, necessary condition is that Vy E + and d R"
(v),, 20
(120)
i.e. ('V d) S
Σ (1 - Vy (
(121)
Training details For N,, use the DRUNet architecture from [Zhang et al., 2021] with 2 residual blocks at each scale and softplus activation functions. We condition the network N, on y similarly to what is done for DRUNet. We stack to 3 color channels of the input image additional channel containing an image with constant pixel value equal to 1/7. We use the same training dataset as [Zhang et al., 2021]. Training is performed with ADAM during 1200 epochs. The learning rate is initialized with learning rate 10-1 and is divided by 2 at epochs 300, 600 and 900.
Validation of the convexity hypothesis Given the trained Bregman Score Denoiser B, = y V (y), we propose to verify validity of the hypothesis of convexity from Proposition 1 i.e. the convexity of i/ x + thy o Vh'( (x) = 1 ( (=1) on inl dom(k)" = IR__. Vi is defined in (11). We showed above that for E R and y = 1/72 Vel IR",
(2 Σ (v2 2y Vg, di (W'V (w),, i=1
(122)
In order to verify this assumption around the image manifold, we represent in Figure 4
c(7,))=
min =E(X.),dE(D.)
V2mm
d.
(123)
y/
min (E.
(vē 2y V Vgr d -VVV' (v
(124)
where the minimum is taken over
[x:): the images from an evaluation dataset of clean images, here the 68 images from CSBD68.
{Di}: 100 random images realizations of a uniform random variable D; 4100
and y is obtained from by interpolating between dE ~ pe(u)a:) noisy version of a and B, (ūe)
with a ~ (1 as
V/ aie + (1 (B)).
(125)
In the figures, we represent c( with respect to:
in the z-axis: the denoiser strength, i.e. the parameter of the denoiser. in the y-axis: the noise level & the input image. Figure 4(a), the noise model is Inverse Gamma noise (e(v)=) = Ig(g - 1, Ex) i.e. the noise model used for training the denoiser. In figure 4(b), the noise model is Poisson Pe (y]xx) = P(fr) i.e. the noise model of the degradation on which our plug-and-play algorithms are evaluated, and in Figure 4(c), the noise model Gaussian ( x = (², Id).
We observe that for a large variety of noise levels 7 and even far away from the image manifold i.e. for large noise, verify the convexity assumption.
24

---

180
140
220
F nille loval
200
80
40
30-
960
20
880
denoiser
deneizer
(a) Inverse Gamma noise
(b) Poisson noise Pe (Mc (z) = Ig(Ex).
(0c (z) I9(E - 1, {x).
70
7700
2000
tolee lesel
40
1600
30-
1400
800
50
100
ydaaaaier
(c) Gaussian noise
Pc (0c xx) Ig( Id).
Figure 4: Plot of c( min V
-1 d, d) w.r.t the denoiser parameter and the noise
level in the input image ç. From {x:} the natural images from the CBSD68 testset and 100 random (0, 1]", the image V/ obtained from {a,} via the interpolation y = age + (1 (1-aa) B, where is noisy version of sampled with (a) Inverse Gamma (b) Poisson (c) Gaussian noise distributions. B. is the trained Bregman Score Denoiser.
E.3 The MMSE denoiser verifies the convexity hypothesis from Proposition
Lemma For Burg's 's entropy, and thus for the noise model (23), with the MMSE denoiser fMMSE(N (v) = E|w wy), we have T4 x + w.y Vh is strictly convex int dom(h"), where thyy is defined in (11) from 9x = = log
Proof. Recall that the MMSE denoiser corresponds to the Bregman Score Denoiser
By (g) (ibb o Vh") Vh(y) Vn, Vh(y)
(126)
in the particular case where we have exactly
9y (v)
logpy (v).
(127)
25

---

Note that in the Euclidean case (i.e. for the L Bregman divergence h = 211-11 and = Po = P 2 Id)), it is shown in Gribonval [2011] to be verified by the MMSE denoiser. We now show the same result for Burg's entropy Bregman potential. For x E R
h(z) = Σ log(z;),
(128)
Recall that, in this case, the Bregman noise model writes
p(y|z) exp(p(x) + mY)
II
exp
(-1
a(x)
II (1)) exp
(-9
(129)
yi
Vi
Vi
where axxx) exp(p(>) + n0)) I x > O. We verify that Vor E int domn(h")
= R
(2)dd, > 0.
(130)
We showed in (119) that, Vy E R
VnH n (h(y)) Diag (1- 2y (v)))
(131)
Similar to Gribonval [2011], for simplicity. we first write the proof in the single variable case (n = 1).
We have (gr(v)
-1 ply (v) py (M)
and (ッ) 4 (w)) (v (v)e and thus
14 (v)
2y (3)
3² (y)22)) p² (v)Px())
(x)(x) y²
1+
(132)
7 py (v)
7
(v(
y²
(ND (v)
3² (2 +
(133)
ApX (v)
Moreover
(v) -f (((2)xx)>)) exp (-rp))
2 dz / u((,, u(z,))dd.
(134)
(v) = _7 pv((u) + 2 y
3/2 /. zu(2z
-7 (v) + y
(2)
(135)
y
where we denote (()()) = S2 zu(z, v)dz.
(v) - (7² +7) _ for(()) 2(72 +7) y3
1 +t + 72 1
13)
(136)
y
where we denote (1² ʃʃ =² y)dz. Eventually, we get, after simplification.
(x))) = 2² 72
(2
(2))
(y)PY (v) (1)12)
(137)
(v) y2
2
+ 2
22
u12, v)(('', ddzde
(138)
(v)
(v) [.]
(2
22)2
u(z, y)u(z', vddzdz
(139)
TV 27
2/)2
2 (-r=+²
(v)
2 exp
(d')d 2 0
(140)
As a(x) 0, the previous term is 0 if and only if px(=)px(=!) when 2', which would imply px (x(=) = 0 Vz. This is impossible since px a proper pdf.
The extension for > 1 is straightforward and (z z²)² becomes (a 22,d)²
26

---

E.4 Convergence of B-RED and B-PnP algorithm for Poisson inverse problems
In this section we verify that the Burg's entropy h in (22) and the Poisson data likelihood f defined in (2) verify the assumptions required for convergence of B-RED and B-PnP. We remind the expression
of f and
h(z)
Σ log(z;), 1m1
(141)
m
f(x) Σ
y: log
y/
+ a(Ax); - Vs:
(142)
a(Aa);
for A E Rmxx, Note that as done in Bauschke et al. [2017]. denoting (a,)ısisn the columns of A, we assume that a; & Om and A1 S j S m. Ex= ai > 0 such that Ar E IRT if .. This is verified for A representing the blur with circular boundary conditions with the (normalized) kernels used in Section 5.2.
We first check Assumptions 1.
(ii) f is real analytic and thus subanalytic on RT+
(iii) It is shown in [Bauschke et al., 2017, Lemma 7] that f verifies the NoLip assumption, i.e. h is convex on +. for L, 2 [|a|l y stands for the Poisson degraded observation appearing in the definition of f.
(iv) First, h is strongly convex everywhere on its domain except in +00. For C' bounded subset of R4 as <²h(z) = उ ld, we have Vrr E C, Vd E R", (V²h(z)d, d} > sup.ce x7 ||d| indicating that h is strongly convex Lipschitz continuous everywhere on R except close to 0. For both B-RED and B-PnP V/ (z) - +oo when x and { / (z) S a}} avoids the case æ 0.
bounded subsets of IR Second, h and are
(v) We remind the parameterizations B, Id -Vgy and 9- ||1 N₂ (2)] with U-Net N,, (with softplus activations). The softplus activation function being real analytic, using the fact that the composition and sum of real analytic functions are real analytic [Krantz and Parks, 2002], Nr and then subsequently 9, are real analytic functions.
For the convergence of B-PnP, we also need to prove that Φ, is real analytic (and thus subanalytic) on its domain. For this assumption to be verified, we require in Theorem 2. bb,y Vh not only to be convex, but to be strictly convex. Indeed, using the expression (64) (v) = (w, Vh*) o Vh(y), we have for y E int dom(h)
JB, (v) V (V tb Vh*) Vh(y))
(143)
v²h(y) ( Vh")(Vh(u))
(144)
By bijectivity of Vh between int doo()h) and int dom(h), we have by Vh strictly convex on int dom(h")
(145)
Vz E int dom(h"), ), Vu R", (V Vh") u) >0
(146)
+ Vy E int dom(h), Vu E R", (V o Vh*) M u) > 0.
(147)
Now, using the fact that hyy is positive definite, it follows
Vy E int dom Vu E R"
(V"h(y) o (V")()) u) u) 0
(148)
Vy E int dom(h), Vu E IR", (JB, (8)u, u) > 0
(149)
Js, positive definite on int (h).
(150)
By the inverse function theorem [Krantz and Parks, 2002], as Vy E int dom(h), J.. (v) > 0.
B,1 is then real analytic on (m.)..) We finally obtain, using the expression (12) of Φ, on Im(,.), again by sum and composition, that Φ. is real analytic (and thus subanalytic) on its domain.
Eventually, 9, is non-negative and we now prove that we have Vy E RR. (v) 2 (1) y & int dom(h), Im(B.) C. int dom(f), dy(v) +oo, this is verified. If dom(h).
27

---

as (v. = 0, we have
Φn(1) - (v()) + (v,
2 δy (V) ()(B,()).))
(151)
= Р (v);
where the inequality comes from (13) and the last equality from (12). Therefore by also lower-bounded.
Second, verify Assumption required for the convergence of B-RED.
(i) [0, R is non-empty closed, bounded, convex and semi-algebraic subset of IR
(ii) With the parametrization gy(y) = _11x - N₁ [2) with a neural network Nyy gy can be shown to have Lipschitz gradient (see Hurault et al. [2021, Appendix B] for proof). We can not show a global NoLip property for 9y- However, as Vgr is Lip(g,,) we have Vx (0, Vd I""
(v (x)d,, d) & Liplg>)lld(& S Lip(g-) RR = Lip(g, ( (152)
which proves that, for Ly Lip(g,) Ligh 9, is convex on (0, R]".
B-PnP additional assumptions We finally discuss the additional assumptions required for the convergence of B-PnP in Theorem 2.
Im (B,, C dom(h)= - We train the denoiser B, to restore images in [e, 1J" (with
3)- the denoiser is thus softly enforeed to have its image in this range. In practice, we empirically verify during the iterations that we always get I >
bb Vh convex on ind dom(h"). This is discussed in detail in Appendix E.2.
Im (vh - (V C dom(Vh") We now show that this condition is true if A||v|l < 1. For x
0
xxV
Vh(z) - AV)(z)
XVf(2)
(153)
Thus we need to verify that V1 I term, using (Ax), Ex=1 uz we have sisn. n.
x > O. For f Poisson data-fidelity
1 + Xx, Vf();;
m
Vf(z),
Σ
aaj,i
(154)
Σ
j=1
A=
and
ag
1+ Vo aA
aj
y) Σ aj,0%%
(155)
j=1
We assumed that has positive entries and Σn Oj,3 = r > Therefore, using 0 S 0g,
< 1, we get
1+ Ar, Vf(x) 2 Allx/11
(156)
which positive Al/n// < 1.
The stepsize condition XL < 1. Using the NoLip constant proposed in [Bauschke et al., 2017, Lemma 7] L, = ||y||ı, the condition boils down to Al|ø/l < 1. The condition A|lvlla < 1 is too restrictive in practice, as ||v|/ could get very big, especially for large images. This is due to the fact that the NoLip constant I. 2 ||y/|| can be largely over- estimated. Indeed, in the proof of [Bauschke al., 2017, Lemma 7] well as in the proof of the previous point, the upper bound
(157)
1= Aj.ATK
can be very loose in practice. For B-RED this not a problem, as we use automatic stepsize backtracking. However, for B-PnP, backtracking is not possible as the stepsize is fixed to т = 1. In order still guarantee the convergence, we adopt the following backtracking- strategy to adjust the regularization parameter X:
28

---

Choose an initial value for À > 0. At each iteration of the B-PnP algorithm, check sufficient decrease of the objective function Xf + 0, i.c.
F(Z2)= F( S 6D, (x)+1,
(158)
some iteration, this condition not satisfied before convergence, we alert the user and restart the algorithm with À <- y7 for some 1j (0, 1). We also let the user know that, for optimal performance, it might be necessary to adjust the regularization parameter 7 of the denoiser, in order compensate for this decrease of X.
In our deblurring experiments, for the proposed value of X, over the variety of blur kernels and noise levels experimented, the sufficient decrease property was always verified and this backtracking algorithm was never activated. This illustrates that ||ll| is rough approximation of the NoLip constant.
E.5 Experiments
We give here more details and results on the evaluation of B-PnP and B-RED algorithms for Poisson image deblurring. We present in Figure the four blur kernels used for evaluation. Initialization is done with zα = y The algorithm terminates when the relative difference between consecutive values of the objective function less than 10- the number of iterations exceeds K 500.
(a)
(b)
(c)
(d)
Figure 5: The 4 blur kernels used for deblurring evaluation. (a) and (b) are real-world camera shake kernels from Levin al. [2009]. (c) is 9 uniform kernel. (d) is a 25 25 Gaussian kernel with standard deviation 1.6.
Choice of hyperparameters For B-RED, stepsize backtracking is performed with y = 0.8 and 17 0.5.
When performing plug-and-play image deblurring with our Bregman Score Denoiser trained with Inverse Gamma noise, for the right choice of hyperparameters À and we may observe the following behavior. The algorithm first converges towards meaningful solution. After hundreds of iterations, it can converge towards a different stationary point that does not correspond to visually good reconstruction. We illustrate this behavior Figure 6 where we plot the evolution the PSNR and of the function values fxxx and g, (II) along the algorithm.
7-471
10700
REE
24000
B-4F3
10400
22500
10100
22500
10:00
21003
(a) (SNR((a)
(b) f(xx
(c) (xt)
Figure Evolution the PSNR, [*x) and 9- (xx) when deblurring with B-RED with the initializa-
tion parameters from Table 2 and without first phase fast decrease of both the data-fidelity term and regularization term values, resulting in fast PSNR increase. After approximately 100 iterations, the regularization continue decreasing and the iterates converge towards a different stationary point with low PSNR.
update after 100 iterations. We observe
This phenomenon can be mitigated by using small stepsize, large and small À values at the expense of slowing down significantly the algorithm. circumvent this issue, we propose to first initialize the algorithm with 100 steps with initial and À values and then to change this parameters for
29

---

the actual algorithm. Note that is possible for B-RED to change the stepsize but not for B-PnI which has fixed stepsize r1 For B-PuP, done in Hurault et [2022] the Euclidean setting. we propose to multiply go by a parameter 0 < < 1 such that the Bregman Score denoiser becomes (y) = y - aa--aVVViia (V"h(y) VV (y)): The convergence B-PnP with this denoiser follows identically, The overall hyperparameters X, ", T and ar for B-PnP and B-RED algorithms for initialization and for the actual algorithm are given in Table 2.
a Initialization
20
40
B-RED
7 50
x= 1.5
X = 2.
2.5
Algorithm r = 0.05 7 500
A=
X = 0.5
À = 0.5
B-PnP
Initialization I 50
=
x=
2.5
Algorithm 0.05 7 = 500
== 0.025
A = 0.025
X = 0.025
Table B-RED and B-PnP hyperparameters
Additional experimental results We provide in Table a quantitative comparison between our algorithms B-RED and B-PnP and other methods. (a) PnP-PGD corresponds to the plug-and and-play proximal gradient descent algorithm 2k+1 = D, o (Id fV with D, the DRUNet denoiser (same architecture than B-RED and B- B-PnP) trained denoiser Gaussian noise. (b) PnP-BPG corresponds to with again the DRUNet denoiser D, trained for Gaussian noise. For both (a) and (b) the parameters and
the B-PnP algorithm 0: 8+1 = D o V (Vh-rV/jj
optimized for each noise level o. (c) ALM Unfolded [Sanghvi et al., 2022] uses the Augmented I Lagrangian Method for deriving 3-operator splitting algorithm that is then trained specifically in an unfolded fashion for image deblurring with a variety of blurs and noise levels a. The publicly available model being trained on grayscale images, for restoring our color images, we treat each color channel independently.
Note that contrary to the proposed B-PnP and B-RED algorithms,, the 3 compared methods do not have any convergence guarantees. We observe that our algorithms performs the best when the Poisson noise is not too intense (a = 40 and a 60) but that PSNR performance decreases for intense noise Inverse Gamma
(a 20). We assume that this is due the fact that the denoising prior trained noise is not powerful enough for such a strong noise. A visual example for a = 20 is given Figure 7. As future direction, we plan investigating how to increase the regularization capacity of the deep inverse gamma noise denoiser to better handle intense noise.
20
40
60
PnP.PGD
23.81
24.41
24.45
PnP-BPG
23.85
24.26
24.71
ALM Unfolded [Sanghvi et al., 2022]
23.39
23.91
24.22
B-RED
23.58
24.54
24.90
B-PnP
23.29
24.54
24.80
Table 3: PSNR (dB) of Poisson deblurring methods the CBSD68 dataset. 4 blur kernels represented Figure for each noise levels a.
PSNR averaged over the
(a) Clean
(b) Observed (14.91dB) (c) B B-RED (22.44dB)
(d) B B-PnP (22.40dB)
Figure 7: Deblurring from the indicated motion kernel and Poisson noise with a = 20.
30


Reviews:
[{"review_id": "kF9gGEBgAq", "reviewer": "Reviewer_6Scw", "summary": "This paper develops a Bregman Plug and Play image restoration algorithm for solving under-determined inverse problems in the presence of Poisson measurement noise. The framework trains the image denoising algorithms used within PnP iterations (the \"Bregman Score Denoiser\") to remove noise with an exponential distribution that depends on the distribution of the measurement noise (as opposed to the Gaussian noise used for training previous PnP denoisers). The paper integrates the proposed denoiser into iterative algorithms to form B-PnP and B-RED. Both algorithms effectively deblur images in the presence of Poisson noise and are provably convergent to fixed points.\n", "questions": "Typos\nLine 90: \"see also a review Kamilov et al.\" --> \"see also a review by Kamilov et al.\"\nLine 197: Puting (i) on a separate line would read cleaner\n", "limitations": "Not discussed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "-Generally well-written\n\n-Well motivated\n\n-I believe technically sound\n\n-Evaluated on several different blur kernels, including those based on camera shake. (Blur kernels were known)\n", "weaknesses": "-Based on table 3 in the appendix, the proposed method does not meaningfully improve performance over existing methods (though it does come with convergence guarantees)\n\n-Comparisons with existing methods are not particularly comprehensive\n\n-In practice, PnP algorithms are quite sensitive to how hyperparameters are chosen and setting them correctly can be a challenge. This algorithm introduces another parameter, gamma.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "- We recognize that in terms of numerical performance, our algorithms compare with existing methods. Yet, as you correctly pointed out, our method stands alone in offering guarantees of convergence. We believe that the combination of comparable performance to these methods along with our convergence guarantees is a relevant novelty. In fact, we have introduced the first convergent and effective technique for Poisson image restoration using a nonconvex and deep prior. Furthermore, although our experiments focus on Poisson Inverse Problems, the core contribution of our work lies in its theoretical advancements. Our theoretical study contains novel convergence findings for the Bregman Proximal Gradient algorithm (see Appendix D.1) and a new characterization of Bregman proximity operators (see Appendix C). These comprehensive theoretical outcomes extend beyond our immediate plug-and-play objectives, holding potential for diverse applications.\n\n- Thanks for your feedback on the comparison section. In the updated version of the publication, we propose to bring more details when decribing the compared methods, as follows: \n\na) PnP-PGD corresponds to the standard plug-and-play proximal gradient descent algorithm in the Euclidean geometry \n$$ x_{k+1} = D_\\sigma \\circ (Id - \\tau \\nabla f)(x_k) $$ \nwhere $f$ is the Poisson data-fidelity term (2). The plugged denoiser $D_\\sigma$ is the DRUNet  network (i.e. the same architecture than the one used for parametrizing our denoiser) but now trained to denoise **Gaussian** noise of std $\\sigma$. We train $D_\\sigma$  simultaneously for all noise levels $\\sigma \\in [0,50]$ with the $L^2$ loss. As $f$ does not have Lipschitz gradient, this PnP-PGD algorithm does not have convergence guarantees.\n\nb) PnP-BPG realizes the same iterations than our B-PnP algorithm (20) \n$$x^{k+1} = D_\\sigma \\circ \\nabla h^*(\\nabla h - \\tau \\nabla f)(x_k)$$ \nbut our Bregman Score Denoiser in (20) is replaced by the more classical Gaussian denoiser $D_\\sigma$ introduced above. This scheme does not have guarantees of convergence as $D_\\sigma$ is not a Bregman proximity operator anymore. \n\nFor both PnP-PGD and  PnP-BPG the parameters $\\sigma$ and $\\tau$ are optimized for each noise level by grid-search.\n\nc) ALM Unfolded [3] uses the Augmented Lagrangian Method for decoupling linear operator $A$ from the data-fidelity term. They derive a 3-step algorithm that is unfolded and trained for specific degradations. In particular, it is trained for image deblurring with a variety of blurs and noise levels. The publicly available model being trained on grayscale images, for restoring our color images, we treat each color channel independently.\n\n- Concerning the hyperparameters involved in our algorithms, the parameter $\\gamma$ is the noise level of the denoiser. $D_\\gamma$ is trained to denoise images degraded with noise parameter $\\gamma$. It is the equivalent of the parameter $\\sigma$ commonly used by PnP method with Gaussian denoisers. Thus, compared to standard PnP algorithms, there is no additional hyperparameter involved. Actually, we even have one hyperparameter less as the stepsize of the algorithm has not to be chosen: for B-RED it is automatically tuned with a backtracking line-search strategy, and for B-PnP it has to be set to $1$. We agree on the fact that hyperparameters selection is a common difficulty in variational image restoration. Given an input image, automatic parameter tuning strategies have been developed and could be added to our framework. For example, [1] proposes a Bayesian method for setting the regularization parameter and [2] employs deep reinforcement learning to train a policy network yielding well-suited parameters for PnP. \n\n- In the global rebuttal to all reviewers, we propose a limitation paragraph which will be added at the end of the main paper.\n\n[1] Vidal et al. \"Maximum likelihood estimation of regularization parameters in high-dimensional inverse problems: An empirical bayesian approach part i: Methodology and experiments.\", 2020\n\n[2] Wei et al. \"Tuning-free plug-and-play proximal algorithm for inverse imaging problems\". In ICML 2020.\n\n[3] Sanghvi et al.  Photon limited non-blind deblurring using algorithm unrolling. IEEE Transactions on Computational Imaging, 2022."}, {"review_id": "GQ2PQZUoVn", "reviewer": "Reviewer_5jWQ", "summary": "This paper studies an extension of the Plug-and-Play (PnP) framework for solving inverse imaging problems by considering descent schemes in metrics different from L2: Motivated by the fact that some data fidelity terms such as the Kullback-Leibler divergence allow for an efficient minimization with the Bregman Proximal Gradient (BPG) method (an extension of proximal gradient descent to arbitrary Bregman distances instead of squared L2), the authors develop a parametrization of a learnable denoiser which can be interpreted as a proximal operator (or descent step) of a cost function w.r.t. the underlying Bregman distance that combines well with a particular data fidelity term. Under some additional assumptions, this allows proving the convergence of the resulting Bregman PnP framework. Numerical experiments illustrate that the resulting scheme can successfully solve deblurring problems with Poisson noise.\n", "questions": "Considering a difficult (2-stage) optimization with different parameters to avoid bad minimizers as shown in Fig. 5 of the supplement, a lack of strict convergence guarantee as the convexity condition cannot be enforced, and negligible difference to a plain PnP or RED approach in terms of the PSNR, what is the advantage of the proposed method?", "limitations": "I don't think there is a potential negative societal impact of this work. In terms of general limitations, I think the authors should be more open about the fact, that the use of the Bregman distance framework did not result in improved results over prior algorithmic schemes. ", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper is technically sound and presents the technical construction of the Bregman PnP approach very well. It is an elegant solution that closes a gap missing in the PnP (and RED) framework with learnable priors. I found it very convincing in terms of its theory and it even derives some general (smaller) missing pieces for convergence beyond the seminal works by Bolte, Bauschke, Teboulle and co-authors. \n\nDespite this strengths section being short in comparison to the weaknesses, I think that the strong theoretical contribution along with an illustration of the implementation outweighs the weaknesses in terms of (benchmark) results, such that I am leaning towards accepting the paper.", "weaknesses": "The numerical experiments/results are not very convincing from a practical point of view. \n- First, the denoising performance of the dedicated network is not better than that of a plain network (might not be very important).\n- Second, the main paper does not compare to the plain PnP approaches with L2 fidelity. The supplementary material reports tiny differences only (with the proposed approach being slightly worse for large noise and slightly better for small noise)\n- According to line 275, the desired condition of $\\phi_\\gamma \\circ \\nabla h^*$ being convex is not enforced but seems to hold empirically when training the network. Thus, any convergence guarantee is lost.\n- Figure 5 in the supplementary material shows a concerning risk of ending up with a bad result. The mitigation strategy of first running 100 iterations with a first set of parameters and then switching to a different set of parameters represents a significant amount of fine-tuning (possibly exceeding the number of hyperparameters and amount of fine-tuning used for the standard L2 PnP approach), such that even the small improvements in table 3 of the supplementary material need to be viewed with care. \n\nMinor aspects:\n- In 225 the author decide to do backtracking line search to avoid estimating the NoLip constant, but would backtracking on the (differentiable but not L-smooth) data fidelity term considered here not work in the L2 case?\n- The authors mention that there is no result on the convergence of the PDHG method for nonconvex regularizers (line 81). As I was curious, I briefly searched online and found \"Precompact convergence of the nonconvex Primal–Dual Hybrid Gradient algorithm\" by Sun et al., Journal of Computational and Applied Mathematics, 2018. Precompactness of the primal variable (the image) would be easy to ensure if one restricts every value to [0,1]. Is their result applicable? (Honestly, I have not read the paper yet). \n- The condition $\\lambda L_f <1 $ in Theorem 2 seems to limit the amount of data fidelity one can use in order to still have a convergent algorithm - is this a limitation?\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:** \n- We do no expect the B-DRUNet denoiser to outperform the DRUNET denoiser. Indeed, both are based on the same architecture but the former is additionally constrained to take the specific form (25). We are thus satisfied with the fact that B-DRUNET almost reach the performance of DRUNET (difference of ~0.05 dB).\n- We recognize that in terms of numerical performance, our algorithms compare with existing methods. Yet, our method stands alone in offering guarantees of convergence. We believe that the combination of comparable performance to these methods along with our convergence guarantees holds great promise. In fact, we have introduced the first convergent and effective technique for Poisson image restoration using a nonconvex and deep prior. \n- Indeed, the convexity condition (26)-(27) is not hard-coded in the network architecture. However, it is naturally verified in practice by the trained denoiser. Please refer to the global rebuttal to all reviewers, where we discuss and plot, in the attached pdf, empirical verifications of the convexity hypothesis. We note a clear confirmation of the condition (26)-(27), even when evaluated far from the image manifold. We also prove in the rebuttal to all reviewers that the assumption is **verified by the true MMSE denoiser**.  By minimizing the $L^2$ cost, the denoiser is actually trained to approximate this MMSE denoiser. This explains why our denoiser naturally satisfies the convexity condition without necessitating supplementary constraint. In comparaison, as explained in the introduction, most of the studies in the literature on the convergence of PnP methods assume *nonexpansive denoisers*. However, the **true MMSE denoiser is not nonexpansive**. Consequently, it is observed in different works [1,2,3] that a network trained to denoise without additional constraints is *not nonexpansive*. Furthermore, ensuring nonexpansivity of a denoiser while training is often achieved through the use of soft penalties or approximations (which also lack explicit guarantees) and significantly degrades denoising performance. \n- On the numerical side, we agree that the adopted two-step process is a limitation of our approach. The learned regularizer being nonconvex, the algorithms can be sensitive to initialization and for strong degradations the algorithm might not converge towards to right critical point if not initialized properly. \n- For Poisson noise, we can show that the data-fidelity term does not verify any NoLip condition (4) w.r.t the L2 metric. However, with Proposition 4 (Appendix D), the backtracking strategy is provably guaranteed to find a stepsize such that the objective function decreases, **provided** that a NoLip is satisfied for some $L>0$. Indeed, Proposition 4 is based on the sufficient decrease property (81) derived from the NoLip condition. Thus, in the L2 case, adding backtracking does not lead to a converging scheme.\n- We are indeed aware of the result from Sun et al. Even though very interesting, the precompactness is however not easy to enforce for the PDHG algorithm. Indeed, hard constraining the values of one variable in [0,1] amounts to add a proximal step (a projection) in the algorithm. In our B-RED algorithm, when adding the hard constraint $i_C$ (18) to the Bregman Gradient Descent algorithm (17), we still fit the general BPG algorithm (14). However, adding a proximal step in the PDHG algorithm does not correspond to some known algorithm. \nNote that one could prove that it is possible to ensure precompactness and thus convergence of nonconvex PDHG if we add the assumption that $f$ has Lipschitz gradient. Such a convergence result would indeed makes sense when compared to the recent convergence results of ADMM/DRS in the nonconvex setting that require one function to have Lipschitz gradient [4]. However, here for Poisson deblurring, $\\nabla f$ is not Lipschitz and this is not applicable. \n- The constraint on $\\lambda$  is indeed a limitation. Thanks for pointing this out. It is due to the fact that the denoiser writes as a proximal step *with stepsize $1$*. We are thus forced to keep a fixed stepsize $\\tau=1$. The BPG condition for convergence $\\tau \\lambda L_f < 1$ then becomes $\\lambda L_f < 1$ i.e. a constraint on $\\lambda$. As a future work, we plan to explore solutions for relaxing this constraint, for instance by using a different algorithm than BPG with a loosen stepsize constraint. \n\n**Questions:** \n\nCompared to other Poisson image restoration methods, the main advantage of our algorithms is their convergence. First and foremost, the convexity condition is required for the convergence of the B-PnP algorithm but not for the convergence of B-RED. B-RED has then strict convergence guarantees. In addition, for the convergence of B-PnP, we refer to our comment above on the practical and semi-theoretical verification of the convexity condition.\n\nFurthermore, although our experiments focus on Poisson Inverse Problems, the core contribution of our work lies in its theoretical advancements. We present the first Bregman extension within the plug-and-play and Regularization-by-Denoising frameworks, supported by a strong theoretical foundation. This encompasses novel convergence findings for the Bregman Proximal Gradient algorithm (see Appendix D.1) and a new characterization of Bregman proximity operators (see Appendix C). These comprehensive theoretical outcomes extend beyond our immediate plug-and-play objectives, holding potential for diverse applications.\n\n[1] Hertrich et al. \"Convolutional proximal neural networks and PnP algorithms\". In Linear Algebra and its App., 2021.\n\n[2] Pesquet et al. \"Learning maximally monotone operators for image recovery\" SIAM J. on Im. Sc., 2022.\n\n[3] Hurault et al. \"Proximal denoiser for convergent PnP optimization with nonconvex regularization\". ICML 2022.\n\n[4] Themelis & Patrinos. \"DRS and ADMM for nonconvex optimization: Tight convergence results.\" SIAM J. on Optim, 2020."}, {"review_id": "bCTuJ4guyr", "reviewer": "Reviewer_oqFb", "summary": "Inspired by the No-Lips literature on the optimization of convex objectives which are not globally L-smooth, this paper adds to the PnP literature the “Bregman score denoiser” which extends the BPG algorithm by a Bregman-based prox-map along with convergence conditions despite NN-parametrized non-convex potentials are involved.", "questions": "Before eq. (14), the convexity assumption regarding F and R, adopted in [Bauschke et al 2017] to which you refer, is missing.\n\nline 92: why is nonexpansiveness considered as an “unrealistic requirement” to ensure convergence pf PnP, in view of your own assumptions which are not easy to check either?\n\nHow is the set C in (18) determined?", "limitations": "This point has not been explicitly addressed in the paper, apparently. Please comment.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "Originality:\nThis is a plausible extension of recent related work on PnP networks based on state-of-the-art theory from the field of numerical optimization.\n\nQuality:\nIt is apparant that the authors do not both fields very well.\n\nSignificance: \nThe paper adds a new concept to the PnP literature.", "weaknesses": "Clarity: \nThe presentation intersperses references, top-level arguments and technical details in a confusing manner. I had to read few times forth and back in order to get an idea what this paper is about.\n\nAuthors criticize “unrealistic assumptions” (e.g. line 92) in related work but have to admit later on that their own assumptions are hard to check as well, even for a simple scenario (lines 275-277). The need for backtracking line search is not convenient either.\n\nSignificance: \nRegarding the theory, I did not get if some `generic’ properties of the Bregman-prox-maps exist, that play the very same role like, say, the firmly-nonexpansiveness of Euclidean prox-maps, to achieve convergence in the considered generalized scenarios. The concrete scenario (Poisson noise) is classical. A comparison to a related approach shows no improvement (Fig.1, (c), (d)). In particular, artificially corrupting images looks like old-style image denoising papers.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:**\n\n- We first bring more details on our assumption of convexity of $\\psi\\_\\gamma \\circ \\nabla h^\\star(x)$, compared to the common assumption of nonexpansivity of the denoiser. \n\nIt is observed in different works [1,2,3] that a state-of-the-art network trained to denoise without additional constraints is **not nonexpansive**, and constraining a network during training to be nonexpansive can severely degrade its denoising performance (see for example [1, Table 1], [2, Figure 1], [3, Table 5&6]). A semi-theoretical reason for this fact is that, by minimizing the $L^2$ cost, the denoiser is trained to approximate the true MMSE denoiser, which is *not nonexpansive*. This is why we call nonexpansivity an \"unrealistic assumption”. Conversely, as proved in the global rebuttal to all reviewers, our assumption (convexity of $\\psi\\_\\gamma \\circ \\nabla h^\\star(x)$) **is verified by the true MMSE denoiser**, and is thus more realistic. This explains why, without any additional constraints, our trained denoiser verifies this assumption.\n\nWe would like to put forward the distinction between *verifying (or checking)* the validity of the assumption, after training, and *enforcing* the assumption, while training. Note that our convexity condition is not *\"hard to check\"* but *hard to enforce explicitly*. Indeed, the convexity condition can be easily verified empirically, after training, with (26)-(27). Please refer to the global rebuttal to all reviewers, where we discuss and plot, in the attached pdf, empirical verifications of this convexity. We observe that the condition (26)-(27) is clearly verified even far way from the image manifold. \n\nIn comparison, coming back to the nonexpansivity  assumption, finding the Lipschitz constant of a neural network being NP-hard, *verifying* the nonexpansivity of a neural network is typically done with a similar empirical verification. Moreover, previous methods *enforcing* nonexpansivity of a denoiser while training do not have better theoretical guarantees. Indeed, existing methods enforce nonexpansivity by either:\n- normalizing each layer by its spectral norm, with the spectral norm *approximated* via the power method [4]. As this is only a local approximation of the Lipschitz constant, there is no strict guarantee to get nonexpansivity.\n- regularizing the training loss [5,3] and thus without strict guarantee either.\n\nLet us also point out that the convexity condition is required for the convergence of the B-PnP algorithm but not for the convergence of B-RED. B-RED has then exact convergence guarantees.\n\n- Concerning the backtracking line search strategy, it allows to keep a fast algorithm along with convergence guarantees as it automatically sets the stepsize to its maximal value for convergence. Backtracking is commonly used when the Lipschitz (or NoLip) constant of the gradient of the smooth potential is unknown, which is typically the case when regularizing inverse problems with deep explicit priors, see for example [6]. \n\n- Regarding the \"generic\" properties of the Bregman-prox-maps, Euclidean proximity operators of *nonconvex potentials* are not firmly-nonexpansive (this is actually why our denoiser is not nonexpansive) . The convergence of proximal algorithms like Proximal Gradient Descent in the nonconvex setting is only due to the *sufficient decrease property* obtained by combining the first order optimality of the prox and the descent lemma on the potential with Lipschitz gradient. This is the same idea for the proof for BPG algorithm with a Bregman prox, where the descent lemma holds thanks to the NoLip (4) property, which replaces the Lipschitz gradient assumption. \n\n- Concerning denoising, Figure 1 is not evaluating Poisson denoising but Inverse Gamma noise denoising. In this experiment, we do not expect the B-DRUNet denoiser to outperform the DRUNET denoiser. Indeed, both are based on the same architecture  but the former is additonally constrained to take the specific form (25). We are thus satisfied with the fact that this specific B-DRUNET almost reach the performance of DRUNET (difference of ~0.05 dB).\n\n- Eventually, regarding the artificial corruption of images, the purpose of training a denoiser using artificial Inverse Gamma noise (or artificial Gaussian noise for standard Plug-and-Play methods) is to establish a robust deep prior that effectively regularizes inverse problems. The selection of the (artificial) Bregman noise model holds significant importance in ensuring that our Bregman Score denoiser offers a meaningful prior. Once this prior is established through training, it becomes applicable in an unsupervised manner for any inverse problem, even for restoring images corrupted with real noise.\n\n**Questions:**\n\n- The NoLip condition will be added before eq (14), thanks for noticing this omission. \n- Concerning the \"nonexpansivity considered as an unrealistic requirement”, please refer to our comment above. \n- In practice, we choose $C = [0,1]^n$.\n\n**Limitations:** \n\nThanks for pointing this out. In the global rebuttal to all reviewers, we provide a limitation paragraph which will be added to the conclusion. \n\n[1] Hertrich et al. \"Convolutional proximal neural networks and PnP algorithms\". In Linear Algebra and its App., 2021.\n\n[2] Bohra et al. \"Learning lipschitz-controlled activation functions in neural networks for PnP im. rec. methods\". NeurIPS Workshop on DL and Inv Prob, 2021.\n\n[3] Hurault et al. \"Proximal denoiser for convergent PnP optimization with nonconvex regul.\". ICML 2022.\n\n[4] Ryu et al. \"PnP methods provably converge with properly trained denoisers\". ICML 2019.\n\n[5] Pesquet et al. \"Learning maximally monotone operators for im. recovery.\" SIAM Journal on Im Sciences, 2021\n\n[6] Romano et al. \"The little engine that could: Regularization by denoising (red)\". SIAM Journal on Im Sciences, 2017."}]

Author rebuttal:
We would like to thank all reviewers for their careful reading of our submission, and their helpful comments and suggestions.  In each individual rebuttal, we tried to answer to all the objections raised by the reviewers. We give here more general remarks about aspects that multiple reviewers have highlighted.

**Regarding the hypothesis of convexity of $\psi\_\gamma \circ \nabla h^\star(x)$ in Proposition 1**

- We first proceed to a detailed experimental validation of this assumption (find the plots in the attached pdf):

After training the Bregman Score Denoiser $\mathcal{B}\_{\gamma}(y) = y - y^2 \nabla g\_\gamma(y)$ to denoise Inverse Gamma noise at different noise levels $\gamma$, we verify the convexity of $\eta_\gamma : x \to \psi_\gamma \circ \nabla h^\star(x) = \psi_\gamma\left(\frac{-1}{x}\right)$ on $int\ dom(h^\star) = \mathbb{R}^n\_{--}$ where $\psi_\gamma$ is defined in (11). 

In appendix E.2, we showed that for $z \in \mathbb{R}^n\_{--}$ and $y = - 1/z$,  $\forall d \in \mathbb{R}^n$,

$\bigl \langle \nabla^2 \eta_{\gamma}(z) d,d \bigr \rangle = \sum_{i=1}^n \left( y^2(1 - 2y \nabla g\_\gamma(y)) \right)_i d_i ^2 - \bigl \langle Diag(y^4)\nabla^2 g\_\gamma(y)d, d \bigr \rangle = \mathcal{C}\_\gamma(y, d)$ 

To confirm the convexity assumption around the image manifold, we need to verify that the above quantity is positive for any image $y$. We represent, **in the attached pdf**

$c(\gamma, \xi) = \min_{y_\xi, d} \biggl \langle \nabla^2 \eta_{\gamma}\left(\frac{-1}{y_\xi}\right) d,d \biggr \rangle =  \min_{y_\xi, d} \mathcal{C}\_\gamma(y_\xi, d)$

The minimum is taken over $\\{y_{\xi}\\}$ and $\\{d\\}$ where
- Given $\\{ x_i \\}$ the CSBD68 testset of clean and natural images, $y_\xi$ is obtained from $x \in \\{x_i\\}$ by interpolating between $\hat y_\xi \sim p_\xi(y|x)$ a noisy version of $x$ and $\mathcal{B}\_{\gamma}(\hat y_\xi)$ via $y\_\xi = \alpha \hat y\_\xi + (1-\alpha) \mathcal{B}\_{\gamma}(\hat y\_\xi)$, where $\alpha \sim \mathcal{U}[0,1]$ . We use different noise models $p_\xi$ and noise levels $\xi$. This enables to explore the space around the image manifold. 
- For each test image, $100$ random vectors $d$ are sampled from $\mathcal{U}[0,1]^n$.

In the figures, we represent $c(\gamma, \xi)$ w.r.t: 
- $x$-axis: the $\gamma$ parameter of the denoiser.
- $y$-axis: the noise level $\xi$ in the input image. Figure (a), the noise model is Inverse Gamma noise $p_\xi(y|x) = \mathcal{IG}(\xi-1, \xi x)$ i.e. the noise model used for training the denoiser. Figure (b), the noise model is Poisson $p_\xi(y|x) = \mathcal{P}(\xi x)$, i.e. the noise model on which our plug-and-play algorithms are evaluated, and Figure (c), the noise model is Gaussian $p_\xi(y|x) = \mathcal{N}(x,\xi^2 Id)$.

We observe that for all  $\gamma$ seen during training, and even far away from the image manifold (i.e. for large noise level) we have $c(\gamma, \xi) > 0$  by a large margin, i.e. the convexity condition of Proposition 1 is clearly verified. These plots will be added to the appendix of the paper.  

- We also provide a theoretical argument supporting the assumption:

Our denoiser is trained by minimizing the $L^2$ cost for which the optimal Bayes estimator is the MMSE. The MMSE is then a theoretical denoiser that our denoiser tries to approximate. **We can show that the proposed assumption ($\psi\_\gamma \circ \nabla h^\star(x)$ convex) is verified by the MMSE denoiser.** This explains why our denoiser, after training, naturally satisfies this condition without necessitating supplementary constraints. This result is proven in [1, Lemma A.1] in the Euclidean case, we here extend the proof for Burg's entropy Bregman potential.

We summarize the proof here without including the details of the calculus. We use the same notations as above. With $h(x) = \sum_{i=1}^n \log(x_i)$ Burg's entropy, it consists in showing that, for 

$g_\gamma(y) = - \frac{1}{\gamma} \log p_Y(y)$ (i.e. when the denoiser is the MMSE) 

we have 

$\forall z \in \mathbb{R}^n_{--}, \langle \eta_\gamma(z) d,d \rangle \geq 0$. 

We showed in Appendix E.2 that  $\forall z \in \mathbb{R}^n_{--}$ and  $y = \nabla h^*(x) = -\frac{1}{z}$,

$\nabla^2 \eta\_\gamma(z) = Diag \left(y^2(1 - 2y \nabla g\_\gamma(y)) \right) - Diag(y^4) \nabla^2 g\_\gamma(y)$

We consider the single variable case $n=1$. After differentiating $g_\gamma$ twice, we get 

$\eta_\gamma''(z) = \frac{y^2}{\gamma p_Y^2(y)} \left( \gamma p_Y^2(y) + 2y p'_Y(y)p_Y(y) - y^2[p'_Y(y)]^2  + y^2 p''_Y(y)p_y(y) \right)$

Recall that for Burg's entropy, the Bregman noise model writes $p(y|x) = \alpha(x) \prod_{i=1}^n (y\_i)^{-\gamma} \exp\left(- \gamma \frac{x\_i}{y\_i}\right)$. 
Differentiating $p_Y(y)= \int\_z p(y|z)p\_X(x)dx$ twice, after simplification we get 

$\eta_\gamma''(z) = \frac{\gamma}{p_Y^2(y)} \int_x \int_{x'} \frac{(x-x')^2}{2} p(y|x)p(y|x')p(x)p(x')dxdx' \geq 0$

The extension to $n>1$ is straightforward, where $(x-x')^2$ becomes $\langle x-x', d \rangle^2$.

The detailed proof will be added to the appendix of the paper.  

**Regarding the limitations, this section will be added in the main paper**:

The central significance of our work stems from its theoretical study but we recognize certain limits within our experimental results. First, when applied to deblurring with Poisson noise, our proposed algorithms do not outperform existing methods in terms of PSNR. Second, while we prove that B-RED is convergent without restriction, the convergence of B-PnP depends on a specific convexity condition. Despite being confirmed with experiments and having robust theoretical foundations, this assumption could potentially be not verified when applied to non-natural images that significantly differ from those in the training dataset. Finally, due to the nonconvex nature of our proposed prior, the practical performance of the algorithms can be sensitive to their initialization."

[1] Gribonval, "Should penalized ...." IEEE Trans on Signal Proc., 2011.

Meta-review: The paper develops plug-n-play methods for inverse problems whose losses are not gradient-Lipschitz. The most significant example of this is inverse problems with Poisson noise, which are prevalent in medical and scientific imaging. The paper develops a framework for this problem, in which one trains a ``Bregman score denoiser’’, which at reconstruction time can be deployed as a component in a Bregman proximal method, with provable convergence. 

Reviewers found the paper’s proposals to be well motivated, addressing a need in the PnP/RED frameworks. Reviewers found the paper to be theoretically solid, well written, and sufficiently novel to merit acceptance. Minor concerns involved the experimental results (where the proposed methods match, but do not exceed, the state of the art), and the existence of prior art on learnable inverse problem solvers with poisson noise. In this regard, the main contributions of the paper are (i) to develop a framework with provable convergence, and (ii) proving convergence. After considering the authors responses, the reviewers converged to a recommendation to accept the paper.
Decision: Accept (poster)