{'fraction_non_alphanumeric': 0.06828041052891454, 'fraction_numerical': 0.021067234102585205, 'mean_word_length': 3.8248440158820194, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Proximal Splitting Algorithms for Convex Optimization: A Tour of Recent Advances, with New Twists

Convex nonsmooth optimization problems, whose solutions live in very high dimensional spaces, have become ubiquitous. To solve them, the class of first-order algorithms known as proximal splitting algorithms is particularly adequate: they consist of simple operations, by handling the terms in the objective function separately. In this overview, we demystify a selection of recent proximal splitting algorithms: we present them within a unified framework, which consists in applying splitting methods for monotone inclusions in primal-dual product spaces, with well-chosen metrics. Along the way, we easily derive new variants of the algorithms and revisit existing convergence results, by extending the parameter ranges in several cases. In particular, we emphasize that when the smooth term in the objective function is quadratic, e.g. for least-squares problems, convergence is guaranteed with larger values of the relaxation parameter than previously known. Such larger values are usually beneficial to the convergence speed in practice.


Introduction
Optimization, also called mathematical programming, is the search of a best object in a set with respect to some criterion. Mathematically, this amounts to characterize a point where a function attains its minimum value, and computationally, this often consists in exhibiting an iterative algorithm, in which a variable converges to such an optimal point. When a function is convex, its minimizers are global: a point which attains a minimum value in its neighborhood is actually optimal in the whole search space. Thus, knowing about the principles behind convex optimization algorithms is important, before considering more general problems. This paper is a constructive and self-contained introduction to the class of proximal splitting algorithms, which are efficient for large-scale convex optimization.
Many problems in statistics, machine learning, signal and image processing, control, and many other fields, can be formulated as convex optimization problems [4,30,100,108,124,125], of the form: for some M ≥ 1, where X is a real Hilbert space, the g m are convex, possibly nonsmooth, functions, and the L m are linear operators. We will see that this generic formulation encompasses the possible presence of constraints. In the age of 'big data', with the explosion in size and complexity of the data to process, it is increasingly important to be able to solve such optimization problems, whose solutions live in very high dimensional spaces [19,64,65,107,123]. For instance, in image processing and computer vision, there are one or several variables to estimate for each pixel, with several millions of pixels, typically. There is an extensive literature about proximal splitting algorithms for solving convex optimization problems, with applications in various fields [8,10,42,45,53,75,85,101,118,123]. They consist of simple, easy to compute, steps that can deal with the terms in the objective function separately: if a function is smooth, its gradient can be used, whereas for a nonsmooth function, the proximity operator will be called instead; hence the wording proximal algorithm. Since the proximity operator of a sum of functions does not have a computable form, in general, it is necessary to handle the proximity operators of the functions, which appear in the problem, individually; this is what splitting refers to, a notion reminiscent to the divide-and-conquer principle, which is ubiquitous in computational sciences. It has been known for several decades that the dual problem 1 , associated to the (primal) optimization problem under consideration, can be easier to solve and, by doing so, may simplify the computation of the primal solution. Yet, the development of primal-dual algorithms, which solve the primal and dual problems simultaneously, in an intertwined way, is more recent. In this tutorial survey, we introduce a selection of primal-dual proximal splitting algorithms, which have been developed over the last decade. We present them in a unified way, which consists in solving a monotone inclusion, expressed in a well-chosen primal-dual product space. Along the way, some new variants are naturally obtained. An important idea for the development of primal-dual algorithms is preconditioning, or equivalently changing the metric of the ambient Hilbert space [13,46,106]. We will make an extensive use of this notion throughout the paper.
In our presentation, we put the emphasis on the possibility to overrelax the algorithms: let us consider an algorithm producing a new estimate z (i+ 1 2 ) of some variable at iteration number i + 1, given the previous estimate z (i) . Instead of setting z (i+1) = z (i+ 1 2 ) before continuing with the next iteration, overrelaxation consists in the extrapolated update z (i+1) = z (i) + ρ (i) (z (i+ 1 2 ) − z (i) ), for some ρ (i) > 1. This step is easy to implement and can improve the convergence speed of the algorithm. Thus, it is important to study the conditions on the parameters ρ (i) , so that convergence is guaranteed. In particular, we show that, if the optimization problem involves a smooth term, which is quadratic, convergence is guaranteed with larger values of the relaxation parameter than previously known. Regularized leastsquares problems are an important class of problems, to which this extension applies.
Also, a desirable property of the proximal splitting algorithms is that they can be parallelized. Thus, when minimizing the sum of a large number of functions on a parallel computing architecture, each computing node can handle one of the functions. Mathematical splitting gracefully blends into the physical splitting of the computing hardware, in that case. The needs for such parallel optimization algorithms are rapidly growing in distributed or federated learning [83,86,92] and artificial intelligence. The last part of the paper is dedicated to parallel versions of the algorithms.
The paper is organized as follows. In section 2, we present the most important notions and principles underlying the class of iterative algorithms for large-scale convex nonsmooth optimization under study. In section 3, we present relaxed versions of the forward-backward splitting algorithm to solve monotone inclusions. This analysis is subsequently applied in section 4 to the Loris-Verhoeven algorithm, which is a primal-dual forward-backward algorithm. In section 5, we analyse the Chambolle-Pock algorithm and its particular case, the Douglas-Rachford algorithm, which is equivalent to the Alternating Direction Method of Multipliers (ADMM). For them too, we derive convergence results for relaxed versions, which extend previously known results. In section 6, we study an algorithm, which we call the Generalized Chambolle-Pock algorithm, since it is applicable to the minimization of two functions composed with two linear operators, and reverts back to the Chambolle-Pock algorithm when one the linear operator is the identity. The algorithm has been presented in the literature as a linearized or preconditioned version of the ADMM, but without relaxation. In section 7, we study an algorithm proposed independently by the first author and by B. C. Vũ, which is also a primal-dual forward-backward algorithm. Section 8 is devoted to recently proposed algorithms based on a three-operator splitting scheme, which can be viewed as the fusion of the forward-backward and Douglas-Rachford two-operator splitting schemes. Along the way, new algorithm variants are proposed, like in subsections 7.1 and 8.1. Finally, in section 9, we propose parallel versions of all these algorithms, by applying them in product spaces.

Notions of convex analysis
In this section, we introduce some notions and notations, which will be used throughout the paper. The reader can find a much more complete account of convex analysis and operator theory in textbooks, e.g. [7,12], or in papers like [40,118]. We consider optimization over real Hilbert spaces: a real Hilbert space is a vector space equipped with a real inner product · , · . When we make use of a norm, without further specification, this is the norm · = · , · 1/2 induced by the inner product. This general formalism allows us to optimize with respect to vectors, matrices, or tensors [70], which are real-valued or complex-valued with some Hermitian properties, but also more complicated objects like quaternions [96]. The Hilbert spaces can be of infinite dimension; this can be useful in control [136] or when solving PDEs, and we will see that most results hold in this general setting. However, there are a few technicalities to keep in mind: 1) Let L : X → U be a linear operator, with X and U two real Hilbert spaces. We define the operator norm of L as L = sup{ Lx : x ∈ X and x ≤ 1}. If L < +∞, L is said to be bounded. L is bounded if and only if it is continuous. In the paper, we will assume all linear operators to be bounded. If X is of finite dimension, L is necessarily bounded.
2) A sequence of points (x (i) ) i∈N in a real Hilbert space X is said to converge weakly to a point x ⋆ ∈ X if, for every y ∈ X , x (i) , y → x ⋆ , y . On the other hand, (x (i) ) i∈N is said to converge strongly to x ⋆ ∈ X if x (i) − x ⋆ → 0. Strong convergence implies weak convergence. If X is of finite dimension, both notions are equivalent and we just say that (x (i) ) i∈N converges to x ⋆ . We will see that most convergence results state weak convergence; this is not a weakness of the proof techniques, and strong convergence in infinite-dimensional spaces cannot be guaranteed without further assumptions [20]. In particular, given a continuous operator T on X , if a sequence (x (i) ) i∈N converges weakly to x ⋆ ∈ X , one cannot deduce that T (x (i) ) i∈N converges weakly to T (x ⋆ ). This makes the convergence proofs of algorithms more difficult to derive.
For the rest of section 2, let X be a real Hilbert space. A subset Ω of X is said to be convex if, for every (x, y) ∈ Ω 2 and a ∈ (0, 1), ax + (1 − a)y ∈ Ω. Note that a convex set can be open, like the interval (0, 1) in R, or closed, like [0, 1]; it can be bounded or unbounded, like [0, +∞). The intersection of convex sets is convex. The union of convex sets is not convex, in general.
A function f : X → R ∪ {+∞} is said to be convex if, for every (x, y) ∈ X 2 and a ∈ (0, 1), We can note that f is allowed to take the value +∞, a distinctive and very useful feature of convex analysis, compared to classical analysis. For calculus, we only need to adopt the following rules: for every t > 0, t(+∞) = +∞, and for every t ∈ R, t + (+∞) = +∞. Note that for every t ∈ R, the intervals (t, +∞] are considered open in R ∪ {+∞}, so that, for instance, the function t ∈ R → (1/t if t > 0, +∞ otherwise) is continuous on R [7, Example 1.22]. The domain of f is the convex set dom f = {x ∈ X : f (x) = +∞}. f is said to be proper if its domain is nonempty. Indicator functions relate the two notions of convexity of a set and of a function: given a set Ω ⊂ X , we define the indicator function of Ω, denoted by ı Ω , as: Then, if Ω is convex, ı Ω is convex; if Ω is nonempty, ı Ω is proper. Also, given two convex sets Ω 1 ⊂ X and Ω 2 ⊂ X , ı Ω 1 ∩Ω 2 = ı Ω 1 + ı Ω 2 .
Indicator functions are important, since they allow to remove constraints and integrate them into the objective function to minimize. Indeed, given a nonempty convex set Ω ⊂ X and a convex function f on X , ≡ minimize We see here the interest of allowing the functions to take the value +∞: it is used to exclude some parts of X from the set of possible solutions. Thus, we will consider the problem of minimizing a sum of convex functions, knowing that thanks to indicator functions, this cover the possible presence of constraints.
A convex optimization problem consists in estimating a minimizer x ⋆ ∈ X , supposed to exist, of a proper convex function f : X → R ∪ {+∞}; that is, a point x ⋆ at which f attains its minimum value min x∈X f (x). The set of minimizers of f , which is convex, is denoted by arg min x∈X f (x). Thus, x ⋆ ∈ arg min x∈X f (x) if and only if, for every x ∈ X , f (x ⋆ ) ≤ f (x). Note that the function t ∈ R → (1/t if t > 0, +∞ otherwise) is proper, convex, bounded from below, but does not have any minimizer.
An optimization algorithm to minimize a proper convex function f : X → R ∪ {+∞}, which has a minimizer, constructs a sequence (x (i) ) i∈N of points, such that f (x (i) ) decreases overall. Thus, for the process to be fruitful, if (x (i) ) i∈N converges to x ⋆ , we must have converges to min f , we cannot deduce that x ⋆ is a minimizer of f . So, we will only consider functions having this lower semicontinuity property. Equivalently, a function f is said to be lower semicontinuous if the convex set {x ∈ X : f (x) ≤ t} is closed, for every t ∈ R. A continuous function is lower semicontinuous. A proper convex lower semicontinuous function is continuous on the interior of its domain [7,Corollary 8.39], so that lower semicontinuity concerns the boundary of the domain, where the function jumps from real values to +∞.
We denote by Γ 0 (X ) the set of convex, proper, lower semicontinuous functions from X to R ∪ {+∞}. We will only consider such functions in optimization problems. If Ω ⊂ X is nonempty and closed, ı Ω ∈ Γ 0 (X ). If f ∈ Γ 0 (X ) and L is a bounded linear operator with values in X , f • L is convex and lower semicontinuous [7,Proposition 9.5], where • denotes the composition of functions; it is proper if and only if dom f ∩ ran L = ∅, where ran L denotes the range of L. Also, if f ∈ Γ 0 (X ) and g ∈ Γ 0 (X ), f + g is convex and lower semicontinuous [7,Corollary 9.4]; it is proper if and only if dom f ∩ dom g = ∅.
The difficulty of convex optimization stems from the fact that even if f and g are two simple convex functions with known minimizers, f + g is difficult to minimize, in general. A well known example is the LASSO problem [129], which consists in minimizing over x ∈ X the function Lx − b 2 + x 1 , where L is a linear operator, b an element, and · 1 is the l 1 norm; note that every norm on X is convex. The minimizer of the least-squares term can be obtained efficiently by solving a linear system, the minimizer of the l 1 norm is the zero element, but there is no straightforward way of finding a minimizer of the sum of the two terms, in general. Algorithms able to solve such problems are precisely the purpose of this paper.
When dealing with vectors, a linear operator can be represented by a matrix, so that the application of a linear operator to a vector can be viewed as a matrix-vector product. But since we place ourselves in the general setting of Hilbert spaces, some definitions are in order. Let L : X → U be a bounded linear operator, for some real Hilbert space U . The adjoint operator L * : U → X of L is the only bounded linear operator such that Lx, u U = x, L * u X , for every x ∈ X and u ∈ U , where we have noted ·, · X and ·, · U the inner products in X and U , respectively, to differentiate them. L is said to be self-adjoint if L * = L. We have L = L * = LL * = L * L . Id denotes the identity operator. Let P be a self-adjoint bounded linear operator on X . P is said to be positive if x, P x ≥ 0, for every x ∈ X , and strongly positive if P − aId is positive, for some real a > 0. If P is positive, for every y ∈ X , the function x → P x + y, x belongs to Γ 0 (X ); for instance, with P = Id and y = 0, we obtain that · 2 ∈ Γ 0 (X ).
In the paper, we consider first-order algorithms; that is, exploiting information about the gradient or subdifferential of the functions. Let us define these notions. A continuous realvalued function f ∈ Γ 0 (X ) is smooth, or differentiable, if at every x ∈ X , there exists an element of X , called the gradient of f at x and denoted by ∇f (x), such that, for every e ∈ X , That is, the affine function y → f (x) + y − x, ∇f (x) is a first-order approximation of f around x. Note that, as a consequence of convexity, this affine function is a minorant of f : Thus, ∇f , the gradient of f , is the operator on X : x → ∇f (x). It is continuous on X [7, Corollary 17.43]. For smooth functions f and g, we have ∇(f +g) = ∇f +∇g, and, for a linear operator L, ∇(f • L) = L * (∇f )L, where AB denotes the composition of the two operators A and B. We can note that the gradient of the quadratic convex function x → P x + y, x , for some positive self-adjoint bounded linear operator P on X and some y ∈ X , is x → 2P x + y. We want to deal with nonsmooth functions; for instance, the absolute value function f : R → R : t → |t|, and by extension the l 1 norm, are convex but not differentiable at 0. Moreover, we allow convex functions to take the value +∞. Thus, we need a more general notion, capturing first-order information, than the gradient. This is where the subdifferential comes to the rescue. Let f ∈ Γ 0 (X ). For every x ∈ X , we define the subdifferential of f at x, denoted by ∂f (x), as the closed convex set [7,Proposition 16.4] ∂f (x) = u ∈ X : ∀y ∈ X , f (x) + y − x, u ≤ f (y) .
The elements of ∂f (x) are called subgradients of f at x. A subgradient u ∈ ∂f (x) can be viewed as the gradient of the affine minorant y → f (x) + y − x, u of f . If f is differentiable at x, there is only one subgradient, which is the gradient: It is important to note that, for any g ∈ Γ 0 (X ) and x ∈ X [7, Proposition 16 where the sum of the two sets is the Minkowski sum, but the inclusion may be strict [7,Remark 16.8]. Similarly, given a bounded linear operator L, for every x ∈ X , We denote by 2 Ω the set of all subsets of a set Ω. Then the subdifferential of f is the set-valued operator ∂f : X → 2 X , x → ∂f (x). We assimilate a set-valued operator M : X → 2 X with exactly one element in each set M (x) as a single-valued operator M : X → X , and conversely. For instance, for a smooth convex function f , we have ∂f = ∇f .
We will consider convex optimization problems of the form: where M ≥ 1, each L m is a bounded linear operator from X to some real Hilbert space U m , possibly the identity if U m = X , and each g m ∈ Γ 0 (U m ). The well known Fermat's rule states that x ⋆ ∈ X is a solution to (11) if and only if This fundamental property is easy to demonstrate: let f ∈ Γ 0 (X ) and x ⋆ ∈ X . Then As already mentioned in (10), for every x ∈ X , we have but the inclusion may be strict. All algorithms considered in the paper actually find a solution if such a solution exists. But in view of (16), the existence of a minimizer of m g m • L m is not sufficient to guarantee the existence of a solution to (17). Therefore, throughout the paper, we assume that for every optimization problem considered, the solution set of the corresponding monotone inclusion (17) is nonempty. Given (16) and Fermat's rule, a solution x ⋆ of (17) is a solution to the optimization problem (11). If the functions satisfy a so-called qualification constraint, we have an equality instead of an inclusion in (16); in that case, it is sufficient to show that a solution to (11) exists. A qualification constraint is [43,Proposition 4.3]: where, for any linear operator L : X → U and Ω ⊂ U , L −1 (Ω) = {x ∈ X : Lx ∈ Ω}, and sri Ω denotes the strong relative interior of a convex set Ω, which is the set of points u ∈ Ω such that the cone generated by −u + Ω is a closed linear subspace. Other qualification constraints can be found in [43,Proposition 4.3]. Note that for M m=1 g m • L m to be proper, it is necessary to have ∩ M m=1 L −1 m (dom g m ) = ∅. Since the strong relative interior is even larger than the interior of a set, the condition (18) is really not a lot to ask.
We have seen that the subdifferential is an important notion, in particular with Fermat's rule. But the algorithms considered in the paper to solve optimization problems involving nonsmooth functions will not make use of subgradients directly. Instead, they will make calls to the proximity operators of the functions: we define the proximity operator [97] of any function f ∈ Γ 0 (X ) as By strong convexity of f + 1 2 · −x 2 , its minimizer always exists and is unique, so that the proximity operator is well defined. For any nonempty closed convex set Ω ⊂ X and τ > 0, is the projection onto Ω, which does not depend on τ . Thus, proximal algorithms can be viewed as generalizing algorithms to find points in intersections of convex sets [38]. Although the definition of the proximity operator is implicit, it has a closed form for many functions of practical interest. For instance, for the absolute value (and by extension for the l 1 norm by elementwise application), this is soft-thresholding: we have, for any τ > 0, There are fast and exact methods to compute the proximity operators for a large class of functions [32,42,48,50,63,109], see also the website http://proximity-operator.net. In this paper, we consider proximal splitting algorithms, which are iterative algorithms producing an estimate of a solution to a convex optimization problem of the form (11), by calling at each iteration the proximity operator or the gradient of each function g m , as well as the operators L m and L * m . An optimization problem can often be written under the form (11) in several equivalent but different ways, with each formulation having appropriate algorithms to solve it. For instance, when a function ı Ω for some convex set Ω appears in the problem, it is implicitly assumed, in general, that one knows how to project onto Ω. For instance, to minimize f ∈ Γ 0 (X ) subject to the constraint L · = b, for some linear operator L and element b, the two following formulations are equivalent: The first formulation is appropriate if one knows how to project on the affine constraint set efficiently, whereas the second formulation is preferable if one does not, and asks for an algorithm making calls to L and L * .
For every f ∈ Γ 0 (X ), we denote by f * the conjugate of f , defined by which belongs to Γ 0 (X ) as well. We have f * * = f . An important property is the Moreau identity, which allows us to compute the proximity operator of f * from the one of f , and conversely: for every τ > 0, We do not define here further notions related to duality, which will be introduced progressively along the presentation of primal-dual algorithms.

Notions of operator theory and fixed-point algorithms
Let M : X → 2 X be a set-valued operator; as mentioned, single-valued operators can be viewed as set-valued operators, so that the subsequent notions apply to them as well. The graph of M is the set {(x, v) ∈ X 2 : v ∈ M x} (we may use the shortened notation M x for M (x), usually adopted for linear operators, with any operator). The inverse operator of M , denoted by M −1 , is the set-valued operator with graph For instance, for any bounded linear operator L, L −1 is well defined everywhere, in this sense. But L −1 x can be the empty set, if x / ∈ ran L, or can contain several elements: L −1 0 is the kernel of L. The inverse relates the subdifferentials of a function f ∈ Γ 0 (X ) and its conjugate f * ∈ Γ 0 (X ) defined in (24), via the simple relation [7,Corollary 16.30]: M is said to be monotone if, for every (x, We can note that for a self-adjoint bounded linear operator, being positive is equivalent to being monotone. A monotone operator M is said to be maximally monotone if there exists no monotone operator whose graph strictly contains the graph of M . For every f ∈ Γ 0 (X ), ∂f is maximally monotone (28) [7, Theorem 20.25]. The first derivative of a smooth convex function f : R → R is nondecreasing; so, monotonicity can be viewed as a generalization of being nondecreasing to multidimensional and set-valued operators. The problem of finding x ⋆ such that 0 ∈ M x ⋆ is called a monotone inclusion. Thus, as we have seen with Fermat's rule (11)- (12), convex optimization problems consist in solving monotone inclusions. That is why monotone operator theory plays an essential role in optimization. If M is monotone, we define its resolvent, denoted by J M , as the operator J M is a well-defined single-valued operator on X . Now, let us apply Fermat's rule to the definition of the proximity operator in (19): if y = prox f (x), for any f ∈ Γ 0 (X ) and x ∈ X , then 0 ∈ ∂f (y) + y − x. Equivalently, y = (∂f + Id) −1 (x). That is, Thus, proximal algorithms handle nonsmooth functions by calls to their proximity operators; this can be viewed as a regularized and implicit way of dealing with their set-valued subdifferentials. But this comes at a price: prox f +g cannot be expressed from prox f and prox g , in general; in rare cases, prox f +g = prox f • prox g [109]. And prox f •L does not have a computable form either, in general. For instance, the proximity operator of the least-squares term τ 2 L · −b 2 , for some bounded linear operator L : X → U , b ∈ U and τ > 0, is which has a closed form, but amounts in practice to solve a linear system, which is generally too costly. Thus, even the apparently simple problem of minimizing the sum of two nonsmooth functions is in fact difficult. Thus, splitting is necessary: we have to call the proximity operator of each simple-enough function individually, when minimizing a sum of functions, and to deal with the linear operators separately.
Let T : X → X be a single-valued operator. A fixed point of T is any x ⋆ ∈ X satisfying T is said to be β-Lipschitz continuous, for some real β ≥ 0, if, for every (x, x ′ ) ∈ X 2 , We can note that a bounded linear operator L is L -Lipschitz continuous. T is said to be nonexpansive if it is 1-Lipschitz continuous. T is said to be α-averaged, for some α ∈ (0, 1), if there exists a nonexpansive operator T ′ : X → X such that If T is α-averaged, for some α ∈ (0, 1), it is α ′ -averaged, for every α ′ ∈ (α, 1), and nonexpansive; in addition, for every ρ ∈ (0, 1 α ), the relaxed operator is ρα-averaged. Firmly nonexpansive means (1/2)-averaged. Finally, T is said to be ξcocoercive, for some ξ > 0, if ξT is firmly nonexpansive; or, equivalently [7,Proposition 4.4 Let Z be a real Hilbert space and T be a nonexpansive operator on Z with nonempty set of fixed points. A fixed-point algorithm, initialized with some z (0) ∈ Z, consists in iterating where i is the iteration counter and z (i) denotes the version of the variable z produced during the i-th iteration. Note that we have changed the name of the variable from x to z, to emphasize that z (i) may not be the variable converging to a solution of the optimization problem at hand: it can be a concatenation of variables, like z (i) = (x (i) , u (i) ), with x (i) converging to a desired solution and u (i) an auxiliary variable, the value of which is not of primary concern. Or it can be that there is an internal variable x (i) , among the set of operations modeled by the abstract operator T , which converges to a solution, whereas z (i) is not interesting in itself. Then, by combining (32) and (33), we have, for every fixed point z ⋆ of T and i ≥ 0, Thus, the variable z (i) gets closer and closer, at every iteration, to the set of fixed points of T . But this is not sufficient to obtain convergence of the sequence (z (i) ) i∈N to some fixed point. A counter-example is the rotation operator of angle θ ∈ [0, 2π) in R 2 , represented by the matrix Then the inequality in (38) is an equality, so that (z (i) ) i∈N generated by (37) turns over and over again in the plane, and does not converge to the fixed point (0, 0) of T . Thus, we have to be slightly more strict in our assumptions: we will suppose that T in (37) is α-averaged, for some α ∈ (0, 1), and not only nonexpansive. And in view of (35), given some sequence (ρ (i) ) i∈N of relaxation parameters, we will consider the relaxed iteration Note that with ρ (i) = 1, there is no relaxation and we recover the iteration (37). While there is no interest in doing underrelaxation with ρ (i) less than 1, overrelaxation with ρ (i) larger than 1 may be beneficial to the convergence speed, as observed in practice in many cases. The idea of using overrelaxation comes from the following intuition: if the sequence (z (i) ) i∈N generated by (37) converges to some fixed point z ⋆ , this means that z (i+1) = T z (i) tends to be closer to z ⋆ than z (i) , on average. Hence, we may want, starting at z (i) , to move further in the direction T z (i) − z (i) , which improves the estimate. This is exactly what the relaxed iteration (40) does, if ρ (i) > 1. Now, we can state the fundamental result, on which we will rely to establish convergence of optimization algorithms: Krasnosel'skiȋ-Mann theorem (algorithm (40)) [7,Proposition 5.16]. Let Z be a real Hilbert space and let T : Z → Z be an α-averaged operator, for some α ∈ (0, 1), with a nonempty set of fixed points. Let (ρ (i) ) i∈N be a sequence in [0, 1/α], such that i∈N ρ (i) (1 − αρ (i) ) = +∞. Let z (0) ∈ Z. Then the sequence (z (i) ) i∈N defined by the iteration (40) converges weakly to some fixed point z ⋆ of T , and (z (i+ 1 2 ) − z (i) ) i∈N converges strongly to 0.
We can note that in this theorem, as a direct implication, (z (i+ 1 2 ) ) i∈N converges weakly to the same fixed point z ⋆ . In practice, one might prefer considering z (i+ 1 2 ) instead of z (i+1) as an estimate of the solution, because it is in the range of T .
In this theorem, the assumption that T has a nonempty set of fixed points is important: consider for instance T = Id + c, for some c = 0, which is α-averaged, for every α ∈ (0, 1). Then T has no fixed point and the iteration (40) diverges.
Finally, we can note that ρ (i) = 1 is allowed in the theorem. In practice, we recommend to set ρ (i) = 1 when running an algorithm for the first time, before trying larger values.
All algorithms presented in the paper are instances of the simple form (40), and their convergence analysis relies on the Krasnosel'skiȋ-Mann theorem.
There remains to see how to construct averaged operators, whose fixed points are solutions to optimization problems. To that aim, we will combine gradients and proximity operators of functions, like Lego bricks. Indeed, the set of averaged operators is stable by convex combinations (a convex combination is a linear combination with nonnegative weights whose sum is one) and compositions: [99] [7, Proposition 4.44] (composition of averaged operators) Let T 1 : X → X be a α 1 -averaged operator and T 2 : X → X be a α 2 -averaged operator, for some α 1 ∈ (0, 1) and α 2 ∈ (0, 1). Set and T = T 2 T 1 . Then α ∈ (0, 1) and T is α-averaged.
Moreover, if M : X → 2 X is a maximally monotone operator, J M is firmly nonexpansive (42) [7, Proposition 23.8]. And since the proximity operator is the resolvent of the subdifferential, it is firmly nonexpansive as well. Finally, in the paper, when we deal with a smooth function using calls to its gradient, we will assume that the gradient is Lipschitz continuous. Indeed, the Baillon-Haddad theorem [7, Corollary 18.17] states the following: let f : X → R be a differentiable convex function, whose gradient ∇f is β-Lipschitz continuous, for some β > 0. Then We now have all the ingredients to start designing algorithms. Let us start with the minimization of a differentiable convex function f : X → R, whose gradient ∇f is β-Lipschitz continuous, for some β > 0, and whose set of minimizers is nonempty. Gradient descent, initialized at some x (0) ∈ X , consists in iterating, for i ≥ 0, for some parameter γ ∈ (0, 2 β ). It is called a descent method, because one can show that f (x (i+1) ) ≤ f (x (i) ), so that the algorithm "goes down" along the function values. To analyze this algorithm, we first observe that x (i+1) = T x (i) , for the operator T = Id − γ∇f . Since 1 β ∇f is firmly nonexpansive, it is easy to show that T is α-averaged, with α = γβ 2 ∈ (0, 1). Moreover, x ⋆ ∈ X is a fixed point of T if and only if ∇f (x ⋆ ) = 0, which is equivalent to x ⋆ being a minimizer of f . Hence, we can invoke the Krasnosel'skiȋ-Mann theorem (without relaxation) to get weak convergence of the sequence (x (i) ) i∈N to a minimizer x ⋆ of f .
We can note that the relaxed version (40) with a fixed parameter ρ of gradient descent simplifies to x (i+1) = x (i) − ργ∇f (x (i) ), with ργ ∈ (0, 2 β ), so that we can just set ρ = 1 and keep γ as the only parameter. How to choose the value of γ then? As γ tends to 2 β , the operator Id − γ∇f is α-averaged, with α tending to 1, and we have seen that convergence is not guaranteed with merely nonexpansive operators. So, we can expect convergence to become slower as γ tends to 2 β ; but this is true in the worst case: if f = 1 2 · 2 , whose gradient ∇f = Id is β = 1-Lipschitz continuous, we have x (i+1) = (1 − γ)x (i) , and γ = 1 β = 1 is the best choice, whereas convergence becomes slower and slower as γ → 2 β = 2. However, a more insightful and realistic example is the least-squares function f : x → 1 2 Ax − b 2 , for some bounded linear operator A with A = 1 and element b. Gradient descent is then equivalent to the classical Richardson iteration to solve the linear system A * Ax = A * b: for some γ ∈ (0, 2). The operator mapping x (i) to x (i+1) is γ/2-averaged; but it is also affine and quasicontractive, and the convergence rate is fully characterized by the eigenvalues of A * A. Let λ min > 0 be its smallest positive eigenvalue; we have assumed that the largest eigenvalue is 1. Then the best convergence rate is obtained for γ = 2/(1 + λ min ), which, for λ min small, is close to 2. In many applications, the spectrum of A * A is concentrated around 0, so that λ min , even if unknown, is likely to be small. In any case, there is no interest in choosing γ < 1. So, as a rule of thumb, in gradient descent, one should try the two values γ = 1 β and γ = 1.9 β . Similarly, in the Krasnosel'skiȋ-Mann theorem, we recommend tuning the parameters in the algorithms with ρ = 1, and then trying overrelaxation with ρ = 1.9 2α . Let us now consider the minimization of a function f ∈ Γ 0 (X ), whose set of minimizers is nonempty. The proximal point algorithm [114], initialized at some x (0) ∈ X , consists in iterating, for i ≥ 0, for some γ > 0. Since, by definition of the proximity operator, , this is also a descent algorithm. We can easily show that so that the proximal point algorithm looks similar to gradient descent, but with an implicit call to the subdifferential via the proximity operator. prox γf is firmly nonexpansive, and x ⋆ ∈ X is a fixed point of prox γf if and only if 0 ∈ ∂f (x ⋆ ), which is equivalent to x ⋆ being a minimizer of f . So, we can again invoke the Krasnosel'skiȋ-Mann theorem (without relaxation) to get weak convergence of the sequence (x (i) ) i∈N to a minimizer x ⋆ of f . Note that if the proximal point algorithm can be applied with every γ > 0, by taking γ very large, it essentially converges in one single iteration. So, it is rare to be able to minimize a function, whose minimizer is unknown, using the proximal point algorithm as is. However, we will see in the rest of the paper that proximal splitting algorithms are able to minimize a sum of functions, using calls to their gradients or proximity operators.

The forward-backward algorithm
Let Z be a real Hilbert space. Let M : Z → 2 Z be a set-valued maximally monotone operator and let C : Z → Z be a ξ-cocoercive operator, for some real ξ > 0. Let us consider the monotone inclusion whose solution set is supposed nonempty. Let z (0) ∈ Z be some initial estimate of a solution and let γ > 0 be some real parameter. The classical forward-backward iteration, to find a solution of (48), consists in iterating This method, proposed by Mercier [95], was further developed by many authors [33,37,47,68,88,102,130].
A not-so-well known extension of the forward-backward iteration consists in relaxing it. Let (ρ (i) ) i∈N be a sequence of relaxation parameters. The iteration becomes: Forward-backward iteration for (48): for i = 0, 1, . . .
To the best of our knowledge, Theorem 25.8 of [6] is the first convergence result about the relaxed forward-backward algorithm, with a smaller relaxation range than in Lemma 3.1 below, though.
We can remark that the explicit mapping from z (i) to z (i+ 1 2 ) in (50) can be equivalently written under the implicit form: The now standard convergence result for the forward-backward iteration is the following [49,Lemma 4.4] [7, Theorem 26.14]: We can allow γ = 2ξ in Lemma 3.1, but in this case δ = 1 and we cannot set ρ (i) = 1.
Since we are interested in overrelaxation, we write all theorems such that the default choice ρ (i) = 1 is allowed.
We can note that in this algorithm, like in every algorithm of the paper, the sequence (z (i+ 1 2 ) ) i∈N converges weakly to the same solution as (z (i) ) i∈N . For any i, z (i+ 1 2 ) may be a better estimate of the solution than z (i) , since it is in the domain of M ; that is, M z (i+ 1 2 ) = ∅. We can also remark that if γ is close to 2ξ, δ is close to 1, so that overrelaxation cannot be used. This explains why the relaxed forward-backward iteration is not so well known. Now, let P be a self-adjoint, strongly positive, bounded linear operator on Z. Clearly, solving (48) is equivalent to solving Let Z P be the Hilbert space obtained by endowing Z with the inner product (x, Proposition 20.24]. However, the cocoercivity of P −1 C in Z P has to be checked on a case-by-case basis. The preconditioned forward-backward iteration to solve (52) is: The corresponding convergence result follows: Let z (0) ∈ Z and let (ρ (i) ) i∈N be a sequence in [0, δ] such that i∈N ρ (i) (δ − ρ (i) ) = +∞. Then the sequence (z (i) ) i∈N defined by the iteration (53) converges weakly to a solution of (48).
Proof. This is a straightforward application of Lemma 3.1 in Z P instead of Z, with γ = 1. Weak convergence in Z P is equivalent to weak convergence in Z, and the solution sets of (48) and (52) are the same.
We can remark that the explicit mapping from z (i) to z (i+ 1 2 ) in (53) can be equivalently written under the implicit form: Also, we can write the iteration as z ( , an iteration studied with nonlinear kernels P in [71] and [22].
If C = 0, the forward-backward iteration reduces to the proximal point algorithm [7,114]; we have discussed this algorithm in section 2 with proximity operators, but it generalizes to any resolvent. Let us give a formal convergence statement for the preconditioned proximal point algorithm. Let M : Z → 2 Z be a maximally monotone operator. The problem is to solve whose solution set is supposed nonempty. Let z (0) ∈ Z and let (ρ (i) ) i∈N be a sequence of relaxation parameters. The relaxed and preconditioned proximal point algorithm consists in the iteration: Proximal point iteration for (55): for i = 0, 1, . . .
The mapping from z (i) to z (i+ 1 2 ) in (56) can be equivalently written as The convergence of the proximal point algorithm can be stated as follows: (56)). Let z (0) ∈ Z and let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2 − ρ (i) ) = +∞. Then the sequence (z (i) ) i∈N defined by the iteration (56) converges weakly to a solution of (55).

The case where C is affine
Let us continue with our analysis of the forward-backward iteration to solve (48). In this section, we suppose that, in addition to being ξ-cocoercive, C is affine; that is, for some self-adjoint, positive, nonzero, bounded linear operator Q on Z and some element c ∈ Z. We have ξ = 1/ Q . Now, we can write (51) as where Therefore, the forward-backward iteration (50) can be interpreted as a preconditioned proximal point iteration (56), applied to find a zero of M + C! Since P must be strongly positive, we must have γ ∈ (0, ξ), so that the admissible range for γ is halved. But in return, we get the larger range (0, 2) for relaxation. Hence, we have the following new convergence result: Then the sequence (z (i) ) i∈N defined by the iteration (50) converges weakly to a solution of (48).
Proof. In view of (57) and (59), this is Theorem 3.3 applied to the problem (55) with M + C as the maximally monotone operator.
Now, let us improve the range (0, ξ) for γ in Theorem 3.4 into (0, ξ]. Let γ ∈ (0, ξ]. We first recall that for every self-adjoint, positive, bounded linear operator P on Z, there exists a unique self-adjoint, positive, bounded linear operator on Z, denoted by √ P and called the positive square root of P , such that √ P √ P = P . So, let A be a bounded linear operator from Z to some real Hilbert space A, such that A * A = Q. A = √ Q is a valid choice, but in some cases, Q = A * A for some A in the first place, for instance in least-squares problems, see eq. (74). In any case, we do not need to exhibit A, the fact that it exists is sufficient here. Furthermore, let B be a bounded linear operator from some real Hilbert space Again, we do not have to exhibit B, the fact that it exists is sufficient. Then we can rewrite the problem (48) as finding a pair w = (z, b) ∈ W = Z × B solution to the monotone inclusion The operators S and N defined in (61) are maximally monotone in W. To solve a monotone inclusion involving two monotone operators, the Douglas-Rachford splitting algorithm [37,60,88,127] is a natural choice. So, let us consider a general monotone inclusion in some arbitrary real Hilbert space W, with nonempty solution set: for some maximally monotone operators S and N on W. Let w (0) ∈ W and y (0) ∈ W, let τ > 0 be a real parameter and let (ρ (i) ) i∈N be a sequence of relaxation parameters. The Douglas-Rachford iteration can be written in several ways; one of them is the following: We have the following convergence theorem [7, Theorem 26.11]: Lemma 3.5 (Douglas-Rachford algorithm (63)). Let w (0) ∈ W and y (0) ∈ W, let τ > 0 and let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2 − ρ (i) ) = +∞. Then the sequence (w (i) ) i∈N defined by the iteration (63) converges weakly to a solution of (62).
Thus, let us apply the Douglas-Rachford algorithm to (61), in the lifted space W = Z × B, with τ = γ. For this, we need to compute the resolvents; we have [98, eq. 15]: Hence, assuming that w (0) = (z (0) , 0) for some z (0) ∈ Z and subsequently keeping only the first part of the variable w as the variable z, the Douglas-Rachford iteration in this setting is: Suppose that (y So, we can remove all variables but z and we recover the forward-backward iteration (50)! It is remarkable that the forward-backward algorithm can be viewed as an instance of the Douglas-Rachford algorithm, since they are both fundamental algorithms to find a zero of a sum of two monotone operators. This is the case only because the cocoercive operator C is supposed here to be affine, though. Hence, as an application of Lemma 3.5, we have the following convergence result, which extends Theorem 3.4 with the possibility of setting γ = ξ: Theorem 3.6 (forward-backward algorithm (50), affine case). Let z (0) ∈ Z, let γ ∈ (0, ξ] and let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2−ρ (i) ) = +∞. Then the sequence (z (i) ) i∈N defined by the iteration (50) converges weakly to a solution of (48).
Thus, if γ = ξ, we can set ρ (i) = 1.4 according to Lemma 3.1; with Theorem 3.6, we can do better and set ρ (i) = 1.9. So, which values of the parameters should be used in practice? As discussed in section 2, the authors' recommendation for a given practical problem is to try the three settings: • γ = ξ and ρ (i) = 1, • γ = 1.9ξ and ρ (i) = 1, • γ = ξ and ρ (i) = 1.9. We illustrate these three choices within an image processing experiment in Figure 1.

Applications to convex optimization
Let X be a real Hilbert space. Let f ∈ Γ 0 (X ) and let h : X → R be a convex and differentiable function, whose gradient ∇h is β-Lipschitz continuous, for some real β > 0. We consider the convex optimization problem: minimize From Fermat's rule, the problem (71) is equivalent to (48) with M = ∂f , which is maximally monotone, and C = ∇h, which is ξ-cocoercive, with ξ = 1/β (and as mentioned in section 2, a zero of ∂f + ∇h is supposed to exist; equivalently here, a minimizer of f + h is supposed to exist). So, it is natural to use the forward-backward iteration to solve (71). Let γ > 0, let x (0) ∈ X and let (ρ (i) ) i∈N be a sequence of relaxation parameters. The forward-backward iteration, also called proximal gradient algorithm, is: Forward-backward iteration for (71): for i = 0, 1, . . .
Now, let us focus on the case where h is quadratic: for some self-adjoint, positive, nonzero, bounded linear operator Q on X , some c ∈ X , and some t ∈ R. We thereafter omit the constant t in the notations of quadratic functions, since it does not play any role, so we can assume that t = 0. A very common example is a least-squares penalty, in particular to solve inverse problems [30]; that is, for some bounded linear operator A from X to some real Hilbert space Y and some y ∈ Y. Clearly, (74) is an instance of (73), with Q = A * A and c = −A * y.
We have, for every x ∈ X , So, β = Q . Suppose that γ ∈ (0, 1/β) and set P = 1 γ Id − Q, which is strongly positive. Then the update step in (72) can be written as where we introduce the norm · P : x → x, P x . So, x (i+ 1 2 ) can be viewed as the result of applying the proximity operator of f + h with the preconditioned norm · P . Thus, convergence follows from Theorem 3.4. Furthermore, as a direct consequence of Theorem 3.6, we have: Suppose that h is quadratic. Let x (0) ∈ X , let γ ∈ (0, 1/β] and let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2 − ρ (i) ) = +∞. Then the sequence (x (i) ) i∈N defined by the iteration (72) converges weakly to a solution of (71).

