Mark Schmidt
email: schmidtm@cs.ubc.ca

Nicolas Le Roux
email: nicolas@le-roux.name

Francis Bach
email: francis.bach@ens.fr

Minimizing Finite Sums with the Stochastic Average Gradient

, which only lead to a faster rate for well-conditioned strongly-convex problems. 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.

Introduction

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, minimize

x∈R p 1 n n i=1 (a T i x -b i ) 2 ,
where the a i ∈ R p and b i ∈ R are the data samples associated with a regression problem. Another important example is logistic regression, minimize

x∈R p 1 n n i=1 log(1 + exp(-b i a i x)),
where the a i ∈ R p and b i ∈ {-1, 1} are the data samples associated with a binary classification problem. A key challenge arising in modern applications is that the number of data points n (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 n is very large are stochastic gradient (SG) methods [START_REF] Robbins | A stochastic approximation method[END_REF]Monro, 1951, Bottou and[START_REF] Bottou | Large scale online learning[END_REF]. 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, minimize

x∈R p g(x) := 1 n n i=1 f i (x). (1)
In this work, we focus on such finite data problems where each f i is smooth and convex.

In addition to this basic setting, we will also be interested in cases where the sum g 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 2 -norm, resulting in problems of the form minimize

x∈R p λ 2 x 2 + 1 n n i=1 l i (x), (2)
where each l i 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

f i (x) := λ 2 x 2 + l i (x).
The resulting function g will be strongly-convex provided that the individual loss functions l i are convex. An extensive list of convex loss functions used in a statistical data-fitting context is given by [START_REF] Teo | A scalable modular convex solver for regularized risk minimization[END_REF], and non-smooth loss functions (or regularizers) can also be put in this framework by using smooth approximations (for example, see [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]).

For optimizing problem (1), the standard deterministic or full gradient (FG) method, which dates back to [START_REF] Cauchy | Méthode générale pour la résolution des systèmes d'équations simultanées[END_REF], uses iterations of the form

x k+1 = x k -α k g (x k ) = x k - α k n n i=1 f i (x k ), (3)
where α k is the step size on iteration k. Assuming that a minimizer x * exists, then under standard assumptions the sub-optimality achieved on iteration k of the FG method with a constant step size is given by g(x k ) -g(x * ) = O(1/k), when g is convex [see Nesterov, 2004, Corollary 2.1.2]. This results in a sublinear convergence rate. When g is strongly-convex, the error also satisfies

g(x k ) -g(x * ) = O(ρ k ),
for some ρ < 1 which depends on the condition number of g [see Nesterov, 2004, Theorem 2.1.5]. 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 n is large because its iteration cost scales linearly in n.

The main appeal of SG methods is that they have an iteration cost which is independent of n, making them suited for modern problems where n may be very large. The basic SG method for optimizing (1) uses iterations of the form

x k+1 = x k -α k f i k (x k ), (4)
where at each iteration an index i k is sampled uniformly from the set {1, . . . , n}. The randomly chosen gradient f i k (x k ) yields an unbiased estimate of the true gradient g (x k ) and one can show under standard assumptions (see [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF]) that, for a suitably chosen decreasing step-size sequence {α k }, the SG iterations have an expected sub-optimality for convex objectives of

E[g(x k )] -g(x * ) = O(1/ √ k),
and an expected sub-optimality for strongly-convex objectives of

E[g(x k )] -g(x * ) = O(1/k).
In these rates, the expectations are taken with respect to the selection of the i k 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 [START_REF] Nemirovski | Problem complexity and method efficiency in optimization[END_REF], [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF], [START_REF] Agarwal | Information-theoretic lower bounds on the oracle complexity of stochastic convex optimization[END_REF]). 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 [START_REF] Blatt | A convergent incremental gradient method with a constant step size[END_REF]. 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

x k+1 = x k - α k n n i=1 y k i , (5)
where at each iteration a random index i k is selected and we set

y k i = f i (x k ) if i = i k , y k-1 i otherwise. ( 6
)
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 n. Despite the low cost of the SAG iterations, we show in this paper that with a constant step-size the SAG iterations have an O(1/k) convergence rate for convex objectives and a linear convergence rate for strongly-convex objectives, like the FG method. That is, by having access to i k and by keeping a memory of the most recent gradient value computed for each index i, 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.

One of the main contexts where minimizing the sum of smooth convex functions arises is machine learning. In this context, g is often an empirical risk (or a regularized empirical risk), which is a sample average approximation to the true risk that we are interested in. It is known that with n training examples the empirical risk minimizer (ERM) has an error for the true risk of O(1/ √ n) in the convex case and O(1/n) in the strongly-convex case. Since these rates are achieved by doing one pass through the data with an SG method, in the worst case the SAG algorithm applied to the empirical risk cannot improve the convergence rate in terms of the true risk over this simple method. Nevertheless, [START_REF] Srebro | Theoretical basis for "more data less work[END_REF] note that "overwhelming empirical evidence shows that for almost all actual data, the ERM is better. However, we have no understanding of why this happens". Although our analysis does not give insight into the better performance of ERM, our analysis shows that the SAG algorithm will be preferable to SG methods for finding the ERM and hence for many machine learning applications.

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 60 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 O(np) to O(n) in the common scenario where each f i only depends on a linear function of x, 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 O(1/k) 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 n 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 n times as fast as the standard gradient method. To prove this stronger result, Theorem 1 employs a Lyapunov function that generalizes the Lyapunov functions used in Propositions 1 and 2 of the previous work. This new Lyapunov function leads to a unified proof for both the convex and the strongly-convex cases, and for both well-conditioned and ill-conditioned problems. However, this more general Lyapunov function leads to a more complicated analysis. To significantly simplify the formal proof, we use a computed-aided strategy to verify the non-negativity of certain polynomials that arise in the proof. 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).

Related Work

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

x k+1 = x k -α k f i k (x k ) + β k (x k -x k-1 ),
see [START_REF] Tseng | An incremental gradient(-projection) method with momentum term and adaptive stepsize rule[END_REF]. It is common to set all β k = β for some constant β, and in this case we can rewrite the SG with momentum method as

x k+1 = x k - k j=1 α j β k-j f ij (x j ).
We can re-write the SAG updates (5) in a similar form as

x k+1 = x k - k j=1 α k S(j, i 1:k )f ij (x j ), (7)
where the selection function S(j, i 1:k ) is equal to 1/n if j is the maximum iteration number where example i j was selected and is set to 0 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,

x k+1 = x k -α k k k j=1 f ij (x j
), 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 [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] 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 [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF], 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 x k 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 [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF]Juditsky, 1992, Bach and[START_REF] Bach | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF]]. Bather'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 O(1/k) rate for strongly-convex optimization even for non-smooth objectives [Hazan andKale, 2011, Rakhlin et al., 2012]. However, none of these methods improve on the O(1/ √ k) and O(1/k) 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 [START_REF] Nesterov | A method for unconstrained convex minimization problem with the rate of convergence O(1/k 2 )[END_REF], as well as classical techniques based on quadratic approximations such as diagonallyscaled FG methods, non-linear conjugate gradient, quasi-Newton, and Hessian-free Newton methods (see [START_REF] Nocedal | Numerical Optimization[END_REF]). There has been a substantial amount of work on developing stochastic variants of these algorithms, with several of the notable recent examples including [START_REF] Bordes | Sgd-qn: Careful quasi-newton stochastic gradient descent[END_REF], [START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF], [START_REF] Sunehag | Variable metric stochastic approximation theory[END_REF], [START_REF] Ghadimi | Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization[END_REF], [START_REF] Martens | Deep learning via Hessian-free optimization[END_REF], [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF]. [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF] 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 [START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF][START_REF] Ghadimi | Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization[END_REF][START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF]. However, we are not aware of any existing method of this flavor that improves on the O(1/ √ k) and O(1/k) 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 f i .

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 k and converges linearly to 0. The second part is independent of k and does not converge to 0. 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 [START_REF] Bach | Non-strongly-convex smooth stochastic approximation with convergence rate o(1/n)[END_REF], 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) [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF], or under extremely strong assumptions about the relationship between the functions f i [START_REF] Solodov | Incremental gradient algorithms with stepsizes bounded away from zero[END_REF]. 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[START_REF] Delyon | Accelerated stochastic approximation[END_REF]. 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. [START_REF] Bertsekas | A new class of incremental gradient methods for least squares problems[END_REF] 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 f i into 'batches' of increasing size and perform SG iterations on the batches. [START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF] give conditions under which this strategy achieves the O(1/k) and O(ρ k ) convergence rates of FG methods. However, in both cases the iterations that achieve the faster rates have a cost that is not independent of n, as opposed to SAG iterations.

Incremental Aggregated Gradient: [START_REF] Blatt | A convergent incremental gradient method with a constant step size[END_REF] 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 i k rather than sampling the i k 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 O(1/k) 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 certain highly-restricted classes of problems, it is possible to show faster convergence rates for methods that only operate on a single function f i . For example, [START_REF] Strohmer | A randomized kaczmarz algorithm with exponential convergence[END_REF] show that the randomized Kaczmarz method with a particular sampling scheme achieves a linear convergence rate for the problem of solving consistent linear systems. It is also known that the SG method with a constant step-size has the O(1/k) and O(ρ k ) convergence rates of FG methods if f i (x) is bounded by a linear function of g (x) for all i and x [START_REF] Schmidt | Fast convergence of stochastic gradient descent under a strong growth condition[END_REF]. This is the strong condition required by [START_REF] Solodov | Incremental gradient algorithms with stepsizes bounded away from zero[END_REF] to show convergence of the SG method with a constant step size. Unlike these previous works, our analysis in the next section applies to general f i that satisfy standard assumptions, and only requires gradient evaluations of the functions f i rather than dual block-coordinate steps.

