Three Operator Splitting with Subgradients, Stochastic Gradients, and Adaptive Learning Rates
Year: 2021

Abstract:
Three Operator Splitting (TOS) (Davis & Yin, 2017) can minimize the sum of multiple convex functions effectively when an efficient gradient oracle or proximal operator is available for each term. This requirement often fails in machine learning applications: (i) instead of full gradients only stochastic gradients may be available; and (ii) instead of proximal operators, using subgradients to handle complex penalty functions may be more efficient and realistic. Motivated by these concerns, we analyze three potentially valuable extensions of TOS. The first two permit using subgradients and stochastic gradients, and are shown to ensure a O(1/ √ t) convergence rate. The third extension ADAPTOS endows TOS with adaptive stepsizes. For the important setting of optimizing a convex loss over the intersection of convex sets ADAPTOS attains universal convergence rates, i.e., the rate adapts to the unknown smoothness degree of the objective function. We compare our proposed methods with competing methods on various applications.

text: √ t)
convergence rate. The third extension ADAPTOS endows TOS with adaptive stepsizes. For the important setting of optimizing a convex loss over the intersection of convex sets ADAPTOS attains universal convergence rates, i.e., the rate adapts to the unknown smoothness degree of the objective function. We compare our proposed methods with competing methods on various applications.

heading: 1 Introduction
text: We study convex optimization problems of the form
min x∈Rn φ(x) := f(x) + g(x) + h(x), (1)
where f : Rn → R and g, h : Rn → R ∪ {+∞} are proper, lower semicontinuous and convex functions. Importantly, this template captures constrained problems via indicator functions. To avoid pathological examples, we assume that the relative interiors of dom(f), dom(g) and dom(h) have a nonempty intersection.
Problem (1) is motivated by a number of applications in machine learning, statistics, and signal processing, where the three functions comprising the objective φ model data fitting, structural priors, or decision constraints. Examples include overlapping group lasso (Yuan et al., 2011), isotonic regression (Tibshirani et al., 2011), dispersive sparsity (El Halabi & Cevher, 2015), graph transduction (Shivanna et al., 2015), learning with correlation matrices (Higham & Strabić, 2016), and multidimensional total variation denoising (Barbero & Sra, 2018).
An important technique for addressing composite problems is operator splitting (Bauschke et al., 2011). However, the basic proximal-(sub)gradient method may be unsuitable for Problem (1) since it
Alp Yurtsever and Alex Gu contributed equally to this paper. The paper is based primarily on the work done while Alp Yurtsever was at Massachusetts Institute of Technology.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
requires the prox-operator of g + h, computing which may be vastly more expensive than individual prox-operators of g and h. An elegant, recent method, Three Operator Splitting (TOS, Davis & Yin (2017), see Algorithm 1) offers a practical choice for solving Problem (1) when f is smooth. Importantly, at each iteration, TOS evaluates the gradient of f and the proximal operators of g and h only once. Moreover, composite problems with more than three functions can be reformulated as an instance of Problem (1) in a product-space and solved by using TOS. This is an effective method as long as each function has an efficient gradient oracle or proximal operator (see Section 2).
Unfortunately, TOS is not readily applicable to many optimization problems that arise in machine learning. Most important among those are problems where only access to stochastic gradients is feasible, e.g., when performing large-scale empirical risk minimization and online learning. Moreover, prox-operators for some complex penalty functions are computationally expensive and it may be more efficient to instead use subgradients. For example, proximal operator for the maximum eigenvalue function that appears in dual-form semidefinite programs (e.g., see Section 6.1 in (Ding et al., 2019)) may require computing a full eigendecomposition. In contrast, we can form a subgradient by computing only the top eigenvector via power method or Lanczos algorithm.
Contributions. With the above motivation, this paper contributes three key extensions of TOS. We tackle nonsmoothness in Section 3 and stochasticity in Section 4. These two extensions enable us to use subgradients and stochastic gradients of f (see Section 2 for a comparison with related work), and satisfy a O(1/ √ T ) error bound in function value after T iterations. The third main contribution is ADAPTOS in Section 5. This extension provides an adaptive step-size rule in the spirit of AdaGrad (Duchi et al., 2011; Levy, 2017) for an important subclass of Problem (1). Notably, for optimizing a convex loss over the intersection of two convex sets, ADAPTOS ensures universal convergence rates. That is, ADAPTOS implicitly adapts to the unknown smoothness degree of the problem, and ensures a Õ(1/ √ t) convergence rate when the problem is nonsmooth but the rate improves to Õ(1/t) if the problem is smooth and a solution lies in the relative interior of the feasible set.
In Section 6, we discuss empirical performance of our methods by comparing them against present established methods on various benchmark problems from COPT Library (Pedregosa et al., 2020) including the overlapping group lasso, total variation deblurring, and sparse and low-rank matrix recovery. We also test our methods on nonconvex optimization by training a neural network model. We present more experiments on isotonic regression and portfolio optimization in the supplements.
Notation. We denote a solution of Problem (1) by x? and φ? := φ(x?). The distance between a point x ∈ Rn and a closed and convex set G ⊆ Rn is dist(x,G) := miny∈G ‖x − y‖; the projection of x onto G is given by projG(x) := arg miny∈G ‖x− y‖. The prox-operator of a function g : Rn → R ∪ {+∞} is defined by proxg(x) := arg miny∈Rn{g(y) + 12‖x− y‖
2}. The indicator function of G gives 0 for all x ∈ G and +∞ otherwise. Clearly, the prox-operator of an indicator function is the projection onto the corresponding set.