The Loris-Verhoeven algorithm
Let X and U be two real Hilbert spaces. Let g ∈ Γ 0 (U ) and let h : X → R be a convex and differentiable function with β-Lipschitz continuous gradient, for some real β > 0. Let L : X → U be a bounded linear operator.
Often, the template problem (71) of minimizing the sum of two functions is too simple and we would like, instead, to minimize We assume that there is no simple way to compute the proximity operator of g • L. The corresponding monotone inclusion, which we will actually solve, is To bypass the annoying operator L, we introduce an auxiliary variable u ∈ ∂g(Lx), which shall be called the dual variable, so that the problem now consists in finding x ∈ X and u ∈ U such that u ∈ ∂g(Lx) The interest in increasing the dimension of the problem is that we obtain a system of two monotone inclusions, which are decoupled: ∂g and ∇h appear separately in the two inclusions. So, equivalently, the problem is to find a pair of objects The operator M : is natural to think of the forward-backward iteration to solve the problem (80). However, to make the resolvent of M computable with the proximity operator of g, we need preconditioning. The solution consists in the iteration, that we first write in implicit form: where τ > 0 and σ > 0 are two real parameters, ). It is not straightforward to see that this yields an explicit iteration. The key is to remark that we have so that we can update the primal variable x (i+ 1 2 ) , once the dual variable u (i+ 1 2 ) is available. So, the first step of the algorithm is to construct u (i+ 1 2 ) . It depends on Lx (i+ 1 2 ) , which is not yet available, but using (82), we can express it using x (i) and LL * u (i+ 1 2 ) . This last term is canceled in the preconditioner P to make the update explicit. That is, We can note that (79) is equivalent to which implies that 0 ∈ (∂g) −1 u − L(∇h) −1 (−L * u); this is the first-order characterization of the convex optimization problem which is called the dual problem to the primal problem (77). So, if a pair (x, u) ∈ X × U is a solution to (80), then x is a solution to (77) and u is a solution to (89). Thus, let τ > 0 and σ > 0, let x (0) ∈ X and u (0) ∈ U , and let (ρ (i) ) i∈N be a sequence of relaxation parameters. The primal-dual forward-backward iteration, which we call the Loris-Verhoeven iteration is: Loris-Verhoeven iteration for (77) and (89): This algorithm was first proposed by Loris and Verhoeven, in the case where h is a leastsquares term [90]. It was then rediscovered several times and named Primal-Dual Fixed-Point algorithm based on the Proximity Operator (PDFP2O) [34] or Proximal Alternating Predictor-Corrector (PAPC) algorithm [57]. The above interpretation of the algorithm as a primal-dual forward-backward iteration has been presented in [41]. We can note that the Loris-Verhoeven iteration can be written in such a way that there is only one call to ∇h and L * per iteration, see the form (206) for instance.
Proof. In view of (81) and (54), this is Theorem 3.2 applied to the problem (80). For this, P must be strongly positive, which is the case if and only if στ The following result makes it possible to have στ L 2 = 1; it is a consequence of the analysis by O'Connor and Vandenberghe [98] of the PD3O algorithm [134], discussed in section 8, of which the Loris-Verhoeven algorithm is a particular case. See also Theorems 3.4 and 3.5 in [34] for the same result, but without relaxation.  (90)). Suppose that X and U are of finite dimension. Let x (0) ∈ X and u (0) ∈ U . Let τ ∈ (0, 2/β) and σ > 0 be such that στ L 2 ≤ 1.
Then the sequences (x (i) ) i∈N and (u (i) ) i∈N defined by the iteration (90) converge to a solution of (77) and to a solution of (89), respectively.
Note that the finite dimension assumption is not necessary in the proof technique [98], so it could be removed.
Thus, in practice, one can keep τ as the only parameter to tune and set σ = 1/(τ L 2 ). For primal-dual algorithms, whose iteration satisfies an inclusion like (81), we want P to be as close to 0 as possible, broadly speaking. We also want the proximal parameters, here σ and τ , to be as large as possible. But these objectives are antagonistic: if τ is large, σ must be small. Nevertheless, for a given value of τ , setting σ to the largest possible value 1/(τ L 2 ) is the best we can do.
If X = U , L = Id, and we set σ = 1/τ , the Loris-Verhoeven iteration becomes: So, we can discard the dual variable and we recover the forward-backward iteration (72). It is interesting that in this primal algorithm, there is an implicit dual variable u ( )/τ , which converges to a solution of the dual problem; that is, to a minimizer of g * (u) + h * (−u).
Again, let us focus on the case where h is quadratic; that is, h : x → 1 2 x, Qx + x, c , for some self-adjoint, positive, nonzero, bounded linear operator Q on X , and some c ∈ X . We have β = Q . We can rewrite the primal-dual inclusion (81), which characterizes the Loris-Verhoeven iteration (90), as We recognize the structure of a primal-dual proximal point algorithm, so we can apply Theorem 3.3. We can do better and apply Theorem 3.6 to obtain convergence under the conditions τ ∈ (0, 1/β] and στ L 2 < 1. Indeed, we want to solve 0 ∈ P −1 M z + P −1 Cz in Z P , with P defined in (81). We have P −1 C : (x, u) → (τ Qx + τ c, 0), which is affine and 1/(τ β)-cocoercive in Z P . Hence, Theorem 3.6 can be applied with γ = 1.
In Figure 1, we illustrate how overrelaxation can accelerate the Loris-Verhoeven algorithm, when h is quadratic. We consider the problem of image deblurring, which consists in restoring an image y corrupted by Gaussian blur and Gaussian noise, shown in Figure 1 (b). We formulate deblurring as an optimization problem of the form (77); that is, we aim at minimizing Ψ = h + g • L, where h : x → 1 2 Ax − y 2 is a least-squares data-fidelity term, with A the convolution with a lowpass filter (β = A = 1), and g•L is the 'isotropic' total variation [115]: L maps the image to the vector field made of all horizontal and vertical finite differences, and g is λ times the l 1,2 norm, see [30,51] for details, where λ = 0.002 is the regularization parameter. So, we use the Loris-Verhoeven algorithm to solve the problem; a solution x ⋆ (which seems to be unique and is estimated by x (i) with i = 10 6 ) is shown in Figure 1 (c). We can see from the plots of x (i) − x ⋆ 2 and Ψ(x (i) ) that the choice τ = 1 β , στ L 2 = 1 and ρ (i) = 1.9, as allowed by Theorem 4.3, speeds up the convergence, in comparison with standard choices allowed by Theorem 4.2.