Subsequent work:

Since the first version of this work was released, there has been an explosion of research into stochastic gradient methods with faster convergence rates. It has been shown that similar rates can be achieved for certain constrained and non-smooth problems, that similar rates can be achieved without the memory requirements, that Newton-like variants of the method may be possible, and that similar rates can be achieved with other algorithms. In Section 6, we survey these recent developments.

Convergence Analysis

In our analysis we assume that each function f i in (1) is convex and differentiable, and that each gradient f i is Lipschitz-continuous with constant L, meaning that for all x and y in R p and each i we have

f i (x) -f i (y) L x -y . (8)
This is a fairly weak assumption on the f i functions, and in cases where the f i are twice-differentiable it is equivalent to saying that the eigenvalues of the Hessians of each f i are bounded above by L. We will also assume the existence of at least one minimizer x * that achieves the optimal function value. We denote the average iterate by xk = 1 k k-1 i=0 x i , and the variance of the gradient norms at the optimum x * by σ 2 = 1 n i f i (x * ) 2 . Our convergence results consider two different initializations for the y 0 i variables: setting y 0 i = 0 for all i, or setting them to the centered gradient at the initial point x 0 given by y 0 i = f i (x 0 ) -g (x 0 ). 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 i k ), 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 g = 1 n n i=1 f i is strongly-convex with constant µ > 0, meaning that the function x → g(x) -µ 2 x 2 is convex. For twice-differentiable g, this is equivalent to requiring that the eigenvalues of the Hessian of g 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 2 -regularization term as in (2) to transform any convex problem into a strongly-convex problem (in this case we have µ ≥ λ). Note that strongconvexity implies the existence of a unique x * that achieves the optimal function value.

Under these standard assumptions, we now state our convergence result.

Theorem 1. With a constant step size of α k = 1 16L , the SAG iterations satisfy for k ≥ 1:

E[g(x k )] -g(x * ) 32n k C 0 ,
where if we initialize with y 0 i = 0 we have

C 0 = g(x 0 ) -g(x * ) + 4L n x 0 -x * 2 + σ 2 16L ,
and if we initialize with y 0 i = f i (x 0 ) -g (x 0 ) we have

C 0 = 3 2 g(x 0 ) -g(x * ) + 4L n x 0 -x * 2 .
Further, if g is µ-strongly convex we have

E[g(x k )] -g(x * ) 1 -min µ 16L , 1 8n k C 0 .
The proof is given in Appendix B of the extended version of this paper, and involves finding a Lyapunov function for a non-linear stochastic dynamical system defined on the y k i and x k variables that converges to zero at the above rates, and showing that this function dominates the expected sub-optimality [E[g(x k )] -g(x * )]. This is the same approach used to show Proposition 1 and 2 in the conference version of the paper [START_REF] Le Roux | A stochastic gradient method with an exponential convergence rate for strongly-convex optimization with finite training sets[END_REF], but in this work we use a more general Lyapunov function that gives a much faster rate for ill-conditioned problems and also allows us to analyze problems that are not strongly-convex. To simplify the analysis of this more complicated Lyapunov function, our new proof verifies positivity of certain polynomials that arise in the bound using a computer-aided approach.

Note that while the first part of Theorem (1) is stated for the average xk , with a trivial change to the proof technique it can be shown to also hold for any iterate x k where g(x k ) is lower than the average function value up to iteration k,

1 k k-1 i=0 g(x i
). Thus, in addition to xk the result also holds for the best iterate. We also note that our bounds are valid for any L greater than or equal to the minimum L satisfying (8), implying an O(1/k) and linear convergence rate for any α k 1/16L (but the bound becomes worse as L grows). Although initializing each y 0 i with the centered gradient may have an additional cost and slightly worsens the dependency on the initial sub-optimality (g(x 0 ) -g(x * )), it removes the dependency on the variance σ 2 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

µ 2 x k -x * 2 g(x k ) -g(x * ).
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 n. We can improve the dependence on n using an appropriate choice of x 0 . In particular, following Le Roux et al. [2012] we can set x 0 to the result of n iterations of an appropriate SG method. In this setting, the expectation of g(x 0 ) -g(x * ) is O(1/ √ n) in the convex setting, while both g(x 0 ) -g(x * ) and x 0 -x * 2 would be in O(1/n) in the strongly-convex setting. If we use this initialization of x 0 and set y 0 i = f i (x 0 ) -g (x 0 ), then in terms of n and k the SAG convergence rates take the form O( √ n/k) and O(ρ k /n) in the convex and stronglyconvex settings, instead of the O(n/k) and O(ρ k ) 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 α k = 1/16L 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 x * , 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 α k = 1 16L 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 1/nL in order to converge. Thus, the SAG iterations can tolerate a step size that is roughly n 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 [START_REF] Nedic | Convergence rate of incremental subgradient algorithms[END_REF] 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 n times larger (see

Algorithm

Step Size Theoretical Rate Rate in Case 1 Rate Case 2 FG

1 L 1 -µ L 2 0.9998 1.000 FG 2 µ+L 1 -2µ L+µ 2 0.9996 1.000 AFG 1 L 1 -µ L 0.9900 0.9990 Lower-Bound - 1 - 2 √ µ √ L+ √ µ 2 0.9608 0.9960 SAG (n iters) 1 16L 1 -min{ µ 16L , 1 8n } n 0.8825 0.9938
Table 1: Comparison of convergence rates of first-order methods to the convergence rates of n iterations of SAG.

In the examples we take n = 100000, L = 100, µ = 0.01 (Case 1), and µ = 0.0001 (Case 2).

the discussion after Proposition 3.4). Related results also appear in [START_REF] Collins | Exponentiated gradient algorithms for conditional random fields and max-margin markov networks[END_REF], 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 n 2L µ , n 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 1/16L, even though it uses iterations which are n 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 n appears in the convergence rate in the wellconditioned setting where n > 2L µ , if we perform n iterations of SAG (i.e., one effective pass through the data), the error is multiplied by (1 -1/8n) n exp(-1/8), which is independent of n. 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 g satisfying our assumptions. In this table, we compare the rate obtained for these FG methods to the rate obtained by running n iterations of SAG, since this requires the same number of evaluations of f i . Case 1 in this table focuses on a well-conditioned case where the rate of SAG is (1 -1/8n), while Case 2 focuses on an ill-conditioned example where the rate is (1 -µ/16L). Note that in the latter case the O(1/k) rate for the method may be faster.

In Table 1 we see that performing n 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 L in the SAG rates is the Lipschitz constant of the f i functions, while in the FG method we only require that L is an upper bound on the Lipschitz continuity of g 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 of the extended version of this paper 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 2 -regularized least-squares regression.

Implementation Details

In Algorithm 1 we give pseudo-code for an implementation of the basic method, where we use a variable d to track the quantity ( n i=1 y i ). 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.

Algorithm 1 Basic SAG method for minimizing

1 n n i=1 f i (x) with step size α. d = 0, y i = 0 for i = 1, 2, . . . , n for k = 0, 1, . . . do Sample i from {1, 2, . . . , n} d = d -y i + f i (x) y i = f i (x) x = x -α n d end for

Structured gradients and just-in-time parameter updates

For many problems the storage cost of O(np) for the y k i vectors is prohibitive, but we can often use the structure of the gradients f i to reduce this cost. For example, a commonly-used specialization of ( 1) is linearly-parameterized models which take form minimize

x∈R p g(x) := 1 n n i=1 f i (a i x). (9)
Since each a i is constant, for these problems we only need to store the scalar

f i k (u k i ) for u k i = a i k x k rather than the full gradient a i f i (u k i )
. This reduces the storage cost from O(np) down to O(n). For problems where the vectors a i are sparse, an individual gradient f i will inherit the sparsity pattern of the corresponding a i . However, the update of x in Algorithm 1 appears unappealing since in general d will be dense, resulting in an iteration cost of O(p). 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 a i k . In particular, we do this by not explicitly storing the full vector x k after each iteration. Instead, on each iteration we only compute the elements x k j corresponding to non-zero elements of a i k , by applying the sequence of updates to each variable x k j since the last iteration where it was non-zero in a i k . This sequence of updates can be applied efficiently since it simply involves changing the step size. For example, if variable j has been zero in a i k for 5 iterations, then we can compute the needed value x k j using

x k j = x k-5 j - 5α n n i=1 (y k i ) j .
This update allows SAG to be efficiently applied to sparse data sets where n and p are both in the millions or higher but the number of non-zeros is much less than np.

Re-weighting on early iterations

In the update of x in Algorithm 1, we normalize the direction d by the total number of data points n. When initializing with y 0 i = 0 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 y i variables contributing to d are set to the uninformative zero-vector. Following [START_REF] Blatt | A convergent incremental gradient method with a constant step size[END_REF], the more logical normalization is to divide d by m, the number of data points that we have seen at least once (which converges to n once we have seen the entire data set), leading to the update x = x -α m d. 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 k iterations, this modified SAG algorithm is also advantageous over hybrid SG/SAG algorithms because it only requires estimating a single constant step size.

Exact and efficient regularization

In the case of regularized objectives like (2), the cost of computing the gradient of the regularizer is independent of n. Thus, we can use the exact gradient of the regularizer in the update of x, and only use d to approximate the sum of the l i functions. By incorporating the gradient of the regularizer explicitly, the update for y i in Algorithm 1 becomes y i = l i (x), and in the case of 2 -regularization the update for x becomes

