Our contributions are the following:

(1) New algorithms: We propose the first distributed algorithms to solve (Eq. 1) in whole generality, with proved convergence to an exact solution, and having the full splitting, or decoupling, property: ∇F _m , p r o x H m , K _m and K m * are applied at the m-th node, and the proximity operator of R is applied at the master node connected to all others. No other more complicated operation, like an inner loop or a linear system to solve, is involved.

Many estimation problems in a wide range of scientific fields can be formulated as large-scale convex optimization problems (Palomar and Eldar, 2009; Sra et al., 2011; Bach et al., 2012; Bubeck, 2015; Polson et al., 2015; Chambolle and Pock, 2016a; Glowinski et al., 2016; Stathopoulos et al., 2016; Condat, 2017a; Condat et al., 2019b). Proximal splitting algorithms (Combettes and Pesquet, 2010; Boţ et al., 2014; Parikh and Boyd, 2014; Komodakis and Pesquet, 2015; Beck, 2017; Condat et al., 2019a) are particularly well suited to solve them; they consist of simple, easy to compute, steps that can deal with the terms in the objective function separately.

These algorithms are generally designed as sequential ones, for M = 1, and then they can be extended by lifting in product space to parallel versions, well suited to minimize F + R + ∑ _m H _m ○ K _m , see for instance Condat et al., 2019a, Section 8. However, it is not straightforward to adapt lifting to the case of a finite-sum F = 1 M ∑ m F m , with each function F _m handled by a different node, which is of primary importance in machine learning. This generalization is one of our contributions.

There is a vast literature on distributed optimization to minimize 1 M ∑ m F m + R , with a focus on strategies based on (block-)coordinate or randomized activation, as well as replacing the gradients by cheaper stochastic estimates (Cevher et al., 2014; Richtárik and Takáč, 2014; Gorbunov et al., 2020; Salim et al., 2020). Replacing the full gradient by a stochastic oracle in the accelerated algorithms with varying stepsizes we propose is not straightforward; we leave this direction for future research. In any case, the generalized setting, with the smooth functions F _m at the nodes supplemented or replaced by nonsmooth functions H _m , possibly composed with linear operators, seems to have received little attention. We want to make up for that. Decentralized optimization over networks is an active research topic (Latafat et al., 2019; Alghunaim et al., 2021). In this paper, we focus on the centralized client–server model, with one master node connected to several client nodes, working in parallel. We leave the study of decentralized algorithms for future work.

When M = 1 and K = I, where I denotes the identity, Davis and Yin (2017) proposed an efficient algorithm, along with an extensive study of its convergence rates and possible accelerations. But the ability to handle a nontrivial K is behind the success of the Chambolle and Pock (2011) or \Condat (2013), Vũ (2013): they are well suited for regularized inverse problems in imaging (Chambolle and Pock, 2016a), for instance with the total variation and its variants (Bredies et al., 2010; Condat, 2014, 2017b; Duran et al., 2016); other examples are computer vision problems (Cremers et al., 2011), overlapping group norms for sparse estimation in data science (Bach et al., 2012), and trend filtering on graphs (Wang et al., 2016). Another prominent case is when H is an indicator function, so that the problem becomes: minimize F(x) + R(x) subject to Kx = b. If K is a gossip matrix like the minus graph Laplacian, decentralized optimization over a network can be tackled (Shi et al., 2015; Scaman et al., 2017; Salim et al., 2021).

When M = 1 and K is arbitrary, there exist algorithms to solve (Eq. 1) in full generality, for example, the Combettes and Pesquet (2012), Condat, (2013), Vũ (2013), PD3O (Yan, 2018) and PDDY (Salim et al., 2020) algorithms. However, their convergence rates and possible accelerations are little understood. Our main contribution is to derive new convergence rates and accelerated versions of the PD3O and PDDY algorithms, and their particular cases, including Chambolle and Pock (2011) and Loris and Verhoeven (2011) algorithms. In order to do this, we show that these two algorithms can be viewed as instances of the Davis–Yin algorithm. This reformulation technique is inspired by the recent one of O’Connor and Vandenberghe (O’Connor and Vandenberghe, 2020); it makes it possible to split the composition H°K and to derive algorithms, which call the operators prox _H , K, K* separately. This technique is fundamentally different from the one in Salim et al. (2020), showing that the PD3O and PDDY algorithms are primal–dual instances of the operator version of Davis–Yin splitting to solve monotone inclusions. Notably, we can derive convergence rates with respect to the objective function and accelerations, which is not possible with the primal–dual reformulation of Salim et al. (2020). On the other hand, the latter encompasses the Condat–Vũ algorithm (Condat, 2013; Vũ, 2013), which is not the case of our approach. So, these are complementary interpretations.

