ACCELERATING SGD WITH MOMENTUM FOR OVER- PARAMETERIZED LEARNING
Year: 2020

Abstract:
Nesterov SGD is widely used for training modern neural networks and other machine learning models. Yet, its advantages over SGD have not been theoretically clarified. Indeed, as we show in this paper, both theoretically and empirically, Nesterov SGD with any parameter selection does not in general provide acceleration over ordinary SGD. Furthermore, Nesterov SGD may diverge for step sizes that ensure convergence of ordinary SGD. This is in contrast to the classical results in the deterministic setting, where the same step size ensures accelerated convergence of the Nesterov’s method over optimal gradient descent. To address the non-acceleration issue, we introduce a compensation term to Nesterov SGD. The resulting algorithm, which we call MaSS, converges for same step sizes as SGD. We prove that MaSS obtains an accelerated convergence rates over SGD for any mini-batch size in the linear setting. For full batch, the convergence rate of MaSS matches the well-known accelerated rate of the Nesterov’s method. We also analyze the practically important question of the dependence of the convergence rate and optimal hyper-parameters on the mini-batch size, demonstrating three distinct regimes: linear scaling, diminishing returns and saturation. Experimental evaluation of MaSS for several standard architectures of deep networks, including ResNet and convolutional networks, shows improved performance over SGD, SGD+Nesterov and Adam.

heading: 1 INTRODUCTION
text: Many modern neural networks and other machine learning models are over-parametrized (5). These models are typically trained to have near zero training loss, known as interpolation and often have strong generalization performance, as indicated by a range of empirical evidence including (23; 3). Due to a key property of interpolation – automatic variance reduction (discussed in Section 2.1), stochastic gradient descent (SGD) with constant step size is shown to converge to the optimum of a convex loss function for a wide range of step sizes (12). Moreover, the optimal choice of step size η∗ for SGD in that setting can be derived analytically.
The goal of this paper is to take a step toward understanding momentum-based SGD in the interpolating setting. Among them, stochastic version of Nesterov’s acceleration method (SGD+Nesterov) is arguably the most widely used to train modern machine learning models in practice. The popularity of SGD+Nesterov is tied to the well-known acceleration of the deterministic Nesterov’s method over gradient descent (15). Yet, has not not theoretically clear whether Nesterov SGD accelerates over SGD.
As we show in this work, both theoretically and empirically, Nesterov SGD with any parameter selection does not in general provide acceleration over ordinary SGD. Furthermore, Nesterov SGD may diverge, even in the linear setting, for step sizes that guarantee convergence of ordinary SGD. Intuitively, the lack of acceleration stems from the fact that, to ensure convergence, the step size of SGD+Nesterov has to be much smaller than the optimal step size for SGD. This is in contrast to the deterministic Nesterov method, which accelerates using the same step size as optimal gradient descent. As we prove rigorously in this paper, the slow-down of convergence caused by the small step
size negates the benefit brought by the momentum term. We note that a similar lack of acceleration for the stochastic Heavy Ball method was analyzed in (9).
To address the non-acceleration of SGD+Nesterov, we introduce an additional compensation term to allow convergence for the same range of step sizes as SGD. The resulting algorithm, MaSS (Momentum-added Stochastic Solver)1 updates the weights w and u using the following rules (with the compensation term underlined):
wt+1 ← ut − η1∇̃f(ut), ut+1 ← (1 + γ)wt+1 − γwt + η2∇̃f(ut).
(1)
Here, ∇̃ represents the stochastic gradient. The step size η1, the momentum parameter γ ∈ (0, 1) and the compensation parameter η2 are independent of t. We proceed to analyze theoretical convergence properties of MaSS in the interpolated regime. Specifically, we show that in the linear setting MaSS converges exponen-
tially for the same range of step sizes as plain SGD, and the optimal choice of step size for MaSS is exactly η∗ which is optimal for SGD. Our key theoretical result shows that MaSS has accelerated convergence rate over SGD. Furthermore, in the full batch (deterministic) scenario, our analysis selects η2 = 0, thus reducing MaSS to the classical Nesterov’s method (15). In this case our convergence rate also matches the well-known convergence rate for the Nesterov’s method (15; 4). This acceleration is illustrated in Figure 1. Note that SGD+Nesterov (as well as Stochastic Heavy Ball) does not converge faster than SGD, in line with our theoretical analysis. We also prove exponential convergence of MaSS in more general convex setting under additional conditions.
We further analyze the dependence of the convergence rate e−s(m)t and optimal hyper-parameters on the mini-batch size m. We identify three distinct regimes of dependence defined by two critical values m∗1 and m ∗ 2: linear scaling, diminishing returns and saturation, as illustrated in Figure 2. The convergence speed per iteration s(m), as well as the optimal hyper-parameters, increase linearly as m in the linear scaling regime, sub-linearly in the diminishing returns regime, and can only increase by a small constant factor in the saturation regime. The critical values m∗1 and m ∗ 2 are derived analytically. We note that the intermediate “dimin-
ishing terurns” regime is new and is not found in SGD (12). To the best of our knowledge, this is the first analysis of mini-batch dependence for accelerated stochastic gradient methods.
We also experimentally evaluate MaSS on deep neural networks, which are non-convex. We show that MaSS outperforms SGD, SGD+Nesterov and Adam (10) both in optimization and generalization, on different architectures of deep neural networks including convolutional networks and ResNet (7).
The paper is organized as follows: In section 2, we introduce notations and preliminary results. In section 3, we discuss the non-acceleration of SGD+Nesterov. In section 4 we introduce MaSS and analyze its convergence and optimal hyper-parameter selection. In section 5, we analyze the mini-batch MaSS. In Section 6, we show experimental results.

heading: 1.1 RELATED WORK
text: Over-parameterized models have drawn increasing attention in the literature as many modern machine learning models, especially neural networks, are over-parameterized (5) and show strong generalization performance (16; 23; 2). Over-parameterized models usually result in nearly perfect fit (or interpolation) of the training data (23; 18; 3). Exponential convergence of SGD with constant step size under interpolation and its dependence on the batch size is analyzed in (12).
1Code url: https://github.com/ts66395/MaSS
There are a few works that show or indicate the non-acceleration of existing stochastic momentum methods. First of all, the work (9) theoretically proves non-acceleration of stochastic Heavy Ball method (SGD+HB) over SGD on certain synthetic data. Furthermore, these authors provide experimental evidence that SGD+Nesterov also converges at the same rate as SGD on the same data. The work (22) theoretically shows that, for sufficiently small step-sizes, SGD+Nesterov and SGD+HB is equivalent to SGD with a larger step size. However, the results in (22) do not exclude the possibility that acceleration is possible when the step size is larger. The work (11) concludes that “momentum hurts the convergence within the neighborhood of global optima”, based on a theoretical analysis of SGD+HB. These results are consistent with our analysis of the standard SGD+Nesterov. However, this conclusion does not apply to all momentum methods. Indeed, we will show that MaSS provably improves convergence over SGD.
There is a large body of work, both practical and theoretical, on SGD with momentum, including (10; 8; 1). Adam (10), and its variant AMSGrad (17), are among the most practically used SGD methods with momentum. Unlike our method, Adam adaptively adjusts the step size according to a weight-decayed accumulation of gradient history. In (8) the authors proposed an accelerated SGD algorithm, which can be written in the form shown on the right hand side in Eq.8, but with different hyper-parameter selection. Their ASGD algorithm also has a tail-averaging step at the final stage. In the interpolated setting (no additive noise) their analysis yields a convergence rate of O(Poly(κ, κ̃) exp(− t
9 √ κ1κ̃ )) compared to O(exp(− t√ κ1κ̃ )) for our algorithm with batch size 1. We provide some experimental comparisons between their ASGD algorithm and MaSS in Fig. 4.
The work (21) proposes and analyzes another first-order momentum algorithm and derives convergence rates under a different set of conditions – the strong growth condition for the loss function in addition to convexity. As shown in Appendix F.3, on the example of a Gaussian distributed data, the rates obtained in (21) can be slower than those for SGD. In contrast, our algorithm is guaranteed to never have a slower convergence rate than SGD. Furthermore, in the same Gaussian setting MaSS matches the optimal accelerated full-gradient Nesterov rate.
Additionally, in our work we consider the practically important dependence of the convergence rate and optimal parameter selection on the mini-batch size, which to the best of our knowledge, has not been analyzed for momentum methods.

