EIGENCURVE: OPTIMAL LEARNING RATE SCHEDULE
Year: 2021

Abstract:
Learning rate schedulers have been widely adopted in training deep neural networks. Despite their practical importance, there is a discrepancy between its practice and its theoretical analysis. For instance, it is not known what schedules of SGD achieve best convergence, even for simple problems such as optimizing quadratic objectives. In this paper, we propose Eigencurve, the first family of learning rate schedules that can achieve minimax optimal convergence rates (up to a constant) for SGD on quadratic objectives when the eigenvalue distribution of the underlying Hessian matrix is skewed. The condition is quite common in practice. Experimental results show that Eigencurve can significantly outperform step decay in image classification tasks on CIFAR-10, especially when the number of epochs is small. Moreover, the theory inspires two simple learning rate schedulers for practical applications that can approximate eigencurve. For some problems, the optimal shape of the proposed schedulers resembles that of cosine decay, which sheds light to the success of cosine decay for such situations. For other situations, the proposed schedulers are superior to cosine decay.

heading: 1 INTRODUCTION
text: Many machine learning models can be represented as the following optimization problem:
min w f(w) ,
1
n n∑ i=1 fi(w), (1.1)
such as logistic regression, deep neural networks. To solve the above problem, stochastic gradient descent (SGD) (Robbins & Monro, 1951) has been widely adopted due to its computation efficiency in large-scale learning problems (Bottou & Bousquet, 2008), especially for training deep neural networks.
Given the popularity of SGD in this field, different learning rate schedules have been proposed to further improve its convergence rates. Among them, the most famous and widely used ones are inverse time decay, step decay (Goffin, 1977), and cosine scheduler (Loshchilov & Hutter, 2017). The learning rates generated by the inverse time decay scheduler depends on the current iteration number inversely. Such a scheduling strategy comes from the theory of SGD on strongly convex functions, and is extended to non-convex objectives like neural networks while still achieving good performance. Step decay scheduler keeps the learning rate piecewise constant and decreases it by a factor after a given amount of epochs. It is theoretically proved in Ge et al. (2019) that when the objective is quadratic, step decay scheduler outperforms inverse time decay. Empirical results are also provided in the same work to demonstrate the better convergence property of step decay in training neural networks when compared with inverse time decay. However, even step decay is proved to be near optimal on quadratic objectives, it is not truly optimal. There still exists a log T gap away from the minimax optimal convergence rate, which turns out to be non-trivial in a wide range of settings and may greatly impact step decay’s empirical performance. Cosine decay scheduler (Loshchilov & Hutter, 2017) generates cosine-like learning rates in the range [0, T ], with T being the maximum iteration. It is a heuristic scheduling strategy which relies on the observation that good performance in practice can be achieved via slowly decreasing the learning rate in the
beginning and “refining” the solution in the end with a very small learning rate. Its convergence property on smooth non-convex functions has been shown in Li et al. (2021), but the provided bound is still not tight enough to explain its success in practice.
Except cosine decay scheduler, all aforementioned learning rate schedulers have (or will have) a tight convergence bound on quadratic objectives. In fact, studying their convergence property on quadratic objective functions is quite important for understanding their behaviors in general nonconvex problems. Recent studies in Neural Tangent Kernel (NTK) (Arora et al., 2019; Jacot et al., 2020) suggest that when neural networks are sufficiently wide, the gradient descent dynamic of neural networks can be approximated by NTK. In particular, when the loss function is least-square loss, neural network’s inference is equivalent to kernel ridge regression with respect to the NTK in expectation. In other words, for regression tasks, the non-convex objective in neural networks resembles quadratic objectives when the network is wide enough.
The existence of log T gap in step decay’s convergence analysis, which will be proven to be tight in a wide range of settings, implies that there is still room for improvement in theory. Meanwhile, the existence of cosine decay scheduler, which has no strong theoretical convergence guarantees but possesses good empirical performance in certain tasks, suggests that its convergence rate may depend on some specific properties of the objective determined by the network and dataset in practice. Hence it is natural to ask what those key properties may be, and whether it is possible to find theoretically-optimal schedulers whose empirical performance are comparable to cosine decay if those properties are available. In this paper, we offer an answer to these questions. We first derive a novel eigenvalue-distribution-based learning rate scheduler called eigencurve for the quadratic function. Combining with the different eigenvalue distributions of different types of networks, new neural-network-based learning rate schedulers can be generated based on our proposed paradigm, which achieve better convergence properties than step decay in Ge et al. (2019). Specifically, eigencurve closes the log T gap in step decay and reaches minimax optimal convergence rates if the Hessian spectrum is skewed. We summarize the main contributions of this paper as follows.
1. To the best of our knowledge, this is the first work that incorporates the eigenvalue distribution of objective function’s Hessian matrix into designing learning rate schedulers. Accordingly, based on the eigenvalue distribution of the Hessian, we propose a novel eigenvalue distribution based learning rate scheduler named eigencurve.
2. Theoretically, eigencurve can achieve optimal convergence rate (up to a constant) for SGD on quadratic objectives when the eigenvalue distribution of the Hessian is skewed. Furthermore, even when the Hessian is not skewed, eigencurve can still achieve no worse convergence rate than the step decay schedule in Ge et al. (2019), whose convergence rate are proven to be sub-optimal in a wide range of settings.
3. Empirically, on image classification tasks, eigencurve achieves optimal convergence rate for several models on CIFAR-10 and ImageNet if the loss can be approximated by quadratic objectives. Moreover, it obtains much better performance than step decay on CIFAR-10, especially when the number of epochs is small.
4. Intuitively, our learning rate scheduler sheds light on the theoretical property of cosine decay and provides a perspective of understanding the reason why it can achieve good performance on image classification tasks. The same idea has been used to inspire and discover several simple families of schedules that works in practice.
Problem Setup For the theoretical analysis and the aim to derive our eigenvalue-dependent learning rate schedulers, we mainly focus on the quadratic function, that is,
min w f(w) , Eξ [f(w, ξ)] , where f(w, ξ) =
1 2 w>H(ξ)w − b(ξ)>w, (1.2)
where ξ denotes the data sample. Hence, the Hessian of f(w) is
H = Eξ [H(ξ)] . (1.3)
Letting us denote b = Eξ[b(ξ)], we can obtain the optima of problem (1.2)
w∗ = H −1b. (1.4)
Given an initial iterate w0 and the learning rate sequence {ηt}, the stochastic gradient update is
wt+1 = wt − ηt∇f(wt, ξ) = wt − ηt(H(ξ)wt − b(ξ)). (1.5)
We denote that
nt = Hwt − b− (H(ξ)wt − b(ξ)) , µ , λmin(H), L , λmax(H), and, κ , L/µ. (1.6)
In this paper, we assume that
Eξ [ ntn > t ] σ2H. (1.7)
The reason for this assumption is presented in Appendix G.5.
Related Work In convergence analysis, one key property that separates SGD from vanilla gradient descent is that in SGD, noise in gradients dominates. In gradient descent (GD), constant learning rate can achieve linear convergence O(cT ) with 0 < c < 1 for strongly convex objectives, i.e. obtaining f(w(t)) − f(w∗) ≤ in O(log( 1 )) iterations. However, in SGD, f(w(t)) cannot even be guaranteed to converge to f(w∗) due to the existence of gradient noise (Bottou et al., 2018). Intuitively, this noise leads to a variance proportional to the learning rate size, so constant learning rate will always introduce a Ω(ηt) = Ω(η0) gap when compared with the con-
vergence rate of GD. Fortunately, inverse time decay scheduler solves the problem by decaying the learning rate inversely proportional to the iteration number t, which achieves O( 1T ) convergence rate for strongly convex objectives, specifically, O(dσ 2
T · κ). However, this is sub-optimal since the minimax optimal rate for SGD is O(dσ 2
T ) (Ge et al., 2019; Jain et al., 2018). Moreover, in practice, κ can be very big for large neural networks, which makes inverse time decay scheduler undesirable for those models. This is when step decay (Goffin, 1977) comes to play. Empirically, it is widely adopted in tasks such as image classification and serves as a baseline for a lot of models. Theoretically, it has been proven that step decay can achieve nearly optimal convergence rateO(dσ 2
T · log T ) for strongly convex least square regression (Ge et al., 2019). Nevertheless, step decay is not always the best choice for image classification tasks. In practice, cosine decay (Loshchilov & Hutter, 2017) can achieve comparable or even better performance, but the reason behind this superior performance is still unknown in theory (Gotmare et al., 2018). All the aforementioned results are summarized in Table 1, along with our results in this paper. It is worth mentioning that the minimax optimal rate O(dσ 2
T ) can be achieved by iterate averaging methods (Jain et al., 2018; Bach & Moulines, 2013; Défossez & Bach, 2015; Frostig et al., 2015; Jain et al., 2016; Neu & Rosasco, 2018), but it is not a common practice to use them in training deep neural networks, so only the final iterate (Shamir, 2012) behavior of SGD is analyzed in this paper, i.e. the point right after the last iteration.
Paper organization: Section 2 describes the motivation of our eigencurve scheduler. Section 3 presents the exact form and convergence rate of the proposed eigencurve scheduler, along with the lower bound of step decay. Section 4 shows the experimental results. Section 5 discusses the discovery and limitation of eigencurve and Section 6 gives our conclusion.

