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

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 Robbins and Monro (1951); Bottou and LeCun (2003). 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.,

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.

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 Teo et al. (2007), 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 Cauchy (1847), uses iterations of the form

Using $x^{\ast}$ to denote the unique minimizer of $g$, the FG method with a constant step size achieves a linear convergence rate:

for some $\rho<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

where $\alpha_{k}$ is a step-size and a training example $i_{k}$ is selected uniformly among the set $\{1,\dots,n\}$. The randomly chosen gradient $f_{i_{k}}^{\prime}(x^{k})$ yields an unbiased estimate of the true gradient $g^{\prime}(x^{k})$, and one can show under standard assumptions that, for a suitably chosen decreasing step-size sequence $\{\alpha_{k}\}$, the SG iterations achieve the sublinear convergence rate

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 Nemirovski and Yudin (1983); Nemirovski et al. (2009); Agarwal et al. (2012)). 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 Blatt et al. (2007), 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

where at each iteration a random training example $i_{k}$ is selected and we set

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

see Tseng (1998). It is common to set all $\beta_{k}=\beta$ for some constant $\beta$, and in this case we can rewrite the SG with momentum method as

We can re-write the SAG updates (5) in a similar form as

where the selection function $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 otherwise. Thus, momentum uses a geometric weighting of previous gradients while the SAG iterations select and average the most recent evaluation of each previous gradient. While momentum can lead to improved practical performance, it still requires the use of a decreasing sequence of step sizes and is not known to lead to a faster convergence rate.

Gradient Averaging: Closely related to momentum is using the sample average of all previous gradients,

which is similar to the SAG iteration in the form (5) but where all previous gradients are used. This approach is used in the dual averaging method Nesterov (2009), and while this averaging procedure leads to convergence for a constant step size and can improve the constants in the convergence rate Xiao (2010), 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 Polyak and Juditsky (1992). 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 Hazan and Kale (2011). 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 Nesterov (1983), 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 Schraudolph (1999); Sunehag et al. (2009); Ghadimi and Lan (2010); Martens (2010); Xiao (2010). Under certain conditions these variants are convergent with an $O(1/k)$ rate Sunehag et al. (2009). Alternately, if we split the convergence rate into a deterministic and stochastic part, these methods can improve the dependency on the deterministic part Ghadimi and Lan (2010); Xiao (2010). 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 and the second part is independent of $k$ but does not converge to . Thus, with a constant step size the SG iterations have a linear convergence rate up to some tolerance, and in general after this point the iterations do not make further progress. Indeed, 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}$ Solodov (1998). This contrasts with the method we present in this work which converges to the optimal solution using a constant step size and does so with a linear rate (without additional assumptions).

Accelerated methods: Accelerated SG methods, which despite their name are not related to the aforementioned AFG method, take advantage of the fast convergence rate of SG methods with a constant step size. In particular, accelerated SG methods use a constant step size by default, and only decrease the step size on iterations where the inner-product between successive gradient estimates is negative Kesten (1958); Delyon and Juditsky (1993). This leads to convergence of the method and allows it to potentially achieve periods of 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 Bertsekas (1997). 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 Friedlander and Schmidt (2012). 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 Blatt et al. (2007). 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

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