heading: 2 Background and related work
text: TOS, proposed recently by Davis & Yin (2017), can be seen as a generic extension of various operator splitting schemes, including the forward-backward splitting, Douglas-Rachford splitting, forward-Douglas-Rachford splitting (Briceño-Arias, 2015), and the generalized forward-backward splitting (Raguet et al., 2013). It covers these aforementioned approaches as special instances when the terms f, g and h in Problem (1) are chosen appropriately. Convergence of TOS is well studied when f has Lipschitz continuous gradients. It ensures O(1/t) convergence rate in this setting, see (Davis & Yin, 2017) and (Pedregosa, 2016) for details.
Other related methods that can be used for Problem (1) when f is smooth are the primal-dual hybrid gradient (PDHG) method (Condat, 2013; Vũ, 2013) and the primal-dual three operator splitting methods in (Yan, 2018) and (Salim et al., 2020). These methods can handle a more general template where g or h is composed with a linear map, however, they require f to be smooth. The convergence rate of PDHG is studied in (Chambolle & Pock, 2016).
Nonsmooth setting. We are unaware of any prior result that permits using subgradients in TOS (or in other methods that can use the prox-operator of g and h separately for Problem (1)). The closest match is the proximal subgradient method which applies when h is removed from Problem (1), and it is covered by our nonsmooth TOS as a special case.
Stochastic setting. There are multiple attempts to devise a stochastic TOS in the literature. Yurtsever et al. (2016) studied Problem (1) under the assumption that f is smooth and strongly convex, and an unbiased gradient estimator with bounded variance is available. Their stochastic TOS has a guaranteed O(1/t) convergence rate. In (Cevher et al., 2018), they drop the strong convexity assumption, instead they assume that the variance is summable. They show asymptotic convergence with no guarantees on the rate. Later, Pedregosa et al. (2019) proposed a stochastic variance-reduced TOS and analyzed its non-asymptotic convergence guarantees. Their method gets O(1/t) convergence rate when f is smooth. The rate becomes linear if f is smooth and strongly convex and g (or h) is also smooth. Recently, Yurtsever et al. (2021) studied TOS on problems where f can be nonconvex and showed that the method finds a first-order stationary point with O(1/ 3 √ t) convergence rate under a diminishing variance assumption. They increase the batch size over the iterations to satisfy this assumption.
None of these prior works cover the broad template we consider: f is smooth or Lipschitz continuous and the stochastic first-order oracle has bounded variance. To our knowledge, our paper gives the first analysis for stochastic TOS without strong convexity assumption or variance reduction.
Other related methods are the stochastic PDHG in (Zhao & Cevher, 2018), the decoupling method in (Mishchenko & Richtárik, 2019), the stochastic primal-dual method in (Zhao et al., 2019), and the stochastic primal-dual three operator splitting in (Salim et al., 2020). The method in (Zhao et al., 2019) can be viewed as an extension of stochastic ADMM (Ouyang et al., 2013; Azadi & Sra, 2014) from the sum of two terms to three terms in the objective. Similar to the existing stochastic TOS variants, these methods either assume strong convexity or require variance-reduction.
Adaptive step-sizes. The standard writings of TOS and PDHG require the knowledge of the smoothness constant of f for the step-size. Backtracking line-search strategies (for finding a suitable step-size when the smoothness constant is unknown) are proposed for PDHG in (Malitsky & Pock, 2018) and for TOS in (Pedregosa & Gidel, 2018). These line-search strategies are significantly different than our adaptive learning rate. Importantly, these methods work only when f is smooth. They require extra function evaluations, and are thus not suitable for stochastic optimization. And their goal is to estimate the smoothness constant. In contrast, our goal is to design an algorithm that adapts to the unknown smoothness degree. Our method does not require function evaluations, and it can be used in smooth, nonsmooth, or stochastic settings.
At the heart of our method lie adaptive online learning algorithms (Duchi et al., 2011; Rakhlin & Sridharan, 2013) together with online to offline conversion techniques (Levy, 2017; Cutkosky, 2019). Similar methods appear in the literature for other problem templates with no constraint or a single constraint in (Levy, 2017; Levy et al., 2018; Kavis et al., 2019; Cutkosky, 2019; Bach & Levy, 2019). Our method extends these results to optimization over the intersection of convex sets. When f is nonsmooth, ADAPTOS ensures a Õ(1/ √ t) rate, whereas the rate improves to Õ(1/t) if f is smooth and there is a solution in the relative interior of the feasible set.
TOS for more than three functions. TOS can be used for solving problems with more than three convex functions by a product-space reformulation technique (Briceño-Arias, 2015). Consider
min x∈Rd q∑ i=1 φi(x), (2)
where each component φi : Rd → R∪{+∞} is a proper, lower semicontinuous and convex function. Without loss of generality, suppose φ1, . . . , φp are prox-friendly. Then, we can reformulate (2) in the product-space Rd×(p+1) as
min (x0,x1,...,xp)∈Rd×(p+1) p∑ i=1 φi(xi) + q∑ i=p+1 φi(x0) subject to x0 = x1 = . . . = xp. (3)
This is an instance of Problem (1) with n = d×(p+1) and x = (x0, x1, . . . , xp). We can choose g(x) as the indicator of the equality constraint, f(x) = ∑q i=p+1 φi(x0), and h(x) = ∑p i=1 φi(xi). Then, the (sub)gradient of f is the sum of (sub)gradients of φp+1, . . . , φq; proxg is a mapping that averages x0, x1, . . . , xp; and proxh is the concatenation of the individual prox-operators of φ1, . . . , φp.
To our knowledge, TOS has been studied only for problems with smooth f , and this forces us to assign all nonsmooth components φi in (2) to the proximal term h in (3). In this work, by enabling subgradient steps for nonsmooth f , we provide the flexibility to choose how to process each nonsmooth component φi in (3), either by its proximal operator through h or by its subgradient via f .
Algorithm 1 Three Operator Splitting (TOS)
Input: Initial point y0 ∈ Rn, step-size sequence {γt}Tt=0 for t = 0, 1, 2, . . . , T do zt = proxγtg(yt) Choose an update direction ut ∈ Rn {ut = ∇f(zt) captures the standard version of TOS} xt = proxγth(2zt − yt − γtut) yt+1 = yt − zt + xt end for Return: Ergodic sequence x̄t and z̄t defined in (5)

