Caroline Chaux
email: caroline.chaux@univ-paris-est.fr.

Jean-Christophe Pesquet
email: jean-christophe.pesquet@univ-paris-est.fr.

Nelly Pustelnik
email: nelly.pustelnik@univ-paris-est.fr.

† C Chaux

J.-C Pesquet

Nested iterative algorithms for convex constrained image recovery problems *

The objective of this paper is to develop methods for solving image recovery problems subject to constraints on the solution. More precisely, we will be interested in problems which can be formulated as the minimization over a closed convex constraint set of the sum of two convex functions f and g, where f may be non-smooth and g is differentiable with a Lipschitz-continuous gradient. To reach this goal, we derive two types of algorithms that combine forward-backward and Douglas-Rachford iterations. The weak convergence of the proposed algorithms is proved. In the case when the Lipschitz-continuity property of the gradient of g is not satisfied, we also show that, under some assumptions, it remains possible to apply these methods to the considered optimization problem by making use of a quadratic extension technique. The effectiveness of the algorithms is demonstrated for two wavelet-based image restoration problems involving a signal-dependent Gaussian noise and a Poisson noise, respectively.

Introduction

Wavelet decompositions [START_REF] Mallat | A wavelet tour of signal processing[END_REF] proved their efficiency in solving many inverse problems. More recently, frame representations such as Bandlets [START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF], Curvelets [START_REF] Candès | Recovering edges in ill-posed inverse problems: Optimality of curvelet frames[END_REF], Grouplets [START_REF]Geometrical grouplets[END_REF] or dual-trees [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF][START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF] have gained much popularity. These linear tools provide geometrical representations of images and they are able to easily incorporate a priori information (e.g. via some simple statistical models) on the data. Variational or Bayesian formulations of inverse problems using such representations often lead to the minimization of convex objective functions including a non-differentiable term having a sparsity promoting role [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF][START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF][START_REF] Antoniadis | Wavelet thresholding for some classes of non-Gaussian noise[END_REF][START_REF] Candès | Sparsity and incoherence in compressive sampling[END_REF][START_REF] Tropp | Just Convex programming methods for identifying sparse signals in noise[END_REF][START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF].

In restoration problems, the observed data are corrupted by a linear operator and a noise which is not necessarily additive. To solve this problem, one can adopt a variational approach, aiming at minimizing the sum of two functions f and g over a convex set C in the transform domain. Throughout the paper, f and g are assumed to be in the class Γ 0 (H) of lower semicontinuous convex functions taking their values in ] -∞, +∞] which are proper (i.e. not identically equal to +∞) and defined on a real separable Hilbert space H. Then, our objective is to solve the following: Problem 1.1 Let C be a nonempty closed convex subset of H. Let f and g be in Γ 0 (H), where g is differentiable on H with a β-Lipschitz continuous gradient for some β ∈ ]0, +∞[. (∀x ∈ H) ι C (x) = 0, if x ∈ C; +∞, otherwise.

Up to now, many authors devoted their works to the unconstrained case, i.e. C = H. So-called thresholded Landweber algorithms belonging to the more general class of forward-backward optimization methods were proposed in [START_REF] Figueiredo | An EM algorithm for wavelet-based image restoration[END_REF][START_REF] Bect | A l 1 -unified variational framework for image restoration[END_REF][START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Bredies | A generalized conditional gradient method and its connection to an iterative shrinkage method[END_REF] in order to solve the problem numerically. Daubechies et al. [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] investigated the convergence of these algorithms in the particular case when g is a quadratic function and f is a weighted ℓ p -norm with p ∈ [START_REF] Alber | On the projected subgradient method for nonsmooth convex optimization in a hilbert space[END_REF][START_REF] Anscombe | The transformation of Poisson, binomial and data[END_REF]. These approaches were put into a more general convex analysis framework in [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] and extended to frame representations in [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]. Attention was also paid to the improvement of the convergence speed of the forward-backward algorithm in [START_REF] Bioucas-Dias | A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF], for some specific choices of f and g. In [START_REF] Vonesch | A fast thresholded Landweber algorithm for wavelet-regularized multidimensional deconvolution[END_REF], an accelerated method was suggested in the specific case of a deconvolution in a Shannon wavelet basis. Then, a Douglas-Rachford algorithm relaxing the assumption of differentiability of g was introduced in [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF]. In recent works [START_REF] Dupé | Deconvolution of confocal microscopy images using proximal iteration and sparse representations[END_REF][START_REF]A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF], a variational approach, which is grounded on a judicious use of the Anscombe transform, was developed for the deconvolution of data contaminated by Poisson noise. A modification of the forward-backward algorithm was subsequently proposed in finite dimension in order to solve the associated optimization problem. Additional comments concerning this approach will be given in Sections 3.2.2 and 5.4. A key tool in the study of the aforementioned methods is the proximity operator introduced by Moreau in 1962 [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF][START_REF]Proximité et dualité dans un espace hilbertien[END_REF]. The proximity operator of f ∈ Γ 0 (H) is prox f : H → H : x → arg min y∈H 1 2 y -x 2 + f (y). We thus see that prox ιC reduces to the projection P C onto the convex set C. The function f in Problem 1.1 may be non-smooth and, actually, it is often chosen as an ℓ 1 -norm, in which case its proximity operator reduces to a componentwise soft-thresholding [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. In [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF], the authors derived the concept of proximal thresholding by considering a larger set of non-differentiable convex functions.

The goal of this paper is to propose iterative algorithms allowing us to solve Problem 1.1 when C = H. The relevance of the proposed methods is shown for image recovery problems where convex constraints on the solution need to be satisfied.

In Section 2, we start by recalling some properties of the proximity operator. Then, in Section 3 we briefly describe the forward-backward and Douglas-Rachford methods. As the proximity operator of the sum of the indicator function of a convex set and a function in Γ 0 (H) cannot be easily expressed in general, we propose two iterative methods to compute this operator: the first one is a forward-backward algorithm, whereas the second one is a Douglas-Rachford algorithm. We also investigate the specific convergence properties of these two algorithms. In Section 4, we derive two iterative methods to solve Problem 1.1 and their convergence behaviours are studied. Finally, in Section 5, these algorithms are applied to a class of image recovery problems. In this case, the Lipschitz-continuity property of the gradient of g is not satisfied in the considered maximum a posteriori criterion. To overcome this difficulty, a quadratic extension technique providing a lower approximation of the objective function is introduced. Numerical results concerning deconvolution problems in the presence of signal-dependent Gaussian noise or Poisson noise are then provided.

Some properties of proximity operators

As already mentioned, the proximity operator of ι C + f plays a key role in our approach. Some useful results for the calculation of prox ιC+f are first recalled. Subsequently, the domain of a function f : Proposition 12] Let f ∈ Γ 0 (H) and let C be a closed convex subset of H such that C ∩ dom f = ∅. Then the following properties hold.

H →]-∞, +∞] is denoted by dom f = {x ∈ H | f (x) < +∞}. Proposition 2.1 [18,
(i) (∀x ∈ H), prox f x ∈ C ⇒ prox ιC+f x = prox f x (ii) Suppose that H = R. Then prox ιC+f = P C • prox f . (1)
Note that, the second part of this proposition can be generalized, yielding the following result which appears also as an extension of [14, Proposition 2.10] when C = H: Proposition 2.2 Let K be a nonempty subset of N, (o k ) k∈K be an orthonormal basis of H and (ϕ k ) k∈K be functions in

Γ 0 (R). Set f : H → ]-∞, +∞] : x → k∈K ϕ k ( x, o k ). ( 2
)
Let C = k∈K {x ∈ H | x, o k ∈ C k } (3)
where

(C k ) k∈K are closed intervals in R such that (∀k ∈ K) C k ∩ dom ϕ k = ∅.
Suppose that either K is finite, or there exists a subset L of K such that:

(i) K L is finite; (ii) (∀k ∈ L) ϕ k ≥ ϕ k (0) = 0 and 0 ∈ C k .
Then,

(∀x ∈ H) prox ιC+f x = k∈K π k o k (4)
where

π k =      inf C k if prox ϕ k x, o k < inf C k sup C k if prox ϕ k x, o k > sup C k prox ϕ k x, o k otherwise. (5)
Proof. Due to the form of f and C, one can write,

(∀x ∈ H) f + ι C (x) = k∈K (ϕ k + ι C k )( x, o k ). For every k ∈ K, ϕ k + ι C k ∈ Γ 0 (R) since ϕ k ∈ Γ 0 (R)
and C k is assumed to be a closed convex set having a nonempty intersection with dom ϕ k . If K is not finite, in view of Assumption (ii), we have (∀k ∈ L)

ϕ k + ι C k ≥ (ϕ k + i C k )(0) = 0.
From [14, Remark 3.2(ii) and Proposition 2.10], it can be deduced that

(∀x ∈ H) prox f +ιC x = k∈K prox ϕ k +ιC k x, o k o k . (6)
On the other hand, since for every

k ∈ K, C k is a closed interval in R such that C k ∩ dom ϕ k = ∅, it follows from Proposition 2.1(ii), that prox ϕ k +ιC k x, o k =(P C k • prox ϕ k )( x, o k ) =      inf C k , if prox ϕ k x, o k < inf C k prox ϕ k x, o k , if prox ϕ k x, o k ∈ C k sup C k , if prox ϕ k x, o k > sup C k . (7)
Combining ( 6) and ( 7) yields ( 4) and [START_REF] Bect | A l 1 -unified variational framework for image restoration[END_REF].

A function f (resp. convex C) satisfying (2) (resp. (3)) will be said separable. Note that ( 4) and (5) imply that (1) holds. However, this relation has been proved under the restrictive assumption that both f and C are separable. In general, when either f or C is not separable, (1) is no longer valid. Let us give two simple counterexamples to illustrate this fact.

Example 2.3 Let H = R 2 and f be the function defined by

(∀x ∈ R 2 ) f (x) = 1 2 x ⊤ Λx with Λ = 1 Λ 1,2 Λ 1,2 Λ 2,2
where

Λ 2,2 ≥ 0 and |Λ 1,2 | ≤ Λ 1/2 2,2 . Let C = [-1, 1] 2 . This convex set is separable w.r.t. the canonical basis of R 2 . Now, set x = 2(Λ 1,2 , 1 + Λ 2,2 ) ⊤ .
After some calculations (see Appendix A), one obtains:

