A Simple Practical Accelerated Method for Finite Sums
Abstract
We describe a novel optimization method for finite sums (such as empirical risk minimization problems) building on the recently introduced SAGA method. Our method achieves an accelerated convergence rate on strongly convex smooth problems. Our method has only one parameter (a step size), and is radically simpler than other accelerated methods for finite sums. Additionally it can be applied when the terms are nonsmooth, yielding a method applicable in many areas where operator splitting methods would traditionally be applied.
Introduction
A large body of recent developments in optimization have focused on minimization of convex finite sums of the form:
a very general class of problems including the empirical risk minimization (ERM) framework as a special case. Any function can be written in this form by setting and for , however when each is sufficiently regular in a way that can be made precise, it is possible to optimize such sums more efficiently than by treating them as black box functions.
In most cases recently developed methods such as SAG (Schmidt et al., 2013) can find an minimum faster than either stochastic gradient descent or accelerated blackbox approaches, both in theory and in practice. We call this class of methods fast incremental gradient methods (FIG).
FIG methods are randomized methods similar to SGD, however unlike SGD they are able to achieve linear convergence rates under Lipschitzsmooth and strong convexity conditions (Mairal, 2014; Defazio et al., 2014b; Johnson and Zhang, 2013; Konečný and Richtárik, 2013). The linear rate in the first wave of FIG methods directly depended on the condition number () of the problem, whereas recently several methods have been developed that depend on the squareroot of the condition number (Lan and Zhou, 2015; Lin et al., 2015; ShalevShwartz and Zhang, 2013c; Nitanda, 2014). Analogous to the blackbox case, these methods are known as accelerated methods.
In this work we develop another accelerated method, which is conceptually simpler and requires less tuning than existing accelerated methods. The method we give is a primal approach, however it makes use of a proximal operator oracle for each instead of a gradient oracle, unlike other primal approaches. The proximal operator is also used by dual methods such as some variants of SDCA (ShalevShwartz and Zhang, 2013a).
1 Algorithm
Pick some starting point and step size . Initialize each where is any gradient/subgradient at .
Then at step :

Pick index from to uniformly at random.

Update :

Update the gradient table: Set , and leave the rest of the entries unchanged ( for ).
Our algorithm’s main step makes use of the proximal operator for a randomly chosen . For convenience, we use the following compact notation:
This proximal operator can be computed efficiently or in closed form in many cases, see Section 4 for details. Like SAGA, we also maintain a table of gradients , one for each function . We denote the state of at the end of step by . The iterate (our guess at the solution) at the end of step is denoted The starting iterate may be chosen arbitrarily.
The full algorithm is given as Algorithm 1. The sum of gradients can be cached and updated efficiently at each step, and in most cases instead of storing a full vector for each , only a single real value needs to be stored. This is the case for linear regression or binary classification with logistic loss or hinge loss, in precisely the same way as for standard SAGA. A discussion of further implementation details is given in Section 4.
With step size
the expected convergence rate in terms of squared distance to the solution is given by:
when each is smooth and strongly convex. See Nesterov (1998) for definitions of these conditions. Using bigO notation, the number of steps required to reduce the distance to the solution by a factor is:
as . This rate matches the lower bound known for this problem (Lan and Zhou, 2015) under the gradient oracle. We conjecture that this rate is optimal under the proximal operator oracle as well. Unlike other accelerated approaches though, we have only a single tunable parameter (the step size ), and the algorithm doesn’t need knowledge of or except for their appearance in the step size.
Compared to the rate for SAGA and other nonaccelerated FIG methods, accelerated FIG methods are significantly faster when is small compared to , however for the performance is essentially the same. All known FIG methods hit a kind of wall at , where they decrease the error at each step by no more than . Indeed, when the problem is so well conditioned so as to be easy for any FIG method to solve it efficiently. This is sometimes called the big data setting (Defazio et al., 2014b).
Our convergence rate can also be compared to that of optimal firstorder black box methods, which have rates of the form per epoch equivalent. We are able to achieve a speedup on a perepoch basis, for not too large. Of course, all of the mentioned rates are significantly better than the rate of gradient descent.
For nonsmooth but strongly convex problems, we prove a type rate under a standard iterate averaging scheme. This rate does not require the use of decreasing step sizes, so our algorithm requires less tuning than other primal approaches on nonsmooth problems.
2 Relation to other approaches
Our method is most closely related to the SAGA method. To make the relation clear, we may write our method’s main step as:
whereas SAGA has a step of the form:
The difference is the point at which the gradient of is evaluated at. The proximal operator has the effect of evaluating the gradient at instead of . While a small difference on the surface, this change has profound effects. It allows the method to be applied directly to nonsmooth problems using fixed step sizes, a property not shared by SAGA or other primal FIG methods. Additionally, it allows for much larger step sizes to be used, which is why the method is able to achieve an accelerated rate.
It is also illustrative to look at how the methods behave at . SAGA degenerates into regular gradient descent, whereas our method becomes the proximalpoint method (Rockafellar, 1976):
The proximal point method has quite remarkable properties. For strongly convex problems, it converges for any at a linear rate. The downside being the inherent difficulty of evaluating the proximal operator. For the case, if each term is an indicator function for a convex set, our algorithm matches Dykstra’s projection algorithm if we take and use cyclic instead of random steps.
Accelerated incremental gradient methods
Several acceleration schemes have been recently developed as extensions of nonaccelerated FIG methods. The earliest approach developed was the ASDCA algorithm (ShalevShwartz and Zhang, 2013b, c). The general approach of applying the proximalpoint method as the outerloop of a doubleloop scheme has been dubbed the Catalyst algorithm Lin et al. (2015). It can be applied to accelerate any FIG method. Recently a very interesting primaldual approach has been proposed by Lan and Zhou (2015). All of the prior accelerated methods are significantly more complex than the approach we propose, and have more complex proofs.
3 Theory
3.1 Proximal operator bounds
In this section we rehash some simple bounds from proximal operator theory that we will use in this work. Define the shorthand , and let , so that . Note that is a subgradient of at the point . This relation is known as the optimality condition of the proximal operator. Note that proofs for the following two propositions are in the supplementary material.
Proposition 1.
(Strengthening firm nonexpansiveness under strong convexity) For any , and any convex function with strong convexity constant ,
In operator theory this property is known as cocoerciveness of .
Proposition 2.
(Moreau decomposition) For any , and any convex function with Fenchel conjugate :
(1) 
Recall our definition of also. After combining, the following relation thus holds between the proximal operator of the conjugate and :
(2) 
Theorem 3.
For any , and any convex smooth function :
Proof.
We will apply cocoerciveness of the proximal operator of as it appears in the decomposition. Note that Lsmoothness of implies strong convexity of . In particular we apply it to the points and :
Pulling from the right side of the inner product out, and plugging in Equation 9, gives the result. ∎
3.2 Notation
Let be the unique minimizer (due to strong convexity) of . In addition to the notation used in the description of the algorithm, we also fix a set of subgradients , one for each of at , chosen such that . We also define Note that at the solution , we want to apply a proximal step for component of the form:
Notation  Description  Additional relation 

