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

Mark Schmidt
email: mark.schmidt@inria.fr

Francis Bach
email: francis.bach@ens.fr

A Stochastic Gradient Method with an Exponential Convergence Rate for Finite Training Sets

We propose a new stochastic gradient method for optimizing the sum of a finite set of smooth functions, where the sum is strongly convex. While standard stochastic gradient methods converge at sublinear rates for this problem, the proposed method incorporates a memory of previous gradient values in order to achieve a linear convergence rate. In a machine learning context, numerical experiments indicate that the new algorithm can dramatically outperform standard algorithms, both in terms of optimizing the training error and reducing the test error quickly.

Introduction

A plethora of the problems arising in machine learning involve computing an approximate minimizer of the sum of a loss function over a large number of training examples, where there is a large amount of redundancy between examples. The most wildly successful class of algorithms for taking advantage of this type of problem structure are stochastic gradient (SG) methods [START_REF] Robbins | A stochastic approximation method[END_REF], [START_REF] Bottou | Large scale online learning[END_REF]. Although the theory behind SG methods allows them to be applied more generally, in the context of machine learning SG methods are typically used to solve the problem of optimizing a sample average over a finite training set, i.e., minimize

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

As an example, in the case of 2 -regularized logistic regression we have f i (x) := λ 2 x 2 + log(1 + exp(-b i a T i x)), where a i ∈ R p and b i ∈ {-1, 1} are the training examples associated with a binary classification problem and λ is a regularization parameter. More generally, any 2 -regularized empirical risk minimization problem of the form minimize

x∈R p λ 2 x 2 + 1 n n i=1 l i (x), (2)
falls in the framework of (1) provided that the loss functions l i are convex and smooth. An extensive list of convex loss functions used in machine learning is given by [START_REF] Teo | A scalable modular convex solver for regularized risk minimization[END_REF], and we can even include non-smooth loss functions (or regularizers) by using smooth approximations.

The standard 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)
Using x * to denote the unique minimizer of g, the FG method with a constant step size achieves a linear convergence rate:

g(x k ) -g(x * ) = O(ρ k ),
for some ρ < 1 which depends on the condition number of g [Nesterov, 2004, Theorem 2.1.15]. Linear convergence is also known as geometric or exponential convergence, because the cost is cut by a fixed fraction on each iteration. Despite the fast convergence rate of the FG method, it can be unappealing when n is large because its iteration cost scales linearly in n. SG methods, on the other hand, have an iteration cost which is independent of n, making them suited for that setting.

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 α k is a step-size and a training example i k is selected uniformly among 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 that, for a suitably chosen decreasing step-size sequence {α k }, the SG iterations achieve the sublinear convergence rate

E[g(x k )] -g(x * ) = O(1/k),
where the expectation is taken with respect to the selection of the i k variables. Under certain assumptions this convergence rate is optimal for strongly-convex optimization 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 cannot hope to obtain a better convergence rate 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 an exponential converengence rate while preserving the iteration cost of SG methods.

The primay 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 [START_REF] Blatt | A convergent incremental gradient method with a constant step size[END_REF], which combines the low iteration cost of SG methods with a linear convergence rate as in FG methods. The SAG method uses iterations of the form

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

y k i = f i (x k ) if i = i k , y k-1 i otherwise.
That is, like the FG method, the step incorporates a gradient with respect to each training example. But, like the SG method, each iteration only computes the gradient with respect to a single training example and the cost of the iterations is independent of n. Despite the low cost of the SAG iterations, in this paper we show that the SAG iterations have a linear convergence rate, 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 training example 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 also show that for certain problems the convergence rate of SAG is faster than is possible for standard FG methods.

In a machine learning context where g(x) is a training cost associated with a predictor parameterized by x, we are often ultimately interested in the testing cost, the expected loss on unseen data points.

Note that a linear convergence rate for the training cost does not translate into a similar rate for the testing cost, and an appealing propertly of SG methods is that they achieve the optimal O(1/k) rate for the testing cost as long as every datapoint is seen only once. However, as is common in machine learning, we assume that we are only given a finite training data set and thus that datapoints are revisited multiple times. In this context, the analysis of SG methods only applies to the training cost and, although our analysis also focuses on the training cost, in our experiments the SAG method typically reached the optimal testing cost faster than both FG and SG methods.

The next section reviews 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, most of the applications focusing on finite datasets, we are not aware of any other SG method that achieves a linear convergence rate while preserving the iteration cost of standard SG methods. Section 3 states the (standard) assumptions underlying our analysis and gives the main technical results; we first give a slow linear convergence rate that applies for any problem, and then give a very fast linear convergence rate that applies when n is sufficiently large. Section 4 discusses practical implementation issues, including how to reduce the storage cost from O(np) to O(n) when each f i only depends on a linear combination of x. Section 5 presents a numerical comparison of an implementation based on SAG to 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.

Related Work

There is 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 ), (6)
where the selection function S(j, i 1:k ) is equal to 1/n if j corresponds to the last iteration where j = i k 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 [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF], and while this averaging procedure leads 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 O(1/k) rate.