• P C (prox f x) = (0, 1) ⊤ • prox ιC+f x = (π, 1) ⊤ where π =      Λ1,2 2 if Λ 1,2 ∈ [-2, 2] 1 if Λ 1,2 > 2 -1 if Λ 1,2 < -2. ( 8
)
We conclude that (1) is not satisfied as soon as Λ 1,2 = 0, that is f is not separable.

Example 2.4 Let H = R 2 . Consider the separable function defined by (∀x = (x (1) , x (2)

) ⊤ ∈ R 2 ) f (x) = (1 + Λ 1,2 )(x (1) ) 2 + (1 -Λ 1,2 )(x (2) ) 2 where 0 < |Λ 1,2 | ≤ 1. Let the nonseparable convex set C be defined by C = {x = (x (1) , x (2) ) ⊤ ∈ R 2 | max(|x (1) -x (2) |, |x (1) + x (2) |) ≤ √ 2}.
In this case, it is shown in Appendix B that (1) does not hold.

In summary, for an arbitrary function in Γ 0 (H) and an arbitrary closed convex set, we cannot trust (1) to determine the proximity operator of the sum of this function and the indicator function of the convex set. In the next section, we will propose efficient approaches to compute the desired proximity operator in a general setting.

Other more classical properties of the proximity operator which will be used in the paper are provided in the sequel.

Proposition 2.5 (i) If f = h + κ •, x where h ∈ Γ 0 (H), x ∈ H and κ ∈ R, then prox f = prox h (• -κx). (ii) If f = h + ϑ • 2 /2 where h ∈ Γ 0 (H) and ϑ ∈ ]0, +∞[, then (a) prox f = prox (1+ϑ) -1 h • /(1 + ϑ) (b) (∀(y, z) ∈ H 2 ) prox f y -prox f z, y -z ≥ (1 + ϑ) prox f y -prox f z 2 (c) prox f is strictly contractive 1 with constant (1 + ϑ) -1 .
Proof. Properties (i) and (ii)(a) result from straightforward calculations [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]Lemma 2.6]. (ii)(b) follows from the fact that prox (1+ϑ) -1 h is firmly nonexpansive [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]Lemma 2.4], i.e.

(∀(y, z) ∈ H 2 ) prox h 1+ϑ y -prox h 1+ϑ z, y -z ≥ prox h 1+ϑ y -prox h 1+ϑ z 2 .
Thus, by using (ii)(a), we have

(∀(y, z) ∈ H 2 ) prox f y -prox f z, y -z =(1 + ϑ) prox h 1+ϑ y 1 + ϑ -prox h 1+ϑ z 1 + ϑ , y 1 + ϑ - z 1 + ϑ ≥(1 + ϑ) prox h 1+ϑ y 1 + ϑ -prox h 1+ϑ z 1 + ϑ 2 =(1 + ϑ) prox f y -prox f z 2 .
Property (ii)(c) can then be deduced, by invoking the Cauchy-Schwarz inequality:

(∀(y, z) ∈ H 2 ) (1 + ϑ) prox f y -prox f z 2 ≤ prox f y -prox f z, y -z ≤ prox f y -prox f z y -z .
Recall that a function f ∈ Γ 0 (H) satisfying the assumptions in (ii) is said to be strongly convex with modulus ϑ. Proposition 2.6 [18, Proposition 11] Let G be a real Hilbert space, let f ∈ Γ 0 (G), and let L : H → G be a bounded linear operator. Suppose that the composition of L and L * satisfies L • L * = ν Id, for some ν ∈ ]0, +∞[. Then f • L ∈ Γ 0 (H) and

prox f •L = Id + ν -1 L * • (prox νf -Id) • L. ( 9
)
3 Iterative solutions to the minimization of a sum of two convex functions

Forward-backward approach

Consider the following optimization problem, which is a specialization of Problem 1.1:

Problem 3.1 Let f 1 and f 2 be two functions in Γ 0 (H) such that Argmin f 1 +f 2 = ∅ and f 2 is differentiable on H with a β-Lipschitz continuous gradient for some β ∈ ]0, +∞[.

Find min x∈H f 1 (x) + f 2 (x).
As mentioned in the introduction, the forward-backward algorithm is an effective method to solve the above problem.

Algorithm [20, Eq.(3.6)]

Let x 0 ∈ H be an initial value. The algorithm constructs a sequence (x n ) n∈N by setting, for every n ∈ N,

x n+1 = x n + λ n prox γnf1 (x n -γ n ∇f 2 (x n ) + b n ) + a n -x n ( 10
)
where γ n > 0 is the algorithm step-size, λ n > 0 is a relaxation parameter and a n ∈ H (resp. b n ∈ H) represents an error allowed in the computation of the proximity operator (resp. the gradient). The weak convergence of (x n ) n∈N to a solution to Problem 3.1 is then guaranteed provided that:

Assumption 3.1 (i) 0 < γ ≤ γ < 2β -1 where γ = inf n∈N γ n and γ = sup n∈N γ n . (ii) (∀n ∈ N) 0 < λ ≤ λ n ≤ 1.
(iii) n∈N a n < +∞ and n∈N b n < +∞.

More details concerning this algorithm can be found in [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF] and conditions for the strong convergence of the algorithm are also given in [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]. An additional result which will be useful in this paper is the following: Lemma 3.2 Suppose that Assumptions 3.1(i) and (ii) as well as the assumptions of Problem 3.1 hold. If f 1 is a strongly convex function with modulus ϑ, then the forward-backward algorithm in [START_REF] Byrne | Iterative image reconstruction algorithms based on cross-entropy minimization[END_REF] with a n ≡ b n ≡ 0 converges linearly to the unique solution x to Problem 3.1. More precisely, we have

(∀n ∈ N) x n -x ≤ 1 - λγϑ 1 + γϑ n x 0 -x . (11)
Proof. Since Argmin f 1 + f 2 = ∅ and f 1 is strongly (thus strictly) convex, there exists a unique minimizer x of f 1 + f 2 . Then, x is a fixed point of the forward-backward algorithm in [START_REF] Byrne | Iterative image reconstruction algorithms based on cross-entropy minimization[END_REF] when a n ≡ b n ≡ 0. Thus, we have, for all n ∈ N,

x n+1 -x = (1 -λ n )(x n -x) + λ n prox γnf1 (x n -γ n ∇f 2 (x n )) -prox γnf1 ( x -γ n ∇f 2 ( x))
which yields

x n+1 -x ≤ (1 -λ n ) x n -x + λ n prox γnf1 (x n -γ n ∇f 2 (x n )) -prox γnf1 ( x -γ n ∇f 2 ( x)) .
Since f 1 has been assumed strongly convex with modulus ϑ, γ n f 1 is strongly convex with modulus γ n ϑ and, according to Assumption 3.1(i), it is also strongly convex with modulus γϑ. We deduce from Proposition 2.5(ii)(c) that prox γnf1 is strictly contractive with constant (1 + γϑ) -1 . Hence, we have

x n+1 -x ≤ (1 -λ n ) x n -x + λ n 1 + γϑ x n -γ n ∇f 2 (x n ) -x + γ n ∇f 2 ( x) .
Recall that an operator R :

H → H is nonexpansive if (∀(y, z) ∈ H 2 ) R(x) -R(y) ≤ x -y . An operator T : H → H is α-averaged with α ∈]0, 1[ if T = (1 -α)Id+ αR where R is a nonexpansive operator.
Since f 2 is a differentiable convex function having a β-Lipschitz continuous gradient with β > 0, we deduce from the Baillon-Haddad theorem [START_REF] Baillon | Quelques propriétés des opérateurs angle-bornés et n-cycliquement monotones[END_REF], that ∇f 2 /β is 1/2-average. As γ n ∈]0, 2/β[ , by using [ This entails that

x n -γ n ∇f 2 (x n ) -x + γ n ∇f 2 ( x) ≤ x n -x
and, consequently,

x n+1 -x ≤ 1 - λ n γϑ 1 + γϑ x n -x ≤ 1 - λγϑ 1 + γϑ x n -x
which results in [START_REF] Candès | Recovering edges in ill-posed inverse problems: Optimality of curvelet frames[END_REF].

The linear convergence of the forward-backward algorithm was also proved in [START_REF] Bredies | Linear convergence of iterative soft-thresholding[END_REF][START_REF] Chen | Convergence rates in forward-backward splitting[END_REF] under different assumptions.

Computation of prox ιC+κg

Let κ > 0 and g be a differentiable function with β-Lipschitz continuous gradient where β ∈ ]0, +∞[. Let C be a closed convex set such that C = ∅. Then, for every x ∈ H, the determination of prox ιC+κg x can be viewed as a minimization problem of the form of Problem 3.1. Indeed, by using the definition of the proximity operator, we have:

(∀x ∈ H) prox κg+ιC x = arg min y∈H 1 2 y -x 2 + κg(y) + ι C (y). Now, we can set f 1 = 1 2 . -x 2 + ι C and f 2 = κg. The proximity operator of γ n f 1 with γ n ∈ ]0, +∞[, is the proximity operator of γn 2 • 2 -γ n •, x + ι C
, which is straightforwardly deduced from Proposition 2.5(i) and (ii)(a):

(∀y ∈ H) prox γnf1 y = P C y + γ n x 1 + γ n . ( 12
)
whereas f 2 has a κβ-Lipschitz continuous gradient. In this case, by setting a n ≡ b n ≡ 0 in Algorithm [START_REF] Byrne | Iterative image reconstruction algorithms based on cross-entropy minimization[END_REF], we get (∀n ∈ N)

x n+1 = x n + λ n P C x n -γ n (κ∇g(x n ) -x) 1 + γ n -x n (13) with 0 < γ ≤ γ n ≤ γ < 2κ -1 β -1 . (14)
The obtained algorithm possesses the following properties:

Proposition 3.3 Suppose that Condition (14) and Assumption 3.1(ii) hold. Consider the algorithm in [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF] where x ∈ H. Then, (i) we have:

(∀n ∈ N) x n -prox ιC+κg x ≤ ρ n x 0 -prox ιC+κg x (15)
where

ρ = 1 - λγ 1 + γ ; ( 16
)
(ii) by setting x 0 = prox κg x, we get:

prox κg x ∈ C ⇒ (∀n ∈ N) x n = prox ιC+κg x. (17)
Proof. (i) : As f 1 is strongly convex with modulus 1, ( 15) is obtained by invoking Lemma 3.2.

(ii) : If x 0 = prox κg x ∈ C, then (15) leads to

(∀n ∈ N) x n -prox ιC+κg x ≤ 1 - λγ 1 + γ n prox κg x -prox ιC+κg x = 0 ( 18
)
where Proposition 2.1(i) has been used in the last equality. This shows that ( 17) is satisfied.

Remark 3.4

(i) Eq. [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF] shows that (x n ) n∈N converges linearly to prox ιC+κg x. Although this equation provides an upper bound, it suggests to choose λ n and γ n as large as possible (i.e. λ n ≡ 1 and γ n close to 2κ -1 β -1 ) to optimize the convergence rate. This fact was confirmed by our simulations.

(ii) Proposition 3.3(ii) may appear as a desirable property since Proposition 2.1(i) states that, when prox κg x ∈ C, prox ιC+κg x takes a trivial form. In this case, the convergence is indeed guaranteed in just one iteration by appropriately initializing the algorithm. Note however that prox κg x may not always be simple to compute, depending on the form of g.

(iii) An alternative numerical method for the computation of prox ιC+κg x would consist of setting

f 1 = ι C and f 2 = 1 2 . -x 2
+ κg in the forward-backward algorithm, so yielding

(∀n ∈ N) x n+1 = x n + λ n P C (x n -γ n (κ∇g(x n ) + x n -x)) -x n with 0 < γ ≤ γ < 2(κβ + 1) -1 .
It can be noticed that the forward-backward algorithm then reduces to a projected gradient algorithm [6, Chap. 3., Sect.

3.3.2][1], when λ n ≡ 1.
In our experiments, it was however observed that the convergence of this algorithm is slower than that in [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF], probably due to the fact that prox γnf1 is no longer strictly contractive for the second choice of f 1 .

Douglas-Rachford approach

Let us relax the Lipschitz continuity assumption in Problem 3.1 and turn our attention to the optimization problem:

Problem 3.2 Let g 1 and g 2 be functions in Γ 0 (H) such that Argmin g 1 + g 2 = ∅. Assume that one of the following three conditions is satisfied:

(i) dom g 2 ∩ int dom g 1 = ∅. 2 (ii) dom g 1 ∩ int dom g 2 = ∅. (iii) H is finite dimensional and rint dom g 1 ∩ rint dom g 2 = ∅. Find min z∈H g 1 (z) + g 2 (z).
In the statement of the above problem, the notation differs from that used in Problem 3.1 to emphasize the difference in the assumptions which have been adopted and facilitate the presentation of the algorithms subsequently presented in Section 4.

The Douglas-Rachford algorithm, proposed in [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF][START_REF] Eckstein | On the Douglas-Rachford splitting methods and the proximal point algorithm for maximal monotone operators[END_REF], provides an appealing numerical solution to Problem 3.2, as described next.

Algorithm [18, Eq.(19)]

Set z 0 ∈ H and compute, for every m ∈ N,

z m+ 1 2 = prox κg2 z m + b m z m+1 = z m + τ m prox κg1 (2z m+ 1 2 -z m ) + a m -z m+ 1 2 (19)
where κ > 0, (τ m ) m∈N is a sequence of positive reals, and (a m ) m∈N (resp. (b m ) m∈N ) is a sequence of errors in H allowed in the computation of the proximity operator of κg 1 (resp. κg 2 ). Then, (z m ) m∈N converges weakly to z ∈ H [17, Corollary 5.2] such that prox κg2 z is a solution to Problem 3.2, provided that:

Assumption 3.5 (i) (∀m ∈ N) τ m ∈]0, 2[ and m∈N τ m (2 -τ m ) = +∞. (ii) m∈N τ m ( a m + b m ) < +∞.
An alternate convergence result is the following:

Proposition 3.6 Suppose that the assumptions of Problem 3.2 hold. If g 2 is a strongly convex function, then the Douglas-Rachford algorithm in [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]

with inf m∈N τ m > 0, sup m∈N τ m ≤ 2 and a m ≡ b m ≡ 0 is such that (z m+1/2
) m∈N converges strongly to the unique solution to Problem 3.2.

Proof. Let the rprox operator be defined, for every f ∈ Γ 0 (H), by

rprox f = 2prox f -Id. ( 20
)
Let us rewrite the Douglas-Rachford iteration in [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF] with a m ≡ b m ≡ 0 as z m+1 = S m z m , where

S m = τ m prox κg1 (rprox κg2 ) + Id -τ m prox κg2 . (21)
For all (y, y ′ ) ∈ H 2 , we have Lemma 2.4] and, the expression in ( 22) can be upper bounded as follows

S m y -S m y ′ 2 = τ 2 m prox κg1 (rprox κg2 y) -prox κg1 (rprox κg2 y ′ ) 2 + 2τ m prox κg1 (rprox κg2 y) -prox κg1 (rprox κg2 y ′ ), y -τ m prox κg2 y -y ′ + τ m prox κg2 y ′ + y -τ m prox κg2 y -y ′ + τ m prox κg2 y ′ 2 . (22) Since κg 1 ∈ Γ 0 (H), prox κg1 is firmly nonexpansive [20,
S m y -S m y ′ 2 ≤ τ 2 m prox κg1 (rprox κg2 y) -prox κg1 (rprox κg2 y ′ ), rprox κg2 y -rprox κg2 y ′ +2τ m prox κg1 (rprox κg2 y) -prox κg1 (rprox κg2 y ′ ), y -τ m prox κg2 y -y ′ + τ m prox κg2 y ′ + y -τ m prox κg2 y -y ′ + τ m prox κg2 y ′ 2
which yields after simplifications:

S m y -S m y ′ 2 ≤ τ m (2 -τ m ) prox κg1 (rprox κg2 y) -prox κg1 (rprox κg2 y ′ ), y -y ′ + y -τ m prox κg2 y -y ′ + τ m prox κg2 y ′ 2 .
Using the definition of the operator S m in (21), we thus obtain, after some simple calculations,

S m y -S m y ′ 2 ≤ (2 -τ m ) S m y -S m y ′ , y -y ′ + (τ m -1) y -y ′ 2 -τ 2 m prox κg2 y -prox κg2 y ′ , y -y ′ -prox κg2 y -prox κg2 y ′ 2 . (23)
Let θ be the modulus of the strongly convex function g 2 . Then κg 2 is strongly convex with modulus κθ and Proposition 2.5(ii)(b) states that the following inequality holds:

prox κg2 y -prox κg2 y ′ , y -y ′ ≥ (κθ + 1) prox κg2 y -prox κg2 y ′ 2 ,
which combined with [START_REF] Dupé | Deconvolution of confocal microscopy images using proximal iteration and sparse representations[END_REF] leads to

S m y -S m y ′ 2 +κθτ 2 m prox κg2 y -prox κg2 y ′ 2 ≤ (2 -τ m ) S m y -S m y ′ , y -y ′ + (τ m -1) y -y ′ 2 . ( 24
)
Now, let z be the unique minimizer of g 1 + g 2 . Hence, z = prox κg2 z where z is a fixed point of S m . Consequently, by setting y = z m and y ′ = z in ( 24), we deduce that

z m+1 -z 2 + κθτ 2 m z m+ 1 2 -z 2 ≤ (2 -τ m ) z m+1 -z, z m -z + (τ m -1) z m -z 2 . ( 25
)
Using the fact that

2 z m+1 -z, z m -z = z m+1 -z 2 + z m -z 2 -z m+1 -z m 2 (25)
can be rewritten as

τ m z m+1 -z 2 + (2 -τ m ) z m+1 -z m 2 + 2κθτ 2 m z m+ 1 2 -z 2 ≤ τ m z m -z 2 . (26)
Considering Assumption 3.5, (2 -τ m ) z m+1 -z m 2 is nonnegative and the left-hand side term of inequality (26) can be lower bounded, so yielding

τ m z m+1 -z 2 + 2κθτ 2 m z m+ 1 2 -z 2 ≤ τ m z m -z 2 .
Finally, by using the assumption that τ = inf m∈N τ m > 0, we obtain

z m+1 -z 2 + 2κθτ z m+ 1 2 -z 2 ≤ z m -z 2 . ( 27
)
This entails that z m+1 -z 2 ≤ z m -z 2 and, the sequence ( z m -z ) m∈N being decreasing, there exists c ∈ ]0, +∞[ such that lim m→+∞ z m -z = c. In turn, from [START_REF] Fessler | Hybrid poisson/polynomial objective functions for tomographic image reconstruction from transmission scans[END_REF], we conclude that lim m→+∞ z m+ 1 2 = z, which shows the strong convergence of (z m+1/2 ) m∈N to the unique minimizer of g 1 + g 2 .

It can be noticed that, although the convergence of the Douglas-Rachford algorithm generally requires that τ m < 2, the strong convergence is obtained under the above assumptions, when τ m = 2. The limit case of the Douglas-Rachford corresponding to τ m ≡ 2 is known as the Peaceman-Rachford algorithm [START_REF] Peaceman | The numerical solution of parabolic and elliptic differential equations[END_REF][START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF].

Computation of prox ιC+γf

Let C be a nonempty closed convex set of H. The Douglas-Rachford algorithm can be used to compute prox ιC+γf where f ∈ Γ 0 (H) and γ is a positive constant, using again the definition of the proximity operator:

(∀x ∈ H) prox ιC+γf x = arg min y∈H 1 2 y -x 2 + ι C (y) + γf (y). ( 28
)
The above minimization problem appears as a specialization of Problem 3.2 by setting g 1 = γf and g 2 = 1 2 • -x 2 + ι C , provided that one of the following three conditions holds:

Assumption 3.7 (i) C ∩ int dom f = ∅. (ii) dom f ∩ int C = ∅. (iii) H is finite dimensional and rint C ∩ rint dom f = ∅.
Subsequently, we propose to use the Douglas-Rachford algorithm in [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF] with a m ≡ b m ≡ 0, to compute the desired proximity operator. Note that both prox κg1 and prox κg2 with κ > 0, have to be calculated to apply this algorithm. In our case, we have prox κg1 = prox κγf and, similarly to [START_REF] Candès | Sparsity and incoherence in compressive sampling[END_REF],

(∀z ∈ H) prox κg2 z = P C z + κx 1 + κ .
The resulting Douglas-Rachford iterations read: for every m ∈ N,

   z m+ 1 2 = P C z m + κx 1 + κ z m+1 = z m + τ m prox κγf (2z m+ 1 2 -z m ) -z m+ 1 2 . ( 29
)
This algorithm enjoys the following properties:

Proposition 3.8 Suppose that one of Assumptions 3.7(i), 3.7(ii) or 3.7(iii) holds. Consider the algorithm in (29) where x ∈ H, inf m∈N τ m > 0 and sup m∈N τ m ≤ 2. Then,

(i) (z m+ 1
2 ) m∈N converges strongly to prox ιC+γf x; (ii) by setting κ = 1 and z 0 = 2prox γf x -x, we get:

prox γf x ∈ C ⇒ (∀m ∈ N) z m+ 1 2 = prox ιC+γf x. (30)
Proof. (i): As g 2 is strongly convex with modulus 1, (i) holds by invoking Proposition 3.6.

(ii): Set κ = 1, z 0 = 2prox γf x -x with prox γf x ∈ C. By considering the first iteration of the Douglas-Rachford algorithm (m = 0), we have z 1 2 = prox γf x and z 1 = z 0 . So, by induction, (∀m ∈ N) z m+ 1 2 = prox γf x, which is also equal to prox ιC+γf x according to Proposition 2.1(i). Remark 3.9 (i) As already observed in Remark 3.4(ii), (30) is a desirable property. It shows that the proposed algorithm converges in one iteration when prox γf x ∈ C, which appears quite consistent in the light of Proposition 2.1(i).

(ii) Other choices can be envisaged for g 1 and g 2 , namely

(a) g 1 = 1 2 • -x 2 + ι C and g 2 = γf (b) g 1 = 1 2 • -x 2 + γf and g 2 = ι C (c) g 1 = ι C and g 2 = 1 2 • -x 2 + γf .
Nevertheless, the strong convergence of (z m+1/2 ) m∈N in virtue of Proposition 3.6 is only guaranteed in the third case, whereas Property (30) holds only in the first case (when κ = 1 and z 0 = x). The second case was investigated in [START_REF]A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF], where the good numerical behaviour of the resulting algorithm was demonstrated.

Discussion

Both Algorithms ( 13) and ( 29) allow us to determine the proximity operator of the sum of the indicator function of a closed convex set and a function in Γ 0 (H). The main difference between the two methods is that, in the former one, the convex function needs to be differentiable with a Lipschitz-continuous gradient, whereas the latter requires that the proximity operator of the convex function is easy to compute. In addition, the forward-backward algorithm converges linearly, while we were only able to prove the strong convergence of the Douglas-Rachford algorithm. As we have shown also, the two algorithms are consistent with Proposition 2.1(i).

4 Proposed algorithms to minimize f + g + ι C

We have presented two approaches to minimize the sum of two functions in Γ 0 (H). We have also seen that these methods can be employed to compute the proximity operator of the sum of the indicator function of a closed convex set C and a function in Γ 0 (H).

We now come back to the more general form of Problem 1.1, for which we will propose two solutions. Both of them correspond to a combination of the forward-backward algorithm and the Douglas-Rachford one.

First method: insertion of a forward-backward step in the Douglas-Rachford algorithm

We propose to apply the Douglas-Rachford algorithm as described in Section 3.2, when g 1 = f and g 2 = ι C + g. If we refer to [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF], we need to determine prox κg1 = prox κf and prox κg2 = prox ιC+κg , where κ > 0.

The main difficulty lies in the computation of the second proximity operator. As proposed in Section 3.1.2, we can use a forward-backward algorithm to achieve this goal. The resulting algorithm is:

Algorithm 4.1 ➀ Set γ ∈]0, 2κ -1 β -1 [, λ ∈]0, 1] and κ ∈ ]0, +∞[. Choose (τ m ) m∈N satisfying Assumption 3.5(i). ➁ Set m = 0, z 0 = z -1/2 ∈ C. ➂ Set x m,0 = z m-1/2 . ➃ For n = 0, . . . , N m -1 a) Choose γ m,n ∈ [γ, 2κ -1 β -1 [ and λ m,n ∈ [λ, 1]. b) Compute x m,n+1 = x m,n + λ m,n P C x m,n -γ m,n (κ∇g(x m,n ) -z m ) 1 + γ m,n -x m,n . ➄ Set z m+ 1 2 = x m,Nm . ➅ Set z m+1 = z m + τ m prox κf (2z m+ 1 2 -z m ) -z m+ 1 2 . ➆ Increment m (m ← m + 1) and goto ➂.
Step ➀ allows us to set the algorithm parameters and Step ➁ corresponds to the initialization of the algorithm. At iteration m ≥ 0, Step ➃ consists of N m ≥ 1 iterations of the forward-backward part of the algorithm, where possibly varying step-sizes (γ m,n ) n and relaxation parameters (λ m,n ) n are used. Finally Steps ➄ and ➅ correspond to the Douglas-Rachford iteration. Here, the error term a m in the computation of prox κf is assumed to be equal to zero but, due to the finite number of iterations N m performed in Step ➃, an error b m = z m+1/2 -prox ιC+κg z m may be introduced in Step ➃.

It can be noticed that the forward-backward algorithm has not been initialized in Step ➂ as suggested by Proposition 3.3(ii). Indeed, as already mentioned, the computation of prox κg z m would be generally costly. Furthermore, the initialization in Step ➂ is useful to guarantee the following properties: Proposition 4.1 Suppose that Problem 1.1 has a solution and that one of Assumptions 3.7(i), 3.7(ii) or 3.7(iii) holds.

(i) Let ξ > 0 and ρ be given by [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF]. If inf g(C) > -∞ and, for every m ∈ N, the positive integer N m is chosen such that

ρ Nm √ 2κ g(z 0 ) -inf g(C) 1/2 ≤ ξ if m = 0 (31a) ρ Nm-1 1 + ξ -1 ρ 1-m z m -z m-1 ≤ 1 if m > 0 (31b)
then, (z m ) m∈N converges weakly to z ∈ H such that prox ιC +κg z is solution to Problem 1.1.

(ii) For every m ∈ N, (x m,n ) 0≤n≤Nm (and thus, z m+1/2 ) lies in C.

Proof. (i): According to Proposition 3.3(i), for every m ∈ N,

(∀n ∈ {0, . . . , N m }) x m,n -prox ιC+κg z m ≤ ρ n x m,0 -prox ιC+κg z m
and, consequently

b m = z m+1/2 -prox ιC+κg z m ≤ ρ Nm z m-1/2 -prox ιC+κg z m . (32)
Let us next show by induction that Conditions (31a) and (31b) allow us to guarantee that

b m ≤ ρ m ξ. (33)
• If m = 0, we deduce from (32) that

b 0 ≤ ρ N0 z 0 -prox ιC+κg z 0 . ( 34
)
From the definition of the proximity operator, we have

(∀x ∈ C) 1 2 z 0 -x 2 + κ g(x) ≥ 1 2 z 0 -prox ιC+κg z 0 2 + κ g prox ιC+κg z 0 ≥ 1 2 z 0 -prox ιC+κg z 0 2 + κ inf g(C)
and, since

z 0 ∈ C, z 0 -prox ιC+κg z 0 2 ≤ 2κ g(z 0 ) -inf g(C) .
By combining the latter inequality with [START_REF] Mallat | A wavelet tour of signal processing[END_REF] and (31a), we conclude that b 0 ≤ ξ.

• Now, let us show that (33) holds for m > 0, by assuming that b m-1 ≤ ρ m-1 ξ. Using (32), we have

b m ≤ ρ Nm z m-1/2 -prox ιC+κg z m-1 + prox ιC+κg z m-1 -prox ιC+κg z m ≤ ρ Nm b m-1 + prox ιC+κg z m-1 -prox ιC+κg z m ≤ ρ Nm b m-1 + z m-1 -z m
where the nonexpansivity of prox ιC+κg has been used in the last inequality. From the induction assumption, we deduce that

b m ≤ ρ Nm (ρ m-1 ξ + z m-1 -z m
which, according to (31b), leads to [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF].

Then, [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF] allows us to claim that Assumption (3.5)(ii) is satisfied since

m∈N τ m ( a m + b m ) ≤ 2ξ(1 -ρ) -1 .
By further noticing that Assumption 3.7 is equivalent to (i) (dom (ι

C + g) ∩ int dom f = ∅), (ii) (dom f ∩ int dom (ι C + g) = ∅) or, (iii) H is finite dimensional and (rint dom f ∩ rint dom (ι C + g) = ∅
), the conditions for the weak convergence of the Douglas-Rachford algorithm are therefore fulfilled.

(ii): The property can be proved by induction by noticing that x 0,0 = z -1/2 ∈ C and that x m,n+1 is a convex combination of x m,n and the projection onto C of an element of H.

Eqs. (31a) and (31b) constitute more a theoretical guaranty for the convergence of the proposed algorithm than a practical guideline for the choice of N m . In our numerical experiments, these conditions were indeed observed to provide overpessimistic values of the number of forward-backward iterations to be applied in Step ➃.

As a consequence of Proposition 4.1(ii), in Step ➃b), the gradient of g is only evaluated on C. This means that the assumption of Lipschitz-continuity on the gradient of g is only required on C and therefore, the algorithm can be applied to the following more general setting: Problem 4.1 Let C be a nonempty closed convex subset of H. Let f and g be in Γ 0 (H), where g is differentiable on C with a β-Lipschitz continuous gradient for some β ∈ ]0, +∞[. 3 Find min x∈C f (x) + g(x).

Note that, in the latter problem, the function g does need to be finite.

Second method: insertion of a Douglas-Rachford step in the forwardbackward algorithm

For this method, a different association between the functions involved in Problem 1.1 is considered by setting f 1 = ι C + f and f 2 = g. Since f 2 has then a β-Lipschitz continuous gradient, we can apply the forwardbackward algorithm presented in Section 3.1.1. This requires however to compute prox γnf1 = prox ιC+γnf , which can be performed with Douglas-Rachford iterations.

Let us summarize the complete form of the second algorithm we propose to solve Problem 1.1.

3 That is there exists an open set containing C on which g is differentiable with a β-Lipschitz continuous gradient.

➁ Set n = 0, x 0 ∈ C. ➂ Set x ′ n = x n -γ n ∇g(x n ). ➃ Set z n,0 = 2prox γnf x ′ n -x ′ n . ➄ For m = 0, . . . , M n -1 a) Compute z n,m+ 1 2 = P C z n,m + x ′ n 2 . b) Choose τ n,m ∈ [τ , 2]. c) Compute z n,m+1 = z n,m + τ n,m prox γnf (2z n,m+ 1 2 -z n,m ) -z n,m+ 1 2 . d) If z n,m+1 = z n,m , then goto ➅. ➅ Set x n+1 = x n + λ n z n,m+ 1 2 -x n . ➆ Increment n (n ← n + 1) and goto ➂.
We see that Step ➄ consists of at most M n ≥ 1 iterations of the Douglas-Rachford algorithm described in Section 3.2.2, which is initialized in accordance with Proposition 3.8(ii). Steps ➂ and ➅ correspond to a forward-backward iteration. Let m n < M n be the iteration number where the Douglas-Rachford algorithm stops. The error terms involved in Step ➅ are a n = z n,mn+ 1 2 -prox ιC+γnf x n and b n = 0. The properties of the algorithm are then the following: (ii) The sequence (x n ) n∈N lies in C.