Current iterate at step  
Solution  
Step size  
Shorthand in results for generic  
Proximal operator of at  
A stored subgradient of as seen at step  
A subgradient of at  
Chosen component index (random variable)  
Lemma 4.
(Technical lemma needed by main proof) Under Algorithm 1, taking the expectation over the random choice of , conditioning on and each , allows us to bound the following inner product at step :
The proof is in the supplementary material.
3.3 Main result
Theorem 5.
(single step Lyapunov descent) We define the Lyapunov function of our algorithm (PointSAGA) at step as:
for . Then using step size
when each is smooth and strongly convex and . This is the same Lyapunov function as used by Hofmann et al. (2015).
Proof.
Term 1 of is straightforward to simplify:
For term of we start by applying cocoerciveness (Theorem 11):
Now we add and subtract
where we have pulled out the quadratic term by using (we can take the expectation since the left hand side of the inner product doesn’t depend on ). We now expand further:
(3) 
We further split the left side of the inner product to give two separate inner products:
(4) 
The first inner product in Equation 4 is the quantity we bounded in Lemma 8 by . The second inner product in Equation 4, can be simplified using Theorem 3 (note the right side of the inner product is equal to ):
Combing these gives the following bound on :
Define , where . Now we multiply the above inequality through by and combine with the rest of the Lyapunov function, giving:
We want an convergence rate, so we pull out the required terms:
Now to complete the proof we note that and
Corollary 6.
Theorem 7.
(Nonsmooth case) Suppose each is strongly convex, and . Then after iterations of PointSAGA with step size :
where The proof of this theorem is included in the supplementary material.
4 Implementation
Care must be taken for efficient implementation, particularly in the sparse gradient case. We discuss the key points below. A fast Cython implementation is available on the author’s website incorporating these techniques.
 Proximal operators

For the most common binary classification and regression methods, implementing the proximal operator is straightforward. We include details of the computation of the proximal operators for the hinge, square and logistic losses in the supplementary material. The logistic loss does not have a closed form proximal operator, however it may be computed very efficiently in practice using Newton’s method on a 1D subproblem. For problems of a nontrivial dimensionality the cost of the dot products in the main step is much greater than the cost of the proximal operator evaluation. We also detail how to handle a quadratic regularizer within each term’s prox operator, which has a closed form in terms of the unregularized prox operator.
 Initialization