Iterate Averaging: Rather than averaging the gradients, some authors use the basic SG iteration but take an average over x k values. 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]. Baher's method [Kushner and Yin, 2003, §1.3.4] combines gradient averaging with online iterate averaging, and also displays appealing asymptotic properties. The epoch SG method uses averaging to obtain the O(1/k) rate even for non-smooth objectives [START_REF] Hazan | Beyond the regret minimization barrier: an optimal algorithm for stochastic strongly-convex optimization[END_REF]. However, the convergence rates of these averaging methods remain sublinear.

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 [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 non-linear conjugate gradient, quasi-Newton, and Hessian-free Newton methods. Several authors have analyzed stochastic variants of these algorithms [START_REF] Schraudolph | Local gain adaptation in stochastic gradient descent[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]. Under certain conditions these variants are convergent with an O(1/k) rate [START_REF] Sunehag | Variable metric stochastic approximation theory[END_REF]. Alternately, if we split the convergence rate into a deterministic and stochastic part, these methods can improve the dependency on the deterministic part [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, as with all other methods we have discussed thus far in this section, we are not aware of any existing method of this flavor that improves on the O(1/k) rate.

Constant step size: If the SG iterations are used with a constant step size (rather than a decreasing sequence), then the convergence rate of the method can be split into two parts [Nedic and Bertsekas, 2000, Proposition 2.4], where the first part depends on k and converges linearly to 0 and the second part is independent of k but 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, convergence of the basic SG method with a constant step size has only been shown 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 [START_REF] Kesten | Accelerated stochastic approximation[END_REF], [START_REF] Delyon | Accelerated stochastic approximation[END_REF]. This leads to convergence of the method and allows it to potentially achieve periods of linear convergence where the step size stays constant. However, the overall convergence rate of the method remains sublinear.

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 linear convergence rate. Bertsekas 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 [START_REF] Bertsekas | A new class of incremental gradient methods for least squares problems[END_REF]. However, the weighting is numerically unstable and the linear convergence rate treats full passes through the data as iterations. A related strategy is to group the f i functions into 'batches' of increasing size and perform SG iterations on the batches [START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF]. In both cases, the iterations that achieve the linear rate have a cost that is not independent of n, as opposed to SAG.

Incremental Aggregated Gradient: Finally, Blatt et al. presents the most closely-related algorithm, the IAG method [START_REF] Blatt | A convergent incremental gradient method with a constant step size[END_REF]. 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 proof technique involving bounding the gradients and iterates simultaneously by a Lyapunov potential function, in this work we give an explicit linear convergence rate for general strongly-convex functions using the SAG iterations that only examine a single training example. Further, as our analysis and experiments show, when the number of training examples is sufficiently large, the SAG iterations achieve a linear convergence rate under a much larger set of step sizes than 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 improved practical performance. We also emphasize that in our experiments IAG and the basic FG method perform similarly, while SAG performs much better, showing that the simple change (random selection vs. cycling) can dramatically improve optimization performance.

Convergence Analysis

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

f i (x) -f i (y) ≤ L x -y .
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. In addition, we also assume that 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. This is a stronger assumption and is not satisfied by all machine learning models. However, note that in machine learning we are typically free to choose the regularizer, and we can always add an 2 -regularization term as in Eq. ( 2) to transform any convex problem into a strongly-convex problem (in this case we have µ ≥ λ). Note that strong-convexity implies that the problem is solvable, meaning that there exists some unique x * that achieves the optimal function value. Our convergence results assume that we initialize y 0 i to a zero vector for all i, and our results depend on the variance of the gradient norms at the optimum x * , denoted by

σ 2 = 1 n i f i (x * ) 2 .
Finally, all our convergence results consider expectations with respect to the internal randomization of the algorithm, and not with respect to the data (which are assumed to be deterministic and fixed).

We first consider the convergence rate of the method when using a constant step size of α k = 1 2nL , which is similar to the step size needed for convergence of the IAG method in practice.

Proposition 1 With a constant step size of α k = 1 2nL , the SAG iterations satisfy for k ≥ 1:

E x k -x * 2 1 - µ 8Ln k 3 x 0 -x * 2 + 9σ 2 4L 2 .
The proof is given in the Appendix. Note that the SAG iterations also trivially obtain the O(1/k) rate achieved by SG methods, since

1 - µ 8Ln k exp - kµ 8Ln 8Ln kµ = O(n/k),
albeit with a constant which is proportional to n. Despite this constant, they are advantageous over SG methods in later iterations because they obtain an exponential convergence rate as in FG methods. We also note that an exponential convergence rate is obtained for any constant step size smaller than1 2nL . In terms of passes through the data, the rate in Proposition 1 is similar to that achieved by the basic FG method. However, our next result shows that, if the number of training examples is slightly larger than L/µ (which will often be the case, as discussed in Section 6), then the SAG iterations can use a larger step size and obtain a better convergence rate that is independent of µ and L (see proof in the Appendix).

Proposition 2 If n 8L
µ , with a step size of α k = 1 2nµ the SAG iterations satisfy for k n:

E g(x k ) -g(x * ) C 1 - 1 8n k , with C = 16L 3n x 0 -x * 2 + 4σ 2 3nµ 8 log 1 + µn 4L + 1 .
We state this result for k n because we assume that the first n iterations of the algorithm use an SG method and that we initialize the subsequent SAG iterations with the average of the iterates, which leads to an O((log n)/k) rate. In contrast, using the SAG iterations from the beginning gives the same rate but with a constant proportional to n. Note that this bound is obtained when initializing all y i to zero after the SG phase. 1 However, in our experiments we do not use the SG initialization but rather use a minor variant of SAG (discussed in the next section), which appears more difficult to analyze but which gives better performance.

It is interesting to compare this convergence rate with the known convergence rates of first-order methods [Nesterov, 2004, see §2]. For example, if we take n = 100000, L = 100, and µ = 0.01 then the basic FG method has a rate of ((L -µ)/(L + µ))2 = 0.9996 and the 'optimal' AFG method has a faster rate of (1 -µ/L) = 0.9900. In contrast, running n iterations of SAG has a much faster rate of (1 -1/8n) n = 0.8825 using the same number of evaluations of f i . Further, the lower-bound for a black-box first-order method is ((

√ L - √ µ)/( √ L + √ µ)) 2 = 0.9608, indicating that SAG can
be substantially faster than any FG method that does not use the structure of the problem. 2 In the Appendix, we compare Propositions 1 and 2 to the rates of primal and dual FG and coordinate-wise methods for the special case of 2 -regularized leasts squares.

Even though n appears in the convergence rate, 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, each pass through the data reduces the excess cost by a constant multiplicative factor that is independent of the problem, as long as n 8L/µ. Further, while the step size in Proposition 2 depends on µ and n, we can obtain the same convergence rate by using a step size as large as α k = 1 16L . This is because the proposition is true for all values of µ satisfying µ L 8 n , so we can choose the smallest possible value of µ = 8L n . We have observed in practice that the IAG method with a step size of α k = 1 2nµ may diverge, even under these assumptions. Thus, for certain problems the SAG iterations can tolerate a much larger step size, 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.

While we have stated Proposition 1 in terms of the iterates and Proposition 2 in terms of the function values, the rates obtained on iterates and function values are equivalent because, by the Lipschitz and strong-convexity assumptions, we have µ

2 x k -x * 2 g(x k ) -g(x * ) L 2 x k -x * 2 .

Implementation Details

In this section we describe modifications that substantially reduce the SAG iteration's memory requirements, as well as modifications that lead to better practical performance.

Structured gradients: For many problems the storage cost of O(np) for the y k i vectors is prohibitive, but we can often use structure in the f i to reduce this cost. For example, many loss functions f i take the form f i (a T i x) for a vector a i . Since 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 T i k x k rather than the full gradient a T i f i (u k i )
, reducing the storage cost to O(n). Further, because of the simple form of the SAG updates, if a i is sparse we can use 'lazy updates' in order to reduce the iteration cost from O(p) down to the sparsity level of a i .

Mini-batches: To employ vectorization and parallelism, practical SG implementations often group training examples into 'mini-batches' and perform SG iterations on the mini-batches. We can also use mini-batches within the SAG iterations, and for problems with dense gradients this decreases the storage requirements of the algorithm since we only need a y k i for each mini-batch. Thus, for example, using mini-batches of size 100 leads to a 100-fold reduction in the storage cost.

Step-size re-weighting: On early iterations of the SAG algorithm, when most y k i are set to the uninformative zero vector, rather than dividing α k in (5) by n we found it was more effective to divide by m, the number of unique i k values that we have sampled so far (which converges to n). This modification appears more difficult to analyze, but with this modification we found that the SAG algorithm outperformed the SG/SAG hybrid algorithm analyzed in Proposition 2.

Exact regularization: For regularized objectives like (2) we can use the exact gradient of the regularizer, rather than approximating it. For example, our experiments on 2 -regularized optimization problems used the recursion

d ← d -y i , y i ← l i (x k ), d ← d + y i , x ← 1 -αλ x - α m d . (7)
This can be implemented efficiently for sparse data sets by using the representation x = κz, where κ is a scalar and z is a vector, since the update based on the regularizer simply updates κ.

Large step sizes: Proposition 1 requires α k 1/2Ln while under an additional assumption Proposition 2 allows α k 1/16L. In practice we observed better performance using step sizes of α k = 1/L and α k = 2/(L + nµ). These step sizes seem to work even when the additional assumption of Proposition 2 is not satisfied, and we conjecture that the convergence rates under these step sizes are much faster than the rate obtained in Proposition 1 for the general case.

Line search: Since L is generally not known, we experimented with a basic line-search, where we start with an initial estimate L 0 , and we double this estimate whenever we do not satisfy the instantiated Lipschitz 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 .
To avoid instability caused by comparing very small numbers, we only do this test when f i k (x k ) 2 > 10 -8 . To allow the algorithm to potentially achieve a faster rate due to a higher degree of local smoothness, we multiply L k by 2 (-1/n) after each iteration.

Experimental Results

Our experiments compared an extensive variety of competitive FG and SG methods. Our first experiments focus on the following methods, which we chose because they have no dataset-dependent tuning parameters:

-Steepest: The full gradient method described by iteration (3), with a line-search that uses cubic Hermite polynomial interpolation to find a step size satisfying the strong Wolfe conditions, and where the parameters of the line-search were tuned for the problems at hand.

-AFG: Nesterov's accelerated full gradient method [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 fixed step size are interleaved with an extrapolation step, and we use an adaptive line-search based on [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.3 This method is by far the most complicated method we considered.

-Pegasos: The state-of-the-art SG method described by iteration (4) with a step size of α k = 1/µk and a projection step onto a norm-ball known to contain the optimal solution Shalev-Shwartz et al. [2007].

-RDA: The regularized dual averaging method [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF], another recent state-of-the-art SG method.

-ESG: The epoch SG method [START_REF] Hazan | Beyond the regret minimization barrier: an optimal algorithm for stochastic strongly-convex optimization[END_REF], which runs SG with a constant step size and averaging in a series of epochs, and is optimal for non-smooth stochastic strongly-convex optimization.

-NOSG: The nearly-optimal SG method Ghadimi and Lan [2010], which combines ideas from SG and AFG methods to obtain a nearly-optimal dependency on a variety of problem-dependent constants.

-SAG: The proposed stochastic average gradient method described by iteration (5) using the modifications discussed in the previous section. We used a step-size of α k = 2/(L k + nλ) where L k is either set constant to the global Lipschitz constant (SAG-C) or set by adaptively estimating the constant with respect to the logistic loss function using the line-search described in the previous section (SAG-LS). The SAG-LS method was initialized with L 0 = 1 .

The theoretical convergence rates suggest the following strategies for deciding on whether to use an FG or an SG method:

1. If we can only afford one pass through the data, then an SG method should be used.

2. 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 L-BFGS. To test whether this is indeed the case on real data sets, we performed experiments on a set of freely available benchmark binary classification data sets. The quantum (p = 50000, p = 78) and protein (n = 145751, p = 74) data set was obtained from the KDD Cup 2004 website,4 the sido data set was obtained from the Causality Workbench website,5 while the rcv1 (n = 20242, p = 47236) and covertype (n = 581012, p = 54) data sets were obtained from the LIBSVM data website. 6 Although our method can be applied to any differentiable function, on these data sets we focus on the 2 -regularized logistic regression problem, with λ = 1/n. We split each dataset in two, training on one half and testing on the other half. 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. We measure the training and testing costs as a function of the number of effective passes through the data, measured as the number of f i evaluations divided by n. These results are thus independent of the practical implementation of the algorithms. We plot the training and testing costs of the different methods for 30 effective passes through the data in Figure 1.

In our second series of experiments, we sought to test whether SG methods (or the IAG method) with a very carefully chosen step size would be competitive with the SAG iterations. In particular, we compared the following variety of basic FG and SG methods.

1. FG: The full gradient method described by iteration (3).

AFG:

The accelerated full gradient method [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) are interleaved with an extrapolation step.

peg:

The pegasos algorithm of Shalev-Shwartz et al. [2007], but where we multiply the step size by a constant.

SG:

The stochastic gradient method described by iteration (4), where we use a constant stepsize.

ASG:

The stochastic gradient method described by iteration (4), where we use a constant step size and average the iterates.7

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) but with a cyclic choice of i k .

7. SAG: The proposed stochastic average gradient method described by iteration (5).

For all of the above methods, we chose the step size that gave the best performance among powers of 10. On the full data sets, we compare these methods to each other and to the L-BFGS and the SAG-LS algorithms from the previous experiment in Figure 2, which also shows the selected step sizes.

We can observe several trends across these experiments:

-FG vs. SG: Although the performance of SG methods can be catastrophic if the step size is not chosen carefully (e.g., the quantum and covertype data), with a carefully-chosen step-size the SG methods always do substantially better than FG methods on the first few passes through the data. In contrast, the adaptive FG methods in the first experiment are not sensitive to the step size and because of its steady progress the best FG method (L-BFGS) always eventually 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, but continue to make steady (linear) progress which leads to better performance than SG methods. The significant speed-up observed for SAG in reaching low training costs often also seems to translate into reaching the optimal testing cost more quickly than the other methods. We also note that the proposed line-search seems to perform as well or better than choosing the optimal fixed step-size in hind sight.

-IAG vs. SAG: The second experiment shows that the IAG method performs similarly to the regular FG method, and they also show the surprising result that the randomized SAG method outperforms the closely-related deterministic IAG method by a very large margin. This is due to the larger step sizes used by the SAG iterations, which would cause the IAG iterations to diverge.

Discussion

Optimal regularization strength: One might wonder if the additional hypothesis in Proposition 2 is satisfied in practice. In a learning context, where each function f i is the loss associated to a single data point, L is equal to the largest value of the loss second derivative ξ (1 for the square loss, 1/4 for the logistic loss) times R 2 , where R is a the uniform bound on the norm of each data point. Thus, the constraint µ L 8 n is satisfied when λ 8ξR 2 n . In low-dimensional settings, the optimal regularization parameter is of the form C/n [START_REF] Liang | Asymptotically optimal regularization in smooth parametric models[END_REF] where C is a scalar constant, and may thus violate the constraint. However, the improvement with respect to regularization parameters of the form λ = C/ √ n is known to be asymptotically negligible, and in any case in such low-dimensional settings, regular stochastic or batch gradient descent may be efficient enough in practice. In the more interesting high-dimensional settings where the dimension p of our covariates is not small compared to the sample size n, then all theoretical analyses we are aware of advocate settings of λ which satisfy this constraint. For example, [START_REF] Sridharan | Fast rates for regularized objectives[END_REF] considers parameters of the form λ = C √ n in the parametric setting, while [START_REF] Eberts | Optimal learning rates for least squares SVMs using Gaussian kernels[END_REF] considers λ = C n β with β < 1 in a non-parametric setting.

Training cost vs. testing cost: The theoretical contribution of this work is limited to the convergence rate of the training cost. Though there are several settings where this is the metric of interest (e.g., variational inference in graphical models), in many cases one will be interested in the convergence speed of the testing cost. Since the O(1/k) convergence rate of the testing cost, achieved by SG methods with decreasing step sizes (and a single pass through the data), is provably optimal when the algorithm only accesses the function through unbiased measurements of the objective and its gradient, it is unlikely that one can obtain a linear convergence rate for the testing cost with the SAG iterations. However, as shown in our experiments, the testing cost of the SAG iterates often reaches its minimum quicker than existing SG methods, and we could expect to improve the constant in the O(1/k) convergence rate, as is the case with online second-order methods [START_REF] Bottou | The tradeoffs of large scale learning[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 algorithm, and (iii) choosing the step size while running the algorithm. This paper showed that the SAG iterations achieve a much faster convergence rate, but the SAG iterations may also be advantageous in terms of tuning step sizes and designing termination criteria. In particular, the SAG iterations suggest a natural termination criterion; since the average of the y k i variables converges to g (x k ) as x k -x k-1 converges to zero, we can use (1/n) i y k i as an approximation of the optimality of x k . Further, while SG methods require specifying a sequence of step sizes and mispecifying this sequence can have a disastrous effect on the convergence rate [Nemirovski et al., 2009, §2.1], our theory shows that the SAG iterations iterations achieve a linear convergence rate for any sufficiently small constant step size and our experiments indicate that a simple line-search gives strong performance.

Appendix

In this Appendix, we first give the proofs of the two propositions. Subsequently, we compare the convergence rates of primal and dual FG and coordinate-wise methods to the rates of SAG for 2 -regularized least squares in terms of effective passes through the data.

A Proofs of the propositions

We present here the proofs of Propositions 1 and 2.

A.1 Problem set-up and notations

We use g = 1 n n i=1 f i to denote a µ-strongly convex function, where the functions f i , i = 1, . . . , n are convex functions from R p to R with L-Lipschitz continuous gradients. Let us denote by x * the unique minimizer of g.

For k 1, the stochastic average gradient algorithm performs the recursion

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

y k i = f i (x k-1 ) if i = i k , y k-1 i otherwise.
Denoting z k i a random variable which takes the value 1 -1 n with probability 1 n and -1 n otherwise (thus with zero expectation), this is equivalent to

y k i = 1 - 1 n y k-1 i + 1 n f i (x k-1 ) + z k i f i (x k-1 ) -y k-1 i x k = x k-1 - α n n i=1 1 - 1 n y k-1 i + 1 n f i (x k-1 ) + z k i f i (x k-1 ) -y k-1 i = x k-1 - α n 1 - 1 n e y k-1 + g (x k-1 ) + (z k ) f (x k-1 ) -y k-1 , with e =    I . . . I    ∈ R np×p , f (x) =    f 1 (x) . . . f n (x)    ∈ R np , z k =    z k 1 I . . . z k n I    ∈ R np×p .
Using this definition of z k , we have E[(z k )(z k ) ] = 1 n I -1 n 2 ee . Note that, for a given k, the variables z k 1 , . . . , z k n are not independent. We also use the notation

θ k =      y k 1 . . . y k n x k      ∈ R (n+1)p , θ * =      f 1 (x * ) . . . f n (x * ) x *      ∈ R (n+1)p .
Finally, if M is a tp × tp matrix and m is a tp × p matrix, then:

• 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.

A.2 Outline of the proofs

Each Proposition will be proved in multiple steps.

1. We shall find a Lyapunov function Q from R (n+1)p to R such that the sequence EQ(θ k ) decreases at a linear rate.

2. We shall prove that Q(θ k ) dominates x k -x * 2 (in the case of Proposition 2) or g(x k ) -g(x * ) (in the case of Proposition 2) by a constant for all k.

3. In the case of Proposition 2, we show how using one pass of stochastic gradient as the initialization provides the desired result.

Throughout the proofs, F k will denote the σ-field of information up to (and including time k), i.e., F k is the σ-field generated by z 1 , . . . , z k .

A.3 Convergence results for stochastic gradient descent

The constant in both our bounds depends on the initialization chosen. While this does not affect the linear convergence of the algorithm, the bound we obtain for the first few passes through the data is the O(1/k) rate one would get using stochastic gradient descent, but with a constant proportional to n. This problem can be alleviated for the second bound by running stochastic gradient descent for a few iterations before running the SAG algorithm. In this section, we provide bounds for the stochastic gradient descent algorithm which will prove useful for the SAG algorithm.

The assumptions made in this section about the functions f i and the function g are the same as the ones used for SAG. To get initial values for x 0 and y 0 , we will do one pass of standard stochastic gradient.

We denote by σ 2 = 1 n n i=1 f i (x * ) 2 the variance of the gradients at the optimum. We will use the following recursion:

xk = xk-1 -γ k f i k xk-1 . Denoting δ k = E xk -x * 2
, we have (following [START_REF] Bach | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF])

δ k δ k-1 -2γ k (1 -γ k L)E g (x k-1 ) (x k-1 -x * ) + 2γ 2 k σ 2 . Indeed, we have xk -x * 2 = xk-1 -x * 2 -2γ k f i k (x k-1 ) (x k-1 -x * ) + γ 2 k f i k (x k-1 ) 2 xk-1 -x * 2 -2γ k f i k (x k-1 ) (x k-1 -x * ) + 2γ 2 k f i k (x * ) 2 + 2γ 2 k f i k (x k-1 ) -f i k (x * ) 2 xk-1 -x * 2 -2γ k f i k (x k-1 ) (x k-1 -x * ) + 2γ 2 k f i k (x * ) 2 + 2Lγ 2 k (f i k (x k-1 ) -f i k (x * )) (x k-1 -x * ) .
By taking expectations, we get

E xk -x * 2 |F k-1 xk-1 -x * 2 -2γ k g (x k-1 ) (x k-1 -x * ) + 2γ 2 k σ 2 + 2Lγ 2 k g (x k-1 ) (x k-1 -x * ) E xk -x * 2 E xk-1 -x * 2 -2γ k (1 -γ k L)E g (x k-1 ) (x k-1 -x * ) + 2γ 2 k σ 2
Thus, if we take

γ k = 1 2L + µ 2 k , we have γ k 2γ k (1 -γ k L) and δ k δ k-1 -γ k E g (x k-1 ) (x k-1 -x * ) + 2γ 2 k σ 2 δ k-1 -γ k E g(x k-1 ) -g(x * ) + µ 2 δ k-1 + 2γ 2 k σ 2 using the strong convexity of g Eg(x k-1 ) -g(x * ) - 1 γ k δ k + 1 γ k - µ 2 δ k-1 + 2γ k σ 2 -2L + µ 2 k δ k + 2L + µ 2 (k -1) δ k-1 + 2γ k σ 2 .
Averaging from i = 0 to k -1 and using the convexity of g, we have

1 k k-1 i=0 Eg(x k-1 ) -g(x * ) 2L k δ 0 + 2σ 2 k k i=1 γ i Eg 1 k k-1 i=0 x i -g(x * ) 2L k δ 0 + 2σ 2 k k i=1 γ i 2L k x 0 -x * 2 + 2σ 2 k k i=1 1 2L + µ 2 i 2L k L x 0 -x * 2 + 2σ 2 k k 0 1 2L + µ 2 t dt 2L k x 0 -x * 2 + 4σ 2 kµ log 1 + µk 4L .

A.4 Important lemma

In both proofs, our Lyapunov function contains a quadratic term

R(θ k ) = (θ k -θ * ) A b b c (θ k -θ * )
for some values of A, b and c. The lemma below computes the value of R(θ k ) in terms of 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 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 .
Note that for square n × n matrix, diag(M ) denotes 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.

Proof Throughout the proof, we will use the equality g (x) = e f (x)/n. Moreover, all conditional expectations of linear functions of z k will be equal to zero.

We have

E (θ k -θ * ) A b b c (θ k -θ * ) F k-1 = E (y k -f (x * )) A(y k -f (x * )) + 2(y k -f (x * )) b(x k -x * ) + (x k -x * ) c(x k -x * )|F k-1 . (8)
The first term (within the expectation) on the right-hand side of Eq. ( 8) is equal to

(y k -f (x * )) A(y k -f (x * )) = 1 - 1 n 2 (y k-1 -f (x * )) A(y k-1 -f (x * )) + 1 n 2 (f (x k-1 ) -f (x * )) A(f (x k-1 ) -f (x * )) + [Diag(z k )(f (x k-1 ) -y k-1 )] A[Diag(z k )(f (x k-1 ) -y k-1 )] + 2 n 1 - 1 n (y k-1 -f (x * )) A(f (x k-1 ) -f (x * )) .
The only random term (given F k-1 ) is the third one whose expectation is equal to

E [Diag(z k )(f (x k-1 ) -y k-1 )] A[Diag(z k )(f (x k-1 ) -y k-1 )]|F k-1 = 1 n (f (x k-1 ) -y k-1 ) Diag(diag(A)) - 1 n A (f (x k-1 ) -y k-1 ) .
A.5 Analysis for α = 1 2nL

We now prove Proposition 1, providing a bound for the convergence rate of the SAG algorithm in the case of a small step size, α = 1 2nL . Proof

Step 1 -Linear convergence of the Lyapunov function In this case, our Lyapunov function is quadratic, i.e.,

Q(θ k ) = (θ k -θ * ) A b b c (θ k -θ * ) .
We consider

A = 3nα 2 I + α 2 n ( 1 n -2)ee b = -α(1 - 1 n )e c = I S = 3nα 2 I b - α n ec = -αe .
The goal will be to prove that

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1
) is negative for some δ > 0. This will be achieved by bounding all the terms by a term depending on g (x k-1 ) (x k-1 -x * ) whose positivity is guaranteed by the convexity of g.

