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

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

A Parallel Inertial Proximal Optimization Method

Keywords: Monotone operators, convex optimization, proximal algorithms, parallel algorithms. Mathematical Subject Classifications: 47H05, 47N10, 65K10, 68W10

The Douglas-Rachford algorithm is a popular iterative method for finding a zero of a sum of two maximally monotone operators defined on a Hilbert space. In this paper, we propose an extension of this algorithm including inertia parameters and develop parallel versions to deal with the case of a sum of an arbitrary number of maximal operators. Based on this algorithm, parallel proximal algorithms are proposed to minimize over a linear subspace of a Hilbert space the sum of a finite number of proper, lower semicontinuous convex functions composed with linear operators. It is shown that particular cases of these methods are the simultaneous direction method of multipliers proposed by Stetzer et al., the parallel proximal algorithm developed by Combettes and Pesquet, and a parallelized version of an algorithm proposed by Attouch and Soueycatt.

Introduction

The use of parallel methods for convex optimization has been an active research area for about two decades [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF][START_REF] Censor | Parallel Optimization: Theory, Algorithms and Applications[END_REF][START_REF] Han | A parallel algorithm for a class of convex programs[END_REF]. With the widespread use of multicore computer architectures, it can be expected that parallel optimization algorithms will play a more prominent role. Recently, a number of proximal parallel convex optimization algorithms have appeared in the literature [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. These algorithms are especially useful for solving large-size optimization problems arising in the fields of inverse problems and imaging.

A splitting algorithm which will be subsequently designated by the simultaneous direction method of multipliers (SDMM) was recently proposed in [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF] to solve the following problem:

minimize y∈G m i=1 f i (L i y) (1.1)
where G is a real Hilbert space, for every i ∈ {1, . . . , m}, f i is a proper lower semicontinuous convex functions from a real Hilbert space H i to ]-∞, +∞], and L i belongs to B (G, H i ), the set of bounded linear operators from G to H i . This algorithm was derived from the Douglas-Rachford algorithm [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF][START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF][START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF] by invoking a duality argument. Such a duality argument can be traced back to the work in [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite elements approximations[END_REF] (see also [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Steidl | Removing multiplicative noise by Douglas-Rachford splitting methods[END_REF]). It was also pointed out in [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF] that the proposed method generalizes to a sum of more than two functions the alternating-direction method of multipliers [START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems[END_REF][START_REF] Glowinski | Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics[END_REF]. Note that parallel alternating-direction algorithms had been previously proposed in [START_REF] Eckstein | Parallel alternating direction multiplier decomposition of convex programs[END_REF] for monotropic convex programming problems and for block-separable convex optimization problems in finite dimensional Hilbert spaces. It is also worth noticing that augmented Lagrangian techniques have become increasingly popular for solving imaging problems. These approaches appeared under different names such as alternating split Bregman algorithm [START_REF] Esser | Applications of Lagrangian-based alternating direction methods and connections to split Bregman[END_REF][START_REF] Goldstein | The split Bregman method for L1-regularized problems[END_REF][START_REF] Yin | Bregman iterative algorithms for l1minimization with applications to compressed sensing[END_REF][START_REF] Zhang | Bregmanized nonlocal regularization for deconvolution and sparse reconstruction[END_REF] or split augmented Lagrangian shrinkage algorithm (SALSA) [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF][START_REF] Afonso | An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems[END_REF] in the recent image processing literature.

On the other hand, another splitting method called the parallel proximal algorithm (PPXA) was proposed in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] to minimize a sum of convex functions in a possibly infinite dimensional Hilbert space (see also extensions in [START_REF] Combettes | Iterative construction of the resolvent of a sum of maximal monotone operators[END_REF] for monotone inclusion problems). Since the splitting algorithms in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF][START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF] were derived from the Douglas-Rachford algorithm by working in a product space, a natural question is to know whether there exist connections between them.

Despite the fact that the work in [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF] puts emphasis on duality issues whereas [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF] relies on results from the theory of proximity operators, this paper shows that both algorithms are particular instances of a more general splitting algorithm. This optimization algorithm itself is a specialization of an algorithm for solving the following problem:

find y ∈ E such that m i=1 ω i L * i A i (L i y) ∩ E ⊥ = ∅ (1.2)
where E is a closed linear subspace of G, E ⊥ is its orthogonal complement and, for every i ∈ {1, . . . , m}, A i : H i → 2 H i is a maximally monotone operator, L * i is the adjoint operator of L i , and ω i ∈ ]0, +∞[. The latter algorithm for finding a zero of a sum of maximally monotone operators is derived from an extension of the Douglas-Rachford algorithm which includes inertia parameters. Note that another framework for splitting sums of maximally monotone operators was developed in [START_REF] Eckstein | General projective splitting methods for sums of maximal monotone operators[END_REF]. Other parallel methods based on a forward-backward approach were also investigated in [START_REF] Attouch | Augmented Lagrangian and proximal alternating direction methods of multipliers in Hilbert spaces -applications to games, PDEs and control[END_REF].

Convergence results concerning various inertial algorithms can be found in [START_REF] Alvarez | Weak convergence of a relaxed and inertial hybrid projection-proximal point algorithm for maximal monotone operators in Hilbert space[END_REF][START_REF] Alvarez | Convergence and asymptotic stabilization for some damped hyperbolic equations with non-isolated equilibria[END_REF][START_REF] Maingé | Inertial iterative process for fixed points of certain quasi-nonexpansive mappings[END_REF][START_REF] Maingé | Convergence theorems for inertial KM-type algorithms[END_REF][START_REF] Maingé | Regularized and inertial algorithms for common fixed points of nonlinear operators[END_REF][START_REF] Moudafi | A hybrid inertial projection-proximal method for variational inequalities[END_REF][START_REF] Moudafi | An approximate inertial proximal method using enlargement of a maximal monotone operator[END_REF][START_REF] Moudafi | Convergence of a splitting inertial proximal method for monotone operators[END_REF]. Recently, a proximal alternating direction of multipliers method was introduced in [START_REF] Attouch | Augmented Lagrangian and proximal alternating direction methods of multipliers in Hilbert spaces -applications to games, PDEs and control[END_REF] which can be viewed as an algorithm within this class. Applications to game theory, PDEs and control were described in [START_REF] Briceño-Arias | Convex variational formulation with smooth coupling for multicomponent signal decomposition and recovery[END_REF].

Conditions for the convergence of inexact forms of the algorithms presented in this paper will be given in the following. These convergence results are valid in any finite or infinite dimensional Hilbert space. First versions of the algorithms are provided, the convergence of which requires that m i=1 ω i L * i L i is an isomorphism. This assumption is relaxed in slightly more complex variants of the algorithms. The resulting unrelaxed optimization method is shown to correspond to a parallelized version of the algorithm in [START_REF] Attouch | Augmented Lagrangian and proximal alternating direction methods of multipliers in Hilbert spaces -applications to games, PDEs and control[END_REF] for equal values of the inertia parameters.

In Section 2, we introduce our notation and recall some useful properties of monotone operators. In Section 3, we propose a variant of the Douglas-Rachford algorithm in a product space, which includes inertia parameters, and study its convergence. In Section 4, we deduce from this algorithm an inertial algorithm for solving Problem (1.2). In Section 5, we consider an application of this algorithm to convex minimization problems and examine the connections of the resulting parallel convex optimization methods with some existing approaches.

Notation

Let (H i , • i ) 1≤i≤m be real Hilbert spaces. We define the product space H = H 1 × • • • × H m endowed with the norm

||| • ||| : (x 1 , . . . , x m ) → m i=1 ω i x i 2 i (2.1)
where (ω i ) 1≤i≤m ∈]0, +∞[ m . The associated scalar product is denoted by • | • whereas, for every i ∈ {1, . . . , m}, • | • i denotes the scalar product of H i . Let k ∈ {0, . . . , m} and let us define the following closed subspace of H:

K k = {(x 1 , . . . , x k , 0, . . . , 0) ∈ H} (2.2)
(with the convention (x 1 , . . . , x k , 0, . . . , 0) = 0 if k = 0). In the following (x 1 , . . . , x m ) ∈ H 1 × • • • × H m will denote the components of a generic element x of H.

An operator

A : H → 2 H is monotone if (∀(x, y) ∈ H 2 )(∀(u, v) ∈ A(x) × A(y)) x -y | u -v ≥ 0. (2.3)
For more details concerning the properties of monotone operators and the definitions recalled below the reader is referred to [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF].

The set of zeros of A :

H → 2 H is zer A = x ∈ H 0 ∈ A(x) and its graph is gra A = (x, u) ∈ H 2 u ∈ Ax . A : H → 2 H is maximally monotone if, for every (x, u) ∈ H 2 , (x, u) ∈ gra A ⇔ (∀(y, v) ∈ gra A) x -y | u -v ≥ 0.
(2.4)

Recall that an operator T :

H → H is • β-cocoercive with β ∈ ]0, +∞[ if (∀(x, y) ∈ H 2 ) β|||T x -T y||| 2 ≤ T x -T y | x -y ; (2.5)
• firmly nonexpansive if it is 1-cocoercive;

• 1/β-Lipschitz continuous with β ∈ ]0, +∞[ if (∀(x, y) ∈ H 2 ) β|||T x -T y||| ≤ |||x -y|||; (2.6) • nonexpansive if it is 1-Lipschitz continuous; • 1/β-strictly contractive if it is 1/β-Lipschitz continuous with β ∈]1, +∞[.
If A is maximally monotone, then its resolvent J A = (Id +A) -1 is a firmly nonexpansive operator from H to H, and its reflection R A = 2J A -Id is a nonexpansive operator.

Let Γ 0 (H) be the set of proper (i.e. not identically equal to +∞) lower semicontinuous convex functions from H to ]-∞, +∞] and let f ∈ Γ 0 (H). Then, its subdifferential ∂f is a maximally monotone operator and the proximity operator of f is prox f = J ∂f [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF][START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF]. The domain of

f is dom f = x ∈ H f (x) < +∞ . The conjugate of f is f * : H → ]-∞, +∞] : x → sup y∈H x | y -f (y).
Let C be a nonempty closed convex subset of H. Its indicator function ι C ∈ Γ 0 (H) is defined as

(∀x ∈ H) ι C (x) = 0 if x ∈ C +∞ otherwise.
(2.7)

The projection onto

C is P C = prox ι C and the normal cone operator to C is N C : H → 2 H defined as (∀x ∈ H) N C (x) = u ∈ H (∀y ∈ C) y -x | u ≤ 0 if x ∈ C ∅ otherwise. (2.8)
The relative interior (resp. strong relative interior) of C is denoted by ri C (resp. sri C).

The weak convergence (resp. strong convergence) is denoted by ⇀ (resp. →).

An extension of the Douglas-Rachford algorithm

We will first consider the following problem:

find x ∈ zer(A + B) = x ∈ H 0 ∈ A(x) + B(x) (3.1)
where B : H → 2 H is a maximally monotone operator,

A : H → 2 H : (x 1 , . . . , x m ) → A 1 (x 1 ) × • • • × A m (x m ) (3.2)
and, for every i ∈ {1, . . . , m}, A i : H i → 2 H i is a maximally monotone operator.

We first state a quasi-Fejérian property which will be useful in the following.

Lemma 3.1 Let S be a nonempty set of H × K k . Let (t n ) n∈N be a sequence in H and let (p n ) n∈N = (p 1,n , . . . , p k,n , 0, . . . , 0) n∈N be a sequence in K k such that ∀( t, p ∈ S)(∀n ∈ N) |||t n+1 -t||| 2 + k i=1 µ i,n+1 p i,n+1 -p i 2 i 1/2 ≤ |||t n -t||| 2 + k i=1 µ i,n p i,n -p i 2 i 1/2 + ∆ n (3.3)
where, for every i ∈ {1, . .

(∀n ∈ N) |||t n -t||| 2 + k i=1 µ i,n p i,n -p i 2 i 1/2 ≤ ρ. (3.4)
Since, for every i ∈ {1, . . . , k}, (µ i,n ) n∈N is a sequence of nonnegative reals converging to a positive limit,

∃µ ∈ ]0, +∞[ , ∃n 0 ∈ N (∀i ∈ {1, . . . , k})(∀n ≥ n 0 ) µ i,n ≥ µ ω i (3.5)
It can be deduced that

(∀n ≥ n 0 ) |||t n -t||| 2 + µ|||p n -p||| 2 1/2 ≤ ρ. (3.6)
This implies that

(∀n ≥ n 0 ) |||t n -t||| ≤ ρ |||p n -p||| ≤ µ -1/2 ρ (3.7)
which shows the boundedness of the sequences (t n ) n∈N and (p n ) n∈N .

(ii) This fact follows from [START_REF] Combettes | Quasi-Fejérian analysis of some optimization algorithms[END_REF]Lemma 3.1 (ii)].

(iii) Since (t n , p n ) n∈N is bounded, it is enough to prove that this sequence cannot have two distinct weak sequential cluster points. Let ( t 1 , p 1 ) ∈ S and ( t 2 , p 2 ) ∈ S be two such cluster points. We have

(∀n ∈ N) 2 t n | t 1 -t 2 = |||t n -t 2 ||| 2 -|||t n -t 1 ||| 2 + ||| t 1 ||| 2 -||| t 2 ||| 2 and (∀i ∈ {1, . . . , k})(∀n ∈ N) 2 p i,n | p i,1 -p i,2 i = p i,n -p i,2 2 i -p i,n -p i,1 2 i + p i,1 2 i -p i,2 2
i . (3.8) Then, it can be deduced from (ii) and the fact that, for every i ∈ {1, . . . , k},

(µ i,n ) n∈N is convergent that ( t n | t 1 -t 2 + k i=1 µ i,n p i,n | p i,1 -p i,2 i
) n∈N converges to some limit ζ. In addition, there exist some subsequences (t n ℓ , p n ℓ ) ℓ∈N and (t n ℓ ′ , p n ℓ ′ ) ℓ ′ ∈N converging weakly to ( t 1 , p 1 ) and ( t 2 , p 2 ), respectively. Passing to the limit, we have thus

t 1 | t 1 -t 2 + k i=1 µ i,∞ p i,1 | p i,1 -p i,2 i = ζ = t 2 | t 1 -t 2 + k i=1 µ i,∞ p i,2 | p i,1 -p i,2 i ⇔ ||| t 1 -t 2 ||| 2 + k i=1 µ i,∞ p i,1 -p i,2 2 i = 0 (3.9)
where, for every i ∈ {1, . . . , k}, µ i,∞ = lim n→+∞ µ i,n > 0. Consequently, t 1 = t 2 and p 1 = p 2 .

We will also need the following property concerning nonexpansive operators.