heading: 3 TOS for Nonsmooth Setting
text: Algorithm 1 presents the generalized TOS for Problem (1). It recovers the standard version in (Davis & Yin, 2017) if we choose ut = ∇f(zt) when f is smooth. For convenience, we define the mapping
TOSγ(y, u) := y − proxγg(y) + proxγh ( 2 · proxγg(y)− y − γu ) (4)
which represents one iteration of Algorithm 1.
The first step of the analysis is the fixed-point characterization of TOS. The following lemma is a straightforward extension of Lemma 2.2 in (Davis & Yin, 2017) to permit subgradients. The proof is similar to (Davis & Yin, 2017), we present it in the supplementary material for completeness. Lemma 1 (Fixed points of TOS). Let γ > 0. Then, there exists a subgradient u ∈ ∂f(proxγg(y)) that satisfies TOSγ(y, u) = y if and only if proxγg(y) is a solution of Problem (1).
When f is Lf -smooth, TOS with ut = ∇f(zt) is known to be an averaged operator1 if γ ∈ (0, 2/Lf ) (see Proposition 2.1 in (Davis & Yin, 2017)) and the analysis in prior work is based on this property. In particular, averagedness implies Fejér monotonicity, i.e., that ‖yt − y?‖ is non-increasing, where y? denotes a fixed point of TOS. However, when f is nonsmooth and ut is replaced with a subgradient, TOS operator is no longer averaged and the standard analysis fails. One of our key observations is that ‖yt − y?‖ remains bounded even-though we loose averagedness and Fejér monotonicity in this setting, see Theorem S.6 in the supplements.
Ergodic sequence. Convergence of operator splitting methods are often given in terms of ergodic (averaged) sequences. This strategy requires maintaining the running averages of zt and xt:
x̄t = 1
t+ 1 t∑ τ=0 xτ and z̄t = 1 t+ 1 t∑ τ=0 zτ . (5)
Clearly, we do not need to store the history of xt and zt to maintain these sequences. In practice, the last iterate often converges faster than the ergodic sequence. We can evaluate the objective function at both points and return the one with the smaller value.
We are ready to present convergence guarantees of TOS for the nonsmooth setting. Theorem 1. Consider Problem (1) and employ TOS (Algorithm 1) with the update directions and step-size chosen as
ut ∈ ∂f(zt) and γt = γ0√ T + 1 for some γ0 > 0, for t = 0, 1, . . . , T . (6)
Assume that ‖ut‖ ≤ Gf for all t. Then, the following guarantees hold:
f(z̄T ) + g(z̄T ) + h(x̄T )− φ? ≤ 1
2 √ T + 1
( D2
γ0 + γ0G
2 f
) (7)
and ‖x̄T − z̄T ‖ ≤ 2
T + 1 (D + γ0Gf ) , where D = max{‖y0 − x?‖, ‖y0 − y?‖}. (8)
1An operator T : Rn → Rn is ω-averaged if ‖Tx− Ty‖2 ≤ ‖x− y‖2 − 1−ω ω ‖(x− Tx)− (y − Ty)‖2 for
some ω ∈ (0, 1) for all x, y ∈ Rn.
Remark 1. The boundedness of subgradients is a standard assumption in nonsmooth optimization. It is equivalent to assuming that f is Gf -Lipschitz continuous on dom(g).
If D and Gf are known, we can optimize the constants in (7) by choosing γ0 = D/Gf . This gives f(z̄T ) + g(z̄T ) + h(x̄T )− φ? ≤ O(DGf/ √ T ) and ‖x̄T − z̄T ‖ ≤ O(D/T ).
Proof sketch. We start by writing the optimality conditions for the proximal steps for zt and xt. Through algebraic modifications and by using convexity of f , g and h, we obtain
f(zt) + g(zt) + h(xt)− φ? ≤ 1
2γ ‖yt − x?‖2 −
1
2γ ‖yt+1 − x?‖2 +
γ 2 ‖ut‖2. (9)
‖ut‖ ≤ Gf by assumption. Then, we average this inequality over t = 0, 1, . . . , T and use Jensen’s inequality to get (7).
The bound in (8) is an immediate consequence of the boundedness of ‖yT+1 − y?‖ that we show in Theorem S.6 in the supplementary material:
‖yT+1 − y?‖ ≤ ‖y0 − y?‖+ 2γ0Gf . (10)
By definition, ‖x̄T − z̄T ‖ = 1T ‖yT+1 − y0‖ ≤ 1 T (‖yT+1 − y?‖+ ‖y? − y0‖).
Theorem 1 does not immediately yield convergence to a solution of Problem (1) because f + g and h are evaluated at different points in (7). Next corollary solves this issue.
Corollary 1. We are interested in two particular cases of Theorem 1:
(i). Suppose h is Gh-Lipschitz continuous. Then,
φ(z̄T )− φ? ≤ 1
2 √ T + 1
( D2
γ0 + γ0G
2 f
) +
2Gh T + 1 (D + γ0Gf ) . (11)
(ii). Suppose h is the indicator function of a convex setH ⊆ Rn. Then,
f(z̄T ) + g(z̄T )− φ? ≤ 1
2 √ T + 1
( D2
γ0 + γ0G
2 f
) (12)
and dist(z̄T ,H) ≤ 2
T + 1 (D + γ0Gf ) . (13)
Proof. (i). Since h is Gh-Lipschitz, φ(z̄T ) ≤ f(z̄T ) + g(z̄T ) + h(x̄T ) +Gh‖x̄T − z̄T ‖. (ii). h(x̄T ) = 0 since x̄T ∈ H. Moreover, dist(z̄T ,H) := infx∈H ‖z̄T − x‖ ≤ ‖z̄T − x̄T ‖.
Remark 2. We fix time horizon T for the ease of analysis and presentation. In practice, we use γt = γ0/ √ t+ 1.
Theorem 1 covers the case in which g is the indicator of a convex set G ⊆ Rn. By definition, z̄T ∈ G and x? ∈ G, hence g(z̄T ) = g(x?) = 0. If both g and h are indicator functions, TOS gives an approximately feasible solution, in G, and close to H. We can also consider a stronger notion of approximate feasibility, measured by dist(z̄T ,G ∩H). However, this requires additional regularity assumptions on G and H to avoid pathological examples, see Lemma 1 in (Hoffmann, 1992) and Definition 2 in (Kundu et al., 2018).
Problem (1) captures unconstrained minimization problems when g = h = 0. Therefore, the convergence rate in Theorem 1 is optimal in the sense that it matches the information theoretical lower bounds for first-order black-box methods, see Section 3.2.1 in (Nesterov, 2003). Remark that the subgradient method can achieve a O(1/t) rate when f is strongly convex. We leave the analysis of TOS for strongly convex nonsmooth f as an open problem.