heading: 2 NOTATIONS AND PRELIMINARIES
text: Given dataset D = {(xi, yi)}ni=1 ⊂ Rd × C, we consider an objective function of the form f(w) = 1 n ∑n i=1 fi(w), where fi only depends on a single data point (xi, yi). Let ∇f denote the exact gradient, and ∇̃mf denote the unbiased stochastic gradient evaluated based on a mini-batch of size m. For simplicity, we also denote ∇̃f(w) := ∇̃1f(w).We use the concepts of strong convexity and smoothness of functions, see definitions in Appendix B.1. For loss function with µ-strong convexity and L-smoothness, the condition number κ is defined as κ = L/µ.
In the case of the square loss, fi(w) = 12 (w Txi − yi)2, and the Hessian matrix is H := 1 n ∑n i=1 xix T i . Let L and µ be the largest and the smallest non-zero eigenvalues of the Hessian respectively. Then the condition number is then κ = L/µ (note that zero eigenvalues can be ignored in our setting, see Section 4).
Stochastic Condition Numbers. For a quadratic loss function, let H̃m := 1m ∑m i=1 x̃ix̃ T i denotes
a mini-batch estimate of H . Define L1 be the smallest positive number such that E [ ‖x̃‖2x̃x̃T ] L1H, and denote Lm := L1/m+ (m− 1)L/m. (2)
Given a mini-batch size m, we define the m-stochastic condition number as κm := Lm/µ.
Following (8), we introduce the quantity κ̃ (called statistical condition number in (8)), which is the smallest positive real number such that E [ ‖x̃‖2H−1 x̃x̃ T ] κ̃H.
Remark 1. We note that L1 ≥ L, since E [ ‖x̃‖2x̃x̃T ] = E[H̃21 ] = H2 + E[(H̃1 − H)2] H2. Consequently, according to the definition of Lm in Eq.2, we have
Lm ≥ L, ∀m ≥ 1; Lm → L, as m→∞. (3)
Hence, the quadratic loss function is also Lm-smooth, for all m ≥ 1. By the definition of κm, we also have
κm ≥ κ, ∀m ≥ 1; κm → κ, as m→∞. (4) Remark 2. It is important to note that κ̃ ≤ κ1, since E [ ‖x̃‖2H−1 x̃x̃ T ] 1µE [ ‖x̃‖2x̃x̃T ] κ1H .

heading: 2.1 CONVERGENCE OF SGD FOR OVER-PARAMETRIZED MODELS AND OPTIMAL STEP SIZE
text: We consider over-parametrized models that have zero training loss solutions on the training data (e.g., (23)). A solution fi(w) which fits the training data perfectly fi(w) = 0, ∀i = 1, 2, · · · , n, is known as interpolating. In the linear setting, interpolation implies that the linear system {xTi w = yi}ni=1 has at least one solution. A key property of interpolation is Automatic Variance Reduction (AVR), where the variance of the stochastic gradient decreases to zero as the weight w approaches the optimal w∗.
Var[∇̃mf(w)] ‖w −w∗‖2E[(H̃m −H)2]. (5) For a detailed discussion of AVR see Appendix B.2.
Thanks to AVR, plain SGD with constant step size can be shown to converge exponentially for strongly convex loss functions (13; 19; 14; 12). The set of acceptable step sizes is (0, 2/Lm), where Lm is defined in Eq.2 and m is the mini-batch size. Moreover, the optimal step size η∗(m) of SGD that induces fastest convergence guarantee is proven to be 1/Lm (12).