The Primal-Dual Fixed Point (PDFP) algorithm
Let us extend the problem (77) to using the same assumptions as before in this section, with an additional function f ∈ Γ 0 (X ). The dual problem is minimize  Can we modify the Loris-Verhoeven algorithm to handle this extended problem? We will see one way to do this in section 8, with the PD3O algorithm. Meanwhile, let us extend the Loris-Verhoeven algorithm, by keeping its primal-dual forward-backward structure. For this, we suppose that, for any γ > 0, the proximity operator of γf is affine; that is, there exist some self-adjoint, positive, bounded linear operator R γ on X , with R γ ≤ 1, and some element a γ ∈ X , such that We can mention two cases of practical interest, where the proximity operator is affine. First, this is the case if f is quadratic: f : x ∈ X → 1 2 x, Q ′ x + x, c ′ + t ′ , for some self-adjoint, positive, bounded linear operator Q ′ , some c ′ ∈ X , and some t ′ ∈ R, with R γ = (γQ ′ + Id) −1 and a γ = −(γQ ′ + Id) −1 c ′ . Second, this is the case if f is the indicator function of an affine subspace: f : x → (0 if x ∈ A, +∞ otherwise), where A is a closed affine subspace of X ; for instance, A = {x ∈ X : Ax = y}, for some bounded linear operator A and some element y in the range of A; then prox γf is the projector onto A, which does not depend on γ.
This corresponds to the iteration: PDFP iteration for (93) and (94): This algorithm has been analyzed by Chen et al. [35] and called the Primal-Dual Fixed Point (PDFP) algorithm, so we give it this name.
It remains an open question whether the PDFP algorithm can be relaxed like in Theorem 4.4, for an arbitrary f .
We can check that when h = β 2 · 2 , it is equivalent to apply the Loris-Verhoeven algorithm to minimize g • L + h or to apply the PDFP to minimize f + g • L, with f = β 2 · 2 . Indeed, in the later case, let τ ′ > 0 and σ > 0 be the parameters of the PDFP algorithm. Set τ = τ ′ 1+τ ′ β .
We have prox τ ′ f : x ∈ X → 1 1+τ ′ β x, which is linear. Then the PDFP algorithm (97) becomes: and we recover the Loris-Verhoeven algorithm (90) applied to minimize g • L + h, with parameters τ and σ. The conditions for convergence from the PDFP interpretation are τ < 1/β, στ L 2 < 1, and δ = 2, so that we do not gain anything in comparison with Theorem 4.3.