heading: 2 MOTIVATION
text: In this section, we will give the main motivation and intuition of our eigencurve learning rate scheduler. We first give the scheduling strategy to achieve the optimal convergence rate in the case that the Hessian is diagonal. Then, we show that the inverse time learning rate is sub-optimal in
most cases. Comparing these two scheduling methods brings up the reason why we should design eigenvalue distribution dependent learning rate scheduler.
Letting H be a diagonal matrix diag(λ1, λ2, . . . , λd) and reformulating Eqn. (1.5), we have
wt+1 − w∗ =wt − w∗ − ηt(H(ξ)wt − b(ξ)) =wt − w∗ − ηt(Hwt − b) + ηt (Hwt − b− (H(ξ)wt − b(ξ))) =wt − w∗ − ηt(Hwt − b− (Hw∗ − b)) + ηt (Hwt − b− (H(ξ)wt − b(ξ))) = (I − ηtH) (wt − w∗) + ηtnt.
It follows, E [ λj (wt+1,j − w∗,j)2 ] E[nt]=0 = λj [ (1− ηtλj)2(wt,j − w∗,j)2 + η2tE ‖[nt]j‖ 2 ]
(1.7) ≤ (1− ηtλj)2 · λj(wt,j − w∗,j)2 + η2t λ2jσ2.
(2.1)
Since H is diagonal, we can set step size scheduling for each dimension separately. Letting us choose step size ηt coordinately with the step size ηt,j = 1λj(t+1) being optimal for the j-th coordinate, then we have the following proposition.
Proposition 1. Assume that H is diagonal matrix with eigenvalues λ1 ≥ λ2 ≥ . . . , λd ≥ 0 and Eqn. (1.7) holds. If we set step size ηt,j = 1λj(t+1) , it holds that
2E [f(wt+1)− f(w∗)] = E  d∑ j=1 λj (wt+1,j − w∗,j)2  ≤ ∑dj=1 λj(w1,j − w∗,j)2 (t+ 1)2 +
t
(t+ 1)2 ·dσ2.
(2.2)
The leading equality here is proved in Appendix G.1, with the followed inequality proved in Appendix E. From Eqn. (2.2), we can observe that choosing proper step sizes coordinately can achieve the optimal convergence rate (Ge et al., 2019; Jain et al., 2018). Instead, if the widely used inverse time scheduler ηt = 1/(L + µt) is chosen, we can show that only a sub-optimal convergence rate can be obtained, especially when λj’s vary from each other.
Proposition 2. If we set the inverse time step size ηt = 1(L+µt) , then we have
2E [f(wt+1)− f(w∗)] = E  d∑ j=1 λj (wt+1,j − w∗,j)2 
≤ ( L+ µ
L+ µt )2 d∑ j=1 λj(w1,j − w∗,j)2 + d∑ j=1 ( λ2j 2λj − µ · 1 L+ µt · σ2 + λ2jσ 2 (L+ µt)2 ) . (2.3)
Remark 1. Eqn. (2.2) shows that if one can choose step size coordinate-wise with step size ηt,j = 1 λj(t+1) , then SGD can achieve a convergence rate
E [f(wT+1)− f(w∗)] ≤ O ( d T · σ2 ) . (2.4)
which matches the lower bound (Ge et al., 2019; Jain et al., 2018). In contrast, replacing L = λ1 and µ = λd in Proposition 2, we can obtain that the convergence rate of SGD being
E [f(wT+1)− f(w∗)] ≤ O  1 T d∑ j=1 λj λd · σ2  . (2.5)
Since it holds that λj ≥ λd, the convergence rate in Eqn. (2.4) is better than the one in Eqn. (2.5), especially when the eigenvalues of the Hessian (H matrix) decay rapidly. In fact, the upper bound in Eqn. (2.5) is tight for the inverse time decay scheduling (Ge et al., 2019).
Main Intuition Though SGD with coordinately proper step sizes can obtain the optimal convergence rate, since H is commonly non-diagonal in practice, one is not able to choose proper step sizes coordinately. However, the diagonal case still provide an important intuition for designing eigenvalue dependent learning rate scheduling. In fact, letting H = UΛU> be the spectral decomposition of the Hessian and setting w′ = U>w, then the Hessian becomes a diagonal matrix from perspective of updating w′, with the variance of the stochastic gradient being unchanged since U is a unitary matrix.
Proposition 1 and 2 imply that a good learning rate scheduler should decrease the error of each coordinate. The inverse time decay scheduler is only optimal for the coordinate related to the smallest eigenvalue. That’s the reason why the inverse time scheduler is sub-optimal. Thus, we should reduce the learning rate gradually such that we can run a optimal learning rate associated to λj to sufficiently drop the error of j-th coordinate. Furthermore, given the total iteration T and the eigenvalue distribution of the Hessian, we should allocate the running time for each optimal learning rate associated to λj .

heading: 3 EIGENVALUE DEPENDENT STEP SCHEDULING
text: Just as discussed in Section 2, to obtain better convergence rate for SGD, we should consider Hessian’s eigenvalue distribution and schedule the learning rate based on the distribution. In this section, we propose a novel learning rate scheduler for this task, which can be regarded as piecewise inverse time decay (see Figure 1). The method is very simple, we group eigenvalues according to their value and denote si to be the number of eigenvalues lie in the rangeRi = [µ · 2i, µ · 2i+1), that is,
si = #λj ∈ [µ · 2i, µ · 2i+1). (3.1)
Then, there are at most Imax = log2 κ such ranges. By the inverse time decay theory, the optimal learning rate associated to eigenvalues in the rangeRi should be
ηt = O (
1
2i−1µ · (t− ti)
) , with 0 = t0 < t1 < t2 < · · · < tImax = T. (3.2)
Our scheduling strategy is to run the optimal learning rate for eigenvalues in eachRi for a period to sufficiently decrease the error associated to eigenvalues inRi. To make the step size sequence {ηt}Tt=1 monotonely decreasing, we define the step sizes as
ηt = 1 L+ µ ∑i−1 j=1 ∆j2 j−1 + 2i−1µ(t− ti−1) if t ∈ [ti−1, ti) (3.3)
where
0 = t0 < t1 < t2 < · · · < tImax = T, ∆i = ti − ti−1, and Imax = log2 κ. (3.4)
To make the total error, that is the sum of error associated with Ri, to be small, we should allocate ∆i according to si−1’s. Intuitively, a large portion of eigenvalues lying in the range Ri should allocate a large portion of iterations. Specifically, we propose to allocate ∆i as follows:
∆i = √ si−1∑Imax−1
i′=0
√ si′ · T, with si = #λj ∈ [µ · 2i, µ · 2i+1). (3.5)
In the rest of this section, we will show that the step size scheduling according to Eqn. (3.3) and (3.5) can achieve better convergence rate than the one in Ge et al. (2019) when si is non-uniformly distributed. In fact, a better ∆i allocation can be calculated using numerical optimization.

heading: 3.1 THEORETICAL ANALYSIS
text: Lemma 1. Let objective function f(x) be quadratic and Assumption (1.7) hold. Running SGD for T -steps starting fromw0 and a learning rate sequence {ηt}Tt=1 defined in Eqn. (3.3), the final iterate
wT+1 satisfies
E [f(wT+1)− f(w∗)] ≤(f(w0)− f(w∗)) · κ2
∆21 +
15
2 · σ2µ Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 .
(3.6)
Since the bias term is a high order term, the variance term in Eqn. (3.6) dominates the error for wT+1. For simplicity, instead of using numerical methods to find the optimal {∆i}, we propose to use ∆i defined in Eqn. (3.5). The value of ∆i is linear to square root of the number of eigenvalues lying in the range [µ · 2i−1, µ · 2i). Using such ∆i, eigencurve has the following convergence property. Theorem 1. Let objective function f(x) be quadratic and Assumption (1.7) hold. Running SGD for T -steps starting from w0, a learning rate sequence {ηt}Tt=1 defined in Eqn. (3.3) and ∆i defined in Eqn. (3.5), the final iterate wT+1 satisfies
E [f(wT+1)− f(w∗)] ≤(f(w0)− f(w∗)) · κ2 ·
(∑Imax−1 i=0 √ si )2 s0T 2 + 15 (∑Imax−1 i=0 √ si )2 T · σ2.
Please refer to Appendix D, F and G for the full proof of Lemma 1 and Theorem 1. The variance term 15( ∑Imax−1 i=0 √ si) 2
T · σ 2 in above theorem shows that when si’s vary largely from each other, then the variance can be close to O ( d T · σ 2 )
which matches the lower bound (Ge et al., 2019). For example, letting Imax = 100, s0 = 0.99d and si = 0.0199 d, we can obtain that(∑99
i=0
√ si )2 T · σ2 = ( √ 0.99 + 99× √ 0.01/99)2 · d T · σ2 < 4d T · σ2.
We can observe that if the variance of si’s is large, the variance term in Theorem 1 can be close to dσ2/T . More generally, as rigorously stated in Corollary 2, eigencurve achieves minimax optimal convergence rate if the Hessian spectrum satisfies an extra assumption of what-we-called “power power law”: the density of eigenvalue λ is exponentially decaying with increasing value of λ in log scale, i.e. ln(λ). This assumption comes from the observation of estimated Hessian spectrums in practice (see Figure 2), which will be further illustrated in Section 4.1. Corollary 2. Given the same setting as in Theorem 1, when Hessian H’s eigenvalue distribution p(λ) satisfies “power power law”, i.e.
p(λ) = 1 Z · exp(−α(ln(λ)− ln(µ))) = 1 Z · (µ λ )α (3.7)
for some α > 1, where Z = ∫ L µ
(µ/λ)αdλ, there exists a constant C(α) which only depends on α, such that the final iterate wT+1 satisfies
E [f(wT+1)− f(w∗)] ≤ ( (f(w0)− f(w∗)) · κ2
T 2 + dσ2 T
) · C(α).
Please refer to Appendix G.3 for the proof. As for the worst-case guarantee, it is easy to check that only when si’s are equal to each other, that is, si = d/Imax = d/ log2(κ), the variance term reaches its maximum. Corollary 3. Given the same setting as in Theorem 1, when si = d/ log2(κ) for all 0 ≤ i ≤ Imax − 1, the variance term in Theorem 1 reaches its maximum and wT+1 satisfies
E [f(wT+1)− f(w∗)] ≤ (f(w0)− f(w∗)) · κ2 log2 κ
T 2 + 15d · log κ T σ2.
Remark 2. When si’s vary from each other, our eigenvalue dependent learning rate scheduler can achieve faster convergence rate than eigenvalue independent scheduler such as step decay which suffers from an extra log(T ) term (Ge et al., 2019). Only when si’s are equal to each other, Corollary 3 shows that the bound of variance matches to lower bound up to log κ which is same to the one in Proposition 3 of Ge et al. (2019).
Furthermore, we show that this log T gap between step decay and eigencurve certainly exists for problem instances of skewed Hessian spectrums. For simplicity, we only discuss the case where H is diagonal. Theorem 4. Let objective function f(x) be quadratic. We run SGD for T -steps starting fromw0 and a step decay learning rate sequence {ηt}Tt=1 defined in Algorithm 1 of Ge et al. (2019) with η1 = 1/L. As long as (1) H is diagonal, (2) The equality in Assumption (1.7) holds, i.e. Eξ [ ntn > t ] = σ2H and (3) T/ log T ≥ max(κ, 216, 16), the final iterate wT+1 satisfies,
E [f(wT+1)− f(w∗)] = Ω ( dσ2
T · log T ) The proof is provided in Appendix G.4. Removing this extra log T term may not seem to be a big deal in theory, but experimental results suggest the opposite.

heading: 4 EXPERIMENTS
text: To demonstrate eigencurve ’s practical value, empirical experiments are conducted on the task of image classification 1. Two well-known dataset are used: CIFAR-10 (Krizhevsky et al., 2009) and ImageNet (Deng et al., 2009). For full experimental results on more datasets, please refer to Appendix A.

heading: 4.1 HESSIAN SPECTRUM’S SKEWNESS IN PRACTICE
text: According to estimated2 eigenvalue distributions of Hessian on CIFAR-10 and ImageNet, as shown in Figure 3, it can be observed that all of them are highly skewed and share a similar tendency: A large portion of small eigenvalues and a tiny portion of large eigenvalues. This phenomenon has also been observed and explained by other researchers in the past(Sagun et al., 2017; Arjevani & Field, 2020). On top of that, when we plot both eigenvalues and density in log scale, the “power power law” arises. Therefore, if the loss surface can be approximated by quadratic objectives, then eigencurve has already achieved optimal convergence rate for those practical settings. The exact values of the extra constant terms are presented in Appendix A.2.
4.2 IMAGE CLASSIFICATION ON CIFAR-10 WITH EIGENCURVE SCHEDULING
This optimality in theory induces eigencurve ’s superior performance in practice, which is demonstrated in Table 2 and Figure 4. The full set of figures are available in Appendix A.7. All models are trained with stochastic gradient descent (SGD), no momentum, batch size 128 and weight decay wd = 0.0005. For full details of the experiment setup, please refer to Appendix B.

heading: 4.3 INSPIRED PRACTICAL SCHEDULES WITH SIMPLE FORMS
text: By simplifying the form of eigencurve and capturing some of its key properties, two simple and practical schedules are proposed: Elastic Step Decay and Cosine-power Decay, whose empirical performance are better than or at least comparable to cosine decay. Due to page limit, we leave all the experimental results in Appendix A.4, A.5, A.6.
Elastic Step Decay: ηt = η0/2k , if t ∈ [ (1− rk)T, (1− rk+1)T ) (4.1)
Cosine-power Decay: ηt = ηmin + (η0 − ηmin) [ 1
2 (1 + cos(
t
tmax π))
]α (4.2)
1Code: https://github.com/opensource12345678/why_cosine_works/tree/main 2Please refer to Appendix B.2 for details of the estimation and preprocessing procedure.

