Emilie Chouzenoux
email: emilie.chouzenoux@univ-paris-est.fr.

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

Audrey Repetti
email: audrey.repetti@univ-paris-est.fr.

J.-C Pesquet

Variable Metric Forward-Backward algorithm for minimizing the sum of a differentiable function and a convex function

Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function

Emilie Chouzenoux, Jean-Christophe Pesquet, Audrey

Introduction

We consider the following problem:

Find x ∈ Argmin G, (1)
where G : R N → (-∞, +∞] is a coercive (i.e. lim x →+∞ G(x) = +∞) function. In addition, we assume that G can be split as

G = F + R, (2)
where F is a differentiable function and R is a proper lower semicontinuous convex function. A standard approach in this context consists of using the proximal Forward-Backward (FB) algorithm [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF][START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF], which generates a sequence (x k ) k∈N by the following iterations:

x 0 ∈ R N For k = 0, 1, . . . y k = prox γ k R (x k -γ k ∇F (x k )), x k+1 = x k + λ k (y k -x k ), (3)
where, for every k ∈ N, (γ k , λ k ) ∈ (0, +∞) 2 , ∇F (x k ) is the gradient of F at x k , and prox γ k R denotes the so-called proximity operator of γ k R. Let us introduce the weighted norm:

∀x ∈ R N x U = x ⊤ U x 1/2 , (4)
where U ∈ R N ×N is some symmetric positive definite matrix. Then, the proximity operator ( [3, Sec. XV.4], [START_REF] Combettes | Variable metric quasi-Fejér monotonicity[END_REF]) is defined as follows:

Definition 1.1. Let ψ : R N → (-∞, +∞] be a proper, lower semicontinuous, convex function, let U ∈ R N ×N be a symmetric positive definite matrix, and let x ∈ R N . The proximity operator of ψ at x relative to the metric induced by U is the unique minimizer of ψ + 1 2 • -x 2 U , and it is denoted by prox U,ψ (x). If U is equal to I N , the identity matrix of R N ×N , then prox ψ ≡ prox I N ,ψ is the proximity operator originally defined in [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF].

When F is a convex function having an L-Lipschitzian gradient with L > 0, the sequence (x k ) k∈N given by (3) converges to a solution to Problem [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF], under the following assumptions on the step lengths (γ k ) k∈N and the relaxation parameters (λ k ) k∈N [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]:

• 0 < inf l∈N γ l ≤ sup l∈N γ l < 2L -1 , • (∀k ∈ N) 0 < inf l∈N λ l ≤ λ k ≤ 1.
The convergence properties of the FB algorithm have been recently extended to the case of non-convex functions F and R in [8,9] when λ k ≡ 1. The convergence results in [8,9] mainly rely on the assumption that the objective function G satisfies the Kurdyka-Lojasiewicz (KL) inequality [START_REF] Lojasiewicz | Une propriété topologique des sous-ensembles analytiques réels[END_REF]. The interesting point is that this inequality holds for a wide class of functions. In particular, it is satisfied by real analytic functions, semi-algebraic functions and many others [START_REF] Lojasiewicz | Une propriété topologique des sous-ensembles analytiques réels[END_REF][START_REF] Bolte | The Lojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems[END_REF][START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF][START_REF] Kurdyka | w f -stratification of subanalytic functions and the Lojasiewicz inequality[END_REF].

In the case of large scale optimization problems such as those encountered in image restoration, one major concern is to find an optimization algorithm able to deliver reliable numerical solutions in a reasonable time. The FB algorithm is characterized by a low computational cost per iteration. However, as many first-order minimization methods, it may suffer from slow convergence [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF]. Two families of acceleration strategies can be distinguished in the literature. The first approach, adopted for example in the FISTA method, relies on subspace acceleration [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF][START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF][START_REF] Zibulevsky | ℓ 2 -ℓ 1 optimization in signal and image processing[END_REF][START_REF] Kowalski | Proximal algorithm meets a conjugate descent[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]. In such methods, the convergence rate is improved by using informations from previous iterates for the construction of the new estimate. Another efficient way to accelerate the convergence of the FB algorithm is based on a variable metric strategy [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF][START_REF] Bonnans | A family of variable metric proximal methods[END_REF][START_REF] Burke | A variable metric proximal point algorithm for monotone operators[END_REF][START_REF] Parente | A class of inexact variable metric proximal point algorithms[END_REF][START_REF] Lotito | A class of variable metric decomposition methods for monotone variational inclusions[END_REF][START_REF] Combettes | Variable metric forward-backward splitting with applications to monotone inclusions in duality[END_REF][START_REF] Becker | A quasi-Newton proximal splitting method[END_REF]. The underlying metric of FB is modified at each iteration, giving rise to the so-called Variable Metric Forward-Backward (VMFB) algorithm:

x 0 ∈ R N For k = 0, 1, . . . y k = prox γ -1 k A k ,R (x k -γ k A -1 k ∇F (x k )), x k+1 = x k + λ k (y k -x k ), (5)
where, for every k ∈ N, A k ∈ R N ×N is a symmetric positive definite matrix. On the one hand, when A k is the identity matrix, the FB algorithm (3) is recovered. On the other hand, when R ≡ 0, Algorithm (5) corresponds to a preconditioned gradient algorithm. If F is a twice differentiable convex function, the preconditioning matrix A k is then usually chosen as an approximation of the Hessian of F at x k . This amounts to performing a change of variables leading to a function whose Hessian has more clustered eigenvalues [START_REF] Bertsekas | Nonlinear Programming, 2nd edn[END_REF]Sec.1.3.], [START_REF] Chen | Matrix Preconditioning Techniques and Applications[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF]). A convergence analysis of Algorithm ( 5) is provided in [START_REF] Combettes | Variable metric forward-backward splitting with applications to monotone inclusions in duality[END_REF], under the assumptions that F and R are convex functions and that there exists a positive bounded sequence (η k ) k∈N such that, for every k ∈ N,