Proof. (i): Set ρ ∈]0, 1[. Let n ∈ N and (z n,m ) m∈N be defined by iterating Steps ➄a), b) and c). By invoking Proposition 3.8(i), we know that (z n,m+ 1 2 ) m∈N converges strongly to prox ιC+γnf x ′ n . This implies that there exists M n ≥ 1 such that

(∀m ∈ N) m ≥ M n -1 ⇒ z n,m+ 1 2 -prox ιC+γnf x ′ n ≤ ρ n . If M n ≥ M n , we deduce that a n = z n,mn+ 1 2 -prox ιC+γnf x ′ n ≤ ρ n
since either m n = M n -1 or the algorithm stops in Step ➄d) (in which case z n,mn is a fixed point of the recursion in Step ➄c) and z n,mn+ 1 2 = prox ιC+γnf x ′ n ). We therefore have n∈N a n < +∞ and the conditions for the weak convergence of the forward-backward algorithm are fulfilled. (ii): We have chosen x 0 in C. In addition, (∀n ∈ N) (z n,m+ 1 2 ) m lies in C and x n+1 is convex combination of x n and z n,m+ 1 2 . Hence, it is easily shown by induction that (∀n ≥ 1) x n ∈ C. Proposition 4.2(i) guarantees that, by choosing M n large enough, the algorithm allows us to solve Problem 1.1. Although this result may appear somehow imprecise regarding the practical choice of M n , it was observed in our simulations that small values of M n are sufficient to ensure the convergence.