In Section 2, we propose new nonstationary versions (i.e. with varying stepsizes) of several algorithms for optimization problems made of three terms, and we analyze their convergence rates. The derivation details are pushed to the end of the paper in Section 5 for ease of reading. In Section 3, we further propose distributed algorithms, which can minimize the sum of an arbitrary number of terms. Again, the derivation details are deferred to Section 6. Numerical experiments illustrating the good match between our theoretical results and practical performance are shown in Section 4.

The main difficulty in (Eq. 2) is the presence of the linear operator K. Indeed, if K = I, the Davis–Yin algorithm (Davis and Yin, 2017) is well suited to minimize F + R + H. Note that there is a minor mistake in the way Algorithm 3 in Davis and Yin (2017) is initialized. This is corrected here. Thus, the Davis–Yin algorithm is as follows:

Let ( γ k ) k ∈ N be a sequence of stepsizes. Let x H 0 ∈ X and u 0 ∈ X . For k = 0, 1, … iterate ⌊ x k + 1 = p r o x γ k R ( x H k + γ k u k ) u k + 1 = u k + 1 γ k ( x H k − x k + 1 ) x H k + 1 = p r o x γ k + 1 H x k + 1 − γ k + 1 u k + 1 − γ k + 1 ∇ F ( x k + 1 ) . (4)

To make this algorithm applicable to K ≠ I, we reformulate the problem (Eq. 2) as follows:

(1) We choose a value η ≥‖K‖²; we recommend to set η = ‖K‖² in practice. Then there exists a real Hilbert space W and a linear operator C : W → U such that KK* + CC* = ηI. C is not unique, for instance, we can set C = ( η I − K K * ) 1 / 2 . We actually don’t need to exhibit C, its existence is sufficient here and there will be no call to C in the algorithms.

We have ∇ F ̃ ( x , w ) = ( ∇ F ( x ) , μ F w ) , p r o x R ̃ ( x , w ) = ( p r o x R ( x ) , 0 ) . Most importantly, for every γ > 0, we have (O’Connor and Vandenberghe, 2020): p r o x H ̃ * / γ ( x , w ) = ( K * u , C * u ) , where u = p r o x H * / ( γ η ) ( K x + C w ) / η . (6)

Note that in O’Connor and Vandenberghe (2020), the authors use F ̃ ( x , w ) = F ( x ) , whereas we add μ F 2 ‖ w ‖ 2 . This difference is essential, so that F ̃ is L _F -smooth and μ _F -strongly convex. Also, R ̃ is μ _R -strongly convex.

Then, we can apply the Davis–Yin algorithm (4) to solve the problem (Eq. 5). We set F, R, H in (Eq. 4) as F ̃ , R ̃ , H ̃ , respectively. The details of the substitutions yielding the algorithms are deferred to Section 5 for the convenience of reading; most notably, whenever CC* appears, it is replaced by ηI − KK*. The obtained algorithms turns out to be a nonstationary version of the PD3O algorithm (Yan, 2018), shown above. On the other hand, if we exchange the two functions and set F, R, H in (Eq. 4) as F ̃ , H ̃ , R ̃ , we obtain a different algorithm. It turns out to be a nonstationary version of the PDDY algorithm proposed recently (Salim et al., 2020), shown above too. With constant stepsizes γ _k ≡ γ ∈ (0, 2/L _F ), for both the PD3O and PDDY algorithms, x ^k and u ^k converge to some solutions x ^* and u ^* of (Eq. 2) and (Eq. 3), respectively; this result was known for η > ‖K‖² (Yan, 2018; Salim et al., 2020) and shown for η = ‖K‖² for the PD3O algorithm in O’Connor and Vandenberghe (2020), but convergence with η = ‖K‖² for the PDDY algorithm, as stated in Theorem 2, is new.

Particular cases of the PD3O and PDDY algorithms, which are shown above, are the following:

(1) If K = I and η = 1, the PD3O algorithm reverts to the Davis–Yin algorithm (Eq. 4); the PDDY algorithm too, but with H and R exchanged in (Eq. 4).