We have, with our definition of A, b and c:

S -Diag(diag(S)) = 3nα 2 I -3nα 2 I = 0 e (f (x k-1 ) -f (x * )) = n[g (x k-1 ) -g (x * )] = ng (x k-1 ) .
The difference is then:

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1 ) (y k-1 -f (x * )) 3nα 2 δ - 1 n I + (1 -δ) α 2 n 2 - 1 n ee (y k-1 -f (x * )) + δ(x k-1 -x * ) (x k-1 -x * ) -(2α -3α 2 nL)(x k-1 -x * ) g (x k-1 ) -2αδ 1 - 1 n (y k-1 -f (x * )) e(x k-1 -x * ).
Note that for any symmetric negative definite matrix M and for any vectors s and t we have

(s + 1 2 M -1 t) M (s + 1 2 M -1 t) 0,
and thus that s M s + s t -

1 4 t M -1 t .
Using this fact with

M = 3nα 2 δ - 1 n I + (1 -δ) α 2 n 2 - 1 n ee = 3nα 2 δ - 1 n I - ee n + α 2 3nδ -1 -2δ + δ -1 n ee n s = y k-1 -f (x * ) t = -2αδ 1 - 1 n e(x k-1 -x * ) ,
we have

(y k-1 -f (x * )) 3nα 2 δ - 1 n I + (1 -δ) α 2 n 2 - 1 n ee (y k-1 -f (x * )) -2αδ 1 - 1 n (y k-1 -f (x * )) e(x k-1 -x * ) -α 2 δ 2 1 - 1 n 2 (x k-1 -x * ) e 3nα 2 δ - 1 n I - ee n +α 2 3nδ -1 -2δ + δ -1 n ee n -1 e(x k-1 -x * ) = - α 2 δ 2 1 -1 n 2 n α 2 3nδ -1 -2δ + δ-1 n x k-1 -x * 2 = - δ 2 1 -1 n 2 n 3nδ -1 -2δ + δ-1 n x k-1 -x * 2 .
A sufficient condition for M to be negative definite is to have δ 1 3n . The bound then becomes

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1 ) -(2α -3α 2 nL)(x k-1 -x * ) g (x k-1 ) + δ - δ 2 1 -1 n 2 3nδ -1 -2δ + δ-1 n n x k-1 -x * 2 .
We shall use the Schur complement condition for positive definiteness. Since A is positive definite, the other condition to verify is 2