The Chambolle-Pock and the Douglas-Rachford algorithms
We now consider a similar problem to (77), but this time we want to make use of the proximity operators of the two functions. So, let X and U be two real Hilbert spaces. Let f ∈ Γ 0 (X ) and g ∈ Γ 0 (U ). Let L : X → U be a nonzero bounded linear operator. We want to The corresponding monotone inclusion is Again, to bypass the annoying operator L, we introduce an auxiliary variable u ∈ ∂g(Lx), which will be called the dual variable, so that the problem now consists in finding x ∈ X and u ∈ U such that u ∈ ∂g(Lx) Let us define the dual convex optimization problem associated to the primal problem (99): If a pair (x, u) ∈ X × U is a solution to (101), then x is a solution to (99) and u is a solution to (102).
To solve the primal and dual problems (99) and (102) jointly, Chambolle and Pock [28] proposed the following algorithms (without relaxation), see also Esser et al. [66] and Zhang et al. [137]: Chambolle-Pock iteration, form I, for (99) and (102): Chambolle-Pock iteration, form II, for (99) and (102): We can note that the Chambolle-Pock algorithm is self-dual: if we apply the Chambolle-Pock iteration form I to the problem (102), to minimizef +g •L withf = g * ,g = f * , L = −L * , and the roles of x and u are switched, as well as the roles of τ and σ, we obtain exactly the Chambolle-Pock iteration form II.
Chambolle and Pock proved the convergence in the finite-dimensional case, assuming that τ σ L 2 < 1 and ρ (i) = 1 [28]. The convergence was proved in a different way by He and Yuan, with a constant relaxation parameter ρ (i) = ρ ∈ (0, 2) and the same other hypotheses [82]; indeed, they observed that the Chambolle-Pock algorithm is a primal-dual proximal point algorithm to find a primal-dual pair z = (x, u) in Z = X × U , solution to the monotone inclusion The operator M : The Chambolle-Pock iteration form II satisfies the same primal-dual inclusion, but with −L replaced by L in P . Thus, the Chambolle-Pock iteration is a preconditioned primal-dual proximal point algorithm and, as a consequence of Theorem 3.3, we have [49, Theorem 3.2]: Theorem 5.1 (Chambolle-Pock algorithm (103) or (104)). Let x (0) ∈ X and u (0) ∈ U . Let τ > 0 and σ > 0 be such that στ L 2 < 1. Let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2 − ρ (i) ) = +∞. Then the sequences (x (i) ) i∈N and (u (i) ) i∈N defined either by the iteration (103) or by the iteration (104) converge weakly to a solution of (99) and to a solution of (102), respectively.
The difference between Theorem 5.1 and Theorem 5.2 is that στ L 2 = 1 is allowed in the latter. This is a significant improvement: in practice, one can set σ = 1/(τ L 2 ) in the algorithms and have only one parameter left, namely τ , to tune.
In Figure 2, we illustrate how overrelaxation can accelerate the Chambolle-Pock algorithm. We consider the problem of image inpainting, which consists in estimating missing pixel values of an image. We formulate inpainting as an optimization problem of the form (99); that is, we aim at minimizing Ψ = f + g • L, where f is the indicator function of the affine subspace of images having the prescribed pixel values, and g•L is the 'isotropic' total variation [30,51,115]. So, we use the Chambolle-Pock algorithm to solve the problem; a solution x ⋆ (which seems to be unique and is estimated by x (i) with i = 10 6 ) is shown in Figure 2 (c). We can see from the plots of x (i) − x ⋆ 2 and Ψ(x (i) ) − Ψ(x ⋆ ) = g(Lx (i) ) − g(Lx ⋆ ) (since f (x (i) ) = 0 for every i), that the best value of τ depends on the desired accuracy or number of iterations: τ = 0.01 is better after 10 3 iterations, whereas a smaller value is better after 10 4 iterations. In any case, ρ (i) = 1.9 significantly speeds up the convergence, in comparison with ρ (i) = 1.