Instead of setting before commencing the algorithm, we recommend using instead. This avoids the cost of a initial pass over the data. In practical effect this is similar to the SDCA initialization of each dual variable to 0.
5 Experiments
We tested our algorithm which we call PointSAGA against SAGA (Defazio et al., 2014a), SDCA (ShalevShwartz and Zhang, 2013a), Pegasos/SGD (ShalevShwartz et al., 2011) and the catalyst acceleration scheme (Lin et al., 2015). SDCA was chosen as the inner algorithm for the catalyst scheme as it doesn’t require a stepsize, making it the most practical of the variants. Catalyst applied to SDCA is essentially the same algorithm as proposed in ShalevShwartz and Zhang (2013c). A single inner epoch was used for each SDCA invocation. Accelerated MISO as well as the primaldual FIG method (Lan and Zhou, 2015) were excluded as we wanted to test on sparse problems and they are not designed to take advantage of sparsity. The stepsize parameter for each method ( for catalystSDCA) was chosen using a grid search of powers of . The step size that gives the lowest error at the final epoch is used for each method.
We selected a set of commonly used datasets from the LIBSVM repository (Chang and Lin, 2011). The prescaled versions were used when available. Logistic regression with regularization was applied to each problem. The regularization constant for each problem was set by hand to ensure was not in the big data regime ; as noted above, all the methods perform essentially the same when . The constant used is noted beneath each plot. Open source code to exactly replicate the experimental results is available at https://github.com/adefazio/pointsaga.
Algorithm scaling with respect to
The key property that distinguishes accelerated FIG methods from their nonaccelerated counterparts is their performance scaling with respect to the dataset size. For large datasets on wellconditioned problems we expect from the theory to see little difference between the methods. To this end, we ran experiments including versions of the datasets subsampled randomly without replacement in 10% and 5% increments, in order to show the scaling with empirically. The same amount of regularization was used for each subset.
Figure 1 shows the function value suboptimality for each datasetsubset combination. We see that in general accelerated methods dominate the performance of their nonaccelerated counterparts. Both SDCA and SAGA are much slower on some datasets comparatively than others. For example, SDCA is very slow on the 5 and 10% COVTYPE datasets, whereas both SAGA and SDCA are much slower than the accelerated methods on the AUSTRALIAN dataset. These differences reflect known properties of the two methods. SAGA is able to adapt to inherent strong convexity while SDCA can be faster on very wellconditioned problems.
There is no clear winner between the two accelerated methods, each gives excellent results on each problem. The Pegasos (stochastic gradient descent) algorithm with its slower than linear rate is a clear loser on each problem, almost appearing as an almost horizontal line on the log scale of these plots.




Nonsmooth problems
We also tested the RCV1 dataset on the hinge loss. In general we did not expect an accelerated rate for this problem, and indeed we observe that PointSAGA is roughly as fast as SDCA across the different dataset sizes.
References
 Chang and Lin [2011] ChihChung Chang and ChihJen Lin. Libsvm : a library for support vector machines. ACM Transactions on Intelligent Systems and Technology, 2:27:1–27:27, 2011.
 Defazio et al. [2014a] Aaron Defazio, Francis Bach, and Simon LacosteJulien. Saga: A fast incremental gradient method with support for nonstrongly convex composite objectives. Advances in Neural Information Processing Systems 27 (NIPS 2014), 2014a.
 Defazio et al. [2014b] Aaron Defazio, Tiberio Caetano, and Justin Domke. Finito: A faster, permutable incremental gradient method for big data problems. Proceedings of the 31st International Conference on Machine Learning, 2014b.
 Hofmann et al. [2015] Thomas Hofmann, Aurelien Lucchi, Simon LacosteJulien, and Brian McWilliams. Variance reduced stochastic gradient descent with neighbors. In C. Cortes, N.D. Lawrence, D.D. Lee, M. Sugiyama, and R. Garnett, editors, Advances in Neural Information Processing Systems 28, pages 2296–2304. Curran Associates, Inc., 2015.
 Johnson and Zhang [2013] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. NIPS, 2013.
 Konečný and Richtárik [2013] Jakub Konečný and Peter Richtárik. SemiStochastic Gradient Descent Methods. ArXiv eprints, December 2013.
 Lan and Zhou [2015] G. Lan and Y. Zhou. An optimal randomized incremental gradient method. ArXiv eprints, July 2015.
 Lin et al. [2015] Hongzhou Lin, Julien Mairal, and Zaid Harchaoui. A universal catalyst for firstorder optimization. In C. Cortes, N.D. Lawrence, D.D. Lee, M. Sugiyama, and R. Garnett, editors, Advances in Neural Information Processing Systems 28, pages 3366–3374. Curran Associates, Inc., 2015.
 Mairal [2014] Julien Mairal. Incremental majorizationminimization optimization with application to largescale machine learning. Technical report, INRIA Grenoble RhôneAlpes / LJK Laboratoire Jean Kuntzmann, 2014.
 Nesterov [1998] Yu. Nesterov. Introductory Lectures On Convex Programming. Springer, 1998.
 Nitanda [2014] Atsushi Nitanda. Stochastic proximal gradient descent with acceleration techniques. In Z. Ghahramani, M. Welling, C. Cortes, N.D. Lawrence, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems 27, pages 1574–1582. Curran Associates, Inc., 2014.
 Rockafellar [1976] R Tyrrell Rockafellar. Monotone operators and the proximal point algorithm. SIAM journal on control and optimization, 14(5):877–898, 1976.
 Schmidt et al. [2013] Mark Schmidt, Nicolas Le Roux, and Francis Bach. Minimizing finite sums with the stochastic average gradient. Technical report, INRIA, 2013.
 ShalevShwartz and Zhang [2013a] Shai ShalevShwartz and Tong Zhang. Stochastic dual coordinate ascent methods for regularized loss minimization. JMLR, 2013a.
 ShalevShwartz and Zhang [2013b] Shai ShalevShwartz and Tong Zhang. Accelerated minibatch stochastic dual coordinate ascent. In C.J.C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems 26, pages 378–385. Curran Associates, Inc., 2013b.
 ShalevShwartz and Zhang [2013c] Shai ShalevShwartz and Tong Zhang. Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization. Technical report, The Hebrew University, Jerusalem and Rutgers University, NJ, USA, 2013c.
 ShalevShwartz et al. [2011] Shai ShalevShwartz, Yoram Singer, Nathan Srebro, and Andrew Cotter. Pegasos: Primal estimated subgradient solver for svm. Mathematical programming, 127(1):3–30, 2011.