3 I -b A -1 b 0. 2 3 I -α 2 1 - 1 n 2 e 3nα 2 + α 2 n -2α 2 ee n -1 e = 2 3 I - n 1 -1 n 2 3n + 1 n -2 ee n 2 3 I - n 3n -2 ee n 0 for n 2 ,
and so P dominates 0 0 0 1 3 I .

This yields

E x k -x * 2 3EQ(θ k ) 3 1 - µ 8nL k Q(θ 0 ) .
We have

Q(θ 0 ) = 3nα 2 i y 0 i -f i (x * ) 2 + (1 -2n)α n 2 i y 0 i 2 -2α 1 - 1 n (x 0 -x * ) i y 0 i + x 0 -x * 2 = 3 4nL 2 i y 0 i -f i (x * ) 2 + (1 -2n) 2n 3 L i y 0 i 2 - n -1 n 2 L (x 0 -x * ) i y 0 i + x 0 -x * 2 .
Initializing all the y 0 i to 0, we get

Q(θ 0 ) = 3σ 2 4L 2 + x 0 -x * 2 ,
and

E x k -x * 2 1 - µ 8nL k 9σ 2 4L 2 + 3 x 0 -x * 2 .
A.6 Analysis for α = 1 2nµ

Step 1 -Linear convergence of the Lyapunov function