In addition, as a direct consequence of Proposition 4.2(ii), in Step ➂, the gradient of g is only evaluated on C. This means that, similarly to Algorithm 4.1, this algorithm is able to solve Problem 4.1. In the next section, we will show that a number of image restoration problems can be formulated as Problem 4.1.

Application to a class of image restoration problems 5.1 Context

We aim at restoring an image y in a real separable Hilbert space G from a degraded observation z ∈ G. Here, digital images of size N 1 × N 2 are considered and thus G = R N with N = N 1 N 2 . Let T be a linear operator from G to G modelling a linear degradation process, e.g. a convolutive blur. The image u = T y (resp.

z = (z (i) ) 1≤i≤N ) is a realization of a real-valued random vector U = (U (i) ) 1≤i≤N (resp. Z = (Z (i) ) 1≤i≤N ).
The image U is contaminated by noise. Conditionally to U = (u (i) ) 1≤i≤N ∈ G, the random vector Z is assumed to have independent components, which are either discrete with conditional probability mass functions (µ Zi|U (i) =u (i) ) 1≤i≤N , or absolutely continuous with conditional probability density functions which also denoted by (µ Zi|U (i) =u (i) ) 1≤i≤N . In this paper, we are interested in probability distributions such that:

∈ {1, . . . , N })(∀υ ∈ R) µ Z (i) |U (i) =υ (z (i) ) ∝ exp -ψ i (υ) (35)
the functions (ψ i ) 1≤i≤N take their values in ] -∞, +∞] and satisfy the following assumption.

Assumption 5.1 There exists a nonempty subset I of {1, . . . , N } and a constant δ ∈ R such that, for all i ∈ {1, . . . , N },