(∀x ∈ R N ) (1 + η k )(x ⊤ A k+1 x) ≥ x ⊤ A k x. (6)
More specific convergence results are available in the literature in the particular case when R is the indicator function of a convex set [START_REF] Bertsekas | Projected Newton methods for optimization problems with simple constraints[END_REF][START_REF] Birgin | Nonmonotone spectral projected gradient methods on convex sets[END_REF][START_REF] Bonettini | A scaled gradient projection method for constrained image deblurring[END_REF]. However, in the aforementioned works, the convergence study is limited to the case of a convex smooth function F . As pointed out in [START_REF] Becker | A quasi-Newton proximal splitting method[END_REF], for an arbitrary matrix A k , the proximal step in ( 5) is not explicit in general, and sub-iterations are thus needed. Our contribution in this paper is to derive an inexact version of the VMFB algorithm, based on majorize-minimize arguments. The convergence of this algorithm is established for a non necessarily convex smooth function F . The rest of the paper is organized as follows: Section 2 introduces the assumptions made in the paper and presents the proposed inexact VMFB strategy. In Section 3, we investigate the convergence properties of the proposed algorithm. Finally, Section 4 provides some numerical results and a discussion of the algorithm performance by means of experiments concerning image recovery problems.

Proposed optimization method 2.1 Background and assumptions

Let us first recall some definitions and notations that will be used throughout the paper. [START_REF] Mordukhovich | Variational Analysis and Generalized Differentiation[END_REF]Sec.1.3] Let ψ : R N → (-∞, +∞] be a proper function and let x ∈ dom ψ. The Fréchet sub-differential of ψ at x is the following set:

Definition 2.1. Let ψ be a function from R N to (-∞, +∞]. The domain of ψ is dom ψ = x ∈ R N | ψ(x) < +∞ . Function ψ is proper if dom ψ is nonempty. The level set of ψ at height δ ∈ R is lev ≤δ ψ = {x ∈ R N | ψ(x) ≤ δ}. Definition 2.2. [31, Def. 8.3],
∂ψ(x) =    t ∈ R N | lim inf y→x y =x 1 x -y ψ(y) -ψ(x) -(y -x) ⊤ t ≥ 0    . If x ∈ dom ψ, then ∂ψ(x) = ∅.
The sub-differential of ψ at x is defined as

∂ψ(x) = t ∈ R N | ∃y k → x, ψ(y k ) → ψ(x), t k ∈ ∂ψ(y k ) → t .
Recall that a necessary condition for x ∈ R N to be a minimizer of ψ is that x is a critical point of ψ, i.e. 0 ∈ ∂ψ(x). Moreover, if ψ is convex, this condition is sufficient.

Remark 2.1. Definition 2.2 implies that ∂ψ is closed [9]. More precisely, we have the following property:

Let (y k , tk ) k∈N be a sequence of Graph ∂ψ = (x, t) ∈ R N × R N | t ∈ ∂ψ(x) . If y k , tk converges to x, t and ψ(y k ) converges to ψ(x), then (x, t) ∈ Graph ∂ψ.
Let us introduce our notation for linear operators. S N denotes the space of symmetric matrices of R N ×N . The Loewner partial ordering on R N ×N is defined as

(∀U 1 ∈ R N ×N )(∀U 2 ∈ R N ×N ) U 1 U 2 ⇔ (∀x ∈ R N ) x ⊤ U 1 x ≥ x ⊤ U 2 x.
In the remainder of this work, we will focus on functions F and R satisfying the following assumptions: Assumption 2.1.

(i) R : R N → (-∞, +∞] is proper, lower semicontinuous and convex, and its restriction to its domain is continuous.

(ii) F : R N → R is differentiable. Moreover, F has an L-Lipschitzian gradient on dom R where L > 0, i.e.

∀(x, y) ∈ (dom R) 2 ∇F (x) -∇F (y) ≤ L x -y .

(iii) G as defined by (2) is coercive. Some comments on these assumptions which will be useful in the rest of the paper are made below.

Remark 2.2.

(i) Assumption 2.1(ii) is weaker than the assumption of Lipschitz differentiability of F usually adopted to prove the convergence of the FB algorithm [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]9]. In particular, if dom R is compact and F is twice continuously differentiable, Assumption 2.1(ii) holds.

(ii) According to Assumption 2.1(ii), dom R ⊂ dom F . Then, as a consequence of Assumption 2.1(i), dom G = dom R is a nonempty convex set.

(iii) Under Assumption 2.1, G is proper and lower semicontinuous, and its restriction to its domain is continuous. Hence, due to the coercivity of G, for every x ∈ dom R, lev ≤G(x) G is a compact set. Moreover, the set of minimizers of G is nonempty and compact.

Assumption 2.2. Function G satisfies the Kurdyka-Lojasiewicz inequality i.e., for every ξ ∈ R, and, for every bounded subset E of R N , there exist three constants κ > 0, ζ > 0 and θ ∈ [0, 1) such that ∀t(x) ∈ ∂G(x) t(x) ≥ κ|G(x) -ξ| θ , for every x ∈ E such that |G(x) -ξ| ≤ ζ (with the convention 0 0 = 0).

Note that other forms of the KL inequality can be found in the literature [START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF][START_REF] Bolte | Characterizations of Lojasiewicz inequalities and applications[END_REF].

Majorize-Minimize metric

Some matrices serving to define some appropriate variable metric will play a central role in the algorithm proposed in this work. More specifically, let (x k ) k∈N be some given sequence of dom R and let (A k ) k∈N be a sequence of matrices of S N that fulfill the following so-called majorization conditions:

Assumption 2.3.

(i) For every k ∈ N, the quadratic function defined as

(∀x ∈ R N ) Q(x, x k ) = F (x k ) + (x -x k ) ⊤ ∇F (x k ) + 1 2 (x -x k ) ⊤ A k (x -x k ), is a majorant function of F at x k on dom R, i.e., (∀x ∈ dom R) F (x) ≤ Q(x, x k ).
(ii) There exists (ν, ν) ∈ (0, +∞) 2 such that