With a constant step size of $\alpha_{k}=\frac{1}{2nL}$, the SAG iterations satisfy for $k\geq 1$:

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

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 than $\frac{1}{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/\mu$ (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 $\mu$ and $L$ (see proof in the Appendix).

If $n\geqslant\frac{8L}{\mu}$, with a step size of $\alpha_{k}=\frac{1}{2n\mu}$ the SAG iterations satisfy for $k\geqslant n$:

We state this result for $k\geqslant 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.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. 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.

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}\leq\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\geqslant 8L/\mu$. Further, while the step size in Proposition 2 depends on $\mu$ and $n$, we can obtain the same convergence rate by using a step size as large as $\alpha_{k}=\frac{1}{16L}$. This is because the proposition is true for all values of $\mu$ satisfying $\frac{\mu}{L}\geqslant\frac{8}{n}$, so we can choose the smallest possible value of $\mu=\frac{8L}{n}$. We have observed in practice that the IAG method with a step size of $\alpha_{k}=\frac{1}{2n\mu}$ 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 $\frac{\mu}{2}\|x^{k}-x^{\ast}\|^{2}\leqslant g(x^{k})-g(x^{\ast})\leqslant\frac{L}{2}\|x^{k}-x^{\ast}\|^{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_{i}^{k}$ vectors is prohibitive, but we can often use structure in the $f_{i}^{\prime}$ to reduce this cost. For example, many loss functions $f_{i}$ take the form $f_{i}(a_{i}^{T}x)$ for a vector $a_{i}$. Since $a_{i}$ is constant, for these problems we only need to store the scalar $f_{i_{k}}^{\prime}(u_{i}^{k})$ for $u_{i}^{k}=a_{i_{k}}^{T}x^{k}$ rather than the full gradient $a_{i}^{T}f_{i}^{\prime}(u_{i}^{k})$, 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_{i}^{k}$ 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_{i}^{k}$ are set to the uninformative zero vector, rather than dividing $\alpha_{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.

This can be implemented efficiently for sparse data sets by using the representation $x=\kappa z$, where $\kappa$ is a scalar and $z$ is a vector, since the update based on the regularizer simply updates $\kappa$.

Large step sizes: Proposition 1 requires $\alpha_{k}\leqslant 1/2Ln$ while under an additional assumption Proposition 2 allows $\alpha_{k}\leqslant 1/16L$. In practice we observed better performance using step sizes of $\alpha_{k}=1/L$ and $\alpha_{k}=2/(L+n\mu)$. 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

To avoid instability caused by comparing very small numbers, we only do this test when $\|f_{i_{k}}^{\prime}(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 Nesterov (1983), where iterations of (3) with a fixed step size are interleaved with an extrapolation step, and we use an adaptive line-search based on Liu et al. (2009).

L-BFGS: A publicly-available limited-memory quasi-Newton method that has been tuned for log-linear models.http://www.di.ens.fr/~mschmidt/Software/minFunc.html 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 $\alpha_{k}=1/\mu 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 Xiao (2010), another recent state-of-the-art SG method.

ESG: The epoch SG method Hazan and Kale (2011), 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 $\alpha_{k}=2/(L_{k}+n\lambda)$ 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:

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.

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.

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

AFG: The accelerated full gradient method Nesterov (1983), 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 step-size.

ASG: The stochastic gradient method described by iteration (4), where we use a constant step size and average the iterates.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.

IAG: The incremental aggregated gradient method of Blatt et al. (2007) described by iteration (5) but with a cyclic choice of $i_{k}$.

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 $\xi$ (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 $\frac{\mu}{L}\geqslant\frac{8}{n}$ is satisfied when $\lambda\geqslant\frac{8\xi R^{2}}{n}$. In low-dimensional settings, the optimal regularization parameter is of the form $C/n$ Liang et al. (2009) where $C$ is a scalar constant, and may thus violate the constraint. However, the improvement with respect to regularization parameters of the form $\lambda=C/\sqrt{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 $\lambda$ which satisfy this constraint. For example, Sridharan et al. (2008) considers parameters of the form $\lambda=\frac{C}{\sqrt{n}}$ in the parametric setting, while Eberts and Steinwart (2011) considers $\lambda=\frac{C}{n^{\beta}}$ with $\beta<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 Bottou and Bousquet (2007).

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_{i}^{k}$ variables converges to $g^{\prime}(x^{k})$ as $\|x^{k}-x^{k-1}\|$ converges to zero, we can use $(1/n)\|\sum_{i}y_{i}^{k}\|$ 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.

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

Appendix

Appendix A Proofs of the propositions

We present here the proofs of Propositions 1 and 2.

For $k\geqslant 1$, the stochastic average gradient algorithm performs the recursion

where an $i_{k}$ is selected in $\{1,\dots,n\}$ uniformly at random and we set

Denoting $z^{k}_{i}$ a random variable which takes the value $1-\frac{1}{n}$ with probability $\frac{1}{n}$ and $-\frac{1}{n}$ otherwise (thus with zero expectation), this is equivalent to

Finally, if $M$ is a $tp\times tp$ matrix and $m$ is a $tp\times p$ matrix, then:

$\operatorname{diag}(M)$ is the $tp\times p$ matrix being the concatenation of the $t$ ($p\times p$)-blocks on the diagonal of $M$;

$\mathop{\rm Diag}(m)$ is the $tp\times tp$ block-diagonal matrix whose ($p\times p$)-blocks on the diagonal are equal to the ($p\times p$)-blocks of $m$.

A.2 Outline of the proofs

Each Proposition will be proved in multiple steps.

We shall prove that $Q(\theta^{k})$ dominates $\|x^{k}-x^{\ast}\|^{2}$ (in the case of Proposition 2) or $g(x^{k})-g(x^{\ast})$ (in the case of Proposition 2) by a constant for all $k$.

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, $\mathcal{F}_{k}$ will denote the $\sigma$-field of information up to (and including time $k$), i.e., $\mathcal{F}_{k}$ is the $\sigma$-field generated by $z^{1},\dots,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 $\sigma^{2}=\frac{1}{n}\sum_{i=1}^{n}\|f_{i}^{\prime}(x^{\ast})\|^{2}$ the variance of the gradients at the optimum. We will use the following recursion:

we have $\gamma_{k}\leqslant 2\gamma_{k}(1-\gamma_{k}L)$ and

Averaging from $i=0$ to $k-1$ and using the convexity of $g$, we have

A.4 Important lemma

In both proofs, our Lyapunov function contains a quadratic term R(\theta^{k})=(\theta^{k}-\theta^{*})^{\top}\left(\begin{array}[]{cc}A&b\\ b^{\top}&c\end{array}\right)(\theta^{k}-\theta^{*}) for some values of $A$, $b$ and $c$. The lemma below computes the value of $R(\theta^{k})$ in terms of elements of $\theta^{k-1}$.

Note that for square $n\times n$ matrix, $\operatorname{diag}(M)$ denotes a vector of size $n$ composed of the diagonal of $M$, while for a vector $m$ of dimension $n$, $\mathop{\rm Diag}(m)$ is the $n\times n$ diagonal matrix with $m$ on its diagonal. Thus $\mathop{\rm Diag}(\operatorname{diag}(M))$ is a diagonal matrix with the diagonal elements of $M$ on its diagonal, and $\operatorname{diag}(\mathop{\rm Diag}(m))=m$.

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

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

The only random term (given $\mathcal{F}_{k-1}$) is the third one whose expectation is equal to

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

The only random term (given $\mathcal{F}_{k-1}$) is the last one whose expectation is equal to

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

The only random term (given $\mathcal{F}_{k-1}$) is the last one whose expectation is equal to

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

with $S=A-\frac{\alpha}{n}be^{\top}-\frac{\alpha}{n}eb^{\top}+\frac{\alpha^{2}}{n^{2}}ece^{\top}=A-bc^{-1}b^{\top}+(b-\frac{\alpha}{n}ec)c^{-1}(b-\frac{\alpha}{n}ec)^{\top}$.

Rewriting $f^{\prime}(x^{k-1})-y^{k-1}=(f^{\prime}(x^{k-1})-f^{\prime}(x^{\ast}))-(y^{k-1}-f^{\prime}(x^{\ast}))$, we have

A.5 Analysis for α=12​n​L𝛼12𝑛𝐿\alpha=\frac{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, $\alpha=\frac{1}{2nL}$.

In this case, our Lyapunov function is quadratic, i.e.,

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

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

Note that for any symmetric negative definite matrix $M$ and for any vectors $s$ and $t$ we have

A sufficient condition for $M$ to be negative definite is to have $\delta\leqslant\frac{1}{3n}$.

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

Since we know that $(x^{k-1}-x^{\ast})^{\top}g^{\prime}(x^{k-1})$ is positive, due to the convexity of $g$, we need to prove that $\displaystyle\left(2\alpha-3\alpha^{2}nL+\frac{\delta^{2}\left(1-\frac{1}{n}\right)^{2}}{\left[3n\delta-1-2\delta+\frac{\delta-1}{n}\right]}\frac{n}{\mu}-\frac{\delta}{\mu}\right)$ is positive.

Using $\delta=\frac{\mu}{8nL}$ and $\alpha=\frac{1}{2nL}$ gives

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

We now need to prove that $Q(\theta^{k})$ dominates $\|x^{k}-x^{\ast}\|^{2}$. If \displaystyle P-\left(\begin{array}[]{cc}0&0\\ 0&\frac{1}{3}I\end{array}\right) is positive definite, then $Q(\theta^{k})\geqslant\frac{1}{3}\|x^{k}-x^{\ast}\|^{2}$.

We shall use the Schur complement condition for positive definiteness. Since $A$ is positive definite, the other condition to verify is $\frac{2}{3}I-b^{\top}A^{-1}b\succ 0$.

and so $P$ dominates \left(\begin{array}[]{cc}0&0\\ 0&\frac{1}{3}I\end{array}\right).

Initializing all the $y_{i}^{0}$ to 0, we get

A.6 Analysis for α=12​n​μ𝛼12𝑛𝜇\alpha=\frac{1}{2n\mu}

We now prove Proposition 2, providing a bound for the convergence rate of the SAG algorithm in the case of a small step size, $\alpha=\frac{1}{2n\mu}$.

We shall use the following Lyapunov function:

Our goal will now be to express all the quantities in terms of $(x^{k-1}-x^{\ast})^{\top}g^{\prime}(x^{k-1})$ whose positivity is guaranteed by the convexity of $g$.

Using the Lipschitz property of the gradients of $f_{i}$, we have

Using $e^{\top}[f^{\prime}(x^{k-1})-f^{\prime}(x^{\ast})]=ng^{\prime}(x^{k-1})$, we have

Reassembling all the terms together, we get

If we regroup all the terms in $[(x^{k-1})-(x^{\ast})]^{\top}g^{\prime}(x^{k-1})$ together, and all the terms in $(y^{k-1}-f^{\prime}(x^{\ast}))^{\top}$ together, we get

Assuming that $\tau_{y,I}$ and $\tau_{y,e}$ are negative, we have by completing the square that

where we used the fact that $(f^{\prime}(x^{k-1})-f^{\prime}(x^{\ast}))^{\top}e=g^{\prime}(x^{k-1})$. After reorganization of the terms, we obtain

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

If we choose $\delta=\frac{\widetilde{\delta}}{n}$ with $\widetilde{\delta}\leqslant\frac{1}{2}$, $\nu=\frac{1}{2n}$, $\eta=2$ and $\alpha=\frac{1}{2n\mu}$, we get

This quantity is negative for $\widetilde{\delta}\leqslant\frac{1}{3}$ and $\frac{\mu}{L}\geqslant\frac{4-6\widetilde{\delta}}{n(1-2\widetilde{\delta})(1-3\widetilde{\delta})}$. If we choose $\widetilde{\delta}=\frac{1}{8}$, then it is sufficient to have $\frac{n\mu}{L}\geqslant 8$.

To finish the proof, we need to prove the positivity of the factor of $\|g^{\prime}(x^{k-1})\|^{2}$.

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

with $\sigma^{2}=\frac{1}{n}\sum_{i}\|f_{i}^{\prime}(x^{\ast})\|^{2}$ the variance of the gradients at the optimum.

We now need to prove that $Q(\theta^{k})$ dominates $g(x^{k})-g(x^{\ast})$.

The quantity on the right-hand side is minimized for $e^{\top}y=\frac{n^{3}\mu}{n+1}\left(\frac{1}{n}(x^{k}-x^{\ast})-\frac{2\alpha}{n}g^{\prime}(x^{k})\right)$. Hence, we have

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

Appendix B Comparison of convergence rates

where to apply SG methods and SAG we can use

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

We can obtain the primal variables from the dual variables by the formula $x=(-1/\lambda)A^{\top}y$. Convergence rates of different primal and dual algorithms are often expressed in terms of the following Lipschitz constants:

Here, we use $M_{\sigma}$ to denote the maximum eigenvalue of $A^{\top}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}\{\sum_{i=1}^{n}(a_{i})^{2}_{j}\}$. We use $g_{j}^{\prime}$ to refer to element of $j$ of $g^{\prime}$, and similarly for $d_{i}^{\prime}$. The convergence rates will also depend on the primal and dual strong-convexity constants:

Here, $m_{\sigma}$ is the minimum eigenvalue of $A^{\top}A$, and $m_{\sigma}^{\prime}$ is the minimum eigenvalue of $AA^{\top}$.

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^{\ast}))$ converges to zero with rate

while a larger step-size of $2/(L_{g}+\mu_{g})$ gives a faster rate of

where the speed improvement is determined by the size of $m_{\sigma}$.

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^{\ast}))$ converges to zero with rate

and with a step-size of $2/(L_{d}+\mu_{d})$ the rate is

Thus, whether we can solve the primal or dual method faster depends on $m_{\sigma}$ and $m_{\sigma}^{\prime}$. 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_{\sigma}=m_{\sigma}^{\prime}$) or it has neither independent columns nor independent rows (so that $m_{\sigma}=m_{\sigma}^{\prime}=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)

and applied to the dual with a step-size of $1/L_{d}$ it has a rate of

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_{g}^{j}$ is (Nesterov, 2010, Theorem 2)

Here, we see that applying a coordinate-descent method can be much more efficient than an FG method if $M_{j}<<M_{\sigma}$. 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_{\sigma}$ 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_{d}^{i}$ this would obtain a rate on the dual objective of

which will be faster than the dual FG method if $M_{i}<<M_{\sigma}$. 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_{\sigma}$ compared to $m_{\sigma}^{\prime}$.

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_{g}$, performing $n$ iterations of the SAG algorithm has a rate of

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 over-determined problems when $m_{\sigma}>>m_{\sigma}^{\prime}$.

Under the condition $n\geqslant 8L_{g}^{i}/\mu_{g}=8(\lambda+M_{i})/(\lambda+m_{\sigma}/n)$, with a step size of $1/2n\mu_{g}$ performing $n$ iterations of the SAG algorithm has a rate of

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_{\sigma}=m_{\sigma}^{\prime}=0$ with $M_{i}=O(p)$ and $M_{j}=O(n)$, then if we have $n=8(\lambda+M_{i})/\lambda$ we obtain

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.

References