heading: 5 DISCUSSION
text: Cosine Decay and Eigencurve For ResNet-18 on CIFAR-10 dataset, eigencurve scheduler presents an extremely similar learning rate curve to cosine decay, especially when the number of training epochs is set to 100, as shown in Figure 5. This directly links cosine decay to our theory: the empirically superior performance of cosine decay is very likely to stem from the utilization of the “skewness” among Hessian matrix’s eigenvalues. For other situations, especially when the number of iterations is small, as shown in Table 2, eigencurve presents a better performance than cosine decay .
Sensitiveness to Hessian’s Eigenvalue Distributions One limitation of eigencurve is that it requires a precomputed eigenvalue distribution of objective functions’s Hessian matrix, which can be time-consuming for large models. This issue can be overcome by reusing the estimated eigenvalue distribution from similar settings. Further experiments on CIFAR-10 suggest the effectiveness of this approach. Please refer to Appendix A.3 for more details. This evidence suggests that eigencurve’s performance is not very sensitive to estimated eigenvalue distributions.
Relationship with Numerically Near-optimal Schedulers In Zhang et al. (2019), a dynamic programming algorithm was proposed to find almost optimal schedulers if the exact loss of the quadratic objective is accessible. While it is certainly the case, eigencurve still possesses several additional advantages over this type of approaches. First, eigencurve can be used to find simpleformed schedulers. Compared with schedulers numerically computed by dynamic programming, eigencurve provides an analytic framework, so it is able to bypass the Hessian spectrum estimation process if some useful assumptions of the Hessian spectrum can be obtained, such as ”power power law”. Second, eigencurve has a clear theoretical convergence guarantee. Dynamic programming can find almost optimal schedulers, but the convergence property of the computed scheduler is still unclear. Our work fills this gap.

heading: 6 CONCLUSION
text: In this paper, a novel learning rate schedule named eigencurve is proposed, which utilizes the “skewness” of objective’s Hessian matrix’s eigenvalue distribution and reaches minimax optimal convergence rates for SGD on quadratic objectives with skewed Hessian spectrums. This condition of skewed Hessian spectrums is observed and indeed satisfied in practical settings of image classification. Theoretically, eigencurve achieves no worse convergence guarantee than step decay for quadratic functions and reaches minimax optimal convergence rate (up to a constant) with skewed Hessian spectrums, e.g. under “power power law”. Empirically, experimental results on CIFAR-10 show that eigencurve significantly outperforms step decay, especially when the number of epochs is small. The idea of eigencurve offers a possible explanation for cosine decay’s effectiveness in practice and inspires two practical families of schedules with simple forms.

heading: A MORE EXPERIMENTAL RESULTS

heading: A.1 RIDGE REGRESSION
text: We compare different types of schedulings on ridge regression
f(w) = 1
n ||Xw − Y ||22 + α||w||22.
This experiment is only an empirical proof of our theory. In fact, the optima of ridge regression has a closed form and can be directly computed with
w∗ = ( X>X + nαI )−1 XTY.
Thus the optimal training loss f(w∗) can be calculated accordingly. In all experiments, we use the loss gap f(wT )− f(w∗) as our performance metric. Experiments are conducted on a4a datasets (Chang & Lin, 2011; Dua & Graff, 2017) (https:// www.csie.ntu.edu.tw/˜cjlin/libsvmtools/datasets/binary.html#a4a/), which contains 4, 781 samples and 123 features. This dataset is chosen majorly because it has a moderate number of samples and features, which enables us to compute the exact Hessian matrix H = 2(X>X/n + αI) and its corresponding eigenvalue distribution in acceptable time and space consumption.
In all of our experiments, we set α = 10−3. The model is optimized via SGD without momentum, with batch size 1, initial learning rate η0 ∈ {0.1, 0.06, 0.03, 0.02, 0.01, 0.006, 0.003, 0.002, 0.001, 0.0006, 0.0003, 0.0002, 0.0001} and learning rate of last iteration ηmin ∈ {0.1, 0.01, 0.001, 0.0001, 0.00001, 0, “UNRESTRICTED”}. Here “UNRESTRICTED” denotes the case where ηmin is not set, which is useful for eigencurve, who can decide the learning rate curve without setting ηmin. Given η0 and ηmin, we adjust all schedulers as follows. For inverse time decay ηt = η0/(1 + γη0t) and exponential decay ηt = γtη0, the hyperparameter γ is computed accordingly based on η0 and ηmin. For cosine decay, η0 and ηmin is directly used, with no restart adopted. For eigencurve, the learning rate curve is linearly scaled to match the given ηmin.
In addition, for eigencurve, we use the eigenvalue distribution of the Hessian matrix, which is directly computed via eigenvalue decomposition, as shown in Figure 6.
All experimental results demonstrate that eigencurve can obtain similar or better training losses when compared with other schedulers, as shown in Table 3.

heading: A.2 EXACT VALUE OF THE EXTRA TERM ON CIFAR-10 EXPERIMENTS
text: In Section 4.1, we have already given the qualitative evidence that shows eigencurve ’s optimality for practical settings on CIFAR-10. Here we strengthen this argument by providing the quantitative evidence as well. The exact value of the extra term is presented in Table 4, where we assume CIFAR-10 has batch size 128, number of epochs 200 and weight decay 5 × 10−4, while ImageNet has batch size 256, number of epochs 90 and weight decay 10−4.
It is worth noticing that the extra term’s value of eigencurve is independent from the number of iterations T , since the value ( ∑Imax−1 i=0 √ si)
2/d only depends on the Hessian spectrum. So basically eigencurve has already achieved the minimax optimal rate (up to a constant) for models and datasets listed in Table 4, if the loss landscape around the optima can be approximated by quadratic functions. For full details of the estimation process, please refer to Appendix B.
A.3 REUSING EIGENCURVE FOR DIFFERENT MODELS ON CIFAR-10
For image classification tasks on CIFAR-10, we check the performance of reusing ResNet-18’s eigenvalue distribution for other models. As shown in Table 5, experimental results demonstrate that Hessian’s eigenvalue distribution of Resnet-18 on CIFAR-10 can be applied to GoogLeNet/VGG16 and still achieves good peformance. Here the experiment settings are exactly the same as Section 4.2 in main paper, except we only run one trial of experiment and use common step decay with three intervals + decay factor 10.
A.4 IMAGENET CLASSIFICATION WITH ELASTIC STEP DECAY
One key observation in CIFAR-10 experiments is the existence of “power power law” shown in Figure 3. Also, notice that in the form of eigencurve , specifically Eqn. (3.5), iteration interval length ∆i is proportional to the square root of eigenvalue density si in range [µ · 2i, µ · 2i+1). Combining those two facts together, it suggests the length of “learning rate interval” should have lengths exponentially decreasing.
Based on this idea, Elastic Step Decay (ESD) is proposed, which has the following form,
ηt = η0/2 k , if t ∈ [ (1− rk)T, (1− rk+1)T ) Compared to general step decay with adjustable interval lengths, elastic step decay does not require manual adjustment for the length of each interval. Instead, they are all controlled by one hyperparameter r ∈ (0, 1), which decides the “shrinking speed” of interval lengths. Experiments on CIFAR-10, CIFAR-100 and ImageNet demonstrate its superiority in practice, as shown in Table 6, Table 7.
For experiments on CIFAR-10/CIFAR-100, we adopt the same settings as eigencurve , except we only run one trial of experiment and use common step decay with three intervals + decay factor 10. For experiments on ImageNet, we use ResNet-50 trained via SGD without momentum, batch size 256 and weight decay wd = 10−4. Since no momentum is used, the initial learning rate is set to η0 = 1.0 instead of η0 = 0.1. The hyperparameters ηmin is set to be 0 for cosine decay scheduler. As for the dataset, we use the common ILSVRC 2012 dataset, which contains 1000 classes, around 1.2M images for training and 50,000 images for validation. For all experiments, we search r ∈ {1/2, 1/ √ 2} for ESD, with other hyperparameter search and selection process being the same as eigencurve .

heading: A.5 LANGUAGE MODELING WITH ELASTIC STEP DECAY
text: More experiments on language modeling are conducted to further demonstrate Elastic Step Decay’s superiority over other schedulers.
For all experiments, we follow almost the same setting in Zaremba et al. (2015), where a large regularized LSTM recurrent neural network (Hochreiter & Schmidhuber, 1997) is trained on Penn Treebank (Marcus et al., 1993) for language modeling task. The Penn Treebank dataset has a training set of 929k words, a validation set of 73k words and a test set of 82k words. SGD without momentum is adopted for training, with batch size 20 and 35 unrolling steps in LSTM.
Other details are exactly the same, except for the number of training epochs. In Zaremba et al. (2015), it uses 55 epochs to train the large regularized LSTM, which is changed to 30 epochs in our setting, since we found that the model starts overfitting after 30 epochs. We conducted hyperparameter search for all schedules, as shown in Table 8.
Experimental results show that Elastic Step Decay significantly outperforms other schedulers, as shown in Table 9.
A.6 IMAGE CLASSIFICATION ON IMAGENET WITH COSINE-POWER SCHEDULING
Another key observation in CIFAR-10 experiments is that eigencurve ’s learning rate curve shape changes in a fixed tendency: more “concave” learning rate curves for less training epochs, which inspire the cosine-power schedule in following form.
Cosine-power : ηt = ηmin + (η0 − ηmin) [ 1
2 (1 + cos(
t
tmax π))
]α
Results in Table 10 show the schedulings’ performance with α = 0.5/1/2, which are denoted as√ Cosine/Cosine/Cosine2 respectively. Notice that the best scheduler gradually moves from small α to larger α when the number of epochs increases. For #epoch=270, since the number of epochs is large enough to make model converge, it is reasonable that the accuracy gap between all schedulers is small.
For experiments on ImageNet, we use ResNet-18 trained via SGD without momentum, batch size 256 and weight decay wd = 10−4. Since no momentum is used, the initial learning rate is set to η0 = 1.0 instead of η0 = 0.1. The hyperparameters ηmin is set to be 0 for all cosine-power scheduler. As for the dataset, we use the common ILSVRC 2012 dataset, which contains 1000 classes, around 1.2M images for training and 50,000 images for validation.

heading: A.7 FULL FIGURES FOR EIGENCURVE EXPERIMENTS IN SECTION 4.2
text: Please refer to Figure 8, 9, 10, 11, 12 and 13.

heading: B DETAILED EXPERIMENTAL SETTINGS FOR IMAGE CLASSIFICATION ON
text: CIFAR-10/CIFAR-100

heading: B.1 BASIC SETTINGS
text: As mentioned in the main paper, all models are trained with stochastic gradient descent (SGD), no momentum, batch size 128 and weight decay wd = 0.0005. Furthermore, we perform a grid search to choose the best hyperparameters of all schedulers, with a validation set created from 5, 000
samples in the training set, i.e. one-tenth of the training set. The remaining 45, 000 samples are then used for training the model. After obtaining hyperparameters with the best validation accuracy, we train the model again with the full training set and test the trained model on test set, where 5 trials of experiments are conducted with different random seeds. The mean and standard deviation of the test results are reported.
Here the grid search explores hyperparameters η0 ∈ {1.0, 0.6, 0.3, 0.2, 0.1} and ηmin ∈ {0.01, 0.001, 0.0001, 0, “UNRESTRICTED”}, where η0 denotes the initial learning rate and ηmin stands for the learning rate of last iteration. “UNRESTRICTED” denotes the case where ηmin is not set, which is useful for eigencurve, who can decide the learning rate curve without setting ηmin. Given η0 and ηmin, we adjust all schedulers as follows. For inverse time decay, the hyperparameter γ is computed accordingly based on η0 and ηmin. For cosine decay, η0 and ηmin is directly used, with no restart adopted. For general step decay, we search the interval number in {3, 4, 5, blog T c, blog T c + 1} and decay factor in {2, 5, 10}. For step decay proposed in Ge et al. (2019), the interval number is fixed to be blog T c, along with a decay factor 2. For eigencurve, two major modifications are made to make it more suitable for practical settings:
ηt = 1/L 1 + 1κ ∑i−1 j=1 ∆j2 j−1 + 2 i−1 κ (t− ti−1) =
η0 1 + 1κ ∑i−1 j=1 ∆jβ j−1 + β i−1 κ (t− ti−1) .
Here we change 1/L to η0 so that it is possible to adjust the initial learning rate of eigencurve. We also change the fixed constant 2 to a general constant β > 1, which is aimed at making the learning rate curve smoother. The learning rate curve of eigencurve is then linearly scaled to match the given ηmin.
Notice that the learning rate η0 can be larger than 1/L, while the loss still does not explode. There are several explanations for this phenomenon. First, in basic non-smooth analysis of GD and SGD with inverse time decay scheduler, the learning rate can be larger than 1/L if the gradient norm is bounded (Shamir & Zhang, 2012). Second, deep learning has a non-convex loss landscape, especially when the parameter is far away from the optima. Hence it is common to use larger learning rate at first. As long as the loss does not explode, it is okay. So we still include large learning rate η0 in our grid search process.
B.2 SETTINGS FOR EIGENCURVE
In addition, for our eigencurve scheduler, we use PyHessian (Yao et al., 2020) to generate Hessian matrix’s eigenvalue distribution for all models. The whole process consists of three phases, which are illustrated as follows.
1) Training the model Almost all CNN models on CIFAR-10 have non-convex objectives, thus the Hessian’s eigenvalue distributions are different for different parameters. Normally, we want the this distribution to reflect the overall tendency of most parts of the training process. According to the phenomenon demonstrated in Appendix E, figure A.11-A.17 of Yao et al. (2020), the eigenvalue distribution of ResNet’s Hessian presents similar tendency after training 30 epochs, which suggests that the Hessian’s eigenvalue distribution can be used after sufficient training.
In all CIFAR-10 experiments, we use the Hessian’s eigenvalue distribution of models after training 180 epochs. Since the goal here is to sufficiently train the model, not to obtain good performance, common baseline settings are adopted for training. For all models used for eigenvalue distribution estimation, we adopt SGD with momentum = 0.9, batch size 128, weight decay wd = 0.0005 and initial learning rate 0.1. On top of that, we use step decay, which decays the learning rate by a factor of 10 at epoch 80 and 120. All of them are default settings of the PyHessian code (https://github.com/amirgholami/PyHessian/blob/master/training. py, commit: f4c3f77).
ImageNet adopts a similar setting, with training epochs being 90, SGD with momentum = 0.9, batch size 256, weight decay wd = 0.0001, inital learning rate 0.1 and step decay schedule decays learning rate by a factor of 10 at epoch 30 and 60.
2) Estimating Hessian matrix’s eigenvalue distribution for the trained model After obtaining the checkpoint of a sufficiently trained model, we then run PyHessian to estimate the Hessian’s
eigenvalue distribution for that checkpoint. The goal here is to obtain the Hessian’s eigenvalue distribution with sufficient precision. To be more specific, the length of intervals around each estimated eigenvalue. PyHessian estimates the eigenvalue spectral density (ESD) of a model’s Hessian, in other words, the output is a list of eigenvalue intervals, along with the density of each interval, where the whole density adds up to 1. Precision means the interval length here.
It is natural that the estimation precision is related to the complexity of the PyHessian algorithm, e.g. the better precision it yields, the more time and space it consumes. More specifically, the algorithm has a time complexity of O(Nn2vd) and space complexity O(Bd+ nvd), where d is the number of model parameters, N is the number of samples used for estimating the ESD, B is the batch size and nv is the iteration number of Stochastic Lanczos Quadrature used in PyHessian, which controls the estimation precision (see Algorithm 1 of Yao et al. (2020)).
In our experiments, we use nv = 5000 for ResNet-18 and nv = 3000 for GoogLeNet/VGG16, which gives an eigenvalue distribution estimation with precision around 10−5 to 10−4. N and B are both set to 200 due to GPU memory constraint, i.e. we use one mini-batch to estimate the eigenvalue distribution. It turns out that this one-batch estimation is good enough and yields similar results to full dataset settings shown in Yao et al. (2020).
However, space complexity is still a bottleneck here. Due to the large number of nv and space complexity O(Bd + nvd) of PyHessian, the value of d cannot be very large. In practice, with a NVIDIA GeForce 2080 Ti GPU, which has around 11GB memory, the maximum acceptable parameter number d is around 200K − 400K. This implies that the model has to be compressed. In our experiments, we reduce the number of channels by a factor of C for all models. For ResNet-18, C = 16. For GoogLeNet, C = 4. For VGG16, C = 8. Notice that those compressed models are only used for eigenvalue distribution estimation. In experiments of comparing different scheduling, we still use the original model with no compression.
One may refer to https://github.com/opensource12345678/why_cosine_ works/tree/main/eigenvalue_distribution for generated eigenvalue distributions.
3) Generating eigencurve scheduler with the estimated eigenvalue distribution After obtaining the eigenvalue distribution, we do a preprocessing before plug it into our eigencurve scheduler.
First, we notice that there are negative eigenvalues in the final distribution. Theoretically, if the parameter is right at the optimal point, no negative eigenvalues should exist for Hessian matrix. Thus we conjecture that those negative eigenvalues are caused by the fact that the model is closed to optima w∗, but not exactly at that point. Furthermore, the estimation precision loss can be another cause. In fact, most of those negative eigenvalues are small, e.g. 98.6% of those negative eigenvalues lie in [−0.1, 0), and can be generally ignored without much loss. In our case, we set them to their absolute values.
Second, for a given weight decay value wd, we need to take the implicit L2 regularization into account, since it affects the Hessian matrix as well. Therefore, for all eigenvalues after the first step, we add wd to them.
After preprocessing, we plug the eigenvalue distribution into our eigencurve scheduler and generates the exact form of eigencurve.
ηt = 1/L 1 + 1κ ∑i−1 j=1 ∆j2 j−1 + 2 i−1 κ (t− ti−1) =
η0 1 + 1κ ∑i−1 j=1 ∆jβ j−1 + β i−1 κ (t− ti−1)
For experiments with 100 epochs, we set β = 1.000005, so that the learning rate curve is much smoother. For experiments with 10 epochs, we set β = 2.0. In our experiments, β serves as a fixed constant, not hyperparameters. So no hyperparameter search is conducted on β. One can do that in practice though, if computation resource allows.