We first give convergence rates for the PD3O algorithm with constant stepsizes.

Theorem 1

(convergence rate of the PD3O algorithm). In the PD3O algorithm, suppose that γ _k ≡ γ ∈ (0, 2/L _F ) and η ≥‖K‖² . Then x ^k and u ^k converge to some solutions x ^* and u ^* of (2) and (3), respectively. In addition, suppose that H is continuous on an open ball centered at Kx ^* . Then the following hold: ( i ) Ψ ( x k ) − Ψ ( x * ) = o ( 1 / k ) .

Define the weighted ergodic iterate x ̄ k = 2 k ( k + 1 ) ∑ i = 1 k i x i , for every k ≥ 1. Then ( i i ) Ψ ( x ̄ k ) − Ψ ( x * ) = O ( 1 / k ) .

Furthermore, if H is L-smooth for some L > 0, we have a faster decay for the best iterate so far: ( i i i ) min i = 1 , … , k Ψ ( x i ) − Ψ ( x * ) = o ( 1 / k ) .

Proof

The convergence of x ^k follows from Davis and Yin, 2017, Theorem 2.1 and the convergence of u ^k follows from the one of the variable u B k = ( z k − x A k ) / γ in the notations of Davis and Yin (2017). (i) follows from Davis and Yin, 2017, Theorem 3.1, using the following facts; first, in this theorem, the function corresponding to H ̃ is supposed to be Lipschitz-continuous on a certain ball, but since the rate is asymptotic and Kx ^k → Kx ^*, it is sufficient to consider the property around Kx ^*; second, it is well known that if a convex real-valued function is continuous on a convex open set, it is Lipschitz-continuous on every compact subset of this set (Unknown author, 1972); third, if H is continuous, H ̃ is continuous too. (ii) follows from Davis and Yin (2017), Theorem 3.2 and (iii) follows from Theorem D.5 in the preprint of Davis and Yin (2017).□

Theorem 1 applies to the particular cases of the PD3O algorithm, like the Loris–Verhoeven, Chambolle–Pock, Douglas–Rachford algorithms. Our results are new even for them.

Remark 1

We can note that the forward–backward algorithm x ^k+1 = prox _γR (x ^k − γ∇F(x ^k )), which is a particular case of the PD3O algorithm when H = 0, is monotonic. So, the best iterate so far is the last iterate. Hence, Theorem 1 (iii) yields Ψ(x ^k ) − Ψ(x ^*) = o(1/k) for the forward–backward algorithm.

For the PDDY algorithm, we cannot derive a similar theorem, since R ̃ is not continuous around (x ^*, 0). Still, we can establish convergence of the variables:

Theorem 2

(convergence of the PDDY algorithm). In the PDDY algorithm, suppose that γ _k ≡ γ ∈ (0, 2/L _F ) and η ≥‖K‖² . Then x ^k and x R k both converge to some solution x ^* of (Eq. 2), and u ^k converges to some solution u ^* of (Eq. 3).

Proof

The convergence of x ^k and x R k to the same solution x ^* of (Eq. 2) follows from Davis and Yin, 2017, Theorem 2.1. The convergence of the variable u B k = ( z k − x A k ) / γ , in the notations of Davis and Yin (2017), implies in our setting, according to (6), that K*u ^k and C*u ^k both converge to some elements. But since ηu ^k = KK*u ^k + CC*u ^k , u ^k converges to some element u * ∈ U . Finally, we have x ^* = prox _γR (x ^* − γ∇F(x ^*) − γK*u ^*), so that 0 ∈ ∂R(x ^*) + ∇F(x ^*) + K*u ^*, and u * = p r o x H * / ( γ η ) ( u * + 1 γ η K x * ) , so that Kx ^* ∈ (∂H)⁻¹(u ^*). Hence, u ^* is a solution to (Eq. 3).□

We now give accelerated convergence results using varying stepsizes, when F or R is strongly convex; that is, μ _F + μ _R > 0. In that case, we denote by x ^* the unique solution to (Eq. 2).

Theorem 3

(convergence rate of the accelerated PD3O algorithm). Suppose that μ _F + μ _R > 0. Let κ ∈ (0, 1) and γ ₀ ∈ (0, 2(1 − κ)/L _F ). Set γ ₁ = γ ₀ and γ k + 1 = − γ k 2 μ F κ + γ k ( γ k μ F κ ) 2 + 1 + 2 γ k μ R 1 + 2 γ k μ R , for every  k ≥ 1 . (7)