(∀k ∈ N) ν I N A k ν I N .
The following lemma points out the existence of such a sequence of matrices:

Lemma 2.1. Suppose that Assumption 2.1 holds and, for every k ∈ N, let A k be equal to L I N , where L > 0 is the Lipschitz constant of ∇F . Then, (A k ) k∈N satisfies Assumption 2.3 with ν = ν = L.

Proof. Under Assumption 2.1, dom R is a convex set and, since F is Lipschitz-differentiable on dom R, the Descent Lemma [25, Prop.A.24] applies, yielding:

∀(x, y) ∈ (dom R) 2 F (x) ≤ F (y) + (x -y) ⊤ ∇F (y) + L 2 x -y 2 .
Consequently, when (∀k ∈ N)

A k = L I N , Assumption 2.3(i) is satisfied while Assumption 2.3(ii) obviously holds..
Although the above lemma provides a simple choice for sequence (A k ) k∈N , it is worth noticing that other choices have been investigated in the literature [START_REF] Hunter | A tutorial on MM algorithms[END_REF][START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF] for some subclasses of functions F .

Inexact Variable Metric Forward-Backward algorithm

In general, the proximity operator relative to an arbitrary metric does not have a closed form expression. To circumvent this difficulty, we propose to solve Problem 1 by introducing the following inexact version of the Variable Metric FB method:

τ ∈ (0, +∞), x 0 ∈ dom R For k = 0, 1, . . .         Find y k ∈ R N and r(y k ) ∈ ∂R(y k ) such that R(y k ) + (y k -x k ) ⊤ ∇F (x k ) + γ -1 k y k -x k 2 A k ≤ R(x k ), ∇F (x k ) + r(y k ) ≤ τ y k -x k A k , x k+1 = (1 -λ k )x k + λ k y k , (7a) (7b) (7c)
where (A k ) k∈N is a sequence of S N associated with (x k ) k∈N for which Assumption 2.3 holds. In addition, (γ k ) k∈N and (λ k ) k∈N are sequences of nonnegative reals satisfying the following two assumptions:

Assumption 2.4. (i) There exists (η, η) ∈ (0, +∞) 2 such that, for every k ∈ N, η ≤ γ k λ k ≤ 2 -η. (ii) There exists λ ∈ (0, +∞) such that, for every k ∈ N, λ ≤ λ k ≤ 1.

Assumption 2.5.

There exists α ∈ (0, 1] such that, for every k ∈ N, (ii) (λ k ) k∈N can always be chosen such that Assumption 2.5 is satisfied (by taking for every k ∈ N, λ k = α = 1).

G(x k+1 ) ≤ (1 -α)G(x k ) + αG(y k ).
(iii) Under Assumption 2.1(i), dom R is convex. Hence, for every k ∈ N, both x k and y k belong to dom R.

(iv) As already mentioned, Algorithm (7) can be viewed as an inexact version of Algorithm [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]. Indeed, let (x k ) k∈N and (y k ) k∈N be sequences generated by Algorithm (5).

On the one hand, due to the variational characterization of the proximity operator and the convexity of R, for every k ∈ N, there exists r(y k ) ∈ ∂R(y k ) such that

r(y k ) = -∇F (x k ) + γ -1 k A k (x k -y k ) (y k -x k ) ⊤ r(y k ) ≥ R(y k ) -R(x k ), (8)
which yields

R(y k ) + (y k -x k ) ⊤ ∇F (x k ) + γ -1 k y k -x k 2 A k ≤ R(x k ).
So the sufficient-decrease condition (7a) holds. On the other hand, let us assume that Assumption 2.4 holds. According to (8), we have

∇F (x k ) + r(y k ) = γ -1 k A k (y k -x k ) ≤ η -1 √ ν y k -x k A k ,
which is the inexact optimality condition (7b) with τ = η -1 √ ν.

3 Convergence analysis

Descent properties

The present section gathers some technical results concerning the behaviour of the sequences G(y k ) k∈N and G(x k ) k∈N generated by Algorithm [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF], which will be used to prove the convergence of the proposed algorithm.

Lemma 3.1. Under Assumptions 2.1, 2.3 and 2.4, there exists µ 1 ∈ (0, +∞) such that for every k ∈ N,

G(x k+1 ) ≤ G(x k ) - µ 1 2 x k+1 -x k 2 (9) ≤ G(x k ) -λ 2 µ 1 2 y k -x k 2 . ( 10
)
Proof. For every k ∈ N, the update equation (7c) yields

G(x k+1 ) = F (x k+1 ) + R ((1 -λ k )x k + λ k y k ) .
The convexity of R and Assumption 2.3(i) allow us to deduce that

G(x k+1 ) ≤ F (x k+1 ) + (1 -λ k )R(x k ) + λ k R(y k ) ≤ F (x k ) + (x k+1 -x k ) ⊤ ∇F (x k ) + 1 2 x k+1 -x k 2 A k + (1 -λ k )R(x k ) + λ k R(y k ). (11)
In addition, according to (7c),

x k+1 -x k = λ k (y k -x k ). (12)
Using (7a) and ( 12) leads to the following inequality:

(x k+1 -x k ) ⊤ ∇F (x k ) ≤ -γ -1 k λ -1 k x k+1 -x k 2 A k + λ k R(x k ) -R(y k ) . (13)
Therefore, by combining ( 11) and ( 13), we obtain

G(x k+1 ) ≤ G(x k ) -(γ -1 k λ -1 k - 1 2 ) x k+1 -x k 2 A k ≤ G(x k ) - 1 2 η 2 -η x k+1 -x k 2 A k ,
where the last inequality follows from Assumption 2.4(i). Then, the lower bound in As- As a consequence of the above lemma, Assumption 2.5 can be reexpressed in a different form:

Corollary 3.1. Let α ∈ (0,
k ∈ [α, 1] such that G(x k+1 ) ≤ (1 -α k )G(x k ) + α k G(y k ). ( 14
)
Proof. Under Assumption 2.5, ( 14) holds if we take, for every k ∈ N, α k = α. Conversely, ( 14) is equivalent to

α k (G(x k ) -G(y k )) ≤ G(x k ) -G(x k+1 ). (15)
If the above inequality holds with α k ∈ [α, 1], then two cases may arise:

(i) Case when G(x k ) ≤ G(y k ). From Lemma 3.1 we have G(x k+1 ) ≤ G(x k ). Thus, α(G(x k ) -G(y k )) ≤ 0 ≤ G(x k ) -G(x k+1 ).
(ii) Case when G(x k ) ≥ G(y k ). Then, [START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF] yields

α(G(x k ) -G(y k )) ≤ α k (G(x k ) -G(y k )) ≤ G(x k ) -G(x k+1 ).
This shows that, if (14) holds, then Assumption 2.5 is satisfied.

When G satisfies some convexity property, we recover standard assumptions on the relaxation parameter as shown below. Proof. According to (7c), we have

(∀k ∈ N) G(x k+1 ) = G (1 -λ k )x k + λ k y k , where λ k ∈ (0, 1]. If G is convex on [x k , y k ] for every k ∈ N, then G(x k+1 ) ≤ (1 -λ k )G(x k ) + λ k G(y k ).
Using Corollary 3.1 and the fact that, for every k ∈ N, λ k is lower-bounded by λ > 0, we conclude that Assumption 2.5 is satisfied.

The next result will allow us to evaluate the variations of G when going from x k to y k at each iteration k of Algorithm [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]. Lemma 3.2. Under Assumptions 2.1, 2.3 and 2.4, there exists µ 2 ∈ R such that

(∀k ∈ N) G(y k ) ≤ G(x k ) -µ 2 y k -x k 2 .
Proof. According to Assumption 2.3(i) and (7a), we have

F (y k ) ≤ Q(y k , x k ) ≤ F (x k ) + R(x k ) -R(y k ) -(γ -1 k - 1 2 ) y k -x k 2 A k ,
which, by using Assumption 2.4, yields,

G(y k ) ≤ G(x k ) -( λ 2 -η - 1 2 ) y k -x k 2 A k .
The result then follows from Assumption 2.3(ii) by setting

µ 2 =      ν( λ 2 -η - 1 2 ) if 2λ + η ≥ 2 ν( λ 2 -η - 1 2
) otherwise.

Convergence result

Our convergence proof hinges upon the following preliminary result:

Lemma 3.3. Let (u k ) k∈N , (g k ) k∈N , (g ′ k )
k∈N and (∆ k ) k∈N be sequences of nonnegative reals and let θ ∈ (0, 1). Assume that

(i) For every k ∈ N, u 2 k ≤ g θ k ∆ k . (ii) (∆ k ) k∈N is summable. (iii) For every k ∈ N, g k+1 ≤ (1 -α)g k + g ′ k where α ∈ (0, 1]. (iv) For every k ≥ k * , (g ′ k ) θ ≤ βu k where β > 0 and k * ∈ N.
Then, (u k ) k∈N is a summable sequence.

Proof. According to (iii), for every k ∈ N,

g θ k+1 ≤ (1 -α) θ g θ k + (g ′ k ) θ .
Assumption (iv) then yields

(∀k ≥ k * ) g θ k+1 ≤ (1 -α) θ g θ k + βu k , which implies that, for every K > k * , K k=k * +1 g θ k ≤ (1 -α) θ K-1 k=k * g θ k + β K-1 k=k * u k ⇔ 1 -(1 -α) θ K-1 k=k * g θ k ≤ g θ k * -g θ K + β K-1 k=k * u k . (16)
On the other hand, (i) can be rewritten as

(∀k ∈ N) u 2 k ≤ β -1 1 -(1 -α) θ g θ k β 1 -(1 -α) θ -1 ∆ k .
By using now the inequality ∀(v, v ′ ) ∈ [0, +∞) 2 √ vv ′ ≤ (v + v ′ )/2 and since, for every k ∈ N, u k ≥ 0, we get

(∀k ∈ N) u k ≤ 1 2 β -1 1 -(1 -α) θ g θ k + β 1 -(1 -α) θ -1 ∆ k . (17)
We deduce from ( 16) and ( 17) that, for every

K > k * , K-1 k=k * u k ≤ 1 2 K-1 k=k * u k + β -1 (g θ k * -g θ K ) + β 1 -(1 -α) θ -1 K-1 k=k * ∆ k ⇒ K-1 k=k * u k ≤ β -1 g θ k * + β 1 -(1 -α) θ -1 K-1 k=k * ∆ k .
The summability of (u k ) k∈N then follows from (ii).

Our main result concerning the asymptotic behaviour of Algorithm ( 7) can now be stated: Theorem 3.1. Under Assumptions 2.1-2.5, the following hold. (iii) G(x k ) k∈N and G(y k ) k∈N are sequences converging to G(x). Moreover, G(x k ) k∈N is a nonincreasing sequence.

Proof. According to Lemma 3.1, we have

(∀k ∈ N) G(x k+1 ) ≤ G(x k ),
thus, (G(x k )) k∈N is a nonincreasing sequence. In addition, by Remark 2.2(iii) and Remark 2.3(iii), the sequence x k k∈N belongs to a compact subset E of lev ≤G(x 0 ) G ⊂ dom R and G is lower bounded. Thus, G(x k ) k∈N converges to a real ξ, and G(x k ) -ξ k∈N is a nonnegative sequence converging to 0. Moreover, by invoking again Lemma 3.1, we have

(∀k ∈ N) λ 2 µ 1 2 y k -x k 2 ≤ (G(x k ) -ξ) -(G(x k+1 ) -ξ). (18)
Hence, the sequence (y k -x k ) k∈N converges to 0. On the other hand, Assumption 2.5 implies that, for every k ∈ N,

G(x k+1 ) -ξ ≤ (1 -α)(G(x k ) -ξ) + α(G(y k ) -ξ).
Then, combining the last inequality with Lemma 3.2, we obtain that, for every k ∈ N,

α -1 G(x k+1 ) -ξ -(1 -α)(G(x k ) -ξ) ≤ G(y k ) -ξ ≤ G(x k ) -ξ -µ 2 y k -x k 2 .
Thus, since (y k -x k ) k∈N and (G(x k ) -ξ) k∈N both converge to 0, the sequence (G(y k )) k∈N converges to ξ.

Let us come back to [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] and let us apply to the convex function ψ : [0, +∞) → [0, +∞) : u → u 1 1-θ , with θ ∈ [0, 1), the gradient inequality

(∀(u, v) ∈ [0, +∞) 2 ) ψ(u) -ψ(v) ≤ ψ(u)(u -v),
which, after a change of variables, can be rewritten as

(∀(u, v) ∈ [0, +∞) 2 ) u -v ≤ (1 -θ) -1 u θ (u 1-θ -v 1-θ ).
Using the latter inequality with u = G(x k ) -ξ and v = G(x k+1 ) -ξ leads to

(∀k ∈ N) G(x k ) -ξ -G(x k+1 ) -ξ ≤ (1 -θ) -1 G(x k ) -ξ θ ∆ ′ k , where (∀k ∈ N) ∆ ′ k = G(x k ) -ξ 1-θ -G(x k+1 ) -ξ 1-θ .
Thus, combining the above inequality with ( 18) yields

(∀k ∈ N) y k -x k 2 ≤ 2λ -2 µ -1 1 (1 -θ) -1 G(x k ) -ξ θ ∆ ′ k . (19)
On the other hand, since E is bounded and Assumption 2.2 holds, there exist constants κ > 0, ζ > 0 and θ ∈ [0, 1) such that Let r(y k ) be defined as in Algorithm [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]. Then, we have

∀r(x) ∈ ∂R(x) κ|G(x) -ξ| θ ≤ ∇F (x) + r(x) , (20)
κ|G(y k ) -ξ| θ ≤ ∇F (y k ) -∇F (x k ) + ∇F (x k ) + r(y k ) ≤ ∇F (y k ) -∇F (x k ) + ∇F (x k ) + r(y k ) ≤ ∇F (y k ) -∇F (x k ) + τ x k -y k A k . (21)
Thus, by using Assumptions 2.1(ii) and 2.3(ii), we get

|G(y k ) -ξ| θ ≤ κ -1 (L + τ √ ν) x k -y k . ( 22
)
In addition, according to Assumption 2.5,

G(x k+1 ) -ξ ≤ (1 -α)(G(x k ) -ξ) + |G(y k ) -ξ|. ( 23
)
Besides, it can be noticed that

+∞ k=k * ∆ ′ k = +∞ k=k * G(x k ) -ξ 1-θ -G(x k+1 ) -ξ 1-θ = G(x k * ) -ξ 1-θ ,
which shows that (∆ ′ k ) k∈N is a summable sequence. From [START_REF] Bonnans | A family of variable metric proximal methods[END_REF], the summability of (∆ ′ k ) k∈N , [START_REF] Combettes | Variable metric forward-backward splitting with applications to monotone inclusions in duality[END_REF], and ( 22), and by setting

(∀k ∈ N)                u k = y k -x k g k = G(x k ) -ξ ≥ 0 g ′ k = |G(y k ) -ξ| ≥ 0 β = κ -1 (L + τ √ ν) > 0 ∆ k = 2λ -2 µ -1 1 (1 -θ) -1 ∆ ′ k , Lemma 3.3 allows us to conclude that ( y k -x k ) k∈N is summable when θ = 0. When θ = 0, as x k -y k → 0, there exists k * * ≥ k * such that (∀k ≥ k * * ) κ -1 (L + τ √ ν) x k -y k < 1.
Hence, according to [START_REF] Lotito | A class of variable metric decomposition methods for monotone variational inclusions[END_REF] (recall that 0 0 = 0), one necessarily has, for every k ≥ k * * , G(y k ) = ξ. Then, according to [START_REF] Bonnans | A family of variable metric proximal methods[END_REF], for every k ≥ k * * , x k = y k , which trivially shows that ( y k -x k ) k∈N is summable. Moreover, according to (7c) and Assumption 2.4(ii), we have

(∀k ∈ N) x k+1 -x k = λ k y k -x k ≤ y k -x k .
Hence, the sequence (x k ) k∈N satisfies the finite length property. In addition, since this latter condition implies that (x k ) k∈N is a Cauchy sequence, it converges towards a point x. As (x k -y k ) k∈N converges to 0, (y k ) k∈N also converges to the same limit x, and (y k ) k∈N satisfies the finite length property since

(∀k ∈ N) y k+1 -y k ≤ y k+1 -x k + x k -y k ≤ y k+1 -x k+1 + x k+1 -x k + x k -y k .
It remains us to show that the limit x is a critical point of G. To this end, let us define

(∀k ∈ N) t(y k ) = ∇F (y k ) + r(y k ),
where r(y k ) is given by ( 7), so that (y k , t(y k )) ∈ Graph ∂G. In addition, by proceeding like in [START_REF] Parente | A class of inexact variable metric proximal point algorithms[END_REF], we obtain

(∀k ∈ N) t(y k ) ≤ (L + τ √ ν) x k -y k .
Since the sequences (x k ) k∈N and (y k ) k∈N both converge to x, (y k , t(y k )) k∈N converges to (x, 0). Furthermore, according to Remark 2.2(iii), the restriction of G to its domain is continuous. Thus, as (∀k ∈ N)

y k ∈ dom G, the sequence (G(y k )) k∈N converges to G(x).
Finally, according to the closedness property of ∂G (see Remark 2.1), (x, 0) ∈ Graph ∂G i.e., x is a critical point of G.

As an offspring of the previous theorem, the proposed algorithm can be shown to locally converge to a global minimizer of G: Corollary 3.3. Suppose that Assumptions 2.1-2.5 hold, and suppose that (x k ) k∈N and (y k ) k∈N are sequences generated by Algorithm [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]. There exists υ > 0 such that, if

G(x 0 ) ≤ inf x∈R N G(x) + υ, (24)
then (x k ) k∈N and (y k ) k∈N both converge to a solution to Problem (1).

Proof. According to Remark 2.2(iii),

ξ = inf x∈R N G(x) < +∞.
Let E = lev ξ+δ G, where δ > 0. As a consequence of Assumption 2.1(iii), E is bounded. In view of Assumption 2.2, there exists constants κ > 0, ζ > 0 and θ ∈ [0, 1) such that [START_REF] Burke | A variable metric proximal point algorithm for monotone operators[END_REF] holds for every x ∈ E such that |G(x) -ξ| ≤ ζ, that is, G(x) ≤ ξ + ζ since, by definition of ξ, we always have G(x) ≥ ξ. Let us now set υ = min{δ, ζ} > 0 and choose x 0 satisfying [START_REF] Becker | A quasi-Newton proximal splitting method[END_REF]. It follows from Theorem 3.1(iii) that, for every k ∈ N,

G(x k ) ≤ ξ + υ.
By continuity of the restriction of G to its domain and Theorem 3.1(i), (x k ) k∈N and (y k ) k∈N converge to x which is such that G(x) ≤ ξ + υ. In other words, the Kurdyka-Lojasiewicz inequality is satisfied at x:

(∀t(x) ∈ ∂G(x)) t(x) ≥ κ|G(x) -ξ| θ .
As x is a critical point of G, 0 belongs to ∂G(x), and then we have

|G(x) -ξ| θ ≤ 0. This shows that G(x) = inf x∈R N G(x).
We now comment on the differences between the results in some related works and our results. (i) In [START_REF] Combettes | Variable metric forward-backward splitting with applications to monotone inclusions in duality[END_REF], the convergence is established under the assumption that G is convex, while our study relies on the fulfillment of Assumption 2.2. Moreover, it can be noticed that Assumption 2.3 on the matrices (A k ) k∈N are less restrictive than Condition (6) considered in [START_REF] Combettes | Variable metric forward-backward splitting with applications to monotone inclusions in duality[END_REF].

(ii) Note that the convergence of (7) was established in [9] in the non-preconditioned case, i.e. A k ≡ L I N , but for a non necessarily convex function R. Constant values of the relaxation parameters (λ k ) k∈N and of the step sizes (γ k ) k∈N were considered.

Application to image reconstruction 4.1 Optimization problem

In this section, we consider an inverse problem where a degraded image z = (z (m) ) 1≤m≤M ∈ R M related to an original image x ∈ [0, +∞) N is observed through the model:

(∀m ∈ {1, . . . , M }) z (m) = [Hx] (m) + σ (m) ([Hx] (m) ) w (m) ,
where H ∈ R M ×N is a matrix with non-negative elements and, for every m ∈ {1, . . . , M }, [Hx] (m) denotes the m-th component of Hx. Moreover, (w (m) ) 1≤m≤M is a realization of a Gaussian random vector with zero-mean and covariance matrix I M , and

(∀m ∈ {1, . . . , M }) σ (m) : [0, +∞) → (0, +∞) u → a (m) u + b (m) (25)
with (a (m) ) 1≤m≤M ∈ [0, +∞) M , (b (m) ) 1≤m≤M ∈ (0, +∞) M . Such a noise model arises in a number of digital imaging devices [START_REF] Healey | Radiometric CCD camera calibration and noise estimation[END_REF][START_REF] Tian | Analysis of temporal noise in CMOS photodiode active pixel sensor[END_REF][START_REF] Janesick | Photon Transfer[END_REF] where the acquired image is contaminated by signal-dependent Photon shot noise and by independent electrical or thermal noise. Signaldependent Gaussian noise can also be viewed as a second-order approximation of Poisson-Gauss noise which is frequently encountered in astronomy, medicine and biology [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Li | A reweighted ℓ 2 method for image restoration with Poisson and mixed Poisson-Gaussian noise[END_REF].

Our objective is to produce an estimate x ∈ [0, +∞) N of the target image x from the observed data z.

The original image can be estimated by solving [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF] where F is a so-called data fidelity term and R is a penalty function serving to incorporate a priori information. In the Bayesian framework, this is equivalent to compute a maximum a posteriori (MAP) estimate [START_REF] Sorenson | Parameter Estimation : Principles and Problems, Control and systems theory[END_REF] of the original image. In this context, a usual choice for the data fidelity term is the neg-log-likelihood of the data which is expressed as

(∀x ∈ R N ) F (x) = F 1 (x) + F 2 (x) if x ∈ [0, +∞) N +∞ otherwise,
where

(∀x ∈ [0, +∞) N ) F 1 (x) = 1 2 M m=1 ρ (m) 1 ([Hx] (m) ), (26)
F 2 (x) = 1 2 M m=1 ρ (m) 2 ([Hx] (m) ), (27)
and

(∀m ∈ {1, . . . , M })(∀u ∈ [0, +∞)) ρ (m) 1 (u) = 1 2 u -z (m) 2 a (m) u + b (m) (28)
ρ (m) 2 (u) = 1 2 log(a (m) u + b (m) ). ( 29
)
In the equation [START_REF] Birgin | Nonmonotone spectral projected gradient methods on convex sets[END_REF], log is a semi-algebraic approximation of the logarithm defined on (0, +∞), which, like the original function, is concave and Lipschitz differentiable on any interval [b, +∞) with b ∈ (0, +∞). Such approximations are commonly used in numerical implementations of the logarithmic function [START_REF] Dahlquist | Numerical Methods[END_REF]Chap.4]. Furthermore, a hybrid penalty function, made up of two terms R = R 1 + R 2 is considered. First, in order to take into account the dynamic range of the target image, we define R 1 = ι C , where C = [x min , x max ] N , x min ∈ [0, +∞) and x max ∈ (x min , +∞) are the minimal and the maximal values of the components of x, respectively, and ι C is the indicator function of C defined as

ι C (x) = 0 if x ∈ C +∞ otherwise.
Secondly, a sparsity prior in an analysis frame [START_REF] Chaâri | Solving inverse problems with overcomplete transforms and convex optimization techniques[END_REF][START_REF] Elad | Analysis versus synthesis in signal priors[END_REF][START_REF] Pustelnik | Relaxing tight frame condition in parallel proximal methods for signal restoration[END_REF] is introduced by taking

(∀x ∈ R N ) R 2 (x) = J j=1 ϑ (j) |[W x] (j) |,
where (ϑ (j) ) 1≤j≤J ∈ [0, +∞) J and W ∈ R J×N is a tight frame operator, i.e. there exists µ W ∈ (0, +∞) such that W ⊤ W = µ W I N .

Backward step

The implementation of the VMFB algorithm (in its exact form) requires to compute, at each iteration k ∈ N, variable y k corresponding to the proximity operator of R at xk = x k -γ k A -1 k ∇F (x k ), relative to the metric induced by γ -1 k A k , i.e.

y k = argmin x∈R N J j=1 ϑ (j) |[W x] (j) | + ι C (x) + 1 2 x -xk 2 γ -1 k A k . ( 45
)
Due to the presence of matrix W , it is not possible to obtain an explicit expression for y k . Sub-iterations are thus needed to compute it. Equation ( 45) is equivalent to

y k = γ 1/2 k A -1/2 k argmin x∈R N 1 2 x -γ -1/2 k A 1/2 k xk 2 + J j=1 ϑ (j) γ 1/2 k |[W A -1/2 k x] (j) | + ι C (γ 1/2 k A -1/2 k x) .
The above optimization problem can be solved by various algorithms. The numerical results provided in the next section have been obtained by using the Dual Forward-Backward algorithm [START_REF] Combettes | Proximity for sums of composite functions[END_REF].

Experimental results

We now demonstrate the practical performance of our algorithm on two image reconstruction scenarios. In the first scenario, the standard Peppers image of size 256 × 256 from http://sipi.usc.edu/database/ is degraded by a blur operator H corresponding to a uniform convolution kernel of size 5 × 5, and further corrupted with the considered signal-dependent additive noise with standard-deviation given by [START_REF] Bertsekas | Nonlinear Programming, 2nd edn[END_REF] where, for every m ∈ {1, . . . , M }, a (m) = 0.5 and b (m) = 1. Here, x min = 0, x max = 252, and the Lipschitz constant of F is equal to L = 1.8 × 10 4 . In the second scenario, x corresponds to one slice of the standard Zubal phantom from [START_REF] Zubal | Computerized three-dimensional segmented human anatomy[END_REF] with dimensions 128 × 128 and H is the Radon matrix modeling M = 16384 parallel projections from 128 acquisition lines and 128 angles. The sinogram Hx is corrupted with signal-dependent noise where, for every m ∈ {1, . . . , M }, a (m) = 0.01 and b (m) = 0.1. In this case, x min = 0, x max = 1, and the Lipschitz constant is L = 6.0 × 10 6 . For both experiments, the employed frame is a redundant (undecimated) wavelet transform using Daubechies' eight-taps filters over three resolution levels. The related tight frame constant is equal to µ W = 64. Parameters ϑ (j) 1≤j≤J are adjusted so as to maximize the signal-to-noise ratio (SNR) between the original image x and the reconstructed one x, which is expressed as SNR = 20 log 10 x x -x .

Figs. 1, 3 and. 4 show the degraded data, and the reconstructed images with VMFB, for the considered deblurring and reconstruction problems. We also present in Fig. 4 the reconstruction result obtained using the standard filtered back-projection approach [START_REF] Slaney | Principles of computerized tomographic imaging[END_REF]. The proposed VMFB algorithm is compared with FB [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF] and FISTA [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] algorithms. The setting λ k ≡ 1 has been chosen for FB and VMFB. Two values of the step-size are tested, namely γ k ≡ 1 and γ k ≡ 1.9. For these values of (λ k ) k∈N and (γ k ) k∈N , Assumption 2.4 is satisfied. Moreover, according to Remark 2.3(ii), Assumption 2.5 also holds. Figs. 2 and 5 illustrate the variations of G(x k ) -G(x) k and x k -x k with respect to the computation time, using the proposed VMFB algorithm, FB and FISTA algorithms, when performing tests on an Intel(R) Xeon(R) CPU X5570 at 2.93GHz, in a single thread, using a Matlab 7 implementation. Note that the optimal solution x has been precomputed for each algorithm, using a large number of iterations. We can observe that, in the deblurring experiment, FISTA converges faster than FB algorithm, while, in the tomography example, the two methods behave similarly. Concerning the choice of the step-size, our results show that, for FB and VMFB algorithms, a faster convergence is obtained for γ k ≡ 1.9 in both experiments.

In conclusion, the variable metric strategy leads to a significant acceleration in terms of decay of both the objective function and the error on the iterates in each experiment.

  (i) If, for every k ∈ N, one chooses x k+1 such that G(x k+1 ) ≤ G(y k ), then Assumption 2.5 holds for α = 1.

sumption 2 . 3 (

 23 ii) allows us to derive (9) by setting µ 1 = νη 2 -η > 0. Inequality (10) results from (12) and Assumption 2.4(ii).

Corollary 3 . 2 .

 32 Under Assumptions 2.1, 2.3 and 2.4, if G is convex on [x k , y k ] for every k ∈ N, then Assumption 2.5 holds.

( i )

 i The sequences (x k ) k∈N and (y k ) k∈N defined by (7) both converge to a critical point x of G. (ii) These sequences have a finite length in the sense that +∞ k=0 x k+1 -x k < +∞ and +∞ k=0 y k+1 -y k < +∞.

  for every x ∈ E such that |G(x) -ξ| ≤ ζ. Since G(y k ) k∈N converges to ξ, there exists k * ∈ N, such that, for every k ≥ k * , |G(y k ) -ξ| < ζ. Hence, we have, for every r(y k ) ∈ ∂R(y k ), (∀k ≥ k * ) κ|G(y k ) -ξ| θ ≤ ∇F (y k ) + r(y k ) .

Figure 1 :

 1 Figure 1: Deblurring: Degraded image, SNR=19.3 dB (left) and restored image (right) with the proposed algorithm, SNR=24.3 dB.

Figure 2 :

 2 Figure 2: Deblurring: Comparison of VMFB algorithm with γ k ≡ 1 (solid thin line) and γ k ≡ 1.9 (solid thick line), FB algorithm with γ k ≡ 1 (dashed thin line) and γ k ≡ 1.9 (dashed thick line) and FISTA algorithm (dotted line).

Figure 3 :

 3 Figure 3: Reconstruction: Original image (left) and degraded sinogram (right).

Figure 4 :Figure 5 :

 45 Figure 4: Reconstruction: Restored images with filtered back-projection, SNR=7 dB (left) and with the proposed algorithm, SNR=18.9 dB (right).

  1] and let k ∈ N. Under Assumptions 2.1, 2.3 and 2.4, Assumption 2.5 is satisfied if and only if there exists α