heading: B.3 COMPUTE RESOURCE AND IMPLEMENTATION DETAILS
text: All the code for results in main paper can be found in https://github.com/ opensource12345678/why_cosine_works/tree/main, which is released under the MIT license.
All experiments on CIFAR-10/CIFAR-100 are conducted on a single NVIDIA GeForce 2080 Ti GPU, where ResNet-18/GoogLeNet/VGG16 takes around 20mins/90mins/40mins to train 100 epochs, respectively. High-precision eigenvalue distribution estimation, e.g. nv ≥ 3000, requires around 1-2 days to complete, but this is no longer necessary given the released results.
The ResNet-18 model is implemented in Tensorflow 2.0. We use tensorflow-gpu 2.3.1 in our code. The GoogLeNet and VGG16 model is implemented in Pytorch, specifically, 1.7.0+cu101.

heading: B.4 LICENSE OF PYHESSIAN
text: According to https://github.com/amirgholami/PyHessian/blob/master/ LICENSE, PyHessian (Yao et al., 2020) is released under the MIT License.

heading: C DETAILED EXPERIMENTAL SETTINGS FOR IMAGE CLASSIFICATION ON IMAGENET
text: One may refer to https://www.image-net.org/download for specific terms of access for ImageNet. The dataset can be downloaded from https://image-net.org/challenges/ LSVRC/2012/2012-downloads.php, with training set being “Training images (Task 1 & 2)” and validation set being “Validation images (all tasks)”. Notice that registration and verification of institute is required for successful download.
All experiments on ImageNet are conducted on two NVIDIA GeForce 2080 Ti GPUs with data parallelism. The ResNet-18 model takes around 2 days to train 90 epochs, about 20mins-30mins per epoch. The ResNet-18 model on ImageNet is implemented in Pytorch, specifically, 1.7.0+cu101.
D IMPORTANT PROPOSITIONS AND LEMMAS
Proposition 3. Letting f(x) be a monotonically increasing function in the range [t0, t̃], then it holds that
t̃−1∑ k=t0 f(k) ≤ ∫ t̃ t0 f(x) dx. (D.1)
If f(x) is monotonically decreasing in the range [t0, t̃], then it holds that∫ t̃ t0 f(x) dx ≤ t̃−1∑ k=t0 f(k) ≤ t̃∑ k=t0 f(k) ≤ f(t0) + ∫ t̃ t0 f(x) dx. (D.2)
Lemma 2. Function f(x) = exp(−αx)x2 with 0 < α and x ∈ (0, 1] is monotone decreasing in the range x ∈ ( 2α ,+∞) and monotone increasing in the range x ∈ [0, 2 α ].
Proof. We can obtain the derivative of f(x) as
∇f(x) = x exp(−αx)(2− αx).
Thus, it holds that∇f(x) ≥ 0 when x ∈ [0, 2α ]. This implies that f(x) is monotone increasing when x ∈ [0, 2α ]. Similarly, we can obtain that f(x) is monotone decreasing when x ∈ ( 2 α ,+∞). Lemma 3. It holds that∫ exp(−αx)x dx = −α−1(x exp(−αx) + α−1 exp(−αx)). (D.3)
Proof.∫ exp(−αx)x dx =− α−1 ∫ x d exp(−αx) = −α−1(x exp(−αx)− ∫ exp(−αx)dx)
=− α−1(x exp(−αx) + α−1 exp(−αx)).