Appendix A Proximal operators
For the most common binary classification and regression methods, implementing the proximal operator is straightforward. In this section let be the label or target for regression, and the data instance vector. We assume for binary classification that .
Hinge loss:
The proximal operator has a closed form expression:
where:
Logistic loss:
There is no closed form expression, however it can be computed very efficiently using Newton iteration, since it can be reduced to a 1D minimization problem. In particular, let , , and . Then iterate until convergence:
The prox operator is then . Three iterations are generally enough, but illconditioned problems or large step sizes may require up to 12. Correct initialization is important, as it will diverge when initialized with a point on the opposite side of 0 from the solution.
Squared loss:
Let and . Define:
Then
L2 regularization
Including a regularizer within each , i.e. can be done using the proximal operator of . Define the scaling factor:
Then .
Appendix B Proofs
Lemma 8.
Under Algorithm 1, taking the expectation over the random choice of , conditioning on and each , allows us to bound the following inner product at step :
Proof.
We start by splitting on the right hand side of the inner product:
(5) 
The first inner product has expectation on the left hand side (Recall that ), so it’s simply 0 in expectation (we may take expectation on the left since the right doesn’t depend on ). The second inner product is the same on both sides, so we may convert it to a normsquared term. So we have:
The inequality used is just an application of the variance formula ∎
Corollary 9.
Chaining the main theorem gives a convergence rate for pointsaga at step under the constants given in of:
if each is smooth and strongly convex.
Proof.
First we simplify using and use Lipschitz smoothness:
Now recall that the main theorem gives a bound where the expectation is conditional on and each from step , taking expectation over the randomness in the choice of . We can further take expectation with respect to and each , giving the unconditional bound:
Chaining over gives the result.∎
Theorem 10.
Suppose each is strongly convex, and . Then after iterations of PointSAGA with step size :
where
Proof.
Recall the bound on the Lyapunov function established in the main theorem:
In the nonsmooth case this holds with . In particular, if we take , then:
Recall that this expectation is (implicitly) conditional on and each from step , Taking expectation over the randomness in the choice of . We can further take expectation with respect to and each , and negate the inequality, giving the unconditional bound:
We now sum this over :
We can drop the since it is always negative. Dividing through by :
Now using Jensen’s inequality on the left gives:
where Now we plug in with :
Now we plug in the bounds in terms of and :
In order to balance the terms on the right, we need:
So we can take , giving a rate of:
∎
Appendix C Proximal operator bounds with proofs
In this section we prove some simple bounds from proximal operator theory that we will use in this work. Define the shorthand , and let , so that . Note that is a subgradient of at the point . This relation is known as the optimality condition of the proximal operator.
We will also use a few standard convexity bounds without proof. Let be a convex function with strong convexity constant and Lipschitz smoothness constant . Let be the minimizer of , then for any :
(6) 
(7) 
Proposition 11.
(Firm nonexpansiveness) For any , and any convex function with strong convexity constant ,
Proof.
Using strong convexity of we apply Equation 6 at the (sub)gradients and , and their corresponding points and :
We now multiply both sides by , then add to both sides: