Hoang Trieu
email: hoang.le@ens-lyon.fr

Vy Le ⋆

Audrey Repetti
email: a.repetti@hw.ac.uk

Nelly Pustelnik
email: nelly.pustelnik@ens-lyon.fr

PNN: From proximal algorithms to robust unfolded image denoising networks and Plug-and-Play methods

Keywords: Image denoising, image restoration, unrolled proximal algorithms, unfolded neural networks, inertial methods

A common approach to solve inverse imaging problems relies on finding a maximum a posteriori (MAP) estimate of the original unknown image, by solving a minimization problem. In this context, iterative proximal algorithms are widely used, enabling to handle non-smooth functions and linear operators. Recently, these algorithms have been paired with deep learning strategies, to further improve the estimate quality. In particular, proximal neural networks (PNNs) have been introduced, obtained by unrolling a proximal algorithm as for finding a MAP estimate, but over a fixed number of iterations, with learned linear operators and parameters. As PNNs are based on optimization theory, they are very flexible, and can be adapted to any image restoration task, as soon as a proximal algorithm can solve it. They further have much lighter architectures than traditional networks. In this article we propose a unified framework to build PNNs for the Gaussian denoising task, based on both the dual-FB and the primal-dual Chambolle-Pock algorithms. We further show that accelerated inertial versions of these algorithms enable skip connections in the associated NN layers. We propose different learning strategies for our PNN framework, and investigate their robustness (Lipschitz property) and denoising efficiency. Finally, we assess the robustness of our PNNs when plugged in a forward-backward algorithm for an image deblurring problem.

Introduction

Image denoising aims to find an estimate of an unknown image x ∈ R N , from noisy measurements z ∈ R N . The present contribution focuses on the Gaussian denoising problem z = x + b, [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] where b ∈ R M models an additive white Gaussian noise with standard deviation δ > 0. A common method to denoise z is to rely on a maximum a posteriori (MAP) approach, and to define the estimate x MAP ∈ R N as a minimizer of a penalized least-squares objective function. A general formulation of this problem is to find

x MAP = argmin x∈R N F(x), (2)
where

F (x) := 1 2 ∥x -z∥ 2 2 + νg(Dx) + ι C (x), (3)
C ⊂ R N is a closed, convex, non-empty constraint set, ν > 0 is a regularization parameter proportional to δ 2 , D : R N → R |F| is a linear operator mapping an image from R N to a feature space R |F| , and g : R |F| → (-∞, +∞] denotes a proper, lower-semicontinous, convex function. The function g and the operator D are chosen according to the type of images of interest. For instance, functions of choice for piece-wise constant images are those in the family of total variation (TV) regularizations [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], which can be expressed as an ℓ 1 (or an ℓ 1,2 ) norm composed with a linear operator performing horizontal and vertical finite differences of the image. More generally, D can be chosen as a sparsifying operator (e.g., wavelet transform [START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF][START_REF] Pustelnik | Wavelet-based image deconvolution and reconstruction[END_REF]), and g as a function promoting sparsity (e.g., ℓ 1 ). For any choice of g and D, the parameter ν > 0 is used to balance the penalization term (i.e., function g • D) with the data-fidelity term (i.e., least-squares function).