heading: E PROOF OF SECTION 2
text: Proof of Proposition 1. By iteratively applying Eqn. (2.1), we can obtain that
E [ λj (wt+1,j − w∗,j)2 ] (2.1) ≤ t∏ i=1 (1− ηi,jλj)2 · λj(w1,j − w∗,j)2
+ λ2jσ 2 · t∑ k=1 t∏ i=k+1 (1− ηi,jλj)2η2k,j
= λj(w1,j − w∗,j)2
(t+ 1)2 + σ2 · t∑ k=1 (k + 1)2 (t+ 1)2 · 1 (k + 1)2
= λj(w1,j − w∗,j)2
(t+ 1)2 +
t
(t+ 1)2 · σ2.
Summing up each coordinate, we can obtain the result.
Proof of Proposition 2. Let us denote σ2j = λ 2 jσ 2. By Eqn. (2.1), we can obtain that E [ λj (wt+1,j − w∗,j)2 ] ≤Πti=1(1− ηi,jλj)2 · λj(w1,j − w∗,j)2 +
t∑ k=1 Πi=k(1− ηi,jλj)η2k,jσ2j
≤ exp ( −2
t∑ i=1 ηi,jλj
) · λj(w1,j − w∗,j)2 +
t∑ k=1 exp
( −2
t∑ i=k ηi,jλj ) η2k,jσ 2 j
= exp ( −2λj
t∑ i=1 1 L+ µi
) · λj(w1,j − w∗,j)2 +
t∑ k=1 exp
( −2λj
t∑ i=k 1 L+ µi ) η2k,jσ 2 j
≤ exp (
2λj µ ln
( L+ µ
L+ µt
)) · λj(w1,j − w∗,j)2 + t∑ k=1 exp ( 2λj µ ln ( L+ µk L+ µt )) · σ2j (L+ µk)2
=
( L+ µ
L+ µt
) 2λj µ · λj(w1,j − w∗,j)2 + t∑
k=1
(L+ µk)
( 2λj µ −2 ) (L+ µt) 2λj µ · σ2j
≤ ( L+ µ
L+ µt
) 2λj µ
· λj(w1,j − w∗,j)2 + σ2j 2λj − µ (L+ µt)
2λj µ −1
(L+ µt) 2λj µ
+ σ2j
(L+ µt)2
=
( L+ µ
L+ µt
) 2λj µ
· λj(w1,j − w∗,j)2 + 1 2λj − µ · 1 L+ µt · σ2j + σ2j (L+ µt)2 .
The third inequality is because function F (x) = 1/(L + µx) is monotone decreasing in the range [1,∞), and it holds that
t∑ i=k 1 L+ µi ≥ ∫ t i=k
1
L+ µi di =
1 µ ln
( L+ µt
L+ µk
) .
The last inequality is because function F (x) = (L + µx)2λj/µ−2 is monotone increasing in the range [0,∞), and it holds that
t∑ k=1 (L+ µk) 2λj µ −2 ≤ ∫ t k=1 (L+ µk) 2λj µ −2dk + (L+ µt) 2λj µ −2
= 1 µ (
2λj µ − 1 ) · ((L+ µt) 2λjµ −1 − (L+ µ) 2λjµ −1)+ (L+ µt) 2λjµ −2 < 1
2λj − µ (L+ µt)
2λj µ −1 + (L+ µt) 2λj µ −2.
By µ ≤ λj , σ2j = λ2jσ2, and summing up from i = 1 to d, we can obtain the result.

heading: F PRELIMINARIES
text: Lemma 4. Let objective function f(x) be quadratic. Running SGD for T -steps starting from w0 and a learning rate sequence {ηt}Tt=1, the final iterate wT+1 satisfies
E [ (wT+1 − w∗)>H(wT+1 − w∗) ] =E [ (w0 − w∗)> · PT . . . P0HP0 . . . PT · (w0 − w∗)
] +
T∑ τ=0 E [ η2τn > τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ] ,
(F.1)
where Pt = I − ηtH .
Proof. Reformulating Eqn. (1.5), we have
wt+1 − w∗ =wt − w∗ − ηt(H(ξ)wt − b(ξ)) =wt − w∗ − ηt(Hwt − b) + ηt (Hwt − b− (H(ξ)wt − b(ξ))) =wt − w∗ − ηt(Hwt − b− (Hw∗ − b)) + ηt (Hwt − b− (H(ξ)wt − b(ξ))) = (I − ηtH) (wt − w∗) + ηtnt =Pt(wt − w∗) + ηtnt.
Thus, we can obtain that
wt+1 − w∗ = Pt . . . P0(w0 − w∗) + t∑
τ=0
Pt . . . Pτ+1ητnτ . (F.2)
We can decompose above stochastic process associated with SGD’s update into two simpler processes as follows:
wbt+1 − w∗ = Pt(wbt − w∗), and wvt+1 − w∗ = Pt(wvt − w∗) + ηtnt, with wv0 = w∗, (F.3)
which entails that
wbt+1 − w∗ = Pt . . . P0(wb0 − w∗) = P0 . . . Pt(wb0 − w∗) (F.4)
wvt+1 − w∗ = t∑
τ=0
Pt . . . Pτ+1ητnτ = t∑ τ=0 Pτ+1 . . . Ptητnτ (F.5)
(F.2)⇒wt+1 − w∗ = ( wbt+1 − w∗ ) + ( wvt+1 − w∗ ) (F.6)
where the last equality in Eqn. (F.4) and Eqn. (F.5) is because the commutative property PtPt′ = (I − ηtH)(I − ηt′H) = (I − ηt′H)(I − ηtH) = Pt′Pt holds for ∀t, t′.
Thus, we have
E [ (wT+1 − w∗)>H(wT+1 − w∗) ] (F.6) = E [ (wbT+1 − w∗)>H(wbT+1 − w∗) + 2(wvT+1 − w∗)>H(wb0 − w∗)
+(wvT+1 − w∗)>H(wvT+1 − w∗) ] .
(F.4) = E [ (wb0 − w∗)>PT . . . P0HP0 . . . PT (wb0 − w∗) + 2(wvT+1 − w∗)>HP0 . . . PT (wb0 − w∗)
+(wvT+1 − w∗)>H(wvT+1 − w∗) ]
(F.5) = E [ (wb0 − w∗)>PT . . . P0HP0 . . . PT (wb0 − w∗)
+ 2 ( T∑ τ=0 PT . . . Pτ+1ητnτ )> HP0 . . . PT (w b 0 − w∗)
+ ( T∑ τ=0 PT . . . Pτ+1ητnτ )> H ( T∑ τ=0 PT . . . Pτ+1ητnτ ) E[nτ ]=0
= E [ (wb0 − w∗)>PT . . . P0HP0 . . . PT (wb0 − w∗) ] + E
 T∑ τ=0,τ ′=0 ητητ ′ · n>τ PT . . . Pτ+1 ·H · Pτ ′+1 . . . PTnτ ′  =E [ (wb0 − w∗)>PT . . . P0HP0 . . . PT (wb0 − w∗)
] +
T∑ τ=0 E [ η2τn > τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ] ,
where the last equality is because when τ and τ ′ are different, it holds that
E[n>τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ′ ] = 0
due to independence between nτ and nτ ′ .
Lemma 5. Given the assumption that Eξ [ ntn > t ] σ2H , then the variance term satisfies that
T∑ τ=0 E [ η2τn > τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ] ≤ σ2 d∑ j=1 λ2j T∑ k=0 η2k T∏ i=k+1 (1− ηiλj)2 , (F.7)
where Pt = I − ηtH .
Proof. Denote Aτ , PT . . . Pτ+1H 1 2 , then
A>τ = ( PT . . . Pτ+1H 1 2 )> = ( H 1 2 )> P>τ+1 . . . P > T = H 1 2Pτ+1 . . . PT , (F.8)
where the second equality is entailed by the fact that H 1 2 , Pτ+1, . . . , PT are symmetric matrices.
Therefore, we have, T∑ τ=0 E [ η2τn > τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ] (F.8) =
T∑ τ=0 E [ η2τn > τ AτA > τ nτ ] = T∑ τ=0 η2τE [ tr ( n>τ AτA > τ nτ )] = T∑ τ=0 η2τE [ tr ( A>τ nτn > τ Aτ )] =
T∑ τ=0 η2τ tr ( E [ A>τ nτn > τ Aτ ]) = T∑ τ=0 η2τ tr ( A>τ E [ nτn > τ ] Aτ )
≤σ2 · T∑ τ=0 η2τ tr ( A>τ HAτ ) = σ2 · T∑ τ=0 η2τ tr ( AτA > τ H ) =σ2 ·
T∑ τ=0 η2τ tr (PT . . . Pτ+1HPτ+1 . . . PTH)
=σ2 d∑ j=1 λ2j T∑ k=0 η2k T∏ i=k+1 (1− ηiλj)2 ,
where the third and sixth equality come from the cyclic property of trace, while the first inequality is because of the condition Eξ [ ntn > t ] σ2H , where
∀x, x>E[nτn>τ ]x ≤ σ2x>Hx ⇒ ∀z, z>A>τ E[nτn>τ ]Aτz = (Aτz)>E[nτn>τ ](Aτz) ≤ σ2(Aτz)>H(Aτz) = σ2z>A>τ HAτz ⇒ A>τ E[nτn>τ ]Aτ σ2A>τ HAτ ⇒ tr ( A>τ E[nτn>τ ]Aτ ) ≤ σ2tr ( A>τ HAτ ) .
Lemma 6. Letting λj̃ be the smallest positive eigenvalue of H , then the bias term satisfies that E [ (w0 − w∗)> · PT . . . P0HP0 . . . PT · (w0 − w∗) ] ≤(w0 − w∗)>H(w0 − w∗) · exp ( −2λj̃
T∑ k=0 ηk
) . (F.9)
Proof. Letting H = UΛU> be the spectral decomposition of H and uj be j-th column of U , we can obtain that
E [ (w0 − w∗)> · PT . . . P0HP0 . . . PT · (w0 − w∗) ] =
d∑ j=1 λj · (u>j (w0 − w∗))2 · T∏ k=0 (1− ηkλj)2
≤ d∑ j=1 λj · (u>j (w0 − w∗))2 · exp
( −2λj
T∑ k=0 ηk
) .
Since λj̃ is the smallest positive eigenvalue of H , it holds that d∑ j=1 λj · (u>j (w0 − w∗))2 · exp ( −2λj T∑ k=0 ηk ) ≤ d∑ j=1 λj · (u>j (w0 − w∗))2 · exp ( −2λj̃ T∑ k=0 ηk )
=(w0 − w∗)>H(w0 − w∗) · exp ( −2λj̃ T∑ k=0 ηk ) .

heading: G PROOF OF THEOREMS
text: Lemma 7. Let learning rate ηt is defined in Eqn. (3.3). Assuming k ∈ [ti′−1, ti′ ] with 1 ≤ i′ ≤ ĩ ≤ T , the sequence {ηt}Tt=0 satisfies that
tĩ+1−1∑ t=k ηt ≥ ĩ+1∑ i=i′+1 1 2i−1µ ln αi αi−1 + 1 2i′−1µ ln
αi′
αi′−1 + 2i ′−1µ(k − ti′−1)
, (G.1)
where αi is defined as
αi , L+ µ i∑
j=1
∆j2 j−1 =
1
ηti . (G.2)
Proof. First, we divide learning rates into two groups: those who are guaranteed to cover a full interval and those who may not.
tĩ+1−1∑ t=k ηt = tĩ+1−1∑ t=ti′ ηt + ti′−1∑ t=k ηt = ĩ+1∑ i=i′+1 ti−1∑ t=ti−1 ηt + ti′−1∑ t=k ηt
Furthermore, because ηt is monotonically decreasing with respect to t, by Proposition 3, we have
tĩ+1−1∑ t=k ηt (D.2) ≥ ĩ+1∑ i=i′+1 ∫ ti ti−1 ηtdt+ ∫ ti′ k ηtdt
(3.3) = ĩ+1∑ i=i′+1 ∫ ti ti−1
1 L+ µ ∑i−1 j=1 ∆j2 j−1 + 2i−1µ(t− ti−1) dt
+ ∫ ti′ k
1 L+ µ ∑i′−1 j=1 ∆j2 j−1 + 2i′−1µ(t− ti′−1) dt
= ĩ+1∑ i=i′+1 1 2i−1µ ln L+ µ
∑i−1 j=1 ∆j2
j−1 + 2i−1µ(ti − ti−1) L+ µ ∑i−1 j=1 ∆j2 j−1
+ 1 2i′−1µ ln L+ µ
∑i′−1 j=1 ∆j2 j−1 + 2i ′−1µ(ti′ − ti′−1)
L+ µ ∑i′−1 j=1 ∆j2 j−1 + 2i′−1µ(k − ti′−1)
= ĩ+1∑ i=i′+1 1 2i−1µ ln L+ µ
∑i j=1 ∆j2 j−1
L+ µ ∑i−1 j=1 ∆j2 j−1
+ 1
2i′−1µ ln
L+ µ ∑i′ j=1 ∆j2 j−1
L+ µ ∑i′−1 j=1 ∆j2 j−1 + 2i′−1µ(k − ti′−1)
= ĩ+1∑ i=i′+1 1 2i−1µ ln αi αi−1 + 1 2i′−1µ ln
αi′
αi′−1 + 2i ′−1µ(k − ti′−1)
.
Lemma 8. Letting sequence {αi} be defined in Eqn. (G.2), given 1 ≤ ĩ, it holds that
ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i ti−1∑ k=ti−1 (αi−1 + 2 i−1µ(k − ti−1))2 ĩ−i+2−2
≤2 · 1 2ĩ+1µ · ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2(α2ĩ−i+2−1i − α2ĩ−i+2−1i−1 )α−2ĩ−i+2i . (G.3)
Proof. Notice that g(k) := (αi−1 + 2i−1µ(k − ti−1))2 ĩ−i+2−2 is a monotonically increasing function, we have,
ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i ti−1∑ k=ti−1 (αi−1 + 2 i−1µ(k − ti−1))2 ĩ−i+2−2
(D.1) ≤ ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i ∫ ti ti−1 (αi−1 + 2 i−1µ(t− ti−1))2 ĩ−i+2−2dt
= ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i ((2ĩ−i+2 − 1) · 2i−1µ)−1 (α2ĩ−i+2−1i − α2ĩ−i+2−1i−1 )
= ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i 1 (2ĩ+1 − 2i−1)µ ( α2 ĩ−i+2−1 i − α 2ĩ−i+2−1 i−1 ) ≤ ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i 1 (2ĩ+1 − 2ĩ)µ ( α2 ĩ−i+2−1 i − α 2ĩ−i+2−1 i−1 )
= ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i 1 2ĩ+1µ 1 1− 12 ( α2 ĩ−i+2−1 i − α 2ĩ−i+2−1 i−1 )
=2 · 1 2ĩ+1µ ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i (α2ĩ−i+2−1i − α2ĩ−i+2−1i−1 )
=2 · 1 2ĩ+1µ ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2(α2ĩ−i+2−1i − α2ĩ−i+2−1i−1 )α−2ĩ−i+2i .
Lemma 9. Letting {αi} be a positive sequence, given 1 ≤ ĩ, it holds that
ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2(α2ĩ−i+2−1i − α2ĩ−i+2−1i−1 )α−2ĩ−i+2i ≤ α−1ĩ+1. (G.4) Proof. First, we have
ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2(α2ĩ−i+2−1i − α2ĩ−i+2−1i−1 )α−2ĩ−i+2i =
ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−1i − α2ĩ−i+2−1i−1 α2 ĩ−i+2 i  =
ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−1i − ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α2ĩ−i+2−1i−1 α2 ĩ−i+2 i  =α−1
ĩ+1 + ĩ∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−1i − ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α2ĩ−i+2−1i−1 α2 ĩ−i+2 i  .
Furthermore, we reformulate the term ∑ĩ i=1 [∏ĩ+1 j=i+1 ( αj−1 αj )2ĩ−j+2] α−1i as follows
ĩ∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−1i = ĩ∑ i=1  ĩ+1∏ j=i+2 ( αj−1 αj )2ĩ−j+2( αi αi+1 )2ĩ−i+1 α−1i (G.5)
= ĩ∑ i=1  ĩ+1∏ j=i+2 ( αj−1 αj )2ĩ−j+2(α2ĩ−i+1−1i α2 ĩ−i+1 i+1 ) (G.6)
i′′=i+1 = ĩ+1∑ i′′=2  ĩ+1∏ j=i′′+1 ( αj−1 αj )2ĩ−j+2α2ĩ−i′′+2−1i′′−1 α2 ĩ−i′′+2 i′′  (G.7) i=i′′ =
ĩ+1∑ i=2  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α2ĩ−i+2−1i−1 α2 ĩ−i+2 i  . (G.8) Combining above results, we can obtain that
ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2(α2ĩ−i+2−1i − α2ĩ−i+2−1i−1 )α−2ĩ−i+2i =α−1
ĩ+1 + ĩ+1∑ i=2  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α2ĩ−i+2−1i−1 α2 ĩ−i+2 i  −
ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α2ĩ−i+2−1i−1 α2 ĩ−i+2 i  =α−1
ĩ+1 −  ĩ+1∏ j=2 ( αj−1 αj )2ĩ−j+2(α2ĩ+1−10 α2 ĩ+1 1 ) ≤α−1
ĩ+1 .
Lemma 10. Letting us denote vt+1,j , ∑t k=0 η 2 k ∏t i=k+1 (1− ηiλj)
2 with ηi defined in Eqn. (3.3), for 1 ≤ t ≤ t′, it holds that
vt′,j ≤ max(vt,j , ηt/λj). (G.9)
Proof. If vt+1,j ≤ max(vt,j , ηt/λj) holds for ∀t ≥ 1, then it naturally follows that
vt′,j ≤max ( vt′−1,j ,
ηt′−1 λj
) ≤ max ( vt′−2,j ,
ηt′−2 λj , ηt′−1 λj ) ) ≤ . . .
≤max ( vt,j ,
ηt λj , . . . ηt′−2 λj , ηt′−1 λj ) = max ( vt,j ,
ηt λj ) where the last equality is entailed by the fact that t ≤ t′ and ηt defined in Eqn. (3.3) is monotonically decreasing. We then prove vt+1,j ≤ max(vt,j , ηt/λj) holds for ∀t ≥ 1.
For ∀t ≥ 1, we have
vt+1,j = t∑ k=0 η2k t∏ i=k+1 (1− ηiλj)2
=η2t + t−1∑ k=0 η2k t∏ i=k+1 (1− ηiλj)2
=η2t + (1− ηtλj)2 t−1∑ k=0 η2k t−1∏ i=k+1 (1− ηiλj)2
=η2t + (1− ηtλj)2vt,j
(G.10)
1) If vt+1,j ≤ vt,j , then it naturally follows vt+1,j ≤ max(vt,j , ηt/λj).
2) If vt+1,j > vt,j , denote a , (1 − ηtλj)2, b , η2t , we have vt+1,j = avt,j + b, where a ∈ [0, 1) and b ≥ 0. It follows,
vt+1,j > vt,j
⇒ avt,j + b > vt,j
⇒ vt,j < b
1− a
⇒ vt+1,j = avt,j + b < a · b
1− a + b =
b
1− a
Therefore,
vt+1,j < b
1− a = η2t 1− (1− ηtλj)2 < η2t 1− (1− ηtλj) = ηt λj ≤ max
( vt,j ,
ηt λj
)
where the second inequality is entailed by the fact that 1− ηtλj ∈ [0, 1).
Lemma 11. Letting vt,j be defined as Lemma 10 and index ĩ satisfy λj ∈ [µ · 2ĩ, µ · 2ĩ+1), then vĩ+1,j has the following property
vtĩ+1,j ≤ 15 · ηtĩ+1 λj . (G.11)
Proof. By the fact that (1− x) ≤ exp(−x), we have
vt+1,j ≤ t∑
k=0
exp ( −2
t∑ t′=k+1 ηt′λj ) η2k.
Setting t = t̃i+1 − 1 in above equation, we have
vtĩ+1,j ≤ tĩ+1−1∑ k=0 exp −2 tĩ+1−1∑ t′=k+1 ηt′λj  η2k. (G.12)
Now we bound the variance term. First, we have tĩ+1−1∑ k=0 exp −2 tĩ+1−1∑ t=k+1 ηtλj  η2k = tĩ+1−1∑ k=0 exp −2 tĩ+1−1∑ t=k ηtλj  exp(2ηkλj)η2k ≤ tĩ+1−1∑ k=0 exp −2 tĩ+1−1∑ t=k ηtλj  exp(2λj L ) η2k
≤ exp(2) · tĩ+1−1∑ k=0 exp −2 tĩ+1−1∑ t=k ηtλj  η2k, where the first inequality is because ηk ≤ 1/L. Hence, we can obtain
tĩ+1−1∑ k=0 exp −2 tĩ+1−1∑ t=k+1 ηtλj  η2k ≤ exp(2) ·
tĩ+1−1∑ k=0 exp −2 tĩ+1−1∑ t=k ηtλj  η2k = exp(2) ·
ĩ+1∑ i=1 ti−1∑ k=ti−1 exp −2 tĩ+1−1∑ t=k ηtλj  η2k. Furthermore, combining with Eqn. (G.1) and the condition λj ∈ [µ · 2ĩ, µ · 2ĩ+1), we can obtain
ĩ+1∑ i=1 ti−1∑ k=ti−1 exp −2 tĩ+1−1∑ t=k ηtλj  η2k ≤ ĩ+1∑ i=1 ti−1∑ k=ti−1 exp −2 tĩ+1−1∑ t=k ηtµ · 2ĩ  η2k
(G.1) ≤ ĩ+1∑ i=1 ti−1∑ k=ti−1 exp −2 ĩ+1∑ j=i+1 2ĩ−j+1 ln αj αj−1 − 2 · 2ĩ−i+1 ln αi αi−1 + 2i−1µ(k − ti−1)  η2k =
ĩ+1∑ i=1 ti−1∑ k=ti−1 exp  ĩ+1∑ j=i+1 2ĩ−j+2 ln αj−1 αj + 2ĩ−i+2 ln αi−1 + 2 i−1µ(k − ti−1) αi  η2k =
ĩ+1∑ i=1 ti−1∑ k=ti−1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2 · (αi−1 + 2i−1µ(k − ti−1) αi )2ĩ−i+2 η2k
= ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2 ti−1∑ k=ti−1 ( αi−1 + 2 i−1µ(k − ti−1) αi )2ĩ−i+2 η2k
(3.3) = ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2 · ti−1∑ k=ti−1 ( αi−1 + 2 i−1µ(k − ti−1) αi )2ĩ−i+2 L+ µ i−1∑ j=1 ∆j2 j−1 + 2i−1µ(k − ti−1) −2
(G.2) = ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2 · ti−1∑ k=ti−1 ( αi−1 + 2 i−1µ(k − ti−1) αi )2ĩ−i+2 ( αi−1 + 2 i−1µ(k − ti−1) )−2
= ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i · ti−1∑ k=ti−1 (αi−1 + 2 i−1µ(k − ti−1))2 ĩ−i+2 ( αi−1 + 2 i−1µ(k − ti−1) )−2
= ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2α−2ĩ−i+2i ti−1∑ k=ti−1 (αi−1 + 2 i−1µ(k − ti−1))2 ĩ−i+2−2
(G.3) ≤ 2 · 1 2ĩ+1µ · ĩ+1∑ i=1  ĩ+1∏ j=i+1 ( αj−1 αj )2ĩ−j+2(α2ĩ−i+2−1i − α2ĩ−i+2−1i−1 )α−2ĩ−i+2i (G.4) ≤ 2 · 1
2ĩ+1µ · α−1 ĩ+1 ≤ 2 · ηtĩ+1 λj ,
where the last inequality is because of the condition λj ∈ [µ · 2ĩ, µ · 2ĩ+1) and the definition of αi. Therefore, we have
vtĩ+1,j ≤ 2 exp(2) · ηtĩ+1 λj ≤ 15 · ηtĩ+1 λj .
Lemma 12. Let objective function f(x) be quadratic and Assumption (1.7) hold. Running SGD for T -steps starting from w0 and a learning rate sequence {ηt}Tt=1 defined in Eqn. (3.3), the final iterate wT+1 satisfies
E [ (wT+1 − w∗)>H(wT+1 − w∗) ] ≤(w0 − w∗)>H(w0 − w∗) · exp ( −2µ
T∑ k=0 ηk
)
+ 15σ2µ Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 .
Proof. The target of this lemma is to obtain the explicit form to bound the variance term. By the definition of vt+1,j in Lemma 10, we can obtain that
T∑ τ=0 E [ η2τn > τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ] (F.7) ≤ σ2
d∑ j=1 λ2j · vT+1,j
(G.9) ≤ σ2 d∑ j=1 λ2j ·max ( vtĩ+1+1,j , ηtĩ+1+1 λj ) (G.11) ≤ σ2
d∑ j=1 λ2j ·max ( 15 · ηtĩ+1+1 λj , ηtĩ+1+1 λj )
=15σ2 d∑ j=1 λj · ηtĩ+1+1 ≤ 15σ 2 d∑ j=1 λj · ηtĩ+1 ≤ 15σ 2µ Imax−1∑ ĩ=0 2ĩ+1sĩ · ηtĩ+1 ,
where the last inequality is because λj ∈ [2ĩµ, 2ĩ+1µ) and there are sĩ such λj’s lie in this range. By Eqn. (3.3), we have
ηtĩ+1 = 1 L+ µ ∑ĩ+1 j=1 ∆j2 j−1 . (G.13)
Therefore, we have
T∑ τ=0 E [ η2τn > τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ] ≤ 15σ2µ Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 . (G.14)
Combining with Lemma 4 and Lemma 6, we can obtain that
E [ (wT+1 − w∗)>H(wT+1 − w∗) ] ≤(w0 − w∗)>H(w0 − w∗) · exp ( −2µ
T∑ k=0 ηk
)
+ 15σ2µ Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 .
Lemma 13. For ∀t ≥ 0, the learning rate sequence {ηt}Tt=1 defined in Eqn. (3.3) satisfies
ηt ≤ 1
L+ µt (G.15)
Proof. For ∀t ≥ 0, there ∃i ≥ 1, where t ∈ [ti−1, ti). Given the form defined in Eqn. (3.3), we have,
ηt = 1 L+ µ ∑i−1 j=1 ∆j2 j−1 + 2i−1µ(t− ti−1)
≤ 1 L+ µ ∑i−1 j=1 ∆j + µ(t− ti−1)
(3.4) =
1 L+ µ ∑i−1 j=1(tj − tj−1) + µ(t− ti−1)
= 1
L+ µ(ti−1 − t0) + µ(t− ti−1)
= 1
L+ µ(t− t0)
= 1
L+ µt
Lemma 14. Let objective function f(x) be quadratic and Assumption (1.7) hold. Running SGD for T -steps starting from w0 and a learning rate sequence {ηt}Tt=1 defined in Eqn. (3.3), the final iterate wT+1 satisfies
E [ (wT+1 − w∗)>H(wT+1 − w∗) ] ≤(w0 − w∗)>H(w0 − w∗) · κ2
∆21
+ 15σ2µ Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 .
Proof. The target of this lemma is to obtain the explicit form to bound the bias term.
First, by Eqn. (G.1) and the condition λj ∈ [µ · 2ĩ, µ · 2ĩ+1), we have
exp ( −2λj
T∑ k=0 ηk
) ≤ exp −2 tĩ+1−1∑ k=0 ηkλj  ≤ exp −2 ĩ+1∑ i=1 1 2i−1µ ln αi αi−1 λj  ≤ exp
− ĩ+1∑ i=1 2ĩ−i+2 ln αi αi−1  = ĩ+1∏ i=1 ( αi−1 αi )2ĩ−i+2
≤ ĩ+1∏ i=1 ( αi−1 αi )2 = ( α1 αĩ+1 )2 = L2 · η2tĩ+1
(G.16)
For λj = µ, since µ ∈ [µ · 2ĩ, µ · 2ĩ+1) for ĩ = 0, it follows,
exp ( −2µ
T∑ k=0 ηk ) ≤ L2 · η2t1 (G.15) ≤ ( L L+ µt1 )2 (3.4) = ( L L+ µ∆1 )2 ≤ ( L µ∆1 )2 = κ2
∆21 (G.17)
Combining with Lemma 12, we obtain that,
E [ (wT+1 − w∗)>H(wT+1 − w∗) ] ≤(w0 − w∗)>H(w0 − w∗) · κ2
∆21
+ 15σ2µ Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 .

heading: G.1 PROOF OF LEMMA 1
text: Lemma 1. Let objective function f(x) be quadratic and Assumption (1.7) hold. Running SGD for T -steps starting fromw0 and a learning rate sequence {ηt}Tt=1 defined in Eqn. (3.3), the final iterate wT+1 satisfies
E [f(wT+1)− f(w∗)] ≤(f(w0)− f(w∗)) · κ2
∆21 +
15
2 · σ2µ Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 .
(3.6)
Proof. For ∀t ≥ 0, we have
f(wt)− f(w∗) (1.2) = E
[ 1
2 w>t H(ξ)wt − b(ξ)>wt
] − E [ 1
2 w>∗ H(ξ)w∗ − b(ξ)>w∗ ] = ( 1
2 w>t E[H(ξ)]wt − E[b(ξ)]>wt
) − ( 1
2 w>∗ E[H(ξ)]w∗ − E[b(ξ)]>w∗ ) = ( 1
2 w>t Hwt − b>wt
) − ( 1
2 w>∗ Hw∗ − b>w∗ ) (1.4) = ( 1
2 w>t Hwt − b>wt
) − ( 1 2 b> ( H> )−1 b− b>H−1b ) = ( 1
2 w>t Hwt − b>wt
) − ( 1
2 b>H−1b− b>H−1b ) = 1
2 w>t Hwt − b>wt +
1 2 b>H−1b
= 1
2 w>t Hwt −
1 2 b>wt − 1 2 b>wt + 1 2 b>H−1b
= 1
2 w>t Hwt −
1 2 w>t b− 1 2 b>wt + 1 2 b>H−1b
= 1
2 w>t Hwt −
1 2 w>t Hw∗ − 1 2 w>∗ Hwt + 1 2 w>∗ Hw∗
= 1
2 (wt − w∗)>H(wt − w∗),
where the 5th equality is entailed by the fact that H> = H is a symmetric matrix, and the 9th equality uses both H> = H and Eqn 1.4.
Combine the above result with Lemma 14, we obtain that
E [f(wT+1)− f(w∗)] ≤(f(w0)− f(w∗)) · κ2
∆21 +
15
2 · σ2µ Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 .

heading: G.2 PROOF OF THEOREM 1
text: Theorem 1. Let objective function f(x) be quadratic and Assumption (1.7) hold. Running SGD for T -steps starting from w0, a learning rate sequence {ηt}Tt=1 defined in Eqn. (3.3) and ∆i defined in Eqn. (3.5), the final iterate wT+1 satisfies
E [f(wT+1)− f(w∗)] ≤(f(w0)− f(w∗)) · κ2 ·
(∑Imax−1 i=0 √ si )2 s0T 2 + 15 (∑Imax−1 i=0 √ si )2 T · σ2.
Proof. We have
µ · Imax−1∑ ĩ=0
2ĩ+1sĩ L+ µ ∑ĩ+1 j=1 ∆j2 j−1 <µ · Imax−1∑ ĩ=0 2ĩ+1sĩ µ2ĩ∆ĩ+1 (3.5) = 2 Imax−1∑ ĩ=0 sĩ√ sĩ∑Imax−1
i=0
√ si · T
= 2 T · Imax−1∑ i=0 √ si · Imax−1∑ ĩ=0 √ sĩ =
2 (∑Imax−1
i=0
√ si )2 T .
Combining with Lemma 1 and the definition of ∆1, we can obtain that
E [f(wT+1)− f(w∗)] ≤(f(w0)− f(w∗)) · κ2 ·
(∑Imax−1 i=0 √ si )2 s0T 2 + 15 (∑Imax−1 i=0 √ si )2 T · σ2.

heading: G.3 PROOF OF COROLLARY 2
text: Corollary 2. Given the same setting as in Theorem 1, when Hessian H’s eigenvalue distribution p(λ) satisfies “power power law”, i.e.
p(λ) = 1 Z · exp(−α(ln(λ)− ln(µ))) = 1 Z · (µ λ )α (3.7)
for some α > 1, where Z = ∫ L µ
(µ/λ)αdλ, there exists a constant C(α) which only depends on α, such that the final iterate wT+1 satisfies
E [f(wT+1)− f(w∗)] ≤ ( (f(w0)− f(w∗)) · κ2
T 2 + dσ2 T
) · C(α).
Proof. According to Theorem 1,
E [f(wT+1)− f(w∗)]
≤(f(w0)− f(w∗)) · κ2 ·
(∑Imax−1 i=0 √ si )2 s0T 2 + 15 (∑Imax−1 i=0 √ si )2 T · σ2
=(f(w0)− f(w∗)) · κ2
T 2 ·
(∑Imax−1 i=0 √ si )2 s0 + dσ2 T · 15 (∑Imax−1 i=0 √ si )2 d .
The key terms here are C1 , (∑Imax−1
i=0
√ si )2 /s0 and C2 , 15 (∑Imax−1 i=0 √ si )2 /d. As long as
we can bound both terms with a constant C(α), the corollary will be directly proved.
1) If κ < 2, then there is only one interval with s0 = d. By setting C(α) = max(C1, C2) = 15, this completes the proof.
2) If κ ≥ 2, then bounding C1 and C2 be done by computing the value of si under power power law. For all interval i except the last interval, we have,
si d (3.1) = #λj ∈ [µ · 2i, µ · 2i+1) = ∫ µ·2i+1 µ·2i p(λ)dλ
(3.7) = ∫ µ·2i+1 µ·2i 1 Z · (µ λ )α dλ = 1 Z · µα · ∫ µ·2i+1 µ·2i λ−αdλ
= (∫ L µ (µ λ )α dλ )−1 · µα · ∫ µ·2i+1 µ·2i λ−αdλ = (∫ L µ λ−αdλ )−1 · ∫ µ·2i+1 µ·2i λ−αdλ
=
( λ1−α
1− α ∣∣∣L µ )−1 · ( λ1−α 1− α ∣∣∣µ·2i+1 µ·2i ) = ( λ1−α ∣∣∣L µ )−1 · ( λ1−α ∣∣∣µ·2i+1 µ·2i ) = ( L1−α − µ1−α
)−1 · (µ1−α · (2i+1)1−α − µ1−α · (2i)1−α) = µ1−α · ( 2i+1 )1−α − µ1−α · (2i)1−α L1−α − µ1−α = ( 2i+1 )1−α − (2i)1−α κ1−α − 1
=2i(1−α) · 2 1−α − 1 κ1−α − 1
Therefore, we have
si = d · 2i(1−α) · 21−α − 1 κ1−α − 1 = d · 2 1−α − 1 κ1−α − 1 · 2i(1−α) (G.18)
holds for all interval i except the last interval i′ = Imax − 1 = log2 κ − 1 > 0. This last interval may not completely covers [µ · 2i′ , µ · 2i′+1) due to the boundary truncated by L, but we still have
si′ ≤ d · 21−α − 1 κ1−α − 1 · 2i ′(1−α)
It follows,( Imax−1∑ i=0 √ si )2 ≤d · 2 1−α − 1 κ1−α − 1 · ( Imax−1∑ i=0 √ 2i(1−α) )2 = d · 2 1−α − 1 κ1−α − 1 · ( Imax−1∑ i=0 2i(1−α)/2 )2
=d · 2 1−α − 1 κ1−α − 1 · ( Imax−1∑ i=0 ( 1 2 )i(α−1)/2)2
≤d · 2 1−α − 1 κ1−α − 1 · ( ∞∑ i=0 ( 1 2 )i(α−1)/2)2 = d · 2 1−α − 1 κ1−α − 1 · ( ∞∑ i=0 ( 1 2(α−1)/2 )i)2
=d · 2 1−α − 1 κ1−α − 1 ·
( 1
1− 1 2(α−1)/2
)2 = d · 2
1−α − 1 κ1−α − 1
· (
1
1− 2(1−α)/2
)2 .
Thus,
C1 =
(∑Imax−1 i=0 √ si )2 s0 ≤ d · 2 1−α−1 κ1−α−1 · ( 1 1−2(1−α)/2 )2 d · 21−α−1κ1−α−1 · 20(1−α) = ( 1 1− 2(1−α)/2 )2
C2 = 15 (∑Imax−1 i=0 √ si )2 d ≤ 15 d · d · 2 1−α − 1 κ1−α − 1 · (
1
1− 2(1−α)/2 )2 =15 · 1− ( 1 2
)α−1 1− ( 1 κ )α−1 · ( 11− 2(1−α)/2 )2
≤15 · (
1
1− 2(1−α)/2
)2 .
Here the last inequality for C2 is entailed by κ ≥ 2 and α > 1. By setting C(α) = max(C1, C2) = 15 · (
1 1−2(1−α)/2
)2 , we obtain
E [f(wT+1)− f(w∗)]
≤(f(w0)− f(w∗)) · κ2
T 2 ·
(∑Imax−1 i=0 √ si )2 s0 + dσ2 T · 15 (∑Imax−1 i=0 √ si )2 d .
=(f(w0)− f(w∗)) · κ2
T 2 · C1 +
dσ2
T · C2 ≤ ( (f(w0)− f(w∗)) · κ2
T 2 + dσ2 T
) · C(α).

heading: G.4 PROOF OF THEOREM 4
text: Theorem 4. Let objective function f(x) be quadratic. We run SGD for T -steps starting fromw0 and a step decay learning rate sequence {ηt}Tt=1 defined in Algorithm 1 of Ge et al. (2019) with η1 = 1/L. As long as (1) H is diagonal, (2) The equality in Assumption (1.7) holds, i.e. Eξ [ ntn > t ] = σ2H and (3) T/ log T ≥ max(κ, 216, 16), the final iterate wT+1 satisfies,
E [f(wT+1)− f(w∗)] = Ω ( dσ2
T · log T
)
Proof. In Ge et al. (2019), step decay has following learning rate sequence:
ηt = 1
2`L if t ∈
[ 1 + T
log T · `, T log T · (`+ 1)
] , (G.19)
where ` = 0, 1, . . . , log T − 1. We first translate the general notations to diagonal cases so that the idea of the proof can be clearer.
Since f(x) is quadratic, according to the proof in Lemma 1,
f(wT+1)− f(w∗) = 1
2 (wT+1 − w∗)>H(wT+1 − w∗)
Furthermore, according to Lemma 4, where Pt = I − ηtH , E [ (wT+1 − w∗)>H(wT+1 − w∗) ] =E [ (w0 − w∗)> · PT . . . P0HP0 . . . PT · (w0 − w∗)
] +
T∑ τ=0 E [ η2τn > τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ] ≥
T∑ τ=0 E [ η2τn > τ · PT . . . Pτ+1HPτ+1 . . . PT · nτ ] =
T∑ τ=0 η2τ T∏ k=τ+1 d∑ j=1 λj(1− ηkλj)2E [ n2τ,j ] =
T∑ τ=0 η2τ T∏ k=τ+1 d∑ j=1 λj(1− ηkλj)2 · λjσ2
=σ2 d∑ j=1 λ2j T∑ τ=0 η2τ T∏ k=τ+1 (1− ηkλj)2
Here the second equality is entailed by the fact that H and Pt are diagonal, and the third equality comes from Eξ [ ntn > t ] = σ2H . Thus, by denoting vj , ∑T τ=0 η 2 τ ∏T k=τ+1(1− ηkλj)2, we have,
E [f(wT+1 − f(w∗)] ≥ σ2
2 d∑ j=1 λ2j T∑ τ=0 η2τ T∏ k=τ+1 (1− ηkλj)2 = σ2 2 d∑ j=1 λ2jvj (G.20)
For each dimension λj , step decay’s bottleneck actually occurs at the first interval ` that satisfies
2` ≥ λj L · 8T log T (G.21)
We first show that interval ` is well-defined for any dimension j.
Since we assume T/ log T ≥ κ, i.e. 1/κ ≥ log T/T , it follows
2` ≥ λj L · 8T log T ≥ µ L · 8T log T = 1 κ · 8T log T ≥ log T T · 8T log T = 8
=⇒ ` ≥ 3
Since interval ` is the first interval satisifies Eqn. (G.21), we also have
2` < λj L · 16T log T , (G.22)
which follows,
2` < λj L · 16T log T ≤ 16T log T ≤ T
=⇒ ` < log T
The last inequality is because T ≥ T/ log T ≥ 216 ⇒ log T ≥ 16. Since T ≥ 16⇒ log T ≥ 4 > 3, thus such ` is guaranteed to be well-defined. Back to our main analysis, by focusing on interval ` only, we have,
vj = T∑ τ=0 η2τ T∏ k=τ+1 (1− ηkλj)2 ≥ (`+1)· Tlog T∑ τ=`· Tlog T +1 η2τ T∏ k=τ+1 (1− ηkλj)2
≥ (`+1)· Tlog T∑ τ=`· Tlog T +1 η2τ T∏ k=`· Tlog T +1 (1− ηkλj)2 = (`+1)· Tlog T∑ τ=`· Tlog T +1 ( 1 L2` )2 T∏ k=`· Tlog T +1 (1− ηkλj)2
= T log T · ( 1 L2` )2 T∏ k=`· Tlog T +1 (1− ηkλj)2
(G.22) > T
log T ·
( 1
L · λjL · 16T log T )2 T∏ k=`· Tlog T +1 (1− ηkλj)2
= 1 256 · log T T · 1 λ2j · T∏ k=`· Tlog T +1 (1− ηkλj)2
≥ 1 256 · log T T · 1 λ2j · 1− T∑ k=`· Tlog T +1 2ηkλj  = 1 256 · log T T · 1 λ2j · 1− 2λj T∑ k=`· Tlog T +1 ηk  = 1
256 · log T T · 1 λ2j · 1− 2λj log T−1∑ i=` (i+1)· Tlog T∑ k=i· Tlog T +1 ηk  = 1
256 · log T T · 1 λ2j · 1− 2λj log T−1∑ i=` (i+1)· Tlog T∑ k=i· Tlog T +1 1 L2i  = 1
256 · log T T · 1 λ2j ·
( 1− 2λj
log T−1∑ i=` T log T · 1 L2i
)
≥ 1 256 · log T T · 1 λ2j · ( 1− 2λj · T log T · 1 L2`−1 ) (G.21) ≥ 1
256 · log T T · 1 λ2j ·
( 1− 4λj · T
log T · 1 L · λjL · 8T log T
)
= 1 256 · log T T · 1 λ2j · 1 2
= 1 512 · log T T · 1 λ2j
Here the first inequality is obtained by focusing variance generated in interval ` only. The second inequality utilizes τ + 1 ≥ ` · T/ log T . The fourth inequality is entailed by (1 − a1)(1 − a2) = 1 − a1 − a2 + a1a2 ≥ 1 − a1 − a2 for ∀a1, a2 ∈ [0, 1], where by mathematical induction, we can extend this inequality for more terms ∏n i=1(1 − ai) ≥ 1 − ∑n i=1 ai as long as ∑n i=1 ai ≤ 1. The
fifth inequality comes from ∑log T−1 i=` 1/2 i ≤ ∑∞ i=` 1/2 i = 1/2`−1.
By combining with Eqn. (G.20), we have
E [f(wT+1 − f(w∗)] ≥ σ2
2 d∑ j=1 λ2jvj ≥ σ2 2 d∑ j=1 λ2j ·
( 1
512 · log T T · 1 λ2j
) = 1
1024 · dσ
2
T · log T
= Ω
( dσ2
T · log T
)

heading: G.5 THE REASON OF USING ASSUMPTION (1.7)
text: In all of our analysis, we employ assumption (1.7) Eξ [ ntn > t ] σ2H where nt = Hwt − b− (H(ξ)wt − b(ξ))
which is the same as the one in Appendix C, Theorem 13 of Ge et al. (2019). This key theorem is the major difference between our work and Ge et al. (2019), which directly entails its main theorem by instantiating σ with specific values in its assumptions.
On the other hand, it is possible to use the assumptions in Ge et al. (2019); Bach & Moulines (2013); Jain et al. (2016) instead of our assumption (1.7) for least square regression:
min w f(w) where f(w) ,
1 2 E(x,y)∼D
[ (y − w>x)2 ] (G.23)
y = w>∗ x+ with satisfying E(x,y)∼D [ 2xx> ] σ2H for ∀(x, y) ∼ D (G.24)
E [ ||x||2xx> ] R2H (G.25)
By combining our Lemma 1 and assumption (1.7) with Lemma 5, Lemma 8 and Lemma 9 in Ge et al. (2019), one can obtain similar results in this paper with their assumptions. For simplicity, we just use assumption (1.7) here.

Review summaries:
The authors propose Eigencurve, a new approach to learning rate scheduling that utilizes information form the eigenvalues of the Hessian. They show that this scheduler obtains the minimax optimal rate on the noisy quadratic problem. Empirically, this scheduler demonstrates faster convergence on CIFAR-10 and ImageNet, especially when the number of epochs is small.

The paper studies convergence rates of SGD with different stepsize schemes, in the context of linear regression. For convergence of the last iterate of SGD, the best known result still misses a $\\log T$ factor compared to the minimax rate. This work aims to fill this gap with an improved stepsize scheme that utilizes the eigenvalue distribution of the Hessian. When the true Hessian is known, the proposed method successfully fills the gap, provably and in the sense of the worst problem instances. When the true Hessian is hard to know, practical variants are also proposed, and are shown to be comparable to the state-of-the-art stepsize schemes in standard deep neural network benchmarks.

This paper studies the convergence of SGD with different learning rate schedules and aims to achieve minimax optimal convergence rates on quadratic objectives. To this end, this work proposes a new learning rate schedule (named Eigencurve) based on the Hessian spectrum and provides an optimal last-iterate convergence rate. The proposed Eigencurve gives rise to slightly improved performance on image classification tasks with deep neural networks. In addition, the introduced schedule is similar to the popular cosine learning rate schedule on some problems, which, to some extent, justifies the effectiveness of the cosine schedule.

# === Update === #  I appreciate the thorough response to my review provided by the authors. I have looked through the newest version of the submission and find it to be much improved. The inclusion of repeats for the train loss / test accuracy comparison is great to see and the updated figures are mostly fixed. The experiments now make a meaningful argument for Eigencurve as a (non)-convex optimization method.  I have increased my score to 6 to reflect the improvements and the author response.   # ====== #  This submission proposes a new step-size schedule for stochastic gradient descent which leverages the eigenvalue spectrum of Hessian to speed-up convergence.  The approach, dubbed Eigencurve, is shown to achieve the minimax optimal convergence rate for stochastic gradient descent on quadratic functions under an additional condition that the eigenspectrum of the Hessian decays according to a power law. When this decay condition is not satisfied, Eigencurve still improves upon the popular step-decay schedule; in this case, it is sub-optimal by a factor of $\\log(\\kappa)$ rather than the $\\log(T)$ factor of step-decay. A lower bound is also provided verifying that the $\\log(T)$ sub-optimality of step-decay is tight and cannot be improved for quadratics. The submission concludes with an empirical investigation of Eigencurve for (non-convex) optimization of several popular neural network architectures on the CIFAR-10 and ImageNet datasets.