The Proximal Method of Multipliers
Let X and U be two real Hilbert spaces. Let g ∈ Γ 0 (U ), c ∈ X and L : X → U be a bounded linear operator. We consider the problem and the dual problem minimize We can view the term x, c as a smooth function h with constant gradient ∇h = c and apply the Loris-Verhoeven algorithm. Thus, let τ > 0 and σ > 0, let x (0) ∈ X and u (0) ∈ U , and let (ρ (i) ) i∈N be a sequence of relaxation parameters. The iteration is (we change the name of the variable x into s): (107) and (108):

Loris-Verhoeven iteration for
Each iteration of the algorithm satisfies the inclusion Since the algorithm has the structure of a proximal point algorithm, we can apply Theorem 3.3 and obtain convergence under the condition στ L 2 < 1.
An alternative is to see the term x, c as a function f , with proximity operator prox τ f : x → x − τ c. So, we can solve (107) and (108) by applying the Chambolle-Pock algorithm (we change the name of its variable x (i+ 1 2 ) into s (i) ): Chambolle-Pock iteration, form I, for (107) and (108): We can discard the variable x and the iteration becomes: and we exactly recover the iteration (109). So, the Loris-Verhoeven algorithm and the Chambolle-Pock algorithm are the same for the specific problems (107) and (108); we call this algorithm the Proximal Method of Multipliers [113]. Furthermore, the form (111) provides us with a compact form of the algorithm: Proximal Method of Multipliers for (107) and (108): and since the two updates are independent, they can be performed in parallel. Theorem 5.2 yields: (113)). Let x (0) ∈ X and u (0) ∈ U . Let τ > 0 and σ > 0 be such that στ L 2 ≤ 1. Let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2 − ρ (i) ) = +∞. Then the sequences (x (i) ) i∈N and (u (i) ) i∈N defined by the iteration (113) converge weakly to a solution of (107) and to a solution of (108), respectively.