It is clear that Assumption 2.1(i) holds. For every m ∈ {1, . . . , M }, the first and second derivatives 1

ρ(m)

Hence, ρ(m)

1 is bounded and, as ρ(m)

is Lipschitzian, Assumption 2.1(ii) is satisfied. Assumption 2.1(iii) follows from the fact that dom R = C is bounded.

Finally, since F 1 , F 2 and R 2 are semi-algebraic functions and C is a semi-algebraic set, G is a semi-algebraic function and Assumption 2.2 holds.

Construction of the majorant

Let us now present a family of diagonal matrices (A k ) k∈N that fulfill Assumption 2.3. First, note that, on its domain, F is the sum of the convex function F 1 and the concave function F 2 , respectively defined by ( 26) and [START_REF] Nocedal | Numerical Optimization[END_REF]. For every k ∈ N, let x k be generated by the k-th iteration of Algorithm [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]. A majorant function of

where ∇F 2 (x k ) is the gradient of F 2 at x k . The next lemma allows us to construct a majorant function of F 1 at x k . Before stating this lemma, we introduce the function

where, for every m ∈ {1, . . . , M },

and ρ 1 is defined by [START_REF] Bertsekas | Projected Newton methods for optimization problems with simple constraints[END_REF].

Lemma 4.1. Let F 1 be defined by [START_REF] Chen | Matrix Preconditioning Techniques and Applications[END_REF] where

