Arnaud Lenoir

Philippe Mahey
email: philippe.mahey@isima.fr

A survey on operator splitting and decomposition of convex programs

Keywords: Operator Splitting, Augmented Lagrangian, Decomposition Methods. Mathematics Subject Classification : 65K13, 90C25, 90C30

Many structured convex minimization problems can be modeled by the search of a zero of the sum of two monotone operators. Operator splitting methods have been designed to decompose and regularize at the same time these kind of models. We review here these models and the classical splitting methods. We focus on the numerical sensitivity of these algorithms with respect to the scaling parameters that drive the regularizing terms, in order to accelerate convergence rates for different classes of models.

Introduction

We survey here some classical monotone operator splitting methods and discuss technical issues which address the difficult question of the acceleration of the numerical performance of these techniques when dealing with the decomposition of large-scale convex programs.

Monotone operator theory has been developed for Hilbert spaces in the seventies by different people, where we first distinguish the monograph by Brezis [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] and the seminal research results of J.J. Moreau [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] on the proximal mapping. The Proximal Point Method (PPM) became popular in the Mathematical Programming community with the seminal papers of Rockafellar [START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF][START_REF] Rockafellar | Augmented lagrangians and applications of the proximal point algorithm in convex programming[END_REF] who showed precisely its link with the Augmented Lagrangian algorithm. Operator splitting is generally referred when dealing with the sum of maximal monotone operators and aiming at decomposing the numerical computations on each operator separately. Early splitting algorithms were analyzed by Lieutaud in his thesis [START_REF] Lieutaud | Approximation d'opérateurs par les méthodes de décomposition[END_REF] using the term Fractional Steps method, borrowed from early works by russian researchers (Denidov, Marchuk, Samarskii and Yanenko, see Temam [START_REF] Temam | Sur la stabilité et la convergence de la méthode des pas fractionnaires[END_REF] for a convergence analysis). The first application of the celebrated Alternate Direction Method of Multipliers (ADMM) stems back to the work of Glowinski and Marocco [START_REF] Glowinski | Sur l'approximation par éléments finis d'ordre 1 et la résolution par pénalisation-dualité d'une classe de problèmes de dirichlet[END_REF] who solved some heat conduction equations, see too Gabay and Mercier [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite-element approximations[END_REF] and the theoretical analysis by Gabay [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF]. In 1979, Lions and Mercier [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF] analyzed the convergence of a family of splitting methods (combining forward and backward steps like in the earlier methods of Douglas-Rachford and Peaceman-Rachford for linear operators), for solving a general monotone inclusion problem involving the sum of two maximal monotone operators. In parallel, Cohen [START_REF] Cohen | Auxiliary problem principle and decomposition of optimization problems[END_REF] introduced the Auxiliary Problem principle to define a very general family of decomposition algorithms, generalizing the role of the Proximity step in Augmented Lagrangian functions. The necessity of reformulation to decompose separable convex programs was the motivation of Spingarn's Partial Inverse method [START_REF] Spingarn | Partial inverse of a monotone operator[END_REF]. The Partial Inverse method was then shown to be closely related to ADMM and to the Douglas-Rachford method by J. Eckstein in his thesis and a collection of related papers (see [START_REF] Eckstein | Splitting methods for monotone operators with applications to parallel optimization[END_REF][START_REF] Eckstein | On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF]). In parallel, Chen and Teboulle [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF] derived a decomposition algorithm adapted from Rockafellar's Proximal Method of Multipliers (PMM, see [START_REF] Rockafellar | Augmented lagrangians and applications of the proximal point algorithm in convex programming[END_REF]) and Mahey et al [START_REF] Mahey | Proximal decomposition on the graph of a maximal monotone operator[END_REF] analyzed the convergence rate of the Proximal Decomposition Algorithm (PDA), generalizing Spingarn's block-decomposition algorithm [START_REF] Spingarn | Applications of the method of partial inverses to convex programming:decomposition[END_REF].

In the recent years, the research on the subject has been mainly focused on the following issues :

• Revisiting forward-backward schemes and extending to some classes of non convex models ( [START_REF] Solodov | A class of decomposition methods for convex optimization and monotone variational inclusions via the hybrid inexact proximal point framework[END_REF][START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems : proximal algorithms, forwardbackward splitting and regularized gauss-seidel methods[END_REF][START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF][START_REF] Shefi | Rate of convergence analysis of decomposition methods based on the proximal method of multipliers for convex minimization[END_REF]);

• Extending splitting methods to composite operators ([13, 9, 1]);

• Exploring new strategies adapted to new models coming from Signal Processing and Classification theory ( [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex programs with applications to imaging[END_REF][START_REF] Goldfarb | Fast alternating linearization methods for minimizing the sum of two convex functions[END_REF][START_REF] Boyd | Distributed optimization and statistical learning with the alternating direction method of multipliers[END_REF]);

• Studying worst-case complexity bounds ( [START_REF] Goldfarb | Fast multiple splitting algorithms for convex optimization[END_REF][START_REF] Hong | Convergence of the alternating method of multipliers[END_REF][START_REF] Shefi | Rate of convergence analysis of decomposition methods based on the proximal method of multipliers for convex minimization[END_REF]);

• Introducing scaling parameters to accelerate convergence with selfadaptive update ( [START_REF] Dussault | Separable augmented lagrangian algorithm with multidimensional scaling for monotropic programming[END_REF][START_REF] He | Alternating direction method with self-adaptive penalty parameters for monotone variational inequalities[END_REF]).

We will focus on the last point in our presentation as it is a critical issue to derive efficient algorithms to decompose large optimization problems. For instance, it is well-known that the rate of convergence of the Douglas-Rachford method is highly sensitive to the choice of the parameter introduced in the proximal step. More worrying is the fact that the control of the rate of the primal and dual sequences are conflicting, thus limiting the best expected linear rate of convergence to the value 0.5 and provoking nasty spiralling effects (see [START_REF] Eckstein | Splitting methods for monotone operators with applications to parallel optimization[END_REF][START_REF] Hamdi | A new decomposition method in nonconvex programming via a separable augmented lagrangian[END_REF]).

On the other hand, worst-case complexity bounds have been studied in the spirit of Nesterov's smoothing techniques [START_REF] Nesterov | A method for unconstrained convex minimization problem with the rate of convergence o(1/k 2 )[END_REF] to better understand the limits of performance of splitting methods for decomposing the sum of general convex functions. Some results may look rather frustrating as commented in a recent study [START_REF] Davis | Convergence rate analysis of several splitting schemes[END_REF] : the global convergence of Douglas-Rachford splitting scheme can be as fast as the Proximal iteration in the ergodic sense and as slow as a subgradient method in the non ergodic sense.

The present study does not pretend to cover all research around operator splitting methods and applications. Since the original work of Gabay thirty years ago [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF], a few interesting tutorials and surveys have been proposed in the literature, see [START_REF] Eckstein | Augmented lagrangian and alternating direction method of multipliers : A tutorial and some illustrative computational examples[END_REF][START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Shefi | Rate of convergence analysis of decomposition methods based on the proximal method of multipliers for convex minimization[END_REF][START_REF] Komodakis | Playing with duality : an overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF].

We first present the main splitting techniques to find a zero of the sum of two maximal monotone operators. We apply the most promising techniques to the decomposition of separable convex programs in section 3, focusing on constrained optimization problems in finite dimension. In section 4, we explore theoretical and practical issues on the convergence of operator splitting methods.

2 Splitting the sum of two monotone operators

Forward and backward steps

Let recall first a few useful definitions about monotone operators. We will use the notation •, • to denote the dot product in a Hilbert space X.

Definition 1 T : X → X is monotone if T (x) -T (x ), x -x ≥ 0, ∀x, x ∈ X
It is maximal if its graph is not strictly contained in the graph of any monotone operator. It is strongly monotone with constant a > 0 if

T (x) -T (x ), x -x ≥ a x -x 2 , ∀x, x and co-coercive with module a > 0 if T (x) -T (x ), x -x ≥ a T (x) -T (x ) 2 , ∀x, x
For example, the subdifferential of a closed convex function on a convex subset of IR n is a maximal monotone operator. A symmetric linear operator which is positive definite (thus the subdifferential of a strongly convex function) is strongly monotone (the constant will be the smallest eigenvalue).

Co-coercivity is the dual property, so that an operator is co-coercive if its inverse is strongly monotone.

Let T be a maximal monotone operator on X. Let consider first the following inclusion :

Find x * such that 0 ∈ T (x * ) (1)
When T is the subdifferential operator of a closed convex function, it corresponds naturally to the global minimization of that function on X. The forward step is the following iteration :

x t+1 ∈ (II -λ t T )(x t )
which is not completely defined unless T (x t ) is a singleton. The parameter λ t > 0 is the stepsize and it is generally expected to decrease at each iteration like in the so-called subgradient algorithm for minimizing nonsmooth convex functions. The backward step is the following iteration :

x t+1 = (II + λT ) -1 (x t )
where J T λ = (II+λT ) -1 is the resolvent operator of T which is indeed defined on the whole space so that the backward iteration is now uniquely defined for any λ > 0. It is too a firmly non expansive operator (or co-coercive with modulus 1, see Minty [START_REF] Minty | Monotone operators in hilbert spaces[END_REF] and the definitions of averaged operators below) which means that :

∀x, x ∈ X, x -x , J T λ (x) -J T λ (x ) ≥ J T λ (x) -J T λ (x ) 2
This suggests that the backward step behaves like a fixed-point iteration to solve (1) which can be derived directly by :

0 ∈ T (x) ⇐⇒ 0 ∈ λT (x) ⇐⇒ x ∈ (II + λT )(x) ⇐⇒ x = (II + λT ) -1 (x)
To say more about the parameter λ > 0, the fixed-point equation above says that we have substituted the operator T by its Moreau-Yosida approximation

T λ = 1 λ (II -J T λ ).
It is shown in Brézis [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] that T λ is maximal monotone and Lipschitz with constant 1/λ. When T is the subdifferential of a closed convex function f , the iteration becomes :

x t+1 = arginf f (x) + 1 2λ x -x t 2 (2)
This implicit step can be viewed as an explicit gradient step to minimize the regularized Moreau-Yoshida function f λ (x) = inf{f (z) + 1 2λ z -x 2 } which is indeed smooth with gradient T λ (x). Thus, (2) is equivalent to

x t+1 = x t -λ∇f λ (x t ).
Implementing iteration [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems : proximal algorithms, forwardbackward splitting and regularized gauss-seidel methods[END_REF] gives the celebrated Proximal Point Method (PPM) the convergence of which was first analyzed by Martinet [START_REF] Martinet | Régularistion d'inéquations variationnelles par approximations successives[END_REF]. Rockafellar in his analysis of (PPM) concluded that the rate of convergence improves when λ increases [START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF].

The former subproblem is generally implementable when f corresponds to the dual function associated with a constrained optimization problem, yielding in the primal space the famous Augmented Lagrangian algorithm (the first detailed analysis of this relation is due to Rockafellar [START_REF] Rockafellar | Augmented lagrangians and applications of the proximal point algorithm in convex programming[END_REF]).

To complete our introduction on the basic iterative schemes, we review the notion of averaged operators, introduced by Baillon et al [START_REF] Baillon | On the asymptotic behavior of nonexpansive mappings and semi-groups[END_REF] (see too [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]). Let S be a non expansive operator on X (such that Sx -Sx ≤ x -x , ∀x, x ). A ρ-averaged operator is, for any ρ ∈ (0, 1), defined by :

S ρ = ρS + (1 -ρ)II
Averaged operators are too non expansive, but they are more amenable to fixed-point iterations than general non expansive ones as they share the following property :

S ρ x -S ρ x 2 ≤ x -x 2 - 1 -ρ ρ (II -S ρ )(x) -(II -S ρ )(x ) 2
S ρ shares the same set of fixed-points with S and the convergence of the general fixed-point iteration x t+1 = S ρt (x t ) (referred as the Krasnosel'ski-Mann algorithm) converges to a solution x * in that set (see [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF]). Moreover, it was shown in [START_REF] Davis | Convergence rate analysis of several splitting schemes[END_REF] that the sequence x t -x * is nonincreasing and that Sx t -x t 2 = o(1/t), assuming only that the sequence τ t = ρ t (1 -ρ t ) is bounded away from 0. Conveniently, compositions of averaged operators are easily seen to be averaged. Indeed, if S a and S b are averaged operators with constant a and b respectively, then

S ρ = S a •S b is averaged too with constant 1-(1-a)(1-b) < 1.
It is immediate to observe that the resolvent operator J T λ of a maximal monotone operator T is 1/2-averaged which is incidentally equivalent to be firmly non expansive. In consequence, the backward iteration in the Proximal Point method converges to a fixed-point of J T λ , i.e. a zero of T . These convergence properties of the averaged operator iterations will apply to most splitting schemes studied in the next section.

Finally, we introduce the reflector operator associated with a monotone operator T : Definition 2 Let T : IR n → IR n maximal monotone, P = (II + T ) -1 and Q = II -P . The operator

N = P -Q = 2P -II = II -2Q
is the generalized reflector associated with T .

That notion truly generalizes a symmetry corresponding to the case when T = N A is the normal cone of a linear subspace A so that Graph(T ) = A×A ⊥ . Generalized reflectors correspond exactly to the set of non expansive operators. The correspondence between the graphs of T and N appears clearly in the following construction:

Proposition 1 Let T : IR n → IR n maximal monotone and N its generalized reflector. Then y ∈ T (x) ⇐⇒ d = N (s) with :

s = x + y d = x -y or x = 1 2 (d + s) y = 1 2 (d -s)
which leads to the following decomposition:

Proposition 2 Let T : IR n → IR n maximal monotone and u ∈ IR n . Then the following assertions are equivalent :

i) y ∈ T x ii) x = P (s) and s = x + y iii) y = Q(s) and s = x + y
In other words, there exists a unique pair (x, y) ∈ Graph(T ) such that x+y = s. That decomposition on the graph of the maximal monotone operator T is called the Moreau-Minty decomposition [START_REF] Minty | Monotone operators in hilbert spaces[END_REF]. It will be at the heart of the Proximal Decomposition method presented in the next section.

Main splitting methods

We consider now the basic model of interest to derive decomposition methods, i.e. the case of T = T 1 + T 2 where T 1 and T 2 are two maximal monotone operators on X. The basic problem is then :

Find x * ∈ X such that 0 ∈ T 1 (x * ) + T 2 (x * ) (P )
One generally defines an operator splitting method as one which combines forward and backward steps applied separately to operators T 1 and T 2 but never to T 1 + T 2 . We will use throughout the following trivial reformulation in the primaldual space X × X :

x * is a solution of (1) if and only if there exists y * ∈ X such that (x * , y * ) solves :

y * ∈ T 1 (x * ) (3) -y * ∈ T 2 (x * ) (4)
Before getting into the description of the main splitting algorithms, let introduce a generic example for (P) that will be extensively used in the remainder :

Find x * ∈ A that minimizes f (x) (P 0)
where f is closed convex on X and A is a subspace of X. Typically, f is a separable function and A is the coupling subspace. The optimality conditions for (P0), assuming the existence of an optimal solution x * , are :

(x * , y * ) ∈ A × A ⊥ Graph (T ) (5)
where T = ∂f .

Double Backward splitting

As suggested by its name, that splitting scheme uses two sequential proximal steps on each operator :

x t+1 = J T 2 λt • J T 1 λt (x t ) (6)
In general, the zeroes of the composed operator J T 2 λ • J T 1 λ do not correspond to the zeroes of T 1 + T 2 . A possibility to solve (1) is to force the scaling parameter λ t to decrease to zero and characterize convergence in the following ergodic sense :

Assuming t λ t = +∞, t λ 2 t < +∞, the sequence {z t = t τ =0 λτ x τ t τ =0 λτ } converges to a zero of T (see Passty [START_REF] Passty | Ergodic convergence to a zero of the sum of monotone operators in hilbert spaces[END_REF]).

Forward-Backward splitting

Here we use a forward step associated with (3), i.e. :

x -λy ∈ (I -λT 1 )(x) composed with a backward step on (4) :

x = J T 2 λ (x -λy)
so that the forward-backward iteration is given by :

x t+1 = J T 2 λ • (I -λT 1 )(x t ) (7)
Observe that the same parameter λ should be used in both forward and backward steps at each iteration.

If T 1 = ∇F where F is smooth convex and T 2 = N C , the normal cone of a closed convex set, this is in fact the Projected Gradient method originally proposed by Goldstein [START_REF] Goldstein | Convex programming in hilbert space[END_REF]. Indeed, the resolvent J T 2 λ is the projection operator on C. Convergence was analyzed first by Passty [START_REF] Passty | Ergodic convergence to a zero of the sum of monotone operators in hilbert spaces[END_REF] who proved ergodic convergence of the sequence {z t } like in the former section. To get the convergence of the whole sequence, one needs additional properties like the co-coercivity of T 1 with modulus a and a controlled stepsize λ in the interval [0, 1/2a](see Mercier [START_REF] Mercier | Topics in finite-element solution of elliptic problems[END_REF] and related extensions in [START_REF] Chen | Convergence rates in forwardbackward splitting[END_REF]).

Figure 1 illustrates the convergence on the example of two operators T 1 , T 2 with values in IR.

Peaceman-and Douglas-Rachford splitting methods

Peaceman-Rachford method was originally applied to linear operators [START_REF] Peaceman | The numerical solution of parabolic elliptic differential equations[END_REF] and we consider here its extension to monotone operators as studied by Lions and Mercier in [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF].

Observe first that the optimality conditions (3),( 4) can be rewritten in the following way :

x * + y * = [N 2 • N 1 ](x * + y * )
where N 1 and N 2 are the reflectors associated with operators T 1 and T 2 respectively. Indeed, using Prop. 1, elementary calculations transform (3) in x * -y * = N 1 (x * + y * ) and ( 4) in x * + y * = N 2 (x * -y * ). The operator N 2 • N 1 is non expansive, inducing the following iteration :

s t+1 = [N 2 • N 1 ](s t )
so that the primal and dual sequences of iterates are obtained using the construction of Prop. [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems : proximal algorithms, forwardbackward splitting and regularized gauss-seidel methods[END_REF]. A scaling parameter can be introduced again using the following change of scale of both operators substituting T 1 by aT 1 (b -1 ) and T 2 by aT 2 (b -1 ) with a, b > 0. It corresponds to a simple change of variables, x ← b -1 x for the primal variables, and y ← a -1 y for the dual variables. The iteration is now detailed using the new primal and dual variables with the single parameter λ = b a :

Algorithm 1 Peaceman Rachford algorithm (PRA) Require: λ > 0, t = 0, choose (x 0 , y 0 )

1: x t+ 1 2 = J T 1 λ (x t + λy t ) 2: y t+ 1 2 = λ -1 (x t -x t+ 1 2 ) + y t 3: x t+1 = J T 2 λ (x t+ 1 2 -λy t+ 1 2 ) 4: y t+1 = λ -1 (x t+1 -x t+ 1 2 ) + y t+ 1 2 5: t ← t + 1; Go To step 1
Convergence of algorithm (PRA) is not guaranteed in the general case but can be stated with some additional hypotheses (the following proposition is proved in [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF]) :

Proposition 3 • If T 1 (respectively T 2 ) is strongly monotone, then x t+ 1 2
(resp. x t ) converges to the unique optimal primal solution x * .

• If T 1 (respectively T 2 ) is co-coercive, then y t+ 1 2 (resp. y t ) converges to the unique optimal dual solution y * .

Many authors, following Varga [START_REF] Varga | Matrix Iterative Analysis[END_REF], have considered a natural underrelaxation of (PRA) associated with a new parameter α > 0 :

s t+1 = [(1 -α t )II + α t (N 2 • N 1 )](s t ) (8)
Its convergence is now guaranteed for any 0 < α t < 1 as shown by Lions and Mercier in [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF]. Indeed, the iteration operator is now averaged and, moreover, the case α = 1/2 is of interest because it is exactly the Douglas-Rachford method for linear inclusions [START_REF] Douglas | On the numerical solution of the heat conduction problem in 2 and 3 space variables[END_REF] which is known to be intimately linked with the Alternate Direction method of Multipliers (ADMM) ( [START_REF] Glowinski | Sur l'approximation par éléments finis d'ordre 1 et la résolution par pénalisation-dualité d'une classe de problèmes de dirichlet[END_REF][START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF]) and the Partial Inverse method of Spingarn [START_REF] Spingarn | Partial inverse of a monotone operator[END_REF]. We will not detail here the fine correspondence between these now classical splitting methods. For a complete overview of this material, Eckstein's PhD thesis [START_REF] Eckstein | Splitting methods for monotone operators with applications to parallel optimization[END_REF] is a very accurate and partly unexploited reading. We give below the classical form of Douglas-Rachford algorithm (DRA) for finding a zero of the sum of two maximal monotone operators.

Algorithm 2 Douglas Rachford algorithm (DRA)

Require: λ > 0, t = 0, choose (x 0 , y 0 )

1: x t+ 1 3 = J T 1 λ (x t + λy t ) 2: y t+ 1 3 = λ -1 (x t -x t+ 1 3 ) + y t 3: x t+ 2 3 = J T 2 λ (x t+ 1 3 -λy t+ 1 3 ) 4: y t+ 2 3 = λ -1 (x t+ 2 3 -x t+ 1 3 ) + y t+ 1 3 5: (x t+1 , y t+1 ) = 1 2 [(x t+ 2 3 , y t+ 2 3 ) + (x t , y t )] 6: t ← t + 1; Go To step 1
The only difference with (PRA) is the addition of step 5 where we substitute the last primal-dual estimate in (PRA) by the mean of the two successive estimates during the iteration. Observe that the second proximal step on T 2 flips the sign of the dual counterpart y. To be more precise as observed in [START_REF] Eckstein | Splitting methods for monotone operators with applications to parallel optimization[END_REF], we have the following proposal : Proposition 4 At each step of algorithm (DRA), the sequence {z t } given by z t = x t + λy t satisfies :

J T 2 λ (z t ) = x t (9) (II -J T 2 λ )(z t ) = λy t (10)
N T 2 λ (z t ) = x t -λy t (11)
so that (x t , y t ) ∈ Graph(T 2 ). The role of the new variable z t = x t +λy t , which combines the primal and dual iterates, is central as the sequence generated by (DRA) satisfies :

z t+1 = [J T 1 λ • N T 2 λ + II -J T 2 λ ](z t )
Lenoir [START_REF] Lenoir | Modèles et algorithmes pour la planification de production à moyen terme en environnement incertain[END_REF] has shown that the averaging step 5 of algorithm 2 could be performed between the proximal steps without changing the convergence properties of the method, thus allowing to play with the modeling of the two operators composing the inclusion problem. The sequence of iterates for (DRA) is illustrated on Figure 2.

Proximal decomposition method

We close now this section with the presentation of the Proximal Decomposition on the graph of a maximal monotone operator, motivated by Prop.

(2). The main idea is that the two proximal steps used in (DRA) can be performed in parallel, followed by the averaging step which is itself a proximal step performed on an appropriate coupling subspace. Having in mind the generic model (P 0), this can be done by duplicating the space X and creating two copies x 1 and x 2 of the original variables x to get an equivalent formulation :

Find (x 1 , x 2 ) ∈ X × X | 0 ∈ T 1 (x 1 ) + T 2 (x 2 ) and x 1 = x 2 Let denote by A = {(x 1 , x 2 ) ∈ X × X | x 1 =
x 2 } the coupling subspace. Observe that the dual variables can again be introduced to write the optimality condition in the following way :

y * 1 ∈ T 1 (x * 1 ) y * 2 ∈ T 2 (x * 2 ) (y * 1 , y * 2 ) ∈ A ⊥
= J T 1 λ (x * 1 + λy * 1 ) and x * 2 = J T 2 λ (x * 2 + λy * 2 )
, inducing the alternate proximal steps of the following algorithm proposed in Mahey et al [START_REF] Mahey | Proximal decomposition on the graph of a maximal monotone operator[END_REF] under the name 'Proximal Decomposition method' : Algorithm 3 Proximal Decomposition algorithm (PDA) Require: λ > 0, t = 0, choose (x 0 , y 0 ) and set x 0 1 = x 0 2 = x 0 and y 0 1 = -y 0 2 = y 0 1: x

t+ 1 2 1 = J T 1 λ (x t 1 + λy t 1 ) 2: y t+ 1 2 1 = λ -1 (x t 1 -x t+ 1 2 1 ) + y t 1 3: x t+ 1 2 2 = J T 2 λ (x t 2 + λy t 2 ) 4: y t+ 1 2 2 = λ -1 (x t 2 -x t+ 1 2 2 ) + y t 1 5: x t+1 1 = x t+1 2 = 1 2 (x t+ 1 2 1 + x t+ 1 2 2 ) 6: y t+1 1 = -y t+1 2 = 1 2 (y t+ 1 2 1 -y t+ 1 2 2 ) 7: t ← t + 1; Go To step 1
To resume, (DRA) alternates proximal steps on T 1 and T 2 in a Gauss-Seidel fashion, whereas (PDA) produces the same proximal steps in parallel Figure 3: Proximal decomposition algorithm in a Jacobi fashion. But (PDA) can also be interpreted as the application of (DRA) to the following inclusion in X × X :

Find (x * 1 , x * 2 ) ∈ X × X such that 0 ∈ A(x * 1 , x * 2 ) + B(x * 1 , x * 2 )
where

A = T 1 0 0 T 2 and Graph(B) = A × A ⊥ .
The sequence of iterates for (PDA) is illustrated on Figure 3. The name 'Proximal Decomposition' is justified by the application of that scheme to model (P0). There is no need to introduce the copies of primal variables, simplifying the above steps in :

Step 1 : (x t+ 1 2 , λy t+ 1 2 ) is the Moreau-Minty Decomposition of x t + λy t on the graph of T .

Step 2 : (x t+1 , y t+1 ) = (x

t+ 1 2 A , y t+ 1 2 A ⊥ )
where x A denotes the projection of x on the subspace A.

(PDA) shares a strong link with a former method introduced by Spingarn, the Partial Inverse method [START_REF] Spingarn | Partial inverse of a monotone operator[END_REF]. Motivated by model (P0), he introduced the Partial Inverse operator associated with the maximal monotone operator T and the subspace A, defined by its graph :

Graph (T A ) = {(x A + y A ⊥ , x A ⊥ + y A ) | y ∈ T (x)}
T A is maximal monotone if and only if T is so, and, moreover, we have

0 ∈ T A (s) ⇐⇒ s A ⊥ ∈ T (s A )
In other words, the projections of a zero of T A on the orthogonal subspaces A and A ⊥ are the primal and dual solutions of (P0) respectively. Thus, Spingarn proposed to apply (PPM) to operator T A to solve (P0). But the backward iteration applied to the Partial Inverse operator reduces to :

s t+1 = (II + T A ) -1 ⇐⇒ s t -s t+1 ∈ T A (s t+1 )
Using the definition of the graph of T A , this implies that there exists (x, y) ∈ Graph (T ) such that s t+1 = x A + y A ⊥ and s t -s t+1 = x A ⊥ + y A , or equivalently : s t = x + y (so that (x, y) is the Proximal Decomposition of s t on the graph of T ) and

s t+1 = x A + y A ⊥ .
Observe that the relationship between (PDA) and Spingarn's Partial Inverse needs to set the λ parameter to 1. Eckstein has derived too the direct relationship between (DRA) and the Partial Inverse method in [START_REF] Eckstein | Splitting methods for monotone operators with applications to parallel optimization[END_REF] using appropriate changes of variables.

The main advantage of (PDA) is that it can easily be extended to inclusions with more than two operators. Again, if the inclusion involves p maximal monotone operators T 1 , . . . , T p on X :

Find x * ∈ X such that 0 ∈ T 1 (x * ) + • • • + T p (x * )
we can create p copies of space X and apply (PDA) to the cross-product T = T 1 × . . . × T p which is indeed a maximal monotone operator on X p over the coupling subspace • To reduce the dimension of large-scale optimization problems with different interconnected subsystems; the challenge here is to identify the coupling variables and/or constraints.

A = {(x 1 , . . . , x p ) ∈ X p | x 1 = • • • = x p } :
• To exploit hidden 'easy' submodels or, equivalently, to isolate the 'hard' features of the model without which the model is solvable by ad hoc software.

Reformulations of the model are frequently necessary to identify these situations. A typical example for the first case is block-angular linear programs which gave rise to the first decomposition algorithms for Operations Research like Dantzig-Wolfe's and Benders' decomposition methods (see Lasdon's textbook [START_REF] Lasdon | Optimization for Large Systems[END_REF] for example). In this case, the separable coupling constraints (or variables) contain the hard features of the model as they prevent to try solving the block subproblems separately. These strategies lead frequently to implicit value functions (like the dual function with Lagrangian Relaxation) which are typically nonsmooth. As was mentioned before, the inherent non smoothness of the dual function is the main motivation of the Proximal Point method, leading to Augmented Lagrangian subproblems. The problem is that separability which allows decomposing the Lagrangian subproblems is destroyed by the quadratic terms introduced in the Augmented Lagrangian function. Operator splitting methods will be able to address that issue.

Separable Augmented Lagrangian

We consider first a general convex minimization problem in IR n :

     Minimize p i=1 f i (x) x ∈ S (P 1)
where the f i are extended real valued convex functions supposed to be proper and lower-semi-continuous on the closed convex set S ⊂ IR n . Additional local constraints may be present in the model, here modeled inside the functions f i . A convenient special case for illustrating the next methods is the problem of finding x ∈ i=1,...,p C i where the C i are closed convex sets. A simple way to decouple the p pieces of the objective function is to introduce p copies of the variable x denoted ξ i , i = 1, . . . , p and reformulate problem (P 1) in the product space (IR n ) p :

         Minimize p i=1 f i (ξ i ) ξ i ∈ S, i = 1, . . . , p ξ = (ξ 1 , . . . , ξ p ) ∈ A (12)
where

A = {ξ ∈ (IR n ) p | ξ 1 = • • • = ξ p } is the coupling subspace.
Thus, (P 1) has been reformulated into the generic model (P 0). That reformulation has been early used by Pierra [START_REF] Pierra | Decomposition through formalization on a product space[END_REF] who applied the doublebackward splitting to [START_REF] Briceno-Arias | Forward douglas-rachford splitting and forward partial inverse method for solving monotone inclusions[END_REF] (recall that this requires the parameter to decrease to zero).

The application of algorithm (PDA) to ( 12) is straightforward, working in IR np with primal ξ = (ξ 1 , . . . , ξ p ) and dual variables ζ = (ζ 1 , . . . , ζ p ). The monotone operator in the product space will be the cartesian product of the subdifferential operators

∂f 1 × • • • × ∂f p : Algorithm 4 PDA-separable Require: t = 0, λ > 0, > 0, ξ 0 ∈ A, ζ 0 ∈ A ⊥ 1: repeat 2:
for all i = 1, . . . , p do 3:

ξ t+ 1 2 i := arg min ξ i ∈S f i (ξ i ) + 1 2λ ξ i -ξ t i -λζ t i 2 4: ζ t+ 1 2 i := λ -1 (ξ t i -ξ t+ 1 2 i ) + ζ t i 5: ξ t+1 i ← 1 p i ξ t+ 1 2 i 6: ζ t+1 i ← ζ t+ 1 2 i -1 p i ζ t+ 1 2 i 7:
end for 8:

t ← t + 1 9: until ξ t+1 -ξ t + ζ t+1 -ζ t <
The algorithm (PDA-separable) can be applied to many structured models, but a typical situation which we will describe now is the case of a separable objective with separable coupling constraints.

We are interested in solving the following convex program, called hereafter the S-Model, defined on the product space

IR n = IR n 1 × . . . × IR np .                Minimize p i=1 f i (x i ) p i=1 g i (x i ) = 0 x i ∈ S i , i = 1, . . . , p (S -M odel)
where f i are extended real valued convex functions on closed convex sets S i , supposed to be proper and lower-semi-continuous (l.s.c) and g i are affine:

g i : IR n i → IR m x i → g i (x i ) = G i x i -b i
where G i are (n i × m) not necessarily full-rank matrices and b i ∈ IR m .

We will now apply (PDA-separable) to the Lagrangian dual of the S-Model, i.e.

Maximize u∈IR m

p i=1 h i (u)
where h i (u) = inf x i ∈S i f i (x i ) + u, g i (x i ) . The concave dual problem is thus in the form of problem (P1), so that it can be reformulated by creating p copies of the dual variables u = ξ 1 = • • • = ξ p to get an equivalent model in

IR mp : Maximize ξ=(ξ 1 ,...,ξp)∈A p i=1 h i (ξ i ) (13)
where

A = {(ξ 1 , . . . , ξ p ) ∈ IR mp | ξ 1 = • • • = ξ p } is the coupling subspace. Let ζ = (ζ 1 , .
. . , ζ p ) be the corresponding variables in duality relation with ξ. The application of (PDA) will give the following algorithm :

Algorithm 5 PDA-dualseparable Require: t = 0, λ > 0, > 0, ξ 0 ∈ A, ζ 0 ∈ A ⊥ 1: repeat 2:
for all i = 1, . . . , p do 3:

ξ t+ 1 2 i := arg max h i (ξ i ) -1 2λ ξ i -ξ t i -λζ t i 2 4: ζ t+ 1 2 i := λ -1 (ξ t i -ξ t+ 1 2 i ) + ζ t i 5: ξ t+1 i ← 1 p i ξ t+ 1 2 i 6: ζ t+1 i ← ζ t+ 1 2 i -1 p i ζ t+ 1 2 i 7:
end for 8:

k ← k + 1 9: until ξ t+1 -ξ t + ζ t+1 -ζ t <
Observe that the quadratic term added in the proximal step 3 of algorithm 5 is here subtracted as the dual problem is a maximization problem.

The algorithm can be developed in the primal setting yielding a separable Augmented Lagrangian algorithm early proposed by Spingarn [START_REF] Spingarn | Applications of the method of partial inverses to convex programming:decomposition[END_REF]. Indeed, the optimality conditions of the proximal step are :

ξ t+ 1 2 i = ξ t i + λζ t i + λg i (x t+ 1 2 i )
where g i (x

t+ 1 2 i ) ∈ ∂h i (ξ t+ 1 2 i
) and x in the ordinary Lagrangian function L i (x i , ξ i ) = f i (x i ) + ξ i , g i (x i ) , which gives the following subproblem :

Minimize x i ∈S i f i (x i ) + ξ t i , g i (x i ) + λ 2 g i (x i ) + ζ t i 2
The presentation of the algorithm can be simplified, avoiding the use of auxiliary variables ξ as ξ t = (u t , . . . , u t ) after the projection on A and integrating the intermediate step t + 1 2 by observing that

u t+1 i = u t + 1 p r(x t+1 ) y t+ 1 2 i = -g i (x t+1 i )
where r(x) = i g i (x i ) is the residual of the relaxed coupling constraints. The complete algorithm, called (SALA) for Separable Augmented Lagrangian Algorithm, is then :

Algorithm 6 SALA : a Separable Augmented Lagrangian Algorithm Require: t = 0, λ > 0, > 0, u 0 ∈ IR m , ζ 0 = (y 0 1 , . . . , y 0 p ) ∈ A ⊥ 1: repeat 2:
for all i = 1, . . . , p do 3:

x t+1 i :∈ arg min x i ∈S i f i (x i ) + u t , g i (x i ) + λ 2 g i (x i ) + y t i 2 4:
end for 5:

r t+1 ← p i=1 g i (x t+1 i ) 6:
for all i = 1, . . . , p do 7:

y t+1 i ← -g i (x t+1 i ) + 1 p r t+1 8:
end for 9:

u t+1 ← u t + λ p r t+1 10:

t ← t + 1 11: until r(x t+1 ) <
Observe that the variables ζ must lie in the orthogonal subspace

A ⊥ = {ζ = (ζ 1 , . . . , ζ p ) | i ζ i = 0}.
They are exactly the right-hand side allocations ζ i = -g i (x i ) used in resource-directive primal decomposition (see Lasdon for example [START_REF] Lasdon | Optimization for Large Systems[END_REF]). Indeed, the primal separable counterpart of ( 13) is :

Minimize ζ=(ζ 1 ,...,ζp)∈A ⊥ p i=1 v i (ζ i ) (14)
where

v i (ζ i ) = inf{f i (x i ) | g i (x i ) = -ζ i , x i ∈ S i }

is the implicit primal block function (convex on a convex domain with the current hypotheses).

In other words, (SALA) is exactly the application of an Augmented Lagrangian algorithm to the resource-directive reformulation [START_REF] Chambolle | A first-order primal-dual algorithm for convex programs with applications to imaging[END_REF] followed by the projection steps on the respective subspaces.

We observe here that Separable Augmented Lagrangian algorithms have been early proposed in the literature (see [START_REF] Stephanopoulos | The use of Hestenes method of multipliers to re-solve dual gaps in engineering system optimization[END_REF][START_REF] Bertsekas | Convexification procedures and decomposition methods for nonconvex optimization problems[END_REF][START_REF] Cohen | Decomposition coordination methods in largescale optimization problems. the nondifferentiable case and the use of augmented lagrangians[END_REF] for instance), but these methods do not rely on splitting schemes, rather linearizing the non separable terms to give rise to three-levels decomposition schemes with extensions to the non convex case.

Alternate direction method of multipliers

We consider now another situation involving two different convex functions. The functions can be smooth or not and are generally composite as they include linear transformations of the variables, so that many authors have analyzed the following M-model :

Minimize f (x) + g(M x) (P 2)
with f : IR n → IR strongly convex, g : IR m → IR simply convex and M is a (generally full-rank) (m×n) matrix. A dual formulation is also convenient as the conjugate function of f (i.e. f * (y) = sup{ y, x -f (x)}) is differentiable and the dual problem associated with (P 2) presents the same nice structure as the primal :

Minimize f * (M T u) + g * (-u) (15)
Even if this model can be reformulated as {min f * (v) + g * (-u) | v + M T u = 0}, thus minimizing the sum of two convex functions on a linear subspace, the direct application of the Forward-Backward splitting is natural using T 1 = ∂f * and T 2 = ∂g * • (-M T ). The resulting algorithm (see [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF] and [START_REF] Tseng | Applications of a splitting algorithm to decomposition in convex programming and variational inequalities[END_REF]) is given below :

Algorithm 7 FB-M Require: t = 0, λ > 0, > 0, x 0 , z 0 , u 0 1: x t+1 := arg min f (x) -u t , M x 2: z t+1 := arg min g(z) + u t , z + λ 2 z -M x t+1 2 3: u t+1 := u t -λ(M x t+1 -z t+1 ) 4: t ← t + 1
Convergence properties of the FB-M scheme have been analyzed by Chen and Rockafellar [START_REF] Chen | Convergence rates in forwardbackward splitting[END_REF].

Application of (DRA) to the M -Model leads naturally to what is generally referred to as the Alternate Direction Method of Multipliers (ADMM) originally studied by Gabay [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF] (see too [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]).

Algorithm 8 ADMM

Require: t = 0, λ > 0, > 0, x 0 , z 0 , u 0 1:

x t+1 := arg min f (x) -u t , M x + λ 2 z t -M x 2 2: z t+1 := arg min g(z) + u t , z + λ 2 z -M x t+1 2 3: u t+1 := u t -λ(M x t+1 -z t+1 ) 4: t ← t + 1
The structural links between (ADMM), (DRA) and the Proximal Point method (PPM) have been extensively detailed by Eckstein (see [START_REF] Eckstein | On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Eckstein | Some reformulations and applications of the alternating direction method of multipliers[END_REF]). One recurrent question which remains partly open with (ADMM) is the extension to more than two blocks of variables. As mentioned early by Lions and Mercier [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF], the Douglas-Rachford splitting method does not naturally generalize to more than two operators. Even if an extension of (DRA) with n operators has been early proposed by Douglas and Gunn [START_REF] Douglas | A general formulation of alternating direction methods[END_REF], the sequential steps in (ADMM) turn the convergence analysis more intricate. A recent study has even exhibited an academic example with three blocks which is indeed divergent [START_REF] Chen | The direct extension of admm for multi block convex minimization problems is not necessarily convergent[END_REF]. Of course, the question can be turned around as commented before, by including the separable functions into a single dual operator and forcing the dual copies in a coupling subspace playing the role of the second operator, as in the (PDA) algorithm. Eckstein and Svaiter [START_REF] Eckstein | General projective splitting methods for sums of maximal monotone operators[END_REF] have introduced a projective splitting method which avoids the introduction of copies for the coupling variables when an arbitrary number of monotone operators are present in the model.

There is still a comment on ADMM assumptions and rates of convergence. Many authors have studied the convergence properties of ADMM (which assumptions on the functions and, with which rates of convergence ?) as shown in the recent monograph [START_REF] Boyd | Distributed optimization and statistical learning with the alternating direction method of multipliers[END_REF], beyond the natural results inherited from (DRA) as early studied by Lions and Mercier [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF]. It appears that linear convergence is possible when at least one of the block function is strongly convex. A recent study by Hong and Luo [START_REF] Hong | Convergence of the alternating method of multipliers[END_REF] made new progress in relaxing these assumptions and extending too to more than two blocks. Their separable model extends the M-model to include the S-model used in the former section. It also uses block functions of the form

f i (x i ) = f 0i (A i x i ) + f 1i (x i )
where f 0i are strictly convex and smooth and f 1i are typically polyhedral. Some steering matrices A i can be zero allowing not strongly convex blocks. The coupling constraints should be linearly independent, i.e. the matrix

[G 1 | • • • |G p ] is full-rank,
but some G i may not. These assumptions allow to prove linear convergence of the primal and dual sequences as well as the sequence of the coupling constraints violation. Similar results in earlier papers exploiting the strongly convex parts in the objective function can be cited here, like [START_REF] Han | A parallel algorithm for a class of convex programs[END_REF][START_REF] Mouallif | A perturbed parallel decomposition method for a class of nonsmooth convex minimization problems[END_REF][START_REF] Kiwiel | Proximal decomposition via alternating linearization[END_REF].

Inspired by multi-block Gauss-Seidel algorithm, Goldfarb and Ma [START_REF] Goldfarb | Fast multiple splitting algorithms for convex optimization[END_REF] proposed a generalized ADMM which converges in O(1/ √ ) iterations to obtain an approximation within of the optimal value but requires that all functions are smooth. Other generalized splitting strategies based on blockcoordinate iterations have been recently proposed by Combettes et al [START_REF] Combettes | Stochastic quasi-fejér blockcoordinate fixed-point iterations with random sweeping[END_REF].

We have seen in this section the impact of problem reformulation to get constructive variants of the basic splitting schemes. The literature has focussed mainly on the M-Model and its dual version [START_REF] Chen | The direct extension of admm for multi block convex minimization problems is not necessarily convergent[END_REF]. An alternative and interesting primal-dual setting has been recently proposed by Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex programs with applications to imaging[END_REF] who used a saddle-point formulation :

inf x sup y f (x) -g * (y) + y, M x
and proved convergence of the following splitting algorithm : Algorithm 9 CP Require: t = 0, λ > 0, > 0, x 0 , y 0 1: x t+1 := arg min f (x) + 1 2λ x -x t + λM T y t 2 2: y t+1/2 := arg min g * (y) + 1 2µ y -y t -µM x t+1 2 3: y t+1 := y t+1/2 + θ(y t+1/2 -y t ) 4: t ← t + 1 where 0 ≤ θ ≤ 1 and scaling parameters chosen such that λµ M 2 < 1.

Convergence to a saddle-point (x * , y * ) was proved in the ergodic sense with rate O(1/t).

Convergence results and complexity issues

We close this survey by inspecting practical issues concerning operator splitting techniques. Before relating the different extensions which have been investigated, and they are quite numerous, it is of interest to overview the areas of applications where splitting methods have been used successfully. Historically, R. Temam [START_REF] Temam | Sur la stabilité et la convergence de la méthode des pas fractionnaires[END_REF], R. Glowinski [START_REF] Glowinski | Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics[END_REF], Gabay and Mercier [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite-element approximations[END_REF] among many others, have used (DRA) to solve evolution equations in Mechanics, referring the approach as a penalty-dualization method. As the formal setting involves inclusion problems with general monotone operators, the application of splitting methods to Variational Inequalities is still an attractive subject for applied mathematicians. Besides these applications, the report of Bensoussan et al [START_REF] Bensoussan | Sur les méthodes de décomposition, décentralisation et de coordination et applications[END_REF] and more recently, the textbook by Bertsekas and Tsitsiklis [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF] have largely contributed to disseminate these techniques to the areas of Mathematical Programming and Operations Research where decomposition techniques are very popular since the sixties. Among many different areas of applications, we can cite Multicommodity Flow problems with convex arc costs ( [START_REF] Mahey | A new proximal decomposition algorithm for routing in telecommunication networks[END_REF], [START_REF] Fukushima | Application of the alternating direction method of multipliers to separable convex programming[END_REF]), Stochastic Programming (adapting (DRA) to two-stage stochastic optimization with recourse leads to the Progressive Hedging method of Rockafellar and Wets [START_REF] Rockafellar | Scenarios and policy aggregation in optimization under uncertainty[END_REF]), Fermat-Weber problems (the Partial Inverse of Spingarn was applied to a polyhedral operator splitting model in [START_REF] Idrissi | Applications and numerical convergence of the partial inverse method[END_REF]). More recently, new models received a lot of interest in the areas of Image Reconstruction and Signal Processing ( [START_REF] Chambolle | A first-order primal-dual algorithm for convex programs with applications to imaging[END_REF][START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF]), with similar models in Classification problems [START_REF] Goldfarb | Fast alternating linearization methods for minimizing the sum of two convex functions[END_REF][START_REF] Boyd | Distributed optimization and statistical learning with the alternating direction method of multipliers[END_REF]. These models involve in general the combination of two norms, one being smooth but not the other one, inducing the use of the M-model. For instance, the lasso or compressed sensing problem considers f (x) = x 1 and g(x) = 1 2µ y -Ax 2 , where A is a huge sparse matrix, with applications in deblurring images [START_REF] Chambolle | A first-order primal-dual algorithm for convex programs with applications to imaging[END_REF], classifying big data [START_REF] Goldfarb | Fast alternating linearization methods for minimizing the sum of two convex functions[END_REF] or matrix rank minimization [START_REF] Fazel | Hankel matrix rank minimization with applications to system identification and realization[END_REF]. More references concerning these applications to Image and Signal Processing can be found in [START_REF] Komodakis | Playing with duality : an overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF].

Motivated by these applications, new versions of the classical splitting schemes have been proposed. As a striking example, we note the recent primal-dual splitting, inspired by the Forward-Backward scheme, on inclusion problems involving composite operators like 0 ∈ Ax + K * BKx where A and B are general maximal monotone operators and K is a linear continuous operator with adjoint K * . The convergence of the explicitly composite case was first analyzed by Briceño-Arias and Combettes in [START_REF] Briceno-Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF] (see too [START_REF] Bot | A primal-dual splitting for finding zeros of sums of maximally monotone operators[END_REF] and [START_REF] Alghamdi | A primal-dual method of partial inverses for composite inclusions[END_REF] for an extension of Spingarn's Partial Inverse method).

Another line of rich theoretical studies concerns the exploitation of smooth-ness in the models. This induces a extension of model (P 0 ) to :

Minimize x∈A f (x) + h(x)
where f is still convex and A a coupling subspace, but h is convex and smooth with a Lipschitzian gradient. Smoothness is better exploited by Forward-Backward schemes and Briceño-Arias has proposed a Forward-DRA scheme in [START_REF] Briceno-Arias | Forward douglas-rachford splitting and forward partial inverse method for solving monotone inclusions[END_REF], further improved by Davis in [START_REF] Davis | Convergence rate analysis of forward-douglas-rachford splitting scheme[END_REF] (see too [START_REF] Tseng | A modified forward-backward plotting method for maximal mono tone mappings[END_REF] or [START_REF] Chouzenoux | Variable metric forwardbackward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF] for similar related algorithms). We will present first some algorithmic enhancements relative to convergence issues and, in the second part, discuss numerical scaling issues.

Algorithmic enhancements

Many variants of the basic schemes have been analyzed in the literature and we focus here on a few important choices that can produce new decomposition methods rather than on the variety of models to which these schemes may be applied. Most of these enhancements correspond to known variants of the Proximal Point method itself.

We will discuss the following issues :

• the use of relaxation parameters;

• the introduction of additional regularizing terms;

• approximate solutions in the proximal steps;

The first idea has been early proposed by Glowinski and Marocco in the original presentation of (ADMM) [START_REF] Glowinski | Sur l'approximation par éléments finis d'ordre 1 et la résolution par pénalisation-dualité d'une classe de problèmes de dirichlet[END_REF], substituting the dual update in Algorithm 8 by :

u t+1 := u t -γλ(M x t+1 -z t+1 )
where γ is a relaxation parameter, indeed corresponding to :

ũt+1 = u t -λ(M x t+1 -z t+1 ) u t+1 := γ ũt+1 + (1 -γ)u t
and they proved that convergence is guaranteed for 0 < γ < 1+ √ 5

2 . Observe that this strategy should be compared to the generalization of (PRA) defined in [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]. On the other hand, most users of (DRA) have fixed γ = 1 to focus on the estimation of the scaling parameter λ as seen below. More insight on relaxed versions of (FB), (DRA) and (PRA) and their theoretical rates of convergence may be found in the recent study by Davis and Yin and companion papers [START_REF] Davis | Convergence rate analysis of several splitting schemes[END_REF].

The second aspect concerns the addition of proximal terms in the primal and/or dual update formulae as early proposed for (PPM) by Rockafellar [START_REF] Rockafellar | Augmented lagrangians and applications of the proximal point algorithm in convex programming[END_REF] in the Proximal method of Multipliers. The corresponding extension of the (FB-M) algorithm was proposed by Chen and Teboulle [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF] and a similar idea was used by Eckstein [START_REF] Eckstein | Some saddle-point splitting method for convex programming[END_REF] to extend (ADMM) (see too [START_REF] Xu | Proximal alternating directions method for structured variational inequalities[END_REF][START_REF] Attouch | Augmented lagrangian and proximal alternating direction method of multipliers in hilbert spaces. applications to games, pde's and control[END_REF]). A recent survey by Shefi and Teboulle [START_REF] Shefi | Rate of convergence analysis of decomposition methods based on the proximal method of multipliers for convex minimization[END_REF] recalled the complexity issues to improve global convergence ratios, in particular for regularized versions of the splitting schemes. The Proximal-ADMM can be sketched as the following algorithm :

Algorithm 10 Prox-ADMM Require: t = 0, λ > 0, > 0, x 0 , z 0 , u 0 1:

x t+1 := arg min f (x) -u t , M x + λ 2 z t -M x 2 + 1 2µ 1 x -x t 2 2: z t+1 := arg min g(z) + u t , z + λ 2 z -M x t+1 2 + 1 2µ 2 z -z t 2 3: u t+1 := u t -λ(M x t+1 -z t+1 ) 4: t ← t + 1
where µ 1 , µ 2 > 0 drive the proximal terms for the {x t } and {z t } sequences.

Global convergence of the Proximal-ADMM algorithm and variants are considered in [START_REF] Shefi | Rate of convergence analysis of decomposition methods based on the proximal method of multipliers for convex minimization[END_REF] with refined convergence rates. In the general case without further assumptions on f and g, typical ergodic convergence is exhibited with O(1/ √ t) global rate. Assuming g is convex Lipschitz continuous, then the whole primal-dual sequence converges with O(1/t) global rate.

Approximate solutions in the proximal steps have been considered early by Rockafellar [START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF] and further introduced in many splitting schemes (see [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF]). A typical inexact version of (DRA) will be :

Algorithm 11 Inexact-ADMM Require: t = 0, λ > 0, 0 > 0, x 0 , z 0 , u 0 1:

x t+1 :≈ 1t arg min f (x) -u t , M x + λ 2 z t -M x 2 2: z t+1 :≈ 2t arg min g(z) + u t , z + λ 2 z -M x t+1 2 3: u t+1 := u t -λ(M x t+1 -z t+1 ) 4: t ← t + 1
where a :≈ b is a shorthand for a -b ≤ . As in the inexact version of the Proximal Point method, convergence is maintained if the errors satisfy ∞ t=0 t < +∞ (see [START_REF] Eckstein | On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF]).

Scaling and numerical enhancements

Many authors have revisited the basic splitting methods discussed in the former sections to improve convergence results and obtain implementable algorithms with optimized performance. As mentioned before, the introduction of a second parameter to generate a family of Peaceman-Rachford pure iteration appeared in Lions and Mercier's paper [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF]. Later, Spingarn in [START_REF] Spingarn | Applications of the method of partial inverses to convex programming:decomposition[END_REF] suggested to use a weighted scalar product to scale the projection step, but the difficulty to adjust many parameters during the convergence process reduced the attempts to implement these ideas. By the way, many authors have considered the introduction of different parameters associated with scaling matrices in the same spirit of what has been investigated for Augmented Lagrangian algorithms. The particular sensitivity to parameter λ of the supposedly most efficient splitting, i.e. (DRA) or (PDA), deserves special attention.

Convergence rates for a single scaling parameter

To analyze the role of the scaling parameter, we use here the Proximal Decomposition setting where T 1 = T is a maximal monotone operator, coercive with constant ρ and Lipschitz with constant L, and T 2 is the subdifferential of the indicator of a linear subspace, i.e. Gr(T 2 ) = {(x, y) ∈ X × X | (x, y) ∈ A × A ⊥ }. Recall that this corresponds to minimizing a strongly convex function (possibly separable) on a coupling subspace. The scaled proximal decomposition on the graph of T with parameter λ > 0 corresponds to the following step : given (x, y) ∈ X × X, find the unique pair (u, v) ∈ Gr(T ) such that u + λv = x + λy = z. The corresponding proximal steps on the primal and dual variables are :

u = (II + λT ) -1 (z) v = λ -1 (z -u))
which are followed by the projection steps on A×A ⊥ . It was proved by Lions and Mercier [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF] for the general case (T 2 is a maximal monotone operator) and revisited by Mahey et al [START_REF] Mahey | Proximal decomposition on the graph of a maximal monotone operator[END_REF] for the proximal decomposition that linear convergence is guaranteed with a convergence rate

r u (λ) = 1 - 2λρ (1 + λL) 2
inducing the best choice for λ = 1/L with optimal rate 1 -ρ 2L . Observe now that the update of the dual variable y is equivalent to :

v = (II + µT -1 ) -1 (µz)
with λµ = 1. We can then derive an upper bound for the rate of the dual sequence as we did for the primal one, noting that the coercivity constant for T -1 is 1/L and its Lipschitz constant is 1/ρ. Straightforward calculations give the convergence rate (of the dual sequence) : with the best choice λ v = 1/ρ giving the optimal rate 1 -ρ 2L . Then the optimal rates are equal but for different values of the parameter. Figure 4 (taken from [START_REF] Mahey | Adaptive scaling and convergence rates of a separable augmented lagrangian algorithm[END_REF]) shows the two rates as function of λ. The best compromise is to minimize max(r u (λ), r v (λ) which yields

r v (λ) = 1 - 2λρ 2 L(1 + λρ)
λ * = 1 √ ρL
Following the analysis in [START_REF] Mahey | Adaptive scaling and convergence rates of a separable augmented lagrangian algorithm[END_REF], we can obtain lower and upper bounds for both primal and dual sequences using the properties of operator T :

1 1 + λL z -z 2 ≤ u -u 2 ≤ 1 1 + λρ z -z 2 λρ 1 + λρ z -z 2 ≤ u -u 2 ≤ λL 1 + λL z -z 2
Setting λ = λ * , we get the same bounds (LB for the lower bound and UB for the upper bound) for both sequences : Ideally, we would expect a rate of 0.5 for both sequences when λ is tuned to its compromise value λ * . Figure 5 compares the behaviour of primal and dual rates (estimated during the tail of iterations) τ u , and τ z respectively, with respect to parameter λ, and the primal-dual rate τ uz = lim

LB = 1 1+ L ρ ≤ 0.5 U B = 1 1+ √ ρ L ≥ 0.5
  u t+1 z t+1   -   u * z *     u t z t   -   u * z *  
, the latter keeping a deceivingly slow pace. This rather frustrating limiting rate of 0.5 induces us to explore the possibility of multidimensional parameters, as discussed below.

Another drawback, early observed by Eckstein [START_REF] Eckstein | Splitting methods for monotone operators with applications to parallel optimization[END_REF] is spiralling which tends to slow down the iterates in the neighborhood of a fixed point (a typical behaviour is shown in Figure 7). This phenomenon was currently observed when splitting is applied to polyhedral (thus not strongly monotone) models.

Multidimensional scaling

It is easy to extend the former splitting methods to a multidimensional scaling strategy. To understand the transformation, for a positive definite matrix M , let consider a variable change z = M x and substitute the monotone operator T by T = M -T • T • M -1 . T is indeed maximal monotone if T is so and their graphs correspond in the following way :

y ∈ T x ⇐⇒ u ∈ T (z) for u = M -T y and z = M x
Consequently, a backward step with the scaled operator T derived back in the original x-space using the inverse transformation x = M -1 z, corresponds to the resolvent operator (II + ΛT ) -1 where Λ = M T M .

Setting the scaling parameter λ = 1, it appears that the matrix Λ plays the role of the scaling parameter in the splitting method. Convergence result are identical but different convergence rates are expected if the scaling matrix Λ is adequately chosen. Most practical applications of multidimensional scaling suggest to use a diagonal scaling to turn the estimation of the parameters easier. The choice of the scaling matrix depends on the model (or its reformulation) to maintain the decomposition features of the method. For instance, considering the S-Model, the scaling corresponds to substituting each resource allocation in (S -M odel) by M g i (x i ) + y i = 0. It is shown on a simple quadratic example with two coupling constraints that two different parameters can drastically improve convergence (see Figure 6 illustrating the rate of convergence of the primal-dual sequences with a two-dimensional scaling).

Many applications of operator splitting techniques include preconditioning techniques which is indeed equivalent to multidimensional scaling, the exhaustive list of references is far too long to be included here (see for instance [START_REF] Boyd | Distributed optimization and statistical learning with the alternating direction method of multipliers[END_REF]).

Spiralling and foldings

Analyzing now the spiralling effect observed in many practical situations where strong monotonicity is not present, it was shown in [START_REF] Lenoir | Accelerating a class of splitting algorithms by iterative foldings[END_REF] an acceleration strategy based on iterative foldings. The idea takes roots in a former study by Lawrence and Spingarn [START_REF] Lawrence | On fixed points of nonexpansive piecewise isometric mappings[END_REF] on folding operators. In the case where T 1 = ∂F and T 2 = ∂C with F and C polyhedral, they introduced a family of averaged operators F, called foldings, which are piecewise isometric, non expansive and positively homogeneous, and showed that the iteration

s t+1 = 1 2 s t + 1 2 F(s t )
produces successive directions which generate linear subspaces, themselves converging to a limit subspace where spiralling is expected. In [START_REF] Lenoir | Accelerating a class of splitting algorithms by iterative foldings[END_REF], the authors used the following hypotheses on the operators :

• The subdifferential operator T 1 = ∂F is proto-differentiable, which means that the graph of the directional variations of T 1 converges graphically to a limit operator, the proto-derivative of T 1 , a weaker notion than differentiability (see [START_REF] Rockafellar | Proto-differentiability of set-valued mappings and its application in optimization[END_REF] for details). In this case, they show how to accelerate convergence by application of averaged sequence of foldings to break spiralling (see Figure 7).

Variable scaling and parameter updates

The high sensitivity to the value of the scaling parameters turns their estimation a difficult issue in practice. So, many researchers have considered the generalization of splitting methods with varying scaling parameters at each iteration. The earliest proposal to our knowledge is a variant of the Forward-Backward splitting proposed by Tseng [START_REF] Tseng | Applications of a splitting algorithm to decomposition in convex programming and variational inequalities[END_REF]. In his analysis, he gave conditions on the scaling parameter to guarantee linear convergence of the algorithm, based on the previous knowledge of the co-coercivity radius of one of the two operators. This is of course difficult to check in practice. The use of adaptive updates of parameters in the splitting algorithms has been early studied by Kontogyorgis and Meyer in [START_REF] Kontogiorgis | A variable-penalty alternating directions method for convex optimization[END_REF] with the additional difficulty of revising the theoretical convergence results. Typically, in the case of a single parameter λ, convergence will be maintained if the sequence {λ t } of parameters converges to a limit value and satisfies +∞ t=0 |λ t+1 -λ t | < +∞ Practical choices for the implementation of these updates appeared in [START_REF] He | Alternating direction method with self-adaptive penalty parameters for monotone variational inequalities[END_REF] for (ADMM) and in [START_REF] Dussault | Separable augmented lagrangian algorithm with multidimensional scaling for monotropic programming[END_REF] for (SALA), taking the following form :

λ t+1 =
θλ t if 0 ≤ t ≤ 100 and t ≡ 0 (mod 10) λ t else where θ ∈ [0. 5 1].

More sophisticated updates that take in consideration the relative behaviour of primal and dual sequences are proposed in [START_REF] Mahey | Adaptive scaling and convergence rates of a separable augmented lagrangian algorithm[END_REF]. Based on the S-Model, their analysis estimates the primal and dual rates by computing τ x = r t+1 r t where r t = i g i (x t i ) is the primal residual at iteration t, and τ u = δ t+1 δ t where δ t = ∇L(x t , u t ) the gradient of the ordinary Lagrangian of the S-Model. The proposed update of the parameter is such that both sequences are kept at a similar pace and is implemented by :

λ t+1 = τ x τ u α λ t
with 0 < α < 1.

Conclusion

We have surveyed the main monotone operator splitting methods and their applications to the decomposition of separable convex problems. This is a still very active research area where recent motivations concerning large-scale problems in signal processing and statistical learning have induced many new adaptations of these relatively old methods which appeared as early as the fifties with Douglas-and Peaceman-Rachford algorithms for linear operators. As these techniques can be interpreted as separable versions of the Augmented Lagrangian dual methods, their main benefit is the regularization effect of the proximal steps which induce numerical stability and smoothness of the implicit primal and dual value functions. We have seen the importance of reformulation to better exploit the decomposition features of each splitting scheme, for example by opposing the S-Model and the M-model. This reveals that the operator splitting techniques are potential candidates for the decomposition of nonconvex problems, even if this has been relatively little explored in the applications (see [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems : proximal algorithms, forwardbackward splitting and regularized gauss-seidel methods[END_REF] for theoretical extensions dealing with semi-algebraic functions). The main drawback which has slowed down their practical use is the difficulty to reach better convergence rates, like superlinear convergence which is unlikely to occur even in the strongly monotone models, neither theoretically nor practically. Nevertheless, we have discussed that question in the last section showing that multidimensional scaling along with adaptive updates of the parameters can significantly improve the speed of convergence in many concrete applications.

Finally, as for most decomposition methods, parallel implementations are natural issues that have been tested by different authors which were not surveyed in the present paper (see Bertsekas and Tsitsiklis [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF] or Eckstein's thesis [START_REF] Eckstein | Splitting methods for monotone operators with applications to parallel optimization[END_REF]).

Figure 1 :

 1 Figure 1: Sequence of iterates for Forward-backward algorithm

Figure 2 :

 2 Figure 2: Douglas-Rachford algorithm

Find

  ξ * = (x * 1 , . . . , x * p ), ζ * = (y * 1 , . . . , y * p ) such that (ξ * , ζ * ) ∈ Graph(T ) A×A ⊥ 3 Application to the decomposition of convex programs Decomposition methods are designed to answer two different objectives :

by substituting ξ i by ξ t+ 1 2 i

 2

Figure 4 :

 4 Figure 4: Primal and dual optimal convergence rates

Figure 5 :

 5 Figure 5: Primal-dual residual rates

Figure 6 :

 6 Figure 6: Multidimensional scaling

Figure 7 :

 7 Figure 7: Breaking spiralling by iterative folding

Acknowledgements : The authors thank J.P. Dussault and M. Minoux for their cooperation and advice, P.L. Combettes, J.C. Pesquet and two anonymous referees for helping to complete the references.

P. Mahey was supported by a grant from PGMO project (Programme Gaspard Monge pour l'Optimisation).