(i) dom ψ i =]δ, +∞[ if i ∈ I and, dom ψ i = [δ, +∞[ if i ∈ I; (ii) if i ∈ I, then ψ i is twice continuously differentiable on ]δ, +∞[ such that inf υ∈]δ,+∞[ ψ i (υ) > -∞ and lim υ→δ υ>δ ψ i (υ) = +∞.
Its second-order derivative ψ ′′ i is decreasing and satisfies

lim υ→+∞ ψ ′′ i (υ) = 0; (iii) if i ∈ I, then there exists α i ∈ [0, +∞[ such that (∀υ ∈ [δ, +∞[) ψ i (υ) = α i υ.
From Assumptions 5.1(ii) and (iii), it is clear that the functions (ψ i ) 1≤i≤N are convex (since (∀i ∈ I)

(∀υ ∈]δ, +∞[) ψ ′′ i (υ) ≥ 0) such that lim υ→δ υ>δ ψ ′′ i (υ) = +∞ (36)
and they are lower semicontinuous (since (∀i ∈ {1, . . . , N }) lim inf υ→δ ψ i (υ) ≥ ψ i (δ)). Examples of such functions will be provided in Sections 5.3 and 5.4.

In addition, a both simple and efficient prior probabilistic model on the unknown image y is adopted by using a representation of this image in a frame [START_REF] Daubechies | Ten lectures on wavelets[END_REF][START_REF] Han | Frames, bases, and group representations[END_REF]. The frame coefficient space is the Euclidean space H = R K (K ≥ N ). We thus use a linear representation of the form:

y = F * x
where F * : H → G is a frame synthesis operator, i.e. ν Id ≤ F * • F ≤ ν Id with (ν, ν) ∈ ]0, +∞[ 2 (which implies that F * is surjective). 4 We then assume that the vector x of frame coefficients is a realization of a random vector X with independent components. Each component X (k) with k ∈ {1, . . . , K} of X, has a probability density exp(-

φ k (•))/ +∞ -∞ exp(-φ k (η)) dη where φ k is a finite function in Γ 0 (R).
Finally, we assume that we have prior information on x which can be expressed by the fact that x belongs to a closed convex set C of H. The constraint set will be assumed to satisfy:

(T C * ) ∩ dom Ψ = ∅ ( 37
)
where

C * = F * C = F * x x ∈ C and ∀u = u (i) 1≤i≤N ∈ G Ψ(u) = N i=1 ψ i u (i) .
With these assumptions, it can be shown (see [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF]) that a Maximum A Posteriori (MAP) estimate of the vector of frame coefficients x be obtained from = z (i) 1≤i≤N by minimizing in the Hilbert space

H the function f + g + ι C where ∀x = x (k) 1≤k≤K ∈ H f (x) = K k=1 φ k x (k) (38)
and

g = Ψ • T • F * . ( 39
)
We consequently have:

Proposition 5.2 Let H = R K and G = R N with K ≥ N .
Let f and g be defined by [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF] and (39), respectively, where T : G → G is a linear operator. Under Assumption 5.1 and Condition (37), then (i) f and g are in Γ 0 (H);

(ii) if f is coercive5 or dom g ∩ C is bounded, then the minimization of f + g + ι C admits a solution. In addition, if f is strictly convex on dom g ∩ C, the solution is unique.

Proof. (i): It is clear that f is a finite convex function of H. As the functions (ψ i ) 1≤i≤N are in Γ 0 (R), Ψ belongs to Γ 0 (G). In addition, by using (37), we have ran (T • F * ) ∩ dom Ψ = ∅. This allows us to deduce that dom g = ∅ and, therefore, g ∈ Γ 0 (H).

(ii): We have dom f ∩ dom g ∩ C = ∅ since dom f = H and [START_REF]Proximité et dualité dans un espace hilbertien[END_REF] shows that dom

g ∩ C = dom (Ψ • T • F * ) ∩ C = ∅.
Since f and g are in Γ 0 (H), we deduce that + g + ι C is in Γ 0 (H).

Suppose now that f is coercive. By Assumption 5.1(ii), (∀i ∈ I) inf υ∈]δ,+∞[ ψ i (υ) > -∞ whereas, due to Assumption 5.1(iii), (∀i

∈ I) inf υ∈[δ,+∞[ ψ i (υ) = α i δ. This implies that inf Ψ(G) > -∞ and, consequently, inf g(H) ≥ inf Ψ(G) > -∞. As a result, f + g + ι C ≥ f + ι C + inf g(H) is coercive. When dom g ∩ C is bounded, f + g + ι C also

Quadratic extension

If we now investigate the Lipschitz-continuity of the gradient of g, it turns out that this property may be violated since Ψ is not finite. Due to [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF], the gradient of g is not even guaranteed to be Lipschitz-continuous on int dom g.

To circumvent this problem, it can be noticed that, because of Assumption 5.1(ii) and ( 36), for all i ∈ I, there exists a decreasing function

υ i : ]0, +∞[ →]δ, +∞[ such that lim υ i (θ) = δ and (∀θ ∈ ]0, +∞[)(∀υ ∈]δ, +∞[) 0 ≤ ψ ′′ i (υ) ≤ θ ⇔ υ ≥ υ i (θ). ( 40
)
Let us now consider the function

g θ = Ψ θ • T • F * with θ ∈ ]0, +∞[, where ∀u = u (i) 1≤i≤N ∈ G Ψ θ (u) = N i=1 ψ θ,i u (i)
and the functions (ψ θ,i ) 1≤i≤N are chosen such that,

(∀υ ∈ R) ψ θ,i (υ) =        θ 2 υ 2 + ζ i,1 (θ) υ + ζ i,0 (θ) if i ∈ I and δ -ǫ(θ) ≤ υ < υ i (θ) α i υ if i ∈ I and δ -ǫ(θ) ≤ υ < δ ψ i (υ) otherwise. ( 41
)
Hereabove, ǫ : ]0, +∞[ → ]0, +∞[ is a decreasing function and,

(∀i ∈ I) ζ i,0 (θ) = ψ i υ i (θ) -υ i (θ)ψ ′ i υ i (θ) + θ 2 υ i (θ) 2 ζ i,1 (θ) = ψ ′ i υ i (θ) -θυ i (θ).
For every i ∈ I, the constants i,0 (θ) and ζ i,1 (θ) have been determined so as to guarantee the continuity of ψ θ,i and of its first order derivative at υ i (θ). Consequently, the following result can be obtained:

Proposition 5.4 Suppose that Assumption 5.1 and Condition (37) hold. Then,

(i) (∀θ ∈ ]0, +∞[) g θ ∈ Γ 0 (H). (ii) ∀(θ 1 , θ 2 ) ∈ ]0, +∞[ 2 , θ 1 < θ 2 ⇒ g θ1 ≤ g θ2 ≤ g. (iii) For every θ ∈ ]0, +∞[, if T C * ⊂]δ -ǫ(θ), +∞[ N , then g θ has a Lipschitz-continuous gradient over C with constant β θ = θ T F * 2 ≤ θν T 2 .
(iv) For every θ ∈ ]0, +∞[, if f is coercive or if dom g θ ∩ C is bounded, then the minimization of f +g θ +ι C admits a solution. In addition, if f is strictly convex on dom g θ ∩ C, then f + g θ + ι C has a unique minimizer x θ .

(v) Assume that

(a) lim θ→+∞ ǫ(θ) = 0, (b) T C * ⊂ [δ, +∞[ N , (c) f is coercive or C is bounded, (d) f is strictly convex on C.
Then, there exists θ ∈ ]0, +∞[ such that, for every θ ∈ [θ, +∞[, the minimizer

x θ of f + g θ + ι C is the minimizer of f + g + ι C .
Proof. (i) Since Ψ θ is defined and continuous on [δ -ǫ(θ), +∞[ N and, (∀i ∈ {1, . . . , N })

(∀υ ∈]δ -ǫ(θ), +∞[) ψ ′′ θ,i (υ) ≥ 0, we have Ψ θ ∈ Γ 0 (G). In addition, dom Ψ θ ∩ ran (T • F * ) ⊃ dom Ψ ∩ ran (T • F * ) = ∅. Thus, g θ ∈ Γ 0 (H).
(ii) As a consequence of ( 41) and ( 40), we have, for every i ∈ I,

(∀υ ∈]δ, υ i (θ 2 )[) ψ ′′ i (υ) > ψ ′′ θ2,i (υ) = θ 2 . So ψ ′ i -ψ ′ θ2,i is a strictly increasing function over ]δ, υ i (θ 2 )] and (∀υ ∈]δ, υ i (θ 2 )[) ψ ′ i (υ) -ψ ′ θ2,i (υ) < ψ ′ i υ i (θ 2 ) -ψ ′ θ2,i υ i (θ 2 )) = 0
which, in turn, shows that ψ i -ψ θ2,i is strictly decreasing on ]δ, υ i (θ 2 )] and

(∀υ ∈]δ, υ i (θ 2 )[) ψ i (υ) -ψ θ2,i (υ) > ψ i υ i (θ 2 ) -ψ θ2,i υ i (θ 2 )) = 0.
In addition, we know that, if (i ∈ I and υ ≤ δ) or (i ∈ I and υ < δ), then ψ i (υ) = +∞ and, if i ∈ I and υ ≥ υ i (θ 2 ) or (i ∈ I and υ ≥ δ), then ψ i (υ) = ψ θ2,i (υ). We deduce that, for all ∈ {1, . . . , N }, ψ i ≥ ψ θ2,i and, therefore g is lower bounded by g θ2 . By proceeding similarly, we have, for every i ∈ I,

(∀υ ∈ [υ i (θ 1 ), +∞[) ψ θ2,i (υ) = ψ i (υ) = ψ θ1,i (υ) (∀υ ∈]δ -ǫ(θ 2 ), υ i (θ 1 )[) ψ ′′ θ2,i (υ) > θ 1 = ψ ′′ θ1,i (υ) ⇒ (∀υ ∈]δ -ǫ(θ 2 ), υ i (θ 1 )[) ψ ′ θ2,i (υ) < ψ ′ θ1,i (υ) ⇒ (∀υ ∈ [δ -ǫ(θ 2 ), υ i (θ 1 )[) ψ θ2,i (υ) > ψ θ1,i (υ).
In addition,

(∀i ∈ {1, . . . , N })(∀υ ∈] -∞, δ -ǫ(θ 2 )[) ψ θ2,i (υ) = +∞ ≥ ψ θ1,i (υ)
and

(∀i ∈ I)(∀υ ∈ [δ -ǫ(θ 2 ), +∞[) ψ θ2,i (υ) = ψ θ1,i (υ).
This shows that Ψ θ2 ≥ Ψ θ1 and, consequently, θ2 ≥ g θ1 .

(iii): As already mentioned, dom

Ψ θ = [δ -ǫ(θ), +∞[ N . Consider O θ = (T F * ) -1 (]δ -ǫ(θ), +∞[ N ) = x ∈ H T F * x ∈]δ -ǫ(θ), +∞[ N .
O θ is an open set and, as T C * ⊂]δ -ǫ(θ), +∞[ N , we have: C ⊂ O θ . In addition, the function g θ is differentiable on O θ and its gradient is [26, Chap. 1, Prop. 5.7]

(∀x ∈ O θ ) ∇g θ (x) = F T * ∇Ψ θ (T F * x) (42)
where

∀u = (u (i) ) 1≤i≤n ∈]δ -ǫ(θ), +∞[ N ∇Ψ θ (u) = ψ ′ θ,i (u (i) ) 1≤i≤N . We have then ∀u = (u (i) ) 1≤i≤n ∈]δ -ǫ(θ), +∞[ N ∀v = (v (i) ) 1≤i≤n ∈]δ -ǫ(θ), +∞[ N ∇Ψ θ (u) -∇Ψ θ (v) = N i=1 ψ ′ θ,i (u (i) ) -ψ ′ θ,i (v (i) ) 2 1/2
and, by the mean value theorem,

(∀i ∈ {1, . . . , N }) ψ ′ θ,i (u (i) ) -ψ ′ θ,i (v (i) ) ≤ |u (i) -v (i) | sup ξ∈]δ-ǫ(θ),+∞[ |ψ ′′ θ,i (ξ)| ≤ θ|u (i) -v (i) |.
This yields

∀u ∈]δ -ǫ(θ), +∞[ N ∀v ∈]δ -ǫ(θ), +∞[ N ∇Ψ θ (u) -∇Ψ θ (v) ≤ θ u -v
and, we deduce from (42) that

∀(x, x ′ ) ∈ O 2 θ ∇g θ (x) -∇g θ (x ′ ) ≤ θ T F * 2 x -x ′ . and T F * 2 ≤ F 2 T 2 ≤ ν T 2 .
(iv): The is similar to that of Proposition 5.2(ii). (v): In the following, we use the notation:

h = f + g + ι C and (∀θ ∈ ]0, +∞[) h θ = f + g θ + ι C .
Let (θ ℓ ) ℓ∈N be an increasing sequence of ]0, +∞[ such that lim ℓ→+∞ θ ℓ = +∞. As a consequence (i) and (ii), (h θ ℓ ) ℓ∈N is an increasing sequence of functions in Γ 0 (H). We deduce from [41, Proposition

] that (h θ ℓ ) ℓ∈N epi-converges to its pointwise limit. By using [START_REF] Rockafellar | Variational analysis[END_REF] in combination with the facts that (∀i ∈ I) lim θ→+∞ υ i (θ) = δ and lim θ→+∞ ǫ(θ) = 0, we see that the pointwise limit is equal to h. Under Assumptions (v)(b) and (v)(c), (∀ℓ ∈ N) h θ ℓ is coercive since C ∩ dom g θ ℓ = C. Equivalently, its level sets lev ≤η h θ ℓ = x ∈ H h θ ℓ (x) ≤ η with η ∈ R, are bounded. (h ℓ ) ℓ∈N being a sequence of increasing functions, ℓ∈N lev ≤η h θ ℓ = lev ≤η h θ0 bounded. As the functions h θ ℓ with ℓ ∈ N and h are lower semicontinuous and proper, [START_REF] Rockafellar | Variational analysis[END_REF]Theorem 7.33] allows us to claim that the sequence ( x θ ℓ ) ℓ∈N converges to the minimizer x of h (by Assumptions (v)(b) and (v)(d), both θ ℓ with ℓ ∈ N and h have a unique minimizer due to the strict convexity of on (C ∩ dom g) ⊂ (C ∩ dom g θ ℓ ) = C and, Propositions 5.2(ii) and 5.4(iv)). As x ∈ dom h, (∀i ∈ I) (T F * x) (i) ∈ dom ψ i =]δ, +∞[, where, for every x ∈ H and i ∈ {1, . . . , N }, (T F * x) (i) denotes the i-th component of vector T F * x. Since lim ℓ→+∞ x θ ℓ = x, we have, for every i ∈ I,

(∀η ∈ ]0, +∞[)(∃ℓ η,i ∈ N) such that ∈ N) ℓ ≥ ℓ η,i ⇒|(T F * x θ ℓ ) (i) -(T F * x) (i) | < η ⇒(T F * x θ ℓ ) (i) > min i∈I (T F * x) (i) -η. By setting = min i∈I (T F * x) (i) -δ 2 > 0 and ℓ η = max i∈I ℓ η,i , we deduce that ∈ N) ℓ ≥ ℓ η ⇒ (T F * x θ ℓ ) (i) ≥ υ (43)
where

υ = δ + min i∈I (T F * x) (i) 2 > δ. In addition, since lim ℓ→+∞ θ ℓ = +∞ ⇒ lim ℓ→+∞ max i∈I υ i (θ ℓ ) = δ,
there exists ℓ ≥ ℓ η such that (∀i ∈ I) υ i (θ ℓ ) ≤ υ. By using [START_REF] Rockafellar | Variational analysis[END_REF], this implies that (∀i

∈ I) (∀υ ∈ [υ, +∞[), ψ θ ℓ ,i (υ) = ψ i (υ). By defining now = x ∈ dom g (∀i ∈ I) (T F * x) (i) ∈ [υ, +∞[ we deduce that (∀x ∈ D) h θ ℓ (x) = h(x). Moreover, according to Assumption (v)(b), for every ℓ ∈ N, if i ∈ I, (T F * x θ ℓ ) (i) ∈ [δ, +∞[. (44)
Altogether, ( 43) and [START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF] show that both x θ ℓ and x belong to D. Consequently, as x θ ℓ = arg min x∈H h θ ℓ (x), we have:

h( x θ ℓ ) = h θ ℓ ( x θ ℓ ) ≤ h θ ℓ ( x) = h( x), which proves that x θ ℓ = x. Considering now θ ∈ [θ ℓ , +∞[, from (ii) we get: h θ ℓ ≤ h θ ≤ h. Thus, h( x) = h θ ℓ ( x) ≤ h θ ( x) ≤ h( x), which results in h θ ( x) = h( x), while (∀x ∈ H) h θ (x) ≥ h θ ℓ (x) ≥ h θ ℓ ( x) = h( x).
This allows us to conclude that x θ = x as soon as θ ≥ θ ℓ = θ.

Remark 5.5 (i) A polynomial approximation of the objective function was considered in [START_REF] Fessler | Hybrid poisson/polynomial objective functions for tomographic image reconstruction from transmission scans[END_REF] which is different from the proposed quadratic extension technique.

(ii) As expressed by Proposition 5.4(ii), g θ (resp. f + g θ + ι C ) with θ > 0 constitutes a lower approximation of g (resp. f + g + ι C ), which becomes closer as θ increases.

(iii) As shown by Proposition 5.4(iii), the main role of parameter θ is to control the Lipschitz constant of the gradient of this approximation of g.

(iv) At the same time, Proposition 5.4(v) indicates that this parameter allows us to control the closeness of the approximation to a minimizer of the original MAP criterion. This approximation becomes perfect when θ becomes greater than some value θ.

Under the assumptions of Proposition 5.4(iii), the minimization of + g θ + ι C with θ ∈ ]0, +∞[ is a problem of the type of Problem 4.1. Therefore, Propositions 4.1 and 4.2 show that, provided that f is coercive or C is bounded, Algorithms 4.1 and 4.2 can be applied in this context. In addition, Proposition 5.4(v) suggests that, by choosing θ large enough, a solution to the original MAP criterion can be found. However, according to Proposition 5.4(iii), a large value of θ induces a large value of the Lipschitz constant β θ . This means that a small value of the step-size parameter must also to be used in the forward iteration of the algorithms, which is detrimental to the convergence speed. In practice, the choice of θ results from a trade-off as will be illustrated by the numerical results.

First example

Model

We want to restore an image y ∈ [0, +∞[ N corrupted by a linear operator T : G → G and an additive noise w ∈ G, having the observation z = T y + w = u + w.

In addition, the linear operator T is assumed to be nonnegative-valued (in the sense that the matrix associated to T has nonnegative elements) and, = (w (i) ) 1≤i≤N is a realization of an independent zeromean Gaussian noise = (W (i) ) 1≤i≤N . The variance of each random variable W (i) with i ∈ {1, . . . , N } is signal-dependent and is equal to σ 2 i (u (i) ) where

(∀υ ∈ [0, +∞[) σ 2 i (υ) = υ 2α i
with α i ∈ ]0, +∞[. So, the functions (ψ i ) 1≤i≤N as defined in [START_REF]Geometrical grouplets[END_REF] are, when z (i) = 0,

(∀υ ∈ R) ψ i (υ) =    α i υ -z (i) 2 υ if υ ∈ ]0, +∞[ +∞ otherwise
and, when z (i) = 0,

(∀υ ∈ R) ψ i (υ) = α i υ if υ ∈ [0, +∞[ +∞ otherwise.
So, provided that z = 0, Assumption 5.1 is satisfied with δ = 0 and I = i ∈ {1, . . . , N } z (i) = 0 since, for all i ∈ I,

(∀υ ∈ ]0, +∞[) ψ ′ i (υ) = α i υ 2 -(z (i) ) 2 υ 2 ψ ′′ i (υ) = 2α i (z (i) ) 2 υ 3 .
We deduce from (40) that, for every i ∈ I,

(∀θ ∈ ]0, +∞[) υ i (θ) = 2α i (z (i) ) 2 θ 1/3 .

Simulation results

Here, T is either a 3 × 3 or a 7 × 7 uniform convolutive blur with T = 1. The 512 × 512 satellite image y (N = 512 2 ) shown in Fig. 1(a) has been degraded by T and a signal-dependent additive noise following the model described in the previous section with α i ≡ 1 or α i ≡ 5. The degraded image z displayed in Fig. 1(b) corresponds to a 7 × 7 uniform blur and α i ≡ 1.

A twice redundant dual-tree tight frame representation [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF] (ν = ν = 2, K = 2N ) using symlet filters of length 6 [START_REF] Daubechies | Ten lectures on wavelets[END_REF] has been employed in this example. The potential functions φ k are taken of the form

χ k | . | + ω k | . | p k where (χ k , ω k ) ∈ ]0, +∞[
2 and p k ∈ {4/3, 3/2, 2} are subband adaptive. These parameters have been determined by a maximum likelihood approach. The function f as defined by [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF] is therefore coercive and strictly convex (see Remark 5.3).

A constraint on the solution is introduced to take into account the range of admissible values in the image by choosing

C * = [0, 255] N . ( 45
)
Due to the form of the operator , T C * = C * and Condition ( 37) is therefore satisfied. Proposition 5.2 thus guarantees that a unique solution x to the MAP estimation problem exists. According to Proposition 5.4(iv), for every θ ∈ ]0, +∞[, a unique minimizer θ of f +g θ +ι C also exists which allows us to approximate x as stated by Proposition 5.4(v).

Since, for every θ ∈ ]0, +∞[, T C * = C * ⊂ [-ǫ(θ), +∞[ N , Proposition 5.4(iii) shows that g θ has a Lipschitz-continuous gradient over C and Algorithms 4.1 and 4.2 can be used to compute x θ . The two algorithms are subsequently tested.

On the one hand, when Algorithm 4.1 is used, the initialization is performed by setting z 0 = P C z and we choose κ ≡ 60 and τ m ≡ 1. The projection onto C = (F * ) -1 C * is P C = prox ι C * •F * which can be computed by using Proposition 2.6 with L = F * . The other parameters have been fixed to λ m,n ≡ 1 and γ m,n ≡ 0.995/(κθ), in compliance with Proposition 5.4(iii). The convergence of the algorithm is secured by Proposition 4.1 since Assumption 3.7(i) trivially holds. However, to improve the convergence profile, the following empirical rule for choosing the number N m of forward-backward iterations has been substituted for the necessary Conditions (31a) and (31b):

N m = inf n ∈ N * x m,n -x m,n-1 ≤ η (46)
with η = 10 -4 .

On the other hand, when Algorithm 4.2 is used, the parameters have been chosen as follows : λ n ≡ 1, τ n,m ≡ 1 and γ n ≡ 0.995/θ. The algorithm has been initialized by setting x 0 = P C z where the projection onto C is computed as described previously. The convergence of the algorithm is ensured by Proposition The number M n of Douglas-Rachford iterations has been fixed as follows:

M n = inf m ∈ N * z n,m -z n,m-1 ≤ η (47)
with the same value of η as for the first algorithm.

The error between an image y and the original image y is evaluated by the signal to noise ratio (SNR) defined as 20 log 10 ( y / y -y ).

Three objectives are targeted in our experiments. First, we want to study the performance of the proposed approach, using the redundant dual-tree transform (DTT). The results presented in Tab. As suggested by Proposition 5.4(v), as θ increases, the image is better restored. The effectiveness of the proposed approach is also demonstrated visually in Fig. 1(c) showing the restored image when T is a 7 × 7 uniform blur, i ≡ 1 and θ = 0.05. It can be observed that the algorithm allows us to recover most of the which were not perceptible due to blur and noise.

Secondly, we aim at comparing the two proposed algorithms in terms of convergence for a given value of θ. In Fig. 2, the MAP criterion value is plotted as a function of the computational time for a 7 × 7 blur, α i ≡ 5 and θ = 0.25. For improved readibility, the criterion has been normalized by subtracting the final value and dividing by the initial one. It can be noticed that Algorithm 4.2 converges faster than Algorithm 4.1. This fact was confirmed by other simulation results performed in various contexts. Finally, Fig. 3 illustrates the influence of the choice of the parameter θ when Algorithm 4.2 is used for a 7 × 7 blur and α i ≡ 5. As expected, the larger θ is, the slower the convergence of the algorithms is. A trade-off has therefore to be made: θ must be chosen large enough to reach a good restoration quality but should not be too large in order to get a fast convergence.

Simulation results

Here, T is a 5 × 5 uniform blur with T = 1. A 256 × 256 (N = 256 2 ) medical image y shown in Fig. 5(a) is degraded by and corrupted by a Poisson noise following the model described in the previous section for various intensity levels. The degraded image z is displayed in Fig. 5(b) when α i ≡ 0.01.

An orthonormal wavelet basis representation has been adopted using symlets of length 6 (ν = ν = 1, K = N ). The potential functions φ k are taken of the same form as in the first example and, the function f is therefore coercive and strictly convex.

The constraint imposed on the solution is given by C = (F * ) -1 C * where C * is defined by [START_REF] Vonesch | A fast thresholded Landweber algorithm for wavelet-regularized multidimensional deconvolution[END_REF]. Since T C * = C * , Proposition 5.4(iv) guarantees that a unique minimizer x θ of f + g θ + ι C exists, which has been computed with Algorithm 4.1. The algorithm has been initialized by setting z 0 = P C z and, we have chosen γ m,n ≡ 1.99/(κθ), κ = 60 and λ m,n ≡ τ m ≡ 1. The number of forward-backward iterations is given by (46) with η = 10 -4 . Note that the convergence rate could be accelerated by using adaptive step-size methods such as the Armijo-Goldstein search [START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF][START_REF]A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF]. However, the computational time of the step-size should be taken into account.

To evaluate the performance of our algorithm we use the Signal to Noise Ratio defined in Section 5.3.2. Tab. 2 shows the values of the SNR obtained for different values of i and θ. As predicted by Proposition 5.4(v), beyond some value of θ, which is dependent of i , the optimal value is found. We also compare our results with those provided by two different approaches. The first one is the regularized Expectation Maximization (EM) approach (also sometimes called SMART) [START_REF] Byrne | Iterative image reconstruction algorithms based on cross-entropy minimization[END_REF][START_REF] Lange | A theoretical study of some maximum likelihood algorithms for emission and transmission tomography[END_REF] where the Poisson anti-likelihood penalized by a term proportional to the Kullback-Leibler divergence between the desired solution and a reference image is minimized. Its weighting factor has been adjusted manually so as to maximize the SNR and, the reference image is a constant image whose pixel values has been set to the mean value of the degraded image. The other approach is the method based on the Anscombe transform proposed in [START_REF]A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF] and discussed in Remark 5.6. For fair comparisons, the method here employs the same orthonormal wavelet representation, the same functions (φ k ) 1≤k≤K as ours and the same constraint set C. It can be observed that the approach we propose gives good results. However, for high intensity levels (α i ≥ 0.1), the method based on the Anscombe transform performs equally well in terms of SNR. The restored images are shown in Fig. 5, when α i ≡ 0.01 and θ ≡ 0.001 after 3000 iterations. In spite of an important degradation of the original image, it can be seen that our approach is able to recover the main features in the image. It can also be noticed that the image restored by the two methods exhibit different visual characteristics.

Regularized Anscombe

Quadratic extension α i EM θ = 0.001 θ = 0.005 θ = 0. The first one concerns the minimization on a convex set C of a sum of two functions, one of which g being smooth while the other may be nonsmooth. Such a constrained minimization has been performed by forward-backward and Douglas-Rachford iterations. Various combinations of these algorithms can be envisaged and the study we made tends to show that Algorithm 4.2 is a good choice. It can be noticed that adding a constraint on the solution for a restoration problem was shown to be useful in another work [START_REF] Pustelnik | A constrained forward-backward algorithm for image recovery problems[END_REF], where it appeared that the visual quality of the restored image can be much improved w.r.t. the unconstrained case, when both restoration approaches are applicable.

The second point concerns the quadratic lower approximation technique we have proposed. This method offers a means of applying the proposed algorithms in cases when g is differentiable on C but the gradient of g is not necessary Lipschitz continuous on C. By quadratically extending the proposed constrained minimization algorithms can be used. This extension depends on a parameter θ which controls the precision (closeness to the solution of the original minimization problem) and the convergence speed of the algorithm. As illustrated by the simulations, the choice of this parameter should result from a trade-off. The numerical results have also shown the efficiency of the proposed methods in deconvolution problems involving a signaldependent Gaussian noise or a Poisson noise.

Finally, it may be interesting to note that nested iterative algorithms similar to those developed in this paper can be used to solve min x∈H f + g + h where H is a real separable Hilbert space, f , g and h are functions in Γ 0 (H) and g is β-Lipschitz differentiable.

A Study of Example 2.3

Let p = prox f x and q = prox f +ιC x where x ∈ H. Let g be the convex function defined by (∀y ∈ H) g(y) = 1 2 y -x 2 + 1 2 y ⊤ Λy. Consequently, p = I + Λ -1 x is the minimizer of g on H, whereas q is the minimizer of g on C. Thus, we can write (∀y ∈ H) g(y) = g(y) + h x where g(y) = 1 2 (y -p) ⊤ (I + Λ)(y -p) and h x is a function of x. Then, q also minimizes g on C.

In the example, we have chosen x = 2(Λ 1,2 , 1 + Λ 2,2 ) ⊤ which yields p = (0, 2) ⊤ and P C (p) = (0, 1) ⊤ . Let q = (π, 1) ⊤ . To show that q = q, have check that q minimizes g on C. A necessary and sufficient condition for the latter property to be satisfied [30, p. 293, Theorem 1.1.1] is that (∀y ∈ C) ∇g(q) ⊤ (y -q) ≥ 0 where ∇g(q) = (I + Λ)(q -p) is the gradient of g at q. This is equivalent to prove that (1) , y (2) ) ⊤ ∈ C) (2π -Λ 1,2 )(y (1) -π) + (Λ 1,2 π -Λ 2,2 -1)(y (2) -1) ≥ 0.

(51) Three cases must be considered:

• when Λ 1,2 < -2, (y (1) , y (2) ) ⊤ ∈ C ⇒ y (1) ≥ -1 = π and y (2) ≤ 1. In addition, we have 2π -Λ 1,2 = -2 -Λ 1,2 > 0 and Λ 2,2 -Λ 2 1,2 ≥ 0 ⇒ Λ 1,2 π -Λ 2,2 -1 ≤ -Λ 2 1,2 -Λ 1,2 -1 < 0. So, (51) holds.

• When Λ 1,2 > 2, similar arguments hold.

• When Λ 1,2 ∈ [-2, 2], 2π -Λ 1,2 = 0 and Λ 1,2 π -Λ 2,2 -1 =

Λ 2 1,2 2 -Λ 2,2 -1 ≤ - Λ 2 1,2
2 -1 ≤ 0, which shows that (51) is satisfied. This leads to the conclusion of Example 2.3.

B Study of Example 2.4

Let f be the function defined in Example 2.4. Defining the rotation matrix = 1 addition,

C = {x ∈ R 2 | Rx ∈ [-1, 1] 2 } = R ⊤ [-1, 1] 2 .
It can be noticed that [-1, 1] 2 is the separable convex set considered in Example 2.3 whereas appears as a particular case in the class of quadratic functions considered in this example (by setting Λ 2,2 = 1). Thus, the proximity operator of is (∀x ∈ H) prox f x = arg min , where the expression of π is given by [START_REF] Bredies | Linear convergence of iterative soft-thresholding[END_REF]. It can be concluded that P C (prox f x) = prox f +ιC x.

Find min x∈C f

 x∈C (x) + g(x).Problem 1.1 is equivalent to minimizing f + g + ι C , where ι C denotes the indicator function of C, i.e.

Algorithm 4. 2 ➀

 2 Choose sequences (γ n ) n∈N and (λ n ) n∈N satisfying Assumptions 3.1(i) and (ii). Set τ ∈]0, 2].

Proposition 4 . 2

 42 Suppose that Problem 1.1 has a solution and one of the Assumptions 3.7(i), 3.7(ii) or 3.7(iii) holds.(i) There exists a sequence of positive integers (M n ) n∈N such that, if (∀n ∈ N) M n ≥ M n then, (x n ) n∈N converges weakly to a solution to Problem 1.1.

Figure 1 :Figure 2 :

 12 Figure 1: Results for a satellite image of the city of Marseille. (a) Original image, (b) degraded image, (c) restored using a DTT.

Figure 3 :Figure 4 :

 34 Figure 3: Normalized MAP criterion (for θ = 0.25 in green and θ = 10 in magenta) versus computational time (in seconds) (Intel Xeon 4 Core, 3.00 GHz).

Figure 5 :

 5 Figure 5: Results on the medical image. (a) Original, (b) degraded, (c) restored with EM, (d) restored with Anscombe transform and (e) restored with quadratic extension.

2 ( 2 + Λ 1 , 2 , 2 - 1 √ 2 ( 1 , 1 √ 2 ( 1 +

 22122121121 Ry 2 + f (Ry) = R ⊤ prox f (Rx). and P C (prox f x) = R ⊤ P [-1,1] 2 (Rprox f x) = R ⊤ P [-1,1] 2 prox f (Rx) . Similarly, we have (∀x ∈ H) prox f +ιC x = R ⊤ prox f +ι [-1,1] 2 (Rx). So, if x = 2R ⊤ (Λ 1,2 , 2) ⊤ = √ Λ 1,2 ) ⊤ , we deduce from Example 2.3 that C (prox f x) = -1) ⊤ and prox f +ιC x = π + 1, 1 -π) ⊤

  is coercive. The existence of a solution to the minimization problem follows from classical results in convex analysis[START_REF] Ekeland | Convex analysis and variational problems[END_REF] Chap. 3, Prop. 1.2]. When f is strictly convex on dom g ∩ C, the uniqueness of the solution follows from the fact that f + g + ι C is strictly convex[START_REF] Ekeland | Convex analysis and variational problems[END_REF] Chap. 3, Prop. 1.2]. The function f is coercive (resp. strictly convex) if and only if the functions (φ k ) 1≤k≤N are coercive [14, Prop. 3.3(iii)(c)] strictly convex).

	Remark 5.3

Table 1 :

 1 1 have been generated by Algorithm 4.1, but Algorithm 4.2 leads to the same results. SNR for the satellite image.

	θ	3 × 3 blur 0.025 0.05 5	7	7 × 7 blur 0.025 0.05 5	7
	α i = 1 SNR 13.9 16.3 16.8 16.8	10.9 11.9 12.1 12.1
	θ	0.15 0.25	10	12	0.15 0.25	10	12
	α i = 5 SNR 15.9 18.0 18.8 18.8	12.6 13.3 13.7 13.7

An operator is strictly contractive with constant β if it is β-Lipschitz continuous and β ∈]0, 1[.

The interior (resp. relative interior) of a set S of H is designated by int S (resp. rint S).

The existence of the lower bound implies the existence of the upper bound in finite dimensional case.

This means that lim x →+∞ f (x) = +∞.

* Part of this work appeared in the conference proceedings of EUSIPCO 2008 [40]. This work was supported by the Agence Nationale de la Recherche under grant ANR-05-MMSA-0014-01.

Second example

Model

In this second scenario, we want to restore an image y ∈ [0, +∞[ N which is corrupted by a linear operator T : G → G, assumed to be nonnegative-valued and, which is embedded in (possibly inhomogeneous) Poisson noise. Thus, the observed image z = (z (i) ) 1≤i≤N ∈ N N is Poisson distributed, its conditional probability mass function being given by

where

Consequently, and (48), for every ∈ {1, . . . , N }, we have, when

and, when

As the functions (ψ i ) 1≤i≤N are defined up to additive constants, these constants have been chosen in (49) so as to obtain the expression of the classical Kullback-Leibler divergence term [START_REF] Byrne | Iterative image reconstruction algorithms based on cross-entropy minimization[END_REF].

In this context, provided that z = 0, Assumption 5.1 holds with = 0 and I = i ∈ {1, . . . , N } z (i) > 0 since, for all i ∈ I,

We deduce from (40) that, for every i ∈ I,

Remark 5.6 At this point, it may be interesting to compare the proposed extension with the approach developed in [START_REF]A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF]. The use of the Anscombe transform [START_REF] Anscombe | The transformation of Poisson, binomial and data[END_REF], in [START_REF]A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF] is actually tantamount to approximating the anti log-likelihood ψ i of the Poisson distribution by

The proposed quadratic extension is illustrated in Fig. 4 where a graphical comparison with the Anscombe approximation is performed.