When C = R N , and for simple choices of g • D, (2) may have a closed form solution (see, e.g., [START_REF] Chierchia | The proximity operator repository[END_REF], and references therein). Otherwise, (3) can be minimized efficiently using proximal splitting methods [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Chambolle | An introduction to continuous optimization for imaging[END_REF]. The most appropriate algorithm will be chosen depending on the properties of g, D and C. For simple sets C and when g • D is proximable, Douglas-Rachford (DR) scheme [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] can be considered, alternating, at each iteration, between a proximity step on the sum of the leastsquare term and the indicator function, and a proximity step on g • D. However, when g • D is not proximable nor differentiable (e.g., TV penalization, or when D is a redundant wavelet transform), more advanced algorithms must be used, splitting all the terms in (3) to handle them separately. Such methods usually rely on the Fenchel-Rockafellar duality [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF]. On the one hand, some algorithms can evolve fully in the dual space, such as, e.g., ADMM [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite elements approximations[END_REF][START_REF] Fortin | Augmented Lagrangian Methods : Applications to the Numerical Solution of Boundary-Value Problems[END_REF][START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] or the dual-Forward-Backward (FB) [START_REF] Combettes | Dualization of signal recovery problems[END_REF][START_REF] Combettes | Dualization of signal recovery problems[END_REF]. Note however that ADMM requires the inversion of D ⊤ D. On the other hand, other algorithms can alternate between the primal and the dual spaces, namely primal-dual algorithms [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms[END_REF][START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF][START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, lipschitzian, and parallel-sum type monotone operators[END_REF].

During the last decade, the performances of proximal algorithms have been pushed to the next level by mixing them with deep learning approaches [START_REF] Venkatakrishnan | Plug-and-play priors for model based reconstruction[END_REF][START_REF] Lecun | Deep learning[END_REF][START_REF] Ongie | Deep learning techniques for inverse problems in imaging[END_REF] leading to PNNs that consists in unrolling optimisation algorithms over a fixed number of iterations [START_REF] Adler | Learned primal-dual reconstruction[END_REF][START_REF] Jiu | A deep primal-dual proximal network for image restoration[END_REF]. In the litterature, two classes of PNNs can be encountered: PNN-LO (i.e. PNN with learned Linear Operators) where the involved linear operators are learned or the PNN-PO (i.e. PNN with learned Proximal Operators) where the proximity operators are replaced by small neural networks, typically a Unet [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF]. PNNs have shown to be very efficient for denoising task [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF] and further good performances for image restoration problems including deconvolution [START_REF] Jiu | A deep primal-dual proximal network for image restoration[END_REF][START_REF] Jiu | Alternative design of deeppdnet in the context of image restoration[END_REF], magnetic resonance imaging [START_REF] Adler | Learned primal-dual reconstruction[END_REF], or computed tomography [START_REF] Savanier | Deep unfolding of the DBFB algorithm with application to ROI CT imaging with limited angular density[END_REF]. Note that in the context of PNN-PO, if only the denoiser is learned, and the algorithm is unrolled until convergence, then the resulting approach boils down to plug-and-play (PnP) or deep equilibrium [START_REF] Bai | Deep equilibrium models[END_REF][START_REF] Gilton | Deep equilibrium architectures for inverse problems in imaging[END_REF][START_REF] Zou | Deep equilibrium learning of explicit regularization functionals for imaging inverse problems[END_REF].

Contributions -We propose a unified framework to design denoising PNN-LO. The resulting PNNs encompass our previous work [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF], with a deeper analysis on performance and stability of the proposed PNNs. The proposed PNN architectures are obtained by enrolling proximal algorithms to obtain a MAP estimate of a denoising problem of the form of ( 2)- [START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF], that is equivalent to computing a proximity operator. Precisely, we introduce a generic framework to build such denoisers derived from two proximal algorithms: the dual-FB iterations and the primal-dual Chambolle-Pock (CP) iterations. The proposed global PNN architectures also includes skip connections either on the primal or on the dual domains, corresponding to inertial parameters for acceleration purpose. We investigate the robustness and the denoising performances of the proposed PNNs, for different training strategies. To evaluate the robustness of our PNNs, we evaluate both their Lipchitz constants and their performances when applied in a different denoising setting than the Gaussian denoising task for which they have been trained (i.e., when evaluated on both Poisson-Gauss and Laplace-Gauss denoising). Further, similarly to our preliminary work [START_REF] Repetti | Dual Forward-Backward unfolded network for flexible Plug-and-Play[END_REF], we inject the resulting denoising PNNs in a FB algorithm for solving an image deblurring problem, to obtain a plug-and-play algorithm.

Outline -The remainder of this paper is organized as follows. Section 2 focuses on MAP denoising estimates, with a recall of the considered iterative schemes. Section 3 is dedicated to the design of the proposed unified PNN architectures, relying on the algorithmic schemes presented in Section 2. We evaluate the performances of the proposed PNNs in Section 4. In particular, we compare the denoising performances of our PNNs to those of an established denoising network in terms of both reconstruction quality and robustness. We further present results when the proposed PNNs are used in a PnP framework for an image deblurring task. Finally, our conclusions are given in Section 5.

Notation -In the remainder of this paper we will use the following notations. An element of R N is denoted by x. For every n ∈ {1, . . . , N }, the n-th coefficient of x is denoted by x (n) . The spectral norm is denoted ∥ • ∥ S . Let C ⊂ R N be a closed, non-empty, convex set. The indicator function of C is denoted by ι C , and is equal to 0 if its argument belongs to C, and +∞ otherwise. Let x ∈ R N . The Euclidean projection of x onto C is denoted by P C (x) = argmin v∈C ∥v -x∥ 2 . Let ψ : R N → (-∞, +∞] be a convex, lower semicontinuous, proper function. The proximity operator of ψ at x is given by prox ψ (x) = argmin

v∈R N ψ(v) + 1 2 ∥v -x∥ 2 . The Fenchel-Legendre conjugate function of ψ is given by ψ * (x) = sup v∈R N v ⊤ x -ψ(v). When ψ = λ∥ • ∥ 1 is the ℓ 1 norm, then ψ * = ι B∞(0,λ) corresponds to the indicator function of the ℓ ∞ -ball centred in 0 with radius λ > 0, i.e., B ∞ (0, λ) = {x ∈ R N | (∀n ∈ {1, . . . , N }) -λ ≤ x (n) ≤ λ}.
Let f K Θ be a feedforward NN, with K layers and learnable parameters Θ. It can be written as a composition of operators (i.e., layers)

f Θ = L Θ K • • • • • L Θ 1 ,
where, for every k ∈ {1, . . . , K}, Θ k are the learnable parameters of the k-th layer L Θ k . The k-th layer is defined as

L Θ k : u ∈ R N k → η k (W k u + b k ), (4)
where η k is an (non-linear) activation function, W k is a linear operator, and b k is a bias. In [START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF], authors shown that most of usual activation function involved in NNs correspond to proximal operators.

2 Denoising (accelerated) proximal schemes

In this section, we introduce the two algorithmic schemes we will focus in this work, i.e., FB and CP. Although problem (2)-( 3) does not have a closed form solution in general, iterative methods can be used to approximate it. Multiple proximal algorithms can be used to solve [START_REF] Mallat | A wavelet tour of signal processing[END_REF] as detailed in the introduction section. In this section, we describe two schemes enabling minimizing (3): the FB algorithm applied to the dual problem of (2), and the primal-dual CP algorithm directly applied to (2). In addition, for both schemes we also investigate their accelerated versions, namely DiFB (i.e., Dual inertial Forward Backward, also known as FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Chambolle | On the Convergence of the Iterates of the Fast Iterative Shrinkage/Thresholding Algorithm[END_REF]), and ScCP (i.e., CP for strongly convex functions [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]). For sake of simplicity, in the remainder of the paper we refer to these schemes as FB, CP, DiFB, and ScCP, respectively. for k = 0, 1, . . .

D(i)FB

u k+1 = prox τ k (νg) * v k + τ k D P C (z -D ⊤ v k ) , v k+1 = (1 + ρ k )u k+1 -ρ k u k , (5)
where (u 0 , v 0 ) ∈ R |F| × R |F| and the step-size parameters, for every k ∈ N, τ k > 0 and ρ k ≥ 0. Note that when, for every k ∈ N, ρ k = 0, then algorithm (5) reduces to DFB.

The following convergence result applies.

Theorem 2.1 ( [START_REF] Combettes | Dualization of signal recovery problems[END_REF][START_REF] Chambolle | On the Convergence of the Iterates of the Fast Iterative Shrinkage/Thresholding Algorithm[END_REF]). Let (u k , v k ) k∈N be generated by [START_REF] Chierchia | The proximity operator repository[END_REF]. Assume that one of the following conditions is satisfied.

1. For every k ∈ N, τ k ∈ (0, 2/∥D∥ 2 S ), and ρ k = 0.

2. For every k ∈ N, τ k ∈ (0, 1/∥D∥ 2 S ), and ρ k = t k -1 t k+1 with t k = k+a-1 a and a > 2.

Then we have

x MAP = lim k→∞ P C (z -D ⊤ u k ), (6)
where x MAP is defined in (2)-( 3).

(Sc)CP -A second strategy to solve (2)-(3) consists in applying the (Sc)CP algorithm to problem [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The data-term being ζ-strongly convex with parameter ζ = 1, the accelerated CP [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], dubbed ScCP, can be employed. This algorithm reads

for k = 0, 1, . . .     x k+1 = P C µ k 1+µ k (z -D ⊤ u k ) + 1 1+µ k x k , u k+1 = prox τ k (νg) * u k + τ k D (1 + α k )x k+1 -α k x k , (7)
where x 0 ∈ R N and u 0 ∈ R |F| . Note that when, for every k ∈ N, α k = 1, then algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] reduces to standard iterations of the primal-dual CP algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], while when α k = 0, its leads to the classical Arrow-Hurwicz algorithm [START_REF] Arrow | Studies in linear and non-linear programming[END_REF].

The following convergence result applies.

Theorem 2.2 ( [16]

). Let (u k , x k ) k∈N be generated by [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Assume that (τ k ) k∈N and (µ k ) k∈N are positive sequences, and that one of the following conditions is satisfied.

1. For every k ∈ N, τ k µ k ∥D∥ 2
S < 1, and α k = 1.

For every

k ∈ N, α k = (1 + 2µ k ) -1/2 , µ k+1 = α k µ k , and τ k+1 = τ k α -1 k with µ 0 τ 0 ∥D∥ 2 S ≤ 1.
Then we have

x MAP = lim k→∞ x k , (8)
where x MAP is defined in (2)-(3).

Proposed unfolded denoising NNs

In this section we aim to design PNNs f Θ such that

f Θ (z) ≈ x, (9)
where x ∈ R N is an estimate of x. As discussed in the introduction, such an estimate can correspond to the penalized least-squares MAP estimate of x, defined as in ( 2)-(3).

Primal-dual building block iteration

The iterations described previously in ( 5) and (7) share a similar framework which yields:

for k = 0, 1, . . .     u k+1 = prox τ k (νg) * u k + τ k Dx k x k+1 = P C µ k 1+µ k (z -D ⊤ u k+1 ) + 1 1+µ k x k . (10)
On the one hand, this scheme is a reformulation of the Arrow-Hurwicz (AH) iterations, i.e., algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] with α k ≡ 0. On the other hand, for the limit case when µ k → +∞, the DFB (5) iterations are recovered. Further, the inertia step is activated either on the dual variable for DiFB

u k+1 ← (1 + ρ k )u k+1 -ρ k u k (11)
or on the primal variable for ScCP

x k+1 ← (1 + α k )x k+1 -α k x k . (12)
Based on these observations we propose a strategy to build denoising PNNs satisfying ( 9), reminiscent to our previous works [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF][START_REF] Repetti | Dual Forward-Backward unfolded network for flexible Plug-and-Play[END_REF], to unroll D(i)FB and (Sc)CP algorithms using a primal-dual perspective.

The result provided below aims to emphasize that each iteration of the primal-building block [START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF] can be viewed as the composition of two layers of feedforward networks, acting either on the image domain (i.e., primal domain R N ) or in the features domain (i.e., dual domain R |F| ). DFB, DiFB, CP, and ScCP, hold the same structure, with extra steps that can be assimilated to skip connections in the specific case of DiFB and ScCP, enabling to keep track of previous layer's outputs (see [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite elements approximations[END_REF] and ( 12)).

Proposition 3.1. Let z ∈ R N , x ∈ R N , u ∈ R |F| , and k ∈ N. Let L ν,Θ k,D ,D : R N ×R |F| → R |F| , defined as L ν,Θ k,D ,D (x, u) = η ν,k,D (W k,D x + V k,D u + b k,D ) , (13)
be a sub-layer acting on both the primal variable x and the dual variable u and returning a dual (D) variable. In (13) η ν,k,D : R |F| → R |F| is a fixed activation function with parameter ν > 0, and Θ k,D is a linear parametrization of the learnable parameters including W k,D :

R N → R |F| , V k,D : R |F| → R |F| and b k,D ∈ R |F| . Let L z,Θ k,P ,P : R N × R |F| → R N , defined as L z,Θ k,P ,P (x, u) = η k,P (W k,P x + V k,P u + b z,k,P ) , (14)
be a sub-layer acting on both the primal variable x and the dual variable u and returning a primal (P) variable. In (14) η k,P : R N → R N is a fixed activation functions, and Θ k,P is a linear parametrization of the learnable parameters including W k,P :

R N → R N , V k,P : R |F| → R N and b z,k,P ∈ R N .
Then, the k-th iteration of the joint formulation (10) can be written as a composition of two layers of the form of (4):

L z,Θ k : R N × R |F| → R N : (x k , u k ) → L z,Θ k,P ,P (x, L Θ k,D ,D (x k , u k )), ( 15
)
where Θ k is the combination of Θ k,P and Θ k,D , i.e., the linear parametrization of all learnable parameters for layer k. Proof. This result is obtained by noticing that, for every k ∈ N, the k-th iteration in algorithm [START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF] can be rewritten as [START_REF] Combettes | Dualization of signal recovery problems[END_REF], where the primal (P) operators are given by W

k,P = 1 1+µ k , V k,P = -µ k 1+µ k D ⊤ , b z,k,P = µ k 1+µ k z,

Arrow-Hurwicz unfolded building block

Our unrolled architectures rely on layer structures introduced in Proposition 3.1 where we allow the linear operator D to be different for each layer. For more flexibility, we also introduce, for every k ∈ {1, . . . , K}, operators D k,D : R N → R |F| and D k,P : R |F| → R N , to replace operators D and D ⊤ , respectively, to allow a possible mismatch between operator D and its adjoint D ⊤ . The resulting unfolding Deep Arrow-Hurwicz (DAH) building block is then given below:

f K,DAH z,Θ (x 0 , u 0 ) = L DAH z,ν,Θ K • • • • • L DAH z,ν,Θ 1 (x 0 , u 0 ), (16)
where, for every k ∈ {1, . . . , K},

u k = L ν,Θ k,D ,D (x k-1 , u k-1 ) x k = L z,Θ k,P ,P (x k-1 , u k ) L DAH z,ν,Θ k (x k-1 , u k-1 ) = (x k , u k ). with              W k,D = τ k D k,D , V k,D = Id, b z,k,D = 0, η ν,k,D = prox τ k (νg) * ,
and

             W k,P = 1 1+µ k , V k,P = -µ k 1+µ k D k,P , b z,k,P = µ k 1+µ k z, η k,P = P C .

Proposed unfolded strategies

In this section we describe four unfolded strategies for building denoising NNs as defined in [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF]. All strategies rely on Arrow-Hurwicz building block presented in Section 3.2.

• DDFB stands for Deep Dual Forward-Backward and it fits DAH when µ k → +∞.

f K,DDFB z,ν,Θ (x 0 , u 0 ) = L DDFB z,ν,Θ K • • • • • L DDFB z,ν,Θ 1 (x 0 , u 0 ), (17)
where, for every k ∈ {1, . . . , K},

u k = L ν,Θ k,D ,D (x k-1 , u k-1 ) x k = L z,Θ k,P ,P (x k-1 , u k ) (x k , u k ) = L DDFB z,ν,Θ k (x k-1 , u k-1 ). with              W k,D = τ k D k,D , V k,D = Id, b k,D = 0, η ν,k,D = prox τ k (νg) * , and
             W k,P = 0, V k,P = -D k,P , b z,k,P = z, η k,P = P C .
• DDiFB: stands for Deep Dual inertial Forward-Backward interpreted as a DDFB with skip connections and defined as:

f K,DDiFB z,ν,Θ (x 0 , u 0 ) = (x K , u K ) ( 18
)
where, for every k ∈ {1, . . . , K},

(x k , u k ) = L DDFB z,ν,Θ k (x k-1 , u k-1 ) u k = (1 + ρ k ) u k -ρ k u k-1 .
• DCP stands for Deep Chambolle-Pock relying on DAH with a special update of the primal variable leading to:

f K,DCP z,ν,Θ (x 0 , u 0 ) = (x K , u K ) (19)
where, for every k ∈ {1, . . . , K},

( x k , u k ) = L DAH z,ν,Θ k (x k-1 , u k-1 ) x k = 2 x k -x k-1 .
• DScCP stands for Deep Strong convexity Chambolle-Pock interpreted as a DAH with skip connections on the primal variable:

f K,DScCP z,ν,Θ (x 0 , u 0 ) = (x K , u K ) ( 20
)
where, for every k ∈ {1, . . . , K},

( x k , u k ) = L DAH z,ν,Θ k (x k-1 , u k-1 ) x k = (1 + α k ) x k -α k x k-1 .
Illustration of a single layer of the resulting DD(i)FB and D(Sc)CP architectures are provided in Figure 1.

Since these architectures are reminiscent of D(i)FB and (Sc)CP, given in Section 2, we can deduce limit cases for the proposed unfolded strategies, when K → +∞ and linear operators are fixed over the layers.

Corollary 1 (Limit case for deep unfolded NNs). We consider the unfolded NNs DD(i)FB and D(Sc)CP defined in Section 3.3. Assume that, for every k ∈ {1, . . . , K}, D k,D = D and D k,P = D ⊤ , for D : R N → R |F| . In addition, for each architecture, we further assume that, for every k ∈ {1, . . . , K},

• DDFB: τ k ∈ (0, 2/∥D∥ 2 S ).
• DDiFB:

τ k ∈ (0, 1/∥D∥ 2 S ) and ρ k = t k -1 t k+1 with t k = k+a-1 a
and a > 2.

• DCP: (τ k , µ k ) ∈ (0, +∞) 2 such that τ k µ k ∥D∥ 2 S < 1.
• DScCP:

α k = (1 + 2µ k ) -1/2 , µ k+1 = α k µ k , and τ k+1 = τ k α -1 k with τ 0 µ 0 ∥D∥ 2 S ≤ 1.
Then, we have x K → x when K → +∞, where x K is the output of either of the unfolded NNs DD(i)FB or D(Sc)CP, and x is a solution to (2)

Table 1: Learnable paramaters of each PNN scheme

Θ k Comments DDFB-LFO D k,P , D k,D absorb τ k in D k,D DDiFB-LFO D k,P , D k,D , α k fix α k , and absorb τ k in D k,D DDFB-LNO D k,P = D ⊤ k,D define τ k = 1.99∥D k ∥ -2 DDiFB-LNO D k,P = D ⊤ k,D fix α k = t k -1 t k+1 , t k+1 = k+a-1 a , a >
µ k+1 = α k µ k DCP-LNO D k,P = D ⊤ k,D , µ learn µ = µ0 = • • • = µK , and fix τ k = 0.99µ -1 ∥D k ∥ -2 DScCP-LNO D k,P = D ⊤ k,D , µ k learn µ k , and fix α k = (1 + 2µ k ) -1/2 , and τ k = 0.99µ -1 k ∥D k ∥ -2

Proposed learning strategies

The proposed DD(i)FB and D(Sc)CP allow for flexibility in the learned parameters. In this work, we propose two learning strategies, either satisfying conditions described in Corollary 1 (LNO), or giving flexibility to the parameters (LFO). These two strategies are described below. The learnable parameters are summarized in Table 1. Training strategy for unfolded denoising networks -The network parameters are optimized by minimizing the ℓ 2 empirical loss between noisy and ground-truth images:

Experiments

Θ ∈ Argmin Θ 1 |I| s∈I L(x s , z s ; Θ) (21)
where

L(x s , z s ; Θ) := 1 2 ∥x s -f K zs,δ 2 ,Θ (z s , D 1,D (z s )) ∥ 2 ,
and f K zs,δ 2 ,Θ is either of the unfolded networks described in Section 3.3. The loss (21) will be optimized in Pytorch with Adam algorithm [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. For the sake of simplicity, in the following we drop the indices in the notation of the network f K zs,δ 2 ,Θ and use the notation f Θ .

Architectures -We will compare the different architectures introduced in Section 3.3, namely DDFB, DDiFB, DDCP, and DScCP, considering both LNO and LFO learning strategies (see Table 1 for details). For each architecture and every layer k ∈ {1, . . . , K}, the weight operator D k,D consists of J convolution filters (features), mapping an image in R N to features in R |F| , with |F| = J N . For LFO strategies, weight operator D k,P consists of J convolution filters mapping from R |F| to R N . All convolution filters considered in our work have the same kernel size of 3 × 3.

We evaluate the performance of the proposed models, varying the numbers of layer K and the number of convolution filter J. In our experiments we consider g = ∥ • ∥ 1 leading to HardTanh activation function as in [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF] and recalled below.

Proposition 4.1. The proximity operator of the conjugate of the ℓ 1 -norm scaled by parameter ν > 0 is equivalent to the HardTanh activation function, i.e., for every x = (x (n) ) 1≤n≤N :

(p (n) ) 1≤n≤N = prox (ν∥•∥ 1 ) * (x) = P ∥•∥∞≤ν (x) = HardTanh ν (x) where p (n) =        -ν if x (n) < -ν, ν if x (n) > ν, x (n) otherwise.
Experimental settings -To evaluate and compare the proposed unfolded architectures, we consider 2 training settings. In both cases, we consider RGB images (i.e., C = 3).

• Training Setting 1 -Fixed noise level: The PNNs are trained with |I| = 200 images extracted from BSDS500 dataset [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF], with a fixed noise level δ = 0.08. The learning rate for ADAM is set to 8 × 10 -5 (all other parameters set as default), we use batches of size of 10 and patches of size 50 × 50 randomly selected.

We train the proposed unfolded NNs considering multiple sizes of convolution filters J ∈ {8, 16, 32, 64} and for multiple numbers of layers K ∈ {5, 10, 15, 20, 25}. All PNNs are trained with (J, K) = (64, 20).

In our experiments, we aim to compare the proposed PNNs for three metrics: (i) architecture complexity, (ii) robustness, and (iii) denoising performance. For sake of completeness, these metrics will also be provided for a state-of-the-art denoising network, namely DRUnet [START_REF] Zhang | Plug-and-play image restoration with deep denoiser prior[END_REF].

For both Settings 1 and 2, our test set J corresponds to randomly selected subsets of images from the BSDS500 validation set. The size of the test set |J| will vary depending on the experiments, hence will be specified for each case.

Architecture comparison

We first compare the proposed unfolded NNs (for both LNO and LFO learning strategies) in terms of runtime, FLOPs and number of learnable parameters (i.e., |Θ|). These values, for (K, J) = (20, 64) are summarized in Table 2, also including the metrics for DRUnet. The experiments are conducted in PyTorch, using an Nvidia Tesla V100 PCIe 16GB.

From the table, it is obvious that the unfolded NNs have much lighter architectures than DRUnet.

Robustness comparison

Multiple works in the literature investigated NN robustness against perturbations [START_REF] Jakubovitz | Improving dnn robustness to adversarial attacks using jacobian regularization[END_REF][START_REF] Hoffman | Robust learning with jacobian regularization[END_REF][START_REF] Pesquet | Learning maximally monotone operators for image recovery[END_REF]. Formally, given an input z and a perturbation ϵ, the error on the output can be upper bounded via the inequality

∥f Θ (z + ϵ) -f Θ (z)∥ ≤ χ∥ϵ∥. (22)
The parameter χ > 0 can then be used as a certificate of the robustness for the network. This analysis is important and complementary to quality recovery performance to choose a reliable model.

According to [START_REF] Combettes | Lipschitz certificates for layered network structures driven by averaged activation operators[END_REF], in the context of feedforward NNs as those proposed in this work and described in Section 3, χ can be upper bounded by looking at the norms of each linear operator, i.e.,

χ ≤ K k=1 ∥W k,P ∥ S × ∥W k,D ∥ S . (23)
Unfortunately, as shown in [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF], this bound can be very pessimistic and not ideal to conclude on the robustness of the network. Instead, a tighter bound can be computed using the definition of Lipschitz continuity. Indeed, by definition the parameter χ in ( 22) is a Lipschitz constant of f Θ . This Lipschitz constant can be computed by noticing that it corresponds to the maximum of ∥ J f Θ (z)∥ S over all possible images z ∈ R N , where J denotes the Jacobian operator. Since such a value is impossible to compute for all images of R N , we can restrict our study on images similar to those used for training the network, i.e., χ ≈ max

(zs) s∈I ∥ J f Θ (z s )∥ S . (24)
Such an approach has been proposed in [START_REF] Pesquet | Learning maximally monotone operators for image recovery[END_REF] for constraining the value of χ during the training process. In practice, the norm is computed using power iterations coupled with automatic differentiation in Pytorch.

Motivated by these facts, we evaluate the robustness of our models by computing an approximation of χ, as described in [START_REF] Jiu | A deep primal-dual proximal network for image restoration[END_REF], considering images in the test set J instead of the training set I. Here, J corresponds to 100 images randomly selected from BSD500 validation set, and for every s ∈ J, z s = x s + w s , where w s ∼ N (0, δ s Id).

Training Setting 1 -Fixed noise level: For this setting we fix δ s ≡ 0.08. Corresponding values of χ for DD(i)FB and D(Sc)CP, with both LNO and LFO learning strategies, are reported in Figure 2. For this setting, we show the evolution of the value of χ for K ∈ {5, 10, 15, 20, 25} and J ∈ {8, 16, 32, 64}. We observe that the value of χ for LFO schemes are higher than their LNO counterparts, i.e., LFO schemes are less robust than LNO according to [START_REF] Ongie | Deep learning techniques for inverse problems in imaging[END_REF]. In addition, DDFB-LNO and D(Sc)CP-LNO seems to be the most robust schemes. Their χ value decreases slightly when K increases, and increases slighly when J increases.

Training Setting 2 -Variable noise level: For this setting, Figure 3 gives the box plots showing the distribution of the 100 values of (∥ J f Θ (z s )∥ S ) s∈J , with δ s ∼ U([0, 0.01]). For the worst median, Q1 and Q3 values.

Denoising performance comparison

Training Setting 1 -Fixed noise level: We evaluate the denoising performance of the four proposed architectures considering either the LNO or LFO learning strategy, varying K and J, on |J| = 100 noisy images obtained from the test set, with noise standard deviation δ = 0.08. The average PSNR value for these test noisy images is 21.88 dB.

In Figure 4, we show the averaged PSNR values obtained with the proposed unfolded networks. We observe a strong improvement of the denoising performances of all the NNs when increasing the size J of convolution filters, as well as a moderate improvement when increasing the depth K of the NNs. All methods have very similar performances, with the exception of ScCP (both LNO and LFO) which has much higher denoising power.

Training Setting 2 -Variable noise level: We evaluate the denoising performances of the proposed unfolded NNs on images degraded by a Gaussian noise with standard deviation δ = 0.05. Figure 5 gives the PSNR values obtained from the proposed unfolded NNs, when applied to a random subset of |J| = 20 images of BSDS500 validation set. Furthermore, Table 3 presents the average PSNR values computed for a separate subset of |J| = 100 images from the BDSD500 validation set.

In this table, we further give the average PSNR value for DRUnet. These results show that for all unfolding strategy, the LNO learning strategy improves the denoising performance over LFO. We further observe that D(Sc)CP outperforms DD(i)FB. DRUnet however outperforms the unfolded NNs on this experiments for PSNR values. For completeness, we give examples of a denoised image in Figure 6 obtained with DRUnet, DDFB and DScCP. On visual inspection, results from DDFB may appear still slightly noisy compared with DRUnet and DScCP. However results from DRUnet and DScCP are very comparable, and DScCP might reproduce slightly better some textures (e.g., the water around the castle is slightly over-smoothed with DRUnet).

Denoising performance versus robustness

According to previous observations on our experiments, in the remainder we will focus on DDFB-LNO and DScCP-LNO, both appearing to have the best compromise in terms of denoising performance and robustness.

To further assess the denoising performances of the proposed DDFB-LNO and DScCP-LNO, we evaluate their performance on denoising tasks that are different from those used during the training process. For the sake of completeness, we also compare our PNNs with DRUnet. All networks have been trained as Gaussian denoisers (Training Setting 2). We observed in the previous section that DRUnet, with 1000 times more parameters than our PNNs, shows impressive performances in terms of Gaussian denoising (in PSNR). However, we also observed that DRUnet has a higher Lipschitz constant than our PNNs. Hence we will further assess the robustness of the networks by uqing them for denoising Laplace-Gauss and Poisson-Gauss noises (i.e., on different noise settings than the training Gaussian setting). For both noise settings, the Gaussian noise level is fixed to δ = 0.05. In addition, for the Laplace noise we fix the scale to 0.05, and for the Poisson noise we fix the noise level to 50/255. The evaluation includes both visual and quantitative analysis, presented in Figure 7 and Table 4, respectively. Remarkably, the proposed PNNs lead to higher denoising performances in these scenarios, validating their higher robustness compared to DRUnet.

Image deblurring comparison with PnP

Deblurring problem Another measure to assess the proposed unfolded NNs is to use them in a Plug-and-Play framework, for image deblurring. In this context, the objective is to find an estimate x † ∈ R N of an original unknown image x ∈ R N , from degraded measurements y ∈ R M obtained through

y = Ax + b, (25)
where A : R N → R M is a linear blurring operator, and b ∈ R M models an additive white Gaussian noise, with standard deviation σ > 0. A common method to solve this inverse problem is then to find the MAP estimate x † of x, defined as a minimizer of a penalized least-squares objective function. A general formulation is given by

find x † MAP ∈ Argmin x∈C 1 2 ∥Ax -y∥ 2 + λg(Dx), (26)
where, for every t ∈ N, f K zt,λγ,Θ is either D(i)FB or D(Sc)CP. In algorithm [START_REF] Jiu | Alternative design of deeppdnet in the context of image restoration[END_REF], parameters (λ, γ) are given as inputs of f K zt,λγ,Θ . Precisely, the regularization parameter ν for the denoising problem ( 2) is chosen to be the product between the regularization parameter λ for the restoration problem [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF] and the stepsize of the algorithm γ, i.e., ν = λγ.

The following result is a direct consequence of Corollary 1 combined with convergence results of the FB algorithm [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF].

(∀t ∈ N * ) ∥x t+1 -x t ∥ ≤ ∥x t -x t-1 ∥. (28)
A few comments can be made on Theorem 4.1 and on the PnP algorithm [START_REF] Jiu | Alternative design of deeppdnet in the context of image restoration[END_REF]. First, [32, Prop. 2] is a particular case of Theorem 4.1 for DDFB-LNO. Second, in practice, only a fixed number of layers K are used in the PNNs, although Theorem 28 holds for K → +∞. In [START_REF] Repetti | Dual Forward-Backward unfolded network for flexible Plug-and-Play[END_REF] the authors studied the behavior of ( 27), using DDFB-LNO. They empirically emphasized that using warm restart for the network (i.e., using both primal and dual outputs from the network of the previous iteration) could add robustness to the PnP-FB algorithm, even when the number of layers K is fixed, due to the monotonic behaviour of the FB algorithm on the dual variable. Finally, the regularization parameter λ in [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF] aims to balance the data-fidelity term and the NN denoising power [START_REF] Repetti | Dual Forward-Backward unfolded network for flexible Plug-and-Play[END_REF]. The proposed PNNs take as an input a parameter ν > 0 that has similar interpretation for the denoising problem (3). In the PnP algorithm, we have ν = λγ, hence δ 2 = β 2 σ 2 γ, where δ 2 is the training noise level, and σ 2 is the noise level of the deblurring problem [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF]. The β 2 allows for flexibility in the choice of λ, to possibly improve the reconstruction quality.

In the remainder of the section, we will focus on unfolded NNs trained using Training setting 2 described in Section 4.1 (i.e., with variable noise level) to better fit the noise level of the inverse problem [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF].

Robustness comparison

In the context of PnP methods, robustness can be measured in terms of convergence of the global algorithm. In this context, it is known that the PnP-FB algorithm convergences if the NN is firmly non-expansive (see, e.g., [START_REF] Pesquet | Learning maximally monotone operators for image recovery[END_REF] for details). According to [43, Prop. 2.1], an operator f Θ is firmly nonexpansive if and only if h

Θ = 2f Θ -Id is a 1-Lipschitz operator, i.e., χ h = max z ∥ J h Θ (z)∥ S < 1.
In the same paper, the authors used this result to develop a training strategy to obtain firmly non-expansive NNs.

Here we propose to use this result as a measure of robustness of the proposed PNNs. Similarly to Section 4.3, we approximate χ h by computing χ h ≈ max s∈J ∥ J h Θ (z s )∥ S , where J contains 100 images randomly selected from BSD500 validation set, and z s are noisy images with standard deviation uniformly distributed in [0, 0.01]. Then, the closer this value is to 1, and the closer the associated NN f Θ is to be firmly non-expansive. Figure 8 gives the box plots showing the distribution of (∥ J h Θ (z s )∥ S ) s∈J . Conclusions on these results are very similar to those in Figure 3, in particular that DDFB-LNO and DScCP-LNO have the smallest χ h values.

Parameter choices

As emphasized previously, the denoising NNs in PnP-FB algorithm [START_REF] Jiu | Alternative design of deeppdnet in the context of image restoration[END_REF] depend on two parameters: the step-size γ and the regularization parameter λ = β 2 σ 2 . In this section we investigate the impact of the step-size on the results, while the regularization parameter is fixed to β = 1.

To ensure the convergence of the PnP-FB algorithm, the stepsize must satisfy γ ∈ (0, 2/∥A∥ 2 S ). Figure 9 aims to evaluate the stability of the proposed unfolded NNs by looking at the convergence of the associated PnP-FB iterations, varying γ ∈ {1.2, 1.4, 1.6, 1.8, 1.99, 2.1} (where ∥A∥ S = 1). For this experiment, we fix λ = σ 2 (i.e., β = 1). Further, we we also consider PnP algorithms using either DRUnet or BM3D as denoisers. The plots show the convergence profiles for the deblurring of one image in terms of PSNR values and relative error norm of consecutive iterates ∥x t+1 -x t ∥/∥x 0 ∥, with respect to the iterations. In theory, (∥x t+1 -x t ∥) t should decrease monotonically. Interestingly, we can draw similar conclusions from the curves in Figure 9 as for Figures 8 and3. Looking at the relative error norms, the most robust NNs are DDFB-LNO and DScCP-LNO, both converging monotonically for any choice of γ ≤ 1.99. DCP-LNO seems to have similar convergence profile, with a slower convergence rate. None of the other PnP schemes seems to be stable for γ = 1.99. For γ ≤ 1.8, DRUnet shows an interesting convergence profile, however the error norm is not decreasing monotonically. The remaining PnP schemes do not seem to converge in iterates, as the error norms reach a plateau. In terms of PSNR values, DScCP-LNO and BM3D have the best performances, followed by DDFB-LNO and DRUnet. For these four schemes, γ = 1.99 leads to the best PSNR values.

Restoration performance comparison

In this section, we perform further comparisons of the different PnP schemes on 12 random images

Conclusion

In this work, we have presented a unified framework for building denoising PNNs with learned linear operators, whose architectures are based on D(i)FB and (Sc)CP algorithms. We show through simulations that the proposed PNNs have similar denoising performances to DRUnet, although being much lighter (∼ 1000 less parameters). In particular, we observed that the use of inertia in the unfolding strategies improves the denoising efficiency, i.e., DDiFB (resp. DScCP) is a better denoiser than DDFB (resp. DCP). Similarly, we observe that adopting a learning strategy close to theoretical convergence conditions leads to higher denoising performance. Beside, we show that the proposed PNNs are generally more robust than DRUnet, in particular for DDFB-LNO, DCP-LNO and DScCP-LNO. Robustness is assessed in terms of the network Lipschitz constants, by applying the networks non-Gaussian noises, and by injecting the networks in PnP algorithms for a deblurring task. In particular, DDFB-LNO and DScCP-LNO exhibit the best trade-off in terms of denoising/reconstruction performances, and robustness.

Annex

The proposed ScCP iterations [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] for minimizing (3) is equivalent to the original ScCP algorithm proposed in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]:

for k = 0, 1, . . .     x k+1 = prox µ k ( 1 2 ∥•-z∥ 2 +ι C ) z -µ k D ⊤ u k , u k+1 = prox τ k (νg) * u k + τ k D (1 + α k )x k+1 -α k x k . ( 29
)
This can be shown by noticing that, for every v ∈ R N ,

prox µ k ( 1 2 ∥•-z∥ 2 +ι C ) (v) = argmin x∈R N 1 + µ k 2µ k ∥x - v + µ k z 1 + µ k ∥ 2 + ι C (x) = P C v + µ k z 1 + µ k .

  -A first strategy to solve (2)-(3) consists in applying (i)FB to the dual formulation of problem (2):

Figure 1 :

 1 Figure 1: Top: Architecture of the proposed DAH-Unified block for the k-th layer. Linearities, biases, and activation functions are shown in blue, green and red, respectively. Bottom: Inertial step for ScCP (top) and DiFB (bottom), for the k-th layer.

  and η k,P = P C , and the dual (D) operators are given by W k,D = τ k D, V k,D = Id, b z,k,P = 0, η ν,k,D = prox τ k (νg) * .

  2, and τ k = 0.99∥D k ∥ -2 DCP-LFO D k,P , D k,D , µ learn µ = µ0 = • • • = µK , and absorb τ k in D k,D DScCP-LFO D k,P , D k,D , µ0 learn µ0, absorb τ k in D k,D , and fix α k = (1 + 2µ k ) -1/2 , and

1 . 2 .

 12 Learned Normalized Operators (LNO): This strategy uses theoretical conditions ensuring convergence of D(i)FB and (Sc)CP, i.e., choosing the stepsizes appearing in deep-D(i)FB and deep-(Sc)CP according to the conditions given in Corollary 1. In this context, for every k ∈ {1, . . . , K}, we choose D k,P to be equal to the adjoint D ⊤ k,D of D k,D . However, unlike in Corollary 1, we allow D k,D to vary for the different layers k ∈ {1, . . . , K}. Learned Flexible Operator (LFO): This strategy, introduced in [26], consists in learning the stepsizes appearing in deep-D(i)FB and deep-(Sc)CP without constraints, as well as allowing a mismatch in learning the adjoint operator D ⊤ of D, i.e., learning D k,D and D k,P independently.

4. 1

 1 Denoising training settingTraining dataset -We consider two sets of images: the training set (x s , z s ) s∈I of size |I| and the test set (x s , z s ) s∈J of size |J|. For both sets, each couple (x s , z s ) consists of a clean multichannel image x s of size N s = C × N s (where C denotes the number of channels, and N s the number of pixels in each channel), and a noisy version of this image given by z s = x s + ε s with ε s ∼ N (0, δ 2 Id) for δ > 0.

Figure 2 :

 2 Figure 2: Training Setting 1: Robustness. Values χ = max s∈J ∥ J f Θ (z s )∥ S (log 2 scale) for the proposed PNNs, with J ∈ {8, 16, 32, 64} and K ∈ {5, 10, 15, 20, 25}. Top row: LNO settings. Bottom row: LFO settings.

Figure 3 :Figure 4 :

 34 Figure 3: Training Setting 2: Robustness. Distribution of (∥ J f Θ (z s )∥ S ) s∈J for 100 images extracted from BSDS500 validation dataset J, for the proposed PNNs and DRUnet.

Figure 5 :

 5 Figure 5: Training Setting 2: Denoising performance. PSNR values obtained with the proposed PNNs (with (K, J) = (20, 64)), for 20 images of BSDS500 validation set, degraded with noise level δ = 0.05.

Figure 6 :

 6 Figure 6: Training Setting 2: Denoising performance on Gaussian noise. Example of denoised images (and PSNR values) for Gaussian noise δ = 0.05 obtained with DRUnet and the proposed DDFB-LNO and DScCP-LNO, with (K, J) = (20, 64).

Theorem 4 . 1 .

 41 Let (x t ) t∈N be a sequence generated by[START_REF] Jiu | Alternative design of deeppdnet in the context of image restoration[END_REF], with f K zt,λγ,Θ being DD(i)FB-LNO or D(Sc)CP-LNO. Assume that γ ∈ (0, 2/∥A∥ 2 S ) and that, for every k ∈ {1, . . . , K}, D k,D = D and D k,P = D ⊤ for D : R N → R |F| . Under the same conditions as Corollary 1, if K → ∞, then (x t ) t∈N converges to a solution x † to problem[START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? the study case of image denoising[END_REF], and

Figure 8 :

 8 Figure 8: Deblurring (Training Setting 2): Robustness. Distribution of ∥ J h Θ (z s )∥ S , where h Θ = 2f Θ -Id for 100 images extracted from BSDS500 validation dataset J, for the proposed PNNs and DRUnet.

  = 0.015) -20.11 dB BM3D -27.10 dB DRUnet -25.09 dB DDFB-LNO -27.23 dB DScCP-LNO -26.48 dB

Figure 11 :

 11 Figure 11: Deblurring (Training Setting 2): Restoration performance. Restoration example for σ = 0.015, with parameters γ = 1.99 and β chosen optimally for each scheme.

Figure 12 :

 12 Figure 12: Deblurring (Training Setting 2): Restoration performance. Restoration example for σ = 0.03, with parameters γ = 1.99 and β chosen optimally for each scheme.

Table 2 :

 2 Architecture comparison. Runtime (in sec.), number of parameters |Θ| and FLOPs (in G) of the denoisers when used on 100 images of size 3 × 481 × 321. Values for the PNNs are given for fixed (K, J) = (20, 64).

			Time (msec)	|Θ|	FLOPs (×10 3 G)
	BM3D	13×10 3 ± 317	-	-
	DRUnet	96 ± 21	32, 640, 960	137.24
		DDFB	3 ± 1.5	34, 560
	LNO	DDiFB DDCP	3 ± 0.5 6 ± 1	34, 560 34, 561	2.26
		DDScCP	7 ± 1	34, 580
		DDFB	4 ± 17	69, 120
	LFO	DDiFB DDCP	5 ± 15 7 ± 14	69, 121 69, 121	2.26
		DDScCP	9 ± 15	69, 160

• Training Setting 2 -Variable noise level: The NNs are trained using the test dataset of ImageNet [39] (|I| = 5 × 10 4 ), with learning rate for ADAM set to 1 × 10 -3 (all other parameters set as default), batch size of 200, and patches of size 50 × 50 randomly selected. Further, we consider a variable noise level, i.e., images are degraded by a Gaussian noise with standard deviation δ i , for i ∈ I, selected randomly with uniform distribution in [0, 0.1].

Table 3 :

 3 Training Setting 2: Denoising performance. Average PSNR (and standard deviation) values (in dB) obtained with the proposed PNNs with (K, J) = (20, 64) and with DRUnet, for 100 noisy images of BSDS500 validation set (δ = 0.05, input PSNR= 25.94dB).

	DRUnet	DDFB	DDiFB	DCP	DScCP
			LNO
	34.7 ±1.89	32.42 ±0.86 33.12 ±1.27 33.53 ±1.46 33.57 ±1.53 LFO
		32.09 ±0.91 32.84 ±1.12 33.32 ±1.35 33.27 ±1.32

Table 4 :

 4 Training Setting 2: Denoising performance vs Robustness. Average PSNR (and standard deviation) values (in dB) obtained with DRUnet and the proposed DDFB-LNO and DScCP-LNO (with (K, J) = (20, 64)), for 12 noisy images of BSDS500 validation set degraded with Poisson-Gauss and Laplace-Gauss noises.

		Noisy	DRUnet	DDFB-LNO	DScCP-LNO
	Poisson-Gauss	19.34 ± 1.68	21.59 ±1.69	23.83 ±1.80	25.44 ±1.25
	Laplace-Gauss	20.72 ± 1.09	24.66 ±0.52	26.26 ±0.51	28.42 ±1.07

This work was partly supported by the ANR (Agence Nationale de la Recherche) from France ANR-19-CE48-0009 Multisc'In, Fondation Simone et Cino Del Duca (Institut de France), the CBP (Centre Blaise Pascal-ENS lyon), the Royal Society of Edinburgh, and the EPSRC grant EP/X028860.

where C ⊂ R N , D : R N → R |F| and g : R |F| → (-∞, +∞] are defined as in (3), and λ ∝ σ 2 (i.e., there exists β > 0 such that λ = (βσ) 2 ) is a regularization parameter.

PnP-FB algorithm

The idea of PnP algorithms is to replace the penalization term (often handled by a proximity operator) by a powerful denoiser. There are multiple choices of denoisers, that can be classified into two main categories: hand-crafted denoisers (e.g. BM3D [START_REF] Dabov | Image denoising by sparse 3D transform-domain collaborative filtering[END_REF]) and learning-based denoisers (e.g., DnCNN [START_REF] Zhang | Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising[END_REF] and UNet [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]). PnP methods with NNs have recently been extensively studied in the literature, and widely used for image restoration (see, e.g., [START_REF] Repetti | Dual Forward-Backward unfolded network for flexible Plug-and-Play[END_REF][START_REF] Pesquet | Learning maximally monotone operators for image recovery[END_REF][START_REF] Kamilov | Plug-and-play methods for integrating physical and learned models in computational imaging: Theory, algorithms, and applications[END_REF][START_REF] Hurault | Gradient step denoiser for convergent plug-and-play[END_REF]).

In this section, as proposed in [START_REF] Repetti | Dual Forward-Backward unfolded network for flexible Plug-and-Play[END_REF], we plug the proposed unfolded NNs in a FB algorithm to solve [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF]. The objective is to further assess the robustness of the proposed unfolded strategies. Following the approach proposed in [START_REF] Repetti | Dual Forward-Backward unfolded network for flexible Plug-and-Play[END_REF], the PnP-FB algorithm is given by selected from BSDS500 validation set, degraded as per model [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF]. In particular, we will run PSNR with optimal β experiments for three different noise levels σ ∈ {0.015, 0.03, 0.05}. Since in the previous sections we observed that DDFB-LNO and DScCP-LNO are the more robust unfolded strategies, in this section we only focus on these two schemes, comparing them to BM3D and DRUnet. For each denoiser, we choose δ = λγ, with γ = 1.99 and λ = β 2 σ 2 . As explained after Theorem 4.1, parameter β 2 allows for flexibility to possibly improve the reconstruction quality. In the remainder we choose it to optimize the reconstruction quality of each image (PSNR value).

In Figure 10, we provide best PSNR values obtained (for optimized values of β) with DDFB-LNO, DScCP-LNO, DRUnet and BM3D, on 12 images from BSDS500 validation set degraded according to [START_REF] Aggarwal | MoDL: Model-based deep learning architecture for inverse problems[END_REF], with σ = 0.03. In addition, the averaged PSNR values obtained with the four different schemes for the different noise levels σ ∈ {0.015, 0.03, 0.05} are given in Table 5. It can be observed that, regardless the noise level σ, DDFB-LNO and BM3D always have the highest PSNR values, outperforming DRUnet and DScCP-LNO. However, DDFB-LNO has a much cheaper computation time than BM3D (see Table 2 for details). For low noise level σ = 0.015, DScCP-LNO also outperforms DRUnet. For highest noise levels σ ∈ {0.03, 0.05}, DScCP-LNO and DRUnet have similar performances.

For visual inspection, we also provide in Figures 11 and12 examples of two different images, for noise levels σ = 0.015, σ = 0.03 and σ = 0.05, respectively. We observe that DRUnet and DScCP-