We now prove Proposition 2, providing a bound for the convergence rate of the SAG algorithm in the case of a small step size, α = 1 2nµ . We shall use the following Lyapunov function:

Q(θ k ) = 2g x k + α n e y k -2g(x * ) + (θ k -θ * ) A b b c (θ k -θ * ) , with A = ηα n I + α n (1 -2ν)ee b = -νe c = 0 .
This yields

S = ηα n I + α n ee Diag(diag(S)) = (1 + η)α n I S -Diag(diag(S)) = α n (ee -I) 1 - 2 n S + 1 n Diag(diag(S)) = 1 - 2 n ηα n I + α n ee + 1 n (1 + η)α n I = 1 - 2 n α n ee + η - η -1 n α n I .
We have

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1 ) = 2g(x k-1 ) -2g(x * ) -2(1 -δ)g x k-1 + α n e y k-1 + 2(1 -δ)g(x * ) + (y k-1 -f (x * )) 1 - 2 n α n ee + η - η -1 n α n I -(1 -δ) ηα n I -(1 -δ) α n (1 -2ν)ee (y k-1 -f (x * )) - 2ν n (x k-1 -x * ) e (f (x k-1 ) -f (x * )) + (1 + η)α n 2 (f (x k-1 ) -f (x * )) (f (x k-1 ) -f (x * )) + 2α n 2 (y k-1 -f (x * )) ee -I (f (x k-1 ) -f (x * )) + 2 1 n -δ ν(y k-1 -f (x * )) e(x k-1 -x * ).
Our goal will now be to express all the quantities in terms of (x k-1 -x * ) g (x k-1 ) whose positivity is guaranteed by the convexity of g.