Suppose that η ≥‖K‖ ² . Then in the PD3O algorithm, there exists c ₀ > 0 (whose expression is given in Section 5 ) such that, for every k ≥ 1, ‖ x k + 1 − x * ‖ 2 ≤ γ k + 1 2 1 − γ k + 1 μ F κ c 0 = O 1 / k 2 .

Proof

This result follows from Davis and Yin, 2017, Theorem 3.3, stated for convenience as Lemma 1 in Section 5.□

Note that with the stepsize rule in (Eq. 7), we have k γ _k → 1/(μ _F κ + μ _R ) as k → + ∞, so that γ _k = O(1/k) and γ _k+1/γ _k → 1. Also, when F = 0, L _F can be taken arbitrarily small, so that we can choose any γ ₀ > 0.

Theorem 3 is new for the PD3O and Loris–Verhoeven algorithms, but has been derived in O’Connor and Vandenberghe (2020) for the Chambolle–Pock algorithm. For the forward–backward algorithm, strong convexity yields linear convergence with constant stepsizes, so this nonstationary version does not seem interesting.

Concerning the PDDY algorithm, H ̃ is not necessarily strongly convex, even if H is. So, we only consider the case where F is strongly convex. As a consequence of Lemma 1, we get:

Theorem 4

(convergence rate of the accelerated PDDY algorithm). Suppose that μ _F > 0. Let κ ∈ (0, 1) and γ ₀ ∈ (0, 2(1 − κ)/L _F ). Set γ ₁ = γ ₀ and γ k + 1 = − γ k 2 μ F κ + γ k ( γ k μ F κ ) 2 + 1 , for every  k ≥ 1 . (8)

Suppose that η ≥‖K‖² . Then in the PDDY algorithm, there exists c ₀ > 0 (whose expression is given in Section 5) such that, for every k ≥ 1, ‖ x k + 1 − x * ‖ 2 ≤ γ k + 1 2 1 − γ k + 1 μ F κ c 0 = O 1 / k 2 .

Moreover, if η > ‖K‖², ‖ x R k − x * ‖ 2 = O ( 1 / k 2 ) as well.

Finally, we consider the case where, in addition to strong convexity of F or R, H is smooth; in that case, the algorithms with constant stepsizes converge linearly; that is, as a consequence of Lemma 2, we have:

Theorem 5

(linear convergence of the PD3O and PDDY algorithms). Suppose that μ _F + μ _R > 0 and that H is L _H -smooth, for some L _H > 0. Let x ^* and u ^* be the unique solutions to (2) and (3), respectively. Suppose that γ _k ≡ γ ∈ (0, 2/L _F ) and η ≥‖K‖² . Then the PD3O algorithm converges linearly: there exists ρ ∈ (0, 1] such that, for every k ∈ N , ‖ x k + 1 − x * ‖ 2 ≤ ( 1 − ρ ) k ‖ γ q 0 − x * + γ ∇ F ( x * ) − γ K * ( u 0 − u * ) ‖ 2 + γ 2 η ‖ u 0 − u * ‖ 2 − γ 2 ‖ K * ( u 0 − u * ) ‖ 2 .

The PDDY algorithm converges linearly too: there exists ρ ∈ (0, 1] such that, for every k ∈ N , ‖ x R k + 1 − x * ‖ 2 ≤ 4 ( 1 − ρ ) k ‖ x R 0 − x * + γ K * ( u 0 − u * ) ‖ 2 + γ 2 η ‖ u 0 − u * ‖ 2 − γ 2 ‖ K * ( u 0 − u * ) ‖ 2 .

Linear convergence of the other variables in the algorithms can be derived as well, see Proposition 1. Lower bounds for ρ can be derived from Theorem D.6 in the preprint version of Davis and Yin (2017). We don’t provide them, since they are not tight, as noticed in Remark D.2 of the same preprint. For instance, for the PDDY or Loris–Verhoeven algorithms with μ _F > 0, ρ = γ μ F ( 2 − γ L F ) ( 1 + γ η L H ) 2 .

If H = 0, by setting L _H = 0, we get ρ = γμ _F (2 − γL _F ). But then the PDDY algorithm reverts to the forward–backward algorithm, for which it is known that 1 − ρ = ( 1 − γ μ F ) 2 whenever γ ≤ 2/(L _F + μ _F ), which corresponds to the larger value ρ = γμ _F (2 − γμ _F ).

We emphasize that linear convergence comes for free with the algorithms, if the conditions are met, without any modification. That is, there is no need to know μ _F , μ _R , L _H , since the conditions on the two parameters γ and η do not depend on these values. For the particular case of the Chambolle–Pock algorithm, as pointed out in O’Connor and Vandenberghe (2020), this is in contrast to existing linear convergence results (Chambolle and Pock, 2016a), derived for a modified version of the algorithm, which depends on these values.

We first consider the non-distributed problem (Eq. 2), for the imaging inverse problem of deblurring, which consists in restoring an image y corrupted by blur and noise (Chambolle and Pock, 2016a). So, we set F : x ↦ 1 2 ‖ A x − y ‖ 2 , where the linear operator A corresponds to a 2-D convolution with a lowpass filter, with L _F = 1. The filter is approximately Gaussian and chosen so that F is μ _F -strongly convex with μ _F = 0.01. y is obtained by applying A to the classical 256 × 256 Shepp–Logan phantom image, with additive Gaussian noise. R = ı ₀ enforces nonnegativity of the pixel values. H°K corresponds to the classical ‘isotropic’ total variation (TV) (Chambolle and Pock, 2016a; Condat, 2017b), with H = 0.6 times the l _1,2 norm and K the concatenation of vertical and horizontal finite differences.

We compare the nonaccelerated, i.e. with constant γ _k , and accelerated versions, with decaying γ _k , of the PD3O, PDDY and Condat–Vũ algorithms. We initialize the dual variables at zero and the estimate of the solution as y. We set γ ₀ = 1.7, κ = 0.15, η = 8 ≥‖K‖² (except for the accelerated Condat–Vũ algorithm proposed in Chambolle and Pock (2016b), for which η = 16 and γ = 0.5).

The results are illustrated in Figure 1 (implementation in Matlab). We observe that the PD3O and PDDY algorithms have almost identical variables: the pink, red, blue curves are superimposed; we know that both algorithms are identical and revert to the Loris–Verhoeven algorithm when R = 0. Here R ≠ 0 but the nonnegativity constraint does not change the solution significantly, which explains the similarity of the two algorithms.

Note that x ^k in the PDDY algorithm is not feasible with respect to nonnegativity, and the red curve actually shows F(x ^k ) + H(Kx ^k ) − Ψ(x ^*). In the nonaccelerated case, Ψ(x ^k ) decays faster than O(1/k) but slower than O(1/k ²), which is coherent with Theorem 1. The same holds for ‖ x k − x * ‖ 2 ≤ 2 μ F ( Ψ ( x k ) − Ψ ( x * ) ) .

The accelerated versions improve the convergence speed significantly: Ψ(x ^k ) and ‖x ^k − x ^*‖² decay even faster than O(1/k ²), in line with Theorem 3 and Theorem 4. In all cases, the Condat–Vũ algorithm is outperformed. Also, there is no interest in considering the ergodic iterate instead of the last iterate, since the former converges at the same asymptotic rate as the latter, but slower.

We consider the same deblurring experiment as before, but we make H smooth by taking the Huber function instead of the l ₁ norm in the total variation; that is, λ|⋅| in the latter is replaced by h : t ∈ R ↦ λ 2 ν t 2 if | t | ≤ ν , λ | t | − ν 2 otherwise , for some ν > 0 and λ > 0 (set here as 0.1 and 0.6, respectively). We can also write h without branching as h ( t ) = λ 2 ν max ( ν − | t | , 0 ) 2 + λ ( | t | − ν 2 ) . It is known that h is L _h -smooth with L _h = λ/ν. For any γ > 0 and t ∈ R , we have p r o x h * / γ ( t ) = t / max ( | t | / λ , 1 + ν λ γ ) . Except for H, everything is unchanged.

The results are illustrated in Figure 2. Again, the PD3O and PDDY algorithms behave very similarly; they converge linearly, as proved in Theorem 5, and achieve machine precision in finite time. x ^k in the PDDY algorithm is not feasible and F(x ^k ) + H(Kx ^k ) − Ψ(x ^*) (red curve) takes negative values (not shown in log scale); so, x R k is the variable to study in this setting. We tested the ‘accelerated’ versions of the algorithms with decaying γ _k , but in this scenario, they are much slower and not suitable. Again, the Condat–Vũ algorithm is outperformed and the ergodic sequences converge much slower. Interestingly, the image x ^* is visually the same with TV and with Huber-TV.

Here we consider Problem (Eq. 1) in the special case with X = R d , for some d ≥ 1, F _m ≡ 0, and K _m ≡ I; that is, the problem of minimizing Ψ ( x ) = 1 M ∑ m = 1 M H m ( x ) + R ( x ) . (9)

In particular, to train a binary classifier, we consider the classical SVM problem with hinge loss, which has the form (Eq. 9) with R ( x ) = α 2 ‖ x ‖ 2 , for some α > 0, and H m ( x ) = max ( 1 − b m a m T x , 0 ) , with data samples a m ∈ R d and b _m ∈ { − 1, 1}.

For any γ > 0 we have prox _γR (x) = x/(1 + γα). We could view the dot product x ↦ b m a m T x as a linear operator K _m , but it is more interesting to integrate it in the function H _m . Indeed, as is perhaps not well known, the proximity operator of H _m has a closed form: for any γ > 0, p r o x γ H m : x ∈ R d ↦ x − b m η m max min ( b m a m T x − 1,0 ) , − η m γ a m , where η m = a m T a m = ‖ a m ‖ 2 . Thus, we use the Distributed Douglas–Rachford algorithm, a particular case of the distributed PD3O and PDDY algorithms. Since R is α-strongly convex, we also use the accelerated version of the algorithm with varying stepsizes, like in Theorem 3. We can note that in the context of Federated learning (Konečný et al., 2016; Malinovsky et al., 2020), where each m corresponds to the smart phone or computer of a different user with its own data (a _m , b _m ) stored locally, the problem is solved in a collaborative way but with preserved privacy, without the users sharing their data.

The method was implemented in Python on a single machine and tested on the dataset ‘australian’ from the LibSVM base (Chang and Lin, 2011), with d = 15 and M = 680. We set ω _m ≡ 1/M, α = 0.1, γ ₀ = 0.1, and we used zero vectors for the initialization. The results are shown in Figure 3. Despite the oscillations, we observe that both the objective suboptimality and the squared distance to the solution converge sublinearly, with rates looking like o ( 1 / k ) and O(1/k ²) for the nonaccelerated and accelerated algorithms, respectively, as guaranteed by Theorem 1 and Theorem 3. The proposed accelerated version of the distributed Douglas–Rachford algorithm yields a significant speedup.

In this section, we state the results on the Davis–Yin algorithm, which we be needed to analyze the PD3O and PPDY algorithms.

The Davis–Yin algorithm to minimize the sum of 3 convex functions F ̃ + G + J over a real Hilbert space Z (assuming that there exists a solution z ^* such that 0 ∈ ∇ F ̃ ( z * ) + ∂ G ( z * ) + ∂ J ( z * ) ) is (Davis and Yin, 2017):

Let z J 0 ∈ Z , u G 0 ∈ Z . For k = 0, 1, … iterate: ⌊ z G k + 1 = p r o x γ k G ( z J k + γ k u G k ) u G k + 1 = u G k + 1 γ k ( z J k − z G k + 1 ) z J k + 1 = p r o x γ k + 1 J z G k + 1 − γ k + 1 u G k + 1 − γ k + 1 ∇ F ̃ ( z G k + 1 ) . (10)

Equivalently, introducing the variable r k : = z J k + γ k u G k : let r 0 ∈ Z . For k = 0, 1, … iterate: ⌊ z G k + 1 = p r o x γ k G ( r k ) z J k + 1 = p r o x γ k + 1 J 1 + γ k + 1 γ k z G k + 1 − γ k + 1 γ k r k − γ k + 1 ∇ F ̃ ( z G k + 1 ) r k + 1 = z J k + 1 + γ k + 1 γ k ( r k − z G k + 1 ) . (11)

Equivalently: let r 0 ∈ Z . For k = 0, 1, … iterate: ⌊ z G k + 1 = p r o x γ k G ( r k ) u J k + 1 = p r o x J * / γ k + 1 1 γ k + 1 + 1 γ k z G k + 1 − 1 γ k r k − ∇ F ̃ ( z G k + 1 ) r k + 1 = z G k + 1 − γ k + 1 ∇ F ̃ ( z G k + 1 ) − γ k + 1 u J k + 1 . (12)

In our notations, Theorem 3.3 of Davis and Yin (2017) translates into Lemma 1 as follows; we assume that F ̃ is L F ̃ -smooth and μ F ̃ -strongly convex and that G is μ _G -strongly convex, for some L F ̃ > 0 , μ F ̃ ≥ 0 , μ _G ≥ 0.

Lemma 1

(accelerated Davis–Yin algorithm). Suppose that μ F ̃ + μ G > 0 . Let z ^* be the unique minimizer of F ̃ + G + J ; that is, 0 ∈ ∇ F ̃ ( z * ) + ∂ G ( z * ) + ∂ J ( z * ) . Let u G * be such that u G * ∈ ∂ G ( z * ) and 0 ∈ ∇ F ̃ ( z * ) + ∂ J ( z * ) + u G * . Let κ ∈ (0, 1) and γ 0 ∈ ( 0,2 ( 1 − κ ) / L F ̃ ) . Set γ ₁ = γ ₀ and γ k + 1 = − γ k 2 μ F ̃ κ + γ k ( γ k μ F ̃ κ ) 2 + 1 + 2 γ k μ G 1 + 2 γ k μ G , for every  k ≥ 1 .

Then, for every k ≥ 1, ‖ z G k + 1 − z * ‖ 2 ≤ γ k + 1 2 1 − γ k + 1 μ F ̃ κ c 0 = O 1 / k 2 , where c 0 = 1 − γ 0 μ F ̃ κ γ 0 2 ‖ z G 1 − z * ‖ 2 + ‖ u G 1 − u G * ‖ 2 .

Note that u G 1 = ( r 0 − z G 1 ) / γ 0 .

Linear convergence occurs in the following conditions, according to Theorem D.6 in the preprint version of Davis and Yin (2017), which translates into Lemma 2 as follows. We assume that F ̃ is L F ̃ -smooth and μ F ̃ -strongly convex, G is μ _G -strongly convex, and J is μ _J -strongly convex, for some L F ̃ > 0 , μ F ̃ ≥ 0 , μ _G ≥ 0, μ _J ≥ 0. We consider constant stepsizes γ _k ≡ γ, for some γ ∈ ( 0,2 / L F ̃ ) .

Lemma 2

(linear convergence of the Davis–Yin algorithm). Suppose that μ F ̃ + μ G + μ J > 0 and that G is L _G -smooth, for some L _G > 0, or J is L _J -smooth, for some L _J > 0. Let z ^* be the unique minimizer of F ̃ + G + J ; that is, 0 ∈ ∇ F ̃ ( z * ) + ∂ G ( z * ) + ∂ J ( z * ) . The dual problem of minimizing ( F ̃ + J ) * ( − u ) + G * ( u ) over u ∈ Z is strongly convex too; let u G * be its unique solution. We have u G * ∈ ∂ G ( z * ) and 0 ∈ ∇ F ̃ ( z * ) + ∂ J ( z * ) + u G * . Set r * = z * + γ u G * . Then, the Davis–Yin algorithm (Eq. 11) converges linearly: there exists ρ ∈ (0, 1] such that, for every k ∈ N , ‖ r k − r * ‖ 2 ≤ ( 1 − ρ ) k ‖ r 0 − r * ‖ 2 . (13)

Loose lower bounds for ρ are given in Davis and Yin, 2017, Theorem D.6.

We have the following corollary of Lemma 2:

Proposition 1

(linear convergence of the other variables in the Davis–Yin algorithm). In the same conditions and notations as in Lemma 2 , we have, for every k ∈ N , ‖ z G k + 1 − z * ‖ 2 ≤ ( 1 − ρ ) k ‖ r 0 − r * ‖ 2 ‖ z J k + 1 − z * ‖ 2 ≤ 4 ( 1 − ρ ) k ‖ r 0 − r * ‖ 2 . (14)

Also, in the form (Eq. 12) of the algorithm, ‖ u J k + 1 + u G * + ∇ F ̃ ( z * ) ‖ 2 ≤ 4 γ 2 ( 1 − ρ ) k ‖ r 0 − r * ‖ 2 and, in the form (Eq. 10) of the algorithm, ‖ u G k + 1 − u G * ‖ 2 ≤ 1 γ 2 ( 1 − ρ ) k ‖ r 0 − r * ‖ 2 .

Proof

Let k ∈ N . By nonexpansiveness of the proximity operator, in view of the first line in (Eq. 11), we have ‖ z G k + 1 − z * ‖ ≤ ‖ r k − r * ‖ , so that (Eq. 14) follows from (Eq. 13). In addition, in view of the second line in (Eq. 11), we have ‖ z J k + 1 − z * ‖ 2 ≤ ‖ 2 ( z G k + 1 − z * ) − ( r k − r * ) − γ ( ∇ F ̃ ( z G k + 1 ) − ∇ F ̃ ( z * ) ) ‖ 2 = ‖ ( z G k + 1 − z * ) − ( r k − r * ) + ( I − γ ∇ F ̃ ) ( z G k + 1 ) − ( I − γ ∇ F ̃ ) ( z * ) ‖ 2 = ‖ ( I − p r o x γ G ) ( r k ) − ( I − p r o x γ G ) ( r * ) + ( I − γ ∇ F ̃ ) ( z G k + 1 ) − ( I − γ ∇ F ̃ ) ( z * ) ‖ 2 and, by nonexpansiveness of I − prox _γG and I − γ ∇ F ̃ , ‖ z J k + 1 − z * ‖ 2 ≤ ‖ r k − r * ‖ + ‖ z G k + 1 − z * ‖ 2 ≤ 4 ‖ r k − r * ‖ 2 .

Using the same arguments, in view of the second line in (Eq. 12), ‖ u J k + 1 + u G * + ∇ F ̃ ( z * ) ‖ 2 ≤ 1 γ 2 ‖ r k − r * ‖ + ‖ z G k + 1 − z * ‖ 2 ≤ 4 γ 2 ‖ r k − r * ‖ 2 .

Finally, as visible in the first line of (Eq. 16), since r k = z J k + γ k u G k , and using the Moreau identity, we have u G k + 1 = p r o x G * / γ ( 1 γ z J k + u G k ) = p r o x G * / γ ( 1 γ r k ) , so that ‖ u G k + 1 − u G * ‖ 2 ≤ 1 γ 2 ‖ r k − r * ‖ 2 . □

We set Z = X × W , F ̃ , G = R ̃ , J = H ̃ , as defined in Section 2. Doing the substitutions in (Eq. 12), we get the algorithm:

Let s 0 ∈ X and r w 0 ∈ W . For k = 0, 1, … iterate: ⌊ x k + 1 = p r o x γ k R ( s k ) u k + 1 = p r o x H * / ( γ k + 1 η ) K 1 γ k + 1 + 1 γ k x k + 1 − 1 γ k s k − ∇ F ( x k + 1 ) / η − C r w k / ( γ k η ) s k + 1 = x k + 1 − γ k + 1 ∇ F ( x k + 1 ) − γ k + 1 K * u k + 1 r w k + 1 = − γ k + 1 C * u k + 1 .

We can remove the variable r _w and the algorithm becomes: Let s 0 ∈ X and u 0 ∈ U . For k = 0, 1, … iterate: ⌊ x k + 1 = p r o x γ k R ( s k ) u k + 1 = p r o x H * / ( γ k + 1 η ) 1 η K 1 γ k + 1 + 1 γ k x k + 1 − 1 γ k s k − ∇ F ( x k + 1 ) + 1 η C C * u k s k + 1 = x k + 1 − γ k + 1 ∇ F ( x k + 1 ) − γ k + 1 K * u k + 1 .

After replacing CC* by ηI − KK*, the iteration becomes: ⌊ x k + 1 = p r o x γ k R ( s k ) u k + 1 = p r o x H * / ( γ k + 1 η ) u k + 1 η K 1 γ k + 1 + 1 γ k x k + 1 − 1 γ k s k − ∇ F ( x k + 1 ) − K * u k s k + 1 = x k + 1 − γ k + 1 ∇ F ( x k + 1 ) − γ k + 1 K * u k + 1 .

We can change the variables, so that only one call to ∇F and K* appears, which yields the algorithm: Let q 0 ∈ X and u 0 ∈ U . For k = 0, 1, … iterate: ⌊ x k + 1 = p r o x γ k R γ k ( q k − K * u k ) q k + 1 = 1 γ k + 1 x k + 1 − ∇ F ( x k + 1 ) u k + 1 = p r o x H * / ( γ k + 1 η ) u k + 1 η K 1 γ k x k + 1 + q k + 1 − q k .