Theorem 5.3 (Proximal Method of Multipliers
We can apply the form II of the Chambolle-Pock algorithm as well, with same convergence properties as in Theorem 5.3: (107) and (108):
Note that the Douglas-Rachford algorithm is not a primal-dual proximal point algorithm, since the operator P in (106) is not strongly positive any more. So, its weak convergence, which is not strong in general [21], is difficult to establish; it was shown only in 2011 by Svaiter [127]. Let us have a quick look at the structure of the Douglas-Rachford algorithm. First, we define the auxiliary variable The Douglas-Rachford iteration only depends on this concatenated variable, not on the full pair (x (i) , u (i) ), and we can rewrite it as: Douglas-Rachford iteration for (115): or, equivalently, (115): for i = 0, 1, . . .
Note that it would be misleading to state Theorem 5.4 by invoking Theorem 5.2, since Theorem 5.2 is itself based on the convergence result of the Douglas-Rachford algorithm [7, Corollary 28.3]. In other words, weak convergence of the Douglas-Rachford algorithm is the fundamental result, from which we deduce weak convergence of the Chambolle-Pock algorithm, which generalizes it.
Let us show that, like the Chambolle-Pock algorithm, the Douglas-Rachford algorithm is self-dual [59, Lemma 3.6] [135]: it is equivalent to apply it to the primal problem (115) or to the dual problem: minimize Using the Moreau identity (25) and starting from the form (116), we can write the Douglas-Rachford iteration as: If we introduce the variables (i) = −1 τ s (i) = u (i) − 1 τ x (i) , we can rewrite the iteration as: which can be simplified as: (121):

Douglas-Rachford iteration for
Thus, we recognize the form (119) of the Douglas-Rachford algorithm applied to the dual problem (121), with parameter 1/τ . All in all, there are only two different instances of the Douglas-Rachford algorithm: the one given here, and the one obtained by switching the roles of f and g (which we could obtain as a particular case of the Chambolle-Pock algorithm, form II). Douglas-Rachford splitting has been shown to be the only way, in some sense, to minimize a sum of two functions, by calling their proximity operators [117].

A slightly modified version of the Douglas-Rachford algorithm
Let us consider a slightly modified version of the Douglas-Rachford algorithm, which we will use in the next section. In addition to f and g as before, let c ∈ X . We consider the problem: The dual problem is The linear term x, c does not add any difficulty: one can apply the Douglas-Rachford algorithm to minimizef + g, wheref = f + ·, c , using the fact that prox τf (x) = prox τ f (x − τ c). So, the modified Douglas-Rachford algorithm is: Douglas-Rachford iteration for (125): for i = 0, 1, . . .
Using the Moreau identity (25) and starting from the form (116), we can write the algorithm as Douglas-Rachford iteration for (125): (note that the variable v here is not the same as the one in (124) or (129); it is its opposite).
Thus, in the Douglas-Rachford algorithm, or equivalently the ADMM, there are two primal variables x and w and two dual variables u and v, as can be seen in (129). But depending on how the algorithm is written, only some of these variables appear. Also, two monotone inclusions are satisfied in parallel at every iteration: Note that in the ADMM form (130), the final extrapolation step from v (i+ 1 2 ) to v (i+1) accounts for the absence of relaxation on the variable w, in such a way that w (i+1) − τ v (i+1) takes the appropriate value. That is, , as if the two variables v and w had been relaxed as usually.
We can remark that τ should not be chosen as large as possible in the Douglas-Rachford algorithm or the ADMM: we see in the primal-dual inclusions (131) and (132) that the antagonistic values τ and 1/τ control the primal and dual updates. So there is a tradeoff to achieve in the choice of τ . We insist on the fact that the Douglas-Rachford algorithm and the ADMM are primal-dual algorithms: it is equivalent to apply them on the primal or on the dual problem. There are only two different algorithms to minimize f + g: the one given here, with prox f or prox f * applied first, and the one obtained by exchanging f and g.
Finally, one can find the following relaxed version of the ADMM in the literature: ADMM iteration for (115) and (121): which is proved to converge for ρ ∈ 0, ( √ 5 + 1)/2 [72]. This is different from the relaxation considered here. Both types of relaxation have been unified in an extended setting in [81].

Infimal postcompositions
Suppose that we want to solve an optimization problem, which involves a term g(Lx). Throughout the paper, we are studying splitting algorithms, which call the proximity operator of g, L, and L * separately. Here we look at an alternative, which consists in a change of variable: we introduce the variable r = Lx and we express the optimization problem with respect to r instead of x. Then an algorithm to solve the reformulated problem will construct a sequence (r (i) ) i∈N converging to r ⋆ = Lx ⋆ , where x ⋆ is a solution to the initial problem, that will be recovered as a byproduct. For instance, suppose that we want to minimize f (x) + g(Lx) over x ∈ X , where f , g, L : X → U are as before. This is the same as Thus, we can rewrite the problem with respect to r only, by 'hiding' x as a sub-variable, as follows: wheref : r ∈ U → inf x∈X f (x) + ı {r} (Lx) .f ∈ Γ 0 (U ) under some mild mathematical safeguards; for instance, we will require that the infimum is attained in the definition off . Clearly, minimizingf (r) + g(r) will force r to be in the range of L, elsef (r) = +∞. That is, a solution will be r ⋆ = Lx ⋆ for some x ⋆ ∈ arg min {x∈X :Lx=r ⋆ } f (x), which is the actual element we are interested in. To apply a proximal splitting algorithm to this problem, we do not need to evaluatef , but we must be able to apply its proximity operator That is, proxf = L(∂f + L * L) −1 L * , which is indeed single-valued: even if x ′ in its definition is not unique, Lx ′ is unique. We havef * : u ∈ U → sup x∈X u, Lx − f (x) = f * (L * u) [7,Proposition 13.24]. Thus, the dual problem to min x f (x) + g(Lx) = min r f (r) + g(r) is min u f * (−u) + g * (u) = min u f * (−L * u) + g * (u) , which is consistent with our formulations of dual problems so far. The functionf is called the infimal postcomposition of f by L, which is thereafter denoted by L ⊲ f [7,Definition 12.34]. The property to remember on it is: Thus, as we have seen, behind this technical notion, there is simply a change of variables, which is helpful if prox L⊲f , or by the Moreau identity prox f * •L * , is easy to compute.
Let us apply this notion of infimal postcomposition. Let X , Y, U be real Hilbert spaces. Let f ∈ Γ 0 (X ), g ∈ Γ 0 (Y), c ∈ U . Let K : Y → U and L : X → U be nonzero bounded linear operators. We consider the problem: The dual problem is minimize The corresponding primal-dual monotone inclusion is: find (x, y, u) ∈ X × Y × U , such that For instance, the minimization of f (x) + g(Lx) is a particular case of this problem, with c = 0, U = Y, K = −Id. We want to solve (139) by applying two infimal postcompositions, so that the problem amounts to find an element r ⋆ = Lx ⋆ = −Ky ⋆ + c, where (x ⋆ , y ⋆ ) is a solution to (139). That is, we rewrite (139) as where For the problem to be well posed, we suppose that the following assumptions hold: (i) The solution set of (141) is nonempty.
(ii) The infimal postcompositions are exact [7,Definition 12.34]; that is, for every r ∈ U , the infimum in (143) and (144) is attained, so it is a minimum (its value can be +∞). In other words, for every r ∈ ran L, the set of minimizers of f over the set L −1 {r} is nonempty. Likewise, for every r ∈ ran K, the set of minimizers of g over K −1 {−r} is nonempty.  (130) become Lx (i+ 1 2 ) and −Ky (i) + c in (142), respectively). Since the ADMM and the Douglas-Rachford algorithm are equivalent, we can write the iteration like in (119) instead: Douglas-Rachford iteration for (139): for i = 0, 1, . . .
In case ρ (i) > 1, the Douglas-Rachford form (146) is to be preferred to the ADMM form (145), since it involves less operations.
Here we used infimal postcompositions in the Douglas-Rachford algorithm, but they could be used in any other splitting algorithm. One should study, on a case-by-case basis, whether the subproblems corresponding to the proximity operators of the infimal postcompositions are easy to solve, and if this is the case, whether it is better to do so or to split the problem with separate calls to the linear operators.