Using the convexity of g, we have

-2(1 -δ)g x k-1 + α n e y k-1 -2(1 -δ) g(x k-1 ) + α n g (x k-1 )e y k-1 .
Using the Lipschitz property of the gradients of f i , we have

(f (x k-1 ) -f (x * )) (f (x k-1 ) -f (x * )) = n i=1 f i (x k-1 ) -f i (x * ) 2 n i=1 L(f i (x k-1 ) -f i (x * )) (x k-1 -x * ) = nL(g (x k-1 ) -g (x * )) (x k-1 -x * ) . Using e [f (x k-1 ) -f (x * )] = ng (x k-1 ), we have - 2ν n (x k-1 -x * ) e (f (x k-1 ) -f (x * )) = -2ν(x k-1 -x * ) g (x k-1 ) 2α n 2 (y k-1 -f (x * )) ee (f (x k-1 ) -f (x * )) = 2α n (y k-1 -f (x * )) eg (x k-1 ) .
Reassembling all the terms together, we get

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1 ) 2δ[g(x k-1 ) -g(x * )] + 2δα n g (x k-1 )e y k-1 + (y k-1 -f (x * )) 1 - 2 n α n ee + η - η -1 n α n I -(1 -δ) ηα n I- (1 -δ) α n (1 -2ν)ee (y k-1 -f (x * )) -2ν - (1 + η)αL n (x k-1 -x * ) g (x k-1 ) - 2α n 2 (y k-1 -f (x * )) f (x k-1 ) -f (x * )) + 2 1 n -δ ν(y k-1 -f (x * )) e(x k-1 -x * ).
Using the convexity of g gives 2δ