Lemma 3.2 Let (γ, λ) ∈ ]0, +∞[ 2 and let (ε 1 , . . . , ε m ) ∈ [0, 1[ m be such that (∀i ∈ {1, . . . , m}) (ε i > 0 ⇔ i ≤ k). (3.10)
Let T and U λ be the operators from H × K k to H × K k defined as

(∀t ∈ H)(∀p ∈ K k ) T (t, p) = (z, P K k p ′ ) (3.11) U λ (t, p) = (1 - λ 2 )t + λ 2 z, P K k p ′ (3.12)
where

p ′ = J γ(1-ε 1 )A 1 (1 -ε 1 )t 1 + ε 1 p 1 , . . . , J γ(1-εm)Am (1 -ε m )t m + ε m p m (3.13) z = R γB 2p ′ -t . (3.14)
We have the following properties:

(i) Fix T = Fix U λ , and

( t, p) ∈ Fix T ⇔ J γA t ∈ zer(A + B) and (∀i ∈ {1, . . . , m}) p i = J γA i t i if i ≤ k 0 otherwise. ( 3.15)
(ii) T is nonexpansive operator in H × K k endowed with the weighted norm

ν : (z, p) → |||z||| 2 + 2 k i=1 ω i ε i 1 -ε i p i 2 i 1/2 . ( 3.16)
(iii) For every (t, s) ∈ H 2 and (p, q) ∈ K 2 k , let (z, p) = T (t, p) and (u, q) = T (s, q). Then,

|||(1 - λ 2 )(t -s) + λ 2 (z -u)||| 2 + λ k i=1 ω i ε i 1 -ε i p i -q i 2 i ≤ |||t -s||| 2 + λ k i=1 ω i ε i 1 -ε i p i -q i 2 i -λ k i=1 ω i ε i 1 -ε i p i -q i -p i + q i 2 i - 1 4 (2 -λ)λ|||t -s -z + u||| 2 . (3.17)
Proof.

(i) It can be first noticed that Fix T = Fix U λ . In addition, ( t, p) is a fixed point of T if and only if

     (∀i ∈ {1, . . . , m}) p i = p ′ i if i ≤ k 0 otherwise t = R γB 2 p ′ -t (3.18)
where

(∀i ∈ {1, . . . , m}) p ′ i = J γ(1-ε i )A i (1 -ε i ) t i + ε i p i . (3.19)
The latter relation yields, for every i ∈ {1, . . . , k}, (ii) For every (t, s) ∈ H 2 and (p, q) ∈ K 2 k , let (z, p) = T (t, p) and (u, q) = T (s, q). Let p ′ be given by (3.13) and q ′ be similarly defined as

(1 -ε i ) t i + ε i p ′ i -p ′ i ∈ γ(1 -ε i )A i ( p ′ i ) ⇔ t i -p ′ i ∈ γA i ( p ′ i ) ⇔ p ′ i = J γA i t i . ( 3
q ′ = J γ(1-ε 1 )A 1 (1 -ε 1 )s 1 + ε 1 q 1 , . . . , J γ(1-εm)Am (1 -ε m )s m + ε m q m .
(3.24)

For every i ∈ {1, . . . , m}, as J γ(1-ε i )A i is firmly nonexpansive, we have

p ′ i -q ′ i 2 i ≤ p ′ i -q ′ i | (1 -ε i )(t i -s i ) + ε i (p i -q i ) i ⇔ p ′ i -q ′ i 2 i ≤ (1 -ε i ) p ′ i -q ′ i | t i -s i i + ε i 2 ( p ′ i -q ′ i 2 i + p i -q i 2 i -p ′ i -q ′ i -p i + q i 2 i ) ⇔ -(1 -ε i ) p ′ i -q ′ i | t i -s i i ≤ -1 - ε i 2 p ′ i -q ′ i 2 i + ε i 2 p i -q i 2 i - ε i 2 p ′ i -q ′ i -p i + q i 2 i . (3.25)
On the other hand, we deduce from the nonexpansivity of R γB that

|||z -u||| 2 ≤ |||2(p ′ -q ′ ) -t + s||| 2 ⇔ |||z -u||| 2 ≤ 4|||p ′ -q ′ ||| 2 + |||t -s||| 2 -4 p ′ -q ′ | t -s . (3.26)
By combining this inequality with (3.25), we get

|||z -u||| 2 + 2 k i=1 ω i ε i 1 -ε i p ′ i -q ′ i 2 i ≤ |||t -s||| 2 + 2 k i=1 ω i ε i 1 -ε i p i -q i 2 i -2 k i=1 ω i ε i 1 -ε i p ′ i -q ′ i -p i + q i 2 i . (3.27)
This leads to

|||z -u||| 2 + 2 k i=1 ω i ε i 1 -ε i p ′ i -q ′ i 2 i ≤ |||t -s||| 2 + 2 k i=1 ω i ε i 1 -ε i p i -q i 2 i . (3.28)
Since p = (p ′ 1 , . . . , p ′ k , 0, . . . , 0) and q = (q ′ 1 , . . . , q ′ k , 0, . . . , 0), the above inequality implies that T is a nonexpansive operator, as stated above.

(iii) In addition, we have

|||(1 - λ 2 )(t -s) + λ 2 (z -u)||| 2 = (1 - λ 2 ) 2 |||t -s||| 2 + λ 2 4 |||z -u||| 2 + 1 2 (2 -λ)λ t -s | z -u = (1 - λ 2 )|||t -s||| 2 + λ 2 |||z -u||| 2 - 1 4 (2 -λ)λ|||t -s -z + u||| 2 (3.29)
which, combined with (3.27), yields (3.17).

We are now able to provide the main result of this section, which concerns the convergence of a generalization of the Douglas-Rachford algorithm.

Proposition 3.3 Let γ ∈ ]0, +∞[, (ε 1 , . . . , ε m ) ∈ [0,
Initialization t 0 ∈ H, p -1 ∈ H For n = 0, 1, . . .       For i = 1, . . . , m p i,n = J γ(1-ε i )A i (1 -ε i )t i,n + ε i p i,n-1 + a i,n t n+1 = t n + λ n J γB 2p n -t n + b n -p n .
(3.30)

Suppose that the following hold.

(i) zer(A + B) = ∅ (ii) There exists λ ∈]0, 2[ such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n < 2. (iii) n∈N |||a n ||| + |||b n ||| < +∞. Then, (t n ) n∈N converges weakly to t and J γA t ∈ zer(A + B).
Proof. Without loss of generality, it can be assumed that (3.10) is satisfied (up to some re-indexing). For every n ∈ N, let

p ′ n = J γ(1-ε 1 )A 1 (1 -ε 1 )t 1,n + ε 1 p 1,n-1 , . . . , J γ(1-εm)Am (1 -ε m )t m,n + ε m p m,n-1 (3.31) z n = R γB 2p ′ n -t n (3.32) t ′ n+1 = (1 - λ n 2 )t n + λ n 2 z n . (3.33)
We have then

(∀n ∈ N) (z n , P K k p ′ n ) = T (t n , P K k p n-1 ) (t ′ n+1 , P K k p ′ n ) = U λn (t n , P K k p n-1 ) (3.34)
where T and U λn are defined by (3.11)- (3.14). According to Lemma 3.2(i) and Assumption (i), Fix U λn = Fix T = ∅. Let ( t, p) be an arbitrary element of Fix T and, for every n ∈ N, define

t n = t n -t, p n-1 = p n-1 -p, t ′ n+1 = t ′ n+1 -t, and p ′ n = p ′ n -p. By applying Lemma 3.2(iii), we get ||| t ′ n+1 ||| 2 + λ n k i=1 ω i ε i 1 -ε i p ′ i,n 2 i + λ n k i=1 ω i ε i 1 -ε i p ′ i,n -p i,n-1 2 i + 1 4 (2 -λ n )λ n |||t n -z n ||| 2 ≤ ||| t n ||| 2 + λ n k i=1 ω i ε i 1 -ε i p i,n-1 2 i . (3.35)
From the triangular inequality, we have, for every n ∈ N,

||| t n+1 ||| 2 + λ n k i=1 ω i ε i 1 -ε i p i,n 2 i + λ n k i=1 ω i ε i 1 -ε i p i,n -p i,n-1 2 i + 1 4 (2 -λ n )λ n |||t n -z n ||| 2 1/2 ≤ ||| t ′ n+1 ||| 2 + λ n k i=1 ω i ε i 1 -ε i p ′ i,n 2 i + λ n k i=1 ω i ε i 1 -ε i p ′ i,n -p i,n-1 2 i + 1 4 (2 -λ n )λ n |||t n -z n ||| 2 1/2 + |||t ′ n+1 -t n+1 ||| + (2λ n ) 1/2 k i=1 ω i ε i 1 -ε i p ′ i,n -p i,n 2 i 1/2 . ( 3.36)
By straightforward calculations, the iterations of Algorithm (3.30) can be expressed as

(∀n ∈ N) p n = p ′ n + a n t n+1 = t ′ n+1 + λn 2 R γB (2p n -t n ) -z n + 2b n .
(3.37)

From the nonexpansivity of the R γB operator, it can be deduced that

(∀n ∈ N) |||t n+1 -t ′ n+1 ||| ≤ λ n 2 |||R γB (2p n -t n ) -z n ||| + λ n |||b n ||| (3.38) |||R γB (2p n -t n ) -z n ||| ≤ 2|||p n -p ′ n ||| = 2|||a n ||| (3.39)
and, by using (3.35), (3.36) becomes

||| t n+1 ||| 2 + λ n k i=1 ω i ε i 1 -ε i p i,n 2 i + λ n k i=1 ω i ε i 1 -ε i p i,n -p i,n-1 2 i + 1 4 (2 -λ n )λ n |||t n -z n ||| 2 1/2 ≤ ||| t n ||| 2 + λ n k i=1 ω i ε i 1 -ε i p i,n-1 2 i 1/2 + ∆ n (3.40)
where

∆ n = λ n (|||a n ||| + |||b n |||) + (2λ n ) 1/2 k i=1 ω i ε i 1 -ε i a i,n 2 i 1/2 . (3.41)
Since, due to Assumption (ii), (λ n ) n∈N is decreasing, (3.40) leads to

||| t n+1 ||| 2 + λ n k i=1 ω i ε i 1 -ε i p i,n 2 i 1/2 ≤ ||| t n ||| 2 + λ n-1 k i=1 ω i ε i 1 -ε i p i,n-1 2 i 1/2 + ∆ n . (3.42)
Since Assumption (iii) holds and (λ n ) n∈N is bounded, the above expression of (∆ n ) n∈N shows that it is a nonnegative summable sequence. Then, as (λ n ) n∈N converges to a positive limit, Lemma 3.1(ii) allows us to claim that (|||

t n ||| 2 + λ n-1 k i=1 ω i ε i 1-ε i p i,n-1 2
i ) n≥1 is a convergent sequence. In turn, as lim n→+∞ ∆ n = 0, we deduce from (3.40) and (3.42) that

λ k i=1 ω i ε i 1 -ε i p i,n -p i,n-1 2 i + 1 4 (2 -λ 0 )λ|||t n -z n ||| 2 ≤ λ n k i=1 ω i ε i 1 -ε i p i,n -p i,n-1 2 i + 1 4 (2 -λ n )λ n |||t n -z n ||| 2 → 0, (3.43)
By using now (3.37), we have

k i=1 ω i ε i 1 -ε i p ′ i,n -p i,n-1 2 i + 1 4 (2 -λ 0 ) |||t n -z n ||| 2 1/2 ≤ k i=1 ω i ε i 1 -ε i p i,n -p i,n-1 2 i + 1 4 (2 -λ 0 )|||t n -z n ||| 2 1/2 + ε 1 -ε 1/2 |||a n ||| → 0 (3.44) where ε = max 1≤i≤k ε i . This entails that ν T (t n , P K k p n-1 ) -(t n , P K k p n-1 ) → 0 (3.45)
where ν is the norm given by (3.16). Appealing to Lemma 3.2(ii), if (t n ℓ , P K k p n ℓ -1 ) ℓ∈N is a weakly converging subsequence of (t n , P K k p n-1 ) n∈N , we deduce from the demiclosedness principle [START_REF] Browder | Convergence theorems for sequences of nonlinear operators in Banach spaces[END_REF] that its limit belongs to Fix T . As a result of Lemma 3.1(iii), we conclude that

(t n , P K k p n-1 ) ⇀ ( t, p) ∈ Fix T . (3.46)
By using Lemma 3.2(i), we have thus J γA t ∈ zer(A + B).

Remark 3.4

(i) In Algorithm (3.30), (λ n ) n∈N and (ε 1 , . . . , ε m ) correspond to relaxation and inertia parameters, respectively. The sequences (a i,n ) n∈N with i ∈ {1, . . . , m} and (b n ) n∈N model possible errors in the computation of J γ(1-ε i )A i and J γB .

(ii) When (∀i ∈ {1, . . . , m}) ε i = 0, weaker conditions than (ii) and (iii) are known [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF] to be sufficient to prove the convergence result, namely

(a) n∈N λ n (2 -λ n ) = +∞ (b) n∈N λ n (|||a n ||| + |||b n |||) < +∞.
Stronger convergence results can be deduced from the previous proposition:

Proposition 3.5 Suppose that the assumptions of Proposition 3.3 hold. Let (t n ) n∈N and (p n ) n≥-1 be sequences generated by Algorithm (3.30). Then, (p n ) n∈N and (y n ) n∈N = J γB (2p nt n ) n∈N both converge weakly to y ∈ zer(A + B).

Proof. We proceed similarly to the proof in [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Theorem 25.6] for the convergence of the standard Douglas-Rachford algorithm, which is inspired from [START_REF] Eckstein | A family of projective splitting methods for the sum of two maximal monotone operators[END_REF][START_REF] Svaiter | On weak convergence of the Douglas-Rachford method[END_REF]. By using the same notation as in the proof of Proposition 3.3, it can be deduced from (3.46) that

(∀i ∈ {1, . . . , m}) (1 -ε i )t i,n + ε i p i,n-1 ⇀ (1 -ε i ) t i + ε i p i (3.47) (since i > k ⇒ ε i = 0), which implies that (1 -ε i )t i,n + ε i p i,
y i = p i . (3.50)
Furthermore, we have

(∀n ∈ N) |||p n -y n ||| = 1 2 |||t n -R γB (2p n -t n )||| ≤ 1 2 |||t n -z n ||| + |||z n -R γB (2p n -t n )||| ≤ 1 2 |||t n -z n ||| + |||a n ||| (3.51)
where (3.39) has been used to get the last inequality. In addition, according to (3.45), |||t n -z n ||| → 0 and we know that a n → 0. Hence, (3.51) yields

p n -y n → 0 (3.52) p ′ n -y n → 0. (3.53)
This implies that

y n ℓ ⇀ y. (3.54)
In turn, (3.31) and the relation defining (y n ) n∈N can be rewritten as

(∀n ∈ N) u n ∈ γA(p ′ n ) (3.55) v n ∈ γB(y n ) (3.56)
where On the other hand,

(∀i ∈ {1, . . . , m}) u i,n = t i,n + ε i 1 -ε i p i,n-1 - 1 1 -ε i p ′ i,n (3.57) v n = 2p n -t n -y n . ( 3
(∀n ∈ N)(∀i ∈ {1, . . . , m}) u i,n + v i,n = p i,n -y i,n + ε i 1 -ε i (p i,n-1 -p ′ i,n ) + a i,n . (3
y -t ∈ γB( y). (3.64)
Summing the two inclusion relations leads to y ∈ zer(A + B).

To end the proof, it is sufficient to note that (3.63) is equivalent to y = J γA t. This shows that (p n ) n∈N cannot have a weak cluster point other than J γA t. We have then p n ⇀ y and the weak convergence of (y n ) n∈N follows from (3.52). Under some restrictive assumptions, a linear convergence property can also be proved:

Proposition 3.6 Let γ ∈ ]0, +∞[, (ε 1 , . . . , ε m ) ∈ [0,
(i) zer(A + B) = ∅ (ii) There exists λ ∈]0, 2] such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n ≤ 2.
(iii) For every i ∈ {1, . . . , m},

(ε i > 0 ⇔ i ≤ k). (iv) For every i ∈ {1, . . . , k}, J γ(1-ε i )A i is (1+τ A i )-cocoercive and that R γB is (1+τ B ) -1/2 -strictly contractive, where (τ A 1 , . . . , τ A k , τ B ) ∈ ]0, +∞[ k+1 and
(∀i ∈ {1, . . . , k}) λ > 2 1 - 2τ A i ε i τ B . (3.65)
Then, there exists ρ ∈]0, 1[ such that

(∀n ∈ N) |||t n+1 -t||| 2 + λ n k i=1 ω i ε i 1 -ε i p i,n -p i 2 i ≤ ρ n |||t 1 -t||| 2 + λ 0 k i=1 ω i ε i 1 -ε i p i,0 -p i 2 i
(3.66) where ( t, p) is the unique fixed point of the operator T defined by (3.11)- (3.14).

Proof. Let ( t, p) ∈ Fix T , and set

(∀i ∈ {1, . . . , m}) p ′ i = p i if i ≤ k J γA i t i otherwise (3.67)
and (∀i ∈ {k + 1, . . . , m}) τ A i = 0. Let n be any positive integer. For every i ∈ {1, . . . , m}, as a consequence of the cocoercivity of J γ(1-ε i )A i , we get

(1 + τ A i ) p i,n -p ′ i 2 i ≤ p i,n -p ′ i | (1 -ε i )(t i,n -t i ) + ε i (p i,n-1 -p i ) i ⇒ -(1 -ε i ) p i,n -p ′ i | t i,n -t i ≤ -1 - ε i 2 + τ A i p i,n -p ′ i || 2 i + ε i 2 p i,n-1 -p i 2 i . (3.68)
Furthermore, since R γB has been assumed strictly contractive, we have

(1 + τ B )|||R γB (2p n -t n ) -t||| 2 = (1 + τ B )|||R γB (2p n -t n ) -R γB (2 p ′ -t)||| 2 ≤ |||2(p n -p ′ ) -t n + t||| 2 = 4|||p n -p ′ ||| 2 + |||t n -t||| 2 -4 p n -p ′ | t n -t . (3.69)
Combined with (3.68), this yields

(1 + τ B )|||R γB (2p n -t n ) -t||| 2 ≤ -2 m i=1 ω i ε i + 2τ A i 1 -ε i p i,n -p ′ i 2 i + |||t n -t||| 2 + 2 k i=1 ω i ε i 1 -ε i p i,n-1 -p i 2 i ≤ -2 k i=1 ω i ε i + 2τ A i 1 -ε i p i,n -p i 2 i + |||t n -t||| 2 + 2 k i=1 ω i ε i 1 -ε i p i,n-1 -p i 2 i . (3.70)
By using (3.37), we have then

(1 + τ B )|||t n+1 -t||| 2 ≤ (1 + τ B ) 1 - λ n 2 |||t n -t||| 2 + λ n 2 |||R γB (2p n -t n ) -t||| 2 ≤ 1 + τ B 1 - λ n 2 |||t n -t||| 2 + λ n k i=1 ω i ε i 1 -ε i p i,n-1 -p i 2 i -λ n k i=1 ω i ε i + 2τ A i 1 -ε i p i,n -p i 2 i . (3.71)
By using Assumption (ii), it can be deduced that

1 + min τ B , 2τ A i ε i 1≤i≤k |||t n+1 -t||| 2 + λ n k i=1 ω i ε i 1 -ε i p i,n -p i 2 i ≤ 1 + τ B 1 - λ 2 |||t n -t||| 2 + λ n-1 k i=1 ω i ε i 1 -ε i p i,n-1 -p i 2 i . (3.72)
By setting now

ρ = 1 + τ B 1 -λ 2 1 + min{τ B , (2τ A i /ε i ) 1≤i≤k } (3.73)
(3.66) is obtained by induction. The fact that ρ ∈]0, 1[ follows from (3.65). The uniqueness of ( t, p) is a straightforward consequence of (3.66) and Lemma 3.2(i).

Remark 3.7 If (∀i ∈ {1, . . . , k}) A i = ∂f i where [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]Proposition 2.5], for example). On the other hand, let B = ∂g where

f i ∈ Γ 0 (H i ) is κ f i -strongly convex with κ f i ∈ ]0, +∞[, then J γ(1-ε i )A i = prox γ(1-ε i )f i is (1 + τ A i )-cocoercive with τ A i = γ (1 -ε i )κ f i . (see
(∀x ∈ H) g(x) = h(x) + κ g 2 |||x||| 2 , (3.74) h ∈ Γ 0 (H) has a β-Lipschitz continuous gradient and (β, κ g ) ∈ ]0, +∞[ 2 . Then, it is readily shown that R γB = rprox γg is (1 + τ B ) -1/2 -strictly contractive with τ B = (1 -γκ g ) 2 + 4γκ g (1 + γβ(1 + γκ g )) -1 -1 (1 + γκ g ) 2 -1 (see the appendix).
4 Zero of a sum of an arbitrary number of maximally monotone operators

Let (G, • ) be a real Hilbert space and let • | • be the scalar product of G. Let E be a closed linear subspace of G and let

L : G → H : y → (L 1 y, . . . , L m y) (4.1)
where, for every i ∈ {1, . . . , m},

L i ∈ B (G, H i ) is such that L i (E) is closed. Thus, ∀(x, y) ∈ H × G x, Ly = m i=1 ω i x i , L i y i = m i=1 ω i L * i x i , y . (4.2)
This shows that the adjoint of L is

L * : H → G : (x 1 , . . . , x m ) → m i=1 ω i L * i x i . (4.3)
In this section, the following problem is considered:

find y ∈ zer m i=1 ω i L * i • A i • L i + N E (4.4)
where, for every i ∈ {1, . . . , m}, A i :

H i → 2 H i is a maximally monotone operator. Since (∀y ∈ G) N E (y) = E ⊥ y ∈ E ∅ otherwise, (4.5
Initialization (t i,0 ) 1≤i≤m ∈ H, (p i,-1 ) 1≤i≤m ∈ H y 0 ∈ Arg min z∈E m i=1 ω i L i z -t i,0 2 i For n = 0, 1, . . .            For i = 1, . . . , m p i,n = J γ(1-ε i )A i (1 -ε i )t i,n + ε i p i,n-1 + a i,n c n ∈ Arg min z∈E m i=1 ω i L i z -p i,n 2 i For i = 1, . . . , m t i,n+1 = t i,n + λ n L i (2c n -y n ) -p i,n y n+1 = y n + λ n (c n -y n ). (4.6)
Suppose that the following assumptions hold.

(i) zer m i=1 ω i L * i • A i • L i + N E = ∅ (ii) There exists λ ∈]0, 2[ such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n < 2.
(iii) (∀i ∈ {1, . . . , m}) n∈N a i,n i < +∞.

Then (L 1 y n , . . . , L m y n ) n∈N , (L 1 c n , . . . , L m c n ) n∈N ,
ω i L * i • A i • L i + N E = zer(L * • A • L + N E ) ⇔ y ∈ E and ∃z ∈ A(L y), L * z ∈ E ⊥ ⇔ y ∈ E and ∃z ∈ A(L y), z ∈ L(E) ⊥ ⇔ y ∈ E and L y ∈ zer(A + N D ) (4.7)
where D = L(E). Hence, Problem (4.4) reduces to finding a zero of A + N D . This latter problem is a specialization of Problem (3.1) when B = N D . We have then

J γB = J γ∂ι D = prox ι D = P D .
According to Proposition 3.3, under Assumptions (i)-(iii), the algorithm:

Initialization t 0 ∈ H, p -1 ∈ H For n = 0, 1, . . .          For i = 1, . . . , m p i,n = J γ(1-ε i )A i (1 -ε i )t i,n + ε i p i,n-1 + a i,n d n = P D p n x n = P D t n t n+1 = t n + λ n 2d n -x n -p n (4.8)
allows us to generate a sequence (t n ) n∈N converging weakly to t such that J γA t ∈ zer(A + N D ). As a consequence, J γA t ∈ D. (4.9)

According to the fixed point properties of the Douglas-Rachford algorithm,

t = (2P D -Id )(2J γA t -t). (4.10)
Due to (4.9), (4.10) is equivalent to

P D t = J γA t. (4.11)
By using the weak continuity of P D , it can be deduced that

x n = P D t n ⇀ P D t = J γA t. (4.12)
In other words, (x n ) n∈N converges weakly to L y = J γA t where y ∈ E is a solution to Problem (4.4). In addition, (x n ) n>0 can be computed in a recursive manner through the relation

(∀n ∈ N) x n+1 = P D t n+1 = P D t n + λ n (2d n -x n -P D p n ) = x n + λ n (d n -x n ) (4.13)
where we have used the fact that (x n ) n≥0 and (d n ) n≥0 are sequences of D. Algorithm (4.8) then becomes Initialization

t 0 ∈ H, p -1 ∈ H x 0 = P D t 0 For n = 0, 1, . . .          For i = 1, . . . , m p i,n = J γ(1-ε i )A i (1 -ε i )t i,n + ε i p i,n-1 + a i,n d n = P D p n t n+1 = t n + λ n 2d n -x n -p n x n+1 = x n + λ n (d n -x n ). (4.14)
According to Assumption (ii), (λ n ) n∈N converges to a positive limit, and we deduce from (4.13) that d n ⇀ L y. In turn, the last update equation in Algorithm (4.8) This shows that this problem can also be formulated in terms of the generalized inverse of L w.r.t. E. In particular, if L * L = m i=1 ω i L * i L i is an isomorphism on G, for every u ∈ H, there exists a unique w ∈ G satisfying (4.16) (since z → |||Lz -u||| 2 is strictly convex).

In this case, a stronger convergence result can be obtained as stated below.

Corollary 4.2 Suppose that the assumptions of Proposition 4.1 hold and that m i=1 ω i L * i L i is an isomorphism. Let (y n ) n∈N be a sequence generated by Algorithm (4.6). Then (y n ) n∈N converges weakly to a solution to Problem (4.4).

Proof. According to Proposition 4.1,

Ly n ⇀ L y (4.17

)
where y is a solution to Problem (4.4). By recalling that any bounded linear operator is weakly continuous, we have therefore

L * Ly n ⇀ L * L y (4.18)
and, consequently,

y n = (L * L) -1 L * Ly n ⇀ (L * L) -1 L * L y = y. (4.19)
The assumption that m i=1 ω i L * i L i is an isomorphism may be restrictive. However, a variant of Algorithm (4.6) allows us to relax this requirement. and(λ n ) n∈N be a sequence of reals. For every i ∈ {1, . . . , m}, let (a i,n ) n∈N be a sequence in

Proposition 4.3 Let γ ∈ ]0, +∞[, (ε 1 , . . . , ε m ) ∈ [0, 1[ m , α ∈ ]0, +∞[,
H i . Let {y n } n∈N ⊂ E, {c n } n∈N ⊂ E, {r n } n∈N ⊂ G, {t n } n∈N ⊂ H,
and {p n } n≥-1 ⊂ H be generated by the routine:

Initialization     (t i,0 ) 1≤i≤m ∈ H, (p i,-1 ) 1≤i≤m ∈ H, r 0 ∈ G y 0 = arg min z∈E m i=1 ω i L i z -t i,0 2 i + α z -r 0 2
For n = 0, 1, . . . Proof. Problem (4.4) can be reformulated as

               For i = 1, . . . , m p i,n = J γ(1-ε i )A i (1 -ε i )t i,n + ε i p i,n-1 + a i,n c n = arg min z∈E m i=1 ω i L i z -p i,n 2 i + α z -r n 2 For i = 1, . . . , m t i,n+1 = t i,n + λ n L i (2c n -y n ) -p i,n r n+1 = r n + λ n (2c n -y n -r n ) y n+1 = y n + λ n (c n -y n ).
find y ∈ zer m+1 i=1 ω i L * i • A i • L i + N E
where A m+1 = 0 and L m+1 = Id . (4.21)

In this case, for every

ω m+1 ∈ ]0, +∞[, the self-adjoint operator m+1 i=1 ω i L * i L i = m i=1 ω i L * i L i + ω m+1
Id is an isomorphism. By applying a version of Algorithm (4.6) to the above problem with m + 1 operators, we get:

Initialization     (t i,0 ) 1≤i≤m+1 ∈ H 1 × • • • × H m × G, (p i,-1 ) 1≤i≤m+1 ∈ H 1 × • • • × H m × G y 0 ∈ Arg min z∈E m i=1 ω i L i z -t i,0 2 i + ω m+1 z -t m+1,0 2
For n = 0, 1, . . .

             For i = 1, . . . , m + 1 p i,n = J γ(1-ε i )A i (1 -ε i )t i,n + ε i p i,n-1 + a i,n c n ∈ Arg min z∈E m i=1 ω i L i z -p i,n 2 i + ω m+1 z -p m+1,n 2
For i = 1, . . . , m + 1 t i,n+1 = t i,n + λ n L i (2c n -y n ) -p i,n y n+1 = y n + λ n (c n -y n ). (4.22)
Setting α = ω m+1 , defining (∀n ∈ N) r n = t m+1,n , and setting also ε m+1 = 0 and a m+1,n ≡ 0 yields p m+1,n ≡ t m+1,n , and Algorithm (4.20) is derived.

Remark 4.4 Another variant of Algorithm (4.6) which allows us to take into account errors in the computation of the projection onto D is the following: the same convergence results as for Proposition 4.1 can be proved. Note however that this algorithm requires two projections onto L(E) at each iteration. Algorithm [START_REF] Han | A parallel algorithm for a class of convex programs[END_REF] in [START_REF] Combettes | Iterative construction of the resolvent of a sum of maximal monotone operators[END_REF] is a special case of this algorithm corresponding to the case when

Initialization (t i,0 ) 1≤i≤m ∈ H, (p i,-1 ) 1≤i≤m ∈ H For n = 0, 1, . . .          For i = 1, . . . , m p i,n = J γ(1-ε i )A i (1 -ε i )t i,n + ε i p i,n-1 + a i,n c n -b c n ∈ Arg min z∈E m i=1 ω i L i z -p i,n 2 i y n -b y n ∈ Arg min z∈E m i=1 ω i L i z -t i,n 2 i t n+1 = t n + λ n L(2c n -y n ) -p n .
ε 1 = . . . = ε m = 0, H 1 = . . . = H m = E = G, L 1 = . . . = L m = Id ,

Application to convex optimization

In this section, the following optimization problem is considered:

minimize y∈E m i=1 f i (L i y) (5.1)
where, for every i ∈ {1, . . . , m}, f i ∈ Γ 0 (H i ). As a preliminary result, it can be noticed that:

Proposition 5.1 If there exists y ∈ E such that

(∀i ∈ {1, . . . , m}) L i y ∈ dom f i (5.2)
and lim An algorithm derived from Algorithm (4.6) can be applied to solve Problem (5.1). In the following, (ω i ) 1≤i≤m are positive constants, as in the previous section.

y∈E, L 1 y 1 + ••• + Lmy m→+∞ f 1 (L 1 y) + • • • + f m (L m y) = +∞,
Proposition 5.2 Let (ε 1 , . . . , ε m ) ∈ [0, 1[ m and (λ n ) n∈N be a sequence of reals. For every i ∈ {1, . . . , m}, let (a i,n ) n∈N be a sequence in H i . Let {y n } n∈N ⊂ E, {c n } n∈N ⊂ E, {t n } n∈N ⊂ H and {p n } n≥-1
⊂ H be generated by the routine:

Initialization (t i,0 ) 1≤i≤m ∈ H, (p i,-1 ) 1≤i≤m ∈ H y 0 ∈ Arg min z∈E m i=1 ω i L i z -t i,0 2 i For n = 0, 1, . . .              For i = 1, . . . , m p i,n = prox (1-ε i )f i ω i (1 -ε i )t i,n + ε i p i,n-1 + a i,n c n ∈ Arg min z∈E m i=1 ω i L i z -p i,n 2 i For i = 1, . . . , m t i,n+1 = t i,n + λ n L i (2c n -y n ) -p i,n y n+1 = y n + λ n (c n -y n ).
(5.5)

Suppose that the following assumptions hold.

(i) 0 ∈ sri (L 1 z -x 1 , . . . , L m z -x m ) z ∈ E, x 1 ∈ dom f 1 , . . . , x m ∈ dom f m . (ii) There exists λ ∈]0, 2[ such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n < 2.
(iii) (∀i ∈ {1, . . . , m}) n∈N a i,n i < +∞.

If the set of solutions to Problem (5.

y ∈ zer ∂(f • L + ι E ). (5.7) Provided that 0 ∈ sri Lz -x z ∈ E, x ∈ dom f (5.8)
we have [START_REF] Zȃlinescu | Convex Analysis in General Vector Spaces[END_REF]Theorem 2.8.3]

∂(f • L + ι E ) = L * • ∂f • L + ∂ι E .
(5.9)

Due Assumption (i), (5.8) is obviously satisfied and (5.7) is therefore equivalent to

y ∈ zer(L * • ∂f • L + ∂ι E ) = zer m i=1 L * i • ∂f i • L i + N E (5.10)
where we have used the fact that ∂f = (ω -1

1 ∂f 1 ) × • • • × (ω -1 m ∂f m )
, which follows from the definition of the norm in (2.1). So, Problem 5.1 appears as a specialisation of Problem 4.4 when (∀i ∈ {1, . . . , m}) A i = ∂f i /ω i . Algorithm (5.5) is then derived from Algorithm (4.6) (with γ = 1) since

J (1-ε i )A i = prox (1-ε i )f i /ω i
and its convergence follows from Proposition 4.1.

When m i=1 ω i L * i L i is an isomorphism on G, a stronger convergence result can be deduced from Corollary 4.2.

Corollary 5.3 Suppose that the assumptions of Proposition 5.2 hold and that m

i=1 ω i L * i L i is an isomorphism. Let (y n ) n∈N be a sequence generated by Algorithm (5.5). If the set of solutions to Problem (5.1) is nonempty, then (y n ) n∈N converges weakly to an element of this set.

A variant of Algorithm (5.5) is applicable when m i=1 ω i L * i L i is not an isomorphism.

Proposition 5.4 Let (ε 1 , . . . , ε m ) ∈ [0, 1[ m , α ∈ ]0,
Initialization     (t i,0 ) 1≤i≤m ∈ H, (p i,-1 ) 1≤i≤m ∈ H, r 0 ∈ G y 0 = arg min z∈E m i=1 ω i L i z -t i,0 2 i + α z -r 0 2 For n = 0, 1, . . .                For i = 1, . . . , m p i,n = prox 1-ε i ω i f i (1 -ε i )t i,n + ε i p i,n-1 + a i,n c n = arg min z∈E m i=1 ω i L i z -p i,n 2 i + α z -r n 2 For i = 1, . . . , m t i,n+1 = t i,n + λ n L i (2c n -y n ) -p i,n r n+1 = r n + λ n (2c n -y n -r n ) y n+1 = y n + λ n (c n -y n ).
(5.11)

Suppose that Assumptions (i)-(iii) in Proposition 5.2 hold. If the set of solutions to Problem (5.1) is nonempty, then (y n ) n∈N converges weakly to an element of this set.

Proof.

By proceeding similarly to the proof of Proposition 5.2, it can be deduced from Proposition 4.3 that the sequence (y n ) n∈N generated by Algorithm (5.11) converges to a solution to Problem (5.1) when such a solution exists. (5.12) Indeed, in this case, we have (see Proposition 3.6 in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]):

0 ∈ sri{D -dom f } = ri{D -dom f } = ri D -ri dom f = D -ri dom f ⇔ D ∩ ri dom f = ∅. ( 5

Connections with some other parallel splitting optimization algorithms

Firstly, is interesting to note that PPXA (Algorithm 3.1 in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]) is a special case of Algorithm (5.5) corresponding to the case when

ε 1 = . . . = ε m = 0, H 1 = . . . = H m = E = G, and L 1 = . . . = L m = Id .
In order to better emphasize the link existing with other algorithms, it will appear useful to rewrite Algorithm (5.5) in a different form. Proposition 6.1 Let (ε 1 , . . . , ε m ) ∈ [0, 1[ m and (λ n ) n≥-2 be a sequence of reals. Algorithm (5.5) with (p i,-1 ) 1≤i≤m = (L 1 y 0 , . . . , L m y 0 ) is equivalent to the following routine:

Initialization             λ -2 = λ -1 = 1 (e i,0 ) 1≤i≤m ∈ H u 0 ∈ Argmin z∈E m i=1 ω i L i z + e i,0 2 i For i = 1, . . . , m     ℓ i,0 = 0 w i,-1 = 0 w i,0 = L i u i,0 For n = 0, 1, . . .                 u n+1 ∈ Argmin z∈E m i=1 ω i L i z -λ n-1 w i,n + e i,n 2 i For i = 1, . . . , m v i,n+1 = L i u n+1 w i,n+1 = prox 1-ε i ω i f i (1 -ε i )(v i,n+1 + e i,n -ℓ i,n ) + ε i w i,n + a i,n k n+1 ∈ Argmin z∈E m i=1 ω i L i z -e i,n + λ n-2 w i,n-1 2 i For i = 1, . . . , m ℓ i,n+1 = (1 -λ n )L i k n+1 e i,n+1 = v i,n+1 + e i,n -λ n w i,n+1 + ℓ i,n+1 . (6.1)
Suppose that Assumptions (i), (ii) and (iii) in Proposition 5.2 hold. If the set of solutions to Problem (5.1) is nonempty, then (v 1,n , . . . , v m,n ) n>0 converges weakly to (L 1 y, . . . , L m y) where y is a solution to Problem (5.1).

Proof. Let us start from Algorithm (4.8) and express it as follows: for every n ∈ N,

             (∀i ∈ {1, . . . , m}) p D i,n + p ⊥ i,n = prox 1-ε i ω i f i (1 -ε i )(t D i,n + t ⊥ i,n ) + ε i (p D i,n-1 + p ⊥ i,n-1 ) + a i,n d n = p D n x n = t D n t D n+1 = t D n + λ n 2d n -x n -p D n t ⊥ n+1 = t ⊥ n -λ n p ⊥ n (6.
2) where p D n (resp p ⊥ n ) denotes the projection of p n onto D (resp. D ⊥ ), and p D i,n (resp. p ⊥ i,n ) is its i-th component in H i , a similar notation being used for the other variables. The above set of equations can be rewritten as

     (∀i ∈ {1, . . . , m}) p D i,n + p ⊥ i,n = prox 1-ε i ω i f i (1 -ε i )(t D i,n + t ⊥ i,n ) + ε i (p D i,n-1 + p ⊥ i,n-1 ) + a i,n t D n+1 = t D n + λ n p D n -t D n t ⊥ n+1 = t ⊥ n -λ n p ⊥ n . (6.
3) Let us now introduce sequences (w n ) n∈N and (e n ) n∈N such that

w n = p n-1 (6.4) e ⊥ n = t ⊥ n . (6.5)
Then, for every n ∈ N, (6.3) becomes

             (∀i ∈ {1, . . . , m}) w D i,n+1 + w ⊥ i,n+1 = prox 1-ε i ω i f i (1 -ε i ) (1 -λ n-1 )t D i,n-1 + λ n-1 w D i,n + e ⊥ i,n + ε i (w D i,n + w ⊥ i,n ) +a i,n t D n+1 = (1 -λ n )t D n + λ n w D n+1 e ⊥ n+1 = e ⊥ n -λ n w ⊥ n+1 (6.6) provided that λ -1 = 1 (thus allowing us to choose t D i,-1 arbitrarily) since p -1 = (L 1 y 0 , . . . , L m y 0 ) = x 0 ∈ D ⇔ w D 0 = p -1 = x 0 = t D 0 .
In addition, for every n ∈ N, set

e D n+1 = λ n-1 p D n-1 -t D n+1 = λ n-1 w D n -t D n+1 . (6.7)
Then, we have: for every n ∈ N,

t D n+1 = (1 -λ n )t D n + λ n w D n+1 ⇔ e D n+1 = (1 -λ n )(e D n -λ n-2 w D n-1 ) + λ n-1 w D n -λ n w D n+1 (6.8)
provided that λ -2 = 1 and w D -1e D 0 = w D 0 . By using the two previous relations, we see that (6.6) is equivalent to

               (∀i ∈ {1, . . . , m}) w D i,n+1 + w ⊥ i,n+1 = prox 1-ε i ω i f i (1 -ε i ) (1 -λ n-1 )(λ n-3 w D i,n-2 -e D i,n-1 ) + λ n-1 w D i,n + e ⊥ i,n +ε i (w D i,n + w ⊥ i,n ) + a i,n e D n+1 = (1 -λ n )(e D n -λ n-2 w D n-1 ) + λ n-1 w D n -λ n w D n+1 e ⊥ n+1 = e ⊥ n -λ n w ⊥ n+1 .
(6.9) By introducing intermediate variables v n+1 and ℓ n+1 , this can be rewritten as

                 v n+1 = λ n-1 w D n -e D n (∀i ∈ {1, . . . , m}) w D i,n+1 + w ⊥ i,n+1 = prox 1-ε i ω i f i (1 -ε i )(v i,n+1 + e D i,n + e ⊥ i,n -ℓ i,n ) + ε i (w D i,n + w ⊥ i,n ) + a i,n ℓ n+1 = (1 -λ n )(e D n -λ n-2 w D n-1 ) e D n+1 = v n+1 + e D n -λ n w D n+1 + ℓ n+1 e ⊥ n+1 = e ⊥ n -λ n w ⊥ n+1
(6.10) where an appropriate initialization is

λ -2 = λ -1 = 1 (6.11)
ℓ 0 = 0 (6.12)

w D -1 = 0, w D 0 = -e D 0 . (6.13)
Eqs (6.10)-(6.13) obviously yield the following algorithm:

Initialization       λ -2 = λ -1 = 1 ℓ 0 = 0 e 0 ∈ H w -1 = 0, w 0 = -P D e 0 For n = 0, 1, . . .           v n+1 = P D (λ n-1 w n -e n ) ℓ n+1 = (1 -λ n )P D (e n -λ n-2 w n-1 ) For i = 1, . . . , m w i,n+1 = prox 1-ε i ω i f i (1 -ε i )(v i,n+1 + e i,n -ℓ i,n ) + ε i w i,n + a i,n e n+1 = v n+1 + e n -λ n w n+1 + ℓ n+1 . (6.14)
By setting now w 0 = Lu 0 and, for every n ∈ N, v n+1 = Lu n+1 and ℓ n+1 = (1 -λ n )Lk n+1 where (u n ) n≥0 and (k n ) n>0 are sequences of E, Algorithm (6.1) is obtained.

Since the assumptions of Proposition 5.2 hold, x n = t D n ⇀ L y where L y is a minimizer of f over D. It can be deduced from (6.8) that λ n (w D n+1t D n ) ⇀ 0. In addition, from the definition of v n+1 and (6.7), for every n ∈ N * ,

v n+1 = λ n-1 w D n -λ n-2 w D n-1 + t D n = λ n-1 (w D n -t D n-1 ) -λ n-2 (w D n-1 -t D n-2 ) + λ n-1 t D n-1 -λ n-2 t D n-2 + t D n . (6.15)
Due to Assumption (ii) in Proposition 5.2, λ n-1 t D n-1 -λ n-2 t D n-2 ⇀ 0, which allows us to conclude that v n+1 ⇀ L y.

Similarly to Corollary 5.3, we have also: Corollary 6.2 Suppose that the assumptions of Proposition 5.2 hold and that m i=1 ω i L * i L i is an isomorphism. Let (u n ) n∈N be a sequence generated by Algorithm (6.1). If the set of solutions to Problem (5.1) is nonempty, then (u n ) n∈N converges weakly to an element of this set.

Note that Algorithm (6.1) may appear somewhat less efficient than Algorithm (5.5) in the sense that it requires to compute two projections on D at each iteration. This disadvantage no longer exists in the unrelaxed case (λ n ≡ 1) where Algorithm (6.1) takes a simplified form. Proposition 6.3 Let (ε 1 , . . . , ε m ) ∈ [0, 1[ m . Algorithm (6.1) with λ n ≡ 1 is equivalent to the following routine:

Initialization       (e i,0 ) 1≤i≤m ∈ H u 0 = argmin z∈E m i=1 ω i L i z + e i,0 2 i For i = 1, . . . , m w i,0 = L i u 0 For n = 0, 1, . . .          u n+1 = argmin z∈E m i=1 ω i L i z -w i,n + e i,n 2
i For i = 1, . . . , m      v i,n+1 = L i u n+1 w i,n+1 = prox 1-ε i ω i f i (1 -ε i )(v i,n+1 + e i,n ) + ε i w i,n + a i,n
e i,n+1 = v i,n+1 + e i,n -w i,n+1 .

(6.16)

Suppose that the following assumptions hold.

(i) (0, . . . , 0) ∈ sri (L 1 z -x 1 , . . . , L m z -x m ) z ∈ E, x 1 ∈ dom f 1 , . . . , x m ∈ dom f m . (ii) (∀i ∈ {1, . . . , m}) n∈N a i,n i < +∞. (iii) m i=1 ω i L * i L i is an isomorphism.
If the set of solutions to Problem (5.1) is nonempty, then (u n ) n∈N converges weakly to an element of this set.

The above algorithm when E = G, ε 1 = . . . = ε m = 0, ω 1 = . . . = ω m and (a i,n ) 1≤i≤m ≡ (0, . . . , 0) was derived in [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF] from the Douglas-Rachford algorithm by invoking a duality argument. We have here obtained this algorithm by following a different way. In addition, the convergence result stated in Proposition 6.3 is not restricted to finite dimensional Hilbert spaces and it allows us to take into account an error term (a i,n ) 1≤i≤m n∈N in the computation of the proximity operators of the functions (f i ) 1≤i≤m and to include inertia parameters. It can also be noticed that Algorithm (5.5) may appear more flexible than Algorithm (6.16) since it offers the ability of using relaxation parameters (λ n ) n∈N .

Similarly to the derivation of Algorithm (5.11), when Assumption (iii) of Proposition 6.3 is not satisfied a variant of Algorithm (6.16)

Initialization        (e i,0 ) 1≤i≤m ∈ H u 0 = argmin z∈E m i=1 ω i L i z + e i,0 2 i + α z 2 For i = 1, . . . , m w i,0 = L i u 0 For n = 0, 1, . . .           u n+1 = argmin z∈E m i=1 ω i L i z -w i,n + e i,n 2 i + α z -u n 2 For i = 1, . . . , m      v i,n+1 = L i u n+1 w i,n+1 = prox 1-ε i ω i f i (1 -ε i )(v i,n+1 + e i,n ) + ε i w i,n + a i,n e i,n+1 = v i,n+1 + e i,
f i (L i y) (6.18)
where f m+1 = 0 and L m+1 = Id . Hence, for every

ω m+1 ∈ ]0, +∞[, m+1 i=1 ω i L * i L i = m i=1 ω i L * i L i + ω m+1
Id is an isomorphism. By applying a version of Algorithm (6.16) to the above problem with m + 1 potentials, we get:

Initialization       (e i,0 ) 1≤i≤m+1 ∈ H 1 × • • • × H m × G u 0 = argmin z∈E m i=1 ω i L i z + e i,0 2 i + ω m+1 z + e m+1,0 2 For i = 1, . . . , m + 1 w i,0 = L i u 0 For n = 0, 1, . . .          u n+1 = argmin z∈E m i=1 ω i L i z -w i,n + e i,n 2 i + ω m+1 z -w m+1,n + e m+1,n 2
For i = 1, . . . , m + 1      v i,n+1 = L i u n+1 w i,n+1 = prox 1-ε i ω i f i (1 -ε i )(v i,n+1 + e i,n ) + ε i w i,n + a i,n e i,n+1 = v i,n+1 + e i,n -w i,n+1 . (6.19)
Then, by setting e m+1,0 = 0, ω m+1 = α, ε m+1 = 0, a m+1,n ≡ 0, and noticing that (∀n ∈ N) e m+1,n+1 = 0, v m+1,n+1 = u n+1 and w m+1,n = u n , Algorithm (6.17) is obtained. In addition, Assumption (i) in Proposition 6.3 is equivalent to assume that

V = λ>0 λ(L 1 z -x 1 , . . . , L m z -x m ) z ∈ E, x 1 ∈ dom f 1 , . . . , x m ∈ dom f m is a closed vector subspace of H. V × G is therefore a closed vector subspace of H × G. But, since dom f m+1 = G, we have V × G = λ>0 {λ(L 1 z -x 1 , . . . , L m z -x m , z -x m+1 ) | z ∈ E, x 1 ∈ dom f 1 , . . . , x m ∈ dom f m , x m+1 ∈ dom f m+1 }.
This shows that

0 ∈ sri (L 1 z -x 1 , . . . , L m z -x m , z -x m+1 ) z ∈ E, x 1 ∈ dom f 1 , . . . , x m+1 ∈ dom f m+1 . ( 6
.20) Since this condition is satisfied and Assumptions (ii) and (iii) hold, the convergence of (u n ) n∈N follows from Proposition 6.3 applied to Problem (6.18). Remark 6.5 (i) In the case when ε 1 = . . . = ε m = 0, Algorithm (6.16) can be derived in a more direct way from Spingarn's method of partial inverses [START_REF] Spingarn | Partial inverse of a monotone operator[END_REF][START_REF] Spingarn | Applications of the method of partial inverses to convex programming: decomposition[END_REF], which is recalled below:

Initialization (s 0 , q 0 ) ∈ D × D ⊥ For n = 0, 1, . . .     (w ′ n+1 , e ′ n+1 ) such that w ′ n+1 + e ′ n+1 = s n + q n and e ′ n+1 ∈ ∂f (w ′ n+1 ) s n+1 = P D w ′ n+1 q n+1 = P D ⊥ e ′ n+1 (6.21)
Indeed at each iteration n ∈ N, w ′ n+1 and e ′ n+1 are then computed as

w ′ n+1 = prox f (s n + q n ) (6.22) e ′ n+1 = s n + q n -w ′ n+1 . (6.23)
By setting w ′ 0 = s 0 and e ′ 0 = q 0s 0 , and by defining, for every n ∈ N,

v ′ n+1 = s n -P D e ′ n = P D (w ′ n -e ′ n ) (6.24)
we have

s n + q n = v ′ n+1 + P D e ′ n + P D ⊥ e ′ n = v ′ n+1 + e ′ n . (6.25)
Altogether, (6.22)-(6.25) yield the following algorithm:

Initialization e ′ 0 ∈ H w ′ 0 = -P D e ′ 0 For n = 0, 1, . . .     v ′ n+1 = P D (w ′ n -e ′ n ) w ′ n+1 = prox f (v ′ n+1 + e ′ n ) e ′ n+1 = v ′ n+1 + e ′ n -w ′ n+1 . (6.26) By setting (w ′ n ) n∈N = (w n ) n∈N , (e ′ n ) n∈N = (e n ) n∈N , and (v ′ n ) n>0 = (v n ) n>0
, the above algorithm leads to Algorithm (6.16) when ε 1 = . . . = ε m = 0, in the absence of error terms.

(ii) When ε 1 = . . . = ε m = α/(1 + α) and a n ≡ 0, the iterations of Algorithm (6.17)

w n+1 = argmin x∈H 1 1 + α f (x) + 1 2 |||x - 1 1 + α (Lu n+1 + e n ) - α 1 + α w n ||| 2 = argmin x∈H f (x) + 1 2 |||x -Lu n+1 -e n ||| 2 + α 2 |||x -w n ||| 2 (6.31) e n+1 = Lu n+1 + e n -w n+1 . (6.32)
By setting γ = α -1/2 and (∀n > 0) e ′ n = γe n , these iterations read also

(∀n ∈ N) u n+1 = argmin z∈G ι E (z) + e ′ n | Lz + γ 2 |||Lz -w n ||| 2 + 1 2γ z -u n 2 (6.33)
w n+1 = argmin x∈H γf (x) -e ′ n | x + γ 2 |||x -Lu n+1 ||| 2 + 1 2γ |||x -w n ||| 2 (6.34)
e ′ n+1 = e ′ n + γ(Lu n+1w n+1 ). (6.35) Therefore, this algorithm takes the form of the one proposed in [START_REF] Attouch | Augmented Lagrangian and proximal alternating direction methods of multipliers in Hilbert spaces -applications to games, PDEs and control[END_REF] when applied to the optimization problem minimize y∈G ι E (y) + γf (Ly). (6.36) In other words, when ε 1 = . . . = ε m = α/(1 + α) and a n ≡ 0, Algorithm (6.17) can be viewed as a parallelized version of the algorithm proposed in [START_REF] Attouch | Augmented Lagrangian and proximal alternating direction methods of multipliers in Hilbert spaces -applications to games, PDEs and control[END_REF].

Numerical experiments

In order to evaluate the performance of the proposed algorithm, we present numerical experiments relative to an image restoration problem. The considered degradation model is z = Ay + ǫ where y ∈ R N denotes the image to be recovered, z ∈ R N are the observed data, A ∈ R N ×N corresponds to a circulant-block circulant matrix associated with a blur operator, and ǫ ∈ R N is a realization of a random vector with independent and identically distributed components following a uniform distribution. The restoration process is based on a variational approach aiming at finding an image y ∈ R N as close as possible to y ∈ R N from the observation z, and prior informations such as the dynamic range of the pixel values and the sparsity of y in a wavelet-frame representation. The considered optimization problem is :

find y ∈ Argmin y∈R N z -Ay p ℓp + ϑ F y ℓ 1 + ι C (y) (7.1)
with p ≥ 1 and ϑ > 0. F ∈ R K×N denotes the matrix associated with the frame analysis operator.

The frame is supposed to be tight, i.e. F * F = ν F Id , for some ν F > 0. C = [0, 255] N models a dynamic range constraint. The notation • ℓp stands for the classical ℓ p norm:

(∀y = (η i ) 1≤i≤N ∈ R N ) y ℓp = N i=1 |η i | p 1/p . ( 7
m = 3, E = R N , f 1 = z-• p ℓp , f 2 = • ℓ 1 , f 3 = ι C , L 1 = A, L 2 = F , and L 3 = Id .
In the following, this approach is used to restore a satellite image of size N = 256×256 corrupted by a Gaussian blur with a large kernel size and a noise uniformly distributed on [-30, 55]. The latter values are not assumed to be known in the recovery process. In order to deal efficiently with the uniformly distributed noise, the data fidelity measure f 1 is chosen to be an ℓ 3 -norm (p = 3). F corresponds to a tight frame version of the dual-tree transform (DTT) proposed in [START_REF] Chaux | Image analysis using a dual-tree M -band wavelet transform[END_REF] (ν F = 2) using symmlets of length 6 over 2 resolution levels.

Algorithm (5.5) is then employed to solve (7.1). The proximity operator of f 1 here takes a closed form expression [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF], whereas the proximity operator of f 2 reduces to a soft-thresholding rule (see [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF] and references therein). In order to efficiently compute c n at iteration n of (5.5), fast Fourier diagonalization techniques have been employed. It can be noticed that the qualification condition (i) in Proposition 5.2 is satisfied since dom f 1 = dom f 2 = R N and int C = ]0, 255[ N = ∅ (see Remark 5.5(i)). The convergence to an optimal solution to Problem 7.1 is then guaranteed by Corollary 5.3.

The original, degraded, and restored images are presented in Figure 1. The value of ϑ has been selected so as to maximize the signal-to-noise-ratio (SNR). The SNR between an image y and the original image y is defined as 20 log 10 ( y ℓ 2 / y -y ℓ 2 ). In our experiments, we also compare the images in terms of structural similarity (SSIM) [START_REF] Wang | Mean squared error: love it or leave it?[END_REF]. SSIM can be seen as a percentage of similarity, the value 1 being achieved for two identical images.

Note that for such an application, the Douglas-Rachford algorithm would not be a proper choice since there is no closed form expression for the proximity operator of f 1 • A. Regarding PPXA [21, Algorithm 3.1], it would be possible to split the data fidelity term as proposed in [START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF], but this solution would not be as efficient as the proposed one due to the large kernel size of the blur.

The effects of the parameters (λ n ) n∈N and (ε i ) 1≤i≤m are evaluated in Figures 2(a) and (b) where the criterion variations ( z -Ay n p ℓp + ϑ F y n ℓ 1 ) n are plotted as a function of the iteration number n. On the one hand, it results from our experiments that for a fixed relaxation parameter λ n ≡ λ, a faster convergence is achieved by choosing λ close to 2. Recall that, according to Proposition 6.3, there exist close connections between the proposed algorithm and SDMM when λ = 1. On the other hand, an identical value ε = 0.4 of the inertial parameters (ε i ) 1≤i≤m also appears to be beneficial to the convergence profile. In particular, the introduction of inertial parameters allows us to reduce the oscillations in the variations of the criterion. Proof. From (7.3) and standard properties of the proximity operator [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]

  and (p n ) n∈N converge weakly to (L 1 y, . . . , L m y) where y is a solution to Problem (4.4). Proof. By using (4.5) and the definition of A in (3.2), we have the following equivalences y ∈ zer m i=1

(4. 20 )

 20 Suppose that Assumptions (i)-(iii) in Proposition 4.1 hold. Then, (y n ) n∈N converges weakly to a solution to Problem (4.4).

(4. 23 )

 23 Under the assumptions of Proposition 4.1 and provided that the error sequences (b y n ) n∈N and (b c n ) n∈N are such that +∞ n=0 b y n + b c n < +∞ (4.24)

and b y n ≡ b c n ≡ 0.

  then the set of solutions to Problem (5.1) is nonempty. Proof. Problem (5.1) is equivalent to minimize x∈D f (x) where D = L(E). (5.3) As a consequence of classical results of convex analysis [28, Proposition II.1.2], the set of solutions to Problem (5.3) is nonempty if dom f ∩ D = dom f ∩ L(E) = ∅ and lim x∈D, |||x|||→+∞ f (x) = +∞ (5.4) which yields the desired result.

Remark 5. 5 (

 5 i) When H is finite dimensional, Assumptions (i) in Proposition 5.2 takes the form∃y ∈ E, L 1 y ∈ ri dom f 1 , . . . , L m y ∈ ri dom f m .

  .13) (ii) We have seen that Problem (5.1) can be put under the form of Problem (5.3). The dual formulation of this problem is minimize x∈H ι * D (-x) + f * (x). (5.14) Recall that ι * D = ι D ⊥ . The dual problem can thus be rewritten as minimize x∈D ⊥ f * (x). (5.15) Algorithm (4.8) can again be used to solve the dual problem. However, since P D ⊥ = Id -P D and, (∀γ > 0) (∀i ∈ {1, . . . , m}) prox γf * i = Id -γ prox f i /γ (•/γ) [23, Lemma 2.10], the resulting algorithm for the dual problem takes a form very similar to the algorithm proposed for the primal one.

. 2 )

 2 Problem (7.1) is a particular case of Problem (5.1) when

Figure 1 :Figure 2 :

 12 Figure 1: Image restoration result.

  If every weak sequential cluster point of (t n , p n ) n∈N belongs to S, then (t n , p n ) n∈N converges weakly to a point in S.

	Proof.

. , k}, (µ i,n ) n∈N is a sequence of nonnegative reals converging to a positive limit, and (∆ n ) n∈N is a summable sequence of nonnegative reals. Then, (i) (t n ) n∈N and (p n ) n∈N are bounded sequences. (ii) For every ( t, p ∈ S, (|||t n -t||| 2 + k i=1 µ i,n p i,n -p i 2 i ) n∈N converges. (iii) (i) Let ( t, p ∈ S. According to [16, Lemma 3.1(i)], there exists ρ ∈ ]0, +∞[ such that

  On the other hand, t satisfies (3.23) if and only if it is a fixed point of the standard Douglas-Rachford iteration, that is if and only if J γA t is a zero of A + B[START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF][START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF].

			.20)
	So, by using (3.10), (3.19) is equivalent to
	p ′ = J γA t.	(3.21)
	and ( t, p) ∈ Fix T if and only if
	(∀i ∈ {1, . . . , m}) p i =	J γA i t i if i ≤ k 0 otherwise	(3.22)
	t = R γB 2 J γA t -t .		(3.23)

  1[ m and (λ n ) n∈N be a sequence of reals. Let (a n ) n∈N and (b n ) n∈N be sequences in H, and let (t n ) n∈N and (p n ) n≥-1 be sequences generated by the following routine.

  n-1 n∈N is a bounded sequence. For every i ∈ {1, . . . , m}, as a consequence of the nonexpansivity of J γ(1-ε i )A i , (p ′ n ) n∈N as defined by (3.31) is bounded. Hence, there exists a subsequence (p ′ n ℓ ) ℓ∈N of (p ′ n ) n∈N and y ∈ H such that

	p ′ n ℓ ⇀ y.	(3.48)
	Due to (3.37) and Assumption (iii) of Proposition 3.3, we have also
	p n ℓ ⇀ y	(3.49)

and, by using (3.46), (∀i ∈ {1, . . . , k})

  Let γ ∈ ]0, +∞[, (ε 1 , . . . , ε m ) ∈ [0, 1[ m , and (λ n ) n∈N be a sequence of reals. For every i ∈ {1, . . . , m}, let (a i,n ) n∈N be a sequence in H i . Let {y n } n∈N ⊂ E, {c n } n∈N ⊂ E, {t n } n∈N ⊂ H and {p n } n≥-1 ⊂ H be generated by the routine:

	Proposition 4.1
	)
	Problem (4.4) is equivalent to Problem (1.2).
	An algorithm derived from Algorithm (3.30) can be applied to solve Problem (4.4).

  can be employed to solve Problem (5.1). Let (ε 1 , . . . , ε m ) ∈ [0, 1[ m and α ∈ ]0, +∞[. For every i ∈ {1, . . . , m}, let (a i,n ) n∈N be a sequence in H i . Let {u n } n∈N ⊂ E, {v n } n∈N ⊂ H {w n } n∈N ⊂ H and {e n } n∈N ⊂ H be generated by the routine:

	Proposition 6.4

  n -w i,n+1 .Suppose that Assumptions (i) and (ii) of Proposition 6.3 hold. If the set of solutions to Problem (5.1) is nonempty, then (u n ) n∈N converges weakly to an element of this set.

		(6.17)
	Proof. Problem (5.1) can be reformulated as
		m+1
	minimize y∈E	i=1

  can be re-expressed as(∀n ∈ N) u n+1 = argmin z∈E |||Lzw n + e n ||| 2 + α z -u n N) u n+1 = argmin z∈G ι E (z) + |||Lzw n + e n ||| 2 + α z -u n

					2	(6.27)
	w n+1 = prox f 1+α	1 1 + α	(Lu n+1 + e n ) +	α 1 + α	w n	(6.28)
	e n+1 = Lu n+1 + e n -w n+1			(6.29)
	which is equivalent to
	(∀n ∈ 2		(6.30)

  ∈ H 2 ) ||| rprox g x -rprox g y||| 2 = 1 (1 + κ) 2 (1 -κ) 2 |||x -y||| 2 -4(1 -κ) xy | prox (1+κ)h * x -prox (1+κ)h * y + 4 ||| prox (1+κ)h * x -prox (1+κ)h * y||| 2 ≤ 1 (1 + κ) 2 (1 -κ) 2 |||x -y||| 2 + 4κ xy | prox (1+κ)h * x -prox (1+κ)h * y (7.5)where the fact that the proximity operator is firmly nonexpansive has been used in the last inequality. On the other hand, according to [6, Theorem 2.1], we know that h has a 1/β-Lipschitz continuous gradient if and only if h * is strongly convex with modulus β, that is ∈ H 2 ) ||| rprox g x -rprox g y||| 2

							, we get
		(∀x ∈ H)	prox g x = prox h 1+κ	x 1 + κ
							=	1 1 + κ	x -prox (1+κ)h * x
			⇔ rprox g x =	1 1 + κ	(1 -κ)x -2 prox (1+κ)h * x .	(7.4)
	Then, we have
	(∀(x, y) h * = ϕ +	β 2	||| • ||| 2 .	(7.6)
	where ϕ ∈ Γ 0 (H). Thus,
		(∀x ∈ H)	prox (1+κ)h * x = prox (1+κ)ϕ 1+β(1+κ)	x 1 + β(1 + κ)	.	(7.7)
	By invoking now the Cauchy-Schwarz inequality and the nonexpansivity of prox (1+κ)ϕ
							1+β(1+κ)
	≤ |||x -y||| prox (1+κ)ϕ 1+β(1+κ)	x 1 + β(1 + κ)	-prox (1+κ)ϕ 1+β(1+κ)	y 1 + β(1 + κ)
	≤	1 1 + β(1 + κ)	|||x -y||| 2 .	(7.8)
	By combining (7.5) and (7.8), we finally obtain
		(∀(x, y) ≤	1 (1 + κ) 2 (1 -κ) 2 +	4κ 1 + β(1 + κ)	|||x -y||| 2 .	(7.9)
	Since	1 (1 + κ) 2 (1 -κ) 2 +	4κ 1 + β(1 + κ)	<	1 (1 + κ) 2 (1 -κ) 2 + 4κ = 1	(7.10)
	rprox

, we deduce that xy | prox (1+κ)h * x -prox (1+κ)h * y g is thus strictly contractive.

Appendix

Proposition 7.1 Let g be the function defined as

where h ∈ Γ 0 (H) has a β-Lipschitz continuous gradient and (β,