where ε ∈ [0, +∞) and P = P (m,n) 1≤m≤M,1≤n≤N is the matrix whose elements are given by

Then, Q 1 defined as

is convex and infinitely derivable on [0, +∞). Let us define

1 We consider right derivatives when u = 0.

where functions ρ(m)

and ω (m) are respectively given by ( 28) and [START_REF] Bolte | Characterizations of Lojasiewicz inequalities and applications[END_REF]. If a (m) = 0, then ρ (m) 1 is a quadratic function and we have

Let us now assume that a (m) ∈ (0, +∞). The third derivative of ρ (m) 1 is given by ...

on [0, +∞), for every u ′ ∈ (0, +∞). The second derivative of the latter function is given by

Moreover, according to second-order Taylor's formula, there exists ũ ∈ (0, u ′ ) such that

hence, using ( 33) and ( 37),

1 is strictly decreasing, and [START_REF] Janesick | Photon Transfer[END_REF]

1 is first strictly decreasing on (0, ũ), then strictly increasing on (ũ, +∞). On the one hand, ( 33) and ( 36) yield

Thus, according to the mean value theorem, there exists u * ∈ (0, u ′ ) such that q(m) (u * , u ′ )ρ(m) 1 (u * ) = 0. On the other hand, according to [START_REF] Healey | Radiometric CCD camera calibration and noise estimation[END_REF],