x = x -α 1 m d + λx = (1 -αλ) x - α m d.
If the loss function gradients l i 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 x again appears to require dense vector operations, we can implement the algorithm efficiently if the a i are sparse. In particular, to allow efficient multiplication of x by the scalar (1-αλ), it is useful to represent x in the form x = κz, where κ is a scalar and z is a vector (as done by [START_REF] Shalev-Shwartz | Pegasos: primal estimated sub-gradient solver for svm[END_REF]). Under this representation, we can multiply x by a scalar in O(1) by simply updating κ (though to prevent κ becoming too large or too small we may need to occasionally renormalize by setting z = κz and κ = 1). 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 pseudocode for a variant of SAG that includes all of these modifications, and thus uses no full-vector operations. This code uses a vector y to keep track of the scalars l i (u k i ), a vector C to determine whether a data point has previously been visited, a vector V to track the last time a variable was updated, and a vector S to keep track of the cumulative sums needed to implement the just-in-time updates.

Warm starting

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 x 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 x, y i , d, m, and C i based on a previous run of

Algorithm 2 SAG variant for minimizing λ 2 x 2 + 1 n n i=1 l i (a i x), with step size α and a i sparse. {Initialization, note that x = κz.} d = 0, y i = 0 for i = 1, 2, . . . , n z = x, κ = 1 m = 0, C i = 0 for i = 1, 2, . . . , n S -1 = 0, V j = 0 for j = 1, 2, . . . , p for k = 0, 1, . . . do Sample i from {1, 2, . . . , n} if C i = 0 then {This is the first time we have sampled this data point.} m = m + 1 C i = 1 end if for j non-zero in a i do {Just-in-time calculation of needed values of z.} z j = z j -(S k-1 -S Vj -1 )d j V j = k end for
{Update the memory y and the direction d.} Let J be the support of

a i d J = d J -a iJ (y i -l i (κa T iJ z J )) y i = l i (κa T iJ z J ) {Update κ and the sum needed for z updates.} κ = κ(1 -αλ) S k = S k-1 + α/(κm) end for {Final x is κ times the just-in-time update of all z.} for j = 1, 2, . . . , p do x j = κ(z j -(S k-1 -S Vj -1 )d j ) end for
the SAG algorithm. In this scenario, Theorem 1 suggests that it may be beneficial in this setting to center the y i variables around d.

Larger step sizes

In our experiments we have observed that utilizing a step size of 1/L, as in standard FG methods, always converged and often performed better than the step size of 1/16L suggested by our analysis. Thus, in our experiments we used α k = 1/L even though we do not have a formal analysis of the method under this step size. We also found that a step size of 2/(L + nµ), which in the stronglyconvex case corresponds to the best fixed step size for the FG method in the special case of n = 1 [see Nesterov, 2004, Theorem 2.1.15], sometimes yields even better performance (though in other cases it performs poorly).

Line-search when L is not known

In general the Lipschitz constant L will not be known, but we may obtain a reasonable approximation of a valid L by evaluating f i values while running the algorithm. In our experiments, we used a basic line-search where we start with an initial estimate L 0 , and double this estimate whenever we do not satisfy the inequality

f i k (x k - 1 L k f i k (x k )) f i k (x k ) - 1 2L k f i k (x k ) 2 ,
which must be true if L k is valid. An important property of this test is that it depends on f i k but not on g, and thus the cost of performing this test is independent of n. To avoid instability caused by comparing very small numbers, we only do this test when f i k (x k ) 2 > 10 -8 . Since L 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 L (and thus a larger step size) can be used as we approach x * . To potentially take advantage of this, we initialize with the slightly smaller L k = (L k-1 2 -1/n ) at each iteration, so that the estimate of L is halved if we do n iterations (an effective pass through the data) and never violate the inequality. Note that in the case of 2 -regularized objectives, we can perform the line-search to find an estimate of the Lipschitz constant of l i rather than f i , and then simply add λ to this value to take into account the effect of the regularizer.

Note that the cost of this line-search is independent of n, making it suitable for large problems. Further, for linearly-parameterized models of the form f i (a T i x), it is also possible to implement the line-search so that its cost is also independent of the number of variables p. To see why, if we use δ k = a T i k x k and the structure in the gradient then the left side is given by

f i k a T i k x k - 1 L k f i k (x k ) = f i k δ k - f i k (δ k ) L k a i k 2 .
Thus, if we pre-compute the squared norms a i 2 and note that δ k and f i k (δ k ) are already needed by the SAG update, then each iteration only involves operations on scalar values and the single-variable function f i k .

Mini-batches for vectorized computation and reduced storage

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 y i for each mini-batch rather than for each example. Thus, for example, using a mini-batch of size 100 leads to a 100-fold reduction in the storage cost.

A subtle issue that arises when using mini-batches is that the value of L in the Lipschitz condition ( 8) is based on the mini-batches instead of the original functions f i . For example, consider the case where we have a batch B and the minimum value of L in (8) for each i is given by L i . In this case, a valid value of L for the function x → 1 |B| i∈B f i (x) would be max i∈B {L i }. We refer to this as L max . But we could also consider using L mean = 1 |B| i∈B L i . The value L mean is still valid and will be smaller than L max unless all L i are equal. We could even consider the minimum possible value of L, which we refer to as L Hessian because (if each f i is twice-differentiable) it is equal to the maximum eigenvalue of 1 |B| i∈B f i (x) across all x. Note that L Hessian ≤ L mean ≤ L max , although L Hessian will typically be more difficult to compute than L mean or L max (although a line-search as discussed in the previous section can reduce this cost). Due to the potential of using a smaller L, 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.

Non-uniform example selection

In standard SG methods, it is crucial to sample the functions f i 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 1/n, regardless of the frequency at which the corresponding function is sampled. We might thus consider sampling the functions f i 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 f i whose gradient changes slowly as often as functions f i 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 f i and we note that in subsequent works this intuition has proved correct for related algorithms [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF]Zhang, 2014, Schmidt et al., 2015].

Let L i again be the Lipschitz constant of f i , and assume that the functions are placed in increasing order of Lipschitz constants, so that L 1 L 2 . . . L n . In the ill-conditioned setting where the convergence rate depends on µ L , a simple way to improve the rate by decreasing L is to replace f n by two functions f n1 and f n2 such that

f n1 (x) = f n2 (x) = f n (x) 2 g(x) = 1 n n-1 i=1 f i (x) + f n1 (x) + f n2 (x) = 1 n + 1 n-1 i=1 n + 1 n f i (x) + n + 1 n f n1 (x) + n + 1 n f n2 (x) .
We have thus replaced the original problem by a new, equivalent problem where:

• n has been replaced by (n + 1),

• L i for i (n -1) is Li(n+1) n ,
• L n and L n+1 are equal to Ln(n+1) 2n .

Hence, if L n-1 < nLn n+1 , this problem has the same µ but a smaller L than the original one, improving the bound on the convergence rate. By duplicating f n , we increase its probability of being sampled from 1 n to 2 n+1 , but we also replace y k n by a noisier version, i.e. y k n1 + y k n2 . Using a noisier version of the gradient appears detrimental, so we assume that the improvement comes from increasing the frequency at which f n is sampled, and that logically we might obtain a better rate by simply sampling f n more often in the original problem and not explicitly duplicating the data.

We now consider the extreme case of duplicating each function f i a number of times equal to the Lipschitz constant of their gradient, assuming that these constants are integers. The new problem becomes

g(x) = 1 n n i=1 f i (x) = 1 n n i=1 Li j=1 f i (x) L i = 1 k L k n i=1 Li j=1 k L k n f i (x) L i .
The function g is now written as the sum of k L k functions, each with a gradient with Lipschitz constant k L k n

. The new problem has the same µ as before, but now has an L equal to the average of the Lipschitz constants across the f i , rather than their maximum, thus improving the bound on the convergence rate. Sampling these functions uniformly is now equivalent to sampling the original f i '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 [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF] in the context of coordinate descent methods, and is also somewhat related to the sampling scheme used by [START_REF] Strohmer | A randomized kaczmarz algorithm with exponential convergence[END_REF] in the context of their randomized Kaczmarz algorithm. Since the first version of this work was released, [START_REF] Needell | Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm[END_REF] have analyzed sampling according to the Lipschitz constant in the context of SG iterations. Such a sampling scheme makes the iteration cost depend on n, due to the need to generate samples from a general discrete distribution over n variables. However, after an initial preprocessing cost of O(n) we can sample from such distributions in O(log n) 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 α k = n i Li does not seem to converge in general. However, by changing the number of times we duplicate each f i , we can interpolate between the Lipschitz sampling and the uniform sampling. In particular, if each function f i is duplicated L i + c times, where L i is the Lipschitz constant of f i and c a positive number, then the new problem becomes

g(x) = 1 n n i=1 f i (x) = 1 n n i=1 Li+c j=1 f i (x) L i + c = 1 k (L k + c) n i=1 Li+c j=1 k (L k + c) n f i (x) L i + c .
Unlike in the previous case, these

k (L k + c)

Experimental Results

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.

Comparison to FG and SG Methods

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 2 -regularized logistic regression problem minimize

x∈R p λ 2 x 2 + 1 n n i=1 log(1 + exp(-b i a i x)), (10)
as a canonical problem satisfying our assumptions. In our experiments we set the regularization parameter λ to 1/n, 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;1 the covertype (based on the datset of Blackard, Jock, and Dean), rcv1, news, and rcv1Full data sets were obtained from the LIBSVM Data website;2 ; the sido data set was obtained from the Causality Workbench website,3 the spam data set was prepared by [see Carbonetto, 2009, §2.6.5] using the TREC 2005 corpus4 ; and the alpha data set was obtained from the Pascal Large Scale Learning Challenge website5 . 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 l i divided by the number of examples n. 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.

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 [START_REF] Nesterov | A method for unconstrained convex minimization problem with the rate of convergence O(1/k 2 )[END_REF], where iterations of (3) with a step size of 1/L k are interleaved with an extrapolation step. We used an adaptive linesearch to estimate a local L based on the variant proposed for 2 -regularized logistic regression by [START_REF] Liu | Large-scale sparse logistic regression[END_REF].