[g(x k-1 ) -g(x * )] 2δ[x k-1 -x * ] g (x k-1 ) ,
and, consequently,

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1 ) 2δ[(x k-1 ) -(x * )] g (x k-1 ) + 2δα n g (x k-1 )e y k-1 + (y k-1 -f (x * )) 1 - 2 n α n ee + η - η -1 n α n I -(1 -δ) ηα n I -(1 -δ) α n (1 -2ν)ee (y k-1 -f (x * )) -2ν - (1 + η)αL n (x k-1 -x * ) g (x k-1 ) - 2α n 2 (y k-1 -f (x * )) f (x k-1 ) -f (x * )) + 2 1 n -δ ν(y k-1 -f (x * )) e(x k-1 -x * ) .
If we regroup all the terms in [(x k-1 )-(x * )] g (x k-1 ) together, and all the terms in (y k-1 -f (x * )) together, we get

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1 ) α n (y k-1 -f (x * )) δη - η -1 n I + δ - 2 n + 2ν(1 -δ) ee (y k-1 -f (x * )) -2ν -2δ - (1 + η)αL n (x k-1 -x * ) g (x k-1 ) + 2(y k-1 -f (x * )) - α n 2 (f (x k-1 ) -f (x * )) + ( 1 n -δ)νe(x k-1 -x * ) + δα n eg (x k-1 ) .
Let us rewrite this as

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1 ) (y k-1 -f (x * )) τ y,I I + τ y,e ee n (y k-1 -f (x * )) + τ x,g (x k-1 -x * ) g (x k-1 ) + (y k-1 -f (x * )) τ y,f (f (x k-1 ) -f (x * )) + τ y,x e(x k-1 -x * ) + τ y,g eg (x k-1 ) with τ y,I = α n δη - η -1 n τ y,e = α δ - 2 n + 2ν(1 -δ) τ x,g = -(2ν -2δ - (1 + η)αL n ) τ y,f = - 2α n 2 τ y,x = 2 1 n -δ ν τ y,g = 2δα n .
Assuming that τ y,I and τ y,e are negative, we have by completing the square that

(y k-1 -f (x * )) τ y,I I + τ y,e ee n (y k-1 -f (x * )) + (y k-1 -f (x * )) τ y,f (f (x k-1 ) -f (x * )) + τ y,x e(x k-1 -x * ) + τ y,g eg (x k-1 ) - 1 4 τ y,f (f (x k-1 ) -f (x * )) + τ y,x e(x k-1 -x * ) + τ y,g eg (x k-1 ) 1 τ y,I I - ee n + 1 τ y,I + τ y,e ee n τ y,f (f (x k-1 ) -f (x * )) + τ y,x e(x k-1 -x * ) + τ y,g eg (x k-1 ) = - 1 4 τ 2 y,f τ y,I f (x k-1 ) -f (x * ) 2 - 1 4 τ 2 y,f n g (x k-1 ) 2 1 τ y,I + τ y,e - 1 τ y,I - 1 4 τ 2 y,x n τ y,I + τ y,e x k-1 -x * 2 - 1 4 τ 2 y,g n τ y,I + τ y,e g (x k-1 ) 2 - 1 2 τ y,f τ y,x n τ y,I + τ y,e (x k-1 -x * ) g (x k-1 ) - 1 2 τ y,f τ y,g n τ y,I + τ y,e g (x k-1 ) 2 - 1 2 τ y,g τ y,x n τ y,I + τ y,e (x k-1 -x * ) g (x k-1 ) ,
where we used the fact that (f (x k-1 ) -f (x * )) e = g (x k-1 ). After reorganization of the terms, we obtain