heading: 4 TOS for Stochastic Setting
text: In this section, we focus on the three-composite stochastic optimization template:
min x∈Rn
φ(x) := f(x) + g(x) + h(x) where f(x) := Eξ f̃(x, ξ) (14)
and ξ is a random variable. The following theorem characterizes the convergence rate of Algorithm 1 for Problem (14). Theorem 2. Consider Problem (14) and employ TOS (Algorithm 1) with a fixed step-size γt = γ = γ0/ √ T + 1 for some γ0 > 0. Suppose we are receiving the update directions ut from an unbiased stochastic first-order oracle with bounded variance, i.e.,
ût := E[ut|zt] ∈ ∂f(zt) and E[‖ut − ût‖2] ≤ σ2 for some σ < +∞. (15) Assume that ‖ût‖ ≤ Gf for all t. Then, the following guarantees hold:
E[f(z̄T ) + g(z̄T ) + h(x̄T )]− φ? ≤ 1
2 √ T + 1
( D2
γ0 + γ0(σ
2 +G2f )
) and (16)
E[‖x̄T − z̄T ‖] ≤ 2
T + 1
( D + γ0 ( Gf + σ
2
)) , where D = max{‖y0 − x?‖, ‖y0 − y?‖}. (17)
Remark 3. Similar rate guarantees hold with some restrictions on the choice of γ0 if we replace bounded subgradients assumption with the smoothness of f . We defer details to the supplements.
If we can estimateD,Gf and σ, then we can optimize the bounds by choosing γ0 ≈ D/max{Gf , σ}. This gives f(z̄T ) + g(z̄T ) + h(x̄T )− φ? ≤ O(Dmax{Gf , σ}/ √ T ) and ‖x̄T − z̄T ‖ ≤ O(D/T ).
Analogous to Corollary 1, from Theorem 2 we can derive convergence guarantees when h is Lipschitz continuous or an indicator function. As in the nonsmooth setting, the rates shown in this section are optimal because Problem (14) covers g(x) = h(x) = 0 as a special case.