• L-BFGS: A publicly-available limited-memory quasi-Newton method that has been tuned for log-linear models such as logistic regression [START_REF] Schmidt | minfunc: unconstrained differentiable multivariate optimization in matlab[END_REF]. 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 10 (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 10.6

• IAG: The incremental aggregated gradient method of [START_REF] Blatt | A convergent incremental gradient method with a constant step size[END_REF] described by iteration (5) with a cyclic choice of i k . 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 10.

• 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 α k = 1/L k where L k is an approximation of the Lipschitz constant for the negative log-likelihoods l i (x) = log(1 + exp(-b i a i x)). Although this Lipschitz constant is given by 0.25 max i { a i 2 }, 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 L 0 = 1.

We plot the results of the different methods for the first 50 effective passes through the data in Figure 1. For the stochastic methods, we plot the mean performance as well as the minimum and maximum function values across 10 choices for the initial random seed. 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 50 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.

Comparison to Coordinate Optimization Methods

For the 2 -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 accelerated SG method of [START_REF] Delyon | Accelerated stochastic approximation[END_REF], SG methods that only average the later iterations as in the 'optimal' SG method for non-smooth optimization of [START_REF] Rakhlin | Making gradient descent optimal for strongly convex stochastic optimization[END_REF], the epoch SG method of [START_REF] Hazan | Beyond the regret minimization barrier: an optimal algorithm for stochastic strongly-convex optimization[END_REF], the 'nearly-optimal' SG method of [START_REF] Ghadimi | Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization[END_REF], 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 [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF]. 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. linear convergence rate by applying coordinate descent to the primal [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF] or coordinateascent 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 [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF], using a step-size of 1/L j , where L j is the Lipschitz-constant with respect to coordinate j of g . 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 [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF].

• 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 n (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 n/p so that 1 effective pass for this method has a cost of O(np) 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 n > p (top and bottom rows of Figure 2), there is little difference between uniform and Lipschitz sampling of the coordinates. For the problems with p > n (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 p > n, DCA and PCD-L have similar performance. For the problems with n > p, 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.

Comparison of Step-Size Strategies

In 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 α k = 1/2nL 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 α k = 1/2nµ performs extremely well on the data sets with p > n (middle row). In contrast, for the data sets with n > p 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 {n, µ, L} hold.

• Theorem 1: Using a step-size of α k = 1/16L performs consistently better than the smaller step size α k = 1/2nL, but in some cases it performs worse than α k = 1/2nµ. However, in contrast to α k = 1/2nµ, the step size α k = 1/16L always has reasonable performance.

• Section 4.5: The step size of α k = 1/L performs performs extremely well on the data sets with p > n, and performs better than the step sizes discussed above on all but one of the remaining data sets. The step size of α k = 2/(L + nµ) seems to perform the same or slightly better than using α k = 1/L 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 100 to 10000 times larger than the step sizes chosen by IAG, and always lead to better performance by several orders of magnitude.

Effect of mini-batches

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 L. 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 1/L max , 1/L mean , and 1/L Hessian 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 nµ 2L without loss of convergence speed, this is a very conservative estimate. In our experiments, the original value of nµ L was on the order of 10 -5 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.

Effect of non-uniform sampling

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 L i + c, where L i is the global Lipschitz constant of the corresponding f i and we set c to the average of these constants, c = L mean = (1/n) i L i . Plugging in these values into the formula at the end of Section 4.8 and using L max to denote the maximum L i value, we set the step-size to α k = 1/2L max + 1/2L mean .

• SAG-LS (Lipschitz): In this method we formed estimates of the quantities L i , L max , and L mean . The estimator L k max is computed in the same way as the SAG-LS method. To estimate each L i , we keep track of an estimate L k i for each i and we set L k mean to the average of the L k i values among the f i that we have sampled at least once. We set L k i = L k-1 i if example i was not selected and otherwise we initialize to

L k i = L k-1 i
/2 and perform the line-search until we have a valid L k i (this means that L k mean 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 (n -m)/n, and otherwise we sample from the previously seen f i in proportion to

L k i + L k mean .
To prevent relying too much on our initially-poor estimate of L mean , we use a step size of α k = n-m n α max + m n α mean , where α max = 1/L k max is the step-size we normally use with uniform sampling and α mean = 1/2L k max + 1/2L k mean is the step-size we use with the non-uniform sampling method, so that the method interpolates between these extremes until the entire data set has been sampled.

We make the following observations from these experiments:

• SAG (1/L) vs. SAG (Lipschitz): With access to global quantities and a constant step size, the difference between uniform and non-uniform sampling was typically quite small. However, in some cases the non-uniform sampling method behaved much better (top row of Figure 5).

• SAG-LS vs. SAG-LS (Lipschitz): When estimating the Lipschitz constants of the individual functions, the non-uniform sampling strategy often gave better performance. Indeed, the adaptive non-uniform sampling strategy gave solutions that are orders of magnitude more accurate than any method we examined for several of the data sets (e.g., the protein, covertype, and sido data sets) In the context of logistic regression, it makes sense that an adaptive sampling scheme could lead to better performance, as many correctly-classified data samples might have a very slowly-changing gradient near the solution, and thus they do not need to be sampled often.

Discussion

Since the first version of this work was published [START_REF] Le Roux | A stochastic gradient method with an exponential convergence rate for strongly-convex optimization with finite training sets[END_REF], there has been an explosion of interest in stochastic methods with improved convergence rates. In this section we first review other algorithms that have been discovered to have this property, and then we discuss the many possible variants on these basic algorithms that have been explored. As this is a very quickly-evolving area there are likely to be many new developments in the near future, but we note that this literature review is up to date as of January, 2015.

Alternative Algorithms

SDCA: The first algorithm that was shown to have a similar convergence rate and iteration cost to SAG was in fact a much older algorithm: coordinate optimization applied to a dual problem with randomized coordinate selection, referred to as stochastic dual coordinate ascent (SDCA). Shalev-Schwartz and Zhang [2013b] consider the problem of minimizing an 2 -regularized finite sum, minimize

x∈R p λ 2 x 2 + 1 n n i=1 f i (a T i x),
where each f i is convex and each f i is Lipschitz-continuous, by optimizing its Fenchel dual, maximize

y∈R n - 1 2λ 1 n n i=1 y i a i 2 - 1 n n i=1 f * i (-y i ).
They consider applying exact coordinate optimization to a randomly selected coordinate. Using the primal-dual relationship x = 1 nλ N i=1 x i a i , they show a linear convergence rate in terms of the duality gap. Since there is one dual variable associated with each example i, the iteration cost is independent of n and thus the strategy has similar convergence properties to SAG (the memory requirements are identical in this context, see Section 4.1). This result is related to the earlier result of [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF], who shows a similar convergence rate for randomized coordinate descent.7 However, the result of Nesterov cannot be directly applied to the dual problem in general since the dual does not necessarily have a Lipschitz-continuous gradient. A similar result was also reported even earlier by [START_REF] Collins | Exponentiated gradient algorithms for conditional random fields and max-margin markov networks[END_REF] without requiring that the gradient of the dual is Lipschitz-continuous, but this result again only applied to the dual objective. More recently, Shalev-Schwartz and Zhang [2013a] show linear convergence of SDCA in the more general setting minimize

x∈R p λr(x) + 1 n n i=1 f i (A T i x),
where A i are matrices and r is 1-strongly convex. They also relax the requirement of exact coordinate optimization, providing a variety of more practical alternatives. Further, in subsequent work they obtain a convergence rate in the convex case by adding an explicit strongly-convex regularizer to the problem [START_REF] Shalev-Schwartz | Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization[END_REF].

A disadvantage of SDCA compared to the SAG algorithm is that the SDCA convergence rates depend on λ rather than the strong-convexity constant µ. In the worst case we have µ = λ, but if µ is much larger then the convergence rate of SAG is much faster. Further, even in cases where µ = λ, the convergence rate of SAG might be much faster if the iterates stay in local region with a higher strong-convexity constant. As an extreme example, due to local strong-convexity SAG might have a linear convergence rate in scenarios where SDCA has a sub-linear convergence. This subtle but practically important issue was a key focus in the recent work of [START_REF] Agarwal | A lower bound for the optimization of finite sums[END_REF], and indeed the performance of SDCA was very poor on three of the test problems in our experiments. See Figure 2 (top left, top right, bottom right).

MISO: [START_REF] Mairal | Optimization with first-order surrogate functions[END_REF] analyzes a very general surrogate optimization framework, that includes a wide vareity of existing algorithms. He also considers incremental algorithms in this framework, and specialized to the smooth and unconstraeind setting (with a 'Lipschitz surrogate') obtains an algorithm (MISO) that is very similar to the SAG algorithm,

x k+1 = 1 n n i=1 x k i - α k n n i=1 y k i ,
where y k i is defined as in the SAG algorithm 6, and x k i is the parameter vector used to compute the corresponding y k i . Thus, instead of applying the SAG step to x k , MISO applies the step to the average of the previous x k i values used to form the current y k i variables. [START_REF] Mairal | Optimization with first-order surrogate functions[END_REF] shows that this algorithm also achieves an O(1/k) rate for convex objectives and a linear convergence rate for strongly-convex objectives.

However, MISO has the disadvantage that it not only requires storing the n gradient values but also storing n previous iterations (which are less likely to have a nice structure). Further, the linear convergence rate shown for MISO is substantially slower than the convergence rate shown in Theorem 1. In particular, the rate is more similar to the substantially slower Proposition 1 in our prior work [START_REF] Le Roux | A stochastic gradient method with an exponential convergence rate for strongly-convex optimization with finite training sets[END_REF]. Subsequent work on the MISO algorithm has shown an analogous result to Proposition 2 in our prior work; if the n is sufficiently larger than Lµ, then using a step-size proportional to 1/µ yields a faster convergence rate [START_REF] Mairal | Incremental majorization-minimization optimization with application to large-scale machine learning[END_REF], Defazio et al., 2014b]. However, using this step-size causes divergence if µ is not sufficiently large. SVRG: Another interesting framework that has been considered is known as 'mixed optimization', 'stochastic variance-reduced gradient' (SVRG), or 'semi-stochastic gradient descent' (S2GD) [START_REF] Mahdavi | Mixedgrad: An o(1/t) convergence rate algorithm for stochastic smooth optimization[END_REF][START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF], Zhang et al., 2013[START_REF] Konečnỳ | Semi-stochastic gradient descent methods[END_REF]. 8Unlike FG methods that utilize the full gradient g and SG methods that consider individual gradients f i , these mixed optimization methods combine the two. In particular, they use a (possibly regularized) iteration of the form

x k+1 = x k -α k (f i (x k ) -f i (x k ) + 1 n n j=1 f j (x k )),
where xk is the last iterate where the full gradient g was evaluated. The algorithm alternates between computing the full gradient at xk , and performing some number m of stochastic gradient iterations. Note this is very similar to the SAG algorithm written in the form

x k+1 = x k -α k ( 1 n f i (x k ) - 1 n f i (x k i ) + 1 n n j=1 f j (x k i )),
where x k i is defined as in the MISO algorithm. These methods are thus similar to SAG in the use of potentially-outdated gradient information, but differ in that the outdated gradients are all computed at the same previous iteration xk and the weighting of terms is changed. If the step-size and parameter m are set appropriately, this algorithm has a linear convergence rate in the stronglyconvex case. It can also achieve an O(1/k) convergence rate in the general convex case, by applying it to a purturbed problem where a strongly-convex regularizer has been added.

A key advantage of this strategy is it only requires storing xk , rather than the n gradient values. However, to obtain this improved memory requirement we must evaluate f i twice on each iteration. Further, the algorithm innefficiently requires full passes through the data to evaluate f j (x k ) for all j.

A disadvantage of the SVRG method is that it requires setting two parameters (rather than one) and the convergence rate depends in a non-trivial way on their interaction both with each other and with both the Lipschitz constant and the strong-convexity. As with SDCA, the dependence on the strong-convexity constant is notable. In particular, the algorithm can't be applied without modification to general convex problems, and (although the effect is not as severe as it is with SDCA) the convergence rate may be substantially slower than SAG if there exists hidden strong convexity.

Self-Concordant Objectives: A surprising recent development is that [START_REF] Bach | Non-strongly-convex smooth stochastic approximation with convergence rate o(1/n)[END_REF] have shown that the finite sum assumption is not required to obtain the O(1/k) convergence rate in the special case of least squares. They show that an averaged SG method achieves this rate, and give a 2-phase Newton-like SG method that achieves this rate under an assumption similar to self-concordance. However, this assumption does not hold in general for the class of problems considered here.

SAGA: One of the most recent developments in this area is the SAGA algorithm Defazio et al. [2014a]. This algorithm intelligently combines the updates used in the SAG and SVRG algorithms. This method maintains the appealing properties of SAG, but yields a simpler proof. However, the proof technique used in that work does not yield a simpler analysis of the original SAG algorithm.

Lower Bounds: There has been recent work on determining a lower bound on the convergence rate that can be expected for minimizing finite sums. Defazio et al. [2014b] show that the rate must be at least (1 -1/n) in the strongly-convex case, while [START_REF] Agarwal | A lower bound for the optimization of finite sums[END_REF] establish a bound that also depends on the condition number L/µ.

Generalizations and Other Issues

Accelerated gradient: AFG methods are variants of the basic FG method that can obtain an O(1/k 2 ) convergence rate for the general convex case, and faster linear convergence rates that depend on the square root of the condition number ( L/µ) rather than the condition number (L/µ) in the strongly-convex case [see Nesterov, 2004, §2.2.1]. For strongly-convex objectives, it has been shown that a mini-batch strategy can obtain a better dependence on the condition number in certain regimes for SDCA [START_REF] Shalev | Accelerated mini-batch stochastic dual coordinate ascent[END_REF] and more recently for SVRG [START_REF] Nitanda | Stochastic proximal gradient descent with acceleration techniques[END_REF]. It is possible that similar arguments could hold for SAG algorithms, which could be advantageous over these methods for reasons discussed in the previous section. [START_REF] Shalev-Schwartz | Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization[END_REF] have also given an accelerated version of the SDCA method for ill-conditioned problems, that uses SDCA to solve a sequence of regularized problems up to a prescribed optimality. Since this procedure ultimately relies on the sequence of primal solutions, we can also accelerate SAG using a procedure like this. However, a difficulty with this procedure (whether SAG or SDCA are used) is the cost of measuring the optimality of the sub-problems.

Rather than using this inner-outer procedure, [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF] show that using a deterministic primal iteration (based on the full data set) allows one to construct a primal solution that has the accelerated rate from the result of an accelerated dual coordinate ascent method. [START_REF] Zhang | Stochastic primal-dual coordinate method for regularized empirical risk minimization[END_REF] give a coordinate-wise variant of the accelerated primal-dual method of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] that also achieves this rate. Based on these results, it is possible that accelerated versions of SAG could be developed that do not rely on an inner-outer procedure.

The convergence rates of these accelerated methods have the same form as the lower bound established by [START_REF] Agarwal | A lower bound for the optimization of finite sums[END_REF]. However, Agarwal and Bottou point out that these accelerated convergence rates depend on λ rather than µ. Thus, they conclude that the basic SAG algorithm may still be much faster on some problems than accelerated SDCA methods, and indeed show how to construct a problem where SAG is arbitrarily faster than accelerated SDCA. The possibility of developing an accelerated SAG algorithm that is adaptive to µ remains open.

Proximal gradient and ADMM: It is becoming increasingly common to address problems of the form minimize

x∈R p r(x) + g(x) := r(x) + 1 n n i=1 f i (x),
where f i and g satisfy our assumptions and r is a general convex function that could be non-smooth or enforce that constraints are satisfied. Proximal-gradient methods for problems with this structure use iterations of the form

x k+1 = prox α k x k - α k n n i=1 f i (x k ) ,
where the prox α k [y] operator is the solution to the proximity problem minimize

x∈R p 1 2 x -y 2 + α k r(x).
Proximal-gradient and accelerated proximal-gradient methods are appealing for solving non-smooth optimization problems because they achieve the same convergence rates as FG methods for smooth optimization problems [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF][START_REF] Schmidt | Convergence rates of inexact proximal-gradient methods for convex optimization[END_REF]. We have explored a variant of the proximal-gradient method where the average over the f i (x k ) values is replaced by the SAG approximation (6). Although our analysis does not directly apply to this scenario, we believe that this proximal-SAG algorithm for composite non-smooth optimization achieves the same convergence rates as the SAG algorithm for smooth optimization; this is supported by the experiments of [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF]. Indeed, there now exist proximal-gradient variants of SDCA [Shalev-Schwartz and Zhang, 2013a], MISO [START_REF] Mairal | Optimization with first-order surrogate functions[END_REF], SVRG [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF], and SAGA [Defazio et al., 2014a].

In cases where r is composed with a linear function, we can consider approaches based on the alternating direction method of multipliers (ADMM). This has been explored for SDCA [START_REF] Suzuki | Stochastic dual coordinate ascent with alternating direction method of multipliers[END_REF] and MISO [START_REF] Zhong | Fast stochastic alternating direction method of multipliers[END_REF], and a variant based on SAG is also likely to be possible.

Coordinate-wise:

The key advantage of SG and SAG methods is that the iteration cost is independent of the number of functions n. However, in many applications we may also be concerned with the dependence of the iteration cost on the number of variables p. Randomized coordinate-wise methods offer linear convergence rates with an iteration cost that is linear in n but independent of p [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF]. We can consider a variant of SAG whose iteration cost is independent of both n and p by using the update

[y k i ] j = [f i (x k )] j if i = i k and j = j k [y k-1 i ] j otherwise,
to each coordinate j of each vector y i in the SAG algorithm, and where j k is a sample from the set {1, 2, . . . , p}. [START_REF] Konečnỳ | Semi-stochastic coordinate descent[END_REF] recently proposed an SVRG algorithm of this flavour.

Newton-like:

In cases where g is twice-differentiable, we can also consider Newton-like variants of the SAG algorithm,

x k+1 = x k - α k n B k [ n i=1 y k i ],
where B k is a positive-definite approximation to the inverse Hessian g (x k ). We would expect to obtain a faster convergence rate by using an appropriate choice of the sequence {B k }. However, in order to not increase the iteration cost these matrices should be designed to allow fast multiplication.

For example, we could choose B k to be diagonal, which would also preserve any sparsity present in the updates. [START_REF] Sohl-Dickstein | Fast large-scale optimization by unifying stochastic gradient and quasi-newton methods[END_REF] propose a quasi-Newton method in this class that shows impressive empirical results, although the iteration cost is much higher.

Relaxing Convexity Assumptions: It is likely that the convexity assumptions made in this work could be relaxed. For example, [START_REF] Gong | Linear convergence of variance-reduced projected stochastic gradient without strong convexity[END_REF] show SVRG can obtain a linear convergence under weaker assumptions. Since SAG is adaptive to hidden strong-convexity, the assumptions needed for its linear convergence are likely even weaker. Further, the SAG algorithm may also be useful even if g is non-convex (which is different than SDCA). In this case we expect that, similar to the IAG method [START_REF] Blatt | A convergent incremental gradient method with a constant step size[END_REF], the algorithm converges to a stationary point under very general conditions.

Non-Uniform Sampling: We have given an argument that non-uniform sampling should benefit the SAG algorithm, and shown empirically that it can lead to a substantial improvement. However, we have not yet given a full analysis of this scheme. Subsequent works have shown that the type of dependency we conjecture here (e.g., dependence on the average Lispschitz constant) can be achieved with non-uniform sampling in the context of SDCA [Qu et al., 2014, Zhao andZhang, 2014], SVRG [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF], and SAGA [START_REF] Schmidt | Non-uniform stochastic average gradient method for training conditional random fields[END_REF] Step-size selection and termination criteria: The three major disadvantages of SG methods are: (i) the slow convergence rate, (ii) deciding when to terminate the algorithms, and (iii) choosing the step size while running the algorithm. This work shows that the SAG iterations achieve a much faster convergence rate, but the SAG iterations may also be advantageous in terms of termination criteria and choosing step sizes. In particular, the SAG iterations suggest a natural termination criterion; since the quantity d in Algorithm 1 converges to f (x k ) as x k -x k-1 converges to zero, we can use d as an approximation of the optimality of x k as the iterates converge. Regarding choosing the step-size, a disadvantage of a constant step-size strategy is that a step-size that is too large may cause divergence. But, we expect that it is possible to design line-search or trust-region strategies that avoid this issue. Such strategies might even lead to faster convergence for functions that are locally well-behaved around their optimum, as indicated in our experiments. Further, while SG methods require specifying a sequence of step sizes and mis-specifying this sequence can have a disastrous effect on the convergence rate [see Nemirovski et al., 2009, §2.1], our theory shows that the SAG iterations achieve a fast convergence rate for any sufficiently small constant step size, and our experiments indicate that a simple line-search gives strong performance.

Appendix A: Comparison of convergence rates

We now present a comparison of the convergence rates of primal and dual FG and coordinate-wise methods to the rate of the SAG method in terms of effective passes through the data for 2 -regularized least-squares regression. In particular, in this appendix we will consider the problem minimize

x∈R p g(x) := λ 2 x 2 + 1 2n n i=1 (a T i x -b i ) 2 ,
where to apply the SG or SAG method we can use

f i (x) := λ 2 x 2 + 1 2 (a T i x -b i ) 2 .
If we use b to denote a vector containing the values b i and A to denote a matrix withs rows a i , we can re-write this problem as

minimize x∈R p λ 2 x 2 + 1 2n Ax -b 2 .
The Fenchel dual of this problem is minimize

y∈R n d(y) := n 2 y 2 + 1 2λ y AA y + y b.
We can obtain the primal variables from the dual variables by the formula x = (-1/λ)A y. Convergence rates of different primal and dual algorithms are often expressed in terms of the following Lipschitz constants:

L g = λ + M σ /n (Lipschitz constant of g ) L i g = λ + M i (Lipschitz constant for all f i ) L j g = λ + M j /n (Lipschitz constant of all g j ) L d = n + M σ /λ (Lipschitz constant of d ) L i d = n + M i /λ (Lipschitz constant of all d i )
Here, we use M σ to denote the maximum eigenvalue of A A, M i to denote the maximum squared row-norm max i { a i 2 }, and M j to denote the maximum squared column-norm max j { n i=1 (a i ) 2 j }. We use g j to refer to element of j of g , and similarly for d i . The convergence rates will also depend on the primal and dual strong-convexity constants:

µ g = λ + m σ /n (Strong-convexity constant of g) µ d = n + m σ /λ (Strong-convexity constant of d)
Here, m σ is the minimum eigenvalue of A A, and m σ is the minimum eigenvalue of AA .

A.1 Full Gradient Methods

Using a similar argument to [Nesterov, 2004, Theorem 2.1.15], if we use the basic FG method with a step size of 1/L g , then (f (x k ) -f (x * )) converges to zero with rate

1 - µ g L g 2 = 1 - λ + m σ /n λ + M σ /n 2 = 1 - nλ + m σ nλ + M σ 2 exp -2 nλ + m σ nλ + M σ ,
while a larger step-size of 2/(L g + µ g ) gives a faster rate of

1 - µ g + µ g L g + µ g 2 = 1 - nλ + m σ nλ + (M σ + m σ )/2 2 exp -2 nλ + m σ nλ + (M σ + m σ )/2
, and we see that the speed improvement is determined by how much smaller m σ is than M σ .

If we use the basic FG method on the dual problem with a step size of 1/L d , then (d(x k ) -d(x * )) converges to zero with rate

1 - µ d L d 2 = 1 - n + m σ /λ n + M σ /λ 2 = 1 - nλ + m σ nλ + M σ 2 exp -2 nλ + m σ nλ + M σ ,
and with a step-size of 2/(L d + µ d ) the rate is

1 - µ d + µ d L d + µ d 2 = 1 - nλ + m σ nλ + (M σ + m σ )/2 2 exp -2 nλ + m σ nλ + (M σ + m σ )/2
.

Thus, whether we can solve the primal or dual method faster depends on m σ and m σ . In the over-determined case where A has independent columns, a primal method should be preferred. In the under-determined case where A has independent rows, we can solve the dual more efficiently. However, we note that a convergence rate on the dual objective does not necessarily yield the same rate in the primal objective. If A is invertible (so that m σ = m σ ) or it has neither independent columns nor independent rows (so that m σ = m σ = 0), then there is no difference between the primal and dual rates.

The AFG method achieves a faster rate. Applied to the primal with a step-size of 1/L g it has a rate of [Nesterov, 2004, Theorem 2.2.2] 1 -

µ g L g = 1 - λ + m σ /n λ + M σ /n = 1 - nλ + m σ nλ + M σ exp - nλ + m σ nλ + M σ ,
and applied to the dual with a step-size of 1/L d it has a rate of

1 - µ d L d = 1 - n + m σ λ n + M σ /λ = 1 - nλ + m σ nλ + M σ exp - nλ + m σ nλ + M σ .

A.2 Coordinate-Descent Methods

The cost of applying one iteration of an FG method is O(np). For this same cost we could apply p iterations of a coordinate descent method to the primal, assuming that selecting the coordinate to update has a cost of O(1). If we select coordinates uniformly at random, then the convergence rate for p iterations of coordinate descent with a step-size of 1/L j g is [Nesterov, 2010, Theorem 2] 1 -

µ g pL j g p = 1 - λ + m σ /n p(λ + M j /n) p = 1 - nλ + m σ p(nλ + M j ) p exp - nλ + m σ nλ + M j .
Here, we see that applying a coordinate-descent method can be much more efficient than an FG method if M j << M σ . This can happen, for example, when the number of variables p is much larger than the number of examples n. Further, it is possible for coordinate descent to be faster than the AFG method if the difference between M σ and M j is sufficiently large.

Recall that the SAG algorithm performs the recursion (for k 1):

x k = x k-1 - α n n i=1 y k i ,
where an integer i k is selected uniformly at random from {1, . . . , n} and we set

y k i = f i (x k-1 ) if i = i k , y k-1 i otherwise.
We will use the notation

y k =    y k 1 . . . y k n    ∈ R np , θ k =      y k 1 . . . y k n x k      ∈ R (n+1)p , θ * =      f 1 (x * ) . . . f n (x * ) x *      ∈ R (n+1)p ,
and we will also find it convenient to use

e =    I . . . I    ∈ R np×p , f (x) =    f 1 (x) . . . f n (x)    ∈ R np . (11)
With this notation, note that g (x) = 1 n e f (x) and x k = x k-1 -α n e y k . For a square n × n matrix M , we use diag(M ) to denote a vector of size n composed of the diagonal of M , while for a vector m of dimension n, Diag(m) is the n × n diagonal matrix with m on its diagonal. Thus Diag(diag(M )) is a diagonal matrix with the diagonal elements of M on its diagonal, and diag(Diag(m)) = m.

In addition, if M is a tp × tp matrix and m is a tp × p matrix, then we will use the convention that:

• diag(M ) is the tp × p matrix being the concatenation of the t (p × p)-blocks on the diagonal of M ;

• Diag(m) is the tp × tp block-diagonal matrix whose (p × p)-blocks on the diagonal are equal to the (p × p)-blocks of m.

Finally, F k will denote the σ-field of information up to (and including) time k. In other words, F k is the σ-field generated by i 1 , . . . , i k . Given F k-1 , we can write the expected values of the y k and x k variables in the SAG algorithm as

E(y k |F k-1 ) = 1 - 1 n y k-1 + 1 n f (x k-1 ), E(x k |F k-1 ) = x k-1 - α n 1 - 1 n e y k-1 - α n 2 e f (x k-1 ).
The proof is based on finding a Lyapunov function L from R (n+1)p to R such that the sequence EL(θ k ) decreases at an appropriate rate, and EL(θ k ) dominates [g(x k ) -g(x * )]. We derive these results in a parameterized way, leaving a variety of coefficients undetermined but tracking the constraints required of the coefficients as we go. Subsequently, we guide the setting of these coefficients by using a second-order cone solver, and verify the validity of the resulting coefficients using a symbolic solver to check positivity of certain polynomials. Finally, the constants in the convergence rate are given by the initial values of the Lyapunov function, based on the choice of y 0 .

The Lyapunov function contains a term of the form (θ

k -θ * ) A B B C (θ k -θ * )
E (θ k -θ * ) A B B C (θ k -θ * ) F k-1 = (y k-1 -f (x * )) 1 - 2 n S + 1 n Diag(diag(S)) (y k-1 -f (x * )) + 1 n (f (x k-1 ) -f (x * )) Diag(diag(S))(f (x k-1 ) -f (x * )) + 2 n (y k-1 -f (x * )) [S -Diag(diag(S))] (f (x k-1 ) -f (x * )) + 2 1 - 1 n (y k-1 -f (x * )) B - α n eC (x k-1 -x * ) + 2 n (f (x k-1 ) -f (x * )) B - α n eC (x k-1 -x * ) + (x k-1 -x * ) C(x k-1 -x * ) , with S = A - α n Be - α n eB + α 2 n 2 eCe ∈ R np×np .
Our proof also uses the following lemma, giving the inverse of a highly-structured matrix that arises in the analysis.

Lemma 2. Let I be the identity in R np×np and e ∈ R np×p be defined by stacking identity matrices in R n×n as in (11). If α and β are non-zero real-valued scalars, then it holds that α I -

1 n ee T + β 1 n ee T -1 = 1 α I - 1 n ee T + 1 β 1 n ee T .
Proof. It is sufficient to verify that the inverse of the left side times the right side is equal to the identity matrix, α I -

1 n ee T + β 1 n ee T 1 α I - 1 n ee T + 1 β 1 n ee T = I - 1 n ee T I - 1 n ee T + α β I - 1 n ee T 1 n ee T + β α 1 n ee T I - 1 n ee T + 1 n ee T 1 n ee T = I - 2 n ee T + 1 n ee T + α β 1 n ee T - 1 n ee T + β α 1 n ee T - 1 n ee T + 1 n ee T = I.
where we use that e T e = nI, giving 1 n ee T 1 n ee T = 1 n ee T .

B.2 General Lyapunov function

For some h 0, we consider a Lyapunov function of the form

L(θ k ) = 2hg(x k + de y k ) -2hg(x * ) + (θ k -θ * ) A B B C (θ k -θ * )
with A = a 1 ee + a 2 I, B = be, C = cI.

To achieve the desired convergence rate, our goal is to show for appropriate values of δ 0 and γ 0 that

B.3 Lyapunov upper bound

To show (a), we derive an upper bound on the quantity

E L(θ k )|F k-1 -(1 -δ)L(θ k-1 ) = 2hE g(x k + de y k ) -2hg(x * ) -(1 -δ)(2hg(x k-1 + de y k-1 ) -2hg(x * )) +E (θ k -θ * ) A B B C (θ k -θ * ) -(1 -δ)(θ k-1 -θ * ) A B B C (θ k-1 -θ * ).
Lemma 1 gives an expression for the expectation over the quadratic term in the last line, and we will have

S = a 2 I + a 1 -2 α n b + α 2 n 2 c ee Diag(diag(S)) = a 2 + a 1 -2 α n b + α 2 n 2 c I S -Diag(diag(S)) = a 1 -2 α n b + α 2 n 2 c (ee -I).
We also use that strong convexity of g implies 2g(x k-1 + de y k-1 ) 2g(x k-1 ) + 2dg (x k-1 ) e y k-1 + µd 2 e y k-1 2 . (12)

Further, using the Lipschitz-continuity of g and the identity x k + de y k = x k-1 + (d -α n )e y k , followed by applying Lemma 1 using A = ee to expand the last term, gives us the bound

E 2g(x k + de y k )|F k-1 2g(x k-1 ) + 2 d - α n g (x k-1 ) E e y k |F k-1 + L d - α n 2 E e y k 2 |F k-1 = 2g(x k-1 ) + 2 d - α n g (x k-1 ) (1 - 1 n )e y k-1 + 1 n e f (x k-1 ) + L d - α n 2 (y k-1 -f (x * )) 1 - 2 n ee + 1 n I (y k-1 -f (x * )) + L d - α n 2 1 n f (x k-1 ) -f (x * ) 2 + 2 n (y k-1 -f (x * )) ee -I (f (x k-1 ) -f (x * )) . (13)
Combining Lemma 1 with Inequalities ( 12) and ( 13) yields an upper bound on E L(θ k )|F k-1 -(1 -δ)L(θ k-1 ) with a large number of terms. To help in the process of simplifying these terms, Table B.3 lists all the terms up to a scalar factor s and possibly a matrix M . The table also gives these scalars and matrices for each term, as well as the source of the term.

We combine the expressions from Table B.3 using the stated groups in the last column. For example, we add together all the terms in group 0 to obtain a scalar coefficient B 0 for terms in this group.

We do this in the straightforward way (using that g (x) = 1 n e f (x) and g (x * ) = 0) for groups 0, 1, 2, 5, 8, and 9. For groups 6 and 7 we split into terms that have an identity matrix M (B 6 ) and terms with an ee term (B 7 ). Similarly, B 3 comes from terms with an identity matrix in M and B 4 correspond to terms with an ee term. Combining expressions in this way (and adding/subtracting (B 3 /n)ee ) gives: E L(θ k )|F k-1 -(1 -δ)L(θ k-1 ) B 0 g(x k-1 ) -g(x * ) + B 9 x k-1 -x * 2 + (x k-1 -x * ) g (x k-1 )B 1 -B 2 g (x k-1 ) 2 ).

B.4 Domination of g(x k ) -g(x * )

By again using the strong convexity of g (as in ( 12)), we have n na 1 + a 2 + nµd 2 2dg (x k ) + 2b(x k -x * ) 2 by minimizing with respect to y k ,

x k -x * 2 L 2 (2h -γ) + c - n na 1 + a 2 + nµd 2 dL + b 2 ,
using the smoothness of g and the assumption that 2h -γ 0.

Thus, in order for L(θ k ) to dominate g(x k ) -g(x * ) we require the additional constraints

L 2 (2h -γ) + c - n na 1 + a 2 + nµd 2 dL + b 2 0 na 1 + a 2 + nµd 2 0 2h -γ 0.

B.5 Finding constants

Given the algorithm and result parameters (α, γ, δ), we have the following constraints: We also require C 1 0, but this will follow from Constraint 8 since we will have C 3 0. All of the constraints above are convex in a 1 , a 2 , b, c, d, h. Thus, given candidate values for the remaining values, the feasibility of these constraints may be checked using a numerical toolbox (as a secondorder cone program). However, these parameters should be functions of n and should be valid for all µ. Thus, given a sampling of values of µ and n, representing these parameters as polynomials in 1/n, the candidate functions may be found through a second-order cone programs.

The other two expressions (Constraints 8 and 9) are more complicated as there are second-order polynomials in y. If n 5, they have positive second derivatives, negative derivatives at y = 0 and y = 1, and positive values at y = 1. They are thus positive, which is verified by the code below.

For the remaining scenarios where n ∈ {2, 3, 4}, we have a fixed x so we can check the positivity of the polynomials in y. A Matlab script that does the steps above in addition to checking these not-particularly-interesting cases is available from the first author's web page.

B.7 Convergence Rate

We have that g(x k ) -g(x * ) L(θ k ),

and that

E(L(θ k )|F k-1 ) -(1 -δ)L(θ k-1 ) - 1 32n (x k-1 -x * ) g (x k-1 ) 0.
In the strongly-convex case (δ > 0), combining these gives us the linear convergence rate E(g(x k )) -g(x * ) (1 -δ) k L(θ 0 ).

In the convex case (µ = 0 and δ = 0), we have by convexity that -1 32n (x k-1 -x * ) g (x k-1 ) 1 32n

[g(x * ) -g(x k-1 )],

We then have 1 32n

[g(x k-1 ) -g(x * )] L(θ k-1 ) -E(L(θ k )|F k-1 ).

  functions have gradients with different Lipschitz constants. Denoting L = max i L i , the maximum Lipschitz constant is equal to k (L k +c)

Figure 1 :

 1 Figure 1: Comparison of different FG and SG optimization strategies. The top row gives results on the quantum (left), protein (center) and covertype (right) datasets. The middle row gives results on the rcv1 (left), news (center) and spam (right) datasets. The bottom row gives results on the rcv1Full (left), sido (center), and alpha (right) datasets. This figure is best viewed in colour.

Figure 2 :

 2 Figure 2: Comparison of optimization different FG and SG methods to coordinate optimization methods.The top row gives results on the quantum (left), protein (center) and covertype (right) datasets. The middle row gives results on the rcv1 (left), news (center) and spam (right) datasets. The bottom row gives results on the rcv1Full (left), sido (center), and alpha (right) datasets. This figure is best viewed in colour.

Figure 3 :

 3 Figure 3: Comparison of step size strategies for the SAG method. The top row gives results on the quantum (left), protein (center) and covertype (right) datasets. The middle row gives results on the rcv1 (left), news (center) and spam (right) datasets. The bottom row gives results on the rcv1Full (left), sido (center), and alpha (right) datasets. This figure is best viewed in colour.

Figure 4 :Figure 5 :

 45 Figure 4: Sub-optimality as a function of the number of effective passes through the data for various datasets, step-size selection schemes and mini-batch sizes. The datasets are quantum (top), covertype (middle) and protein (bottom). Left: the step-size is 1/L with L the maximum Lipschitz constant of the individual gradients. It is thus the same for all mini-batch sizes. Center: the step-size is 1/L where L is obtained by taking the maximum among the averages of the Lipschitz constants within mini-batches. Right: the step-size is 1/L where L is obtained by computing the maximum eigenvalue of the Hessian for each mini-batch, then taking the maximum of these quantities across mini-batches.

  for some values of A, B and C. Our analysis makes use of the following lemma, derived in[Le Roux et al., 2012, Appendix A.4], showing how this quadratic form evolves through the SAG recursion in terms of the elements of θ k-1 .Lemma 1. If P = A B B C , for A ∈ R np×np , B ∈ R np×p and C ∈ R p×p , then

  k ) γ g(x k ) -g(x * ) , almost surely. Thus, in addition to the algorithm parameter α, there are 2 parameters of the result {γ, δ} and 6 parameters of the Lyapunov function {a 1 , a 2 , b, c, d, h}.

-

  (y k-1 -y * ) B 3 (I -1 n ee ) + B 4 1 n ee (y k-1 -y * ) +(y k-1 -y * ) B 5 e(x k-1 -x * ) + B 6 (f (x k-1 ) -f (x * )) + B 7 eg (x k-1 ) +B 8 f (x k-1 ) -f (x * ) 2 , -(1 -δ)L(θ k-1

L

  (θ k ) -γ g(x k ) -g(x * ) (2h -γ) g(x k ) -g(x * ) + c x k -x * 2 +(y k -y * ) a 2 (I -1 + a 2 + nµd 2 ) (y k -y * ) +(y k -y * ) e 2dg (x k ) + 2b(x k -x * ) (2h -γ) g(x k ) -g(x * ) + c x k -x * 2 -1 4

Table 2 :

 2 Binary data sets used in the experiments.

	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	[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]
	alpha	500 000	500	Synthetic

http://osmot.cs.cornell.edu/kddcup

http://www.csie.ntu.edu.tw/ ~cjlin/libsvmtools/datasets

http://www.causality.inf.ethz.ch/home.php

http://plg.uwaterloo.ca/ ~gvcormac/treccorpus

http://largescale.ml.tu-berlin.de

Note that we also compared to a variety of other SG methods including the popular Pegasos SG method of[START_REF] Shalev-Shwartz | Pegasos: primal estimated sub-gradient solver for svm[END_REF], SG with momentum, SG with gradient averaging, the regularized dual averaging method of[START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF] (a stochastic variant of the primal-dual subgradient method of[START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] for regularized objectives), the

Local linear convergence rates of deterministic coordinate descent methods had been established much earlier[START_REF] Luo | On the convergence of the coordinate descent method for convex differentiable minimization[END_REF].

[START_REF] Wang | Variance reduction for stochastic gradient optimization[END_REF] consider a related algorithm that maintains an easily-computable approximation to the current gradient g (x k ). Although this can improve the constants in the sublinear SG convergence rates, it does not improve the rates.

Acknowledgements

We would like to thank the anonymous reviewers for their many useful comments. This work was partially supported by the European Research Council (SIERRA-ERC-239993) and a Google Research Award. Mark Schmidt is also supported by a postdoctoral fellowship from the Natural Sciences and Engineering Research Council of Canada.

For the O(np) cost of one iteration of the FG method, we could alternately perform n iterations of coordinate descent on the dual problem. With a step size of 1/L i d this would obtain a rate on the dual objective of

which will be faster than the dual FG method if M i << M σ . This can happen, for example, when the number of examples n is much larger than the number of variables p. The difference between the primal and dual coordinate methods depends on M i compared to M j and m σ compared to m σ .

The analysis of SDCA gives a convergence rate in terms of the duality gap (and hence the primal) rather than simply in terms of the dual Shalev-Schwartz and Zhang [2013b]. Using the SDCA analysis, we obtain a rate in the duality gap of

which is the same as the dual rate given above when m σ = 0, and is slower otherwise.

A.3 Stochastic Average Gradient

For the O(np) cost of one iteration of the FG method, we can perform n iterations of SAG. With a step size of 1/16L, performing n iterations of the SAG algorithm has a rate of

In the case where n 2L i g /µ g , this is most similar to the rate obtained with the dual coordinate descent method, though there is a constant factor of 16 and the rate depends on the primal strong convexity constant m σ rather than the dual m σ . However, in this case the SAG rate will often be faster because the term nλ in the denominator is replaced by λ. An interesting aspect of the SAG rate is that unlike other methods the convergence rate of SAG reaches a limit: the convergence rate improves as n grows and as the condition number L i g /µ g decreases but no further improvement in the convergence rate is obtained beyond the point where n = 2L i g /µ g . Thus, while SAG may not be the method of choice for very well-conditioned problems, it is a robust choice as it will always tend to be among the best methods in most situations.

Appendix B: Proof of the theorem

In this Appendix, we give the proof of Theorem 1.

B.1 Problem set-up and notation

Recall that we use g = 1 n n i=1 f i to denote a µ-strongly convex function, where the functions f i for i = 1, . . . , n are convex functions from R p to R with L-Lipschitz continuous gradients. In this appendix we will use the convention that µ 0 so that the regular convex case (where µ = 0) is allowed. We assume that a minimizer of g is attained by some parameter x * (such a value always exists and is unique when µ > 0). with

In this expression, y k-1 -y * only appears through a quadratic form, with a quadratic term -(y k-1 -

0 and B 4 0, then this quadratic term is non-positive and we may maximize in closed form with respect to y k-1 to obtain an upper bound. Assuming that B 3 = 0 and B 4 = 0, we use Lemma 2, B 3 (I -

In order to show the decrease of the Lyapunov function, we need to show that

If we assume that C 1 0 and C 2 0, then we have by strong-convexity

and thus the condition is that

And if we want to show that

By experimenting with this strategy, we have come up with the following values for the constants:

.

We will assume n > 1, since the result for n = 1 follows because SAG is equivalent to gradient descent in this case. Under this assumption, we see that Constraints 1 and 2 are satisfied by the above parameterizaiton. Below we verify that Constraints 3-9 are also satisfied using symbolic computations in Matlab (the code used in this section is available on the first author's webpage). Specifically, we parameterize the constraints as polynomials and then verify that the polynomials are positive over an appropriate interval using the function below (which simply checks that no roots of the polynomial P lie in the interval (x 1 , x 2 ), and that P is positive at the mid-point of the interval).

B.6 Verifying the result

To verify the result under these constants, we consider whether δ = 1/8n or µ/16L. In B 4 , we discard the term of the form (1-δ)µhd 2 , which does not impact the validity of our result. Moreover, without loss of generality, we may assume L = 1 and µ ∈ [0, 1].

B.6.1 Well-conditioned problems (µ 2/n)

In this situation, it suffices to show the result for µ = 2/n since, (a) if a function is µ-strongly convex, then it is µ -strongly convex for all smaller µ , and (b) the final inequality does not involve µ. All constraints (when properly multiplied where appropriate by B 3 B 4 ) can be written as rational functions of x = 1/n:

We only need to check the positivity of these polynomials in x. We do this using the function above, via the script below. Note that we also have B 3 > 0 and B 4 > 0 for n > 1.

B.6.2 Ill-conditioned problems (µ 2/n)

We consider the variables x = 1/n ∈ [0, 1/2] and y = nµ/2 ∈ [0, 1], so that µ = 2y/n. We may express all quantities using x and y. The difficulty here is that we have two variables. We first check the dependency in terms of y of the expressions B 3 , B 4 , and B 6 + B 7 . These are univariate polynomials with an affine dependence on y, so their positivity can be checked by checking positivity for y = 1 or y = 0. Again making use of symbolic computation, we can deduce that B 3 is non-negative and decreasing in y, (Constraint 3)

B 4 is non-negative and decreasing in y, (Constraint 4)

B 6 + B 7 is non-negative and increasing in y,

We include below the additional commands to verify these properties:

Given the monotonicity of the bounds in B 3 and B 4 , we only need to check our results for the smaller values of B 3 and B 4 , i.e., for y = 1. Similarly, because of the monotonicity in the term (B 6 + B 7 ) 2 , we may replace (B 6 + B 7 ) 2 by (B 6 + B 7 ) times its upper-bound (i.e., its value at y = 1). Also note that B 5 is divisble by y, and we have B 3 > 0 and B 4 > 0 for n > 1. Using this, we can show that Constraints 5-9 are satisfied by checking the positivity of the following polynomials:

Constraints 5-7 are affine in y and we thus only need to check positivity for y = 0 and y = 1.

Summing this up to iteration k yields

Finally, by Jensen's inequality note that

B.8 Intial values of Lyapunov function

To obtain the results of Theorem 1, all that remains is computing the initial value of the Lyapunov function for the two initializations.

B.8.1 Initialization with the zero vector

If we initialize with y 0 = 0 we have

Plugging in the values of our parameters,

and using σ 2 = 1 n i f i (x * ) 2 we obtain (noting that e f (x * ) = 0)

B.8.2 Initialization with average gradient

If we initialize with y 0 i = y 0 i = f i (x 0 ) -(1/n) i f i (x 0 ) we have, by noting that we still have (y 0 ) e = 0,

By observing that y 0 = f (x 0 ) -eg (x 0 ) and using [Nesterov, 2004, Equations 2.17], we can bound the norm in the second term as follows:

y 0 -f (x * ) 2 = (f (x 0 ) -f (x * )) -e(g (x 0 ) -g (x * ) 2 2 f (x 0 ) -f (x * ) 2 + 2 e(g (x 0 ) -g (x * )) 2 = 2 n i=1 f i (x 0 ) -f i (x * ) 2 + 2n g (x 0 ) -g (x * ) 2

4L

n i=1

[f i (x 0 ) -f i (x * ) -f i (x * ) (x 0 -x * )] + 4nL(g(x 0 ) -g(x * ))

= 8nL(g(x 0 ) -g(x * )).

Using this in the Lyapunov function we obtain L(θ 0 ) 3 2 (g(x 0 ) -g(x * )) + 4L n x 0 -x * 2 .