The Generalized Chambolle-Pock algorithm
Let X , U , V be real Hil-bert spaces. Let f ∈ Γ 0 (V), g ∈ Γ 0 (U ), c ∈ X . Let K : X → V and L : X → U be nonzero bounded linear operators. We consider the following problem, generalizing (99): minimize The corresponding monotone inclusion, which we will actually solve, is We introduce two dual variables u ∈ U and v ∈ V, so that the problem is to find Accordingly, the dual problem is Clearly, if K = Id and L = Id, (147) and (150) revert to (125) and (126), and we can use the Douglas-Rachford algorithm (128) to solve these problems. On the other hand, if K = Id and c = 0, (147) and (150) revert to (99) and (102), and we can use the Chambolle-Pock algorithm (103) or (104).
The Generalized Chambolle-Pock algorithm (151) has appeared in the literature under different names, like Alternating Proximal Gradient Method [91], Generalized Alternating Direction Method of Multipliers [56], or Preconditioned ADMM [14,126]. The convergence results derived in this section generalize previously known results.
If K = Id and we set η = 1/τ , we recover the Chambolle-Pock algorithm form I, to minimizef + g • L, withf = f + ·, c . Indeed, in that case, the first two updates become: Therefore, the Generalized Chambolle-Pock algorithm indeed generalizes the Chambolle-Pock algorithm to any linear operator K. But we can remark that setting η = 1/τ is not allowed in Theorem 6.1. So, to extend the range of parameters to τ σ L 2 ≤ 1 and τ η K 2 ≤ 1, we now analyze the algorithm from another point of view: we show that it is a Douglas-Rachford algorithm in an augmented space. This analysis is inspired and generalizes the analysis of the Chambolle-Pock algorithm by O'Connor and Vandenberghe [98].
We can write the Generalized Chambolle-Pock algorithm in a different form, with only one call to K, K * , L, L * per iteration. For this, we introduce the scaled variablex = x/τ and auxiliary variables b = K * v and r = L * u+K * v+c. Set b (0) = K * v (0) and l (0) = L * u (0) +b (0) +c. This yields the iteration: Furthermore, let us show that we recover as a particular case the Loris-Verhoeven algorithm applied to the minimization of f • K + h, when h : x → 1 2 x, Qx + x, c is quadratic, with Q = L * L (again, given Q, such an operator L exists). Indeed, minimizing f • K + h is equivalent to minimizing f •K +g•L+ ·, c , with g = g * = 1 2 · 2 . Let us apply the generalized Chambolle-Pock algorithm (151) on this latter problem, with σ = 1. Since prox σg * = 1 2 Id, we can write the update of u as: so that, if u (0) = Lx (0) , we have u (i) = Lx (i) , for every i ∈ N. Hence, we can remove the variable u and rewrite the iteration as: which is exactly the Loris-Verhoeven iteration (90). Thus, the Loris-Verhoeven algorithm can be viewed as a primal-dual forward-backward algorithm, but also as a primal-dual Douglas-Rachford algorithm, when h is quadratic.

The Condat-Vũ algorithm
Let us consider the primal optimization problem: where f ∈ Γ 0 (X ), g ∈ Γ 0 (U ), h : X → R is a convex and differentiable function with β-Lipschitz continuous gradient ∇h, for some real β > 0, and L : X → U is a bounded linear operator. The corresponding monotone inclusion is Again, we introduce a dual variable u, so that we can rewrite the problem (164) with respect to a pair of objects z = (x, u) in Z = X × U : If (x, u) ∈ X × U is a solution to (165), then x is a solution to (163) and u ∈ ∂g(Lx) is a solution to the dual problem The operator M : [7,Proposition 26.32 (iii)] and C : Z → Z, (x, u) → (∇h(x), 0) is ξ-cocoercive, with ξ = 1/β. Thus, it is again natural to think of the forward-backward iteration, with preconditioning. The difference with the construction in section 4 is the presence of the nonlinear operator ∂f , which prevents us to express x (i+ 1 2 ) in terms of x (i) and u (i+ 1 2 ) . Instead, the iteration is made explicit by canceling the dependence of x (i+ 1 2 ) on u (i+ 1 2 ) in P . That is, the iteration, written in implicit form, is where τ > 0 and σ > 0 are two real parameters, z (i) = (x (i) , u (i) ) and z (i+ 1 2 ) = (x (i+ 1 2 ) , u (i+ 1 2 ) ). Thus, the primal-dual forward-backward iteration is: Condat-Vũ iteration, form I, for (163) and (166): for i = 0, 1, . . .
This algorithm was proposed independently by the first author [49] and by B. C. Vũ [132]. An alternative is to update u before x, instead of the converse. This yields a different algorithm, characterized by the primal-dual inclusion .
We can observe that if h = 0, the Condat-Vũ iteration reverts to the Chambolle-Pock iteration, so the former can be viewed as a generalization of the latter. Accordingly, if we set β = 0 in Theorem 7.1, we recover Theorem 5.1.
The Condat-Vũ algorithm and the Loris-Verhoeven algorithms are both primal-dual forward-backward algorithms, but they are different. When f = 0, larger values of τ and σ are allowed in the latter than in the former; this may be beneficial to the convergence speed in practice.
For the Condat-Vũ algorithm too, let us focus on the case where h is quadratic; that is, for some self-adjoint, positive, nonzero, bounded linear operator Q on X and some c ∈ X . We have β = Q . We can rewrite the primal-dual inclusion (167), which characterizes the Condat-Vũ iteration (168), as Similarly, we can rewrite the primal-dual inclusion (169), which characterizes the second form of the Condat-Vũ iteration (170), as (178), with −L replaced by L in P . In both cases, using the properties of the Schur complement, P is strongly positive if and only if τ Q + σL * L < 1 (which implies that τ < 1/β). A sufficient condition for this inequality to hold is τ (σ L 2 + β) < 1. However, in some applications, Q + σL * L may be smaller than σ L 2 + β, so that larger stepsizes τ and σ may be used when h is quadratic, for the benefit of convergence speed. Thus, when h is quadratic, the Condat-Vũ iteration can be viewed as a preconditioned primal-dual proximal point algorithm, just like the Chambolle-Pock iteration. Accordingly, we can apply Theorem 3.3 and obtain convergence, under the condition τ Q + σL * L < 1. Instead, let us apply the stronger Theorem 3.6 to allow τ Q + σL * L = 1. For this, let us go back to the forward-backward analysis in (167) and (169). We suppose that τ σ L 2 < 1. Then we can strengthen the analysis in (171)-(176): Hence, we have: Theorem 7.2 (Condat-Vũ algorithm (168) or (170), quadratic case). Suppose that h is quadratic. Let x (0) ∈ X and u (0) ∈ U . Let τ > 0 and σ > 0 be such that τ σ L 2 < 1 and τ Q + σL * L ≤ 1. Let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2 − ρ (i) ) = +∞. Then the sequences (x (i) ) i∈N and (u (i) ) i∈N defined either by the iteration (168) or by the iteration (170) converge weakly to a solution of (163) and to a solution of (166), respectively.

