Large Scale Machine learning Optimization through Stochastic Average Gradient
We propose the stochastic average gradient (SAG) method for optimizing the sum of a finite number of smooth convex functions. Like stochastic gradient (SG) methods, the SAG method's iteration cost is independent of the number of terms in the sum. However, by incorporating a memory of previous gradient values the SAG method achieves a faster convergence rate than black-box SG methods. The convergence rate is improved from O(1/k^1/2) to O(1/k) in general, and when the sum is strongly-convex the convergence rate is improved from the sub-linear O(1/k) to a linear convergence rate of the form O(p^k) for p 1. Further, in many cases the convergence rate of the new method is also faster than black-box deterministic gradient methods, in terms of the number of gradient evaluations. Numerical experiments indicate that the new algorithm often dramatically outperforms existing SG and deterministic gradient methods, and that the performance may be further improved through the use of non-uniform sampling strategies.READ FULL TEXT VIEW PDF
In this paper, we propose a couple of new Stochastic Strictly Contractiv...
Recent advances in optimization theory have shown that smooth strongly c...
Black-box optimization is primarily important for many compute-intensive...
We propose and analyze a new stochastic gradient method, which we call
We address the problem of minimizing a smooth function f^0(x) over a
Structured problems arise in many applications. To solve these problems,...
We apply stochastic average gradient (SAG) algorithms for training
Large Scale Machine learning Optimization through Stochastic Average Gradient
A plethora of the optimization problems arising in practice involve computing a minimizer of a finite sum of functions measuring misfit over a large number of data points. A classical example is least-squares regression,
where the and
are the data samples associated with a regression problem. Another important example is logistic regression,
where the and are the data samples associated with a binary classification problem. A key challenge arising in modern applications is that the number of data points (also known as training examples) can be extremely large, while there is often a large amount of redundancy between examples. The most wildly successful class of algorithms for taking advantage of the sum structure for problems where is very large are stochastic gradient (SG) methods (Robbins and Monro, 1951; Bottou and LeCun, 2003). Although the theory behind SG methods allows them to be applied more generally, SG methods are often used to solve the problem of optimizing a finite sample average,
In this work, we focus on such finite data problems where each is smooth and convex.
In addition to this basic setting, we will also be interested in cases where the sum has the additional property that it is strongly-convex. This often arises due to the use of a strongly-convex regularizer such as the squared -norm, resulting in problems of the form
where each is a data-misfit function (as in least-squares and logistic regression) and the positive scalar controls the strength of the regularization. These problems can be put in the framework of (1) by using the choice
The resulting function
will be strongly-convex provided that the individual loss functionsare convex. An extensive list of convex loss functions used in a statistical data-fitting context is given by Teo et al. (2007), and non-smooth loss functions (or regularizers) can also be put in this framework by using smooth approximations (for example, see Nesterov (2005)).
where is the step size on iteration . Assuming that a minimizer exists, then under standard assumptions the sub-optimality achieved on iteration of the FG method with a constant step size is given by
when is convex (see Nesterov, 2004, Corollary 2.1.2). This results in a sublinear convergence rate. When is strongly-convex, the error also satisfies
for some which depends on the condition number of (see Nesterov, 2004, Theorem 2.1.15). This results in a linear convergence rate, which is also known as a geometric or exponential rate because the error is cut by a fixed fraction on each iteration. Unfortunately, the FG method can be unappealing when is large because its iteration cost scales linearly in .
The main appeal of SG methods is that they have an iteration cost which is independent of , making them suited for modern problems where may be very large. The basic SG method for optimizing (1) uses iterations of the form
where at each iteration an index is sampled uniformly from the set . The randomly chosen gradient
yields an unbiased estimate of the true gradientand one can show under standard assumptions (see Nemirovski et al., 2009) that, for a suitably chosen decreasing step-size sequence , the SG iterations have an expected sub-optimality for convex objectives of
and an expected sub-optimality for strongly-convex objectives of
In these rates, the expectations are taken with respect to the selection of the variables. These sublinear rates are slower than the corresponding rates for the FG method, and under certain assumptions these convergence rates are optimal in a model of computation where the algorithm only accesses the function through unbiased measurements of its objective and gradient (see Nemirovski and Yudin (1983); Nemirovski et al. (2009); Agarwal et al. (2012)). Thus, we should not expect to be able to obtain the convergence rates of the FG method if the algorithm only relies on unbiased gradient measurements. Nevertheless, by using the stronger assumption that the functions are sampled from a finite dataset, in this paper we show that we can achieve the convergence rates of FG methods while preserving the iteration complexity of SG methods.
The primary contribution of this work is the analysis of a new algorithm that we call the stochastic average gradient (SAG) method, a randomized variant of the incremental aggregated gradient (IAG) method of Blatt et al. (2007). The SAG method has the low iteration cost of SG methods, but achieves the convergence rates stated above for the FG method. The SAG iterations take the form
where at each iteration a random index is selected and we set
That is, like the FG method, the step incorporates a gradient with respect to each function. But, like the SG method, each iteration only computes the gradient with respect to a single example and the cost of the iterations is independent of . Despite the low cost of the SAG iterations, we show in this paper that with a constant step-size the SAG iterations have an convergence rate for convex objectives and a linear convergence rate for strongly-convex objectives, like the FG method. That is, by having access to and by keeping a memory of the most recent gradient value computed for each index , this iteration achieves a faster convergence rate than is possible for standard SG methods. Further, in terms of effective passes through the data, we will also see that for many problems the convergence rate of the SAG method is also faster than is possible for standard FG methods.
The next section reviews several closely-related algorithms from the literature, including previous attempts to combine the appealing aspects of FG and SG methods. However, despite years of extensive research on SG methods, with a significant portion of the applications focusing on finite datasets, we believe that this is the first general method that achieves the convergence rates of FG methods while preserving the iteration cost of standard SG methods. Section 3 states the (standard) assumptions underlying our analysis and gives our convergence rate results. Section 4 discusses practical implementation issues including how we adaptively set the step size and how we can reduce the storage cost needed by the algorithm. For example, we can reduce the memory requirements from to in the common scenario where each only depends on a linear function of , as in least-squares and logistic regression. Section 5 presents a numerical comparison of an implementation based on SAG to competitive SG and FG methods, indicating that the method may be very useful for problems where we can only afford to do a few passes through a data set.
A preliminary conference version of this work appears in Le Roux et al. (2012), and we extend this work in various ways. Most notably, the analysis in the prior work focuses only on showing linear convergence rates in the strongly-convex case while the present work also gives an convergence rate for the general convex case. In the prior work we show (Proposition 1) that a small step-size gives a slow linear convergence rate (comparable to the rate of FG methods in terms of effective passes through the data), while we also show (Proposition 2) that a much larger step-size yields a much faster convergence rate, but this requires that is sufficiently large compared to the condition number of the problem. In the present work (Section 3) our analysis yields a very fast convergence rate using a large step-size (Theorem 1), even when this condition required by the prior work is not satisfied. Surprisingly, for ill-conditioned problems our new analysis shows that using SAG iterations can be nearly times as fast as the standard gradient method. Beyond this significantly strengthened result, in this work we also argue that yet-faster convergence rates may be achieved by non-uniform sampling (Section 4.8) and present numerical results showing that this can lead to drastically improved performance (Section 5.5).
There are a large variety of approaches available to accelerate the convergence of SG methods, and a full review of this immense literature would be outside the scope of this work. Below, we comment on the relationships between the new method and several of the most closely-related ideas.
Momentum: SG methods that incorporate a momentum term use iterations of the form
see Tseng (1998). It is common to set all for some constant , and in this case we can rewrite the SG with momentum method as
We can re-write the SAG updates (5) in a similar form as
where the selection function is equal to if corresponds to the last iteration where and is set to otherwise. Thus, momentum uses a geometric weighting of previous gradients while the SAG iterations select and average the most recent evaluation of each previous gradient. While momentum can lead to improved practical performance, it still requires the use of a decreasing sequence of step sizes and is not known to lead to a faster convergence rate.
Gradient Averaging: Closely related to momentum is using the sample average of all previous gradients,
which is similar to the SAG iteration in the form (5) but where all previous gradients are used. This approach is used in the dual averaging method of Nesterov (2009) and, while this averaging procedure and its variants lead to convergence for a constant step size and can improve the constants in the convergence rate (Xiao, 2010), it does not improve on the sublinear convergence rates for SG methods.
Iterate Averaging: Rather than averaging the gradients, some authors propose to perform the basic SG iteration but use an average over certain values as the final estimator. With a suitable choice of step-sizes, this gives the same asymptotic efficiency as Newton-like second-order SG methods and also leads to increased robustness of the convergence rate to the exact sequence of step sizes (Polyak and Juditsky, 1992; Bach and Moulines, 2011). Baher’s method (Kushner and Yin, 2003, §1.3.4) combines gradient averaging with online iterate averaging and also displays appealing asymptotic properties. Several authors have recently shown that suitable iterate averaging schemes obtain an rate for strongly-convex optimization even for non-smooth objectives (Hazan and Kale, 2011; Rakhlin et al., 2012). However, none of these methods improve on the and rates for SG methods.
Stochastic versions of FG methods: Various options are available to accelerate the convergence of the FG method for smooth functions, such as the accelerated full gradient (AFG) method of Nesterov (1983), as well as classical techniques based on quadratic approximations such as diagonally-scaled FG methods, non-linear conjugate gradient, quasi-Newton, and Hessian-free Newton methods (see Nocedal and Wright, 2006). There has been a substantial amount of work on developing stochastic variants of these algorithms, with several of the notable recent examples including Bordes et al. (2009); Hu et al. (2009); Sunehag et al. (2009); Ghadimi and Lan (2010); Martens (2010) and Xiao (2010). Duchi et al. (2011) have recently shown an improved regret bound using a diagonal scaling that takes into account previous gradient magnitudes. Alternately, if we split the convergence rate into a deterministic and stochastic part, these methods can improve the dependency of the convergence rate of the deterministic part (Hu et al., 2009; Ghadimi and Lan, 2010; Xiao, 2010). However, we are not aware of any existing method of this flavor that improves on the and dependencies on the stochastic part. Further, many of these methods typically require carefully setting parameters (beyond the step size) and often aren’t able to take advantage of sparsity in the gradients .
Constant step size: If the SG iterations are used for strongly-convex optimization with a constant step size (rather than a decreasing sequence), then Nedic and Bertsekas (2000, Proposition 3.4) showed that the convergence rate of the method can be split into two parts. The first part depends on and converges linearly to . The second part is independent of and does not converge to . Thus, with a constant step size, the SG iterations have a linear convergence rate up to some tolerance, and in general after this point the iterations do not make further progress. Indeed, up until the recent work of Bach and Moulines (2013), convergence of the basic SG method with a constant step size had only been shown for the strongly-convex quadratic case (with averaging of the iterates) (Polyak and Juditsky, 1992), or under extremely strong assumptions about the relationship between the functions (Solodov, 1998). This contrasts with the method we present in this work which converges to the optimal solution using a constant step size and does so with a linear rate (without additional assumptions).
Accelerated methods: Accelerated SG methods, which despite their name are not related to the aforementioned AFG method, take advantage of the fast convergence rate of SG methods with a constant step size. In particular, accelerated SG methods use a constant step size by default, and only decrease the step size on iterations where the inner-product between successive gradient estimates is negative (Kesten, 1958; Delyon and Juditsky, 1993). This leads to convergence of the method and allows it to potentially achieve periods of faster convergence where the step size stays constant. However, the overall convergence rate of the method is not improved.
Hybrid Methods: Some authors have proposed variants of the SG method for problems of the form (1) that seek to gradually transform the iterates into the FG method in order to achieve a faster convergence rate. Bertsekas (1997) proposes to go through the data cyclically with a specialized weighting that allows the method to achieve a linear convergence rate for strongly-convex quadratic functions. However, the weighting is numerically unstable and the linear convergence rate presented treats full passes through the data as iterations. A related strategy is to group the functions into ‘batches’ of increasing size and perform SG iterations on the batches. Friedlander and Schmidt (2012) give conditions under which this strategy achieves the and convergence rates of FG methods. However, in both cases the iterations that achieve the faster rates have a cost that is not independent of , as opposed to SAG iterations.
Incremental Aggregated Gradient: Blatt et al. (2007) present the most closely-related algorithm to the SAG algorithm, the IAG method. This method is identical to the SAG iteration (5), but uses a cyclic choice of rather than sampling the values. This distinction has several important consequences. In particular, Blatt et al. are only able to show that the convergence rate is linear for strongly-convex quadratic functions (without deriving an explicit rate), and their analysis treats full passes through the data as iterations. Using a non-trivial extension of their analysis and a novel proof technique involving bounding the gradient and iterates simultaneously by a Lyapunov potential function, in this work we give an rate for general convex functions and an explicit linear convergence rate for general strongly-convex functions using the SAG iterations that only examine a single function. Further, as our analysis and experiments show, the SAG iterations allow a much larger step size than is required for convergence of the IAG method. This leads to more robustness to the selection of the step size and also, if suitably chosen, leads to a faster convergence rate and substantially improved practical performance. This shows that the simple change (random selection vs. cycling) can dramatically improve optimization performance.
Special Problem Classes: For restricted classes of problems, several recent works have shown faster convergence rates of methods that operate only on a single function . For example, Strohmer and Vershynin (2009) show that the randomized Kaczmarz method with a particular sampling scheme achieves a linear convergence rate for the problem of solving consistent linear systems. We show in Schmidt and Le Roux (2013) that the SG method with a constant step-size has the and convergence rates of FG methods if is bounded by a linear function of for all and . This is the strong condition required by Solodov (1998) to show convergence of the SG method with a constant step size. Since the first version of this work was released, Shalev-Schwartz and Zhang (2013b) have shown that a linear convergence rate can be obtained for linearly-parameterized loss functions with -regularization using a stochastic (Fenchel) dual coordinate optimization method with an exact line-search. Note that in this setting, there is one dual variable associated with each example so the iteration cost is independent of . Our experiments show that although this method performs as well as the SAG method for certain problems, on others it performs poorly. The related work of Collins et al. (2008) shows, again for a particular -regularized class of models, that when the dual problem satisfies certain properties a linear convergence rate in the dual objective value can be obtained using a dual block-coordinate exponentiated-gradient algorithm, where each block corresponds to the dual variables associated with a given function . More recently, Shalev-Schwartz and Zhang (2013a) have more recently shown that a linear convergence rate can be obtained under less-stringent conditions on the step size and under general strongly-convex regularizers, but it is unclear whether these dual ascent methods can be used to obtain a faster convergence rate without an explicit strongly-convex regularizer. Recently, Bach and Moulines (2013) have shown that for least-squares regression an convergence rate can be achieved by an averaged SG method without strong-convexity. Bach and Moulines (2013) also give a correction strategy that allows an rate to be achieved for general convex functions satisfying a self-concordant property (which includes generalized linear models like logistic regression). Unlike these previous works, our analysis in the next section applies to general that satisfy standard assumptions, and only requires gradient evaluations of the functions rather than dual block-coordinate steps.
Incremental Surrogate Optimization: Since the first version of this work was published, Mairal (2013) has considered an algorithm that is closely-related to SAG, but in the more general proximal-gradient setting that allows simple constraints and simple non-smooth terms. Specialized to the smooth and unconstrained setting, the incremental variant of this algorithm can be written in the form
where is defined as in the SAG algorithm 6, and
is the parameter vector used to compute the corresponding. Thus, instead of applying the SAG step to it applies the step to the average of the previous values used to form the current variables. Mairal (2013) shows that this algorithm also achieves an and a linear convergence rate. However, the linear convergence rate obtained is substantially slower than the convergence rate obtained for SAG. In particular, the rate is more comparable to the rate of the basic FG method in terms of passes through the data. More recently, Zhong and Kwok (2013) consider an alternating direction method of multipliers (ADMM) variant of this strategy that may be beneficial for problems with more complicated structures.
Mixed Optimization: Another interesting related framework that has been considered since the first version of this work was published is mixed optimization Mahdavi and Jin (2013). Unlike FG methods that evaluate on each iteration or SG methods that evaluate an individual on each iteration, in mixed optimization both of these operations are considered. In this framework, Mahdavi and Jin (2013) show that an convergence rate can be achieved with SG iterations and only FG evaluations. This algorithm also requires a memory and we have found it to be slower than SAG in practice, but an advantage of this method over SAG is that instead of storing all of the variables it only requires storing a single previous parameter value (along with the corresponding value ), since all gradients in the memory can be re-computed given this .
In our analysis we assume that each function in (1) is convex and differentiable, and that each gradient is Lipschitz-continuous with constant , meaning that for all and in and each we have
This is a fairly weak assumption on the functions, and in cases where the
are twice-differentiable it is equivalent to saying that the eigenvalues of the Hessians of eachare bounded above by . We will also assume the existence of at least one minimizer that achieves the optimal function value. We denote the average iterate by
, and the variance of the gradient norms at the optimumby . Our convergence results consider two different initializations for the variables: setting for all , or setting them to the centered gradient at the initial point given by
. We note that all our convergence results are expressed in terms of expectations with respect to the internal randomization of the algorithm (the selection of the random variables), and not with respect to the data which is assumed to be deterministic and fixed.
In addition to this basic convex case discussed above, we will also consider the case where the average function is strongly-convex with constant , meaning that the function is convex. For twice-differentiable , this is equivalent to requiring that the eigenvalues of the Hessian of are bounded below by . This is a stronger assumption that is often not satisfied in practical applications. Nevertheless, in many applications we are free to choose a regularizer of the parameters, and thus we can add an -regularization term as in (2) to transform any convex problem into a strongly-convex problem (in this case we have ). Note that strong-convexity implies the existence of a unique that achieves the optimal function value.
Under these standard assumptions, we now state our convergence result.
With a constant step size of , the SAG iterations satisfy for :
where if we initialize with we have
and if we initialize with we have
Further, if is -strongly convex we have
The proof is given in Appendix B, and involves finding a Lyapounov function for a non-linear stochastic dynamical system defined on the and variables that converges to zero at the above rates, and showing that this function dominates the expected sub-optimality . Note that while the first part is stated for the average , with a trivial change to the proof technique it can be shown to also hold for any iterate where is lower than the average function value up to iteration , . Thus, in addition to the result also holds for the best iterate. We also note that our bounds are valid for any greater than or equal to the minimum satisfying (8), implying an and linear convergence rate for any (but the bound becomes worse as grows). Although initializing each with the centered gradient may have an additional cost and slightly worsens the dependency on the initial sub-optimality , it removes the dependency on the variance of the gradients at the optimum. While we have stated Theorem 1 in terms of the function values, in the strongly-convex case we also obtain a convergence rate on the iterates because we have
Theorem 1 shows that the SAG iterations are advantageous over SG methods in later iterations because they obtain a faster convergence rate. However, the SAG iterations have a worse constant factor because of the dependence on . We can improve the dependence on using an appropriate choice of . In particular, following Le Roux et al. (2012) we can set to the result of iterations of an appropriate SG method. In this setting, the expectation of is in the convex setting, while both and would be in in the strongly-convex setting. If we use this initialization of and set , then in terms of and the SAG convergence rates take the form and in the convex and strongly-convex settings, instead of the and rates implied by Theorem 1. However, in our experiments we do not use an SG initialization but rather use a minor variant of SAG in the early iterations (discussed in the next section), which appears more difficult to analyze but which gives better empirical performance.
An interesting consequence of using a step-size of is that it makes the method adaptive to the strong-convexity constant . That is, for problems with a higher degree of local strong-convexity around the solution , the algorithm will automatically take advantage of this and yield a faster local rate. This can even lead to a local linear convergence rate if the problem is strongly-convex near the optimum but not globally strongly-convex. This adaptivity to the problem difficulty is in contrast to SG methods whose sequence of step sizes typically depend on global constants and thus do not adapt to local strong-convexity.
We have observed in practice that the IAG method with a step size of may diverge. While the step-size needed for convergence of the IAG iterations is not precisely characterized, we have observed that it requires a step-size of approximately in order to converge. Thus, the SAG iterations can tolerate a step size that is roughly times larger, which leads to increased robustness to the selection of the step size. Further, as our analysis and experiments indicate, the ability to use a large step size leads to improved performance of the SAG iterations. Note that using randomized selection with a larger step-size leading to vastly improved performance is not an unprecedented phenomenon; the analysis of Nedic and Bertsekas (2000) shows that the iterations of the basic stochastic gradient method with a constant step-size can achieve the same error bound as full cycles through the data of the cyclic variant of the method by using steps that are times larger (see the discussion after Proposition 3.4). Related results also appear in Collins et al. (2008); Lacoste-Julien et al. (2013) showing the advantage of stochastic optimization strategies over deterministic optimization strategies in the context of certain dual optimization problems.
The convergence rate of the SAG iterations in the strongly-convex case takes a somewhat surprising form. For ill-conditioned problems where , does not appear in the convergence rate and the SAG algorithm has nearly the same convergence rate as the FG method with a step size of , even though it uses iterations which are times cheaper. This indicates that the basic gradient method (under a slightly sub-optimal step-size) is not slowed down by using out-of-date gradient measurements for ill-conditioned problems. Although appears in the convergence rate in the well-conditioned setting where , if we perform iterations of SAG (i.e., one effective pass through the data), the error is multiplied by , which is independent of . Thus, in this setting each pass through the data reduces the excess objective by a constant multiplicative factor that is independent of the problem.
It is interesting to compare the convergence rate of SAG in the strongly-convex case with the known convergence rates for first-order methods (see Nesterov, 2004, §2). In Table 1, we use two examples to compare the convergence rate of SAG to the convergence rates of the standard FG method, the faster AFG method, and the lower-bound for any first-order strategy (under certain dimensionality assumptions) for optimizing a function satisfying our assumptions. In this table, we compare the rate obtained for these FG methods to the rate obtained by running iterations of SAG, since this requires the same number of evaluations of . We also include the rate obtained by running iterations of the recently-introduced ‘minimization by incremental surrogate optimization’ (MISO) method of Mairal (2013), another general method that achieves a linear convergence rate with an iteration cost that is independent of . Example 1 in this table focuses on a well-conditioned case where the rate of SAG is , while Example 2 focuses on an ill-conditioned example where the rate is . Note that in the latter case the rate for the method may be faster.
In Table 1 we see that performing iterations of SAG can actually lead to a rate that is faster than the lower bound for FG methods. Thus, for certain problems SAG can be substantially faster than any FG method that does not use the structure of the problem. However, we note that the comparison is somewhat problematic because in the SAG (and MISO) rates is the Lipschitz constant of the functions, while in the FG method we only require that is an upper bound on the Lipschitz continuity of so it may be much smaller. To give a concrete example that takes this into account and also considers the rates of dual methods and coordinate-wise methods, in Appendix A we attempt to more carefully compare the rates obtained for SAG with the rates of primal and dual FG and coordinate-wise methods for the special case of -regularized least-squares regression.
|Algorithm||Step Size||Theoretical Rate||Rate in Example 1||Rate Example 2|
|MISO ( iters)|
|SAG ( iters)|
After the first version of this work was released, a convergence rate with a similar form was shown for a restricted class of problems by Shalev-Schwartz and Zhang (2013b). Their work considered a stochastic dual coordinate ascent method with exact coordinate optimization for -regularized linearly-parameterized models (a problem class that is discussed further in the next section). For this class of problems, they show that a sub-optimality of can be achieved in iterations (where is the Lipschitz constant of the loss functions and is a positive parameter of the strongly-convex regularizer). If we apply the SAG algorithm to their special case, Theorem 1 implies that the SAG algorithm requires iterations, and we note that . Shalev-Schwartz and Zhang (2013b) require that be positive and specified in advance, so unlike SAG their method is not adaptive to the strong convexity of the problem.
In Algorithm 1 we give pseudo-code for an implementation of the basic method, where we use a variable to track the quantity . This section focuses on further implementation details that are useful in practice. In particular, we discuss modifications that lead to better practical performance than the basic Algorithm 1, including ways to reduce the storage cost, how to handle regularization, how to set the step size, using mini-batches, and using non-uniform sampling. Note that an implementation of the algorithm that incorporates many of these aspects is available from the first author’s webpage.
For many problems the storage cost of for the vectors is prohibitive, but we can often use the structure of the gradients to reduce this cost. For example, a commonly-used specialization of (1) is linearly-parameterized models which take form
Since each is constant, for these problems we only need to store the scalar for rather than the full gradient . This reduces the storage cost from down to .
For problems where the vectors are sparse, an individual gradient will inherit the sparsity pattern of the corresponding . However, the update of in Algorithm 1 appears unappealing since in general will be dense, resulting in an iteration cost of . Nevertheless, we can take advantage of the simple form of the SAG updates to implement a ‘just-in-time’ variant of the SAG algorithm where the iteration cost is proportional to the number of non-zeroes in . In particular, we do this by not explicitly storing the full vector after each iteration. Instead, on each iteration we only compute the elements corresponding to non-zero elements of , by applying the sequence of updates to each variable since the last iteration where it was non-zero in . This sequence of updates can be applied efficiently since it simply involves changing the step size. For example, if variable has been zero in for iterations, then we can compute the needed value using
This update allows SAG to be efficiently applied to sparse data sets where and are both in the millions or higher but the number of non-zeros is much less than .
In the update of in Algorithm 1, we normalize the direction by the total number of data points . When initializing with we believe this leads to steps that are too small on early iterations of the algorithm where we have only seen a fraction of the data points, because many variables contributing to are set to the uninformative zero-vector. Following Blatt et al. (2007), the more logical normalization is to divide by , the number of data points that we have seen at least once (which converges to once we have seen the entire data set), leading to the update . Although this modified SAG method appears more difficult to analyze, in our experiments we found that running the basic SAG algorithm from the very beginning with this modification outperformed the basic SAG algorithm as well as the SG/SAG hybrid algorithm mentioned in the Section 3. In addition to using the gradient information collected during the first iterations, this modified SAG algorithm is also advantageous over hybrid SG/SAG algorithms because it only requires estimating a single constant step size.
In the case of regularized objectives like (2), the cost of computing the gradient of the regularizer is independent of . Thus, we can use the exact gradient of the regularizer in the update of , and only use to approximate the sum of the functions. By incorporating the gradient of the regularizer explicitly, the update for in Algorithm 1 becomes , and in the case of -regularization the update for becomes
If the loss function gradients are sparse as in Section 4.1, then these modifications lead to a reduced storage requirement even though the gradient of the regularizer is dense. Further, although the update of again appears to require dense vector operations, we can implement the algorithm efficiently if the are sparse. In particular, to allow efficient multiplication of by the scalar , it is useful to represent in the form , where is a scalar and is a vector (as done in Shalev-Shwartz et al., 2011). Under this representation, we can multiply by a scalar in by simply updating (though to prevent becoming too large or too small we may need to occasionally re-normalize by setting and ). To efficiently implement the vector subtraction operation, we can use a variant of the just-in-time updates from Section 4.1. In Algorithm 2, we give pseudo-code for a variant of SAG that includes all of these modifications, and thus uses no full-vector operations. This code uses a vector to keep track of the scalars , a vector to determine whether a data point has previously been visited, a vector to track the last time a variable was updated, and a vector to keep track of the cumulative sums needed to implement the just-in-time updates.
In many scenarios we may need to solve a set of closely-related optimization problems. For example, we may want to apply Algorithm 2 to a regularized objective of the form (2) for several values of the regularization parameter . Rather than solving these problems independently, we might expect to obtain better performance by warm-starting the algorithm. Although initializing with the solution of a related problem can improve performance, we can expect an even larger performance improvement if we also use the gradient information collected from a run of SAG for a close value of . For example, in Algorithm 2 we could initialize , , , , and based on a previous run of the SAG algorithm. In this scenario, Theorem 1 suggests that it may be beneficial in this setting to center the variables around .
In our experiments we have observed that utilizing a step size of , as in standard FG methods, always converged and often performed better than the step size of suggested by our analysis. Thus, in our experiments we used even though we do not have a formal analysis of the method under this step size. We also found that a step size of , which in the strongly-convex case corresponds to the best fixed step size for the FG method in the special case of (Nesterov, 2004, Theorem 2.1.15), sometimes yields even better performance (though in other cases it performs poorly).
In general the Lipschitz constant will not be known, but we may obtain a reasonable approximation of a valid by evaluating values while running the algorithm. In our experiments, we used a basic line-search where we start with an initial estimate , and double this estimate whenever we do not satisfy the inequality
which must be true if is valid. An important property of this test is that it depends on but not on , and thus the cost of performing this test is independent of . To avoid instability caused by comparing very small numbers, we only do this test when . Since is a global quantity but the algorithm will eventually remain within a neighbourhood of the optimal solution, it is possible that a smaller estimate of (and thus a larger step size) can be used as we approach . To potentially take advantage of this, we initialize with the slightly smaller at each iteration, so that the estimate of is halved if we do iterations (an effective pass through the data) and never violate the inequality. Note that in the case of -regularized objectives, we can perform the line-search to find an estimate of the Lipschitz constant of rather than , and then simply add to this value to take into account the effect of the regularizer.
Because of the use of vectorization and parallelism in modern architectures, practical SG implementations often group functions into ‘mini-batches’ and perform SG iterations on the mini-batches. We can also use mini-batches within the SAG iterations to take advantage of the same vectorization and parallelism. Additionally, for problems with dense gradients mini-batches can dramatically decrease the storage requirements of the algorithm, since we only need to store a vector for each mini-batch rather than for each example. Thus, for example, using a mini-batch of size leads to a -fold reduction in the storage cost.
A subtle issue that arises when using mini-batches is that the value of in the Lipschitz condition (8) is based on the mini-batches instead of the original functions . For example, consider the case where we have a batch and the minimum value of in (8) for each is given by . In this case, a valid value of for the function would be . We refer to this as . But we could also consider using . The value is still valid and will be smaller than unless all are equal. We could even consider the minimum possible value of , which we refer to as because (if each is twice-differentiable) it is equal to the maximum eigenvalue of across all . Note that , although will typically be more difficult to compute than or (although a line-search as discussed in the previous section can reduce this cost). Due to the potential of using a smaller , we may obtain a faster convergence rate by using larger mini-batches. However, in terms of passes through the data this faster convergence may be offset by the higher iteration cost associated with using mini-batches.
In standard SG methods, it is crucial to sample the functions uniformly, at least asymptotically, in order to yield an unbiased gradient estimate and subsequently achieve convergence to the optimal value (alternately, the bias induced by non-uniform sampling would need to be asymptotically corrected). In SAG iterations, however, the weight of each gradient is constant and equal to , regardless of the frequency at which the corresponding function is sampled. We might thus consider sampling the functions non-uniformly, without needing to correct for this bias. Though we do not yet have any theoretical proof as to why a non-uniform sampling might be beneficial, intuitively we would expect that we do not need to sample functions whose gradient changes slowly as often as functions whose gradient changes more quickly. Indeed, we provide here an argument to justify a non-uniform sampling strategy based on the Lipschitz constants of the individual gradients .
Let again be the Lipschitz constant of , and assume that the functions are placed in increasing order of Lipschitz constants, so that . In the ill-conditioned setting where the convergence rate depends on , a simple way to improve the rate by decreasing is to replace by two functions and such that
We have thus replaced the original problem by a new, equivalent problem where:
has been replaced by ,
for is ,
and are equal to .
Hence, if , this problem has the same but a smaller than the original one, improving the bound on the convergence rate. By duplicating
, we increase its probability of being sampled fromto , but we also replace by a noisier version, i.e. . Using a noisier version of the gradient appears detrimental, so we assume that the improvement comes from increasing the frequency at which is sampled, and that logically we might obtain a better rate by simply sampling more often in the original problem and not explicitly duplicating the data.
We now consider the extreme case of duplicating each function a number of times equal to the Lipschitz constant of their gradient, assuming that these constants are integers. The new problem becomes
The function is now written as the sum of functions, each with a gradient with Lipschitz constant . The new problem has the same as before, but now has an equal to the average of the Lipschitz constants across the , rather than their maximum, thus improving the bound on the convergence rate. Sampling these functions uniformly is now equivalent to sampling the original ’s according to their Lipschitz constant. Thus, we might expect to obtain better performance by, instead of creating a larger problem by duplicating the functions in proportion to their Lipschitz constant, simply sampling the functions from the original problem in proportion to their Lipschitz constants.
Sampling in proportion to the Lipschitz constants of the gradients was explored by Nesterov (2010) in the context of coordinate descent methods, and is also somewhat related to the sampling scheme used by Strohmer and Vershynin (2009) in the context of their randomized Kaczmarz algorithm. Such a sampling scheme makes the iteration cost depend on , due to the need to generate samples from a general discrete distribution over variables. However, after an initial preprocessing cost of we can sample from such distributions in using a simple binary search (see Robert and Casella, 2004, Example 2.10).
Unfortunately, sampling the functions according to the Lipschitz constants and using a step size of does not seem to converge in general. However, by changing the number of times we duplicate each
, we can interpolate between the Lipschitz sampling and the uniform sampling. In particular, if each functionis duplicated times, where is the Lipschitz constant of and a positive number, then the new problem becomes
Unlike in the previous case, these functions have gradients with different Lipschitz constants. Denoting , the maximum Lipschitz constant is equal to and we must thus use the step size .
In this section we perform empirical evaluations of the SAG iterations. We first compare the convergence of an implementation of the SAG iterations to a variety of competing methods available. We then seek to evaluate the effect of different algorithmic choices such as the step size, mini-batches, and non-uniform sampling.
The theoretical convergence rates suggest the following strategies for deciding on whether to use an FG or an SG method:
If we can only afford one pass through the data, then an SG method should be used.
If we can afford to do many passes through the data (say, several hundred), then an FG method should be used.
We expect that the SAG iterations will be most useful between these two extremes, where we can afford to do more than one pass through the data but cannot afford to do enough passes to warrant using FG algorithms like the AFG or L-BFGS methods. To test whether this is indeed the case in practice, we perform a variety of experiments evaluating the performance of the SAG algorithm in this scenario.
Although the SAG algorithm can be applied more generally, in our experiments we focus on the important and widely-used -regularized logistic regression problem
as a canonical problem satisfying our assumptions. In our experiments we set the regularization parameter to , which is in the range of the smallest values that would typically be used in practice, and thus which results in the most ill-conditioned problems of this form that would be encountered. Our experiments focus on the freely-available benchmark binary classification data sets listed in Table 2. The quantum and protein data set was obtained from the KDD Cup 2004 website;111http://osmot.cs.cornell.edu/kddcup the covertype, rcv1, news, and rcv1Full data sets were obtained from the LIBSVM Data website; 222http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets; the sido data set was obtained from the Causality Workbench website,333 http://www.causality.inf.ethz.ch/home.php the spam data set was prepared by Carbonetto (2009, §2.6.5) using the TREC 2005 corpus444http://plg.uwaterloo.ca/~gvcormac/treccorpus; and the alpha data set was obtained from the Pascal Large Scale Learning Challenge website555http://largescale.ml.tu-berlin.de. We added a (regularized) bias term to all data sets, and for dense features we standardized so that they would have a mean of zero and a variance of one. To obtain results that are independent of the practical implementation of the algorithm, we measure the objective as a function of the number of effective passes through the data, measured as the number of times we evaluate divided by the number of examples . If they are implemented to take advantage of the sparsity present in the data sets, the runtimes of all algorithms examined in this section differ by at most a constant times this measure.
|Data set||Data Points||Variables||Reference|
|quantum||50 000||78||(Caruana et al., 2004)|
|protein||145 751||74||(Caruana et al., 2004)|
|covertype||581 012||54||Blackard, Jock, and Dean (Frank and Asuncion, 2010)|
|rcv1||20 242||47 236||(Lewis et al., 2004)|
|news||19 996||1 355 191||(Keerthi and DeCoste, 2005)|
|spam||92 189||823 470||(Cormack and Lynam, 2005; Carbonetto, 2009)|
|rcv1Full||697 641||47 236||(Lewis et al., 2004)|
|sido||12 678||4 932||(Guyon, 2008)|
In our first experiment we compared the following variety of competitive FG and SG methods:
AFG: A variant of the accelerated full gradient method of Nesterov (1983), where iterations of (3) with a step size of are interleaved with an extrapolation step. We used an adaptive line-search to estimate a local based on the variant proposed for -regularized logistic regression by Liu et al. (2009).
L-BFGS: A publicly-available limited-memory quasi-Newton method that has been tuned for log-linear models such as logistic regression.666http://www.di.ens.fr/~mschmidt/Software/minFunc.html This method is the most complicated method we considered.
SG: The stochastic gradient method described by iteration (4). Since setting the step-size in this method is a tenuous issue, we chose the constant step size that gave the best performance (in hindsight) among all powers of (we found that this constant step-size strategies gave better performance than the variety of decreasing step-size strategies that we experimented with).
ASG: The average of the iterations generated by the SG method above, where again we choose the best step size among all powers of .777Note that we also compared to a variety of other SG methods including the popular Pegasos SG method of Shalev-Shwartz et al. (2011), SG with momentum, SG with gradient averaging, the regularized dual averaging method of Xiao (2010) (a stochastic variant of the primal-dual subgradient method of Nesterov (2009) for regularized objectives), the accelerated SG method of Delyon and Juditsky (1993), SG methods that only average the later iterations as in the ‘optimal’ SG method for non-smooth optimization of Rakhlin et al. (2012) , the epoch SG method of
, the epoch SG method ofHazan and Kale (2011), the ‘nearly-optimal’ SG method of Ghadimi and Lan (2010), a diagonally-scaled SG method using the inverse of the coordinate-wise Lipshitz constants as the diagonal terms, and the adaptive diagonally-scaled AdaGrad method of Duchi et al. (2011). However, we omit results obtained using these algorithms since they never performed substantially better than the minimum between the SG and ASG methods when their step-size was chosen in hindsight.
IAG: The incremental aggregated gradient method of Blatt et al. (2007) described by iteration (5) with a cyclic choice of . We used the re-weighting described in Section 4.2, we used the exact regularization as described in Section 4.3, and we chose the step-size that gave the best performance among all powers of .
SAG-LS: The proposed stochastic average gradient method described by iteration (5). We used the re-weighting described in Section 4.2, the exact regularization as described in Section 4.3, and we used a step size of where is an approximation of the Lipschitz constant for the negative log-likelihoods . Although this Lipschitz constant is given by , we used the line-search described in Section 4.6 to estimate it, to test the ability to use SAG as a black-box algorithm (in addition to avoiding this calculation and potentially obtaining a faster convergence rate by obtaining an estimate that could be smaller than the global value). To initialize the line-search we set .
We plot the results of the different methods for the first effective passes through the data in Figure 1. We can observe several trends across the experiments:
FG vs. SG: Although the performance of SG methods is known to be catastrophic if the step size is not chosen carefully, after giving the SG methods (SG and ASG) an unfair advantage (by allowing them to choose the best step-size in hindsight), the SG methods always do substantially better than the FG methods (AFG and L-BFGS) on the first few passes through the data. However, the SG methods typically make little progress after the first few passes. In contrast, the FG methods make steady progress and eventually the faster FG method (L-BFGS) typically passes the SG methods.
(FG and SG) vs. SAG: The SAG iterations seem to achieve the best of both worlds. They start out substantially better than FG methods, often obtaining similar performance to an SG method with the best step-size chosen in hindsight. But the SAG iterations continue to make steady progress even after the first few passes through the data. This leads to better performance than SG methods on later iterations, and on most data sets the sophisticated FG methods do not catch up to the SAG method even after passes through the data.
IAG vs. SAG: Even though these two algorithms differ in only the seemingly-minor detail of selecting data points at random (SAG) compared to cycling through the data (IAG), the performance improvement of SAG over its deterministic counterpart IAG is striking (even though the IAG method was allowed to choose the best step-size in hindsight). We believe this is due to the larger step sizes allowed by the SAG iterations, which would cause the IAG iterations to diverge.
For the -regularized logistic regression problem, an alternative means to take advantage of the structure of the problem and achieve a linear convergence rate with a cheaper iteration cost than FG methods is to use randomized coordinate optimization methods. In particular, we can achieve a linear convergence rate by applying coordinate descent to the primal (Nesterov, 2010) or coordinate-ascent to the dual (Shalev-Schwartz and Zhang, 2013b). In our second experiment, we included the following additional methods in this comparison:
PCD: The randomized primal coordinate-descent method of Nesterov (2010), using a step-size of , where is the Lipschitz-constant with respect to coordinate of . Here, we sampled the coordinates uniformly.
PCD-L: The same as above, but sampling coordinates according to their Lipschitz constant, which can lead to an improved convergence rate (Nesterov, 2010).
DCA: Applying randomized coordinate ascent to the dual, with uniform example selection and an exact line-search (Shalev-Schwartz and Zhang, 2013b).
As with comparing FG and SG methods, it is difficult to compare coordinate-wise methods to FG and SG methods in an implementation-independent way. To do this in a way that we believe is fair (when discussing convergence rates), we measure the number of effective passes of the DCA method as the number of iterations of the method divided by (since each iteration accesses a single example as in SG and SAG iterations). We measure the number of effective passes for the PCD and PCD-L methods as the number of iterations multiplied by so that effective pass for this method has a cost of as in FG and SG methods. We ignore the additional cost associated with the Lipschitz sampling for the PCD-L method (as well as the expense incurred because the PCD-L method tended to favour updating the bias variable for sparse data sets) and we also ignore the cost of numerically computing the optimal step-size for the DCA method.
We compare the performance of the randomized coordinate optimization methods to several of the best methods from the previous experiment in Figure 2. In these experiments we observe the following trends:
PCD vs. PCD-L: For the problems with (top and bottom rows of Figure 2), there is little difference between uniform and Lipschitz sampling of the coordinates. For the problems with (middle row of Figure 2), sampling according to the Lipschitz constant leads to a large performance improvement over uniform sampling.
PCD vs. DCA: For the problems with , DCA and PCD-L have similar performance. For the problems with , the performance of the methods typically differed but neither strategy tended to dominate the other.
(PCD and DCA) vs. (SAG): For some problems, the PCD and DCA methods have performance that is similar to SAG-LS and the DCA method even gives better performance than SAG-LS on one data set. However, for many data sets either the PCD-L or the DCA method (or both) perform poorly while the SAG-LS iterations are among the best or substantially better than all other methods on every data set. This suggests that, while coordinate optimization methods are clearly extremely effective for some problems, the SAG method tends to be a more robust choice across problems.
In our prior work we analyzed the step-sizes and (Le Roux et al., 2012), while Section 3 considers the choice and Section 4.5 discusses the choices and as in FG methods. In Figure 3 we compare the performance of these various strategies to the performance of the SAG algorithm with our proposed line-search as well as the IAG and SAG algorithms when the best step-size is chosen in hindsight. In this plot we see the following trends:
Proposition 1 of Le Roux et al. (2012): Using a step-size of performs poorly, and makes little progress compared to the other methods. This makes sense because Proposition 1 in Le Roux et al. (2012) implies that the convergence rate (in terms of effective passes through the data) under this step size will be similar to the basic gradient method, which is known to perform poorly unless the problem is very well-conditioned.
Proposition 2 of Le Roux et al. (2012): Using a step-size of performs extremely well on the data sets with (middle row). In contrast, for the data sets with it often performs very poorly, and in some cases appears to diverge. This is consistent with Proposition 2 in Le Roux et al. (2012), which shows a fast convergence rate under this step size only if certain conditions on hold.
Theorem 1: Using a step-size of performs consistently better than the smaller step size , but in some cases it performs worse than . However, in contrast to , the step size always has reasonable performance.
Section 4.5: The step size of performs performs extremely well on the data sets with , and performs better than the step sizes discussed above on all but one of the remaining data sets. The step size of seems to perform the same or slightly better than using except on one data set where it performs poorly.
Line-Search: Using the line-search from Section 4.6 tends to perform as well or better than the various constant step size strategies, and tends to have similar performance to choosing the best step size in hindsight.
IAG vs. SAG: When choosing the best step size in hindsight, the SAG iterations tend to choose a much larger step size than the IAG iterations. The step sizes chosen for SAG were to times larger than the step sizes chosen by IAG, and always lead to better performance by several orders of magnitude.
As we discuss in Section 4.7, when using mini-batches within the SAG iterations there is a trade-off between the higher iteration cost of using mini-batches and the faster convergence rate obtained using mini-batches due to the possibility of using a smaller value of . In Figure 4, we compare (on the dense data sets) the excess sub-optimality as a function of the number of examples seen for various mini-batch sizes and the three step-size strategies , , and discussed in Section 4.7.
Several conclusions may be drawn from these experiments:
Even though Theorem 1 hints at a maximum mini-batch size of without loss of convergence speed, this is a very conservative estimate. In our experiments, the original value of was on the order of and mini-batch sizes of up to 500 could be used without a loss in performance. Not only does this yield large memory storage gains, it would increase the computational efficiency of the algorithm when taking into account vectorization.
To achieve fast convergence, it is essential to use a larger step-size when larger mini-batches are used. For instance, in the case of the quantum dataset with a mini-batch size of , we have , and . In the case of the covertype dataset and for a mini-batch size of , we have , and .
In our final experiment, we explored the effect of using the non-uniform sampling strategy discussed in Section 4.8. In Figure 5, we compare several of the SAG variants with uniform sampling to the following two methods based on non-uniform sampling:
SAG (Lipschitz): In this method we sample the functions in proporition to , where is the global Lipschitz constant of the corresponding and we set to the average of these constants, . Plugging in these values into the formula at the end of Section 4.8 and using to denote the maximum value, we set the step-size to .
SAG-LS (Lipschitz): In this method we formed estimates of the quantities , , and . The estimator is computed in the same way as the SAG-LS method. To estimate each , we keep track of an estimate for each and we set to the average of the values among the that we have sampled at least once. We set if example was not selected and otherwise we initialize to and perform the line-search until we have a valid (this means that will be approximately halved if we perform a full pass through the data set and never violate the inequality). To ensure that we do not ignore important unseen data points for too long, in this method we sample a previously unseen function with probability , and otherwise we sample from the previously seen in proportion to . To prevent relying too much on our initially-poor estimate of , we use a step size of