heading: 5 TOS with Adaptive Learning Rates
text: In this section, we focus on an important subclass of Problem (1) where g and h are indicator functions of some closed and convex sets:
min x∈Rn
f(x) subject to x ∈ G ∩H. (18)
TOS is effective for Problem (18) when projections onto G andH are easy but the projection onto their intersection is challenging. Particular examples include transportation polytopes, doubly nonnegative matrices, and isotonic regression, among many others.
We propose ADAPTOS with an adaptive step-size in the spirit of adaptive online learning algorithms and online to batch conversion techniques, see (Duchi et al., 2011; Rakhlin & Sridharan, 2013; Levy, 2017; Levy et al., 2018; Cutkosky, 2019; Kavis et al., 2019; Bach & Levy, 2019) and the references therein. ADAPTOS employs the following step-size rule:
γt = α√
β + ∑t−1 τ=0 ‖uτ‖2 for some α, β > 0. (19)
β in the denominator prevents γt to become undefined. If D := ‖y0− x?‖ and Gf are known, theory suggests choosing α = D and β = G2f for a tight upper bound, however, this choice affects only the constants and not the rate of convergence as we demonstrate in the rest of this section. Importantly, we do not assume any prior knowledge on D or Gf . In practice, we often discard β and use γ0 = α at the first iteration.
For ADAPTOS, in addition to (5), we will also use a second ergodic sequence with weighted averaging:
x̃t = 1∑t
τ=0 γτ t∑ τ=0 γτxτ and z̃t = 1∑t τ=0 γτ t∑ τ=0 γτzτ . (20)
This sequence was also considered for TOS with line-search in (Pedregosa & Gidel, 2018).
Theorem 3. Consider Problem (18) and TOS (Algorithm 1) with the update directions ut ∈ ∂f(zt) and the adaptive step-size (19). Assume that ‖ut‖ ≤ Gf for all t. Then, the estimates generated by TOS satisfy
f(z̃t)− f? ≤ Õ (
2αGf√ t+ 1
( D2
4α2 + 1 + Gf√ β
)) and (21)
dist(z̄t,H) ≤ Õ (
2α√ t+ 1
( 1 +
Gf√ β
)) where D = ‖y0 − x?‖. (22)
If D and Gf are known, we can choose α = D and β = G2f . This gives f(z̃t)− f? ≤ Õ(GfD/ √ t)
and dist(z̄t,H) ≤ Õ(D/ √ t).
The next theorem establishes a faster rate for the same algorithm when f is smooth and a solution lies in the interior of the feasible set. Theorem 4. Consider Problem (18) and suppose f is Lf -smooth on G. Use TOS (Algorithm 1) with the update directions ut = ∇f(zt) and the adaptive step-size (19). Assume that ‖ut‖ ≤ Gf for all t. Suppose Problem (18) has a solution in the interior of the feasible set. Then, the estimates generated by TOS satisfy
f(z̄t)− f? ≤ Õ ( 2
t+ 1
( 4α2Lf ( D2 4α2 + 1 + G2f β )2 + α √ β ( D2 4α2 + 1 + G2f β ))) and (23)
dist(z̄t,H) ≤ Õ ( 2α
t+ 1
( D α + 1 + Gf√ β )) where D = ‖y0 − x?‖. (24)
If D and Gf are known, we can choose α = D and β = G2f . This gives f(z̄t)− f? ≤ Õ((LfD2 +GfD)/t) and dist(z̄t,H) ≤ Õ(D/t). Remark 4. When f is smooth, the boundedness assumption ‖ut‖ ≤ Gf holds automatically with Gf ≤ LfDG if G has a bounded diameter DG .
We believe the assumption on the location of the solution is a limitation of the analysis and that the method can achieve fast rates when f is smooth regardless of where the solution lies. Remark that this assumption also appears in (Levy, 2017; Levy et al., 2018).
Following the definition in (Nesterov, 2015), we say that an algorithm is universal if it does not require to know whether the objective is smooth or not yet it implicitly adapts to the smoothness of the objective. ADAPTOS attains universal convergence rates for Problem (18). It converges to a solution with Õ(1/ √ t) rate (in function value) when f is nonsmooth. The rate becomes Õ(1/t) if f is smooth and the solution is in the interior of the feasible set.
Finally, the next theorem shows that ADAPTOS can successfully handle stochastic (sub)gradients. Theorem 5. Consider Problem (18). Use TOS (Algorithm 1) with the update directions ut from an unbiased stochastic subgradient oracle such that E[ut|zt] ∈ ∂f(zt) almost surely. Assume that ‖ut‖ ≤ Gf for all t. Suppose Problem (18) has a solution in the interior of the feasible set. Then, the estimates generated by TOS satisfy
E [ f(z̃t)− f? ] ≤ Õ ( 2αGf√ t+ 1 ( D2 4α2 + 1 + G2f β )) and (25)
E[dist(z̄t,H)] ≤ Õ ( 2α
t+ 1
( D α + 1 + Gf√ β )) where D = ‖y0 − x?‖. (26)

heading: 6 Numerical Experiments
text: This section demonstrates empirical performance of the proposed method on a number of convex optimization problems. We also present an experiment on neural networks. Our experiments are performed in Python 3.7 with Intel Core i9-9820X CPU @ 3.30GHz. We present more experiments on isotonic regression and portfolio optimization in the supplementary materials. The source code for the experiments is available in the supplements.

heading: 6.1 Experiments on Convex Optimization with Smooth f
text: In this subsection, we compare ADAPTOS with TOS, PDHG and their line-search variants TOS-LS and PDHG-LS. Our experiments are based on the benchmarks described in (Pedregosa & Gidel, 2018) and their source code available in COPT Library (Pedregosa et al., 2020) under the new BSD License. We implement ADAPTOS and investigate its performance on three different problems:
. Logistic regression with overlapping group lasso penalty:
min x∈Rn
1
N N∑ i=1 log(1 + exp(−bi〈ai, x〉)) + λ ∑ G∈G √ |G|‖xG‖+ λ ∑ H∈H √ |H|‖xH‖, (27)
where {(a1, b1), . . . , (aN , bN )} is a given set of training examples, G andH are the sets of distinct groups and | · | denotes the cardinality. The model we use (from COPT) considers groups of size 10 with 2 overlapping coefficients. In this experiment, we use the benchmarks on synthetic data (dimensions n = 1002, N = 100) and real-sim dataset (Chang & Lin, 2011) (n = 20958, N = 72309).
. Image recovery with total variation penalty:
min X∈Rm×n ‖Y −A(X)‖2F + λ m∑ i=1 n−1∑ j=1 |Xi,j+1 −Xi,j |+ λ n∑ j=1 m−1∑ i=1 |Xi+1,j −Xi,j |, (28)
where Y is a given blurred image and A : Rm×n → Rm×n is a linear operator (blur kernel). The benchmark in COPT solves this problem for an image of size 153× 115 with a provided blur kernel. . Sparse and low-rank matrix recovery via `1 and nuclear-norm regularizations:
min X∈Rn×n
1
N N∑ i=1 huber(bi − 〈Ai, X〉) + λ‖X‖∗ + λ‖X‖1. (29)
We use huber loss. {(A1, b1), . . . , (AN , bN )} is a given set of measurements and ‖X‖1 is the vector `1-norm of X . The benchmark in COPT considers a symmetric ground truth matrix X\ ∈ R20×20 and noisy synthetic measurements (N = 100) where Ai has Gaussian iid entries. bi = 〈Ai, X\〉+ωi where ωi is generated from a zero-mean unit variance Gaussian distribution.
At each problem, we consider two different values for the regularization parameter λ. We use all methods with their default parameters in the benchmark. For ADAPTOS, we discard β and tune α by trying the powers of 10. See the supplementary material for the behavior of the algorithm with different values of α. Figure 1 shows the results of this experiment. In most cases, the performance of ADAPTOS is between TOS-LS and PDHG-LS. Remark that TOS-LS is using the extra knowledge of the Lipschitz constant of h.

heading: 6.2 Experiments on Convex Optimization with Nonsmooth f
text: We examine the empirical performance of ADAPTOS for nonsmooth problems on an image impainting and denoising task from (Zeng & So, 2018; Yurtsever et al., 2018). We are given an occluded image (i.e., missing some pixels) of size 517× 493, contaminated with salt and pepper noise of 10% density. We use the following template where data fitting is measured in terms of vector `p-norm:
min X∈Rm×n
‖A(X)− Y ‖p subject to ‖X‖∗ ≤ λ, 0 ≤ X ≤ 1, (30)
where Y is the observed noisy image with missing pixels. This is essentially a matrix completion problem, A : Rm×n → Rm×n is a linear map that samples the observed pixels in Y . In particular, we consider (30) with p = 1 and p = 2. The `2-loss is common in practice for matrix completion (often in the least-squares form) but it is not robust against the outliers induced by the salt and pepper noise. `1-loss is known to be more reliable for this task.
The subgradients in both cases have a fixed norm at all points (note that the subgradients are binary valued for `1-loss and unit-norm for `2-loss), hence the analytical and the adaptive step-sizes are same up to a constant factor.
Figure 2 shows the results. The empirical rates for p = 1 roughly match our guarantees in Theorem 1. We observe a locally linear convergence rate when `2-loss is used. Interestingly, the ergodic sequence converges faster than the last iterate for p = 1 but significantly slower for p = 2. The runtime of the two settings are approximately the same, with 67 msec per iteration on average. Despite the slower rates, we found `1-loss more practical on this problem. A low-accuracy solution obtained by 1000 iterations on `1-loss yields a high quality recovery with PSNR 26.21 dB, whereas the PSNR saturates at 21.15 dB for the `2-formulation. See the supplements for the recovered images and more details.

heading: 6.3 An Experiment on Neural Networks
text: In this section, we train a regularized deep neural network to test our methods on nonconvex optimization. We consider a regularized neural network problem formulation in (Scardapane et al., 2017). This problem involves a fully connected neural network with the standard cross-entropy loss function, a ReLu activation for the hidden layers, and the softmax activation for the output layer. Two regularizers are added to this loss function: The first one is the standard `1 regularizer, and the second is the group sparse regularizer where the outgoing connections of each neuron is considered as a group. The goal is to force all outgoing connections from the same neurons to be simultaneously zero, so that we can safely remove the neurons from the network. This is shown as an effective way to obtain compact networks (Scardapane et al., 2017), which is crucial for the deployment of the learned parameters on resource-constrained devices such as smartphones (Blalock et al., 2020).
We reuse the open source implementation (built with Lasagne framework based on Theano) published in (Scardapane et al., 2017) under BSD-2 License. We follow their experimental setup and instructions with MNIST database (LeCun, 1998) containing 70k grayscale images (28×28) of handwritten digits (split 75/25 into train and test partitions). We train a fully connected neural network with 784 input features, three hidden layers (400/300/100) and 10-dimensional output layer. Interested readers can find more details on the implementation in the supplementary material or in (Scardapane et al., 2017).
Scardapane et al. (2017) use SGD and Adam with the subgradient of the overall objective. In contrast, our methods can leverage the prox-operators for the regularizers. Figure 3 compares the performance in terms of two measures: the sparsity of the parameters and the accuracy. On the left side, we see the spectrum of weight and neuron magnitudes. The advantage of using prox-operators is outstanding: More than 93% of the weights are zero and 68% of neurons are inactive when trained with ADAPTOS. In contrast, subgradient based methods can achieve only approximately sparse solutions.
The third and the fourth subplots present the training and test accuracies. Remarkably, ADAPTOS performs better than the state-of-the-art (both in train and test). Unfortunately, we could not achieve the same performance gain in preliminary experiments with more complex models like ResNet (He et al., 2016), where SGD with momentum shines. Interested readers can find the code for these preliminary experiments in the supplements. We leave the technical analysis and a comprehensive examination of ADAPTOS for nonconvex problems to a future work.

heading: 7 Conclusions
text: We studied an extension of TOS that permits subgradients and stochastic gradients instead of the gradient step and established convergence guarantees for this extension. Moreover, we proposed an adaptive step-size rule (ADAPTOS) for the minimization of a convex function over the intersection of two convex sets. ADAPTOS guarantees a nearly optimal Õ(1/ √ t) rate on the baseline setting,
and it enjoys the faster Õ(1/t) rate when the problem is smooth and the solution is in the interior of feasible set. We present numerical experiments on various benchmark problems. The empirical performance of the method is promising.
We conclude with a short list of open questions and follow-up directions: (i) In parallel to the subgradient method, we believe TOS can achieve O(1/t) rate guarantees in the nonsmooth setting if f is strongly convex. The analysis remains open. (ii) The faster rate for ADAPTOS on smooth f requires an extra assumption on the location of the solution. We believe this assumption can be removed, and leave this as an open problem. (iii) We analyzed ADAPTOS only for a specific subclass of Problem (1) in which g and h are indicator functions. Extending this result for the whole class is a valuable question for future study.

heading: Acknowledgements
text: The authors would like to thank Ahmet Alacaoglu for carefully reading and reporting an error in the preprint of this paper.
Alp Yurtsever received support from the Swiss National Science Foundation Early Postdoc.Mobility Fellowship P2ELP2_187955, from the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation, and partial postdoctoral support from the NSF-CAREER grant IIS-1846088. Suvrit Sra acknowledges support from an NSF BIGDATA grant (1741341) and an NSF CAREER grant (1846088).

Review summaries:
This paper provides an extension of the Three Operator Splitting methods for composite optimization, of the form $f + g + h$. Previous work assumed $f$ to be smooth. This work tackles two cases: 1) the case where $f$ is non-smooth, and the case where additionally, only a stochastic estimator of the (sub)gradient of $f$ is available. It can be seen as an extension of (Yurtsever et al 2018) in the non-smooth, non-strongly convex case. The paper provides a theoretical worst-case analysis of convergence, and gives non-asymptotic rates for these two cases. It also proposes an analyses an adaptive variant, where step size depends on the previous observed subgradient magnitudes, a la AdaGrad. It provides experimental validation of the method on a non-smooth problem, and on a non-convex non-smooth problem.