An Extended Generalized Chambolle-Pock algorithm
Now that we have seen the different ways to construct primal-dual algorithms based on the forward-backward or proximal point algorithms, we can imagine different variations. For instance, we can extend the Generalized Chambolle-Pock algorithm in section 6 to deal with additional smooth terms on u and v in the dual problem. Let us look at one particular case of this extension: we add a smooth quadratic term on v. The problem is: for some self-adjoint, positive, nonzero, bounded linear operator Q on V and some element t ∈ V. The problem (180) has applications, for instance, to solve inverse problems in imaging regularized with the total variation, in the way it is defined by the first author [51]. Again, we design a primal-dual forward-backward algorithm, whose iteration satisfies: Accordingly, the algorithm is: Primal-dual iteration for (180): for i = 0, 1, . . .

The Primal-Dual Three-Operator Splitting (PD3O) and the Davis-Yin algorithms
We consider the same problem as in the previous section, which involves three functions and a linear operator: minimize The dual problem is, again: Let τ > 0 and σ > 0, let s (0) ∈ X and u (0) ∈ U , and let (ρ (i) ) i∈N be a sequence of relaxation parameters. The Primal-Dual Three-Operator Splitting (PD3O) algorithm proposed by Yan [134] is: (see (216) for a more compact form with one call to ∇h and L * per iteration).
If h is quadratic (∇h is affine) or if ρ (i) = 1, we have s (i) = x (i) − τ ∇h(x (i) ) − τ L * u (i) and the PD3O iteration can be rewritten as Thus, we recover the Chambolle-Pock algorithm when h = 0. When f = 0, the PD3O algorithm reverts to the Loris-Verhoeven algorithm (with s (i) = x (i+ 1 2 ) in the former playing the role of x (i) in the latter).
O'Connor and Vandenberghe showed that the PD3O algorithm can be obtained from the Davis-Yin algorithm, shown below, by a reformulation in an augmented space, of the same type as the one showed for the Generalized Chambolle-Pock algorithm in section 6. Their analysis makes it possible to have στ L 2 = 1 [98]. So, from the results on the Davis-Yin algorithm stated in Theorem 8.3, we can extend Theorem 8.1 to στ L 2 ≤ 1 as follows: PD3O algorithm (185)). Suppose that X and U are of finite dimension. Let s (0) ∈ X and u (0) ∈ U . Let τ ∈ (0, 2/β) and σ > 0 be such that στ L 2 ≤ 1. Set δ = 2 − τ β/2. Let (ρ (i) ) i∈N be a sequence in [0, δ] such that i∈N ρ (i) (δ − ρ (i) ) = +∞. Then the sequences (x (i+ 1 2 ) ) i∈N and (u (i) ) i∈N defined by the iteration (185) converge to a solution of (183) and to a solution of (184), respectively.
Note that the finite dimension assumption is not necessary in the proof of this result [98], so it could be removed.
Thus, the conditions for convergence of the PD3O algorithm are the same as for the Loris-Verhoeven algorithm, so that we get the additional function f for free.
If f = 0, the Davis-Yin algorithm reverts to the forward-backward algorithm. If h = 0, it reverts to the Douglas-Rachford algorithm. The Davis-Yin algorithm generalizes the forward-Douglas-Rachford algorithm and the Generalized forward-backward algorithm [16,110,111].

A new algorithm for the quadratic case
Let us continue with the minimization of f + g • L + h, in the quadratic case where h(x) = 1 2 x, Qx + x, c , for some self-adjoint, positive, nonzero, bounded linear operator Q on X and some c ∈ X . Set β = Q . There seems to be no way to interpret the PD3O algorithm in a different way, to enlarge the relaxation range to (0, 2). However, we can propose a new algorithm, which has the widest range of parameters seen so far. Let τ ∈ (0, 2/β) and σ > 0 be such that στ L 2 < 1. We introduce the self-adjoint, strongly positive, bounded linear operator P in Z: Note that contrary to what we have done so far, only one half of Q appears in P . We define the real Hilbert space Z P as Z endowed with the inner product ·, · P : (z, z ′ ) → z, P z ′ . We will solve the following monotone inclusion in Z P : Let s (0) ∈ Z. The operators M and N are maximally monotone in Z P , so we will use the Douglas-Rachford algorithm to solve (191) (with stepsize γ = 1); it can be written as We have J M : (x, u) → prox τ f (x − τ 2 Qx − τ c), u , and the resolvent J N amounts to one iteration of the Loris-Verhoeven algorithm (92), so it maps (x, u) to (x ′ , u ′ ), with: Substituting these two expressions of the resolvents in (192), we obtain the new primal-dual Douglas-Rachford iteration: Primal-dual Douglas-Rachford iteration, form I, for (183): If we switch the roles of M and N in (192), we obtain the iteration: Primal-dual Douglas-Rachford iteration, form II, for (183): Let s (0) ∈ X and u (0) ∈ U . Let τ ∈ (0, 2/β) and σ > 0 be such that στ L 2 < 1. Let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2 − ρ (i) ) = +∞. Then the sequences (x (i+ 1 2 ) ) i∈N and (u (i) ) i∈N defined either by the iteration (194) or by the iteration (195) converge weakly to a solution of (183) and to a solution of (184), respectively.
In comparison with the Condat-Vũ and the PD3O algorithms, with the new algorithm, τ is allowed to be in (0, 2/β), the condition on τ and σ does not depend on Q, and the relaxation range is (0, 2). Like them, the new algorithm reverts to the Chambolle-Pock algorithm when Q = 0. The price to pay is two evaluations of Q per iteration.

Parallel versions of the algorithms
In this section, we generalize some of the algorithms seen so far, to make it possible to solve optimization problems involving an arbitrary number of functions g m , m = 1, . . . , M , instead of only one function g. The derived algorithms are parallel, in the sense that the proximity operators of the functions g m or g * m are independent of each other and can be applied in parallel. We present two techniques for this generalization, which consist in reformulating the problems in product spaces, an idea which has been known for a long time [105]. When M = 1, we recover the initial algorithms. But we should keep in mind that the variables, which are updated in parallel with respect to antagonistic functions, are then essentially averaged to form the update of the variable x of interest. This means that the larger M , the slower the convergence to a solution. For instance, it should be faster to use the Douglas-Rachford algorithm to minimize f +g than to apply any parallel algorithm to minimize g 1 +g 2 . Similarly, the Generalized Chambolle-Pock algorithm in section 6 to minimize f • K + g • L should be faster than the parallel version of the Proximal Method of Multipliers in subsection 5.1 to minimize g 1 • L 1 + g 2 • L 2 .
Technique 1: Let M ≥ 1 be an integer. Let X be a real Hilbert space. We consider the minimization of f (x) + M m=1 g m (x) + h(x) over x ∈ X , for some f ∈ Γ 0 (X ), g m ∈ Γ 0 (X ), differentiable and convex function h with β-Lipschitz continuous gradient, for some β ≥ 0. Let (ω m ) M m=1 be a sequence of positive weights, whose sum is 1. We introduce the real Hilbert space X = X × · · · × X (M times), endowed with the inner product · , · X : (x, x ′ ) → M m=1 ω m x m , x ′ m . We introduce the function ı = : If we switch the roles of f and g, we obtain the iteration: Douglas-Rachford iteration, form II, for (196): for i = 0, 1, . . .
As an application of Theorem 5.2, we have:
When h is quadratic, as an application of Theorem 7.2, we have: Let (ρ (i) ) i∈N be a sequence in [0, 2] such that i∈N ρ (i) (2 − ρ (i) ) = +∞. Then the sequences (x (i) ) i∈N and (u (i) 1 , . . . , u (i) M ) i∈N defined either by the iteration (209) or by the iteration (210) converge weakly to a solution of (207) and to a solution of (208), respectively.

Conclusion
We have made a small tour of proximal splitting algorithms and we have shown how a principled analysis in terms of primal-dual monotone inclusions makes it possible to derive existing and new algorithms, with general convergence guarantees. Several connections between apparently distinct algorithms have been established, and we have unleashed the relaxation potential of a large class of proximal splitting algorithms 2 . Finally, we have derived parallel versions of the algorithms to minimize the sum of an arbitrary number of functions. This selected overview of proximal splitting methods is by no way exhaustive; for instance, we did not cover the forward-backward-forward splitting method and its applications [18,43,131] [7, section 26.6], or Dykstra-like algorithms [5,29] [7, section 30.2]. Moreover, research around proximal splitting techniques is still an active topic, with recent contributions like the asymmetric forward-backward-adjoint algorithm [87], the forward-backward-half forward algorithm [17], the forward-reflected-backward algorithm [94], or the Forward-Douglas-Rachford-Forward algorithm [119]. Moreover, there are many avenues for extensions, including the use of variable stepsizes or metrics [41,46,93,104], the study of convergence rates and development of accelerated variants [8,9,26,31,52,54,58,84,89,116,128], or the design of randomized versions [25,27,44,80,103,112,120]. It will be very interesting to see proximal splitting methods being enriched by these modern notions.