 There is widespread sentiment that it is not possible to effectively utilize fast gradient methods (e.g. Nesterov's acceleration, conjugate gradient, heavy ball) for the purposes of stochastic optimization due to their instability and error accumulation, a notion made precise in d'Aspremont 2008 and Devolder, Glineur, and Nesterov 2014. This work considers these issues for the special case of stochastic approximation for the least squares regression problem, and our main result refutes the conventional wisdom by showing that acceleration can be made robust to statistical errors. In particular, this work introduces an accelerated stochastic gradient method that provably achieves the minimax optimal statistical risk faster than stochastic gradient descent. Critical to the analysis is a sharp characterization of accelerated stochastic gradient descent as a stochastic process. We hope this characterization gives insights towards the broader question of designing simple and effective accelerated stochastic methods for more general convex and non-convex optimization problems.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Stochastic gradient descent (SGD) is the workhorse algorithm for optimization in machine learning and stochastic approximation problems; improving its runtime dependencies is a central issue in large scale stochastic optimization that often arise in machine learning problems at scale

(Bottou and Bousquet, 2007), where one can only resort to streaming algorithms.

This work examines these broader runtime issues for the special case of stochastic approximation in the following least squares regression problem:

 minx∈RdP(x),where% , P(x)def=12⋅E(a,b)∼D[(b−⟨x,a⟩)2], (1)

where we have access to a stochastic first order oracle, which, when provided with as an input, returns a noisy unbiased stochastic gradient using a tuple sampled from , with being the dimension of the problem. A query to the stochastic first-order oracle at produces:

 ˆ∇P(x)= −(b−⟨a,x⟩)⋅a. (2)

Note (i.e. eq(2

) is an unbiased estimate). Note that nearly all practical stochastic algorithms use sampled gradients of the specific form as in equation

2. We discuss differences to the more general stochastic first order oracle (Nemirovsky and Yudin, 1983) in section 1.4.

Let be a population risk minimizer. Given any estimation procedure which returns using samples, define the excess risk (which we also refer to as the generalization error or the error) of as . Now, equipped a stochastic first-order oracle (equation (2)), our goal is to provide a computationally efficient (and streaming) estimation method whose excess risk is comparable to the optimal statistical minimax rate.

In the limit of large , this minimax rate is achieved by the empirical risk minimizer (ERM), which is defined as follows. Given i.i.d. samples drawn from , define

 ˆxERMndef=argminxPn(x), where Pn(x)def=1nn∑i=112(bi−a⊤ix)2,

where denotes the ERM over the samples . For the case of additive noise models (i.e. where , with being independent of ), the minimax estimation rate is  (Kushner and Clark, 1978; Polyak and Juditsky, 1992; Lehmann and Casella, 1998; van der Vaart, 2000), i.e.:

 limn→∞ESn[P(ˆxERMn)]−P(x∗)dσ2/n =1, (3)

where

is the variance of the additive noise and the expectation is over the samples

drawn from . The seminal works of Ruppert (1988); Polyak and Juditsky (1992) proved that a certain averaged stochastic gradient method enjoys this minimax rate, in the limit. The question we seek to address is: how fast (in a non-asymptotic sense) can we achieve the minimax rate of ?

### 1.1 Review: Acceleration with Exact Gradients

Let us review results in convex optimization in the exact first-order oracle model. Running steps of gradient descent (Cauchy, 1847) with an exact first-order oracle yields the following guarantee:

 P(xt)−P(x∗)≤exp(−t/κo)⋅(P(x0)−P(x∗)),

where is the starting iterate, is the condition number of , where,

are the largest and smallest eigenvalue of the hessian

. Thus gradient descent requires oracle calls to solve the problem to a given target accuracy, which is sub-optimal amongst the class of methods with access to an exact first-order oracle (Nesterov, 2004). This sub-optimality can be addressed through Nesterov’s Accelerated Gradient Descent (Nesterov, 1983), which when run for t-steps, yields the following guarantee:

 P(xt)−P(x∗)≤exp(−t/√κo)⋅(P(x0)−P(x∗)),

which implies that oracle calls are sufficient to achieve a given target accuracy. This matches the oracle lower bounds (Nesterov, 2004) that state that calls to the exact first order oracle are necessary to achieve a given target accuracy. The conjugate gradient method (Hestenes and Stiefel, 1952) and heavy ball method (Polyak, 1964) are also known to obtain this convergence rate for solving a system of linear equations and for quadratic functions. These methods are termed fast gradient methods owing to the improvements offered by these methods over Gradient Descent.

This paper seeks to address the question: “Can we accelerate stochastic approximation in a manner similar to what has been achieved with the exact first order oracle model?”

### 1.2 A thought experiment: Is Accelerating Stochastic Approximation possible?

Let us recollect known results in stochastic approximation for the least squares regression problem (in equation 1). Running -steps of tail-averaged SGD (Jain et al., 2016) (or, streaming SVRG (Frostig et al., 2015b)222Streaming SVRG does not function in the stochastic first order oracle model (Frostig et al., 2015b)) provides an output that satisfies the following excess risk bound:

 E[P(ˆxn)]−P(x∗)≤exp(−n/κ)⋅(P(x0)−P(x∗))+2σ2d/n, (4)

where is the condition number of the distribution, which can be upper bounded as , assuming that

with probability one (refer to section

2 for a precise definition of ). Under appropriate assumptions, these are the best known rates under the stochastic first order oracle model (see section 1.4 for further discussion). A natural implication of the bound implied by averaged SGD is that with oracle calls (Jain et al., 2016) (where, hides factors in ), the excess risk attains (up to constants) the (asymptotic) minimax statistical rate. Note that the excess risk bounds in stochastic approximation consist of two terms: (a) bias: which represents the dependence of the generalization error on the initial excess risk , and (b) the variance: which represents the dependence of the generalization error on the noise level in the problem.

A precise question regarding accelerating stochastic approximation is: “is it possible to improve the rate of decay of the bias term, while retaining (up to constants) the statistical minimax rate?” The key technical challenge in answering this question is in sharply characterizing the error accumulation of fast gradient methods in the stochastic approximation setting. Common folklore and prior work suggest otherwise: several efforts have attempted to quantify instabilities in the face of statistical or non-statistical errors (Paige, 1971; Proakis, 1974; Polyak, 1987; Greenbaum, 1989; Roy and Shynk, 1990; Sharma et al., 1998; d’Aspremont, 2008; Devolder et al., 2013, 2014; Yuan et al., 2016). Refer to section 1.4 for a discussion on robustness of acceleration to error accumulation.

Optimistically, as suggested by the gains enjoyed by accelerated methods in the exact first order oracle model, we may hope to replace the oracle calls achieved by averaged SGD to . We now provide a counter example, showing that such an improvement is not possible. Consider a (discrete) distribution where the input is the

standard basis vector with probability

, . The covariance of in this case is a diagonal matrix with diagonal entries . The condition number of this distribution is . In this case, it is impossible to make non-trivial reduction in error by observing fewer than samples, since with constant probability, we would not have seen the vector corresponding to the smallest probability.

On the other hand, consider a case where the distribution is a Gaussian with a large condition number . Matrix concentration informs us that (with high probability and irrespective of how large is) after observing samples, the empirical covariance matrix will be a spectral approximation to the true covariance matrix, i.e. for some constant , . Here, we may hope to achieve a faster convergence rate, as information theoretically samples suffice to obtain a non-trivial statistical estimate (see Hsu et al. (2014) for further discussion).

Figure 1 shows the behavior of SGD in these cases; both are synthetic examples in dimensions, with a condition number and noise level . See the figure caption for more details.

These examples suggest that if acceleration is indeed possible, then the degree of improvement (say, over averaged SGD) must depend on distributional quantities that go beyond the condition number . A natural conjecture is that this improvement must depend on the number of samples required to spectrally approximate the covariance matrix of the distribution; below this sample size it is not possible to obtain any non-trivial statistical estimate due to information theoretic reasons. This sample size is quantified by a notion which we refer to as the statistical condition number (see section 2 for a precise definition and for further discussion about ). As we will see in section 2, we have , is affine invariant, unlike (i.e.

is invariant to linear transformations over

).

### 1.3 Contributions

This paper introduces an accelerated stochastic gradient descent scheme, which can be viewed as a stochastic variant of Nesterov’s accelerated gradient method (Nesterov, 2012). As pointed out in Section 1.2, the excess risk of this algorithm can be decomposed into two parts namely, bias and variance. For the stochastic approximation problem of least squares regression, this paper establishes bias contraction at a geometric rate of , improving over prior results (Frostig et al., 2015b; Jain et al., 2016),which prove a geometric rate of , while retaining statistical minimax rates (up to constants) for the variance. Here is the condition number and is the statistical condition number of the distribution, and a rate of is an improvement over since (see Section 2 for definitions and a short proof of ).

See Table 1 for a theoretical comparison. Figure 2 provides an empirical comparison of the proposed (tail-averaged) accelerated algorithm to (tail-averaged) SGD (Jain et al., 2016) on our two running examples. Our result gives improvement over SGD even in the noiseless (i.e. realizable) case where ; this case is equivalent to the setting where we have a distribution over a (possibly infinite) set of consistent linear equations. See Figure 3 for a comparison on the case where .

On a more technical note, this paper introduces two new techniques in order to analyze the proposed accelerated stochastic gradient method: (a) the paper introduces a new potential function in order to show faster rates of decaying the bias, and (b) the paper provides a sharp understanding of the behavior of the proposed accelerated stochastic gradient descent updates as a stochastic process and utilizes this in providing a near-exact estimate of the covariance of its iterates. This viewpoint is critical in order to prove that the algorithm achieves the statistical minimax rate.

We use the operator viewpoint for analyzing stochastic gradient methods, introduced in Défossez and Bach (2015). This viewpoint was also used in Dieuleveut and Bach (2015); Jain et al. (2016).

### 1.4 Related Work

##### Non-asymptotic Stochastic Approximation:

Stochastic gradient descent (SGD) and its variants are by far the most widely studied algorithms for the stochastic approximation problem. While initial works (Robbins and Monro, 1951) considered the final iterate of SGD, later works (Ruppert, 1988; Polyak and Juditsky, 1992) demonstrated that averaged SGD obtains statistically optimal estimation rates. Several works provide non-asymptotic analyses for averaged SGD and variants (Bach and Moulines, 2011; Bach, 2014; Frostig et al., 2015b) for various stochastic approximation problems. For stochastic approximation with least squares regression Bach and Moulines (2013); Défossez and Bach (2015); Needell et al. (2016); Frostig et al. (2015b); Jain et al. (2016)

provide non-asymptotic analysis of the behavior of SGD and its variants.

Défossez and Bach (2015); Dieuleveut and Bach (2015) provide non-asymptotic results which achieve the minimax rate on the variance (where the bias is lower order, not geometric). Needell et al. (2016) achieves a geometric rate on the bias (and where the variance is not minimax). Frostig et al. (2015b); Jain et al. (2016) obtain both the minimax rate on the variance and a geometric rate on the bias, as seen in equation 4.

##### Acceleration and Noise Stability:

While there have been several attempts at understanding if it is possible to accelerate SGD , the results have been largely negative. With regards to acceleration with adversarial (non-statistical) errors in the exact first order oracle model, d’Aspremont (2008) provide negative results and Devolder et al. (2013, 2014) provide lower bounds showing that fast gradient methods do not improve upon standard gradient methods. There is also a series of works considering statistical errors. Polyak (1987) suggests that the relative merits of heavy ball (HB) method (Polyak, 1964) in the noiseless case vanish with noise unless strong assumptions on the noise model are considered; an instance of this is when the noise variance decays as the iterates approach the minimizer. The Conjugate Gradient (CG) method (Hestenes and Stiefel, 1952) is suggested to face similar robustness issues in the face of statistical errors (Polyak, 1987); this is in addition to the issues that CG is known to suffer from owing to roundoff errors (due to finite precision arithmetic) (Paige, 1971; Greenbaum, 1989). In the signal processing literature, where SGD goes by Least Mean Squares (LMS) (Widrow and Stearns, 1985), there have been efforts that date to several decades (Proakis, 1974; Roy and Shynk, 1990; Sharma et al., 1998) which study accelerated LMS methods (stochastic variants of CG/HB) in the same oracle model as the one considered by this paper (equation 2). These efforts consider the final iterate (i.e. no iterate averaging) of accelerated LMS methods with a fixed step-size and conclude that while it allows for a faster decay of the initial error (bias) (which is unquantified), their steady state behavior (i.e. variance) is worse compared to that of LMS. Yuan et al. (2016) considered a constant step size accelerated scheme with no iterate averaging in the same oracle model as this paper, and conclude that these do not offer any improvement over standard SGD. More concretely, Yuan et al. (2016) show that the variance of their accelerated SGD method with a sufficiently small constant step size is the same as that of SGD with a significantly larger step size. Note that none of the these efforts (Proakis, 1974; Roy and Shynk, 1990; Sharma et al., 1998; Yuan et al., 2016) achieve minimax error rates or quantify (any improvement whatsoever on the) rate of bias decay.

##### Oracle models and optimality:

With regards to notions of optimality, there are (at least) two lines of thought: one is a statistical objective where the goal is (on every problem instance) to match the rate of the statistically optimal estimator  (Anbar, 1971; Fabian, 1973; Kushner and Clark, 1978; Polyak and Juditsky, 1992); another is on obtaining algorithms whose worst case upper bounds (under various assumptions such as bounded noise) match the lower bounds provided in  Nemirovsky and Yudin (1983). The work of Polyak and Juditsky (1992) are in the former model, where they show that the distribution of the averaged SGD estimator matches, on every problem, that of the statistically optimal estimator, in the limit (under appropriate regularization conditions standard in the statistics literature, where the optimal estimator is often referred to as the maximum likelihood estimator/the empirical risk minimizer/an -estimator (Lehmann and Casella, 1998; van der Vaart, 2000)). Along these lines, non-asymptotic rates towards statistically optimal estimators are given by Bach and Moulines (2013); Bach (2014); Défossez and Bach (2015); Dieuleveut and Bach (2015); Needell et al. (2016); Frostig et al. (2015b); Jain et al. (2016). This work can be seen as improving this non-asymptotic rate (to the statistically optimal estimation rate) using an accelerated method. As to the latter (i.e. matching the worst-case lower bounds in  Nemirovsky and Yudin (1983)), there are a number of positive results on using accelerated stochastic optimization procedures; the works of Lan (2008); Hu et al. (2009); Ghadimi and Lan (2012, 2013); Dieuleveut et al. (2016) match the lower bounds provided in Nemirovsky and Yudin (1983). We compare these assumptions and works in more detail.

In stochastic first order oracle models (see  Kushner and Clark (1978); Kushner and Yin (2003)), one typically has access to sampled gradients of the form:

 ˆ∇P(x)=∇P(x)+η, (5)

where varying assumptions are made on the noise . The worst-case lower bounds in Nemirovsky and Yudin (1983) are based on that is bounded; the accelerated methods in Lan (2008); Hu et al. (2009); Ghadimi and Lan (2012, 2013); Dieuleveut et al. (2016) which match these lower bounds in various cases, all assume either bounded noise or, at least is finite. In the least squares setting (such as the one often considered in practice and also considered in Polyak and Juditsky (1992); Bach and Moulines (2013); Défossez and Bach (2015); Dieuleveut and Bach (2015); Frostig et al. (2015b); Jain et al. (2016)), this assumption does not hold, since is not bounded. To see this, in our oracle model (equation 2) is:

 η=ˆ∇P(x)−∇P(x)=(aa⊤−H)(x−x∗)−ϵ⋅a (6)

which implies that is not uniformly bounded (unless additional assumptions are enforced to ensure that the algorithm’s iterates lie within a compact set). Hence, the assumptions made in Hu et al. (2009); Ghadimi and Lan (2012, 2013); Dieuleveut et al. (2016) do not permit one to obtain finite -sample bounds on the excess risk. Suppose we consider the case of , i.e. where the additive noise is zero and . For this case, this paper provides a geometric rate of convergence to the minimizer , while the results of Ghadimi and Lan (2012, 2013); Dieuleveut et al. (2016) at best indicate a rate. Finally, in contrast to all other existing work, our result is the first to provide finer distribution dependent characteristics of the improvements offered by accelerating SGD (e.g. refer to the Gaussian and discrete examples in section 1.2).

##### Acceleration and Finite Sums:

As a final remark, there have been results (Shalev-Shwartz and Zhang, 2014; Frostig et al., 2015a; Lin et al., 2015; Lan and Zhou, 2015; Allen-Zhu, 2016) that provide accelerated rates for offline stochastic optimization which deal with minimizing sums of convex functions; these results are almost tight due to matching lower bounds (Lan and Zhou, 2015; Woodworth and Srebro, 2016). These results do not immediately translate into rates on the generalization error. Furthermore, these algorithms are not streaming, as they require making multiple passes over a dataset stored in memory. Refer to Frostig et al. (2015b) for more details.

## 2 Main Results

We now provide our assumptions and main result, before which, we have some notation. For a vector and a positive semi-definite matrix (i.e. ), denote .

### 2.1 Assumptions and Definitions

Let

denote the second moment matrix of the input, which is also the hessian

of (1):

 Hdef=E(a,b)∼D[a⊗a]=∇2P(x).

Furthermore, let the fourth moment tensor

of the inputs is defined as:

 M=E(a,b)∼D[a⊗a⊗a⊗a].
1. [leftmargin=*,label=]

2. Finite second and fourth moment: The second moment matrix and the fourth moment tensor exist and are finite.

3. Positive Definiteness: The second moment matrix is strictly positive definite, i.e. .

We assume 1 and 22 implies that is strongly convex and admits a unique minimizer . Denote the noise in a sample as: . First order optimality conditions of imply

 ∇P(x∗)=E[ϵ⋅a]=0.

Let denote the covariance of gradient at optimum (or noise covariance matrix),

 Σdef=E(a,b)∼D[ˆ∇P(x∗)⊗ˆ∇P(x∗)]=E(a,b)∼D[ϵ2⋅a⊗a].

We define the noise level , condition number , statistical condition number below.
Noise level: The noise level is defined to be the smallest positive number such that

 Σ⪯σ2H.

The noise level quantifies the amount of noise in the stochastic gradient oracle and has been utilized in previous work (e.g., see Bach and Moulines (2011, 2013)) for providing non-asymptotic bounds for the stochastic approximation problem. In the homoscedastic (additive noise/well specified) case, where is independent of the input , this condition is satisfied with equality, i.e. with .
Condition number: Let

 μdef=λmin(H).

by 2. Now, let be the smallest positive number such that

 E[∥a∥2 aa⊤]⪯R2 H.

. The condition number of the distribution  (Défossez and Bach, 2015; Jain et al., 2016) is

 κdef=R2/μ.

Statistical condition number: The statistical condition number is defined as the smallest positive number such that

 E[∥a∥2H−1aa⊤]⪯˜κ H.

Remarks on and : Unlike , it is straightforward to see that is affine invariant (i.e. unchanged with linear transformations over ). Since , we note . For the discrete case (from Section 1.2), it is straightforward to see that both and are equal to . In contrast, for the Gaussian case (from Section 1.2), is , while is which may be arbitrarily large (based on choice of the coordinate system).

governs how many samples require to be drawn from so that the empirical covariance is spectrally close to , i.e. for some constant , . In comparison to the matrix Bernstein inequality where stronger (yet related) moment conditions are assumed in order to obtain high probability results, our results hold only in expectation (refer to Hsu et al. (2014) for this definition, wherein is referred to as bounded statistical leverage in theorem and remark ).

### 2.2 Algorithm and Main Theorem

Algorithm 1 presents the pseudo code of the proposed algorithm. ASGD can be viewed as a variant of Nesterov’s accelerated gradient method (Nesterov, 2012), working with a stochastic gradient oracle (equation 2) and with tail-averaging the final iterates. The main result now follows:

###### Theorem 1.

Suppose  1 and  2 hold. Set . After calls to the stochastic first order oracle (equation 2), ASGD outputs satisfying:

 E[P(¯xt,n)]−P(x∗)≤C⋅(κ˜κ)9/4dκ(n−t)2⋅exp(−t9√κ˜κ)⋅(P(x0)−P(x∗))Leading order % bias error+5σ2dn−tLeading order % variance error+ C⋅(κ˜κ)5/4dκ⋅exp(−n9√κ˜κ)(P(x0)−P(x∗))Exponentially vanishing lower% order bias term+C⋅σ2d(n−t)2√κ˜κLower order variance error term+ C⋅exp(−n9√κ˜κ)⋅(σ2d⋅(κ˜κ)7/4+σ2d(n−t)2⋅(κ˜κ)7/2˜κ)+C⋅σ2dn−t(κ˜κ)11/4exp(−(n−t−1)30√κ˜κ)Exponentially vanishing lower order % variance error terms,

where is a universal constant, , and are the noise level, condition number and statistical condition number respectively.

The following corollary holds if the iterates are tail-averaged over the last samples and . The second condition lets us absorb lower order terms into leading order terms.

###### Corollary 2.

Assume the parameter settings of theorem 1 and let and (for an appropriate universal constants ). We have that with calls to the stochastic first order oracle, ASGD outputs a vector satisfying:

 E[P(¯xt,n)]−P(x∗)≤C⋅exp(−n20√κ˜κ)⋅(P(x0)−P(x∗))+11σ2dn.

A few remarks about the result of theorem 1 are due: (i) ASGD decays the initial error at a geometric rate of during the unaveraged phase of iterations, which presents the first improvement over the rate offered by SGD (Robbins and Monro, 1951)/averaged SGD (Polyak and Juditsky, 1992; Jain et al., 2016) for the least squares stochastic approximation problem, (ii) the second term in the error bound indicates that ASGD obtains (up to constants) the minimax rate once . Note that this implies that Theorem 1 provides a sharp non-asymptotic analysis (up to factors) of the behavior of Algorithm 1.

### 2.3 Discussion and Open Problems

A challenging problem in this context is in formalizing a finite sample size lower bound in the oracle model considered in this work. Lower bounds in stochastic oracle models have been considered in the literature (see Nemirovsky and Yudin (1983); Raginsky and Rakhlin (2011); Agarwal et al. (2012)), though it is not evident these oracle models and lower bounds are sharp enough to imply statements in our setting (see section 1.4 for a discussion of these oracles).

Let us now understand theorem 1 in the broader context of stochastic approximation. Under certain regularity conditions, it is known that (Lehmann and Casella, 1998; van der Vaart, 2000) that the rate described in equation 3

for the homoscedastic case holds for a broader set of misspecified models (i.e., heteroscedastic noise case), with an appropriate definition of the noise variance. By defining

, the rate of the ERM is guaranteed to approach  (Lehmann and Casella, 1998; van der Vaart, 2000) in the limit of large , i.e.:

 limn→∞ESn[Pn(ˆxERMn)]−P(x∗)σ2ERM/n =1, (7)

where is the ERM over samples . Averaged SGD (Jain et al., 2016) and streaming SVRG (Frostig et al., 2015b) are known to achieve these rates for the heteroscedastic case. Refer to Frostig et al. (2015b) for more details.Neglecting constants, Theorem 1 is guaranteed to achieve the rate of the ERM for the homoscedastic case (where ) and is tight when the bound is nearly tight (upto constants). We conjecture ASGD achieves the rate of the ERM in the heteroscedastic case by appealing to a more refined analysis as is the case for averaged SGD (see Jain et al. (2016)). It is also an open question to understand acceleration for smooth stochastic approximation (beyond least squares), in situations where the rate represented by equation 7 holds (Polyak and Juditsky, 1992).

## 3 Proof Outline

We now present a brief outline of the proof of Theorem 1. Recall the variables in Algorithm 1. Before presenting the proof outline we require some definitions. We begin by defining the centered estimate as:

 θjdef=[xj−x∗yj−x∗]∈R2d.

Recall that the stepsizes in Algorithm 1 are . The accelerated SGD updates of Algorithm 1 can be written in terms of as:

 θj=ˆAjθj−1+ζj, where, ˆAjdef=[0(I−δaja⊤j)−α(1−β) I(1+α(1−β))I−(αδ+(1−α)γ)aja⊤j] ,ζjdef=[δ⋅ϵjaj(αδ+(1−α)γ)⋅ϵjaj],

where . The tail-averaged iterate is associated with its own centered estimate . Let , where is a filtration generated by . Let be linear operators acting on a matrix so that , , . Denote and matrices as:

Bias-variance decomposition: The proof of theorem 1 employs the bias-variance decomposition, which is well known in the context of stochastic approximation (see Bach and Moulines (2011); Frostig et al. (2015b); Jain et al. (2016)) and is re-derived in the appendix. The bias-variance decomposition allows for the generalization error to be upper-bounded by analyzing two sub-problems: (a) bias, analyzing the algorithm’s behavior on the noiseless problem (i.e. a.s.) while starting at and (b) variance, analyzing the algorithm’s behavior by starting at the solution (i.e. ) and allowing the noise to drive the process. In a similar manner as , the bias and variance sub-problems are associated with and , and these are related as:

 E[¯θt,n⊗¯θt,n]⪯2⋅(E[¯θbiast,n⊗¯θbiast,n]+E[¯θvariancet,n⊗¯θvariancet,n]). (8)

Since we deal with the square loss, the generalization error of the output of algorithm 1 is:

 E[P(¯xt,n)]−P(x∗)=12⋅⟨[H000],E[¯θt,n⊗¯θt,n]⟩, (9)

indicating that the generalization error can be bounded by analyzing the bias and variance sub-problem. We now present the lemmas that bound the bias error.

###### Lemma 3.

The covariance of the bias part of averaged iterate satisfies:

 E[¯θbiast,n⊗¯θbiast,n] =1(n−t)2(I+(I−AL)−1AL+(I−A⊤R)−1A⊤R)(I−B)−1(Bt+1−Bn+1)(θ0⊗θ0) −1(n−t)2n∑j=t+1((I−AL)−1An+1−jL+(I−A⊤R)−1(A⊤R)n+1−j)Bj(θ0⊗θ0).

The quantity that needs to be bounded in the term above is . Lemma 4 presents a result that can be applied recursively to bound ( since ).

###### Lemma 4 (Bias contraction).

For any two vectors , let . We have:

 ⟨G,B(θθ⊤)⟩≤(1−19√κ˜κ)⟨G,θθ⊤⟩
##### Remarks:

(i) the matrices and appearing in are due to the fact that we prove contraction using the variables and instead of and , as used in defining . (ii) The key novelty in lemma 4 is that while standard analyses of accelerated gradient descent (in the exact first order oracle) use the potential function (e.g. Wilson et al. (2016)), we consider it crucial for employing the potential function (this potential function is captured using the matrix ) to prove accelerated rates (of ) for bias decay.

We now present the lemmas associated with bounding the variance error:

###### Lemma 5.

The covariance of the variance error satisfies:

 E[¯θvariancet,n⊗¯θvariancet,n]=1n−t(I+(I−AL)−1AL+(I−A⊤R)−1A⊤R)(I−B)−1ˆΣ −1(n−t)2((I−AL)−2(AL−An+1−tL)+(I−A⊤R)−2(A⊤R−(A⊤R)n+1−t))(I−B)−1ˆΣ +1(n−t)2n∑j=t+1((I−AL)−1An+1−jL+(I−A⊤R)−1(A⊤R)n+1−j)(I−B)−1BjˆΣ.

The covariance of the stationary distribution of the iterates i.e., requires a precise bound to obtain statistically optimal error rates. Lemma 6 presents a bound on this quantity.

###### Lemma 6 (Stationary covariance).

The covariance of limiting distribution of satisfies:

 E[θvariance∞⊗θvariance∞]=(I−B)−1ˆΣ ⪯5σ2((2/3)⋅(1˜κH−1)+(5/6)⋅(δI))⊗.

A crucial implication of this lemma is that the limiting final iterate has an excess risk . This result naturally lends itself to the (tail-)averaged iterate achieving the minimax optimal rate of . Refer to the appendix E and lemma 17 for more details in this regard.

## 4 Conclusion

This paper introduces an accelerated stochastic gradient method, which presents the first improvement in achieving minimax rates faster than averaged SGD (Robbins and Monro, 1951; Ruppert, 1988; Polyak and Juditsky, 1992; Jain et al., 2016)/Streaming SVRG (Frostig et al., 2015b) for the stochastic approximation problem of least squares regression. To obtain this result, the paper presented the need to rethink what acceleration has to offer when working with a stochastic gradient oracle: these thought experiments indicated a need to consider a quantity that captured more fine grained problem characteristics. The statistical condition number (an affine invariant distributional quantity) is shown to characterize the improvements that acceleration offers in the stochastic first order oracle model.

In essence, this paper presents the first known provable analysis of the claim that fast gradient methods are stable when dealing with statistical errors, in contrast to negative results in statistical and non-statistical settings (Paige, 1971; Proakis, 1974; Polyak, 1987; Greenbaum, 1989; Roy and Shynk, 1990; Sharma et al., 1998; d’Aspremont, 2008; Devolder et al., 2013, 2014; Yuan et al., 2016). We hope that this paper provides insights towards developing simple and effective accelerated stochastic gradient schemes for general convex and non-convex optimization problems.

##### Acknowledgments:

Sham Kakade acknowledges funding from Washington Research Foundation Fund for Innovation in Data-Intensive Discovery and the NSF through awards CCF-, CCF- and CCF-.

## References

• Agarwal et al. (2012) A. Agarwal, P. L. Bartlett, P. Ravikumar, and M. J. Wainwright. Information-theoretic lower bounds on the oracle complexity of stochastic convex optimization. IEEE Transactions on Information Theory, 2012.
• Allen-Zhu (2016) Z. Allen-Zhu. Katyusha: The first direct acceleration of stochastic gradient methods. CoRR, abs/1603.05953, 2016.
• Anbar (1971) D. Anbar. On Optimal Estimation Methods Using Stochastic Approximation Procedures. University of California, 1971.
• Bach (2014) F. R. Bach.

Adaptivity of averaged stochastic gradient descent to local strong convexity for logistic regression.

Journal of Machine Learning Research (JMLR), volume 15, 2014.
• Bach and Moulines (2011) F. R. Bach and E. Moulines. Non-asymptotic analysis of stochastic approximation algorithms for machine learning. In NIPS 24, 2011.
• Bach and Moulines (2013) F. R. Bach and E. Moulines. Non-strongly-convex smooth stochastic approximation with convergence rate O(1/n). In NIPS 26, 2013.
• Bottou and Bousquet (2007) L. Bottou and O. Bousquet. The tradeoffs of large scale learning. In NIPS 20, 2007.
• Cauchy (1847) L. A. Cauchy. Méthode générale pour la résolution des systémes d’équations simultanees. C. R. Acad. Sci. Paris, 1847.
• d’Aspremont (2008) A. d’Aspremont. Smooth optimization with approximate gradient. SIAM Journal on Optimization, 19(3):1171–1183, 2008.
• Défossez and Bach (2015) A. Défossez and F. R. Bach. Averaged least-mean-squares: Bias-variance trade-offs and optimal sampling distributions. In AISTATS, volume 38, 2015.
• Devolder et al. (2013) O. Devolder, F. Glineur, and Y. E. Nesterov. First-order methods with inexact oracle: the strongly convex case. CORE Discussion Papers 2013016, 2013.
• Devolder et al. (2014) O. Devolder, F. Glineur, and Y. E. Nesterov. First-order methods of smooth convex optimization with inexact oracle. Mathematical Programming, 146:37–75, 2014.
• Dieuleveut and Bach (2015) A. Dieuleveut and F. R. Bach. Non-parametric stochastic approximation with large step sizes. The Annals of Statistics, 2015.
• Dieuleveut et al. (2016) A. Dieuleveut, N. Flammarion, and F. R. Bach. Harder, better, faster, stronger convergence rates for least-squares regression. CoRR, abs/1602.05419, 2016.
• Fabian (1973) V. Fabian. Asymptotically efficient stochastic approximation; the RM case. Annals of Statistics, 1(3), 1973.
• Frostig et al. (2015a) R. Frostig, R. Ge, S. Kakade, and A. Sidford. Un-regularizing: approximate proximal point and faster stochastic algorithms for empirical risk minimization. In ICML, 2015a.
• Frostig et al. (2015b) R. Frostig, R. Ge, S. M. Kakade, and A. Sidford. Competing with the empirical risk minimizer in a single pass. In COLT, 2015b.
• Ghadimi and Lan (2012) S. Ghadimi and G. Lan. Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization i: A generic algorithmic framework. SIAM Journal on Optimization, 2012.
• Ghadimi and Lan (2013) S. Ghadimi and G. Lan. Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization, ii: shrinking procedures and optimal algorithms. SIAM Journal on Optimization, 2013.
• Greenbaum (1989) A. Greenbaum. Behavior of slightly perturbed lanczos and conjugate-gradient recurrences. Linear Algebra and its Applications, 1989.
• Hestenes and Stiefel (1952) M. R. Hestenes and E. Stiefel. Methods of conjuate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 1952.
• Hsu et al. (2014) D. J. Hsu, S. M. Kakade, and T. Zhang.

Random design analysis of ridge regression.

Foundations of Computational Mathematics, 14(3):569–600, 2014.
• Hu et al. (2009) C. Hu, J. T. Kwok, and W. Pan. Accelerated gradient methods for stochastic optimization and online learning. In NIPS 22, 2009.
• Jain et al. (2016) P. Jain, S. M. Kakade, R. Kidambi, P. Netrapalli, and A. Sidford. Parallelizing stochastic approximation through mini-batching and tail-averaging. CoRR, abs/1610.03774, 2016.
• Kushner and Clark (1978) H. J. Kushner and D. S. Clark. Stochastic Approximation Methods for Constrained and Unconstrained Systems. Springer-Verlag, 1978.
• Kushner and Yin (2003) H. J. Kushner and G. Yin. Stochastic approximation and recursive algorithms and applications. Springer-Verlag, 2003.
• Lan (2008) G. Lan. An optimal method for stochastic composite optimization. Tech. Report, IE, Georgia Tech., 2008.
• Lan and Zhou (2015) G. Lan and Y. Zhou. An optimal randomized incremental gradient method. CoRR, abs/1507.02000, 2015.
• Lehmann and Casella (1998) E. L. Lehmann and G. Casella. Theory of Point Estimation. Springer Texts in Statistics. Springer, 1998. ISBN 9780387985022.
• Lin et al. (2015) H. Lin, J. Mairal, and Z. Harchaoui. A universal catalyst for first-order optimization. In NIPS, 2015.
• Needell et al. (2016) D. Needell, N. Srebro, and R. Ward. Stochastic gradient descent, weighted sampling, and the randomized kaczmarz algorithm. Mathematical Programming, 2016.
• Nemirovsky and Yudin (1983) A. S. Nemirovsky and D. B. Yudin. Problem Complexity and Method Efficiency in Optimization. John Wiley, 1983.
• Nesterov (1983) Y. E. Nesterov. A method for unconstrained convex minimization problem with the rate of convergence . Doklady AN SSSR, 269, 1983.
• Nesterov (2004) Y. E. Nesterov. Introductory lectures on convex optimization: A basic course, volume 87 of Applied Optimization. Kluwer Academic Publishers, 2004.
• Nesterov (2012) Y. E. Nesterov. Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization, 22(2):341–362, 2012.
• Paige (1971) C. C. Paige.

The computation of eigenvalues and eigenvectors of very large sparse matrices.

PhD Thesis, University of London, 1971.
• Polyak (1964) B. T. Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4, 1964.
• Polyak (1987) B. T. Polyak. Introduction to Optimization. Optimization Software, 1987.
• Polyak and Juditsky (1992) B. T. Polyak and A. B. Juditsky. Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, volume 30, 1992.
• Proakis (1974) J. G. Proakis. Channel identification for high speed digital communications. IEEE Transactions on Automatic Control, 1974.
• Raginsky and Rakhlin (2011) M. Raginsky and A. Rakhlin. Information-based complexity, feedback and dynamics in convex programming. IEEE Transactions on Information Theory, 2011.
• Robbins and Monro (1951) H. Robbins and S. Monro. A stochastic approximation method. The Annals of Mathematical Statistics, vol. 22, 1951.
• Roy and Shynk (1990) S. Roy and J. J. Shynk. Analysis of the momentum lms algorithm. IEEE Transactions on Acoustics, Speech and Signal Processing, 1990.
• Ruppert (1988) D. Ruppert. Efficient estimations from a slowly convergent robbins-monro process. Tech. Report, ORIE, Cornell University, 1988.
• Shalev-Shwartz and Zhang (2014) S. Shalev-Shwartz and T. Zhang. Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization. In ICML, 2014.
• Sharma et al. (1998) R. Sharma, W. A. Sethares, and J. A. Bucklew. Analysis of momentum adaptive filtering algorithms. IEEE Transactions on Signal Processing, 1998.
• van der Vaart (2000) A. W. van der Vaart. Asymptotic Statistics. Cambridge University Publishers, 2000.
• Widrow and Stearns (1985) B. Widrow and S. D. Stearns. Adaptive Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1985.
• Wilson et al. (2016) A. C. Wilson, B. Recht, and M. I. Jordan. A lyapunov analysis of momentum methods in optimization. CoRR, abs/1611.02635, 2016.
• Woodworth and Srebro (2016) B. Woodworth and N. Srebro. Tight complexity bounds for optimizing composite objectives. CoRR, abs/1605.08003, 2016.
• Yuan et al. (2016) K. Yuan, B. Ying, and A. H. Sayed. On the influence of momentum acceleration on online learning. Journal of Machine Learning Research (JMLR), volume 17, 2016.

## Appendix A Appendix setup

We will first provide a note on the organization of the appendix and follow that up with introducing the notations.

### a.1 Organization

• In subsection A.2, we will recall notation from the main paper and introduce some new notation that will be used across the appendix.

• In section B, we will write out expressions that characterize the generalization error of the proposed accelerated SGD method. In order to bound the generalization error, we require developing an understanding of two terms namely the bias error and the variance error.

• In section C, we prove lemmas that will be used in subsequent sections to prove bounds on the bias and variance error.

• In section D, we will bound the bias error of the proposed accelerated stochastic gradient method. In particular, lemma 4 is the key lemma that provides a new potential function with which this paper achieves acceleration. Further, lemma 16 is the lemma that bounds all the terms of the bias error.

• In section E, we will bound the variance error of the proposed accelerated stochastic gradient method. In particular, lemma 6 is the key lemma that considers a stochastic process view of the proposed accelerated stochastic gradient method and provides a sharp bound on the covariance of the stationary distribution of the iterates. Furthermore, lemma 20 bounds all terms of the variance error.

• Section F presents the proof of Theorem 1. In particular, this section aggregates the result of lemma 16 (which bounds all terms of the bias error) and lemma 20 (which bounds all terms of the variance error) to present the guarantees of Algorithm 1.

### a.2 Notations

We begin by introducing , which is the fourth moment tensor of the input , i.e.:

 Mdef=E(a,b)∼D[a⊗a⊗a⊗a]

Applying the fourth moment tensor to any matrix produces another matrix in that is expressed as:

 MSdef=E[(a⊤Sa)aa⊤].

With this definition in place, we recall as the smallest number, such that

applied to the identity matrix

satisfies:

 MI=E[∥a∥22aa⊤]⪯R2 H

Moreover, we recall that the condition number of the distribution , where is the smallest eigenvalue of . Furthermore, the definition of the statistical condition number of the distribution follows by applying the fourth moment tensor to , i.e.:

 MH−1 =E[(a⊤H−1a)⋅aa⊤]⪯˜κ H

We denote by and the left and right multiplication operator of any matrix , i.e. for any matrix , and .

Parameter choices: In all of appendix we choose the parameters in Algorithm 1 as

 α=√κ˜κc2√2c1−c21+√κ˜κ,  β=c3c2√2c1−c21√κ˜κ,  γ=c2√2c1−c21μ√κ˜κ,  δ=c1R2

where is an arbitrary constant satisfying . Furthermore, we note that , and . Note that we recover Theorem 1 by choosing . We denote

 cdef=α(1−β)% and, q