heading: 3 NON-ACCELERATION OF SGD+NESTEROV
text: In this section we prove that SGD+Nesterov, with any constant hyper-parameter setting, does not generally improve convergence over optimal SGD. Specifically, we demonstrate a setting where SGD+Nesterov can be proved to have convergence rate of (1 − O(1/κ))t, which is same (up to a constant factor) as SGD. In contrast, the classical accelerated rate for the deterministic Nesterov’s method is (1− 1/ √ κ)t.
We will consider the following two-dimensional data-generating component decoupled model. Fix an arbitrary w∗ ∈ R2 and randomly sample z from N (0, 2). The data points (x, y) are constructed as follow:
x = { σ1 · z · e1, w.p. 0.5, σ2 · z · e2, w.p. 0.5, and y = 〈w ∗,x〉, (6)
where e1, e2 ∈ R2 are canonical basis vectors, σ1 > σ2 > 0. It can be seen that L1 = σ21 , κ1 = 6σ 2 1/σ 2 2 and κ̃ = 6 (See Appendix F.1). This model is similar to that used to analyze the stochastic Heavy Ball method in (9).
The following theorem gives a lower bound for the convergence of SGD+Nesterov, regarding the linear regression problem on the component decoupled data model. See Appendix C for the proof. Theorem 1 (Non-acceleration of SGD+Nesterov). Let {(xi, yi)}ni=1 be a dataset generated according to the component decoupled model. Consider the optimization problem of minimizing quadratic function 12n ∑ i(x T i w − yi)2. For any step size η > 0 and momentum parameter γ ∈ (0, 1) of SGD+Nesterov with random initialization, with probability one, there exists a T ∈ N such that ∀t > T ,
E [f(wt)]− f(w∗) ≥ (
1− C κ
)t [ f(w0)− f(w∗) ] , (7)
where C > 0 is a constant.
Compared with the convergence rate (1−1/κ)t of SGD (12), this theorem shows that SGD+Nesterov does not accelerate over SGD. This result is very different from that in the deterministic gradient scenario, where the classical Nesterov’s method has a strictly faster convergence guarantee than gradient descent (15).
Intuitively, the key reason for the non-acceleration of SGD+Nesterov is a condition on the step size η required for non-divergence of the algorithm. Specifically, when momentum parameter γ is close to
1, η is required to be less than 2(1−γ)/3+O((1−γ)2) (precise formulation is given in Lemma 1 in Appendix C). The slow-down resulting from the small step size necessary to satisfy that condition cannot be compensated by the benefit of the momentum term. In particular, the condition on the step-size of SGD+Nesterov excludes η∗ that achieves fastest convergence for SGD. We show in the following corollary that, with the step size η∗, SGD+Nesterov diverges. This is different from the deterministic scenario, where the Nesterov method accelerates using the same step size as gradient descent. Corollary 1. Consider the same optimization problem as in Theorem 1. Let step-size η = 1L1 = 1 6σ21 and acceleration parameter γ ∈ [0.6, 1], then SGD+Nesterov, with random initialization, diverges with probability 1.
We empirically verify the non-acceleration of SGD+Nesterov as well as Corollary 1, in Section 6 and Appendix F.2.

heading: 4 MASS: ACCELERATING SGD
text: In this section, we propose MaSS, which introduces a compensation term (see Eq.1) onto SGD+Nesterov. We show that MaSS can converge exponentially for all the step sizes that result in convergence of SGD, i.e., η ∈ (0, 2/Lm). Importantly, we derive a convergence rate exp(−t/ √ κ1κ̃), where κ̃ ≤ κ1, for MaSS which is faster than the convergence rate for SGD exp(−t/κ1). Moreover, we give an analytical expression for the optimal hyper-parameter setting. For ease of analysis, we rewrite update rules of MaSS in Eq.1 in the following equivalent form (introducing an additional variable v):{
wt+1 ← ut − η1∇̃f(ut), ut+1 ← (1 + γ)wt+1 − γwt + η2∇̃f(ut) ⇐⇒  wt+1 ← ut − η∇̃f(ut), vt+1 ← (1− α)vt + αut − δ∇̃f(ut), ut+1 ← α1+αvt+1 + 1 1+αwt+1. (8)
There is a bijection between the hyper-parameters (η1, η2, γ) and (η, α, δ), which is given by: γ = (1− α)/(1 + α), η1 = η, η2 = (η − αδ)/(1 + α). (9)
Remark 3 (SGD+Nesterov). In the literature, the Nesterov’s method is sometimes written in a similar form as the R.H.S. of Eq.8. Since SGD+Nesterov has no compensation term, δ has to be fixed as η/α, which is consistent with the parameter setting in (15).
Assumptions. We first assume square loss function, and later extend the analysis to general convex loss functions under additional conditions. For square loss function, the solution setW∗ := {w ∈ Rd|f(w) = 0} is an affine subspace in the parameter space Rd. Given any w, we denote its closest solution as w∗ := arg minv∈W∗ ‖w−v‖, and define the error = w−w∗. Be aware that different w may correspond to different w∗, and that and (stochastic) gradients are always perpendicular to W∗ (see discussion in Appendix B.3). Hence, no actual update happens alongW∗. For this reason, we can ignore zero eigenvalues of H and restrict our analysis to the span of the eigenvectors of the Hessian with non-zero eigenvalues.
Based on the equivalent form of MaSS in Eq.8, the following theorem shows that, for square loss function in the interpolation setting, MaSS is guaranteed to have exponential convergence when hyper-parameters satisfy certain conditions. Theorem 2 (Convergence of MaSS). Consider minimizing a quadratic loss function in the interpolation setting. Let µ be the smallest non-zero eigenvalue of the Hessian matrix H . Let Lm be as defined in Eq.2. Denote κ̃m := κ̃/m + (m − 1)/m. In MaSS with mini batch of size m, if the positive hyper-parameters η, α, δ satisfy the following two conditions:
α/δ ≤ µ, αδκ̃m + η(ηLm − 2) ≤ 0, (10) then, after t iterations,
E [ ‖vt −w∗‖2H−1 + δ
α ‖wt −w∗‖2
] ≤ (1− α)t ( ‖v0 −w∗‖2H−1 + δ
α ‖w0 −w∗‖2
) .
Consequently, ‖wt −w∗‖2 ≤ C · (1− α)t, for some constant C > 0 which depends on the initialization.
Remark 4. By condition Eq.10, the admissible step size η is (0, 2/Lm), exactly the same as SGD for interpolated setting (12). Remark 5. One can easily check that the hyper-parameter setting of SGD+Nesterov does not satisfy the conditions in Eq.10. Proof sketch for Theorem 2. Denote Ft := E [ ‖vt+1 −w∗‖2H−1 + δ α‖wt+1 −w ∗‖2 ] , we show that, under the update rules of MaSS in Eq.8,
Ft+1 ≤ (1− α)Ft + + (α/µ− δ)︸ ︷︷ ︸ c1
‖ut −w∗‖2 + ( δ2κ̃m + δη 2Lm/α− 2ηδ/α )︸ ︷︷ ︸
c2
‖ut −w∗‖2H .
By the condition in Eq.10, c1 ≤ 0, c2 ≤ 0, then the last two terms are non-positive. Hence,
Ft+1 ≤ (1− α)Ft ≤ (1− α)t+1F0. (11)
Using that ‖wt−w∗‖2 ≤ αFt/δ, we get the final conclusion. See detailed proof in Appendix D.
Hyper-parameter Selection. From Theorem 2, we observe that the convergence rate is determined by (1−α)t. Therefore, larger α is preferred for faster convergence. Combining the conditions in Eq.10, we have α2 ≤ η(2− ηLm)µ/κ̃m. (12) By setting η∗ = 1/Lm, which maximizes the right hand side of the inequality, we obtain the optimal selection α∗ = 1/ √ κmκ̃m. Note that this setting of η∗ and α∗ determines a unique δ∗ = α∗/µ by the conditions in Eq.10. In summary,
η∗(m) = 1
Lm , α∗(m) = 1√ κmκ̃m
, δ∗(m) = η∗
α∗κ̃m . (13)
By Eq.9, the optimal selection of (η1, η2, γ) would be:
η∗1(m) = 1
Lm , η∗2(m) =
η∗1 √ κmκ̃m√
κmκ̃m + 1
( 1− 1
κ̃m
) , γ∗(m) =
√ κmκ̃m − 1√ κmκ̃m + 1 . (14)
κ̃m is usually larger than 1, which implies that the coefficient η∗2 of the compensation term is non-negative. The non-negative coefficient η2 indicates that the weight ut is “over-descended” in SGD+Nesterov and needs to be compensated along the gradient direction.
It is important to note that the optimal step size for MaSS as in Eq.13 is exactly the same as the optimal one for SGD (12). With such hyper-parameter selection given in Eq.14, we have the following theorem for optimal convergence: Theorem 3 (Acceleration of MaSS). Under the same assumptions as in Theorem 2, if we set hyperparameters in MaSS as in Eq.13, then after t iteration of MaSS with mini batch of size m,
‖wt −w∗‖2 ≤ C · (1− 1/ √ κmκ̃m) t, (15)
for some constant C > 0 which depends on the initialization.
Remark 6. With the optimal hyper-parameters in Eq.13, the asymptotic convergence rate of MaSS is
O(e−t/ √ κmκ̃m), (16)
which is faster than the rate O(e−t/κm) of SGD (see (12)), since κm ≥ κ̃m. Remark 7 (MaSS Reduces to the Nesterov’s method for full batch). In the limit of full batch m→ ∞, we have κm → κ, κ̃m → 1, the optimal parameter selection in Eq.14 reduces to
η∗1 = 1/L, γ ∗ = (
√ κ− 1)/( √ κ+ 1), η∗2 = 0. (17)
It is interesting to observe that, in the full batch (deterministic) scenario, the compensation term vanishes and η∗1 and γ
∗ are the same as those in Nesterov’s method. Hence MaSS with the optimal hyper-parameter selection reduces to Nesterov’s method in the limit of full batch. Moreover, the convergence rate in Theorem 3 reduces toO(e−t/ √ κ), which is exactly the well-known convergence rate of Nesterov’s method (15; 4).
Extension to Convex Case. First, we extend the definition of L1 to convex functions, L1 := inf{c ∈ R | E[‖∇̃f(w)‖2] ≤ 2c (f(w)− f∗)}, and keep the definition of Lm the same as Eq.2. It can be shown that these definitions of Lm are consistent with those in the quadratic setting. We assume that the smallest positive eigenvalue of Hessian H(x) is lower bounded by µ > 0, for all x.
Theorem 4. Suppose there exists a 1L -strongly convex and 1 µ -smooth non-negative function g : Rd → R such that g(w∗) = 0 and 〈∇g(x),∇f(z)〉 ≥ (1 − )〈x −w∗, z −w∗〉,∀x, z ∈ Rd, for some > 0. In MaSS, if the hyper-parameters are set to be:
η = 1/(2Lm), α = (1− )/(2κm), δ = 1/(2Lm), (18) then after t iterations, there exists a constant C such that, E[f(wt)] ≤ C · ( 1− 1− 2κm )t .

heading: 5 LINEAR SCALING RULE AND THE DIMINISHING RETURNS
text: Based on our analysis, we discuss the effect of selection of mini-batch size m. We show that the domain of mini-batch size m can be partitioned into three intervals by two critical points: m∗1 = min(L1/L, κ̃),m ∗ 2 = max(L1/L, κ̃). The three intervals/regimes are depicted in Figure 2, and the detailed analysis is in Appendix G.
Linear Scaling: m < m∗1. In this regime, we have Lm ≈ L1/m, κm ≈ κ1/m and κ̃m ≈ κ̃/m. The optimal selection of hyper-parameters is approximated by:
η∗(m) ≈ m · η∗(1), α∗(m) ≈ m · α∗(1), δ∗(m) ≈ m · δ∗(1),
and the convergence rate in Eq.16 is approximately O(e−m·t/ √ κ1κ̃). This indicates linear gains in convergence when m increases.
In the linear scaling regime, the hyper-parameter selections follow a Linear Scaling Rule (LSR): When the mini-batch size is multiplied by k, multiply all hyperparameters (η, α, δ) by k. This parallels the linear scaling rule for SGD which is an accepted practice for training neural networks (6).
Moreover, increasing m results in linear gains in the convergence speed, i.e., one MaSS iteration with mini-batch size m is almost as effective as m MaSS iterations with mini-batch size 1.
Diminishing Returns: m ∈ [m∗1,m∗2) In this regime, increasing m results in sublinear gains in convergence speed. One MaSS iteration with mini-batch size m is less effective than m MaSS iterations with mini-batch size 1.
Saturation: m ≥ m∗2. One MaSS iteration with mini-batch size m is nearly as effective (up to a multiplicative factor of 2) as one iteration with full gradient.
This three regimes partition is different from that for SGD (12), where only linear scaling and saturation regimes present. An empirical verification of the dependence of the convergence speed on m is shown in Figure 3. See the setup in Appendix G.

heading: 6 EMPIRICAL EVALUATION
text: Synthetic Data. We empirically verify the non-acceleration of SGD+Nesterov and the fast convergence of MaSS on synthetic data. Specifically, we optimize the quadratic function 12n ∑ i(x T i w − yi) 2, where the dataset {(xi, yi)}ni=1 is generated by the component decoupled model described in Section 3. We compare the convergence behavior of SGD+Nesterov with SGD, as well as our proposed method, MaSS, and several other methods: SGD+HB, ASGD (8). We select the best hyper-parameters from dense grid search for SGD+Nesterov (step-size and momentum parameter),
SGD+HB (step-size and momentum parameter) and SGD (step-size). For MaSS, we do not tune the hyper-parameters but use the hyper-parameter setting suggested by our theoretical analysis in Section 4; For ASGD, we use the setting provided by (8).
Fig. 4 shows the convergence behaviors of these algorithms on the setting of σ21 = 1, σ 2 2 = 1/2 9. We observe that the fastest convergence of SGD+Nesterov is almost identical to that of SGD, indicating the non-acceleration of SGD+Nesterov. We also observe that our proposed method, MaSS, clearly outperforms the others. In Appendix F.2, we provide additional experiments on more settings of the component decoupled data, and Gaussian distributed data. We also show the divergence of SGD+Nesterov with the same step size as SGD and MaSS in Appendix F.2.
Real data: MNIST and CIFAR-10. We compare the optimization performance of SGD, SGD+Nesterov and MaSS on the following tasks: classification of MNIST with a fully-
connected network (FCN), classification of CIFAR-10 with a convolutional neural network (CNN) and Gaussian kernel regression on MNIST. See detailed description of the architectures in Appendix H.1. In all the tasks and for all the algorithms, we select the best hyper-parameter setting over dense grid search, except that we fix the momentum parameter γ = 0.9 for both SGD+Nesterov and MaSS, which is typically used in practice. All algorithms are implemented with mini batches of size 64 for neural network training.
Fig. 5 shows the training curves of MaSS, SGD+Nesterov and SGD, which indicate the fast convergence of MaSS on real tasks, including the non-convex optimization problems on neural networks.
Figure 5: Comparison of SGD, SGD+Nesterov and MaSS on (left) fully-connected neural network, (middle) convolutional neural network, and (right) kernel regression.
Test Performance. We show that the solutions found by MaSS have good generalization performance. We evaluate the classification accuracy of MaSS, and compare with SGD, SGD+Nesterov and Adam, on different modern neural networks: CNN and ResNet (7). See description of the architectures in Appendix H.1. In the training processes, we follow the standard protocol of data augmentation and reduction of learning rate, which are typically used to achieve state-of-the-art results in neural networks. In each task, we use the same initial learning rate for MaSS, SGD and SGD+Nesterov, and run the same number of epochs (150 epochs for CNN and 300 epochs for ResNet-32). Detailed experimental settings are deferred to Appendix H.2.
Table 6 compares the classification accuracy of these algorithms on the test set of CIFAR-10 (average of 3 independent runs).
We observe that MaSS produces the best test performance. We also note that increasing initial learning rate may improves performance of MaSS and SGD, but degrades that of
SGD+Nesterov. Moreover, in our experiment, SGD+Nesterov with large step size η = 0.3 diverges in 5 out of 8 runs on CNN and 2 out of 5 runs on ResNet-32 (for random initialization), while MaSS and SGD converge on every run.

heading: ACKNOWLEDGEMENTS
text: This research was in part supported by NSF funding and a Google Faculty Research Award. GPUs donated by Nvidia were used for the experiments. We thank Ruoyu Sun for helpful comments concerning convergence rates. We thank Xiao Liu for helping with the empirical evaluation of our proposed method.

heading: A PSEUDOCODE FOR MASS
text: Algorithm 1 : MaSS–Momentum-added Stochastic Solver Require: Step-size η1, secondary step-size η2, acceleration parameter γ ∈ (0, 1). Initialize: u0 = w0. while not converged do wt+1 ← ut − η1∇̃f(ut), ut+1 ← (1 + γ)wt+1 − γwt − η2∇̃f(ut).
end while Output: weight wt.
Note that the proposed algorithm initializes the variables w0 and u0 with the same vector, which could be randomly generated.
As discussed in section 4, MaSS can be equivalently implemented using the following update rules:
wt+1 ← ut − η∇̃f(ut), (19a) vt+1 ← (1− α)vt + αut − δ∇̃f(ut), (19b)
ut+1 ← α
1 + α vt+1 +
1
1 + α wt+1. (19c)
In this case, variables u0, v0 and w0 should be initialized with the same vector.
There is a bijection between the hyper-parameters (η1, η2, γ) and (η, α, δ), which is given by:
γ = 1− α 1 + α , η1 = η, η2 = η − αδ 1 + α . (20)

heading: B ADDITIONAL PRELIMINARIES

heading: B.1 STRONG CONVEXITY AND SMOOTHNESS OF FUNCTIONS
text: Definition 1 (Strong Convexity). A differentiable function f : Rd → R is µ-strongly convex (µ > 0), if
f(x) ≥ f(z) + 〈∇f(z),x− z〉+ µ 2 ‖x− z‖2, ∀x, z ∈ Rd. (21)
Definition 2 (Smoothness). A differentiable function f : Rd → R is L-smooth (L > 0), if
f(x) ≤ f(z) + 〈∇f(z),x− z〉+ L 2 ‖x− z‖2, ∀x, z ∈ Rd. (22)

heading: B.2 AUTOMATIC VARIANCE REDUCTION
text: In the interpolation setting, one can write the square loss as
f(w) = 1 2 (w −w∗)TH(w −w∗) = 1 2 ‖w −w∗‖2H . (23)
A key property of interpolation is that the variance of the stochastic gradient of decreases to zero as the weight w approaches an optimal solution w∗. Proposition 1 (Automatic Variance Reduction). For the square loss function f in the interpolation setting, the stochastic gradient at an arbitrary point w can be written as
∇̃mf(w) = H̃m(w −w∗) = H̃m . (24) Moreover, the variance of the stochastic gradient
Var[∇̃mf(w)] ‖ ‖2E[(H̃m −H)2]. (25)
Since E[(H̃m − H)2] is independent of w, the above proposition unveils a linear dependence of variance of stochastic gradient on the norm square of error . This observation underlies exponential convergence of SGD in certain convex settings (20; 13; 19; 14; 12).

heading: B.3 ON IGNORING ZERO EIGENVALUES OF THE HESSIAN
text: Consider the square loss function, f(w) = 12n ∑ i(x T i w − yi)2. The (stochastic) gradient is
∇̃mf(w) = 1
m m∑ i=1 xi(x T i w − yi) = 1 m m∑ i=1 xix T i (w −w∗) = Hm(w −w∗),
where the stochastic gradient is computed based on a randomly sampled batch of size m.
Recall that the solution set W∗ := {w ∈ Rd|f(w) = 0} is an affine subspace in the parameter space, and that w∗ is the solution in W∗ that is closest to w. Hence, w − w∗ is perpendicular to W∗, i.e., w − w∗ ∈ (W∗)⊥. Note that (W∗)⊥ is an invariant space of Hm, hence ∇̃mf(w) = Hm(w −w∗) ∈ (W∗)⊥. Hence the (stochastic) gradient is perpendicular toW∗.

heading: C PROOF OF THEOREM 1
text: The key proof technique is to consider the asymptotic behavior of SGD+Nesterov in the decoupled model of data when the condition number becomes large.
Notations and proof setup. Recall that the square loss function based on the component decoupled dataD, define in Eq.6, is in the interpolation regime, then for SGD+Nesterov, we have the recurrence relation
wt+1 = (1 + γ)(1− ηH̃)wt − γ(1− ηH̃)wt−1, (26)
where H̃ = diag((x[1])2, (x[2])2). It is important to note that each component of wt evolves independently, due to the fact that H̃ is diagonal.
With := w −w∗, we define for each component j = 1, 2 that
Φ [j] t+1 := E
[( [j] t+1
[j] t
) ⊗ ( [j] t+1
[j] t
)] , j = 1, 2, (27)
where [j] is the j-th component of vector .
The recurrence relation in Eq.26 can be rewritten as
Φ [j] t+1 = B[j]Φ [j] t , (28)
with
B[j] =
E  (1 + γ)2(1− η(x̃[j])2)2 −γ(1 + γ)(1− η(x̃[j])2)2 −γ(1 + γ)(1− η(x̃[j])2)2 γ2(1− η(x̃[j])2)2 (1 + γ)(1− η(x̃[j])2) 0 −γ(1− η(x̃[j])2) 0 (1 + γ)(1− η(x̃[j])2) −γ(1− η(x̃[j])2) 0 0
1 0 0 0

=  (1 + γ)2A[j] −γ(1 + γ)A[j] −γ(1 + γ)A[j] γ2A[j]
(1 + γ)(1− ησ2j ) 0 −γ(1− ησ2j ) 0 (1 + γ)(1− ησ2j ) −γ(1− ησ2j ) 0 0
1 0 0 0  where A[j] := (1− ησ2j )2 + 5(ησ2j )2.
For the ease of analysis, we define u := 1 − γ ∈ (0, 1] and tj := ησ2j , j = 1, 2. Without loss of generality, we assume σ21 = 1 in this section. In this case, t1 = η and t2 = η/κ, where κ is the condition number.
Elementary analysis gives the eigenvalues of B[j]:
λ1 = (1− u)(1− tj),
λ2 = T0 − 21/3
3 T1( T2 + √ T 22 + 4T 3 1
)1/3 + 13 · 21/3 ( T2 + √ T 22 + 4T 3 1 )1/3 ,
λ3 = T0 + (1− i √ 3) 21/3
6 T1( T2 + √ T 22 + 4T 3 1
)1/3 + (1 + i√3) 16 · 21/3 ( T2 + √ T 22 + 4T 3 1 )1/3 ,
λ4 = T0 + (1 + i √ 3) 21/3
6 T1( T2 + √ T 22 + 4T 3 1
)1/3 + (1− i√3) 16 · 21/3 ( T2 + √ T 22 + 4T 3 1 )1/3 ,
where
T0 = 1− u+ u2 3 − 7 3 tj + 7 3 utj − 2 3 u2tj , T1 = −(−3 + 3u− u2 + 7tj − 7utj + 2u2tj)2 + 3(3− 6u+ 4u2 − u3 − 10tj + 20utj − 13u2tj + 3u3tj), T2 = 9u
4 − 9u5 + 2u6 − 72u2tj + 144u3tj − 141u4tj + 69u5tj − 12u6tj + 252t2j − 756ut2j +1185u2t2j − 1110u3t2j + 615u4t2j − 186u5t2j + 24u6t2j − 686t3j + 2058ut3j − 2646u2t3j +1862u3t3j − 756u4t3j + 168u5t3j − 16u6t3j .
Proof idea. For the two-dimensional component decoupled data, we have E [f(wt)]− f(w∗) ≥ σ22E [ ‖wt −w∗‖2 ] = σ22E [ ‖ t‖2 ] . (30)
By definition of Φ in Eq.27, we can see that the convergence rate is lower bounded by the convergence rates of the sequences {‖Φ[j]t ‖}t. By the relation Eq.28, we have that the convergence rate of the sequence {‖Φt‖}t is controlled by the magnitude of the top eigenvalue λmax of B, if Φt has nonzero component along the eigenvector of B with eigenvalue λmax(B). Specifically, if |λmax| > 1, ‖Φ[j]t ‖ grows at a rate of |λmax|t, indicating the divergence of SGD+Nesterov; if |λmax| < 1, then ‖Φ[j]t ‖ converges at a rate of |λmax|t.
In the following, We use the eigen-systems of matrices B[j], especially the top eigenvalue, to analyze the convergence behavior of SGD+Nesterov with any hyper-parameter setting. We show that, for any choice of hyper-parameters (i.e., step-size and momentum parameter), at least one of the following statements must holds:
• B[1] has an eigenvalue larger than 1. • B[2] has an eigenvalue of magnitude 1−O (1/κ).
This is formalized in the following two lemmas.
Lemma 1. For any u ∈ (0, 1], if step size
η > η0(u) := −3− 2u+ 2u2 +
√ 9 + 84u− 164u2 + 100u3 − 20u4 2(9− 15u+ 6u2) , (31)
then, B[1] has an eigenvalue larger than 1.
We will analyze the dependence of the eigenvalues on κ, when κ is large to obtain
Lemma 2. For any u ∈ (0, 1], if step size
η ≤ η0(u) := −3− 2u+ 2u2 +
√ 9 + 84u− 164u2 + 100u3 − 20u4 2(9− 15u+ 6u2) , (32)
then, B[2] has an eigenvalue of magnitude 1−O (1/κ).
Finally, we show that Φt has non-zero component along the eigenvector of B with eigenvalue λmax, hence the convergence of SGD+Nesterov is controlled by the eigenvalue of B[j] with the largest magnitude.
Lemma 3. Assume SGD+Nesterov is initialized with w0 such that both components 〈w0−w∗, e1〉 and 〈w0 − w∗, e2〉 are non-zero. Then, for all t > 2, Φ[j]t has a non-zero component in the eigen direction of B[j] that corresponds to the eigenvalue with largest magnitude. Remark 8. When w is randomly initialized, the conditions 〈w0−w∗, e1〉 6= 0 and 〈w0−w∗, e2〉 6= 0 are satisfied with probability 1, since complementary cases form a lower dimensional manifold which has measure 0.
By combining Lemma 1, 2 and 3, we have that SGD+Nesterov either diverges or converges at a rate of (1−O(1/κ))t, and hence, we conclude the non-acceleration of SGD+Nesterov. In addition, Corollary 1 is a special case of Theorem 1 and is proven by combining Lemma 1 and 3.
In high level, the proof ideas of Lemma 1 and 3 is analogous to those of (9), which proves the non-acceleration of stochastic Heavy Ball method over SGD. But the proof idea of Lemma 2 is unique.
Proof of Lemma 1. The characteristic polynomial of Bj , j = 1, 2, are:
Dj(λ) = λ 4 − (1 + γ)2Ajλ3 + [ 2γ(1 + γ)2(1− ησ2j )Aj − γ2(1− ησ2j )2 − γ2Aj ] λ2
−γ2(1 + γ)2(1− ησ2j )2Ajλ+ γ4(1− ησ2j )2Aj (33)
First note that limλ→∞Dj(λ) = +∞ > 0. In order to show B1 has an eigenvalue larger than 1, it suffices to verity that D1(λ = 1) < 0, because Dj(λ) is continuous.
Replacing γ by 1− u and ησ21 by t1, we have
D1(1) = 4u 2t1 − 2u3t1 − 2ut21 − 10u2t21 + 6u3t21 − 6t31 − 16ut31 + 38u2t31 − 16u3t31 − 18t41
+48ut41 − 42u2t41 + 12u3t41. (34)
Solving for the inequality D1(1) < 0, we have
η = t1 > −3− 2u+ 2u2 +
√ 9 + 84u− 164u2 + 100u3 − 20u4 2(9− 15u+ 6u2) ,
for positive step size η.
Proof of Lemma 2. We will show that at least one of the eigenvalues of B[2] is 1 − O(1/κ), under the condition in Eq.32.
First, we note that t2 = t1/κ = η/κ, which is O(1/κ). We consider the following cases separately: 1) u is Θ(t0.52 ); 2) u is o(t 0.5 2 ) and ω(t2); 3) u is O(t2); 4) u is ω(t 0.5 2 ) and o(1); and 5) u is Θ(1), the last of which includes the case where momentum parameter is a constant.
Note that, for cases 1-4, u is o(1). In such cases, the step size condition Eq.32 can be written as
η < 2
3 u+ o(u). (35)
It is interesting to note that η must be o(1) to not diverge, when u is o(1). This is very different to SGD where a constant step size can result in convergence, see (12).
Case 1: u is Θ(t0.52 ). In this case, the terms u 6, u4t2, u2t22 and t 3 2 are of the same order.
We find that
T0 = 1− u+O(t2), T1 = −3u2 + 12t2 +O(t3/22 ),
T2 = 9u 4 − 72u2t2 + 252t22 +O(t 5/2 2 ).
Hence,
λ1 = 1− u− t2 + ut2
λ2 = 1− u− 21/3
3
12t2 − 3u2
(108(4t2 − u2)3)1/6 +
1
3 · 21/3 (108(4t2 − u2)3)1/6 +O(t2)
λ3 = 1− u+ (1− i √ 3) 21/3
6
12t2 − 3u2
(108(4t2 − u2)3)1/6 + (1 + i
√ 3) 1
6 · 21/3 (108(4t2 − u2)3)1/6 +O(t2)
λ4 = 1− u+ (1 + i √ 3) 21/3
6
12t2 − 3u2
(108(4t2 − u2)3)1/6 + (1− i
√ 3) 1
6 · 21/3 (108(4t2 − u2)3)1/6 +O(t2)
Write t2 = cu2 asymptotically for some constant c. If 4t2 ≥ u2, i.e., 4c− 1 ≥ 0, then
λ2 = 1− u+O(t2), λ3 = 1− u+ √
4c− 1√ 3 u+O(t2),
λ4 = 1− u+ √
4c− 1√ 3 u+O(t2).
If 4t2 ≤ u2, i.e., 4c− 1 ≤ 0, then
λ2 = 1− u+O(t2), λ3 = 1− u+ i √
1− 4c√ 3 u+O(t2),
λ4 = 1− u+ i √
1− 4c√ 3 u+O(t2).
In either case, the first-order term is of order u.
Recall that t2 = η/κ, then we have
u2 = ct2 = cη/κ < 2cu 3κ ≤ cu κ (36)
hence, u < c/κ = O(1/κ). Therefore, λi, i = 1, 2, 3, 4 are 1−O(1/κ). Case 2: u is o(t0.52 ) and ω(t2). In this case,
T0 = 1− u− 7t2/3 + o(t2), T1 = 12t2 + o(t2),
T2 = 252t 2 2 − 72u2t2 + o(t22).
Hence,
λ1 = 1− u− t2 + ut2, λ2 = 1− u+O(t2), λ3 = 1− u+ 2i √ t2 +O(t2), λ4 = 1− u− 2i √ t2 +O(t2).
Assume u ∼ tα2 . When α ∈ (0.5, 1), note that t2 = η/κ ≤ 2u3κ = O(t α 2 /κ), we have t 1−α 2 = O(1/κ), hence t2 = o(1/κ2). This implies that t0.52 = o(1/κ) and u = o(1/κ). Therefore, all the eigenvalues λi are of order 1−O(1/κ). Case 3: u is O(t2). This case if forbidden by the assumption of this lemma. This is because, t2 = η/κ ≤ 2u3κ = o(u) (1/κ is o(1)). This is contradictory to u is O(t2).
Case 4: u is Ω(t0.52 ) and o(1). In this case, we first consider the terms independent of t2, i.e., constant term and u-only terms. These terms can be obtained by putting t2 = 0. In such a setting, the eigenvalues are simplified as:
λ1|t2=0 = 1− u, λ2|t2=0 = 1, λ3|t2=0 = 1− 2u+ u2, λ4|t2=0 = 1− u. (37)
Note that the u-only terms cancel in λ2, so the first order after the constant term must be of t2 (could be t2/u2, t2/u etc.). In the following we are going to analyze the t2 terms.
Since u is Ω(t0.52 ), u 2 has lower order than t2, and t2/u2 is o(1). This allows us to do Taylor expansion:
T1 = −3u2(1− 4 t2 u2 ) + f(u) + higher order terms, T2 = 9u 4(1− 8 t2
u2 ) + g(u) + higher order terms,
where f(u) and g(u) are u-terms only, which, by the above analysis in Eq.37, are shown to contribute nothing to λ2. Hence, we use the first terms of T1 and T2 above to analyze the first order term of λ2. Plugging in these term to the expression of λ2, and keeping the lowest order of t2, we find a zero coefficient of the lowest order t2-term: t2/u2.
Hence, λ2 can be written as:
λ2 = 1− c t2 u +O(t2), (38)
where c is the coefficient.
On the other hand, by Eq.35 and t2 = η/κ, we have
t2 u = O(1/κ). (39)
Therefore, we can write Eq.38 as λ2 = 1−O(1/κ). Case 5: u is O(1). This is the case where the momentum parameter is κ-independent. Using the same argument as in case 4, we have zero u-only terms. Then, directly taking Taylor expansion with respect to t2 results in:
λ2 = 1−O(t2) = 1−O(1/κ). (40)
Proof of Lemma 3. This proof follows the idea of the proof for stochastic Heavy Ball method in (9). The idea is to examine the subspace spanned by Φ[j]t , t = 0, 1, 2, · · · , j = 1, 2, and to prove that the eigenvector of B[j] corresponding to top eigenvalue (i.e., eigenvalue with largest magnitude) is not orthogonal to this spanned subspace. This in turn implies that there exists a non-zero component of Φ
[j] t in the eigen direction of B[j] corresponding to top eigenvalue, and this decays/grows at a rate of
λtmax(B[j]).
Recall that Φ[j]t+1 = B[j]Φ [j] t . Hence, if Φ [j] t has non-zero component in the eigen direction of B[j] with top eigenvalue, then Φ[j]t should also have non-zero component in the same direction. Thus, it suffices to show that at least one of Φ[j]0 ,Φ [j] 1 ,Φ [j] 2 and Φ [j] 3 has non-zero component in the eigen direction with top eigenvalue.
Since H̃ is diagonal for this two-dimensional decoupled data, w[1] and w[2] evolves independently, and we can analyze each component separately. In addition, it can be seen that each of the initial values w[j]0 − (w∗)[j], which is non-zero by the assumption of this lemma, just acts as a scale factor during the training. Hence, without loss of generality, we can assume w[j]0 − (w∗)[j] = 1, for each j. Then, according to the recurrence relation of SGD+Nesterov in Eq.26,(
[j] 0 [j] −1
) = ( 1 1 ) , ( [j] 1
[j] 0
) = ( s [j] 1 1 ) , ( [j] 2
[j] 1
) = ( (1 + γ)s [j] 1 s [j] 2 − γs [j] 2
s [j] 1 ) ,(
[j] 3 [j] 2
) = ( s [j] 3 [ (1 + γ)2s [j] 1 s [j] 2 − γ(1 + γ)s [j] 2 ] − γs[j]3 [ (1 + γ)s [j] 1 s [j] 2 − γs [j] 2 ] (1 + γ)s
[j] 1 s [j] 2 − γs [j] 2
) , (41)
where s[j]t = 1− η(x̃ [j] t ) 2, t = 1, 2, 3, · · · .
Denote the vectorized form of Φ[j]t as vec(Φ [j] t ), which is a 4 × 1 column vector. We stack the vectorized forms of Φ[j]0 ,Φ [j] 1 ,Φ [j] 2 ,Φ [j] 3 to make a 4× 4 matrix, denoted asM[j]:
M[j] = [vec(Φ[j]0 ) vec(Φ [j] 1 ) vec(Φ [j] 2 ) vec(Φ [j] 3 )]. (42)
Note that Φ[j]t , t = 0, 1, 2, 3, are symmetric tensors, which implies thatM[j] contains two identical rows. Specifically, the second and third row of M[j] are identical. Therefore, the vector v = (0,−1/ √ 2, 1/ √
2, 0)T is an eigenvector ofM[j] with eigenvalue 0. In fact, v is also an eigenvector of B[j] with eigenvalue γ(1− ησ2j ). Note that
det(B[j]) = γ4(1− ησ2j )2 ( (1− ησ2j )2 + 5(ησ2j )2 ) ≥ γ4(1− ησ2j )4.
Hence, v is not the eigenvector along top eigenvalue, and therefore, is orthogonal to the eigen space with top eigenvalue.
In order to prove at least one of Φ[j]t , t = 0, 1, 2, 3, has a non-zero component along the eigen direction of top eigenvalue, it suffices to verify thatM[j] is rank 3, i.e., spans a three-dimensional space. Equivalently, we consider the following matrix
M′ := E  20,1 21,1 22,1 0,1 −1,1 0,1 1,1 1,1 2,1 2−1,1 2 0,1 2 1,1  , (43) where we omitted the superscript [j] for simplicity of the expression. If the determinant ofM′[j] is not zero, then it is full rank, and henceM[j] spans a three-dimensional space. Plug in the expressions in Eq.41, then we have
det(M′[j]) = (1− u)t3j
2 (36t2j − 18ut2j + 6utj + 3tj − 3u+ 1), (44)
where tj = ησ2j and u = 1 − γ. det(M′[j]) = 0 if and only if the polynomial 36t2j − 18ut2j + 6utj + 3tj − 3u+ 1 = 0. Solving for this equation, we have
tj = −1− 2u±
√ (1 + 2u)2 + 8(3− u)(u− 2)
12(2− u) . (45)
We note that, for all u ∈ [0, 1) both tj are not positive. This means that, for all u and positive tj , the determinant det(M′[j]) can never be zero.
Therefore, for each j = 1, 2,M′[j] is full rank, andM[j] spans a three-dimensional space, which includes the eigenvector with the top eigenvalue of B[j]. Hence, at least one of Φ[j]t , t ∈ {0, 1, 2, 3} has non-zero component in the eigen direction with top eigenvalue. By Φ[j]t+1 = B[j]Φ [j] t , all succeeding Φ [j] t also have non-zero component in the eigen direction with top eigenvalue of B[j].
Proof of Theorem 1. Lemma 1 and 2 show that, for any hyper-parameter setting (η, γ) with η > 0 and γ ∈ (0, 1), either top eigenvalue of B[1] is larger than 1 or top eigenvalue of B[2] is 1−O(1/κ). Hence, |λmax| is either greater than 1 or is 1 − O(1/κ). Lemma 3 shows that Φt has non-zero component along the eigenvector of B with eigenvalue λmax(B). By Eq.28 and Lemma 1 and 2, the sequence {‖Φt‖}t either diverges or converges at a rate of (1 − O(1/κ)t. By definition of Φ in Eq.27, we have that {‖ t‖}t either diverges or converges at a rate of (1−O(1/κ))t. Note that, for the two-dimensional component decoupled data, we have
E [f(wt)]− f(w∗) ≥ σ22E [ ‖wt −w∗‖2 ] = σ22E [ ‖ t‖2 ] . (46)
Therefore, the convergence rate of SGD+Nesterov is lower bounded by (1−O(1/κ))t. Note that the convergence rate of SGD on this data is (1 − O(1/κ))t, hence SGD+Nesterov does not accelerate SGD on the two-dimensional component decoupled dataset.
Proof of Corollary 1. When η = 1/L1 and γ ∈ [0.6, 1], the condition in Eq.31 is satisfied. By Lemma 1, the top eigenvalue λ[1]max of B[1] is larger than 1. By Lemma 3, Φ[1] has non-zero component along the eigenvector with this top eigenvalue λ[1]max. Hence, |〈wt−w∗, e1〉| grows at a rate of (λ [1] max)t.

heading: D PROOF OF THEOREM 2
text: We first give a lemma that is useful for dealing with the mini-batch scenario:
Lemma 4. If square loss f is in the interpolated setting, i.e., there exists w∗ such that f(w∗) = 0, then E[H̃mH−1H̃m]−H 1m (κ̃− 1)H.

heading: Proof.
text: E[H̃mH−1H̃m] = 1 m2 E [ m∑ i=1 x̃ix̃i TH−1x̃ix̃i T + ∑ i6=j x̃ix̃i TH−1x̃jx̃j T ] 1 m κ̃H + m− 1 m H.
Proof of Theorem 2. By Eq.8, we have
‖vt+1 −w∗‖2H−1 = ‖(1− α)vt + αut −w ∗‖2H−1︸ ︷︷ ︸
Ã
+ δ2‖H̃m(ut −w∗)‖2H−1︸ ︷︷ ︸ B̃
−2δ〈H̃m(ut −w∗), (1− α)vt + αut −w∗〉H−1︸ ︷︷ ︸ C̃
Using the convexity of the norm ‖ · ‖H−1 and the fact that µ is the smallest non-zero eigenvalue of the Hessian, we get
E[Ã] ≤ (1− α)‖vt −w∗‖2H−1 + α‖ut −w ∗‖2H−1 ≤ (1− α)‖vt −w ∗‖2H−1 + α µ ‖ut −w∗‖2.
Applying Lemma 4 on the term B̃, we have
E[B̃] ≤ δ2κ̃m‖ut −w∗‖2H . (47)
Note that E[H̃m] = H , then
E[C̃] = −2δ〈ut −w∗, (1− α)vt + αut −w∗〉
= −2δ〈ut −w∗,ut −w∗ + 1− α α (ut −wt)〉 = −2δ‖ut −w∗‖2 + 1− α α δ ( ‖wt −w∗‖2 − ‖ut −w∗‖2 − ‖wt − ut‖2
) ≤ 1− α
α δ‖wt −w∗‖2 −
1 + α
α δ‖ut −w∗‖2. (48)
Therefore, E [ ‖vt+1 −w∗‖2H−1 ] ≤ (1−α)‖vt−w∗‖2H−1+ 1− α α δ‖wt−w∗‖2+ ( α µ − 1 + α α δ ) ‖ut−w∗‖2+δ2κ̃m‖ut−w∗‖2H .
On the other hand, by the fact that E[H̃2m] 1mE[H̃ 2 1 ] + m−1 m H 2, (see (12)),
E[‖wt+1 −w∗‖2] = E [ ‖ut −w∗ − ηH̃m(ut −w∗)‖2 ] ≤ ‖ut −w∗‖2 − 2η‖ut −w∗‖2H + η2Lm‖ut −w∗‖2H . (49)
Hence,
E [ ‖vt+1 −w∗‖2H−1 + δ
α ‖wt+1 −w∗‖2 ] ≤ (1− α)‖vt −w∗‖2H−1 + 1− α α
δ‖wt −w∗‖2 + (α/µ− δ)︸ ︷︷ ︸ c1 ‖ut −w∗‖2
+ ( δ2κ̃m + δη 2Lm/α− 2ηδ/α )︸ ︷︷ ︸
c2
‖ut −w∗‖2H .
By the condition in Eq.10, c1 ≤ 0, c2 ≤ 0, then the last two terms are non-positive. Hence, E [ ‖vt+1 − x∗‖2H−1 + δ
α ‖wt+1 −w∗‖2
] ≤ (1− α) ( ‖vt −w∗‖2H−1 + δ
α ‖wt −w∗‖2 ) ≤ (1− α)t+1 ( ‖v0 −w∗‖2H−1 + δ
α ‖w0 −w∗‖2
) .

heading: E PROOF OF THEOREM 4
text: Proof. The update rule for variable v is, as in Eq.8:
vt+1 = (1− α)vt + αut − δ∇̃mf(ut). (50)
By 1/µ-strong convexity of g, we have
g(vt+1) ≤ g ((1− α)vt + αut) + 〈∇g ((1− α)vt + αut) ,−δ∇̃mf(ut)〉+ δ2 1
2µ ‖∇̃mf(ut)‖2.
(51) Taking expectation on both sides, we get
E[g(vt+1)] = g ((1− α)vt + αut) + 〈∇g ((1− α)vt + αut) ,−δ∇f(ut)〉+ δ2 1
2µ E[‖∇̃mf(ut)‖2]
≤ (1− α)g(vt) + αg(ut)− δ(1− )〈(1− α)vt + αut −w∗,ut −w∗〉+ δ2 2Lm 2µ f(ut),
where in the last inequality we used the convexity of g and the assumption in the Theorem.
By Eq.48,
−δ(1− )〈(1−α)vt+αut−w∗,ut−w∗〉 ≤ δ(1− )
2
( 1− α α ‖wt −w∗‖2 − 1 + α α ‖ut −w∗‖2 ) .
By the strong convexity of g,
αg(ut) ≤ α
2µ ‖ut −w∗‖2.
Hence,
E[g(vt+1)] ≤ (1−α)g(vt)+(1−α) δ(1− )
2α ‖wt−w∗‖2+
( α
2µ − δ(1 + α)(1− ) 2α
) ‖ut−w∗‖2+δ2κmf(ut).
(52) On the other hand,
E[‖wt+1 −w∗‖2] = ‖ut −w∗‖2 − 2η〈ut −w∗,∇f(ut)〉+ η2E[‖∇̃mf(ut)‖2] ≤ ‖ut −w∗‖2 − 2ηf(ut) + η2 · 2Lmf(ut), (53)
where in the last inequality we used the convexity of f .
Multiply a factor of δ(1− )2α onto Eq.53 and add it to Eq.52, then
E [ g(vt+1) +
δ(1− ) 2α
‖wt+1 −w∗‖2 ] ≤ (1− α)g(vt) + (1− α)
δ(1− ) 2α ‖wt −w∗‖2
+ 1
2
( α
µ − δ(1− )
) ‖ut −w∗‖2
+
( δ
α (1− )(η2Lm − η) + δ2κm
) f(ut). (54)
If the hyper-parameters are selected as:
η = 1
2Lm , α = 1− 2κm , δ = 1 2Lm , (55)
then the last two terms in Eq.54 are non-positive. Hence, E [ g(vt+1) +
δ(1− ) 2α
‖wt+1 −w∗‖2 ] ≤ (1− α)E [ g(vt) +
δ(1− ) 2α
‖wt −w∗‖2 ] , (56)
which implies E [ g(vt) +
δ(1− ) 2α
‖wt −w∗‖2 ] ≤ (1− 1− 2κm )tE [ g(v0) + δ(1− ) 2α ‖w0 −w∗‖2 ] . (57)
Since f(wt) ≤ L/2 · ‖wt −w∗‖2 (by smoothness), then we have the final conclusion
E[f(wt)] ≤ C · (1− 1− 2κm )t, (58)
with C being a constant.

heading: F EVALUATION ON SYNTHETIC DATA

heading: F.1 DISCUSSION ON SYNTHETIC DATASETS
text: Component Decoupled Data. The data is defined as follows (also defined in section 3):
Fix an arbitrary w∗ ∈ R2 and let z be randomly drawn from the zero-mean Gaussian distribution with variance E[z2] = 2, i.e. z ∼ N (0, 2). The data points (x, y) ∈ D are constructed as follow:
x = { σ1 · z · e1 w.p. 0.5 σ2 · z · e2 w.p. 0.5, and y = 〈w ∗,x〉, (59)
where e1, e2 ∈ R2 are canonical basis vectors, σ1 > σ2 > 0. Note that the corresponding square loss function on D is in the interpolation regime, since f(w∗) = 0. The Hessian and stochastic Hessian matrices turn out to be
H = [ σ21 0 0 σ22 ] , H̃ = [ (x[1])2 0 0 (x[2])2 ] . (60)
Note that H̃ is diagonal, which implies that stochastic gradient based algorithms applied on this data evolve independently in each coordinate. This allows a simplified directional analysis of the algorithms applied.
Here we list some useful results for our analysis. The fourth-moment of Gaussian variable E[z4] = 3E[z2]2 = 12. Hence, E[(x[j])2] = σ2j and E[(x[j])4] = 6σ4j , where superscript j = 1, 2 is the index for coordinates in R2. It is easy to check that κ1 = 6σ21/σ22 and κ̃ = 6.
Gaussian Data. Suppose the data feature vectors {xi} are zero-mean Gaussian distributed, and yi = 〈w∗,xi〉, ∀i, where w∗ is fixed but unknown. Then, by the fact that E[z1z2z3z4] = E[z1z2]E[z3z4] + E[z1z3]E[z2z4] + E[z1z4]E[z2z3] for zero-mean Gaussian random variables z1, z2, z3 and z4, we have
E[H̃H̃] = (2H + tr(H))H, E[H̃H−1H̃] = E[x̃x̃TH−1x̃x̃T ] = (2 + d)H,
where d is the dimension of the feature vectors. Hence L1 = 2λmax(H) + tr(H) and κ̃ = 2 + d, and κ1 = L1/µ. This implies a convergence rate of O(e−t/ √ (2+d)κ1) of MaSS when batch size is 1. Particularly, if the feature vectors are n-dimensional, e.g., as in kernel learning, then MaSS with mini batches of size 1 has a convergence rate of O(e−t/ √ nκ1).

heading: F.2 EVALUATION OF FAST CONVERGENCE OF MASS AND NON-ACCELERATION OF SGD+NESTEROV
text: In this subsection, we show additional empirical verification for the fast convergence of MaSS, as well as the non-acceleration of SGD+Nesterov, on synthetic data. In addition, we show the divergence of SGD+Nesterov when using the same step size as SGD and MaSS, as indicated by Corollary 1.
We consider two families of synthetic datasets:
• Component decoupled: (as defined in Section 3). Fix an arbitrary w∗ ∈ R2 with all components non-zero. xi is drawn fromN (0, diag(2σ21 , 0)) orN (0, diag(0, 2σ22)) with probability 0.5 each. yi = 〈w∗,xi〉 for all i. • 3-d Gaussian: Fix an arbitrary w∗ ∈ R3 with all components non-zero. xi are independently drawn from N (0, diag(σ21 , σ21 , σ22)), and yi = 〈w∗,xi〉 for all i.
Non-acceleration of SGD+Nesterov and accelerated convergence of MaSS. We compare MaSS with SGD and SGD+Nesterov on linear regression with the above datasets. Each comparison is
performed on either 3-d Gaussian or component decoupled data with fixed σ1 and σ2. For each setting of (σ1, σ2), we randomly select w∗, and generate 2000 samples for the dataset. Batch sizes for all algorithms are set to be 1. We report the performances of SGD, SGD+Nesterov and SGD+HB using their best hyper-parameter setting selected from dense grid search. On the other hand, we do not tune hyper-parameters of MaSS, but use the suggested setting by our theoretical analysis, Eq. 14. Specifically, we use
Component decoupled: η∗1 = 1
6 , η∗2 =
5
36 + 6σ2 , γ∗ = 6− σ2 6 + σ2 ; (61a)
3-d Gaussian: η∗1 = 1
4 , η∗2 =
1
5
√ 20√
20 + σ2 , γ∗ = √ 20− σ2√ 20 + σ2 . (61b)
For ASGD, we use the setting suggested by (8).
Figure 6 (in addition to Fig 4) and Figure 7 show the curves of the compared algorithms under various data settings. We observe that: 1) SGD+Nesterov with its best hyper-parameters is almost identical to the optimal SGD; 2) MaSS, with the suggested hyper-parameter selections, converges faster than all of the other algorithms, especially SGD. These observations are consistent with our theoretical results about non-acceleration of SGD+Nesterov, Theorem 1, and accelerated convergence of MaSS, Theorem 3.
Recall that MaSS differs from SGD+Nesterov by only a compensation term, this experiment illustrates the importance of this term. Note that the vertical axis is log scaled. Then the linear decrease of log losses in the plots implies an exponential loss decrease, and the slopes correspond to the coefficients in the exponents.
Divergence of SGD+Nesterov with large step size. As discussed in Corollary 1, SGD+Nesterov diverges with step size η∗ = 1/L1 (when γ ∈ [0.6, 1]), which is the optimal choice of step size for both SGD and MaSS. We run SGD+Nesterov, with step size η∗ = 1/L1, to optimize the square loss function on component decoupled data mentioned above. Figure 8 shows the divergence of SGD+Nesterov with two common choices of momentum parameter γ: 0.9 and 0.99.

heading: F.3 COMPARISON OF VASWANI’S METHOD (21) WITH SGD AND MASS ON GAUSSIAN-DISTRIBUTED DATA
text: The analysis in Vaswani et. al. (21) is based on the strong growth condition (SGC). Assuming SGC with the parameter ρ on the loss function they prove convergence rate ( 1− √ 1/ρ2κ )t
of their method (called SGD with Nesterov acceleration in their paper), where t is the iteration number. In the following, we show that, on a simple (zero-mean) Gaussian distributed data, this rate is slower than that of SGD, which has a rate of (1− 1/κ)t. On the other hand, MaSS achieves the accelerated rate ( 1− 1/ √ (2 + d)κ )t .
Consider the problem of minimizing the squared loss, f(w) = ∑ i fi(w) = 1 2 ∑ i(w
Txi − yi)2, over a zero-mean Gaussian distributed dataset, as defined in F.1. Then,
Ei‖∇fi(w)‖2 = (w −w∗)TEi [ xix T i xix T i ] (w −w∗)
= (w −w∗)T (2H2 + tr(H)H)(w −w∗), (62)
whereas,
‖∇f(w)‖2 = (w −w∗)TH2(w −w∗). (63)
According to the definition of SGC,
Ei‖∇fi(w)‖2 ≤ ρ‖∇f(w)‖2. (64)
Hence the SGC parameter ρ must satisfy
ρ ≥ 2 + tr(H) λmin(H) > κ, (65)
where κ is the condition number. In this case, the convergence rate ( 1− √ 1/ρ2κ )t
of Vaswani’s method would be slower than( 1− 1/κ3/2 )t , which is slower than SGD.
On the other hand, MaSS accelerates over SGD on this dataset. Recall from Section F.1 that, for this Gaussian data, κ̃ = 2 + d, where d is the dimension of the the data. According to Theorem 3, the
convergence rate of MaSS is ( 1− 1/ √ (2 + d)κ )t , which is faster than that of SGD.

heading: G MINI-BATCH DEPENDENCE OF CONVERGENCE SPEED AND OPTIMAL HYPER-PARAMETERS

heading: G.1 ANALYSIS
text: The two critical points are defined as follow:
m∗1 = min(L1/L, κ̃),m ∗ 2 = max(L1/L, κ̃).
When m < m∗1, we have m < L1/L and m < κ̃, which implies Lm ≈ L1/m, κm ≈ κ1/m and κ̃m ≈ κ̃/m. Plugging into Eq.13, we find that the optimal selection of hyper-parameters is approximated by:
η∗(m) ≈ m · η∗(1), α∗(m) ≈ m · α∗(1), δ∗(m) ≈ m · δ∗(1),
and the convergence rate in Eq.16 is approximately O(e−m·t/ √ κ1κ̃), the latter of which indicates a linear gains in convergence when m increases.
Now consider the following two cases: (i) For m ≥ L1/L, we have κm ≤ 2κ following the definition of Lm in Eq.2; (ii) For m ≥ κ̃, κ̃m = κ̃/m + (m − 1)/m ≤ 2. When either case holds m ≥ m∗1, the convergence rate (Eq.16) becomes
O ( max { e − √ m·t√ 2κ(κ̃+m−1) , e − √ m·t√ 2κ1+2(m−1)κ }) ,
indicating that increasing the mini-batch size results in sublinear gain in the convergence speed, of at most √ m. Moreover, in the saturation regime m > m∗2 (i.e., when both conditions (i) and (ii) hold), the convergence rate is then upper-bounded by O(exp(− t 2 √ κ
)). That implies that a single iteration of MaSS with mini-batch size m is equivalent up to a multiplicative factor of 2 to an iteration with full gradient.

heading: G.2 EMPIRICAL VERIFICATION
text: We empirically verify the three-regime partition observed in section 5 using zero-mean Gaussian data.
In this evaluation, we set the covariant matrix of the (zero-mean) Gaussian to be:
Σ = diag(1, · · · , 1︸ ︷︷ ︸ 8 , 2−10, · · · , 2−10︸ ︷︷ ︸ 40 ). (66)
and generate 2000 samples following the Gaussian distribution N (0,Σ). Recall from Example 1 that, for zero-mean Gaussian data, L1 = 2λmax(H)+tr(H) and κ̃ = 2+d. In this case, L1/L ≈ 10 and κ̃ ≈ 50. The two critical values m∗1,m∗2 of mini-batch size are then:
m∗1 ≈ 10, m∗2 ≈ 50.
In the experiments, we run MaSS with a variaty of mini-batch size m, ranging from 1 to 160, on this Gaussian dataset. For each training process, we compute the convergence speed s(m), which is defined to be the inverse of the number of iterations needed to achieve a training error of ε. Running 1 iteration of MaSS with minibatch of sizem is almost as effective as running s(m)/s(1) iterations of MaSS with mini-batch of size 1.
Figure 9 demonstrates the convergence speed s(m) as a function of the mini-batch size m. We see that the three regimes defined by m∗1 and m∗2 coincide with the analysis in subsection 5: left (m < m∗1), almost linear; middle (m ∈ [m∗1,m∗2]), sublinear; right (m > m∗2), saturation.

heading: H EXPERIMENTAL SETUP

heading: H.1 NEURAL NETWORK ARCHITECTURES
text: Fully-connected Network. The fully-connected neural network has 3 hidden layers, with 100 ReLU-activated neurons in each layer. After each hidden layer, there is a dropout layer with keep probability 0.5. This network takes 784-dimensional vectors as input, and has 10 softmax-activated output neurons. It has ≈99k trainable parameters in total. Convolutional Neural Network (CNN). The CNN we considered has three convolutional layers with kernel size of 5 × 5 and without padding. The first two convolutional layers have 64 channels each, while the last one has 128 channels. Each convolutional layer is followed by a 2 × 2 max pooling layer with stride of 2. On top of the last max pooling layer, there is a fully-connected ReLU-activated layer of size 128 followed by the output layer of size 10 with softmax non-linearity. A dropout layer with keep probability 0.5 is applied after the full-connected layer. The CNN has ≈576k trainable parameters in total. Residual Network (ResNet). We train a ResNet (7) with 32 convolutional layers. The ResNet-32 has a sequence of 15 residual blocks: the first 5 blocks have an output of shape 32 × 32 × 16, the following 5 blocks have an output of shape 16×16×32 and the last 5 blocks have an output of shape 8×8×64. On top of these blocks, there is a 2×2 average pooling layer with stride of 2, followed by a output layer of size 10 with softmax non-linearity. The ResNet-32 has≈467k trainable parameters in total.
We use the fully-connected network to classify the MNIST dataset, and use CNN and ResNet to classify the CIFAR-10 dataset.

heading: H.2 EXPERIMENTAL SETUP FOR TEST PERFORMANCE OF MASS
text: The 3-layer CNN we use for this experiment is slightly different from the aforementioned one. We insert a batch normalization (BN) layer after each convolution computation, and remove the dropout layer in the fully-connected phase.
Reduction of Learning Rate. On both CNN and ResNet-32, we initialize the learning rate of SGD, SGD+Nesterov and Mass using the same value, while that of Adam is set to 0.001, which is common default setting. On CNN, the learning rates for all algorithms drop by 0.1 at 60 and 120 epochs, and we train for total 150 epochs. On ResNet-32, the learning rates drop by 0.1 at 150 and 225 epochs, and we train for total 300 epochs.
Whenever the learning rate of MaSS reduced, we restart the MaSS with the latest learned weights. The reason for restarting is to avoid the mismatch between the large momentum term and small gradient descent update after the reduction of learning rate.
Data Augmentation. We augment the Cifar-10 training data by enabling random horizontal flip, random horizontal shift and random vertical shift of the images. For the random shift of images, we allow a shift range of 10% of the image width or height.
Hyper-parameter Selection for MaSS. Since the reduction of learning rate may affect the suggested value of δ (c.f. Eq.13), we consider (η, α, κ̃m) as independent hyper-parameters instead, where κ̃m = η/(αδ). Observing the fact that κ̃m ≤ κm, we select κ̃m in the range (1, 1/α). In our experiments, we set α = 0.05, corresponding to the “momentum parameter” γ being 0.90, which is commonly used in Nesterov’s method. We find that good choices of κ̃m often locate in the interval [2, 25] for mini-batch size m = 64 in our experiments.
The classification results of MaSS in Table 6 are obtained under the following hyper-parameter settings:
• CNN: η = 0.01 (initial), α = 0.05, κ̃m = 3; η = 0.3 (initial), α = 0.05, κ̃m = 6.
• ResNet-32: η = 0.1 (initial), α = 0.05, κ̃m = 2; η = 0.3 (initial), α = 0.05, κ̃m = 24.

Review summaries:
The authors present a new first order optimization method that adds a corrective term to Nesterov SGD. They demonstrate that this adjustment is necessary and sufficient to benefit from the faster convergence of Nesterov gradient descent in the stochastic case. In the full-batch (non deterministic) setting, their algorithm boils down to the classical formulation of Nesterov GD. Their approach is justified by a well conducted theoretical analysis and some empirical work on toy datasets.