Year: 2022

Abstract:
Plug-and-Play (PnP) methods constitute a class of iterative algorithms for imaging problems where regularization is performed by an off-the-shelf denoiser. Although PnP methods can lead to tremendous visual performance for various image problems, the few existing convergence guarantees are based on unrealistic (or suboptimal) hypotheses on the denoiser, or limited to strongly convex data-fidelity terms. We propose a new type of PnP method, based on half-quadratic splitting, for which the denoiser is realized as a gradient descent step on a functional parameterized by a deep neural network. Exploiting convergence results for proximal gradient descent algorithms in the nonconvex setting, we show that the proposed PnP algorithm is a convergent iterative scheme that targets stationary points of an explicit global functional. Besides, experiments show that it is possible to learn such a deep denoiser while not compromising the performance in comparison to other state-of-the-art deep denoisers used in PnP schemes. We apply our proximal gradient algorithm to various ill-posed inverse problems, e.g. deblurring, superresolution and inpainting. For all these applications, numerical results empirically confirm the convergence results. Experiments also show that this new algorithm reaches state-of-the-art performance, both quantitatively and qualitatively.

heading: 1 INTRODUCTION
text: Image restoration (IR) problems can be formulated as inverse problems of the form
x∗ ∈ arg min x f(x) + λg(x) (1)
where f is a term measuring the fidelity to a degraded observation y, and g is a regularization term weighted by a parameter λ ≥ 0. Generally, the degradation of a clean image x̂ can be modeled by a linear operation y = Ax̂+ ξ, where A is a degradation matrix and ξ a white Gaussian noise. In this context, the maximum a posteriori (MAP) derivation relates the data-fidelity term to the likelihood f(x) = − log p(y|x) = 12σ2 ||Ax− y|| 2, while the regularization term is related to the chosen prior.
Regularization is crucial since it tackles the ill-posedness of the IR task by bringing a priori knowledge on the solution. A lot of research has been dedicated to designing accurate priors g. Among the most classical priors, one can single out total variation (Rudin et al., 1992), wavelet sparsity (Mallat, 2009) or patch-based Gaussian mixtures (Zoran & Weiss, 2011). Designing a relevant prior g is a difficult task and recent approaches rather apply deep learning techniques to directly learn a prior from a database of clean images (Lunz et al., 2018; Prost et al., 2021; González et al., 2021).
Generally, the problem (1) does not have a closed-form solution, and an optimization algorithm is required. First-order proximal splitting algorithms (Combettes & Pesquet, 2011) operate individually on f and g via the proximity operator
Proxf (x) = arg min z
1 2 ||x− z||2 + f(z). (2)
Among them, half-quadratic splitting (HQS) (Geman & Yang, 1995) alternately applies the proximal operators of f and g. Proximal methods are particularly useful when either f or g is nonsmooth.
Plug-and-Play (PnP) methods (Venkatakrishnan et al., 2013) build on proximal splitting algorithms by replacing the proximity operator of g with a generic denoiser, e.g. a pretrained deep network.
∗Corresponding author: samuel.hurault@math.u-bordeaux.fr
These methods achieve state-of-the-art results (Buzzard et al., 2018; Ahmad et al., 2020; Yuan et al., 2020; Zhang et al., 2021) in various IR problems. However, since a generic denoiser cannot generally be expressed as a proximal mapping (Moreau, 1965), convergence results, which stem from the properties of the proximal operator, are difficult to obtain. Moreover, the regularizer g is only made implicit via the denoising operation. Therefore, PnP algorithms do not seek the minimization of an explicit objective functional which strongly limits their interpretation and numerical control.
In order to keep tractability of a minimization problem, Romano et al. (2017) proposed, with regularization by denoising (RED), an explicit prior g that exploits a given generic denoiser D in the form g(x) = 12 〈x, x −D(x)〉. With strong assumptions on the denoiser (in particular a symmetric Jacobian assumption), they show that it verifies
∇xg(x) = x−D(x). (3) Such a denoiser is then plugged in gradient-based minimization schemes. Despite having shown very good results on various image restoration tasks, as later pointed out by Reehorst & Schniter (2018) or Saremi (2019), existing deep denoisers lack Jacobian symmetry. Hence, RED does not minimize an explicit functional and is not guaranteed to converge.
Contributions. In this work, we develop a PnP scheme with novel theoretical convergence guarantees and state-of-the-art IR performance. Departing from the PnP-HQS framework, we plug a denoiser that inherently satisfies equation (3) without sacrificing the denoising performance. The resulting fixed-point algorithm is guaranteed to converge to a stationary point of an explicit functional. This convergence guarantee does not require strong convexity of the data-fidelity term, thus encompassing ill-posed IR tasks like deblurring, super-resolution or inpainting.

heading: 2 RELATED WORKS
text: PnP methods have been successfully applied in the literature with various splitting schemes: HQS (Zhang et al., 2017b; 2021), ADMM (Romano et al., 2017; Ryu et al., 2019), Proximal Gradient Descent (PGD) (Terris et al., 2020). First used with classical non deep denoisers such as BM3D (Chan et al., 2016) and pseudo-linear denoisers (Nair et al., 2021; Gavaskar et al., 2021), more recent PnP approaches (Meinhardt et al., 2017; Ryu et al., 2019) rely on efficient off-theshelf deep denoisers such as DnCNN (Zhang et al., 2017a). State-of-the-art IR results are currently obtained with denoisers that are specifically designed to be integrated in PnP schemes, like IRCNN (Zhang et al., 2017b) or DRUNET (Zhang et al., 2021). Though providing excellent restorations, such schemes are not guaranteed to converge for all kinds of denoisers or IR tasks.
Designing convergence proofs for PnP algorithms is an active research topic. Sreehari et al. (2016) used the proximal theorem of Moreau (Moreau, 1965) to give sufficient conditions for the denoiser to be an explicit proximal map, which are applied to a pseudo-linear denoiser. The convergence with pseudo-linear denoisers have been extensively studied (Gavaskar & Chaudhury, 2020; Nair et al., 2021; Chan, 2019). However, state-of-the-art PnP results are obtained with deep denoisers. Various assumptions have been made to ensure the convergence of the related PnP schemes. With a “bounded denoiser” assumption, Chan et al. (2016); Gavaskar & Chaudhury (2019) showed convergence of PnP-ADMM with stepsizes decreasing to 0. RED (Romano et al., 2017) and RED-PRO (Cohen et al., 2021) respectively consider the classes of denoisers with symmetric Jacobian or demicontractive mappings, but these conditions are either too restrictive or hard to verify in practice. In Appendix A.3, more details are given on RED-based methods. Many works focus on Lipschitz properties of PnP operators. Depending on the splitting algorithm in use, convergence can be obtained by assuming the denoiser averaged (Sun et al., 2019b), firmly nonexpansive (Sun et al., 2021; Terris et al., 2020) or simply nonexpansive (Reehorst & Schniter, 2018; Liu et al., 2021). These settings are unrealistic as deep denoisers do not generally satisfy such properties. Ryu et al. (2019); Terris et al. (2020) propose different ways to train deep denoisers with constrained Lipschitz constants, in order to fit the technical properties required for convergence. But imposing hard Lipschitz constraints on the network alters its denoising performance (Bohra et al., 2021; Hertrich et al., 2021). Yet, Ryu et al. (2019) manages to get a convergent PnP scheme without assuming the nonexpansiveness ofD. This comes at the cost of imposing strong convexity on the data-fidelity term f , which excludes many IR tasks like deblurring, super-resolution or inpainting. Hence, given the ill-posedness of IR problems, looking for a unique solution via contractive operators is a restrictive assumption. In this work, we do not impose contractiveness, but still obtain convergence results with realistic hypotheses.
One can relate the ideal deep denoiser to the “true” natural image prior p via Tweedie’s Identity. In (Efron, 2011), it is indeed shown that the Minimum Mean Square Error (MMSE) denoiser D∗σ (at noise level σ) verifies Dσ(x) = x+ σ2∇x log pσ(x) where pσ is the convolution of p with the density ofN (0, σ2 Id). In a recent line of research (Bigdeli et al., 2017; Xu et al., 2020; Laumont et al., 2021; Kadkhodaie & Simoncelli, 2020), this relation is used to plug a denoiser in gradient-based dynamics. In practice, the MMSE denoiser cannot be computed explicitly and Tweedie’s Identity does not hold for deep approximations of the MMSE. In order to be as exhaustive as possible, we detailed the addressed limitations of existing PnP methods in Appendix A.1.