E[Q(θ k )|F k-1 ] -(1 -δ)Q(θ k-1 ) τ x,g - nτ y,x 2(τ y,I + τ y,e ) (τ y,f + τ y,g ) (x k-1 -x * ) g (x k-1 ) - 1 4 τ 2 y,f n 1 τ y,I + τ y,e - 1 τ y,I + 1 4 τ 2 y,g n τ y,I + τ y,e + 1 2 τ y,f τ y,g n τ y,I + τ y,e g (x k-1 ) 2 - 1 4 τ 2 y,f τ y,I f (x k-1 ) -f (x * ) 2 - 1 4 τ 2 y,x n τ y,I + τ y,e x k-1 -x * 2 .
We thus get

E g(x k ) -g(x * ) 2EQ(θ k ) = 1 - 1 8n k 7 3 (g(x 0 ) -g(x * )) + 7σ 2 6nµ .
Step 3 -Initialization of x 0 using stochastic gradient descent

During the first few iterations, we obtain the O(1/k) rate obtained using stochastic gradient descent, but with a constant which is proportional to n. To circumvent this problem, we will first do n iterations of stochastic gradient descent to initialize x 0 , which will be renamed x n to truly reflect the number of iterations done.

Using the bound from section A.3, we have Eg

1 n n-1 i=0 xi -g(x * ) 2L n x 0 -x * 2 + 4σ 2 nµ log 1 + µn 4L .
And so, using

x n = 1 n n-1 i=0 xi , we have for k n E g(x k ) -g(x * ) 1 - 1 8n k-n 14L 3n x 0 -x * 2 + 28σ 2 3nµ log 1 + µn 4L + 7σ 2 6nµ . Since 1 - 1 8n -n 8 7 ,
we get

E g(x k ) -g(x * ) 1 - 1 8n k 16L 3n x 0 -x * 2 + 32σ 2 3nµ log 1 + µn 4L + 4σ 2 3nµ .

B Comparison of convergence rates

We consider the 2 -regularized least squares 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 SG methods and SAG 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 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 .

B.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 ,
where the speed improvement is determined by the size of 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 σ .

B.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 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.

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 1 -

µ d nL i d n = 1 - n + m σ /λ n(n + M i /λ) n = 1 - nλ + m σ n(nλ + M i ) n ≤ exp - nλ + m σ nλ + M i ,
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 σ .

B.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/2nL

Figure 2 :

 2 Figure 2: Comparison of optimization strategies that choose the best step-size in hindsight. In the top row are the quantum (left), protein (center), and sido (right) data sets. In the bottom row are the rcv1 and covertype data sets. This figure is best viewed in colour.

  λ + m σ /n p(λ + M j /n) p = 1 -nλ + m σ p(nλ + M j ) p ≤ exp -nλ + m σ nλ + M j .

While it may appear suboptimal to not use the gradients computed during the n iterations of stochastic gradient descent, using them only improves the bound by a constant.

Note that L in the SAG rates is based on the f i functions, while in the FG methods it is based on g which can be much smaller.

http://www.di.ens.fr/ ~mschmidt/Software/minFunc.html

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

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

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

We have also compared to a variety of other SG methods, such as SG with momentum, SG with gradient averaging, accelerated SG, and using SG but delaying averaging until after the first effective pass. However, none of these SG methods performed better than the ASG method above so we omit them to keep the plots simple.

Acknowledgements

Nicolas Le Roux, Mark Schmidt, and Francis Bach are supported by the European Research Council (SIERRA-ERC-239993). Mark Schmidt is also supported by a postdoctoral fellowship from the Natural Sciences and Engineering Research Council of Canada (NSERC).

The second term (within the expectation) on the right-hand side of Eq. ( 8) is equal to

The only random term (given F k-1 ) is the last one whose expectation is equal to

) .

The last term on the right-hand side of Eq. ( 8) is equal to

The only random term (given F k-1 ) is the last one whose expectation is equal to

Summing all these terms together, we get the following result:

Hence, the sum may be rewritten as

This concludes the proof.

This leads to (using the lemma of the previous section):

The third line is obtained using the Lipschitz property of the gradient, that is

where the inequality in the second line stems from [Nesterov, 2004, Theorem 2.1.5].

We have

We now use the strong convexity of g to get the inequality

This yields the final bound

Since we know that (x k-1 -x * ) g (x k-1 ) is positive, due to the convexity of g, we need to prove

We can then take a full expectation on both sides to obtain:

Since Q is a non-negative function (we show below that it dominates a non-negative function), this results proves the linear convergence of the sequence EQ(θ k ) with rate 1 -δ. We have

Step 2 -Domination of x k -x * 2 by Q(θ k )

We now need to prove that

We now use the strong convexity of the function to get the following inequalities:

.

Finally, we have

τ y,f τ y,g n τ y,I + τ y,e g (x k-1 ) 2 .

If we choose δ = δ n with δ 1 2 , ν = 1 2n , η = 2 and α = 1 2nµ , we get

This quantity is negative for δ

To finish the proof, we need to prove the positivity of the factor of g (x k-1 ) 2 . (τ y,f + τ y,g ) 2 -n 4

Then, following the same argument as in the previous section, we have

) 2 the variance of the gradients at the optimum.

Step 2 -Domination of g(x k ) -g(x * ) by Q(θ k )

We now need to prove that Q(θ k ) dominates g(x k ) -g(x * ).

using the convexity of g and the fact that

by dropping some terms.

The quantity on the right-hand side is minimized for e y = n 3 µ n+1

) . Hence, we have

This is most similar to the rate obtained with the dual coordinate descent method, but is likely to be slower because of the n term scaling M i . However, the difference will be decreased for overdetermined problems when m σ >> m σ .

Under the condition n 8L i g /µ g = 8(λ + M i )/(λ + m σ /n), with a step size of 1/2nµ g performing n iterations of the SAG algorithm has a rate of 1 -

Note that depending on the constants this may or may not not be faster than coordinate descent methods. However, if we consider the typical case where m σ = m σ = 0 with M i = O(p) and M j = O(n), then if we have n = 8(λ + M i )/λ we obtain 1 -

Despite the constant of 64 (which is likely to be highly sub-optimal), from these rates we see that SAG outperforms coordinate descent methods when n is sufficiently large.