Therefore, from the monotonicity properties of q(m) (•, u ′ ) -ρ(m) 1 , we deduce that u * is the unique zero of this function on (0, u ′ ), and

Equation ( 40) implies that q

1 (•) is strictly increasing on [0, u * ), strictly decreasing on (u * , u ′ ) and strictly increasing on (u ′ , +∞). Thus, given [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF],

Moreover, from the expression of ω (m) (0) in [START_REF] Bolte | Characterizations of Lojasiewicz inequalities and applications[END_REF] and from [START_REF] Rockafellar | Variational Analysis[END_REF], it is easy to show that

Therefore, by gathering ( 41) and ( 42), we obtain

which, as pointed out before, is still valid when a (m) = 0. Majoration [START_REF] Chaâri | Solving inverse problems with overcomplete transforms and convex optimization techniques[END_REF] implies that, for every k ∈ N,

According to Jensen's inequality and [START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF], for every m ∈ {1, . . . , M },

Since the convexity of ρ (m) 1

for m ∈ {1, . . . , M } implies the positivity of ω (m) on [0, +∞), we deduce that

which yields the result.

It can be deduced from [START_REF] Mordukhovich | Variational Analysis and Generalized Differentiation[END_REF] and Lemma 4.1 that Assumption 2.3(i) is satisfied for

It can be further noticed that, for every m ∈ {1, . . . , M }, the derivative of ω

where ũ ∈ (0, u ′ ). Since the third-order derivative of ρ (m) 1

is negative, ω (m) is a strictly decreasing positive function. According to the expression of R, dom R = [x min , x max ] N , so that its image {Hx | x ∈ dom R} is a compact set. Thus, the function x → ω (m) ([Hx] (m) ) admits a minimum value ω

Note that, if each column of H is non-zero, we can choose ε = 0 in [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]. Otherwise, we must choose ε > 0.