heading: 3 THE GRADIENT STEP PLUG-AND-PLAY
text: The proposed method is based on the PnP version of half-quadratic-splitting (PnP-HQS) that amounts to replacing the proximity operator of the prior g with an off-the-shelf denoiser Dσ . In order to define a convergent PnP scheme, we first set up in Section 3.1 a Gradient Step (GS) denoiser. We then introduce the Gradient Step PnP (GS-PnP) algorithm in Section 3.2.

heading: 3.1 GRADIENT STEP DENOISER
text: We propose to plug a denoising operator Dσ that takes the form of a gradient descent step
Dσ = Id−∇gσ, (4) with gσ : Rn → R. Contrary to Romano et al. (2017), our denoiser exactly represents a conservative vector field. The choice of the parameterization of gσ is fundamental for the denoising performance. As already noticed in Salimans & Ho (2021), we experimentally found that directly modeling gσ as a neural network (e.g. a standard network used for classification) leads to poor denoising performance. In order to keep the strength of state-of-the-art unconstrained denoisers, we rather use
gσ(x) = 1
2 ||x−Nσ(x)||2, (5)
which leads to Dσ(x) = x−∇gσ(x) = Nσ(x) + JNσ (x)T (x−Nσ(x)), (6) where Nσ : Rn → Rn is parameterized by a neural network and JNσ (x) is the Jacobian of Nσ at point x. As discussed in Appendix A.2, the formulation (5) for gσ has been proposed in (Romano et al., 2017, Section 5.2) and (Bigdeli & Zwicker, 2017) for a distinct but related purpose, and not exploited for convergence analysis. Thanks to our definition (6) for Dσ , we can parameterize Nσ with any differentiable neural network architecture Rn → Rn that has proven efficient for image denoising. Although the representation power of the denoiser is limited by the particular form (6), we show (see Section 5.1) that such parameterization still yields state-of-the-art denoising results. We train the denoiser Dσ for Gaussian noise by minimizing the MSE loss function
L(Dσ) = Ex∼p,ξσ∼N (0,σ2I)[||Dσ(x+ ξσ)− x|| 2], (7)
or L(gσ) = Ex∼p,ξσ∼N (0,σ2I)[||∇gσ(x+ ξσ)− ξσ|| 2], (8)
when written in terms of gσ using equation (4). Remark 1. By definition, the optimal solution g∗σ ∈ arg minL is related to the MMSE denoiserD∗σ , that is, the best non-linear predictor of x given x + ξσ . Therefore, it satisfies Tweedie’s formula and ∇g∗σ = −σ2∇ log pσ (Efron, 2011) i.e. g∗σ = −σ2 log pσ + C, for some C ∈ R. Hence approximating the MMSE denoiser with a denoiser parameterized as (4) is related to approximating the logarithm of the smoothed image prior of pσ with − 1σ2 gσ . This relation was used for image generation with “Denoising Score Matching” by Saremi & Hyvarinen (2019); Bigdeli et al. (2020).

heading: 3.2 A PLUG-AND-PLAY METHOD FOR EXPLICIT MINIMIZATION
text: The standard PnP-HQS operator is TPnP-HQS = Dσ◦Proxτf , i.e. (Id−∇gσ)◦Proxτf when using the GS denoiser as Dσ . For convergence analysis, we wish to fit the proximal gradient descent (PGD) algorithm. We thus propose to switch the proximal and gradient steps and to relax the denoising step with a parameter λ ≥ 0. Our PnP algorithm with GS denoiser (GS-PnP) then writes
xk+1 = T τ,λ GS-PnP(xk) with T τ,λ GS-PnP = Proxτf ◦(τλDσ + (1− τλ) Id),
= Proxτf ◦(Id−τλ∇gσ). (9)
Under suitable conditions on f and gσ (see Lemma 1 in Appendix C), fixed points of the PGD operator T τ,λGS-PnP correspond to critical points of a classical objective function in IR problems
F (x) = f(x) + λgσ(x). (10)
Therefore, using the GS denoiser from equation (4) is equivalent to include an explicit regularization and thus leads to a tractable global optimization problem solved by the PnP algorithm. Our complete PnP scheme is presented in Algorithm 1. It includes a backtracking procedure on the stepsize τ that will be detailed in Section 4.2. Also, after convergence, we found it useful to apply an extra gradient step Id−λτ∇gσ in order to discard the residual noise brought by the last proximal step Proxτf .
4 CONVERGENCE ANALYSIS Algorithm 1: Plug-and-Play image restoration Param.: init. z0, λ > 0, σ ≥ 0, > 0, τ0 > 0,
K ∈ N∗, η ∈ (0, 1), γ ∈ (0, 1/2). Input : degraded image y. Output: restored image x̂. k = 0; x0 = Proxτf (z0); τ = τ0/η; ∆ > ; while k < K and ∆ > do
zk = λτDσ(xk) + (1− λτ)xk; xk+1 = Proxτf (zk); if F (xk)− F (xk+1) < γτ ||xk − xk+1||
2 ; then τ = ητ ; else ∆ = F (xk)−F (xk+1)F (x0) ; k = k + 1 ;
end x̂ = λτDσ(xK) + (1− λτ)xK ;
In this section, we introduce conditions on f and Dσ that will ensure the convergence of the PnP iterations (9) towards a solution of (10). For that purpose, we make use of the literature of convergence analysis (Attouch et al., 2013; Beck, 2017; Beck & Teboulle, 2009) of the PGD algorithm in the nonconvex setting.

heading: 4.1 CONVERGENCE RESULTS
text: A common setting for image restoration is the convex smooth L2 data-fidelity f(x) = ||Ax− y||2. In order to cover the noiseless case or to deal with a broader range of common degradation models, like Laplace or Poisson noise, we only assume f to be proper, lower semicontinous and convex. Next, the regularizer gσ is assumed to be differentiable with Lipschitz gradient, but not necessarily convex. This assumption on gσ is reasonable from a practical perspective. Indeed, using a network Nσ with differentiable activation functions, the function gσ introduced in Section 3 is differentiable with Lipschitz gradient (details and proof are given in Appendix B). Without further assumptions on f and gσ , the following theorem establishes the convergence of both the objective function values and the residual for a large variety of IR tasks. Theorem 1 (Proof in Appendix C). Let f : Rn → R ∪ {+∞} and gσ : Rn → R be proper lower semicontinous functions with f convex and gσ differentiable with L-Lipschitz gradient. Let λ > 0, F = f + λgσ and assume that F is bounded from below. Then, for τ < 1λL , the iterates xk given by the iterative scheme (9) verify
(i) (F (xk)) is non-increasing and converges.
(ii) The residual ||xk+1 − xk|| converges to 0.
(iii) All cluster points of the sequence (xk) are stationary points of (10). Remark 2. In the nonconvex setting, the quantity γk = min0≤i≤k ||xi+1 − xi||2 is commonly used to analyze the convergence rate of the algorithm (Beck & Teboulle, 2009; Ochs et al., 2014). Following the proof of Theorem 1 in Appendix C, we can obtain γk ≤ 1k F (x0)−limF (xk) 1 2τ− L 2 that is to say a O( 1k ) convergence rate for the squared L 2 residual. Remark 3. Even if most data-fidelity terms f encountered in image restoration are convex, Theorem 1 can be extended to nonconvex f . The proof of (iii) requires technical adaptations that can be found in Li & Lin (2015, Theorem 1) (as the 1-Lipschitz property of Proxτf does not hold anymore). Such nonconvex f appear for example in the context of phase retrieval (Metzler et al., 2018).
We can further obtain convergence of the iterates by assuming that the generated sequence (xk) is bounded and that f and gσ verify the Kurdyka-Lojasiewicz (KL) property. The boundedness of (xk) is discussed in Appendix D. The KL property (defined in Appendix E) has been widely used to study the convergence of optimization algorithms in the nonconvex setting (Attouch et al., 2010; 2013; Ochs et al., 2014). Very large classes of functions, in particular all the semi-algebraic functions, satisfy this property. In practice, in the extent of our analysis, the KL property is always satisfied.
Theorem 2 (Proof in Attouch et al. (2013), Theorem 5.1). Let f : Rn → R ∪ {+∞} and gσ : Rn → R be proper lower semicontinous functions with f convex and gσ differentiable with L-Lipschitz gradient. Let λ > 0, F = f + λgσ and assume that F is bounded from below. Assume that F verify the KL property. Suppose that τ < 1λL . If the sequence (xk) given by the iterative scheme (9) is bounded, then it converges, with finite length, to a critical point x∗ of F .
Remark 4. As explained in Attouch et al. (2013, Remark 5.2), the continuity of f is not required since we use the “exact forward-backward splitting algorithm”. We can thus deal with non-continuous data-fidelity terms, as it is the case with the inpainting application in Appendix J.3.
A more detailed description of all the assumptions of Theorems 1 and 2 is given in Appendix F.
4.2 BACKTRACKING TO HANDLE THE LIPSCHITZ CONSTANT OF ∇gσ
The convergence of Algorithm 1 actually requires to control the Lipschitz constant of ∇gσ only on a small subset of images related to {xk}. Therefore, estimating L for all images and setting the maximum stepsize τ accordingly will lead to sub-optimal convergence speed. In order to avoid small stepsizes, we use the backtracking strategy of Beck (2017, Chapter 10) and Ochs et al. (2014).
The convergence study in the proof of Theorem 1 is based on the sufficient decrease property of F established in equation (33). Without knowing the exact Lipschitz constant L, backtracking aims at finding the maximal stepsize τ yielding the sufficient decrease property. Given γ ∈ (0, 1/2), η ∈ [0, 1) and an initial stepsize τ0 > 0, the following update rule on τ is applied at each iteration k:
while F (xk)− F (T τ,λGS-PnP(xk)) < γ
τ ||T τ,λGS-PnP(xk)− xk|| 2, τ ←− ητ. (11)
Proposition 1 (Proof in Appendix G). Under the assumptions of Theorem 1, at each iteration of the algorithm, the backtracking procedure (11) is finite (i.e. a stepsize satifying F (xk)− F (T τ,λGS-PnP(xk)) ≥ γ τ ||T τ,λ GS-PnP(xk)− xk||2 is found in a finite number if iterations), and with backtracking, the convergence results of Theorem 1 and Theorem 2 still hold.
Remark 5. In practice, as explained in Section 5, we choose to initialize the stepsize as λτ0 = 1 and backtracking hardly ever activates. In this particular case, our results prove convergence of the standard PnP-HQS scheme xk+1 = Proxτf ◦Dσ(xk) applied with our specific denoiser.

heading: 5 EXPERIMENTS
text: In this section, we first show the performance of the GS denoiser. Next we empirically confirm that our GS-PnP method is convergent while providing state-of-the art results for different IR tasks.

heading: 5.1 GRADIENT-DESCENT-BASED DENOISER
text: Denoising Network Architecture We choose to parameterize Nσ with the architecture DRUNet (Zhang et al., 2021)) (represented in Appendix H), a U-Net in which residual blocks are integrated. One first benefit of DRUNet is that it is built to take the noise level σ as input, which is consistent with our formulation. Also, the U-Net models have previously offered good results in the context of prior approximation via Denoising Score Matching (Ho et al., 2020). Furthermore, Zhang et al. (2021) showed that DRUNet yields state-of-the-art results for denoising but also for PnP image restoration. In order to ensure differentiability w.r.t the input, we change RELU activations to ELU. We also limit the number of residual blocks to 2 at each scale to lower the computational burden.
Training details We use the color image training dataset proposed in Zhang et al. (2021) i.e. a combination of the Berkeley segmentation dataset (CBSD) (Martin et al., 2001), Waterloo Exploration Database (Ma et al., 2017), DIV2K dataset (Agustsson & Timofte, 2017) and Flick2K dataset (Lim et al., 2017). During training, the input images are corrupted with a white Gaussian noise ξσ with standard deviation σ randomly chosen in [0, 50/255]. With our parameterization (6) of Dσ , the network Nσ is trained to minimize the L2 loss (7). While the original DRUNet is trained with L1 loss, we stick to the L2 loss to keep the interpretability of gσ as an approximation of the log prior (see Remark 1). For each batch, the gradient∇gσ is computed with PyTorch differentiation tools.
We train the model on 128 × 128 patches randomly sampled from the training images, with batch size 16, during 1500 epochs. We use the ADAM optimizer with learning rate 10−4, divided by 2 every 300 epochs. It takes around one week to train the model on a single Tesla P100 GPU.
Denoising results We evaluate the PSNR performance of the proposed GS denoiser (GS-DRUNet) on 256× 256 color images center-cropped from the original CBSD68 dataset images (Martin et al., 2001). In Table 1, we compare, for various noise levels σ, our model with the simplified DRUNet (called “DRUNet light”) that has the same architecture as our GS-DRUNet (2 residual blocks) and that is trained (with L2 loss) without the conservative field constraint. We also provide comparisons with the original DRUNet (Zhang et al., 2021) (with 4 residual blocks at each scale and trained with L1 loss) and two state-of-the-art denoisers encountered in the PnP literature: FFDNet (Zhang et al., 2018) and DnCNN (Zhang et al., 2017a). For each network, we indicate in Table 1 the average runtime while processing a 256× 256 color image on one Tesla P100 GPU.
Our GS-DRUNet denoiser, despite being constrained to be an exact conservative field, reaches the performance of (and even slightly outperforms) its unconstrained counterpart DRUNet light. Second, departing from the latter, we are able to reduce the processing time by a large margin (÷7) while keeping close PSNR to the original DRUNet (around -0.05dB) and maintaining a significant PSNR gap (around +0.5dB) with other deep denoisers like DnCNN and FFDNet. Note that the time difference between GS-DRUNet and DRUNet light is due to the computation of ∇gσ via backpropagation. These results indicate that GS-DRUNet is likely to yield a competitive and fast PnP algorithm.

heading: 5.2 PLUG-AND-PLAY IMAGE RESTORATION
text: We show in this section that, in addition to being convergent, our PnP Alg. 1 reaches state-of-the-art performance in deblurring and super-resolution (Sections 5.2.1 and 5.2.2) and realizes relevant inpainting (Appendix J.3). In all cases, we seek an estimate x of a clean image x̂ ∈ Rn, from an observation obtained as y = Ax̂+ ξν ∈ Rm, with A a m×n degradation matrix and ξν a white Gaussian noise with zero mean and standard deviation ν. The objective function minimized by Alg. 1 is
F (x) = 1 2ν2 ||Ax− y||2 + λ 2 ||Nσ(x)− x||2. (12)
In practice, for ν > 0, we multiply F in equation (12) by ν2 and consider F = f + λνgσ with f(x) = 12 ||Ax− y|| 2, gσ(x) = 12 ||Nσ(x)− x|| 2 and λν = λν2. With this formulation, the convergences of iterates and objective values are guaranteed by Theorems 1 and 2. We also demonstrate in Appendix J.3 that our framework can be extended to other kinds of objective functions. For example, inpainting noise-free input images leads to a non differentiable data-fidelity term f .
Due to the large computational time of some compared methods, we use for evaluation and comparison a subset of 10 color images taken from the CBSD68 dataset (CBSD10) together with the 3 famous set3C images (butterfly, leaves and starfish). Quantitative results run on the full CBSD68 dataset are given in Appendix J. All images are center-cropped to the size 256× 256. For each IR problem, we provide default values for the parameters σ and λ that can be used to treat sucessfully a large class of images and degradations. The influence of both parameters is analyzed in Appendix J.5. Performance can be marginally improved by tuning λ for each image, for example with the method of Wei et al. (2020) based on reinforcement learning. In our experiments, backtracking is performed with η = 0.9 and γ = 0.1. We observe (see Appendix B) that on a majority of images, the Lipschitz constant L of ∇gσ is slightly larger than 1. As convergence is ensured for λντ = ν2λτ < 1L , we set the initial stepsize to τ0 = (ν2λ)−1. At the first iteration, the gradient step in equation (9) is thus exactly Dσ . In the majority of our experiments, backtracking is never activated. The algorithm is initialized with a proximal step and terminates when the relative difference between consecutive values of the objective function is less than or the number of iterations exceeds K.

heading: 5.2.1 DEBLURRING
text: For image deblurring, the degradation operator A = H is a convolution performed with circular boundary conditions. Therefore, H = F∗ΛF , where F is the orthogonal matrix of the discrete Fourier transform (and F∗ its inverse), and Λ is a diagonal matrix. The proximal operator of the data-fidelity term f(x) = 12 ||Hx− y|| 2 involves only element-wise inversion and writes
Proxτf (z) = F∗(In + τΛ∗Λ)−1F(τHT y + z). (13) We demonstrate the effectiveness of our method on a large variety of blur kernels (represented in Table 2) and noise levels. As in (Zhang et al., 2017b; Pesquet et al., 2021; Zhang et al., 2021), we use the 8 real-world camera shake kernels proposed in Levin et al. (2009) as well as the 9 × 9 uniform kernel and the 25×25 Gaussian kernel with standard deviation 1.6 (as in (Romano et al., 2017)). We consider Gaussian noise with 3 noise levels ν ∈ {2.55, 7.65, 12.75}/255 i.e. ν ∈ {0.01, 0.03, 0.05}. For all noise levels, we set σ = 1.8ν, λν = ν2λ = 0.1 for motion blur (kernels (a) to (h)) and λν = 0.075 for static blur (kernels (i) and (j)). Initialization is done with z0 = y but we show in Appendix J.6 the robustness to the initialization. The stopping criteria are = 10−5 and K = 400.
We compare in Table 2 our method (GS-PnP) against the patch-based method EPLL (Zoran & Weiss, 2011; Hurault et al., 2018), the deep PnP methods IRCNN (Zhang et al., 2017b), DPIR (Zhang et al., 2021), MMO (Pesquet et al., 2021), and the “RED-FP” algorithm (Romano et al., 2017) (with TNRD denoiser (Chen & Pock, 2016)) referred to as RED. Both IRCNN and DPIR use PnP-HQS with a fast decrease of τ and σ in a few iterations (8 iterations for DPIR) without guarantee of convergence. DPIR uses the DRUNet denoiser from Table 1. MMO is the only compared method that guarantees convergence by plugging a DnCNN denoiser trained with Lipschitz constraints (but the only given network was trained for very low noise level). Finally, as RED only treats the Y channel in the YCbCr color space, we also indicate in Appendix J.1, for RED and the proposed method, the PSNR evaluated on the Y channel only.
Among all methods, GS-PnP closely follows DPIR in terms of PSNR for low noise level but performs equally or better for higher noise levels. Other comparisons are conducted in Appendix J.1 on the Set3c and the full CBSD68 datasets (Tables 4 and 5). These results exhibit that GS-PnP reaches state-of-the-art in PnP deblurring for a variety of kernels and noise levels. We underline that the convergence of GS-PnP is guaranteed, whereas DPIR can asymptotically diverge (see Appendix J.7).
For qualitative comparison, we show in Figure 1(c-f) the deblurring obtained with various methods on the image “starfish” (from set3C). Note that our algorithm, compared to competing methods, can recover the sharpest edges. We also give convergence curves that empirically confirm the convergence of the values F (xk) (g) and of the residual min0≤i≤k ||xi+1 − xi||2/||x0||2 (h). These observations are supported by the additional experiment shown in Appendix J.1, Figure 6.

heading: 5.2.2 SUPER-RESOLUTION
text: For single image super-resolution, the low-resolution image y ∈ Rm is obtained from the highresolution one x ∈ Rn via y = SHx + ξν where H ∈ Rn×n is the convolution with anti-aliasing kernel. The matrix S is the standard s-fold downsampling matrix of size m × n and n = s2 ×m. In this context, we make use of the closed-form calculation of the proximal map for the data-fidelity term f(x) = 12 ||SHx− y|| 2, given by Zhao et al. (2016):
Proxτf (z) = ẑτ − 1 s2 F∗Λ∗
( Im + τ
s2 ΛΛ∗
)−1 ΛF ẑτ , (14)
where ẑτ = τHTST y + z and Λ = [Λ1, . . . ,Λs2 ] ∈ Rm×n, with Λ = diag(Λ1, . . . ,Λs2) a blockdiagonal decomposition according to a s × s paving of the Fourier domain. Note that Im + τdΛΛ ∗ is a m × m diagonal matrix and its inverse is computed element-wise. As expected, with s = 1, equation (14) comes down to equation (13).
As in Zhang et al. (2021), we evaluate super-resolution performance on 8 Gaussian blur kernels represented in Table 3: 4 isotropic kernels with different standard deviations (0.7, 1.2, 1.6 and 2.0) and 4 anisotropic kernels. Results are averaged between isotropic and anisotropic. We consider downsampled images at scale s = 2 and s = 3 and Gaussian noise with 3 different noise levels ν ∈ {0.01, 0.03, 0.05}. Our method (GS-PnP) is compared against bicubic upsampling, RED, IRCNN (“IRCNN+” from (Zhang et al., 2021)) and DPIR. We give again in Appendix J.2 the results obtained on the Set3C dataset (Table 7). All our results are obtained with λν = ν2λ = 0.065 and σ = 2ν. Initialization z0 is done with a bicubic interpolation of y (with a shift correction (Zhang et al., 2021)) and the stopping criteria are = 10−6 and K = 400.
Besides being the only compared PnP method with convergence guarantee, GS-PnP outperforms in PSNR all other PnP algorithms over the considered range of blur kernels, noise levels and scale factors. We show in Figure 2 the super-resolution of the image “leaves” downsampled by 2, with an isotropic kernel and noise level ν = 0.03. GS-PnP (f) recovers more accurately structures and color details than competing approaches (c-e), while converging in terms of function values (g) and residual (h). Additional visual comparisons are presented in Appendix J.2.

heading: 6 CONCLUSION
text: In this work, we introduce a new PnP algorithm with convergence guarantees. A denoiser is trained to realize an exact gradient step on a regularization function that is formulated through a neural network. This denoiser is plugged in an iterative scheme closely related to PnP-HQS, which is proved to converge towards a stationary point of an explicit functional. One strength of this approach is to simultaneously allow for a non strongly convex (and non smooth) data-fidelity term with a denoiser that may not be nonexpansive. Experiments conducted on ill-posed imaging problems (deblurring, super-resolution, inpainting) confirm the convergence results and show that the proposed PnP algorithm reaches state-of-the-art image restoration performance. This work also opens several research perspectives. One could first examine which information is encoded in the proposed prior. For example, based on the sharp visual results, one can question if a relation can be drawn between this prior and the gradient energy or sparsity. Also, it would be interesting to see if the recovery analysis developed in (Liu et al., 2021) can adapt to the proposed framework.

heading: 7 REPRODUCIBILITY STATEMENT
text: Anonymous source code is given in supplementary material. It contains a README.md file that explains step by step how to run the algorithm and replicate the results of the paper. Moreover, the pseudocode of our algorithm is given Algorithm 1. In Section 5 it is precisely detailed how all the hyper-parameters are chosen and, for each experiment, which dataset is used. As for the theoretical results presented in Section 4, complete proofs are given in the appendixes.

heading: 8 ETHICS STATEMENT
text: We believe that this work does not raise potential ethical concerns.

heading: ACKNOWLEDGEMENTS
text: This work was funded by the French ministry of research through a CDSN grant of ENS ParisSaclay. This study has also been carried out with financial support from the French Research Agency through the PostProdLEAP and Mistic projects (ANR-19-CE23-0027-01 and ANR-19-CE40-005).

heading: A COMPARISON WITH THE PNP LITERATURE
text: A.1 LIMITATIONS OF PREVIOUS PNP METHODS
In the existing literature, PnP approaches have one of the following limitations:
• They are not able to provide proof of convergence when non strongly convex data-fidelity terms are involved (Ryu et al., 2019), which is the case of some classical IR problems such as deblurring, super-resolution or inpainting.
• They are restricted to (nearly) nonexpansive denoisers (Reehorst & Schniter, 2018; Ryu et al., 2019; Sun et al., 2019b; Xu et al., 2020) or denoisers with a symmetric Jacobian (Romano et al., 2017). But it has already been shown that imposing symmetric Jacobian or Lipschitz property on a deep denoiser network alters its denoising performance (Bohra et al., 2021; Hertrich et al., 2021). We highlight that it was already empirically observed (Romano et al., 2017; Zhang et al., 2021) that the performance of the denoiser directly impacts the performance of the corresponding PnP scheme for IR.
• They show convergence of iterates thanks to decreasing time steps (Chan et al., 2016), but there is no characterization of the obtained solution (it is not a minima or a critical point of any functional).
On the other hand, our method is proved to converge to a stationary point of an explicit functional including a non strongly convex data-fidelity term. It also relies on a (possibly expansive) denoiser that, although being constrained to be a conservative vector field, allows to produce state-of-the-art results for various ill-posed IR problems.
A.2 ON THE REGULARIZATION gσ .
We first underline that the main point of our method is to define the denoiser as Dσ = Id−∇gσ . The choice for gσ is important for the denoising performance. With respect to the convergence properties of GS-PnP, this is nevertheless a secondary issue, as our method would converge for other differentiable regularizers gσ .
The proposed regularization gσ(x) = 12 ||x−Nσ(x)|| 2 was previously mentioned in the RED original paper (Romano et al., 2017) (but explicitly left aside) and used in the DAEP paper (Bigdeli & Zwicker, 2017). The main difference between our regularizer and the one alternately proposed in RED and DAEP is the following:
• RED and DAEP both consider a generic given pretrained denoiser Dσ : Rn → Rn, which is then associated with the regularizer gσ(x) = 12 ||x − Dσ(x)||
2 and used as such in IR problems.
• In our method, we set gσ(x) = 12 ||x−Nσ(x)|| 2 (with Nσ : Rn → Rn differentiable) and
then we train the denoiser as Dσ = Id−∇gσ with the loss function ||Dσ(x+ ξ)− x||2 for clean images x and additive white Gaussian noise (AWGN) ξ.
With this new formulation, we are ensured that Dσ = Id−∇gσ is inherently a conservative vector field, without further assumptions on Nσ . Thanks to this relation, the (slightly modified) PnP-HQS given in relation (9) becomes a proximal gradient descent (PGD). We can then make use of convergence results of the PGD algorithm in the nonconvex setting to show the convergence of PnP-HQS.
In contrast to the original RED paper, we aimed at finding one setting of Plug-and-Play image restoration that allows for a convergence proof with sufficiently general hypotheses. For this purpose, we had to consider this very particular form of regularization.
A.3 RECENT LITERATURE ON REGULARIZATION BY DENOISING
We here provide a more detailed discussion on the follow-up literature on RED.
In parallel to the RED method (Romano et al., 2017), the authors of Bigdeli & Zwicker (2017) propose to use the regularization, mentioned but not exploited in RED, g(x) = ||Dσ(x) − x||2, where Dσ is a pretrained denoising autoencoder. Next, Bigdeli et al. (2017) extended this work with a new prior, which is the Gaussian-smoothed version of the natural image prior. Inspired by Tweedie’s formula, they approximate the gradient of this log smoothed prior with the residual of a pretrained denoising autoencoder. With this new formulation, it is possible to optimize on the restored image but also on other parameters (e.g. the noise level and the used blur kernel).
Initially designed in the context of convex data-fidelity term, RED Romano et al. (2017) has been applied in the nonconvex setting for phase retrieval problems in prDeep Metzler et al. (2018).
Regularization by Artifact-Removal (RARE) (Liu et al., 2020) extends the RED framework by replacing the denoiser by a more general artifact-removal operator. The main advantage of this operator is that it can be trained without groundtruth data, but only by mapping pairs of artifact and noise contaminated images obtained directly from undersampled measurements. This is particularly useful for medical imaging applications where it is difficult to acquire fully-sampled training data.
The convergence of the original RED algorithm is discussed in Reehorst & Schniter (2018). The authors provide a convergence proof for RED-PGD which requires the denoiser to be nonexpansive, which, as detailed in the previous sections, is a restrictive hypothesis.
The authors of Liu et al. (2021) provide a recovery guarantee for the PnP framework, meaning convergence to a x∗ that satisfies y = Ax∗ while being in the set Fix(D) of the fixed points of D. More precisely, they show the convergence of the PnP-PGD method towards such a true solution x∗ under the assumptions that the denoiser residual R = Id−D is bounded and Lipschitz, and that the measurement operator satisfies a “set-restricted eigenvalue condition” (S-REC, which can be understood as strong convexity on the image of the denoiser Im(D)). Under the additional assumptions that the denoiser is nonexpansive and that there exists x ∈ Fix(D) that is also critical for the regularizer g, they show that PnP and RED have the same solutions. As mentioned by the authors, it is nevertheless difficult to verify the S-REC condition for a given measurement operator: since Im(D) is not explicit, it is not clear how much S-REC relaxes the strong convexity. As explained in Sections 2 and 3, our results do not require strong convexity of the data-fidelity term.
Instead of including an explicit regularization in the functional, RED-PRO (Cohen et al., 2021) aims at minimizing the data-fidelity term on the set Fix(D) of fixed points of a generic denoiser D. The study is conducted under the hypothesis that the denoiser is demicontractive, which implies that Fix(D) is convex, thus leading to a convex optimization problem. However, this assumption seems difficult to verify in practice and the existence of fixed points for the RED-PRO operator does not appear straightforward. In contrast, the fixed points of the GS-PnP operator are directly related to the stationary points of the global functional F = f + λgσ (Lemma 1 in Appendix C), whose existence is guaranteed as soon as F is coercive (see the discussion in Appendix D).
ASYNC-RED (Sun et al., 2020) enables faster computation of RED by decomposing the inference into a sequence of partial (block-coordinate) updates on x which can be executed asynchronously in parallel over a multicore system. As in their previous work BC-RED (Sun et al., 2019a), the authors propose to further reduce the computational time by using only a random subset of measurements at every iteration. Convergence of ASYNC-RED is shown, provided the denoiser is nonexpansive. A possible future extension of our work is the integration of the ASYNC framework to accelerate GS-PnP for large scale imaging inverse problems. As our GS-PnP converges without assuming nonexpansiveness of the denoising operation, it would be interesting to see if one can adapt the GS-PnP convergence properties to an ASYNC-GSPnP algorithm.
B LIPSCHITZ CONSTANT OF ∇gσ
First, let us give a result which ensures that a large class of neural networks trained with differentiable activation functions have Lipschitz gradients with respect to the input image.
Proposition 2. LetH = hp◦ . . .◦h1 be a composition of differentiable functions hi : Rdi−1 → Rdi . Let us assume that for any i the differential map h′i is bounded and Lipschitz. Then H ′ is Lipschitz.
Proof. Let us denote Hi = hi ◦ . . . ◦ h1 (and by convention, H0 = Id). Let ‖h′i‖∞ be the best uniform bound on the operator norms ‖h′i(x)‖, x ∈ Rdi−1 (which is also the best Lipschitz constant of hi). Let us also denote ‖h′i‖Lip the Lipschitz constant of h′i. The chain rule gives that for any x, H ′(x) can be expressed as a composition of linear maps
H ′(x) = h′p(Hp−1(x))h ′ p−1(Hp−2(x)) . . . h ′ 1(x) (15)
Therefore, for any x, y,
H ′(x)−H ′(y) = p−1∑ i=0 h′p(Hp−1(x)) . . . h ′ i+2(Hi+1(x))h ′ i+1(Hi(x))h ′ i(Hi−1(y)) . . . h ′ 1(y) (16)
− h′p(Hp−1(x)) . . . h′i+2(Hi+1(x))h′i+1(Hi(y))h′i(Hi−1(y)) . . . h′1(y). (17)
We can thus bound the operator norms
‖H ′(x)−H ′(y)‖ ≤ p−1∑ i=0 ( ‖h′p(Hp−1(x)) . . . h′i+2(Hi+1(x))‖ (18)
‖h′i+1(Hi(x))− h′i+1(Hi(y))‖‖h′i(Hi−1(y)) . . . h′1(y)‖ ) . (19)
and thus
‖H ′(x)−H ′(y)‖ ≤ p−1∑ i=0 ( ∏ j 6=i+1 ‖h′j‖∞ ) ‖h′i+1‖Lip‖H ′i‖∞‖x− y‖ (20)
which concludes because the chain-rule ensures that ‖H ′i‖∞ ≤ ‖h′i‖∞ . . . ‖h′1‖∞.
Proposition 2 applies for a neural network obtained as a composition of fully-connected layers with ELU activation functions, that is, by composing functions of the form
h(x) = E(Ax+ b) (21)
where A is a matrix, b a vector and E is the element-wise ELU defined by
E(x)i = { xi if xi ≥ 0 exi − 1 if xi < 0 .
(22)
It is easy to see that E is differentiable and that E′ is 1-Lipschitz with ‖E′‖∞ ≤ 1. Therefore h′(x) = E′(Ax+ b)A (23)
is also bounded and Lipschitz.
Let us also mention that this proposition encompasses the case of U-nets which, in addition to composing fully-connected layers, also integrates skip-connections. For example, taking a skipconnection on a composition h3 ◦ h2 ◦ h1 amounts to defining
H(x) = h3 ( h2(h1(x)) , h1(x) ) . (24)
This can be simply rewritten H = h3 ◦ h̃2 ◦ h1 where h̃2(x) = ( h2(h1(x)), h1(x) ) . (25)
It is then clear that h̃2 has bounded Lipschitz differential as soon as h1 and h2 do.
The bound obtained in the proof of Proposition 2 is exponential in the depth of the neural network. We now provide some experiments showing that, in practice, the Lipschitz constant of ∇gσ does not explode. We show in Figure 3, for various noise levels σ, the distribution of the spectral norms ||∇2gσ(x)||S on the training image set X , estimated with power iterations. The computed value varies a lot across images. Hence approximating the Lipschitz constant of ∇gσ with L = maxx∈X ||∇2gσ(x)||S would lead to under-estimated stepsizes and slow convergence on most images. Backtracking solves this issue by finding at each iteration the optimal stepsize allowing sufficient decrease of the objective function.

heading: C PROOF OF THEOREM 1
text: We first remind that a function f : Rn −→ R ∪+∞ is proper if its domain dom(f) = {x ∈ R, f(x) < +∞} (26)
is non empty. Also, recall that f is lower semicontinuous at x∗ if lim infx→x∗ f(x) ≥ f(x∗).
Proof.
(i) For ease of notation, we consider λ = 1. The generalisation for any λ > 0 is straightforward by rescaling gσ (and L) accordingly. We denote the proximal gradient fixed point operator Tτ = Proxτf ◦(Id−τ∇xgσ), the objective function F = f + gσ and we introduce
Qτ (x, y) = gσ(y) + 〈x− y,∇gσ(y)〉+ 1
2τ ||x− y||2 + f(x). (27)
We have
arg min x Qτ (x, y) = arg min x
gσ(y) + 〈x− y,∇gσ(y)〉+ 1
2τ ||x− y||2 + f(x)
= arg min x
f(x) + 1
2τ ||x− (y − τ∇gσ(y)||2
= Proxτf ◦(Id−τ∇xgσ)(y) = Tτ (y).
(28)
By definition for the arg min, xk+1 = Tτ (xk)⇒ Qτ (xk+1, xk) ≤ Qτ (xk, xk). Moreover, with gσ being L-smooth, we have by the descent lemma, for any τ ≤ 1L and any x, y ∈ R n,
gσ(x) ≤ gσ(y) + 〈x− y,∇gσ(y)〉+ 1
2τ ||x− y||2, (29)
so that for every x, y ∈ Rn, Qτ (x, x) = F (x) and Qτ (x, y) ≥ F (x). (30)
Therefore, at iteration k,
F (xk+1) ≤ Qτ (xk+1, xk) ≤ Qτ (xk, xk) = F (xk). (31) (F (xk)) is thus non-increasing. Since F is lower-bounded, (F (xk)) thus converges to a limit F ∗.
(ii) Note that Qτ (xk+1, xk) ≤ Qτ (xk, xk) implies
f(xk+1) ≤ f(xk)− 〈xk+1 − xk,∇gσ(xk)〉 − 1
2τ ||xk+1 − xk||2. (32)
Using also (29) with stepsize 1L , we get
F (xk+1) = f(xk+1) + gσ(xk+1)
≤ f(xk)− 〈xk+1 − xk,∇gσ(xk)〉 − 1
2τ ||xk+1 − xk||2
+ gσ(xk) + 〈xk+1 − xk,∇gσ(xk)〉+ L
2 ||xk+1 − xk||2 = F (xk)− ( 1
2τ − L 2
) ||xk+1 − xk||2.
(33)
Summing over k = 0, 1, ...,m gives
m∑ k=0 ||xk+1 − xk||2 ≤ 1 1 2τ − L 2 (F (x0)− F (xm+1))
≤ 1 1 2τ − L 2
(F (x0)− F ∗) . (34)
Therefore, limk→∞ ||xk+1 − xk|| = 0.
(iii) We begin by the two following lemmas characterizing the proximal gradient descent operator Tτ = Proxτf ◦(Id−τ∇xgσ).
Lemma 1. With the assumptions of Theorem 1, for x∗ ∈ Rn, x∗ is a fixed point of the proximal gradient descent operator Tτ = Proxτf ◦(Id−τ∇xgσ), i.e. Tτ (x∗) = x∗, if and only if x∗ is a stationary point of problem (10), i.e. −∇gσ(x∗) ∈ ∂f(x∗).
Proof. By definition of the proximal operator, we have
Tτ (x ∗) = x∗ ⇔ x∗ = Proxτf ◦(Id−τ∇xgσ)(x∗)
⇔ x∗ − τ∇xgσ(x∗)− x∗ ∈ τ∂f(x∗) ⇔ −∇xgσ(x∗) ∈ ∂f(x∗).
(35)
Lemma 2. With the assumptions of Theorem 1, Tτ is 1 + τL Lipschitz.
Proof. Using the fact that for f convex, Proxτf is 1-Lipschitz (Bauschke & Combettes, 2011, Proposition 12.28), and by the Lipschitz property of∇xgσ ,
||Tτ (x)− Tτ (y)|| = ||Proxτf ◦(Id−τ∇xgσ)(x)− Proxτf ◦(Id−τ∇xgσ)(y)|| ≤ ||(Id−τ∇xgσ)(x)− (Id−τ∇xgσ)(y)|| ≤ (1 + τL)||x− y||.
(36)
Note that, by nonconvexity of gσ , the fixed point operator Tτ is not necessarily nonexpansive, but we can still show the convergence of the fixed-point algorithm towards a critical point of the objective function. We can now turn to the proof of (iii). Let x∗ be a cluster point of (xk)k≥0. Then there exists a subsequence (xkj )j≥0 converging to x ∗. We have ∀j ≥ 0,
||x∗ − Tτ (x∗)|| ≤ ||x∗ − xkj ||+ ||xkj − Tτ (xkj )||+ ||Tτ (xkj )− Tτ (x∗)|| ≤ (2 + τL)||x∗ − xkj ||+ ||xkj − Tτ (xkj )|| by Lemma 2.
(37)
Using (ii), the right-hand side of the inequality tends to 0 as j →∞. Thus ||x∗ − Tτ (x∗)|| = 0 and x∗ = Tτ (x ∗), which by Lemma 1 means that x∗ is a stationary point of problem (10).
D ON THE BOUNDEDNESS OF (xk)
In order to obtain convergence of the iterates, in Theorem 2 the generated sequence (xk) is assumed to be bounded. In the experiments (Section 5), we observed that under the rest of assumptions of Theorem 2, boundedness was always verified. A sufficient condition for the boundedness of the iterates is the coercivity of the objective function, that is, lim|x|→∞ F (x) = +∞ (because the non-increasing property gives F (xk) ≤ F (x0)). Similar to Laumont et al. (2021), we can constrain F to be coercive by choosing a convex compact set C ⊂ Rn where the iterates should stay and by adding an extra term to the regularization gσ:
ĝσ(x) = gσ(x) + 1
2 ||x−ΠC(x)||2 =
1 2 ||x−Nσ(x)||2 + 1 2 ||x−ΠC(x)||2 (38)
with ΠC the Euclidian projection on C. As gσ is differentiable, the gradient step becomes
(Id−τλ∇xĝσ)(x) = (Id−τλ∇xgσ) + τλ(x−ΠC(x)). (39)
In our experiments, we choose the compact set C as C = [−1, 2]n. In practice we observe that all the iterates always remain in C and that the extra regularization term ||x−ΠC(x)||2 is never activated. Therefore, we don’t present this technical adaptation in Algorithm 1 but we let the reader aware that boundedness of (xk) is not a limiting assumption.

heading: E KL PROPERTY
text: Definition 1. Kurdyka-Lojasiewicz (KL) property (taken from Attouch et al. (2010))
(a) A function f : Rn −→ R ∪ +∞ is said to have the Kurdyka-Lojasiewicz property at x∗ ∈ dom(f) if there exists η ∈ (0,+∞), a neighborhood U of x∗ and a continuous concave function ψ : [0, η) −→ R+ such that ψ(0) = 0, ψ is C1 on (0, η), ψ′ > 0 on (0, η) and ∀x ∈ U ∩ [f(x∗) < f < f(x∗) + η], the Kurdyka-Lojasiewicz inequality holds:
ψ′(f(x)− f(x∗))dist(0, ∂f(x)) ≥ 1. (40)
(b) Proper lower semicontinuous functions which satisfy the Kurdyka-Lojasiewicz inequality at each point of dom(∂f) are called KL functions.
This condition can be interpreted as the fact that, up to a reparameterization, the function is sharp i.e. we can bound its subgradients away from 0. A big class of functions that have the KL-property is given by real semi-algebraic functions. For more details and interpretations, we refer to Attouch et al. (2010) and Bolte et al. (2010).

heading: F ON THE ASSUMPTIONS OF THEOREMS 1 AND 2
text: In this section, we explicitly list and comment all the assumptions required by Theorems 1 and 2. These assumptions are standard in nonconvex optimization. We now detail why each assumption is verified for our plug-and-play image restoration algorithm.
Assumptions of Theorem 1:
• Data-fidelity term f : Rn → R ∪ {+∞} proper lower semicontinous and convex. This is a general assumption that includes most of the data-fidelity terms classically used in IR problems. Note that we do not require differentiability of f . Degradations with Gaussian, Poisson or Laplacian noise models fall into this hypothesis. As noticed in Remark 3, our results can even be easily extended to a nonconvex data-fidelity term f , which encompasses applications like phase retrieval Metzler et al. (2018). In practice, it is helpful to have f proximable, i.e. Proxf with closed-form formula. Otherwise, Proxf needs to be calculated at each iteration with an optimization algorithm.
• Regularization function gσ : Rn → R proper lower semicontinous and differentiable with L-Lipschitz gradient. We parametrize as gσ(x) = 12 ||x−Nσ(x)||
2 with Nσ a differentiable neural network. This assumption on gσ is thus reasonable from a practical perspective. Indeed, using a network Nσ with differentiable activation functions, our function gσ is differentiable with Lipschitz gradient (details and proof are given in Appendix B).
• Functional F = f + λgσ bounded from below. This is straightforward as all the terms are positive.
• The stepsize τ < 1λL . This is handled by backtracking (see Section 4.2).
Assumptions of Theorem 2:
• Assumptions of Theorem 1 • F verify the KL property. The KL property (defined in Appendix E) has been widely used to
study the convergence of optimization algorithms in the nonconvex setting (Attouch et al., 2010; 2013; Ochs et al., 2014). Very large classes of functions, in particular all the semialgebraic functions, satisfy this technical property. It encompasses all the data-fidelity and regularization terms encountered in inverse problems.
• The sequence (xk) given by the iterative scheme (9) is bounded. As discussed in Appendix D, the boundedness can be ensured with a potential additional projection at each iteration. This is just a theoretical guarantee, as we observed that such a projection is never activated in practice.

heading: G BACKTRACKING AND PROOF OF PROPOSITION 1
text: Before giving the proof of Proposition 1, we first point out that our backtracking line search is a classical Armijo-type backtracking strategy, already used for nonconvex optimization in (Beck, 2017, Chapter 10) or Ochs et al. (2014). Other procedures could be investigated in future work. For instance, Li & Lin (2015) uses a Barzilai-Borwein rule to initialize the backtracking line search. Scheinberg et al. (2014) and Calatroni & Chambolle (2019) have also proposed a backtracking strategy that allows for both decreasing and increasing of the stepsize.
We now give the proof of Proposition 1.
Proof. For a given stepsize τ , we showed in Appendix C, equation (33) that
F (xk)− F (Tτ (xk)) ≥ 1
2
( 1
τ − L
) ||Tτ (xk)− xk||2. (41)
Taking τ < 1−2γL , we get 1 2 ( 1 τ − L) > γ τ so that
F (xk)− F (Tτ (xk)) > γ
τ ||Tτ (xk)− xk||2. (42)
Hence, when τ < 1−2γL , the sufficient decrease condition equation (42) is satisfied and the backtracking procedure (τ ←− ητ ) must end. In the proof of Theorem 1, we can replace the sufficient decrease (33) by (42) and finish the proof with the same arguments. In the same way, in the proof of Theorem 2 given in (Attouch et al., 2013, Theorem 5.1), our sufficient decrease (42) replaces (Attouch et al., 2013, Equation (52)).
H DRUNET light ARCHITECTURE
The architecture of the DRUNet light denoiser of (Zhang et al. (2021)) is given in Figure 4.

heading: I EXPANSIVENESS OF THE DENOISER
text: As gσ is not necessarily convex, our GS-DRUNet denoiser Dσ = Id−∇gσ is not necessarily nonexpansive and neither is the gradient step Id−λτ∇gσ . This is not an issue as, unlike previous theoretical PnP studies (Terris et al., 2020; Reehorst & Schniter, 2018), our convergence results do not require a nonexpansive denoising step. To advocate that our method converges without this assumption, we show in Figure 5 the evolution of ||Dσ(xk+1)−Dσ(xk)||||xk+1−xk|| along the algorithm that was run to obtain the super-resolution results of Figure 7. In this experiment, backtracking did not get activated and stayed fixed at λτ = 1. The gradient step in the PGD algorithm was thus simply a denoising step Dσ = Id−λτ∇gσ . Note that the Lipschitz constant of Dσ goes above 1 but convergence is still observed as shown by the two convergence curves in Figure 7.

heading: J ADDITIONAL EXPERIMENTS
text: J.1 DEBLURRING
We give here additional image deblurring experiments. We first present the PSNR performance comparison on the Set3c dataset in Table 4. We also provide an evaluation of the 3 best methods (GS-PnP, DPIR and IRCNN) on the full CBSD68 dataset in Table 5. For fair comparison with RED, we also display in Table 6 the PSNR calculated on the Y channel only. An additional visual
comparison is finally shown in Figure 6. Details and comments are given in the corresponding captions.
J.2 SUPER-RESOLUTION
We also present additional super-resolution experiments. We realize a full PSNR performance comparison on the Set3c dataset Table 7. We show additional visual comparisons between methods Figure 7 and Figure 8. Details and comments are given in the corresponding captions.
J.3 INPAINTING (WITH NON-DIFFERENTIABLE DATA-FIDELITY TERM)
We now propose to apply our PnP scheme to image inpainting with the degradation model
y = Ax (43)
whereA is a diagonal matrix with values in {0, 1}. For inpainting, no noise is added to the degraded image. In this context, the data-fidelity term is the indicator function of A−1({y}) = {x | Ax = y}: f(x) = ıA−1({y}) (which, by definition, equals 0 on A−1({y}) and +∞ elsewhere). Despite being non differentiable, f still verifies the assumptions of Theorems 1 and 2 and convergence is theoretically ensured. The proximal map becomes the orthogonal projection ΠA−1({y})
Proxτf (x) = ΠA−1({y})(x) = Ay −Ax+ x (44)
In our experiments, the diagonal of A is filled with Bernoulli random variables with parameter p = 0.5. We run our PnP algorithm with σ = 10/255. Given the form of f , we do not use the backtracking strategy and keep a fixed stepsize. Even if we do not exactly know the Lipschitz constant of ∇gσ , we observed in Figure 5 that, for small noise, it was almost always estimated as slightly larger than 1. We thus choose λτ = 1 and empirically confirm convergence with this choice in follow-up experiments (see Figure 9). The algorithm is initialized with x0 = y + 0.5(Id−A)y (masked pixels with value 0.5) and terminates when the number of iterations exceeds K = 100. We found it useful to run the first 10 iterations of the algorithm at larger noise level σ = 50/255. As y does not have noise, we found preferable not to run the last extra gradient pass from Algorithm 1.
We show inpainting results on set3C images Figure 9. Our PnP restores the input images with high accuracy, including its small details. Furthermore, convergence of the residual at rate O( 1k ) is empirically confirmed.
J.4 PSNR CONVERGENCE CURVES
We first plot Figure 10 the evolution of the PSNR along the iterations of the PnP algorithm during the experiments of Figure 1 and Figure 2. This illustrates that the minimization of F coincides with the maximization of the PSNR, which supports the interest of the optimized functional F = f+λgσ .
J.5 INFLUENCE OF THE PARAMETERS
In this section we study more deeply the influence of the parameters involved in the GS-PnP algorithm. Three parameters are involved: the stepsize τ , the denoiser level σ and the regularization parameter λ.
• The stepsize τ is automatically tuned with backtracking and is not tweaked heuristically, contrary to other competing methods based on PnP-HQS.
• The first regularization parameter σ is linked to the used denoiser. • The second regularization parameter λ is introduced so as to target the objective function f + λgσ , which is the main purpose of our method. It is a classical formulation of inverse problems, and the trade-off parameter λ is usually tuned manually.
Thus, like PnP-HQS, we have two parameters that we are free to tune manually. One additional motivation for keeping both λ and σ as regularization parameters is to be able to use our PnP algorithm with noise-blind denoisers like DnCNN that are independent on σ. In practice, in our experiments, we first roughly estimated σ proportionally to the input noise level ν and tweaked λ more precisely. Note that, for each inverse problem, our parameters λ and σ are fixed for a large variety of kernels, images and noise levels ν. The parameters are not optimized for each image.
Figure 11 and Figure 12 respectively plot the average PSNR when deblurring the CBSD10 images with different values λν and σ/ν, and fixed ν = 0.03. Both parameters control the strength of the regularization. We observe that λν and σ have a similar influence on the output: for small λν or small σ, the regularization involved by the denoising pass is not sufficient to counteract the noise amplification done by the proximal steps with large τ (recall that when τ → ∞, Proxτf tends to the pseudo-inverse of A). On the contrary, as expected, increasing λ and σ tends to over-smooth the output result.
On a single image, we also vary the main parameters of both GS-PnP and RED (Figure 13). For fair comparison, the PSNR is computed on the luminance channel only. This experiment confirms that, when manually optimizing the parameters for both methods, the PSNR results obtained with GS-PnP and RED remain close, as already observed in Table 6.
J.6 INFLUENCE OF THE INITIALIZATION
In Figure 14, we examine the robustness of the method to the initialization. As can be seen on this experiment, the output image does not change much even for relatively large perturbation of the initialization. We thus observe a robustness to the initialization, both in terms of visual aspect and PSNR. We also observed that initializing with a uniform image does not change the output of the algorithm. We suggest that this robustness comes from the first proximal steps on the data-fidelity term (with a large τ ), which prevent the algorithm to be stuck in a poor local minimum. Note that the use of large τ in the beginning of the algorithm is possible thanks to the backtracking procedure.
J.7 CONVERGENCE OF DPIR (ZHANG ET AL. (2021))
In this section, we illustrate that, contrary to our method, the DPIR algorithm is not guaranteed to converge and can even easily diverge. In Figure 15, we plot the convergence curves of both DPIR and our GS-PnP when deblurring the “starfish” image degraded with a motion kernel and ν = 0.01. In the original DPIR paper Zhang et al. (2021), only 8 iterations are used with decreasing τ and σ. More precisely, σ decreases uniformly in log-scale from 49 to the input noise level ν, and τ is set proportional to σ2. In order to study the asymptotic behaviour of the method, we propose two strategies to run DPIR with 1000 iterations:
(i) Decreasing σ from 49 to ν over 1000 iterations instead of 8 (Figure 15, row 1). (ii) Decreasing σ from 49 to ν in 8 iterations, and then keep the last values of σ and τ for the
the remaining iterations (Figure 15, row 2). As illustrated by the plot of ∑ i≤k ||xi+1 − xi||2 (third column), DPIR fails to converge with both strategies, even if the residual ||xk+1 − xk||2 tends to decrease with the second strategy. This divergence also involves a loss of restoration performance in terms of PSNR (first column). On the other hand, as theoretically shown in this paper, the residual ||xk+1 − xk||2 with GS-PnP tends to 0 (reaches ∼ 10−13 before the activation of backtracking, versus ∼ 10−4 for DPIR) and its series converges.

Review summaries:
This paper makes an extension of the plug-and-play framework by formulating an explicit regularizer $g(x)=||x - N_\sigma(x)||_2^2$ whose gradient $\nabla g$ corresponds to the noise residual $x-D_\sigma(x)$. By replacing the proximal of regularizer with this gradient step denoiser, the authors proposed the GS-PnP algorithm based on the half quadratic splitting algorithm. Since the explicit regularizer is known, a convergence analysis is hence established by assuming the Lipschitz continuity of $\nabla g$.  This paper is closely related to the recent trend of learning a regularizer functional by using deep neural networks. The difference between the existing literature and this paper is that the former trains a deep neural network to directly output a scalar value, while the latter trains a neural network to output an image-size vector and then envelopes it with a $\ell_2$-norm. In fact, the proposed regularizer shares the same formulation as the one stated in Sec 5.2 *"An Alternative Prior"* in Romano's RED paper. However, the difference between the two is not clearly stated in the paper.  **Strength** 1. An extension of the PnP with explicit regularizer formulated. 2. (Non-convex) convergence analysis under the common assumptions of the Lipschitz continuity of $\nabla g$  3. Extensive validation on image denoising, deblurring, and super-resolution.  

This paper presents a new form of the plug-and-play (PnP) half-quadratic splitting algorithm with provable convergence.  In contrast to directly formulating the denoiser as previous works, they propose to parameterize the denoiser with a learnable score-based function. This builds a bridge between the recently emerged score-based generative model and plug-and-play methods. And even more surprisingly, they show such a new formulation within the PnP framework leads to a very strong theoretical convergence guarantee under mild realistic assumptions. They also demonstrate good empirical results on three image restoration tasks i.e., deblurring, super-resolution and inpainting, verifying the empirical convergence rate is usually faster than $\mathcal{O}(\frac{1}{k})$ -- the worst-case rate established theoretically. 

The paper considers a novel and interesting idea: designing a deep neural network denoiser that makes Plug-and-play priors (PnP) convergence analysis clean and simple, motivated by PnP-HQS [Zhang et al., 2017b] and regularization by denoising (RED) [Romano et al., 2017]. Existing works have proved the convergence of PnP and RED with contractive and nonexpansiven denoisers. It is an activate area for designing deep neural net denoisers that ensure PnP/RED convergence. Specifically, this paper focuses on designing denoisers that correspond to the gradient descent step in terms of the $L$-smooth regularizer function $g_\sigma (x)$. By using such denoising step within PnP followed by the proximal of data fidelity $f$, one can guarantee convergence via traditional non-convex optimization given the objective function $F(x)$.  Since the handy $F(x)$ is available, a standard backtracking scheme is used in this work to ensure convergence without needs for tracking the exact Lipschitz constant of $\nabla g_\sigma$. Finally, the performance and stability of the proposed method is evaluated over three image inverse problems such as deblurring, super-resolution and inpainting, with satisfactory results compared to existing methods based on PnP and RED.