Year: 2020

Abstract:
The Heavy Ball Method (Polyak, 1964), proposed by Polyak over five decades ago, is a first-order method for optimizing continuous functions. While its stochastic counterpart has proven extremely popular in training deep networks, there are almost no known functions where deterministic Heavy Ball is provably faster than the simple and classical gradient descent algorithm in non-convex optimization. The success of Heavy Ball has thus far eluded theoretical understanding. Our goal is to address this gap, and in the present work we identify two non-convex problems where we provably show that the Heavy Ball momentum helps the iterate to enter a benign region that contains a global optimal point faster. We show that Heavy Ball exhibits simple dynamics that clearly reveal the benefit of using a larger value of momentum parameter for the problems. The first of these optimization problems is the phase retrieval problem, which has useful applications in physical science. The second of these optimization problems is the cubic-regularized minimization, a critical subroutine required by Nesterov-Polyak cubic-regularized method (Nesterov & Polyak (2006)) to find second-order stationary points in general smooth non-convex problems.

heading: 1 INTRODUCTION
text: Poylak’s Heavy Ball method (Polyak (1964)) has been very popular in modern non-convex optimization and deep learning, and the stochastic version (a.k.a. SGD with momentum) has become the de facto algorithm for training neural nets. Many empirical results show that the algorithm is better than the standard SGD in deep learning (see e.g. Hoffer et al. (2017); Loshchilov & Hutter (2019); Wilson et al. (2017); Sutskever et al. (2013)), but there are almost no corresponding mathematical results that show a benefit relative to the more standard (stochastic) gradient descent. Despite its popularity, we still have a very poor justification theoretically for its success in non-convex optimization tasks, and Kidambi et al. (2018) were able to establish a negative result, showing that Heavy Ball momentum cannot outperform other methods in certain problems. Furthermore, even for convex problems it appears that strongly convex, smooth, and twice differentiable functions (e.g. strongly convex quadratic functions) are one of just a handful examples for which a provable speedup over standard gradient descent can be shown (e.g (Lessard et al., 2016; Goh, 2017; Ghadimi et al., 2015; Gitman et al., 2019; Loizou & Richtárik, 2017; 2018; Gadat et al., 2016; Scieur & Pedregosa, 2020; Sun et al., 2019; Yang et al., 2018a; Can et al., 2019; Liu et al., 2020; Sebbouh et al., 2020; Flammarion & Bach, 2015)). There are even some negative results when the function is strongly convex but not twice differentiable. That is, Heavy Ball momentum might lead to a divergence in convex optimization (see e.g. (Ghadimi et al., 2015; Lessard et al., 2016)). The algorithm’s apparent success in modern non-convex optimization has remained quite mysterious.
In this paper, we identify two non-convex optimization problems for which the use of Heavy Ball method has a provable advantage over vanilla gradient descent. The first problem is phase retrieval. It has some useful applications in physical science such as microscopy or astronomy (see e.g. (Candés et al., 2013), (Fannjiang & Strohmer, 2020), and (Shechtman et al., 2015)). The objective is
minw∈Rd f(w) := 1
4n ∑n i=1 ( (x>i w) 2 − yi )2 , (1)
where xi ∈ Rd is the design vector and yi = (x>i w∗)2 is the label of sample i. The goal is to recover w∗ up to the sign that is not recoverable (Candés et al., 2013). Under the Gaussian design
setting (i.e. xi ∼ N(0, Id)), it is known that the empirical risk minimizer (1) is w∗ or −w∗, as long as the number of samples n exceeds the order of the dimension d (see e.g. Bandeira et al. (2014)). Therefore, solving (1) allows one to recover the desired vectorw∗ ∈ Rd up to the sign. Unfortunately the problem is non-convex which limits our ability to efficiently find a minimizer. For this problem, there are many specialized algorithms that aim at achieving a better computational complexity and/or sample complexity to recoverw∗ modulo the unrecoverable sign (e.g. (Cai et al., 2016; Candés & Li, 2014; Candés et al., 2015; 2013; Chen & Candés, 2017; Duchi & Ruan, 2018; Ma et al., 2017; 2018; Netrapalli et al., 2013; Qu et al., 2017; Tu et al., 2016; Wang et al., 2017a;b; Yang et al., 2018b; Zhang et al., 2017a;b; Zheng & Lafferty, 2015)). Our goal is not about providing a state-of-the-art algorithm for solving (1). Instead, we treat this problem as a starting point of understanding Heavy Ball momentum in non-convex optimization and hope for getting some insights on why Heavy Ball (Algorithm 1 and 2) can be faster than the vanilla gradient descent in non-convex optimization and deep learning in practice. If we want to understand why Heavy Ball momentum leads to acceleration for a complicated non-convex problem, we should first understand it in the simplest possible setting.
We provably show that Heavy Ball recovers the desired vector w∗, up to a sign flip, given a random isotropic initialization. Our analysis divides the execution of the algorithm into two stages. In the first stage, the ratio of the projection of the current iterate wt on w∗ to the projection of wt on the perpendicular component keeps growing, which makes the iterate eventually enter a benign region which is strongly convex, smooth, twice differentiable, and contains a global optimal point. Therefore, in the second stage, Heavy Ball has a linear convergence rate. Furthermore, up to a value, a larger value of the momentum parameter has a faster linear convergence than the vanilla gradient descent in the second stage. Yet, most importantly, we show that Heavy Ball momentum also has an important role in reducing the number of iterations in the first stage, which is when the iterate might be in a non-convex region. We show that the higher the momentum parameter β, the fewer the iterations spent in the first stage (see also Figure 1). Namely, momentum helps the iterate to enter a benign region faster. Consequently, using a non-zero momentum parameter leads to a speedup over the standard gradient descent (β = 0). Therefore, our result shows a provable acceleration relative to the vanilla gradient descent, for computing a global optimal solution in non-convex optimization.
The second of these is solving a class of cubic-regularized problems,
minw f(w) := 1 2w >Aw + b>w + ρ3‖w‖ 3, (2)
where the matrixA ∈ Rd×d is symmetric and possibly indefinite. Problem (2) is a sub-routine of the Nesterov-Polyak cubic-regularized method (Nesterov & Polyak (2006)), which aims to minimize a non-convex objective F (·) by iteratively solving
wt+1 = arg minw∈Rd{∇F (wt)>(w − wt) + 12 (w − wt) >∇2F (w)(w − wt) + ρ3‖w − wt‖ 3}, (3) With some additional post-processing, the iterate wt converges to an ( g, h) second order stationary point, defined as {w : ‖∇f(w)‖ ≤ g and ∇2f(w) − hId} for any small g, h > 0. However, their algorithm needs to compute a matrix inverse to solve (2), which is computationally expensive when the dimension is high. A very recent result due to Carmon & Duchi (2019) shows that vanilla gradient descent approximately finds the global minimum of (2) under mild conditions, which only needs a Hessian-vector product and can be computed in the same computational complexity as computing gradients (Pearlmutter, 1994), and hence is computationally cheaper than the matrix inversion of the Hessian. Our result shows that, similar to the case of phase retrieval, the use of Heavy Ball momentum helps the iterate to enter a benign region of (3) that contains a global optimal solution faster, compared to vanilla gradient descent. For certain non-convex problems, e.g. dictionary learning (Sun et al., 2015), matrix completion (Chi et al., 2019), robust PCA (Ge et al., 2017), and learning a neural network (Ge et al. (2019); Bai & Lee (2020)), where it suffices to find a second-order stationary point, our result consequently might have application.
To summarize, our theoretical results of the two non-convex problems provably show the benefit of using Heavy Ball momentum. Compared to the vanilla gradient descent, the use of momentum helps to accelerate the optimization process. The key to showing the acceleration in getting into benign regions of these problems is a family of simple dynamics due to Heavy Ball momentum. We will argue that the simple dynamics are not restricted to the two main problems considered in this paper. Specifically, the dynamics also naturally arise when solving the problem of top eigenvector computation (Golub & Loan, 1996) and the problem of saddle points escape (e.g. (Jin et al., 2017;
Wang et al., 2020)), which might imply the broad applicability of the dynamics for analyzing Heavy Ball in non-convex optimization.
Algorithm 1: Heavy Ball method (Polyak, 1964) (Equivalent version 1)
1: Required: step size η and momentum parameter β ∈ [0, 1]. 2: Init: w0 = w−1 ∈ Rd 3: for t = 0 to T do 4: Update iterate wt+1 = wt− η∇f(wt)+β(wt−wt−1). 5: end for
Algorithm 2: Heavy Ball method (Polyak, 1964) (Equivalent version 2)
1: Required: step size η and momentum parameter β ∈ [0, 1]. 2: Init: w0 ∈ Rd and m−1 = 0. 3: for t = 0 to T do 4: Update momentum mt := βmt−1 +∇f(wt). 5: Update iterate wt+1 := wt − ηmt. 6: end for

heading: 2 MORE RELATED WORKS
text: Heavy Ball (HB): HB has two exactly equivalent presentations in the literature (see Algorithm 1 and 2). Given the same initialization, both algorithms generate the same sequence of {wt}. In Algorithm 2, we note that the momentum mt can be written as mt = ∑t s=0 β
t−s∇f(ws) and can be viewed as a weighted sum of gradients. As we described in the opening paragraph, there is little theory of showing a provable acceleration of the method in non-convex optimization. The only exception that we are aware of is (Wang et al., 2020). They show that HB momentum can help to escape saddle points faster and find a second-order stationary point faster for smooth non-convex optimization. They also observed that stochastic HB solves (1) and that using higher values of the momentum parameter β leads to faster convergence. However, while their work focused on the stochastic setting, their main result required some assumptions on the statistical properties of the sequence of observed gradients; it is not clear whether these would hold in general. In appendix A, we provide a more detailed literature review of HB. To summarize, current results in the literature imply that we are still very far from understanding deterministic HB in non-convex optimization, let alone understanding the success of stochastic HB in deep learning. Hence, this work aims to make progress on a simple question: can we give a precise advantage argument for the acceleration effect of Heavy Ball in the deterministic setting?
Phase retrieval: The optimization landscape of problem (1) and its variants has been studied by (Davis et al., 2018; Soltanolkotabi, 2014; Sun et al., 2016; White et al., 2016), which shows that as long as the number of samples is sufficiently large, it has no spurious local optima. We note that the problem can also be viewed as a special case of matrix sensing (e.g. Li et al. (2018); Gunasekar et al. (2017); Li & Lin (2020); Li et al. (2019); Gidel et al. (2019); You et al. (2020)); in Appendix A, we provide a brief summary of matrix sensing. For solving phase retrieval, Mannellia et al. (2020) study gradient flow, while Chen et al. (2018) show that the standard gradient descent with a random initialization like Gaussian initialization solves (1) and recovers w∗ up to the sign. Tan & Vershynin (2019) show that online gradient descent with a simple random initialization can converge to a global optimal point in an online setting where fresh samples are required for each step. In this paper, we show that Heavy Ball converges even faster than the vanilla gradient descent. Zhou et al. (2016) propose leveraging Nesterov’s momentum to solve phase retrieval. However, their approach requires delicate and computationally expensive initialization like spectral initialization so that the initial point is already within the neighborhood of a minimizer. Similarly, Xiong et al. (2018; 2020) show local convergence of Nesterov’s momentem and Heavy Ball momentum for phase retrieval, but require the initial point to be in the neighborhood of an optimal point. Jin et al. (2018) propose an algorithm that uses Nesterov’s momentum together as a subroutine with perturbation for finding a second-order stationary point, which could be applied for solving phase retrieval. Compared to (Zhou et al., 2016; Jin et al., 2018; Xiong et al., 2018; 2020), we consider directly applying gradient descent with Heavy Ball momentum (i.e. HB method) to the objective function with simple random initialization, e.g. Gaussian initialization, which is what people do in practice and is what we want to understand. The goals of the works are different. Finally, we note that there are some efforts in integrating the technique of generative models and phase retrieval, which could help the task of image recovery (e.g. Hand et al. (2018)). Phase retrieval might also be a good entry point of understanding some observations in optimization and neural net training (e.g. Mannellia et al. (2020)).
f(w1)
≥ 0.5}, i.e. if the relative change of objective value
compared to the initial value has been increased to 50%. Algorithm 3 and Algorithm 4 in Appendix H describe the procedures. All the lines are obtained by initializing the iterate at the same pointw0 ∼ N (0, Id/(10000d)) and using the same step size η = 5× 10−4. Here we set w∗ = e1 and sample xi ∼ N (0, Id) with dimension d = 10 and number of samples n = 200. We see that the higher the momentum parameter β, the faster the algorithm enters the linear convergence regime. (a): Objective value (1) vs. iteration t. We see that the higher the momentum parameter β, the faster the algorithm enters the linear convergence regime. (b): The size of projection of wt on w∗ over iterations (i.e. |w‖t | vs. t), which is non-decreasing throughout the iterations until reaching an optimal point (here, ‖w∗‖ = 1). (c): The size of the perpendicular component over iterations (i.e. ‖w⊥t ‖ vs. t), which is increasing in the beginning and then it is decreasing towards zero after some point. We see that the slope of the curve corresponding to a larger momentum parameter β is steeper than that of a smaller one, which confirms Lemma 1 and Lemma 3.

heading: 3 PHASE RETRIEVAL

heading: 3.1 PRELIMINARIES
text: Following the works of Candés et al. (2013); Chen et al. (2018), we assume that the design vectors {xi} (which are known a priori) are from Gaussian distribution xi ∼ N(0, Id). Furthermore, without loss of generality, we assume that w∗ = e1 (so that ‖w∗‖ = 1), where e1 is the standard unit vector whose first element is 1. We also denote w‖t := wt[1] and wt,⊥ := [wt[2], . . . , wt[d]] > . That is, w ‖ t is the projection of the current iterate wt on w∗, while w ⊥ t is the perpendicular component. Throughout the paper, the subscript t is an index of the iterations while the subscript i is an index of the samples.
Before describing the main results, we would like to provide a preliminary analysis to show how momentum helps. Applying gradient descent with Heavy Ball momentum (Algorithm 1) to objective (1), we see that the iterate is generated according to wt+1 = wt − η 1n ∑n i=1 ( (x>i wt)
3 − (x>i wt)yi ) xi +β(wt−wt−1). On the other hand, the population counterpart (i.e. when the number of samples n is infinite) of the update rule turns out to be the key to understanding momentum. The population gradient∇F (w) := Ex∼N(0,In)[∇f(w)] is (proof is available in appendix B)
∇F (w) = ( 3‖w‖2 − 1 ) w − 2(w>∗ w)w∗. (4)
Using the population gradient (4), we have the population update, wt+1 = wt−η∇F (wt)+β(wt− wt−1), which can be decomposed as follows:
w ‖ t+1 = ( 1 + 3η(1− ‖wt‖2) ) w ‖ t + β(w ‖ t − w ‖ t−1)
w⊥t+1 = ( 1 + η(1− 3‖wt‖2) ) w⊥t + β(w ⊥ t − w⊥t−1).
(5)
Assume that the random initialization satisfies ‖w0‖2 13 . From the population recursive system (5), both the magnitude of the signal component w‖t and the perpendicular component w ⊥ t grow exponentially in the first few iterations.
Lemma 1. For a positive number θ and the momentum parameter β ∈ [0, 1], if a non-negative sequence {at} satisfies a0 ≥ a−1 > 0 and that for all t ≤ T ,
at+1 ≥ (1 + θ)at + β(at − at−1), (6)
then {at} satisfies at+1 ≥ ( 1 + ( 1 + β1+θ ) θ ) at, (7)
for every t = 1, . . . , T + 1. Similarly, if a non-positive sequence {at} satisfies a0 ≤ a−1 < 0 and that for all t ≤ T , at+1 ≤ (1 + θ)at + β(at − at−1), then {at} satisfies
at+1 ≤ ( 1 + ( 1 + β1+θ ) θ ) at, (8)
for every t = 1, . . . , T + 1.
One can view at in Lemma 1 as the projection of the current iterate wt onto a vector of interest. The lemma says that with a larger value of the momentum parameter β, the magnitude of at is increasing faster. It also implies that if the projection due to vanilla gradient descent satisfies at+1 ≥ (1+θ)at, then the magnitude of the projection only grows faster with the use of Heavy Ball momentum. The dynamics in the lemma are the keys to showing that Heavy Ball momentum accelerates the process of entering a benign (convex) region. The factor β1+θ θ in (7) and (8) represents the contribution due to the use of momentum, and the contribution is larger with a larger value of momentum parameter β. Now let us apply Lemma 1 to the recursive system (5) and pretend that the magnitude of ‖wt‖ was a constant for a moment. Denote θt := 3η(1−‖wt‖2) and θ̃t := η(1− 3‖wt‖2) and notice that θt > θ̃t > 0 when ‖wt‖2 < 13 . We can rewrite the recursive system as
w ‖ t+1 = (1 + θt)w ‖ t + β(w ‖ t − w ‖ t−1) w⊥t+1 = (1 + θ̃t)w ⊥ t + β(w ⊥ t − w⊥t−1).
(9)
Since the above system is in the form of (6), the dynamics (7) and (8) in Lemma 1 suggest that the larger the momentum parameter β, the faster the growth rate of the magnitude of the signal componentw‖t and the perpendicular componentw ⊥ t . Moreover, the magnitude of the signal componentw ‖ t grows faster than that of the perpendicular componentw⊥t . Both components will grow until the size of iterate ‖wt‖2 is sufficiently large (i.e. ‖wt‖2 > 13 ). After that, the magnitude of the perpendicular component w⊥t starts decaying, while |w ‖ t | keeps growing until it approaches 1. Furthermore, we have that the larger the momentum parameter β, the faster the decay rate of the (magnitude of the) perpendicular component w⊥t . In other words, |w ‖ t | converges to 1 and w⊥t converges to 0 quickly. Lemma 3 in Appendix C, which is a counterpart of Lemma 1, can be used to explain the faster decay of the magnitude due to a larger value of the momentum parameter β.
By using the population recursive system (5) and Lemma 1 & 3 as the tool, we obtain a high-level insight on how momentum helps (see also Figure 1). The momentum helps to drive the iterate to enter the neighborhood of a w∗ (or −w∗) faster.

heading: 3.2 MAIN RESULTS
text: We denote dist(wt, w∗) := min{‖wt−w∗‖2, ‖wt+w∗‖2} as the distance between the current iterate wt and w∗, modulo the unrecoverable sign. Note that both±w∗ achieve zero testing errors. Furthermore, as long as the number of samples is sufficiently large, i.e. n & d log d, there exists a constant ϑ > 0 so that the Hessian satisfies ∇2f(w) ϑId for all w ∈ Bζ(±w∗) with high probability, where Bζ(±w∗) represents the balls centered at ±w∗ with a radius ζ (e.g. (Ma et al., 2017)). So in this paper we consider the case that the local strong convexity holds in the neighborhoodBζ(±w∗). We will divide the iterations into two stages. The first stage consists of those iterations that satisfy 0 ≤ t ≤ Tζ , where Tζ is defined as
Tζ := min{t : ||w‖t | − 1| ≤ ζ 2 and ‖w ⊥ t ‖ ≤ ζ 2}, (10)
and ζ > 0 is sufficiently small so that it makes wTζ be in the neighborhood of w∗ or −w∗ which is smooth, twice differentiable, strongly convex, see e.g. (Ma et al., 2017; Soltanolkotabi, 2014).
Observe that if ||w‖t | − 1| ≤ ζ 2 and ‖w ⊥ t ‖ ≤ ζ 2 , we have that dist(wt, w∗) ≤ ‖wt − w∗‖ ≤ ||w‖t | − 1| + ‖w⊥t ‖ ≤ ζ. The second stage consists of those iterations that satisfy t ≥ Tζ . Given that dist(wTζ , w∗) ≤ ζ and that the local strong convexity holds in Bζ(±w∗), the iterate wt would be in a benign region at the start of this stage, which allows linear convergence to a global optimal point. That is, we have that dist(wt, w∗) ≤ (1 − ν)t−Tζζ for all t ≥ Tζ , where 1 > ν > 0 is some number. Since the behavior of the momentum method in the second stage can be explained by the existing results (e.g. Section 3 of Saunders (2018) or Xiong et al. (2020)), the goal is to understand why momentum helps to drive the iterate into the benign region faster.
To deal with the case that only finite samples are available in practice, we will consider some perturbations from the population dynamics (5). In particular, we consider
w ‖ t+1 = ( 1 + 3η(1− ‖wt‖2) + ηξt ) w ‖ t + β(w ‖ t − w ‖ t−1)
w⊥t+1[j] = ( 1 + η(1− 3‖wt‖2) + ηρt,j ) w⊥t [j] + β(w ⊥ t [j]− w⊥t−1[j]),
(11)
where {ξt} and {ρt,j} for 1 ≤ j ≤ d−1 are the perturbation terms. The perturbation terms are used to model the deviation from the population dynamics. In this paper, we assume that there exists a small number cn > 0 such that for all iterations t ≤ Tζ and all j ∈ [d− 1], max{|ξt|, |ρt,j |} ≤ cn, where the value cn should decay when the number of samples n is increasing and cn u 0 when there are sufficiently large number of samples n. Theorem 1. Suppose that the approximated dynamics (11) holds with max{|ξt|, |ρt,j |} ≤ cn for all iterations t ≤ Tζ and all dimensions j ∈ [d − 1], where cn ≤ ζ. Suppose that the local strong convexity holds in Bζ(±w∗). Assume that the initial point w0 satisfies |w‖0 | & 1√d log d and ‖w0‖ < 13 . Set the momentum parameter β ∈ [0, 1]. Assume that the norm of the momentum ‖mt‖ is bounded for all t ≤ Tζ , i.e. ‖mt‖ ≤ cm for some constant cm > 0. If the step size η satisfies η ≤ 1
36 ( 1+ζ ) max{cm,1} , then Heavy Ball (Algorithm 1 & 2) takes at most
Tζ . log d η ( 1 + cηβ )
number of iterations to enter the benign region Bζ(w∗) or Bζ(−w∗), where cη := 11+η/2 ≤ 1. Furthermore, for all t ≥ Tζ , the distance is shrinking linearly for some values of η and β < 1. That is, we have that dist(wt, w∗) ≤ (1− ν)t−Tζζ, for some number 1 > ν > 0.
The theorem states that the number of iterations required for gradient descent to enter the linear convergence is reduced by a factor of (1 + cηβ), which clearly demonstrates that momentum helps to drive the iterate into a benign region faster. The constant cη suggests that the smaller the step size η, the acceleration due to the use of momentum is more evident. The reduction can be about u 1+β for a small η. After Tζ , Heavy Ball has a locally linear convergence to w∗ or −w∗. Specifically, if w ‖ 0 > 0, then we will have that w ‖ t > 0 for all t and that the iterate will converge to w∗; otherwise, we will have w‖t < 0 for all t and that the iterate will converge to −w∗. The proof of Theorem 1 is in Appendix D.
Remark 1: The initial point is required to satisfy 〈w0, w∗〉 & 1√d log d and ‖w0‖ < 1 3 . The first condition can be achieved by generating w0 from a Gaussian distribution or uniformly sampling from a sphere with high probability (see e.g. Chapter 2 of (Blum et al., 2018)). The second condition can then be satisfied by scaling the size appropriately. We note that the condition that the norm of the momentum is bounded is also assumed in (Wang et al., 2020).
Remark 2: Our theorem indicates that in the early stage of the optimization process, the momentum parameter β can be as large as 1, which is also verified in the experiment (Figure 1). However, to guarantee convergence after the iterate is in the neighborhood of a global optimal solution, the parameter β must satisfy β < 1 (Polyak, 1964; Lessard et al., 2016).
Remark 3: The number ν of the local linear convergence rate due to Heavy Ball momentum actually depends on the smoothness constant L and strongly convexity constant µ in the neighborhood of a global solution, as well as the step size η and the momentum parameter β (see e.g. Section 3 of (Saunders, 2018) or the original paper (Polyak, 1964)). By setting the step size η = 4/( √ L+ √ µ)2
and the momentum parameter β = max{|1 − √ ηL|, |1 − √ηµ|}2, ν will depend on the squared root of the condition number √ κ := √ L/µ instead of κ := L/µ, which means that an optimal local convergence rate is achieved (e.g. (Bubeck, 2014)). In general, up to a certain threshold, a larger value of β leads to a faster rate than that of standard gradient descent.

heading: 4 CUBIC-REGULARIZED PROBLEM

heading: 4.1 NOTATIONS
text: We begin by introducing the notations used in this section. For the symmetric but possibly indefinite matrixA, we denote its eigenvalue in the increasing order, where any λ(i) might be negative. We denote the eigen-decomposition ofA asA := ∑d i=1 λ (i)(A)viv > i , where each vi ∈ Rd is orthonormal. We also denote γ := −λ(1)(A), and γ+ = max{γ, 0}, and ‖A‖2 := max{|λ(1)(A)|, |λ(d)(A)|}. For any vectorw ∈ Rd, we denotew(i) as the projection on the eigenvector ofA,w(i) = 〈w, vi〉. Denote w∗ as a global minimizer of the cubic-regularized problem (2) and denoteA∗ := A+ρ‖w∗‖Id. Previous works of Nesterov & Polyak (2006); Carmon & Duchi (2019) show that the minimizer w∗ has a characterization; it satisfies ρ‖w∗‖ ≥ γ and ∇f(w∗) = A∗w∗ + b = 0. Furthermore, the minimizer w∗ is unique if ρ‖w∗‖ > γ. In this paper, we assume that the problem has a unique minimizer so that ρ‖w∗‖ > γ. The gradient of the cubic-regularized problem (2) is
∇f(w) = Aw + b+ ρ‖w‖w = A∗(w − w∗)− ρ ( ‖w∗‖ − ‖w‖ ) w. (12)
By applying the Heavy Ball algorithm (Algorithm 1) to the cubic-regularized problem (2), we see that it generates the iterates via wt+1 = wt − η∇f(wt) + β(wt − wt−1) = (Id − ηA− ρη‖wt‖Id)wt − ηb+ β(wt − wt−1).
(13)

heading: 4.2 ENTERING A BENIGN REGION FASTER
text: For the cubic-regularized problem, we define a different notion of a benign region from that for the phase retrieval. The benign region here is smooth, contains the unique global optimal point w∗, and satisfies a notion of one-point strong convexity (Kleinberg et al., 2018; Li & Yuan, 2017; Safran et al., 2020), w ∈ Rd : 〈w − w∗,∇f(w)〉 ≥ ϑ‖w − w∗‖2, where ϑ > 0. (14) We note that the standard strong convexity used in the definition of a benign region for phase retrieval could imply the one-point strong convexity here, but not vice versa (Hinder et al., 2020).
Previous work of Carmon & Duchi (2019) shows that if the norm of the iterate is sufficiently large, i.e. ρ‖w‖ ≥ γ − δ for any sufficiently small δ > 0, the iterate is in the benign region that contains the global minimizer w∗. To see this, by using the gradient expression of∇f(w), we have that
〈w − w∗,∇f(w)〉 =(w − w∗)> ( A∗ + ρ 2 (‖w‖ − ‖w∗‖)Id ) (w − w∗)
+ ρ2 ( ‖w∗‖ − ‖w‖ )2 (‖w‖+ ‖w∗‖).
(15)
The first term on the r.h.s. of the equality becomes nonnegative if the matrixA∗+ ρ2 (‖w‖−‖w∗‖)Id becomes PSD. Since A∗ (−γ+ ρ‖w∗‖)Id, it means that if ρ‖w‖ ≥ γ− ( ρ‖w∗‖− γ ) , the matrix becomes PSD (note that by the characterization, ρ‖w∗‖− γ > 0). Furthermore, if the size of iterate ρ‖w‖ satisfies ρ‖w‖ > γ − ( ρ‖w∗‖ − γ ) , the matrix becomes positive definite and consequently, we have that A∗ + ρ2 (‖w‖ − ‖w∗‖)Id ϑId for a number ϑ > 0. Therefore, (15) becomes 〈w − w∗,∇f(w)〉 ≥ ϑ‖w − w∗‖2, ϑ > 0. (16)
Therefore, the benign region of the cubic regularized problem can be characterized as B := {w ∈ Rd : ρ‖w‖ > γ − ( ρ‖w∗‖ − γ ) }. (17)
What we are going to show is that HB with a larger value of the momentum parameter β enters the benign region B faster. We have the following theorem, which shows that the size of the iterate ρ‖wt‖ will grow very fast to exceed any level below γ. Furthermore, the larger the momentum parameter β, the faster the growth, which shows the advantage of Heavy Ball over vanilla gradient descent.
Theorem 2. Fix any number δ that satisfies δ > 0. Define Tδ := min{t : ρ‖wt+1‖ ≥ γ − δ}. Suppose that the initialization satisfies w(1)0 b
(1) ≤ 0. Set the momentum parameter β ∈ [0, 1]. If the step size η satisfies η ≤ 1‖A‖2+ρ‖w∗‖ , then Heavy Ball (Algorithm 1 & 2) takes at most
Tδ ≤ 2 ηδ(1 + β/(1 + ηδ) ) log ( 1 + γ2+(1 + β/(1 + ηδ)) 4ρ|b(1)| )
number of iterations required to enter the benign region B.
Note that the case of β = 0 in Theorem 2 reduces to the result of Carmon & Duchi (2019)which analyzes vanilla gradient descent. The lemma implies that the higher the momentum parameter β, the faster that the iterate enters the benign region for which a linear convergence is possible; see also Figure 2 for the empirical results. Specifically, β reduces the number of iterations Tδ by a factor of (1 +β/(1 + ηδ)) (ignoring the β in the log factor as its effect is small), which also implies that for a smaller step size η, the acceleration effect due to the use of momentum is more evident. The factor can be approximately 1 + β for a small η. Lastly, the condition w(1)0 b
(1) ≤ 0 in Theorem 2 can be satisfied by w0 = w−1 = −r b‖b‖ for any r > 0.
Proof. (sketch; detailed proof is available in Appendix E) The theorem holds trivially when γ ≤ 0, so let us assume γ > 0. Recall the notation that w(1) represents the projection of w on the eigenvector v1 of the least eigenvalue λ(1)(A), i.e. w(1) = 〈w, v1〉. From the update, we have that
w (1) t+1
−ηb(1) = (Id + ηγ − ρη‖wt‖) w
(1) t −ηb(1) + 1 + β(w (1) t −w (1) t−1) −ηb(1) . (18)
Denote at := w
(1) t
−ηb(1) and in the detailed proof we will show that at ≥ 0 for all t ≤ Tδ . We can rewrite (18) as at+1 = (1 + ηγ− ρη‖wt‖)at + 1 +β(at− at−1) ≥ (1 + ηδ)at + 1 +β(at− at−1), where the inequality is due to that ρ‖wt‖ ≤ γ − δ for t ≤ Tδ . So we can now see that the dynamics is essentially in the form of (6) except that there is an additional 1 on the r.h.s of the inequality. Therefore, we can invoke Lemma 1 to show that the higher the momentum, the faster the iterate enters the benign region. In Appendix E, we consider the presence of 1 on the r.h.s and obtain a tighter bound than what Lemma 1 can provide.
We have shown that the simple dynamic results in entering a benign region that is one-point strongly convex to w∗ faster. However, different from the case of phase retrieval, we are not aware of any prior results of showing that the iterate generated by Heavy Ball keeps staying in a region that has the property (16) once the iterate is in the region. Carmon & Duchi (2019) show that for the cubic regularized problem, the iterate generated by vanilla gradient descent stays in the region under certain conditions, which leads to a linear convergence rate after it enters the benign region. Showing that the property holds for HB is not in the scope of this paper, but we empirically observe that Heavy Ball stays in the region. Subfigure (a) on Figure 2 shows that the norm ‖wt‖ is monotone increasing for a wide range of β, which means that the iterate stays in the benign region according to (17). Assuming that the iterate stays in the region, in Appendix F, we show a locally linear convergence of HB for which up to a certain threshold of β, the larger the β the better the convergence rate.

heading: 5 DISCUSSION AND CONCLUSION
text: Let us conclude by a discussion about the applicability of the simple dynamics to other non-convex optimization problems. Let A be a positive semi-definite matrix. Consider applying HB to top eigenvector computations, i.e. solving minw∈Rd:‖w‖≤1− 12w
>Aw, which is a non-convex optimization problem as it is about maximizing a convex function. The update of HB for this objective is wt+1 = (Id + ηA)wt + β(wt − wt−1). By projecting the iterate wt+1 on an eigenvector ui of the matrix A, we have that
〈wt+1, ui〉 = (1 + ηλi)〈wt, ui〉+ β(〈wt, ui〉 − 〈wt−1, ui〉). (19)
We see that this is in the form of the simple dynamics in Lemma 1 again. So one might be able to show that the larger the momentum parameter β, the faster the top eigenvector computation. This
‖w̃‖ w̃ and b = −(A + ρ‖w∗‖Id)w∗. The procedure makes
w∗ the global minimizer of problem instance (A, b, ρ). Patterns shown on this figure exhibit for other random problem instances as well.
connection might be used to show that the dynamics of HB momentum implicitly helps fast saddle points escape. In Appendix G, we provide further discussion and show some empirical evidence. We conjecture that if a non-convex optimization problem has an underlying structure like the ones in this paper, then HB might be able to exploit the structure and hence makes progress faster than vanilla gradient descent.

heading: A RELATED WORKS

heading: A.1 HEAVY BALL
text: We first note that Algorithm 1 and Algorithm 2 generate the same sequence of the iterates {wt} given the same initialization w0, the same step size η, and the same momentum parameter β.
Lessard et al. (2016) analyze the Heavy Ball algorithm for strongly convex quadratic functions by using tools from dynamical systems and prove its accelerated linear rate. Ghadimi et al. (2015) also show an O(1/T ) ergodic convergence rate for general smooth convex problems, while Sun et al. (2019) show the last iterate convergence on some classes of convex problems. Nevertheless, the convergence rate of both results are not better than gradient descent. Maddison et al. (2018) and Diakonikolas & Jordan (2019) study a class of momentum methods which includes Heavy Ball by a continuous time analysis. Can et al. (2019) prove an accelerated linear convergence to a stationary distribution for strongly convex quadratic functions under Wasserstein distance. Gitman et al. (2019) analyze the stationary distribution of the iterate of a class of momentum methods that includes SGD with momentum for a quadratic function with noise, as well as studying the condition of its asymptotic convergence. Loizou & Richtárik (2017) show linear convergence results of the Heavy Ball method for a broad class of least-squares problems. Loizou & Richtárik (2018) study solving a average consensus problem by HB, which could be viewed as a strongly convex quadratic function. Sebbouh et al. (2020) show a convergence result of stochastic HB under a smooth convex setting and show that it can outperform SGD under the assumption that the data is interpolated. Chen & Kolar (2020) study stochastic HB under a growth condition. Yang et al. (2018a) show an O(1/ √ T ) rate of convergence in expected gradient norm for smooth non-convex problems, but the rate is not better than SGD. Liu et al. (2020) provide an improved analysis of SGD with momentum. They show that SGD with momentum can converge as fast as SGD for smooth nonconvex settings in terms of the expected gradient norm. Krichene et al. (2020) show that in the continuous time regime, i.e. infinitesimal step size is used, stochastic HB converges to a stationary solution asymptotically for training a one-hidden-layer network with infinite number of neurons, but the result does not show a clear advantage compared to standard SGD. Lastly, Wang et al. (2020) show that the Heavy Ball’s momentum can help to escape saddle points faster and find a second order stationary point faster for smooth non-convex optimization. However, while their work focused on the stochastic setting, their main result required two assumptions on the statistical properties of the sequence of observed gradients; it is not clear whether these would hold in general. Specifically, they make an assumption called APAG (Almost Positively Aligned with Gradient), i.e. Et[〈∇f(wt),mt − gt〉] ≥ − 12‖∇f(wt)‖
2,where∇f(wt) is the deterministic gradient, gt is the stochastic gradient, and mt is the stochastic momentum. They also make an assumption called APCG (Almost Positively Correlated with Gradient), i.e. Et[〈∇f(wt),Mtmt〉] ≥ −c′ησmax(Mt)‖∇f(wt)‖2, where Mt is a PSD matrix that is related to a local optimization landscape.
We also note that there are negative results regarding Heavy Ball (see e.g. Lessard et al. (2016); Ghadimi et al. (2015); Kidambi et al. (2018)).

heading: A.2 MATRIX SENSING
text: Problem (1) can also be viewed as a special case of matrix factorization or matrix sensing. To see this, one can rewrite (1) as minU∈Rd×1 14n ∑n i=1 ( yi−〈Ai, UU>〉 )2 , whereAi = xix>i ∈ Rd×d and the dot product represents the matrix trace. Li et al. (2018) show that when the matrices {Ai} satisfy restricted isometry property (RIP) and U ∈ Rd×d, gradient descent can converge to a global solution with a close-to-zero random initialization. Yet, if the matrix Ai is in the form of a rank-one matrix product, the matrix might not satisfy RIP and a modification of the algorithm might be required (Li et al. (2018)). Li et al. (2019), a different group of authors, show that with a carefully-designed initialization (e.g. spectral initialization), gradient descent will be in a benign region in the beginning and will converge to a global optimal point. In our work, we do not assume that xix>i satisfies RIP neither do we assume a carefully-designed initialization like spectral initialization is available. Li & Lin (2020) show a local convergence to an optimal solution by Nesterov’s momentum for a matrix factorization problem; the initial point needs to be in the neighborhood of an optimal solution. In contrast, we study Polyak’s momentum and are able to establish global convergence to an optimal solution with a simple random initialization. Gunasekar et al. (2017); Gidel et al. (2019); You et al.
(2020) study implicit regularization of gradient descent for the matrix sensing/matrix factorization problem. The directions are different from ours.

heading: B POPULATION GRADIENT
text: Lemma 2. Assume that ‖w∗‖ = 1. Ex∼N(0,Id) [ (x>w)3x− (x>w∗)2(x>w)x ] = ( 3‖w‖2 − 1 ) w − 2(w>∗ w)w∗.
Proof. In the following, denote
Q := E[(x>w)3x] R := E[(x>w∗)2(x>w)x].
Now define h(w) := E[(x>w)4]. We have that h(w) = 3‖w‖4 as x>w ∼ N(0, ‖w‖2) (i.e. the fourth order moment). Then, ∇h(w) = 4E[(x>w)3x] = 4Q. So Q = 14∇h(w) = 3‖w‖ 2w.
For the other term, define g(w) := E[(x>w∗)2(x>w)2]. Given that[ x>w∗ x>w ] ∼ N ( [0 0 ] , [ ‖w∗‖2 w>∗ w w>∗ w ‖w‖2 ] ) , (20)
we can write x>w∗ d∼ θ1x>w + θ2z, where z ∼ N(0, 1), θ1 := w > ∗ w ‖w‖2 , and θ 2 2 := ‖w∗‖2 − (w>∗ w) 2 ‖w‖2 . Then,
g(w) = θ21E[(x>w)4] + 2θ1θ2E[(x>w)3z] + θ22E[z2(x>w)2] = 3(w>∗ w)
2 + θ22E[z2]E[(x>w)2] = 2(w>∗ w)
2 + ‖w‖2‖w∗‖2. So we have that ∇g(w) = 2E[(x>w∗)2(x>w)x] = 2R, which in turn implies that
R = 1 2 ∇g(w) = 1 2
( 4(w>∗ w)w∗ + 2‖w∗‖w ) = 2(w>∗ w)w∗ + ‖w∗‖2w.
Combining the above results, we have that ∇F (w) = Q−R = 3‖w‖2w − ‖w∗‖2w − 2(w>∗ w)w∗ = ( 3‖w‖2 − 1 ) w − 2(w>∗ w)w∗. (21)

heading: C SIMPLE LEMMAS
text: Lemma 1: For a positive number θ and the momentum parameter β ∈ [0, 1], if a non-negative sequence {at} satisfies a0 ≥ a−1 > 0 and that for all t ≤ T ,
at+1 ≥ (1 + θ)at + β(at − at−1), (22) then {at} satisfies
at+1 ≥ ( 1 + ( 1 + β1+θ ) θ ) at, (23)
for every t = 1, . . . , T + 1. Similarly, if a non-positive sequence {at} satisfies a0 ≤ a−1 < 0 and that for all t ≤ T , at+1 ≤ (1 + θ)at + β(at − at−1), then {at} satisfies
at+1 ≤ ( 1 + ( 1 + β1+θ ) θ ) at, (24)
for every t = 1, . . . , T + 1.
Proof. Let us first prove the first part of the statement. In the following, we denote c := 11+θ . For the base case t = 1, we have that
a2 ≥ (1 + θ)a1 + β(a1 − a0) ≥ (1 + θ)a1 + βcθa1, (25)
where the last inequality holds because a1 − a0 ≥ θca1 ⇐⇒ a1 ≥ 11−θca0, as a1 ≥ (1 + θ)a0 = 1 1−θca0. Now suppose that it holds at iteration t, at+1 ≥ (1+(1+βc)θ)at. Consider iteration t+1, we have that
at+2 ≥ (1 + θ)at+1 + β(at+1 − at) ≥ (1 + θ)at+1 + θcβat+1, (26)
where the last inequality holds because at+1 − at ≥ θcat+1 ⇐⇒ at+1 ≥ 11−θcat as at+1 ≥ (1 + (1 + βc)θ)at ≥ 11−θcat, given the assumption at t and that c := 1 1+θ .
The second part of the statement can be proved similarly.
Lemma 3. For a positive number θ < 1 and the momentum parameter β ∈ [0, 1] that satisfy (1 + β1−θ )θ < 1, if a non-negative sequence {bt} satisfies b0 ≤ b−1 and that for all t ≤ T ,
bt+1 ≤ (1− θ)bt + β(bt − bt−1), (27)
then {bt} satisfies
bt+1 ≤ ( 1− ( 1 + β1−θ ) θ ) bt, (28)
for every t = 1, . . . , T + 1. Similarly, if a non-positive sequence {bt} satisfies b0 ≥ b−1 and that for all t ≤ T , bt+1 ≥ (1− θ)bt + β(bt − bt−1), then {bt} satisfies
bt+1 ≥ ( 1− ( 1 + β1−θ ) θ ) bt, (29)
for every t = 1, . . . , T + 1.
Proof. Let us first prove the first part of the statement. In the following, we denote c := 11−θ . For the base case t = 1, we have that
b2 ≤ (1− θ)b1 + β(b1 − b0) ≤ (1− θ)b1 − βcθb1, (30)
where the last inequality holds because
b1 − b0 ≤ −θcb1 ⇐⇒ b1 ≤ 1
1 + θc b0, (31)
as b1 ≤ (1− θ)b0 = 11+θcb0 due to that c = 1 1−θ . Now suppose that it holds at iteration t.
bt+1 ≤ (1− (1 + βc)θ)bt.
Consider iteration t+ 1, we have that
bt+2 ≤ (1− θ)bt+1 + β(bt+1 − bt) ≤ (1− θ)bt+1 − θcβbt+1, (32)
where the last inequality holds because
bt+1 − bt ≤ −θcbt+1 ⇐⇒ bt+1 ≤ 1
1 + θc bt (33)
as bt+1 ≤ (1− (1 + βc)θ)bt ≤ 11+θcbt, due to the induction at t and that c = 1 1−θ .
The second part of the statement can be proved similarly.

heading: D PROOF OF THEOREM 1
text: Recall the recursive system (11).
w ‖ t+1 = ( 1 + 3η(1− ‖wt‖2) + ηξt ) w ‖ t + β(w ‖ t − w ‖ t−1)
w⊥t+1[j] = ( 1 + η(1− 3‖wt‖2) + ηρt,j ) w⊥t [j] + β(w ⊥ t [j]− w⊥t−1[j]),
(34)
where {ξt} and {ρt,j} are the perturbation terms. In the analysis, we will show that the sign of w‖t never changes during the execution of the algorithm. Our analysis divides the iterations of the first stage to several sub-stages. We assume that the size of the initial point w0 is small so that it begins in Stage 1.1.
• (Stage 1.1) considers the duration when ‖wt‖ ≤ √ 4 9 + cn, which lasts for at most T0
iterations, where T0 is defined in Lemma 4. A by-product of our analysis shows that |w‖T0+1| ≥ √ 4 9 + cn.
• (Stage 1.2) considers the duration when the perpendicular component ‖w⊥t ‖ is decreasing and eventually falls below ζ/2, which consists of all iterations T0 ≤ t ≤ T0 + Tb, where Tb is defined in Lemma 5.
• (Stage 1.3) considers the duration when |w‖t | is converging to the interval [1− ζ 2 , 1 + ζ 2 ], if
it was outside the interval, which consists of all iterations T0 + Tb ≤ t ≤ T0 + Tb + Ta, where Ta is defined in Lemma 6.
In stage 1.1, both the signal component |w‖t | and ‖w⊥t ‖ grow in the beginning (see also Figure 1). The signal component grows exponentially and consequently it only takes a logarithm number of
iterations T0 + 1 to reach |w‖T0+1| > √ 4 9 + cn, which also means that ‖wT0+1‖ > √ 4 9 + cn. Moreover, the larger the momentum parameter β, the smaller the number of T0. After ‖wt‖ passing the threshold, the iterate enters Stage 1.2. Lemma 4. (Stage 1.1) Denote ca := 11+ η2 . There will be at most T0 := log( √ 4 9 +cn/|w ‖ 0 |) log (
1+(1+caβ)η 5 3 ) iterations such that w‖t ≤ √ 4 9 + cn. Furthermore, for all t in this stage, ‖wt‖ ≤ 5 6 − cn 3 .
In stage 1.2, we have that ‖wt‖2 ≥ |w‖t |2 ≥ 49 + cn > 1 3 so that the perpendicular component ‖w ⊥ t ‖ is decaying while the signal component |w‖t | keeps growing before reaching 1 (see also Figure 1). In particular, the perpendicular component decays exponentially so that at most additional Tb iterations is needed to fall below ζ2 (i.e. ‖w ⊥ t ‖ ≤ ζ 2 ), with larger β leading to a smaller number of Tb. Notice that ifw‖t satisfies ||w ‖ t |−1| ≤ ζ 2 when ‖w ⊥ t ‖ falls below ζ 2 , we immediately have that Tζ ≤ T0+Tb; otherwise, the iterate enters the next stage, Stage 1.3.
Lemma 5. (Stage 1.2) Denote cb := 11− η3 and ω̄ := maxt ‖w ⊥ t ‖ ≤ √ 4 9 + cn. There will be at most Tb ≤ log( ζ2ω̄ )
log(1− η3 (1+cbβ)) iterations such that
√ 4 9 + cn ≤ |w ‖ t | and ‖w⊥t ‖ ≥ ζ 2 .
In stage 1.3, we show that |w‖t | converges towards 1 linearly, given that ‖w⊥t ‖ ≤ ζ 2 . Specifically, after at most additional Ta iterations, we have that ||w‖t | − 1| ≤ ζ 2 and that larger β reduces the number of iterations Ta.
Lemma 6. (Stage 1.3) Denote cζ,g := 11+η ζ2 , cζ,d := 11−ηζ , and ω := maxt |w
‖ t | ≤ √ 10 9 + 1 3cn.
There will be at most Ta := max ( log((1−ζ/2)/√ 49 +cn)
log(1+η ζ2 (1+cζ,gβ)) ,
log( 1+ζ/2 ω ) log(1−ηζ(1+cζ,dβ)) ) iterations such that
||w‖t | − 1| ≥ ζ 2 and ‖w ⊥ t ‖ ≤ ζ 2 .
By combining the result of Lemma 4, Lemma 5 and Lemma 6, we have that
T ≤ T0 + Tb + Ta = log(
√ 4 9 +cn
|w‖0 | ) log ( 1 + η 53 (1 + caβ) ) + log( ζ2ω̄ ) log(1− η3 (1 + cbβ))
+ max ( log( 1−ζ/2√ 4 9 +cn )
log(1 + η ζ2 (1 + cζ,gβ)) ,
log( 1+ζ/2ω ) log(1− ηζ(1 + cζ,dβ)) )
. 6 log(
√ 4 9 +cn
|w‖0 | )
5η(1 + caβ) +
6 log( 2ω̄ζ )
η(1 + cbβ) + max
( 2 log( 1−ζ/2√ 4 9 +cn )
ηζ(1 + cζ,gβ)) ,
2 log( ω1+ζ/2 )
ηζ(1 + cζ,dβ)) ) . log d
η (1 + cηβ) ,
(35)
where the last inequality uses that |w‖0 | & 1√d log d due to the random isotropic initialization. The detailed proof of Lemma 4-6 is available in the following subsections.
After t ≥ Tζ , the iterate enters a benign region that is locally strong convex, smooth, twice differentiable, and contains w∗ (or −w∗) (Ma et al., 2017; White et al., 2016), which allows us to use the existing result of gradient descent with Heavy Ball momentum for showing its linear convergence. In particular, the result of local landscape (e.g. Ma et al. (2017); White et al. (2016)) and the known convergence result of gradient descent with Heavy Ball momentum (e.g. Saunders (2018); Polyak (1964); Lessard et al. (2016); Xiong et al. (2020)) can be used to show that for all t > Tζ , dist(wt, w∗) ≤ (1− ν)t−Tζdist(wTζ , w∗) ≤ (1− ν)t−Tζζ, for some number 1 > ν > 0.
D.1 STAGE 1.1 Lemma: 4 (Stage 1.1) Denote ca := 11+ η2 . There will be at most T0 := d log( √ 4 9 +cn/|w ‖ 0 |) log (
1+(1+caβ)η 5 3 )e iterations such that ‖wt‖ ≤ √ 4 9 + cn. Furthermore, we have that |w ‖ T0+1 | > √ 4 9 + cn.
Proof. Let us first assume that w‖t > 0 and denote at := w ‖ t . By using that ‖wt‖ ≤ √ 4 9 + cn in this stage, we can lower-bound the growth rate of at as
at+1 ≥ ( 1 + 3η(1− ‖wt‖2)− η|ξt| ) at + β(at − at−1)
≥ ( 1 + 3η(1− 4
9 + cn)− cn
) at + β(at − at−1)
≥ ( 1 + η 5
3
) at + β(at − at−1)
≥ ( 1 + (1 + caβ)η 5
3
) at,
(36)
where in the last inequality we use Lemma 1 and that
ca = 1
1 + η 53 . (37)
Observe that ( 1 + (1 + caβ)η 5 3 ) > 1. Consequently, the sign of w‖t never change in this stage. So
for w‖t ≥ √ 4 9 + cn it takes number of iterations at most
T0 := d log(
√ 4 9 +cn
|w‖0 | ) log ( 1 + (1 + caβ)η 5 3
)e ≤ 2 log( √ 4 9 +cn |w‖0 | )
(1 + caβ)η 5 3
. (38)
Similar, when w‖t < 0, we can show that after at most T0 iterations, w ‖ t falls below − √ 4 9 + cn.
Since ‖wt‖ > |w‖t |, it means that there will be at most T0 iterations such that ‖wt‖ ≤ √ 4 9 + cn.
D.2 STAGE 1.2
Lemma 5: (Stage 1.2) Denote cb := 11− η3 and ω̄ := maxt ‖w ⊥ t ‖ ≤ √ 4 9 + cn. There will be at most Tb := d log( ζ2ω̄ )
log(1− η3 (1+cbβ)) e iterations such that √ 4 9 + cn ≤ |w ‖ t | and ‖w⊥t ‖ ≥ ζ 2 .
Proof. Let t′ be the last iteration of the previous stage. We have that ‖wt′‖ ≥ ‖w‖t′‖ ≥ √ 4 9 + cn. Denote at := |w‖t |. In this stage, we have that at keeps increasing until ‖wt‖2 & 1. Moreover, w ‖ t remains the same sign as the previous stage. Now fix an element j 6= 1 and denote bt := |w⊥t [j]|. From Lemma 7, we know that the magnitude bt := |w⊥t [j]| is non-increasing in this stage. Furthermore, we can show the decay of bt as follows. If wt[j], wt−1[j] > 0,
wt+1[j] ≤ ( 1 + η(1− 3‖wt‖2) + η|ρt,j | ) wt[j] + β(wt[j]− wt−1[j])
≤ ( 1 + η ( 1− 3(4
9 + cn)
) + ηcn ) wt[j] + β(wt[j]− wt−1[j])
≤ ( 1− η
3
) wt[j] + β(wt[j]− wt−1[j])
≤ ( 1− η
3 (1 + cbβ)
) wt[j],
(39)
where in the last inequality we used Lemma 3, as the condition (1 + β1− η3 ) η 3 < 1 is satisfied, and we denote that
cb := 1
1− η/3 . (40)
On the other hand, if wt[j], wt−1[j] < 0, wt+1[j] ≥ ( 1 + η(1− 3‖wt‖2) + η|ρt,j | ) wt[j] + β(wt[j]− wt−1[j])
≥ ( 1 + η ( 1− 3(4
9 + cn)
) + ηcn ) wt[j] + β(wt[j]− wt−1[j])
≥ ( 1− η
3
) wt[j] + β(wt[j]− wt−1[j])
≥ ( 1− η
3 (1 + cbβ)
) wt[j],
(41)
where in the last inequality we used Lemma 3, as the condition (1 + β1− η3 ) η 3 < 1 is satisfied, and we denote that
cb := 1
1− η/3 . (42)
The inequalities of (39) and (41) allow us to write |wt+1[j]| ≤ ( 1− η3 (1 + cbβ) ) |wt[j]|. Taking the square of both sides and summing all dimension j 6= 1, we have that
‖w⊥t+1‖2 ≤ ( 1− η
3 (1 + cbβ) )2‖w⊥t ‖2. (43) Consequently, for ‖w⊥t ‖ to fall below ζ/2, it takes at most
Tb := d log( ζ/2‖w⊥ t′‖ )
log(1− η3 (1 + cbβ)) e ≤ d
log( ζ2ω̄ )
log(1− η3 (1 + cbβ)) e (44)
iterations.
Lastly, Lemma 7 implies that at the time that the magnitude of w⊥t [j] starts decreasing, ‖wt‖2 ≤ 4 9 + cn, which in turn implies that ω̄ := maxt ‖w ⊥ t ‖ ≤ √ 4 9 + cn.
D.3 STAGE 1.3
Lemma 6: (Stage 1.3) Denote cζ,g := 11+η ζ2 , cζ,d := 11−ηζ , and ω := maxt |w
‖ t | ≤ √ 10 9 + 1 3cn.
There will be at most Ta := max ( log((1−ζ/2)/√ 49 +cn)
log(1+η ζ2 (1+cζ,gβ)) ,
log( 1+ζ/2 ω ) log(1−ηζ(1+cζ,dβ)) ) iterations such that
||w‖t | − 1| ≥ ζ 2 and ‖w ⊥ t ‖ ≤ ζ 2 .
Proof. Denote t′ the last iteration of the previous stage. We have that t′ ≤ T0 + Tb. Since w‖t does not change the sign in stage 1.1 and 1.2, w.l.o.g, we assume that w‖t > 0. Denote at := w ‖ t . We consider at′ in two cases: at′ ≤ 1− ζ2 and at′ ≥ 1 + ζ 2 . If at′ ≤ 1− ζ 2 , then we have that for all t in this stage, ‖wt‖2 ≤ a2t + ‖w⊥t ‖2 ≤ (1− ζ 2 ) 2 + ( ζ2 ) 2 ≤ 1− ζ2 , for any sufficiently small ζ. So at grows as follows.
at+1 ≥ ( 1 + 3η(1− ‖wt‖2)− η|ξt| ) at + β(at − at−1)
≥ ( 1 + 3η ζ
2 − ηcn
) at + β(at − at−1)
(a) ≥ ( 1 + η ζ
2
) at + β(at − at−1)
(b) ≥ (1 + η ζ 2 (1 + cζ,gβ))at,
(45)
where (a) is by ζ ≥ cn and (b) is due to that Lemma 1 and that
cζ,g = 1
1 + η ζ2 . (46)
Consequently, it takes at most
d log 1−ζ/2√ 4 9 +cn
log(1 + η ζ2 (1 + cζ,gβ)) e (47)
number of iterations in this stage for at to rise above 1− ζ2 . On the other hand, if at′ ≥ 1 + ζ 2 , then we can lower bound ‖wt‖2 in this stage as ‖wt‖2 ≥ a2t ≥ (1 + ζ 2 ) 2. We have that
at+1 ≤ ( 1 + 3η(1− ‖wt‖2) + η|ξt| ) at + β(at − at−1)
≤ ( 1− 3ηζ + cn ) at + β(at − at−1)
(a) ≤ ( 1− ηζ ) at + β(at − at−1)
(b) ≤ (1− ηζ(1 + cζ,dβ))at,
(48)
where (a) uses ζ ≥ cn and (b) uses Lemma 3
cζ,d := 1
1− ηζ . (49)
That is, at is decreasing towards 1 + ζ/2. Denote ω := maxt |w‖t |. we see that it takes at most
d log 1+ζ/2ω
log(1− ηζ(1 + cζ,dβ)) e (50)
number of iterations in this stage for at to fall below 1 + ζ/2. Lastly, Lemma 8 implies that at the time that the magnitude of w‖t starts decreasing, ‖wt‖2 ≤ 109 + 1 3cn, which in turn implies that
ω := maxt |w‖t | ≤ √ 10 9 + 1 3cn.
On the other hand, by Lemma 10, the magnitude of the perpendicular component is non-increasing in this stage, and hence ‖w⊥t ‖ keeps staying below ζ/2.
Similar analysis holds for w‖t < 0, hence we omitted the details.

heading: D.4 SOME SUPPORTING LEMMAS
text: Lemma 7. Suppose that η satisfies η ≤ 1 36 (
1 3 +cn ) max{cm,1} . Let t0 be the first time such that
‖wt0‖2 ≥ 13 + cn. Then, there exists a time τ ≤ t0 such that w ⊥ τ+1[j] ≤ w⊥τ [j], if w⊥τ [j] ≥ 0; Similarly, w⊥τ+1[j] ≥ w⊥τ [j], if w⊥τ [j] ≤ 0. Furthermore, we have that ‖wt0‖2 ≤ 49 + cn.
Proof. Recall that w⊥t+1[j] = ( 1 + η(1 − 3‖wt‖2) + ηρt,j ) w⊥t [j] + β(w ⊥ t [j] − w⊥t−1[j]). W.l.o.g, let us consider w⊥0 [j] > 0. Assume that w ⊥ t [j] ≥ w⊥t−1[j] for all t ≤ t0; otherwise, there exists a time τ < t0 such that w⊥τ [j] ≤ w⊥τ−1[j].
Denote λt0,j := η(1− 3‖wt0‖2) + ηρt0,j . Since ‖wt0‖2 ≥ 13 + cn, we have that λt0,j ≤ 0. We can rewrite the dynamics as [
w⊥t0+1[j] w⊥t0 [j]
] = [ 1 + λt0,j + β −β
1 0 ] [ w⊥t0 [j] w⊥t0−1[j] ] (51)
We have that
‖ [ w⊥t0+1[j] w⊥t0 [j] ] ‖ ≤ ‖ [ 1 + λt0,j + β −β 1 0 ] ‖2 · ‖ [ w⊥t0 [j] w⊥t0−1[j] ] ‖. (52)
We will show that w⊥t0+1[j] ≤ w ⊥ t0−1[j] if w ⊥ t0−1[j] > 0, which means that the magnitude of w⊥[j] has stopped increasing. It suffices to show that the spectral norm of the matrix
‖ [ 1 + λt0,j + β −β
1 0
] ‖2 is not greater than 1.
Note that the roots of the characteristic equation of the matrix, z2 − (1 + λt0,j + β)z + β are( 1+λt0,j+β ) ± √ (1+λt0,j+β) 2−4β
2 .
If the roots are complex conjugate, then the magnitude of the roots is at most √ β; consequently the spectral norm is at most √ β ≤ 1. On the other hand, if the roots are real, to show that the larger
root is not larger than 1, it suffices to show that √
(1 + λt0,j + β) 2 − 4β ≤ 1− λt0,j − β, which is
guaranteed if λt0,j ≤ 0. To show that the smaller root is not greater than −1, we need to show that√ (1 + λt0,j + β) 2 − 4β ≤ 3 + λt0,j + β, which is guaranteed if λt0,j ≥ −1.
By definition, ‖wt0−1‖2 < 13 + cn. If η ≤ min{ 1 36‖wt0−1‖cm , 118cm }, then by invoking Lemma 9, we have that ‖wt0‖2 − ‖wt0−1‖2 ≤ 19 and consequently, we have that ‖wt0‖ 2 ≤ 49 + cn. Thus, by choosing η satisfies η ≤ 136(1/3+cn)cm , we have that ‖wt0‖ 2 ≤ 49 + cn. Using the upper-bound of ‖wt0‖2 and the constraint of η, we have that λt0,j ≥ −η( 13 + cn) ≥ −1.
Similar analysis when w⊥t [j] is negative; and hence omitted.
Lemma 8. Suppose that η satisfies η ≤ 1 36 (
1+ 13 cn ) max{cm,1} . Let t1 be the first time (if exist) such
that ‖wt1‖2 ≥ 1+ 13cn. Then, there exists a time τ ≤ t1 such that w ‖ τ+1 ≤ w ‖ τ , if w ‖ τ ≥ 0. Similarly, w ‖ τ+1 ≥ w ‖ τ , if w ‖ τ ≤ 0. Furthermore, we have that ‖wt1‖2 ≤ 109 + 1 3cn. Proof. Recall that w‖t+1 = ( 1 + 3η(1−‖wt‖2) + ηξt ) w ‖ t + β(w ‖ t −w ‖ t−1). W.l.o.g, let us consider w ‖ 0 > 0. Assume that w ‖ t ≥ w ‖ t−1 for all t ≤ t1; otherwise, there exists a time τ < t1 such that w ‖ τ ≤ w‖τ−1.
Denote λt1 := 3η(1 − ‖wt1‖2) + ηξt1 . Since ‖wt1‖2 ≥ 1 + 13cn, we have that λt1 ≤ 0. We can rewrite the dynamics as [
w ‖ t1+1
w ‖ t1
] = [ 1 + λt1 + β −β
1 0
][ w ‖ t1
w ‖ t1−1
] (53)
We have that
‖
[ w ‖ t1+1
w ‖ t1
] ‖ ≤ ‖ [ 1 + λt1 + β −β
1 0
] ‖2 · ‖ [ w ‖ t1
w ‖ t1−1
] ‖. (54)
The analysis essentially follows the same lines as Lemma 7. Specifically, to show that w‖t1+1 ≤ w ‖ t1−1, it suffices to ensure that λt1 := 3η(1− ‖wt1‖
2) + ηξt1 ∈ [−1, 0]. We have that λt1 ≤ 0 by the definition of t1. Furthermore, by the definition of t1 − 1, we have that ‖wt1−1‖2 < 1 + 13cn. So if η ≤ min{ 136‖wt1−1‖cm , 1 18cm }, then by invoking Lemma 9, we have that ‖wt1‖2−‖wt1−1‖2 ≤ 19 and consequently, we have that ‖wt1‖2 ≤ 109 + 1 3cn. Thus, by choosing η satisfies η ≤ 1 36 (
1+ 13 cn ) cm ,
we have that ‖wt1‖2 ≤ 109 + 1 3cn. Using this upper-bound of ‖wt1‖ 2 and the constraint of η, we have that λt1 ≥ −η( 13 + 2cn) ≥ −1. Therefore, we have completed the proof.
Lemma 9. Assume that the norm of the momentum is bounded for all t ≤ Tζ , i.e. ‖mt‖ ≤ cm,∀t ≤ Tζ . Set the step size η satisfies η ≤ min{ 136‖wt−1‖cm , 1 18cm }. Then, we have that
‖wt‖2 − ‖wt−1‖2 ≤ 1
9 .
Proof. To see this, we will use the alternative presentation Algorithm 2, which shows that wt = wt−1 − ηmt−1, where the momentum mt−1 stands for the weighted sum of gradients up to (and including) iteration t− 1, i.e. mt−1 = ∑t−1 s=0 β
t−1−s∇f(ws). Using the expression, we can expand ‖wt‖2 − ‖wt−1‖2 as
‖wt‖2 − ‖wt−1‖2 = ‖wt−1 − ηmt−1‖2 − ‖wt−1‖2 = −2η〈wt−1,mt−1〉+ η2‖mt−1‖2
≤ 2η‖wt−1‖‖mt−1‖+ η2‖mt−1‖2 ≤ 1
9 .
(55)
where the last inequality holds if η ≤ min{ 136‖wt−1‖cm , 1 18cm }.
Lemma 10. Fix an index j. Set η ≤ 1 36(1+ 13 cn) max{cm,1} and β ≤ 1. Suppose that ‖wt‖2 ≥ 13 +cn. If for a number R > 0, we have that ‖ [ w⊥t [j] w⊥t−1[j] ] ‖ ≤ R, then ‖ [ w⊥t+1[j] w⊥t [j] ] ‖ ≤ R.
Proof. The proof is similar to that of Lemma 7. We have that ‖ [ w⊥t+1[j] w⊥t [j] ] ‖ ≤ ‖ [ 1 + λt,j + β −β 1 0 ] ‖2 · ‖ [ w⊥t [j] w⊥t−1[j] ] ‖. (56)
Denote λt,j := η(1 − 3‖wt‖2) + ηρt,j . As the proof of Lemma 7, to show that the spectral norm of the matrix ‖ [ 1 + λt,j + β −β
1 0
] ‖2 is not greater than one. It suffices to have β ≤ 1 and
that λt,j ∈ [−1, 0]. By the assumption, it holds that ‖wt‖2 ≥ 13 + cn, so we have that λt,j ≤ 0. Furthermore, by Lemma 8, if the step size η satisfies η ≤ 1 36 (
1+ 13 cn ) max{cm,1} , we have that
‖wt‖2 ≤ 109 + 1 3cn. Therefore, using the upper-bound of the step size and the norm, we can obtain that λt,j := η(1− 3‖wt‖2) + ηρt,j ≥ −1. Hence, we have completed the proof.

heading: E PROOF OF THEOREM 2
text: To prove Theorem 2, we will need the following lemma. Lemma 11. Fix any number δ > 0. Define Tδ := min{t : ρ‖wt+1‖ ≥ γ − δ}. Assume that η ≤ 1‖A‖2+ρ‖w∗‖ . Suppose that w (1) 0 b (1) ≤ 0. Then, we have that w(1)t b(1) ≤ 0, for all 0 ≤ t ≤ Tδ .
Proof. The lemma holds trivially when γ ≤ 0, so let us assume γ > 0. Recall the Heavy Ball generates the iterates as
w (1) t+1 = (1− ηλ(1)(A)− ρη‖wt‖)w (1) t − ηb(1) + β(w (1) t − w (1) t−1).
We are going to show that for all t,
w (1) t b (1) ≤ 0 and (w(1)t − w (1) t−1)b (1) ≤ −cbww(1)t b(1), (57)
for any constant cbw ≥ 0. The initialization guarantees that w(1)0 b(1) ≤ 0 and that (w (1) 0 − w (1) −1)b (1) = 0 ≤ −cbww(1)0 b(1). Suppose that (57) is true at iteration t. Consider iteration t+ 1.
w (1) t+1b (1) = (1− ηλ(1)(A)− ρη‖wt‖)w(1)t b(1) − η(b(1))2 + β(w (1) t − w (1) t−1)b (1)
≤ (1− ηλ(1)(A)− ρη‖wt‖ − βcbw)w(1)t b(1) − η(b(1))2 ≤ 0, (58)
where the first inequality is by induction at iteration t and the second one is true if (1 − ηλ(1)A − ρη‖wt‖ − βcbw) ≥ 0, which gives a constraints about η,
1− ηλ(1)(A)− ρη‖wt‖ ≥ cbw. (59)
Now let us switch to show that (w(1)t+1 −w (1) t )b (1) ≤ −cbww(1)t+1b(1), which is equivalent to showing that w(1)t+1b (1) ≤ 11+cbww (1) t b (1). From (58), it suffices to show that
(1− ηλ(1)(A)− ρη‖wt‖ − βcbw)w(1)t b(1) − η(b(1))2 ≤ 1
1 + cbw w
(1) t b (1). (60)
Since w(1)t b (1) ≤ 0, a sufficient condition of the above inequality is
1− ηλ(1)(A)− ρη‖wt‖ − βcbw − 1
1 + cbw ≥ 0. (61)
Now using that ρ‖wt‖ ≤ γ − δ and that 11+x ≤ 1− 1 2x for x ∈ [0, 1]. It suffices to have that
1− ηλ(1)(A)− η(γ− δ)−βcbw− 1
1 + cbw ≥ 1 + ηδ−βcbw− 1 +
1 2 cbw ≥ ηδ− 1 2 cbw ≥ 0. (62)
By setting cbw = 0, we have that the inequality is satisfied. Substituting cbw = 0 to (59), we have that
η ≤ 1 λ(1)(A) + ρ‖wt‖ ,
Recall that ρ‖wt‖ ≤ γ − δ for all t ≤ Tδ . So we have that η ≤ 1 ‖A‖2+ ( γ−δ ) ·1γ−δ≥0 . Furthermore, by using that ρ‖w∗‖ > γ, it suffices to have that η ≤ 1‖A‖2+ρ‖w∗‖ . We have completed the proof.
Given Lemma 11, we are ready for proving Theorem 2.
Proof. (of Theorem 2) The lemma holds trivially when γ ≤ 0, so let us assume γ > 0. Recall the notation that w(1) represents the projection of w on the eigenvector v1 of the least eigenvalue λ(1)(A), i.e. w(1) = 〈w, v1〉. From the update rule, we have that
w (1) t+1
−ηb(1) = (Id + ηγ − ρη‖wt‖) w
(1) t −ηb(1) + 1 + β(w (1) t −w (1) t−1) −ηb(1) . (63)
Denote at := w
(1) t
−ηb(1) . We can rewrite (63) as
at+1 = (1 + ηγ − ρη‖wt‖)at + 1 + β(at − at−1) ≥ (1 + ηδ)at + 1 + β(at − at−1), (64)
where the inequality is due to that ρ‖wt‖ ≤ γ − δ for t ≤ Tδ . Now we are going to show that, at+1 ≥ (1+ηδ+ ηδ1+ηδβ)at+1. For the above inequality to hold, it suffices to show that at−at−1 ≥ ηδ 1+ηδat. That is, at ≥ 1
1−ηδ/(1+ηδ)at−1. The base case t = 1 holds because a1 ≥ (1 + ηδ)a0 + 1 ≥ 1
1−ηδ/(1+ηδ)a0. Suppose that at iteration t, we have that at ≥ 1 1−ηδ/(1+ηδ)at−1. Consider t + 1, we have that at+1 ≥ (1 + ηδ)at + 1 + β(at − at−1) ≥ (1 + ηδ)at + 1 ≥ 11−ηδ/(1+ηδ)at, where the second to last inequality is because at ≥ 11−ηδ/(1+ηδ)at−1 implies at ≥ at−1. Therefore, we
have completed the induction. So we have shown that w
(1) t+1
−ηb(1) = (Id + ηγ − ρη‖wt‖) w
(1) t
−ηb(1) + 1 +
β(w (1) t −w (1) t−1)
−ηb(1) ≥ (1 + ηδ + ηδ 1+ηδβ) w
(1) t
−ηb(1) + 1. Recursively expanding the inequality, we have that
w (1) t+1
−ηb(1) ≥(1 + ηδ + ηδ 1+ηδβ) w
(1) t
−ηb(1) + 1 ≥ (1 + ηδ + ηδ 1+ηδβ) 2 w
(1) t−1
−ηb(1) + (1 + ηδ + ηδ 1+ηδβ) + 1
≥ . . . ≥ 1 ηδ(1+ β1+ηδ )
( (1 + ηδ + ηδ1+ηδβ) t − 1 ) .
(65)
Therefore, γ−δρ (a) ≥ ‖wTδ‖ ≥ |w (1) Tδ | (b) ≥ |b (1)| δ+δβ/(1+ηδ) ( (1 + ηδ + ηδ1+ηδβ) Tδ − 1 ) where (a) uses that for t ≤ Tδ , ρ‖wt‖ ≤ γ − δ, and (b) uses (65) and Lemma 11 that w(1)t b(1) ≤ 0. Consequently,
Tδ ≤ log ( 1+ (γ−δ)(δ+δβ/(1+ηδ)) ρ|b(1)| ) log(1+ηδ+ ηδ1+ηδ β) (a) ≤ 2ηδ(1+β/(1+ηδ)) log ( 1 + (γ−δ)δ(1+β/(1+ηδ)) ρ|b(1)|
) (b)
≤ 2ηδ(1+β/(1+ηδ)) log ( 1 + γ2+(1+β/(1+ηδ)) 4ρ|b(1)| ) ,
(66)
where (a) uses that log(1 + x) ≥ x2 for any x ∈ [0,∼ 2.51], and (b) uses γδ − δ 2 ≤ γ
2 + 4 .

heading: F CONVERGENCE OF HB FOR THE CUBIC-REGULARIZED PROBLEM
text: Theorem 3. Assume that the iterate stays in the benign region that exhibits one-point strong convexity to w∗, i.e. B := {w ∈ Rd : ρ‖w‖ ≥ γ − δ} for a number δ > 0, once it enters the benign region. Denote cconverge := 1 − 2c̃0
1−2ηcδ ( ρ‖w∗‖−γ ) , where cδ and c̃0 are some constants that satisfy 0.5 > cδ > 0 and c̃0 > 0. Suppose that there is a number R, such that the size of the iterate satisfies ‖wt‖ ≤ R for all t during the execution of the algorithm and that ‖w∗‖ ≤ R. Denote L := ‖A‖2 + 2ρR. Also, suppose that the momentum parameter β ∈ [0, 0.65] and that
the step size η satisfies η ≤ min (
1 4(‖A∗‖+26ρR) ,
c̃0cδ ( ρ‖w∗‖−γ ) (50L+2(26)2)(‖A∗‖+ρR)(1+‖A∗‖+ρR)+1.3L ) . Set
w0 = w−1 = −r b‖b‖ for any sufficiently small r > 0. Then, in the benign region it takes at most
t . T̂ := 1 ηcδ(ρ‖w∗‖−γ)(1+ βcconverge ) log( (‖A‖2+2ρR)c̃2 )
number of iterations to reach an -approximate error, where c̃ := 4R2(1 + ηC̃) with C̃ = 4 3 c̃0(ρ‖w∗‖ − γ) + 3 ( c̃0cδ(ρ‖w∗‖ − γ) + 10 max{0, γ} ) .
Note that the constraint of η ensures that ηcδ ( ρ‖w∗‖ − γ ) ≤ 18 ; as a consequence, c
converge ≤ 1. Theorem 3 indicates that up to an upper-threshold, a larger value of β reduces the number of iterations to linearly converge to an -optimal point, and hence leads to a faster convergence.
Let us now make a few remarks. First, we want to emphasize that in the linear convergence rate regime, a constant factor improvement of the convergence rate (i.e. of T̂ here) means that the slope of the curve in the log plot of optimization value vs. iteration is steeper. Our experimental result
(Figure 2) confirms this. In this figure, we can see that the curve corresponds to a larger momentum parameter β has a steeper slope than that of the smaller ones. The slope is steeper as β increases, which justifies the effectiveness of the momentum in the linear convergence regime. Our theoretical result also indicates that the acceleration due to the use of momentum is more evident for a small step size η. When η is sufficiently small, the number cconverge are close to 1, which means that the number of iterations can be reduced approximately by a factor of 1 + β.
Secondly, from the theorem, one will need |b(1)| be non-zero, which can be guaranteed with a high probability by adding some Gaussian perturbation on b. We refer the readers to Section 4.2 of Carmon & Duchi (2019) for the technique.
To prove Theorem 3, we will need a series of lemmas which is given in the following subsection.

heading: F.1 SOME SUPPORTING LEMMAS
text: Lemma 12. Denote sequences dt := ρ2 (‖w∗‖ − ‖wt‖)‖wt − w∗‖ 2, et := −(wt − w∗)>(A∗ − 2ηc̃wA 2 ∗)(wt − w∗), gt := − ρ 2 (‖w∗‖ − ‖wt‖) 2 ( ‖w∗‖ + ‖wt‖ − 4ηρc̃w‖wt‖2 ) , where cw := 2β1−β and c̃w := (1 + cw) + (1 + cw)2. For all t, we have that
〈wt − w∗, wt − wt−1〉+ cw‖wt − wt−1‖2 ≤ η t−1∑ s=1 βt−1−s(ds + es + gs).
Proof. We use induction for the proof. The base case t = 0 holds, because w0 = w−1 by initialization and both sides of the inequality is 0. Let us assume that it holds at iteration t. That is, 〈wt−w∗, wt−wt−1〉+cw‖wt−wt−1‖2 ≤ η ∑t−1 s=1 β
t−1−s(ds+es+gs). Consider iteration t+1. We want to prove that
〈wt+1 − w∗, wt+1 − wt〉+ cw‖wt+1 − wt‖2 ≤ η t∑
s=1
βt−s(ds + es + gs).
Denote ∆ := wt+1 − wt. It is equivalent to showing that 〈∆, wt + ∆ − w∗〉 + cw‖∆‖2 ≤ η ∑t s=1 β t−s(ds + es + gs), or
〈−η∇f(wt) + β(wt − wt−1), wt − w∗ − η∇f(wt) + β(wt − wt−1)〉
+ cw‖ − η∇f(wt) + β(wt − wt−1)‖2 ≤ η t∑
s=1
βt−s(ds + es + gs). (67)
which is in turn equivalent to showing that
− η 〈∇f(wt), wt − w∗〉︸ ︷︷ ︸ (a) +η2(1 + cw) ‖∇f(wt)‖2︸ ︷︷ ︸ (b) −2ηβ(1 + cw)〈∇f(wt), wt − wt−1〉︸ ︷︷ ︸ (c)
+ β〈wt − wt−1, wt − w∗〉+ β2‖wt − wt−1‖2 + cwβ2‖wt − wt−1‖2
≤ η t∑
s=1
βt−s(ds + es + gs).
(68)
For term (a), we have that 〈wt − w∗,∇f(wt)〉 = (wt − w∗)>A∗(wt − w∗) + ρ(‖wt‖ − ‖w∗‖)(‖wt‖2 − w>∗ wt)
= (wt − w∗)> ( A∗ + ρ
2 (‖wt‖ − ‖w∗‖)Id
) (wt − w∗) + ρ
2
( ‖w∗‖ − ‖wt‖ )2 (‖wt‖+ ‖w∗‖).
= (wt − w∗)>A∗(wt − w∗) + ρ
2 (‖wt‖ − ‖w∗‖)‖wt − w∗‖2 +
ρ
2
( ‖w∗‖ − ‖wt‖ )2 (‖wt‖+ ‖w∗‖).
(69)
Notice that A∗ −γ+ ρ‖w∗‖Id 0Id, as ρ‖w∗‖ ≥ γ. On the other hand, for term (b), we get that ‖∇f(wt)‖2 = ‖A∗(wt − w∗)− ρ(‖w∗‖ − ‖wt‖)wt‖2
≤ 2(wt − w∗)>A2∗(wt − w∗) + 2ρ2 ( ‖w∗‖ − ‖wt‖ )2‖wt‖2. (70)
For term (c), we can bound it as
− 2ηβ(1 + cw)〈∇f(wt), wt − wt−1〉 ≤ ‖ √ 2η(1 + cw)∇f(wt)‖‖ √ 2β(wt − wt−1)‖
≤ 1 2
( 2η2(1 + cw) 2‖∇f(wt)‖2 + 2β2‖wt − wt−1‖2 )
= η2(1 + cw) 2‖∇f(wt)‖2 + β2‖wt − wt−1‖2.
(71)
Combining the above, we have that
− η〈∇f(wt), wt − w∗〉+ η2(1 + cw)‖∇f(wt)‖2 − 2ηβ(1 + cw)〈∇f(wt), wt − wt−1〉 + β〈wt − wt−1, wt − w∗〉+ β2‖wt − wt−1‖2 + cwβ2‖wt − wt−1‖2
≤ −η〈∇f(wt), wt − w∗〉+ η2c̃w‖∇f(wt)‖2 + β〈wt − wt−1, wt − w∗〉+ β2(2 + cw)‖wt − wt−1‖2 ≤ −η(wt − w∗)> ( A∗ − 2ηc̃wA2∗ ) (wt − w∗) − η ρ 2 (‖w∗‖ − ‖wt‖)2 ( ‖w∗‖+ ‖wt‖ − 4ηρc̃w‖wt‖2 ) + η ρ 2 (‖w∗‖ − ‖wt‖)‖wt − w∗‖2 + β〈wt − wt−1, wt − w∗〉+ β2(2 + cw)‖wt − wt−1‖2︸ ︷︷ ︸ =β〈wt−wt−1,wt−w∗〉+βcw‖wt−wt−1‖2
≤ η(dt + et + gt) + ηβ t−1∑ s=1 βt−1−s(ds + es + gs) := η t∑ s=1 βt−s(ds + es + gs).
(72)
where in the second to last inequality we used β(2 + cw) = cw and the last inequality we used the assumption at iteration t. So we have completed the induction.
Lemma 13. Assume that for all t, ‖wt‖ ≤ R for some number R. Following the notations used in Lemma 12, we denote sequences dt := ρ2 (‖w∗‖ − ‖wt‖)‖wt − w∗‖
2, et := −(wt − w∗)>(A∗ − 2ηc̃wA 2 ∗)(wt − w∗), gt := − ρ 2 (‖w∗‖ − ‖wt‖) 2 ( ‖w∗‖ + ‖wt‖ − 4ηρc̃w‖wt‖2 ) , where cw := 2β1−β and c̃w := (1 + cw) + (1 + cw)2. Let us also denote cβ := (2β2 + 4β + Lβ), L := ‖A‖2 + 2ρR, and zt := ∑t s=0 β t−s+1∇f(ws). If η satisfies: (1) η ≤ 1+β2ρR , and (2) η ≤ 1+β
4‖A∗‖2 , then we have that for all t,
‖wt+1 − w∗‖2
≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2
+ η2(wt−1 − w∗)>(2cβA2∗ + 2β2LId + 2cβρ2R2Id)(wt−1 − w∗) + η2(cβ − 2cw)‖zt−2‖2
+ 2ηβ2 t−2∑ s=0 βt−2−s(ds + es + gs).
(73)
Proof. Recall that the update of heavy ball is
wt+1 = wt − η∇f(wt) + β(wt − wt−1). (74)
So the distance term can be decomposed as ‖wt+1 − w∗‖2 =‖wt − η∇f(wt) + β(wt − wt−1)− w∗‖2
=‖wt − w∗‖2 − 2η〈wt − w∗,∇f(wt)〉+ η2‖∇f(wt)‖2
+ β2‖wt − wt−1‖2 − 2ηβ(wt − wt−1)>∇f(wt) + 2β(wt − w∗)>(wt − wt−1) =‖wt − w∗‖2 − 2η〈wt − w∗,∇f(wt)〉+ η2‖∇f(wt)‖2
+ β2‖wt − wt−1‖2 − 2ηβ〈wt − w∗,∇f(wt)〉 − 2ηβ〈w∗ − wt−1,∇f(wt)〉 + 2β(wt − w∗)>(wt − wt−1)
=‖wt − w∗‖2 − 2η(1 + β) 〈wt − w∗,∇f(wt)〉︸ ︷︷ ︸ (a) +η2 ‖∇f(wt)‖2︸ ︷︷ ︸ (b)
+ (β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈w∗ − wt−1,∇f(wt)〉+ 2β(wt−1 − w∗)>(wt − wt−1). (75)
For term (a), 〈wt−w∗,∇f(wt)〉, by using that∇f(w) = A∗(w−w∗)− ρ(‖w∗‖−‖w‖)w, we can bound it as 〈wt − w∗,∇f(wt)〉 = (wt − w∗)>A∗(wt − w∗) + ρ(‖wt‖ − ‖w∗‖)(‖wt‖2 − w>∗ wt)
= (wt − w∗)> ( A∗ + ρ
2 (‖wt‖ − ‖w∗‖)Id
) (wt − w∗) + ρ
2
( ‖w∗‖ − ‖wt‖ )2 (‖wt‖+ ‖w∗‖).
(76) On the other hand, for term (b), ‖∇f(wt)‖2, we get that
‖∇f(wt)‖2 = ‖A∗(wt − w∗)− ρ(‖w∗‖ − ‖wt‖)wt‖2 ≤ 2(wt − w∗)>A2∗(wt − w∗) + 2ρ2 ( ‖w∗‖ − ‖wt‖ )2‖wt‖2. (77) By combining (75),(76), and (77), we can bound the distance term as ‖wt+1 − w∗‖2 = ‖wt − w∗‖2 − 2η(1 + β)〈wt − w∗,∇f(wt)〉+ η2‖∇f(wt)‖ + (β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈w∗ − wt−1,∇f(wt)〉+ 2β(wt−1 − w∗)>(wt − wt−1) ≤ ‖wt − w∗‖2
− 2η(1 + β) ( (wt − w∗)> ( A∗ + ρ
2 (‖wt‖ − ‖w∗‖)Id
) (wt − w∗) ) − ρη(1 + β) ( ‖w∗‖ − ‖wt‖ )2 (‖wt‖+ ‖w∗‖)
+ η2 ( 2(wt − w∗)>A2∗(wt − w∗) + 2ρ2 ( ‖w∗‖ − ‖wt‖ )2‖wt‖2) + (β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈w∗ − wt−1,∇f(wt)〉+ 2β(wt−1 − w∗)>(wt − wt−1).
(78)
Now let us bound the terms on the last line of (78). For the second to the last term, it is equal to − 2ηβ〈w∗ − wt−1,∇f(wt)〉 = −2ηβ〈w∗ − wt−1,∇f(wt−1)〉 − 2ηβ〈w∗ − wt−1,∇f(wt)−∇f(wt−1)〉, (79)
On the other hand, the last term of (78) is equal to 2β(wt−1 − w∗)>(wt − wt−1) = −2βη〈wt−1 − w∗,∇f(wt−1)〉+ 2β2〈wt−1 − w∗, wt−1 − wt−2〉. (80) By Lemma 12, for all t, we have that
〈wt−1 − w∗, wt−1 − wt−2〉 ≤ −cw‖wt−1 − wt−2‖2 + η t−2∑ s=1 βt−2−s(ds + es + gs). (81)
Therefore, by combining (79), (80), and (81), we have that (β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈w∗ − wt−1,∇f(wt)〉+ 2β(wt−1 − w∗)>(wt − wt−1) ≤ (β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈wt−1 − w∗,∇f(wt−1)−∇f(wt)〉 − 2cwβ2‖wt−1 − wt−2‖2
+ 2ηβ2 t−2∑ s=1 βt−2−s(ds + es + gs).
(82)
Combining (78) and (82) leads to the following, ‖wt+1 − w∗‖2 ≤ ‖wt − w∗‖2 − 2η(1 + β) ( (wt − w∗)> ( A∗ + ρ
2 (‖wt‖ − ‖w∗‖)Id
) (wt − w∗) ) − ρη(1 + β) ( ‖w∗‖ − ‖wt‖ )2 (‖wt‖+ ‖w∗‖)
+ η2 ( 2(wt − w∗)>A2∗(wt − w∗) + 2ρ2 ( ‖w∗‖ − ‖wt‖ )2‖wt‖2) + (β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈wt−1 − w∗,∇f(wt−1)−∇f(wt)〉.
− 2cwβ2‖wt−1 − wt−2‖2 + 2ηβ2 t−2∑ s=1 βt−2−s(ds + es + gs).
(83)
Let us now bound the terms on the second to the last line (83) above,
(β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈wt−1 − w∗,∇f(wt−1)−∇f(wt)〉. (84) First, note that ‖∇2f(w)‖ ≤ ‖A‖2 + 2ρ‖w‖. So we know that f is L := ‖A‖2 + 2ρR smooth on {w : ‖w‖ ≤ R}. Second, denote
zt := t∑ s=0 βt−s+1∇f(ws), (85)
we can bound ‖wt − wt−1‖ as
‖wt − wt−1‖ = ‖ − η t−1∑ s=0 βt−1−s∇f(ws)‖ ≤ η‖∇f(wt−1)‖+ η‖ t−2∑ s=0 βt−1−s∇f(ws)‖
:= η‖∇f(wt−1)‖+ η‖zt−2‖,
‖wt − wt−1‖2 = ‖ − η t−1∑ s=0 βt−1−s∇f(ws)‖2 ≤ 2η2‖∇f(wt−1)‖2 + 2η2‖ t−2∑ s=0 βt−1−s∇f(ws)‖2
:= 2η2‖∇f(wt−1)‖2 + 2η2‖zt−2‖2, (86)
Using the results, we can bound (84) as
(β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈wt−1 − w∗,∇f(wt−1)−∇f(wt)〉 (86)
≤ 2η2(β2 + 2β)(‖∇f(wt−1)‖2 + ‖zt−2‖2) + 2η2βL‖wt−1 − w∗‖(‖∇f(wt−1)‖+ ‖zt−2‖) ≤ 2η2(β2 + 2β)(‖∇f(wt−1)‖2 + ‖zt−2‖2) + η2βL(‖wt−1 − w∗‖2 + ‖∇f(wt−1)‖2) + η2βL(‖wt−1 − w∗‖2 + ‖zt−2‖2) = 2η2(β2 + 2β)(‖∇f(wt−1)‖2 + ‖zt−2‖2) + 2η2βL‖wt−1 − w∗‖2 + η2βL‖∇f(wt−1)‖2
+ η2βL‖zt−2‖2. (87)
Note that
‖∇f(wt−1)‖2 = ‖A∗(wt−1 − w∗)− ρ(‖w∗‖ − ‖wt−1‖)wt−1‖2 ≤ 2(wt−1 − w∗)>A2∗(wt−1 − w∗) + 2ρ2 ( ‖w∗‖ − ‖wt−1‖ )2‖wt−1‖2. (88) Denote
cβ := (2β 2 + 4β + Lβ). (89)
By (87), (88), we have that
(β2 + 2β)‖wt − wt−1‖2 − 2ηβ〈wt−1 − w∗,∇f(wt−1)−∇f(wt)〉. ≤ 2η2cβ(wt−1 − w∗)>A2∗(wt−1 − w∗) + 2η2cβρ2 ( ‖w∗‖ − ‖wt−1‖ )2‖wt−1‖2 + 2η2β2L‖wt−1 − w∗‖2 + η2cβ‖zt−2‖2.
(90)
Let us summarize the results so far, by (83) and (90), we have that
‖wt+1 − w∗‖2 ≤ ‖wt − w∗‖2 − 2η(1 + β) ( (wt − w∗)> ( A∗ + ρ
2 (‖wt‖ − ‖w∗‖)Id
) (wt − w∗) ) − ρη(1 + β) ( ‖w∗‖ − ‖wt‖ )2 (‖wt‖+ ‖w∗‖)
+ η2 ( 2(wt − w∗)>A2∗(wt − w∗) + 2ρ2 ( ‖w∗‖ − ‖wt‖ )2‖wt‖2) + 2η2cβ(wt−1 − w∗)>A2∗(wt−1 − w∗) + 2η2cβρ2 ( ‖w∗‖ − ‖wt−1‖
)2‖wt−1‖2 + 2η2β2L‖wt−1 − w∗‖2 + η2(cβ − 2cw)‖zt−2‖2 + 2ηβ2
t−2∑ s=1 βt−2−s(ds + es + gs),
(91)
where we also used that wt−1 − wt−2 = − ηβ zt−2. We can rewrite the inequality above further as
‖wt+1 − w∗‖2 ≤ (wt − w∗)> ( Id − 2η(1 + β)A∗(Id − η
1 + β A∗)− ηρ(1 + β)(‖wt‖ − ‖w∗‖)Id
) (wt − w∗)
−ηρ(1 + β)(‖w∗‖ − ‖wt‖)2 ( ‖wt‖(1−
2ηρ‖wt‖ 1 + β
) + ‖w∗‖ )
︸ ︷︷ ︸ ≤0
+2η2cβρ 2 ( ‖w∗‖ − ‖wt−1‖ )2‖wt−1‖2 + η2(wt−1 − w∗)>(2cβA2∗ + 2β2LId)(wt−1 − w∗)
+ η2(cβ − 2cw)‖zt−2‖2 + 2ηβ2 t−2∑ s=1 βt−2−s(ds + es + gs),
(92) where we used that −ηρ(1 + β)(‖w∗‖ − ‖wt‖)2 ( ‖wt‖(1 − 2ηρ‖wt‖1+β ) + ‖w∗‖ ) ≤ 0 for any η that satisfies η ≤ 1 + β
2ρR , (93)
as ‖wt‖ ≤ R for all t. Let us simplify the inequality (92) further by writing it as ‖wt+1 − w∗‖2
(a) ≤ (wt − w∗)> ( Id − 2η(1 + β)A∗(Id − η
1 + β A∗)− ηρ(1 + β)(‖wt‖ − ‖w∗‖)Id
) (wt − w∗)
+ η2(wt−1 − w∗)>(2cβA2∗ + 2β2LId + 2cβρ2R2Id)(wt−1 − w∗) + η2(cβ − 2cw)‖zt−2‖2
+ 2ηβ2 t−2∑ s=1 βt−2−s(ds + es + gs).
(b) ≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2
+ η2(wt−1 − w∗)>(2cβA2∗ + 2β2LId + 2cβρ2R2Id)(wt−1 − w∗) + η2(cβ − 2cw)‖zt−2‖2
+ 2ηβ2 t−2∑ s=1 βt−2−s(ds + es + gs),
(94)
where (a) is because that (‖w∗‖−‖wt−1‖)2‖wt−1‖2 ≤ ‖w∗−wt−1‖2‖wt−1‖2 ≤ ‖w∗−wt−1‖2R2 as ‖wt‖ ≤ R for all t, while (b) is by another constraint of η,
η ≤ 1 + β 4‖A∗‖2 , (95)
so that 2ηA∗(Id − η1+βA∗) 3 2ηA∗ 3 2η(−γ + ρ‖w∗‖)Id.
Lemma 14. Fix some numbers c0, c1 > 0. Denote ωβt−2 := ∑t−2 s=0 β. Following the notations and assumptions used in Lemma 12 and Lemma 13. If η satisfies: (1) η ≤ c0β 2cβ‖A∗‖22+2β2L+2cβρ2R2 , (2)
η ≤ c1 β(2c̃w+Lβω β t−2)‖A∗‖ , (3) η ≤ c1 ρβ2Lωβt−2R , (4) η ≤ 14ρc̃wR , and (5) η ≤ c1
( ‖A∗‖+ρR ) Lβ2ωβt−2(‖A∗‖2+ρ2R2) for
all t, then we have that for all t, ‖wt+1 − w∗‖2 ≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2
+ 2ηβc1 ( ‖A∗‖+ ρR ) t−2∑ s=0 βt−2−s‖ws − w∗‖2
− 2ηβ2 t−2∑ s=1 βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖wt‖ ) Id ) (ws − w∗).
Proof. From Lemma 13, we have that
‖wt+1 − w∗‖2
≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2
+ η2(wt−1 − w∗)>(2cβA2∗ + 2β2LId + 2cβρ2R2Id)(wt−1 − w∗) + η2(cβ − 2cw)‖zt−2‖2
+ 2ηβ2 t−2∑ s=1 βt−2−s(ds + es + gs),
(a) ≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2 + η2(cβ − 2cw)‖zt−2‖2
+ 2ηβ2 t−2∑ s=1 βt−2−s(ds + es + gs),
(b) ≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2
+ 2η2(cβ − 2cw)β2ωβt−2 t−2∑ s=0 βt−2−s(ws − w∗)>(A2∗ + ρ2R2Id)(ws − w∗)
+ 2ηβ2 t−2∑ s=1 βt−2−s(ds + es + gs),
(96)
where (a) is by a constraint of η so that η(2cβ‖A∗‖22 + 2β2L+ 2cβρ2R2) ≤ c0β, and (b) is due to the following, denote ωβt := ∑t s=0 β, we have that
‖zt‖2 := ‖ t∑
s=0
βt−s+1∇f(ws)‖2 = β2(ωβt )2‖ t∑
s=0
βt−s
ωβt ∇f(ws)‖2
≤ β2(ωβt )2 ( t∑ s=0 βt−s ωβt ‖∇f(ws)‖2 ) ≤ 2β2ωβt
( t∑ s=0 βt−s(ws − w∗)>(A2∗ + ρ2R2Id)(ws − w∗) ) .
(97)
where the first inequality of (97) is due to Jensen’s inequality and the second inequality of (97) is due to the following upper-bound of the gradient norm, ‖∇f(ws)‖2 = ‖A∗(ws −w∗)− ρ(‖w∗‖ − ‖ws‖)ws‖2 ≤ 2(ws − w∗)>A2∗(ws − w∗) + 2ρ2(‖w∗‖ − ‖ws‖)2‖ws‖2 ≤ 2(ws − w∗)>A2∗(ws −
w∗) + 2ρ 2‖w∗ − ws‖2R2. By the definitions of ds, es, gs, we further have that
‖wt+1 − w∗‖2 ≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖ws − w∗‖2 + ηβc0‖wt−1 − w∗‖2
− 2ηβ2 t−2∑ s=1 βt−2−s(ws − w∗)> ( A∗ − ρ 2 (‖w∗‖ − ‖ws‖)Id ) (ws − w∗)
+ 2ηβ2 t−2∑ s=1 βt−2−s(ws − w∗)> ( η(2c̃w + (cβ − 2cw)ωβt−2)A2∗ ) (ws − w∗)
+ 2η2ρ2R2β2ωβt−2(cβ − 2cw) t−2∑ s=1 βt−2−s‖ws − w∗‖2
− ηβ2ρ t−2∑ s=1 βt−2−s(‖w∗‖ − ‖ws‖)2 ( ‖w∗‖+ ‖ws‖ − 4ηρc̃w‖ws‖2 ) + η2(cβ − 2cw)2βtωβt−2 ( ‖A∗‖22 + ρ2R2)‖w0 − w∗‖2
≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2
− 2ηβ2 t−2∑ s=1 βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖ws‖ ) Id ) (ws − w∗)
+ 2ηβc1 t−2∑ s=1 βt−2−s(ws − w∗)>(A∗ + ρRId)(ws − w∗)
− ηβ2ρ t−2∑ s=1 βt−2−s(‖w∗‖ − ‖ws‖)2 ( ‖w∗‖+ ‖ws‖ − 4ηρc̃w‖ws‖2 ) + η2(cβ − 2cw)2βtωβt−2 ( ‖A∗‖22 + ρ2R2)‖w0 − w∗‖2,
(98)
where the last inequality is due to (1): 2η2β2(2c̃w + (cβ − 2cw)ωβt−2)A2∗ 2η2β2(2c̃w + Lβωβt−2)A
2 ∗ 2ηβc1A∗, as cβ − 2cw ≤ Lβ and that η ≤ c1β(2c̃w+Lβωβt−2)‖A∗‖ , and (2) that
2η2ρ2R2β2ωβt−2(cβ − 2cw) ≤ 2η2ρ2R2Lβ3ω β t−2 ≤ 2ηρβc1R, as η ≤ c1ρβ2Lωβt−2R .
To continue, let us bound the last two terms on (98). For the second to the last term, by using that ‖wt‖ ≤ R for all t and that η ≤ 14ρc̃wR , we have that ( ‖w∗‖ + ‖ws‖ − 4ηρc̃w‖ws‖2 ) ≥ ‖w∗‖. Therefore, we have that the second to last term on (98) is non-positive, namely, −ηβ2ρ ∑t−2 s=1 β t−2−s(‖w∗‖−‖ws‖)2 ( ‖w∗‖+‖ws‖−4ηρc̃w‖ws‖2 ) ≤ 0. For the last term on (98),
by using that η ≤ c1 ( ‖A∗‖+ρR ) Lβ2ωβt−2(‖A∗‖2+ρ2R2) , we have that η2(cβ−2cw)2βtωβt−2 ( ‖A∗‖22 +ρ2R2)‖w0−
w∗‖2 ≤ η22Lβt+1ωβt−2 ( ‖A∗‖22 + ρ2R2)‖w0 − w∗‖2 ≤ 2ηβt−1c1 ( ‖A∗‖+ ρR ) ‖w0 − w∗‖2.
Combining the above results, we have that ‖wt+1 − w∗‖2 ≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2
+ 2ηβc1 ( ‖A∗‖+ ρR ) t−2∑ s=0 βt−2−s‖ws − w∗‖2
− 2ηβ2 t−2∑ s=1 βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖ws‖ ) Id ) (ws − w∗).
(99)
The following lemma will be used for getting the iteration complexity from Lemma 14. Lemma 15. For a non-negative sequence {yt}t≥0, suppose that it satisfies
yt+1 ≤ pyt + qyt−1 + rβȳt−2 + xβ t−2∑ s=τ βt−2−sys + zβ t−τ+1, for all t ≥ τ, (100)
for non-negative numbers p, q, r, x, z ≥ 0 and β ∈ [0, 1), where we denote
ȳt := t∑ s=0 βt−sys. (101)
Fix some numbers φ ∈ (0, 1) and ψ > 0. Define θ = −p+ √ p2+4(q+φ)
2 . Suppose that
β ≤ p+ √ p2+4(q+φ) 2 φ
r + φ+ x and β ≤
p+ √ p2 + 4(q + φ)
2 − 1 ψ .
then we have that for all t ≥ τ , yt ≤ (p+√p2 + 4(q + φ)
2
)t−τ cτ,β , (102)
where cτ,β is an upper bound that satisfies, for any t ≤ τ , yt + θyt−1 + φȳt−2 + βψz ≤ cτ,β , (103)
Proof. For a non-negative sequence {yt}t≥0, suppose that it satisfies
yt+1 ≤ pyt + qyt−1 + rβȳt−2 + xβ t−2∑ s=τ βt−2−sys + zβ t−τ+1, for all t ≥ τ, (104)
for non-negative numbers p, q, r, x, z and β ∈ [0, 1), where we denote
ȳt := t∑ s=0 βt−sys. (105)
yt+1 + θyt + φȳt−1 + ψzβ t−τ+2
(a) ≤ (p+ θ)yt + qyt−1 + rβȳt−2 + φȳt−1 + xβ t−2∑ s=τ βt−2−sys + (ψβ + 1)zβ t−τ+1
(b) = (p+ θ)yt + (q + φ)yt−1 + (rβ + φβ)ȳt−2 + xβ t−2∑ s=τ βt−2−sys + (ψβ + 1)zβ t−τ+1
≤ (p+ θ)(yt + (q + φ)
p+ θ yt−1) + β(r + φ+ x)ȳt−2 + (ψβ + 1)zβ
t−τ+1
(c) ≤ (p+ θ)(yt + θyt−1) + β(r + φ+ x)ȳt−2 + (ψβ + 1)zβt−τ+1
(d) ≤ (p+ θ)(yt + θyt−1 + φȳt−2 + ψzβt−τ+1) ≤ (p+ θ)2(yt−1 + θyt−2 + φȳt−3 + ψzβt−τ−1) ≤ . . . := (p+ θ)t−τ+1cτ,β .
(106)
where (a) is due to the dynamics (100) , (b) is because ȳt−1 = yt−1 + βȳt−2, (c) is by,
q + φ p+ θ ≤ θ, (107)
and (d) is by
β ≤ (p+ θ)φ r + φ+ x and β ≤ p+ θ − 1 ψ . (108)
Note that (107) holds if,
θ ≥ −p+
√ p2 + 4(q + φ)
2 . (109)
Let us choose the minimal θ = −p+ √ p2+4(q+φ) 2 . So we have p + θ = p+ √ p2+4(q+φ)
2 , which completes the proof.
Lemma 16. Assume that for all t, ‖wt‖ ≤ R for some number R. Fix the numbers c0, c1 > 0 in Lemma 14 so that c1 ≤ c040(‖A∗‖+ρR) . Suppose that the step size η satisfies (1) η ≤ c0β 2cβ‖A∗‖22+2β2L+2cβρ2R2 , (2) η ≤ c1 β(2c̃w+Lβω β t−2)‖A∗‖ , (3) η ≤ c1 ρβ2Lωβt−2R , (4) η ≤ 14ρc̃wR , (5)
η ≤ c1 ( ‖A∗‖+ρR ) Lβ2ωβt−2(‖A∗‖2+ρ2R2) , (6) η ≤ 1+β2ρR , and (7) η ≤ 1+β 4‖A∗‖2 for all t, where cβ := (2β 2+4β+Lβ), L := ‖A‖2+2ρR, cw := 2β1−β , c̃w := (1+cw)+(1+cw) 2 and ωβt−2 := ∑t−2 s=0 β. Furthermore, sup-
pose that the momentum parameter β satisfies β ≤ min{ 1011 ( 1− η(1 + β)δ ) , 1− η(1 + β)δ− 0.1 } . Fix numbers δ, δ′ so that 12 (ρ‖w∗‖ − γ) > δ > 0 and c0 20 ≥ δ
′ > 0. Assume that ‖wt‖ is nondecreasing. If ρ‖wt‖ ≥ γ − 12 (ρ‖w∗‖ − γ) + δ and ρ‖wt‖ ≥ γ − δ
′ for some t ≥ t∗, then we have that
‖wt − w∗‖2 ≤ ( 1− η(1 + β)δ + 2ηβc0 1− η(1 + β)δ )t−t∗ cβ,t∗
where cβ,t∗ is a number that satisfies for any t ≤ t∗, ‖wt − w∗‖2 + 2ηc0β 1−η(1+β)δ‖wt−1 − w∗‖ 2 + ηβc0 ∑t−2 s=0 β t−2−s‖ws−w∗‖2 +max{0, β20η ∑t∗−1 s=1 β t∗−1−s ( γ 2 − ρ 2‖ws‖ ) ‖ws−w∗‖2} ≤ cβ,t∗ .
Proof. From Lemma 14, we have that ‖wt+1 − w∗‖2 ≤ ( 1− η(1 + β) [ ρ‖wt‖ − (γ −
ρ‖w∗‖ − γ 2
) ]) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2
+ 2ηβc1(‖A‖∗ + ρR) t−2∑ s=0 βt−2−s‖ws − w∗‖2
− 2ηβ2 t−2∑ s=1 βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖ws‖ ) Id ) (ws − w∗).
Using that for all t ≥ t∗, ρ‖wt‖ ≥ γ − 12 (ρ‖w∗‖ − γ) + δ, we have that ‖wt+1 − w∗‖2 ≤ ( 1− η(1 + β)δ ) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2
+ 2ηβc1(‖A‖∗ + ρR) t−2∑ s=0 βt−2−s‖ws − w∗‖2
− 2ηβ2 t−2∑ s=1 βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖ws‖ ) Id ) (ws − w∗)
= ( 1− η(1 + β)δ ) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2
+ 2ηβc1(‖A‖∗ + ρR) t−2∑ s=0 βt−2−s‖ws − w∗‖2
− 2ηβ2 t−2∑ s=t∗ βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖ws‖ ) Id ) (ws − w∗)
− 2ηβ2 t∗−1∑ s=1 βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖ws‖ ) Id ) (ws − w∗).
(110)
We bound the second to last term of (110) as follows. Note that for s ≥ t∗, we have that ρ‖ws‖ ≥ γ − δ′ by the assumption that ‖wt‖ is non-decreasing. Therefore, once ρ‖wt‖ exceeds any level
blow γ − δ′, it will not fall below γ − δ′. So we have that
A∗ − ρ
2
( ‖w∗‖ − ‖ws‖ ) Id (a) A∗ − ρ
2 ‖w∗‖Id +
γ − δ′
2 Id
(b)
(−γ + ρ‖w∗‖)Id − ρ
2 ‖w∗‖Id +
γ − δ′
2 Id
(c) −δ ′
2 Id
(111)
where (a) uses that ρ‖ws‖ ≥ γ − δ′ for some number δ′ > 0, (b) uses the fact that A∗ (−γ + ρ‖w∗‖)Id, and (c) uses the fact that ρ‖w∗‖ ≥ γ. Using the result, we can bound the second to the last term of (110) as
− 2ηβ2 t−2∑ s=t∗ βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖ws‖ ) Id ) (ws − w∗)
≤ ηβ2δ′ t−2∑ s=t∗ βt−2−s‖ws − w∗‖2.
(112)
Now let us switch to the last term of (110). Using the fact that A∗ (−γ + ρ‖w∗‖)Id, and that ρ‖w∗‖ ≥ γ by the characterization of the optimizer ‖w∗‖, we have that
A∗ − ρ
2
( ‖w∗‖ − ‖ws‖ ) Id ρ
2
( ‖w∗‖+ ‖ws‖ ) Id − γ (ρ 2 ‖ws‖ − γ 2 ) Id. (113)
So we can bound the last term as
− 2ηβ2 t∗−1∑ s=1 βt−2−s(ws − w∗)> ( A∗ − ρ 2 ( ‖w∗‖ − ‖ws‖ ) Id ) (ws − w∗)
≤ 2ηβt−t∗+1 t∗−1∑ s=1 βt∗−1−s(ws − w∗)> (γ 2 − ρ 2 ‖ws‖ ) Id(ws − w∗) := 2ηβt−t∗+1Dt∗, (114)
where we denote Dt∗ := ∑t∗−1 s=1 β t∗−1−s(ws−w∗)> ( γ 2 − ρ 2‖ws‖ ) Id(ws−w∗). Combining (110), (112), (114), we have that
‖wt+1 − w∗‖2 ≤ ( 1− η(1 + β)δ ) ‖wt − w∗‖2 + ηβc0‖wt−1 − w∗‖2
+ 2ηβc1(‖A‖∗ + ρR) t−2∑ s=0 βt−2−s‖ws − w∗‖2
+ ηβ2δ′ t−2∑ s=t∗ βt−2−s‖ws − w∗‖2 + 2ηβt−t∗+1Dt∗.
(115)
Now we are ready to use Lemma 15. Set p, q, r, x, z in Lemma 15 as follows.
• p = 1− η(1 + β)δ
• q = ηβc0
• r = 2η(‖A∗‖+ ρR)c1
• x = ηβδ′
• z = 2ηDt∗1{Dt∗ ≥ 0}.
• φ = q = ηβc0
• ψ = 10.
we have that ‖wt − w∗‖2 ≤ (1− η(1 + β)δ +√(1− η(1 + β)δ)2 + 4(2ηβc0)
2
)t−t∗ cβ,t∗
≤ ( 1− η(1 + β)δ + 2ηβc0 1− η(1 + β)δ )t−t∗ cβ,t∗
(116)
where we use that a+a √ 1+4b/a2
2 ≤ a+a(1+ 2b a2 ) 2 = a + b a for number a, b that satisfy 4b/a 2 ≥
−1. Now let us check the conditions of Lemma 15 to see if β ≤ p+ √ p2+4(q+φ) 2 φ
r+φ+x and β ≤ p+ √ p2+4(q+φ)
2 − 1 ψ . For the first condition, β(r + φ + x) ≤
p+ √ p2+4(q+φ)
2 φ is equivalent to β ( 2η(‖A∗‖ + ρR)c1 + ηβc0 + ηβδ′ ) ≤ p+ √ p2+4(q+φ)
2 ηβc0, which can be satisfied by
c1 ≤ c040(‖A∗‖+ρR) and δ ′ ≤ c020 , leading to an upper bound of β, which is β ≤ 10 11
p+ √ p2+4(q+φ)
2 . Using the expression of p, q, φ, and ψ, it suffices to have that β ≤ 1011 ( 1− η(1 + β)δ ) . On the other
hand, for the second condition, β ≤ p+ √ p2+4(q+φ)
2 − 1 ψ , by using the expression of p, q, φ, and ψ,
it suffices to have that β ≤ 1− η(1 + β)δ − 0.1.

heading: F.2 PROOF OF THEOREM 3
text: Proof. Let δ = cδ ( ρ‖w∗‖ − γ ) for some number cδ < 0.5. Denote c converge ηδ := 1 − 2c̃0 1−2ηδ for some number c̃0 > 0. By Lemma 16, we have that
‖wt − w∗‖2 ≤ ( 1− η(1 + β)δ + 2ηβc̃0δ 1− η(1 + β)δ )t−t∗ cβ,t∗
≤ ( 1− ηδ(1 + βcconvergeηδ ) )t−t∗ cβ,t∗
(117)
where cβ,t∗ is a number that satisfies for any t ≤ t∗, ‖wt − w∗‖2 + 2ηc0β 1−η(1+β)δ‖wt−1 − w∗‖ 2 + ηβc0 ∑t−2 s=0 β t−2−s‖ws−w∗‖2 +max{0, β20η ∑t∗−1 s=1 β t∗−1−s ( γ 2 − ρ 2‖ws‖ ) ‖ws−w∗‖2} ≤ cβ,t∗ . Note that we can obtain a trivial upper-bound cβ,t for any t as follows. Using that ‖wt−w∗‖2 ≤ 4R2 and that c0 ← c̃0δ and that δ ← cδ ( ρ‖w∗‖ − γ ) , we can upper-bound the term as cβ,t ≤ 4R2 ( 1 + 2ηc̃0cδ ( ρ‖w∗‖ − γ ) β
1− η(1 + β)cδ ( ρ‖w∗‖ − γ ) + ηβc̃0cδ(ρ‖w∗‖ − γ)1− βt 1− β + β10ηmax{0, γ}1− β t 1− β )
≤ 4R2(1 + ηβC̃β) := c̃β .
(118)
where we define C̃β := 2c̃0cδ
( ρ‖w∗‖−γ ) 1−η(1+β)cδ ( ρ‖w∗‖−γ ) + 11−β (c̃0cδ(ρ‖w∗‖−γ)+10 max{0, γ}) and c̃β := 4R2(1 + ηβC̃β).
Now denote cconverge := 1− 2c̃0 1−2ηcδ ( ρ‖w∗‖−γ ) . We have that f(wt∗+t)− f(w∗) (a)
≤ ‖A‖2 + 2ρR 2 ‖wt∗+t − w∗‖2
(b) ≤ ‖A‖2 + 2ρR 2
c̃β exp ( − ηcδ(ρ‖w∗‖ − γ)(1 + βcconverge)t ) (119) where (a) uses the (‖A‖2 + 2ρR)-smoothness of function f(·) in the region of {w : ‖w‖ ≤ R}, and (b) uses (117), (118) and that δ := cδ ( ρ‖w∗‖ − γ ) . So we see that the number of iterations in the linear convergence regime is at most
t ≤ T̂ := 1 ηcδ(ρ‖w∗‖ − γ)(1 + βcconverge) log( (‖A‖2 + 2ρR)c̃β 2 ). (120)
Lastly, let us check if the step size η satisfies the constraints of Lemma 16. Recall the notations that cw := 2β 1−β , c̃w := (1+cw)+(1+cw)
2, cβ := (2β2 +4β+Lβ), and L := ‖A‖2 +2ρR. Lemma 16 has the following constraints, (1) η ≤ c0β
2cβ‖A∗‖22+2β2L+2cβρ2R2 , (2) η ≤ c1β(2c̃w+Lβ/(1−β))‖A∗‖ ,
(3) η ≤ c1ρβ2LR/(1−β) , (4) η ≤ 1 4ρc̃wR , (5) η ≤ c1
( ‖A∗‖+ρR ) Lβ2(‖A∗‖2+ρ2R2)/(1−β) , (6) η ≤ 1+β 2ρR , and (7)
η ≤ 1+β4‖A∗‖2 . For the constraints of (1), using c0 = c̃0δ and that δ = cδ ( ρ‖w∗‖ − γ ) , it can be rewritten as
η ≤ c̃0cδ
( ρ‖w∗‖ − γ ) 2(2β + 4 + L) ( ‖A∗‖22 + ρ2R2 ) + 2βL . (121)
For the constraints of (2), using c1 ≤ c̃0δ40(‖A∗‖+ρR) and that δ = cδ ( ρ‖w∗‖ − γ ) , it can be rewritten as
η ≤ c̃0cδ
( ρ‖w∗‖ − γ ) 40β(‖A∗‖2 + ρR‖A∗‖)(2c̃w + Lβ/(1− β)) . (122)
The constraints of (3) can be written as, using c1 ≤ c̃0δ40(‖A∗‖+ρR) and that δ = cδ ( ρ‖w∗‖ − γ ) ,
η ≤ c̃0cδ
( ρ‖w∗‖ − γ ) 40ρR(‖A∗‖+ ρR)β2L/(1− β) . (123)
The constraints of (5) translates into, using c1 ≤ c̃0δ40(‖A∗‖+ρR) and that δ = cδ ( ρ‖w∗‖ − γ ) ,
η ≤ c̃0cδ
( ρ‖w∗‖ − γ ) Lβ2(‖A∗‖2 + ρ2R2)/(1− β) . (124)
Considering all the above constraints, it suffices to let η satisfies
η ≤ min ( 1
4(‖A∗‖+ c̃wρR) ,
c̃0cδ ( ρ‖w∗‖ − γ ) Cβ(‖A∗‖+ ρR)(1 + ‖A∗‖+ ρR) + 2βL ) , (125)
where Cβ := max ( 4β + 8 + 2L, 40β 2L 1−β + 80βc̃w, 4c̃w ) .
Note that the constraint of η satisfies ηcδ ( ρ‖w∗‖ − γ ) ≤ 14cδ ≤ 1 8 . Using this inequality, we can simplify the constraint regarding the parameter β in Lemma 16, which leads to β ∈ [0, 0.65]. Consequently, we can simplify and upper bound the constants c̃w, Cβ and C̃β , which leads to the theorem statement.
Thus, we have completed the proof.

heading: G MORE DISCUSSIONS
text: Recall the discussion in the main text, we showed that the iterate wt+1 generated by HB satisfies
〈wt+1, ui〉 = (1 + ηλi)〈wt, ui〉+ β(〈wt, ui〉 − 〈wt−1, ui〉), (126)
which is in the form of dynamics shown in Lemma 1. Hence, one might be able to show that with the use of the momentum, the growth rate of the projection on the eigenvector ui (i.e. |〈wt+1, ui〉|) is faster as the momentum parameter β increases. Furthermore, the top eigenvector projection |〈wt+1, u1〉| is the one that grows at the fastest rate. As the result, after normalization (i.e. wT|wT | ), the normalized solution will converge to the top eigenvector after a few iterations T .
However, we know that power iteration or Lanczos method are the standard, specialized, state-ofthe-art algorithms for computing the top eigenvector. It is true that HB is outperformed by these methods. But in the next subsection, we will show an implication of the acceleration result, compared to vanilla gradient descent, of top eigenvector computations.
G.1 IMPLICATION: ESCAPE SADDLE POINTS FASTER
Recent years there is a growing trend in designing algorithms to quickly find a second order stationary point in non-convex optimization (e.g. Carmon et al. (2018); Agarwal et al. (2017); AllenZhu & Li (2018); Xu et al. (2018); Ge et al. (2015); Levy (2016); Fang et al. (2019); Jin et al. (2017; 2018; 2019); Daneshmand et al. (2018); Staib et al. (2019); Wang et al. (2020)). The common assumptions are that the gradient is L-Lipschitz: ‖∇f(x) − ∇f(y)‖ ≤ L‖x − y‖ and that the Hessian is ρ-Lipschitz: ‖∇2f(x) − ∇2f(y)‖ ≤ ρ‖x − y‖, while some related works make some additional assumptions. All the related works agree that if the current iterate wt0 is in the region of strict saddle points, defined as that the gradient is small (i.e. ‖∇f(wt0)‖ ≤ g ) but the least eigenvalue of the Hessian is strictly negative (i.e. λmin(∇2f(wt0)) − hId), then the eigenvector corresponding the least eigenvalue λmin(∇2f(wt0)) is the escape direction. To elaborate, by ρ-Lipschitzness of the Hessian, f(wt0+t) − f(wt0) ≤ 〈∇f(wt0), wt0+t − wt0〉 + 1 2 (wt0+t − wt0)>∇2f(wt0)(wt0+t − wt0)︸ ︷︷ ︸
exhibit negative curvature
+ρ3‖wt0+t−wt0‖ 3. So if wt0+t−wt0 is in the direction
of the bottom eigenvector of ∇2f(wt0), then 12 (wt0+t − wt0) >∇2f(wt0)(wt0+t − wt0) ≤ −c′ h for some c′ > 0. Together with the fact that the gradient is small when in the region of saddle points and the use of a sufficiently small step size can guarantee that the function value decreases sufficiently (i.e.f(wt0+t)− f(wt0) ≤ −c h for some c > 0). Therefore, many related works design fast algorithms by leveraging the problem structure to quickly compute the bottom eigenvector of the Hessian (see e.g. Carmon et al. (2018); Agarwal et al. (2017); Allen-Zhu & Li (2018); Xu et al. (2018)).
An interesting question is as follows “If the Heavy Ball algorithm is used directly to solve a nonconvex optimization problem, can it escape possible saddle points faster than gradient descent?” Before answering the question, let us first conduct an experiment to see if the Heavy Ball algorithm can accelerate the process of escaping saddle points. Specifically, we consider a problem that was consider by Staib et al. (2019); Reddi et al. (2018); Wang et al. (2020) for the challenge of escaping saddle points. The problem is
minw f(w) := 1 n ∑n i=1 ( 1 2w >Hw + x>i w + ‖w‖1010 ) (127)
with H := [ 1 0 0 −0.1 ] . where xi ∼ N (0, diag([0.1, 0.001])) and the small variance in the second
component will provide smaller component of gradient in the escape direction. At the origin, we have that the gradient is small but that the Hessian exhibits a negative curvature. For this problem, Wang et al. (2020) observe that SGD with momentum escapes the saddle points faster but they make strong assumptions in their analysis. We instead consider the Heavy Ball algorithm (i.e. Algorithm 1, the deterministic version of SGD with momentum). Figure 4 shows the result and we see that the higher the momentum parameter β, the faster the process of escaping the saddle points.
We are going to argue that the observation can be explained by our theoretical result that the Heavy Ball algorithm computes the top eigenvector faster than gradient descent. Let us denote h(w) := 1 n ∑n i=1 ( x>i w + ‖w‖1010 ) . We can rewrite the objective (127) as f(w) := 12w >Hw + h(w). Then,
the Heavy Ball algorithm generates the iterate according to
wt+1 =
[ 1− η 0
0 1 + η10 ] ︸ ︷︷ ︸
:=A
wt + β(wt − wt−1)− η∇h(wt). (128)
By setting η ≤ 1, we have that the top eigenvector of A is u1 = e2, which is the escape direction. Now observe the similarity (if one ignores the term η∇h(wt)) between the update (128) and 19. It suggests that a similar analysis could be used for explaining the escape process. In Appendix G.2, we provide a detailed analysis of the observation. However, problem (127) has a fixed Hessian and is a synthetic objective function. For general smooth non-convex optimization, can we have a similar explanation? Consider applying the Heavy Ball algorithm to minw f(w) and suppose that at time t0, the iterate is in the region of strict saddle points. We have that wt+1+t0 − wt0 = ( Id − η∇2f(wt0) )( wt+t0 −
wt0 ) + β ( (wt+t0 − wt0) − (wt+t0−1 − wt0) ) + η ( ∇2f(wt0)(wt+t0 − wt0)−∇f(wt+t0)︸ ︷︷ ︸
deviation
) .
By setting η ≤ 1L with L being the smoothness constant of the problem, we have that the top eigenvector ( Id − η∇2f(wt0) ) is the eigenvector that corresponds to the smallest eigenvalue of the Hessian ∇2f(wt0), which is an escape direction. Therefore, if the deviation term can be controlled, the dynamics of the Heavy Ball algorithm can be viewed as implicitly and approximately computing the eigenvector, and hence we should expect that higher values of momentum parameter accelerate the process. To control ∇2f(wt0)(wt+t0−wt0)−∇f(wt+t0), one might want to exploit the ρ-Lipschitzness assumption of the Hessian and might need further mild assumptions, as Du et al. (2017) provide examples showing that gradient descent can take exponential time T = Ω(exp(d)) to escape saddles points. Previous work of Wang et al. (2020) makes strong assumptions to avoid the result of exponential time to escape.
On the other hand, Lee et al. (2019) show that first order methods can escape strict saddle points almost surely. So we conjecture that under additional mild conditions, if gradient descent can escape in polynomial time, then using Heavy Ball momentum (β 6= 0) will accelerate the process of escape. We leave it as a future work.

heading: G.2 ESCAPING SADDLE POINTS OF (127)
text: Recall the objective is
minw f(w) := 1 n ∑n i=1 ( 1 2w >Hw + x>i w + ‖w‖1010 ) (129)
with
H := [ 1 0 0 −0.1 ] . (130)
where xi ∼ N (0, diag([0.1, 0.001])) and the small variance in the second component will provide smaller component of gradient in the escape direction. At the origin, we have that the gradient is small but that the Hessian exhibits a negative curvature. Let us denote h(w) := 1n ∑n i=1 ( x>i w +
‖w‖1010 ) . We can rewrite the objective as f(w) := 12w
>Hw+h(w). Then, the Heavy Ball algorithm generates the iterate according to
wt+1 = wt − ηHwt + β(wt − wt−1)− η∇h(wt)
=
[ 1− η 0
0 1 + 0.1η
] wt + β(wt − wt−1)− η∇h(wt)
:= Awt + β(wt − wt−1)− η∇h(wt),
(131)
where the matrix A is defined as A := [ 1− η 0
0 1 + 0.1η
] . By setting η ≤ 1, we have that the
top eigen-vector of A is u1 = e2 = [ 0 1 ] , which is the escape direction. In the following, let us
denote u2 := e1 = [ 1 0 ] and denote x̄ = 1n ∑n i=1 xi. Since ‖w‖pp := (∑d i=1 |wi|p ) and that ∂‖w‖pp ∂wj := p|wj |p−1 ∂|wj |∂wj = p|wj | p−2wj , we have that
∇h(w) = x̄+ 10abs(w)8 ◦ w, (132) where ◦ denotes element-wise product and abs(·) denotes the absolution value of its argument in the element-wise way.
By initialization w0 = w−1 = 0, we have that (w>0 u1) = (w > −1u1) = (w > 0 u2) = (w > −1u2) = 0 and the dynamics
w>t+1u2 = (1− η)w>t u2 + β(w>t u2 − w>t−1u2)− ηu>2 ∇h(wt) w>t+1u1 = (1 + η
10 )w>t u1 + β(w > t u1 − w>t−1u1)− ηu>1 ∇h(wt),
(133)
while we also have the initial condition that w>1 u2 = −ηu>2 ∇h(w0) = −ηu>2 ∇h( [ 0 0 ] ) = −ηu>2 x̄ = −ηx̄[1]
w>1 u1 = −ηu>1 ∇h(w0) = −ηu>1 ∇h( [ 0 0 ] ) = −ηu>1 x̄ = −ηx̄[2].
(134)
That is, w1 = −ηx̄. (135)
Using (132), we can rewrite (133) as
w>t+1u2 = (1− η)w>t u2 + β(w>t u2 − w>t−1u2)− ηx̄[1]− 10η(w>t+1u2)9 w>t+1u1 = (1 + η
10 )w>t u1 + β(w > t u1 − w>t−1u1)− ηx̄[2]− 10η(w>t+1u1)9.
(136)
Note that we have that ∇f(w) := Hw +∇h(w) = Hw + x̄+ 10abs(w)8 ◦ w
=
[ 1 + 10|w[1]|8 0 0 −0.1 + 10|w[2]|8 ] w + x̄
(137)
So the stationary points of the objective function satisfy
(1 + 10|w[1]|8)w[1] + x̄[1] = 0 (−0.1 + 10|w[2]|8)w[2] + x̄[2] = 0
(138)
To check if the stationary point is a local minimum, we can use the expression of the Hessian
∇2f(w) := H +∇2h(w)
=
[ 1 + 90|w[1]|8 0 0 −0.1 + 90|w[2]|8 ] .
(139)
Therefore, the Hessian at a stationary point is positive semi-definite as long as |w[2]|8 ≥ 1900 , which can be guaranteed by some realizations of x̄[2] according to (138).
Now let us illustrate why higher momentum β leads to the faster convergence in a high level way, From (138), we see that one can specify the stationary point by determining x̄[1] and x̄[2]. Furthermore, from (139), once the iterate wt satisfies |wt[2]|8 > 1900 , the iterate enters the locally strongly convex and smooth region, for which the local convergence of the Heavy Ball algorithm is known (Ghadimi et al. (2015)). W.l.o.g, let us assume that x̄[2] is negative. From (134), we have that w>1 u1 > 0 when x̄[2] is negative. Moreover, from (136), if, before the iterate satisfies |wt[2]|8 > 1900 , we have that 1 10 (w > t u1) − 10(w>t u1)9 − x̄[2] > 0, then the contribution of the projection on the escape direction (i.e. w>t+1u1) due to the momentum term β(w > t u1 − w>t−1u1) is positive, which also implies that the larger β, the larger the contribution and hence the faster the growing rate of |wt[2]|. Now let us check the condition, 110 (w > t u1)− 10(w>t u1)9− x̄[2] > 0 before |wt[2]|8 > 1900 . A sufficient condition is that ( 1 10 − 1 90 )(w > t u1)− x̄[2] > 0, which we immediately see that it is true given that x̄[2] is negative and that w>1 u1 > 0 and that the magnitude of w > t u1 is increasing. A similar reasoning can be conducted on the other coordinate wt+1[1] := w>t+1u2.

heading: H EMPIRICAL RESULTS

heading: H.1 PHASE RETRIEVAL
text: All the lines are obtained by initializing the iterate at the same point w0 ∼ N (0, Id/(10000d)) and using the same step size η = 5 × 10−4. Here we set w∗ = e1 and sample xi ∼ N (0, Id) with dimension d = 10 and number of samples n = 200. We see that the higher the momentum parameter β, the faster the algorithm enters the linear convergence regime. For the line represented by HB+(β = 1.0 → 0.9), it means switching to the use of β = 0.9 from β = 1 after some iterations. Below, Algorithm 3 and Algorithm 4, we show two equivalent presentations of this
practice. In our experiment, for the ease of implementation, we let the criteria of the switch be 1{ f(w1)−f(wt)f(w1) ≥ 0.5}, i.e. if the relative change of objective value compared to the initial value has been increased to 50%.
Algorithm 3: Switching β = 1 to β < 1 1: Required: step size η and momentum parameter β ∈ [0, 1). 2: Init: u = v = w0 ∈ Rd and β̂ = 1. 3: for t = 0 to T do 4: Update iterate wt+1 = wt − η∇f(wt) + β̂(u− v). 5: Update auxiliary iterate u = v 6: Update auxiliary iterate v = wt. 7: If { Criteria is met } 8: β̂ = β. 9: u = v = wt+1 # (reset momentum) 10: end 11: end for
Algorithm 4: Switching β = 1 to β < 1. 1: Required: step size η and momentum parameter β ∈ [0, 1). 2: Init: w0 ∈ Rd, m−1 = 0d, β̂ = 1. 3: for t = 0 to T do 4: Update momentum mt := β̂mt−1 +∇f(wt). 5: Update iterate wt+1 := wt − ηmt. 6: If { Criteria is met } 7: β̂ = β. 8: mt = 0. # (reset momentum) 9: end 10: end for

heading: H.2 CUBIC-REGULARIZED PROBLEM
text: Figure 2 shows empirical results of solving the cubic-regularized problem by Heavy Ball with different values of momentum parameter β. Subfigure (a) shows that larger momentum parameter β results in a faster growth rate of ‖wt‖, which confirms Lemma 2 and shows that it enters the benign region B faster with larger β. Note that here we have that ‖w∗‖ = 1. It suggests that the norm is non-decreasing during the execution of the algorithm for a wide range of β except very large β For β = 0.9, the norm starts decreasing only after it arises above ‖w∗‖. Subfigure (b) show that higher β also accelerates the linear convergence, as one can see that the slope of a line that corresponds to a higher β is steeper than that of the lower one (e.g. compared to β = 0), which verifies Theorem 3. We also observe a very interesting phenomenon: when β is set to a very large value (e.g. 0.9 here), the pattern is intrinsically different from the smaller ones. The convergence is not monotone and its behavior (bump and overshoots when decreasing); furthermore, the norm of ‖wt‖ generated by the high β is larger than ‖w∗‖ of the minimizer at some time during the execution of the algorithm, which is different from the behavior due to using smaller values of β (i.e. non-decreasing of the norm until the convergence). Our theoretical results cannot explain the behavior of such high β, as such value of β exceeds the upper-threshold required by the theorem. An investigation and understanding of the observation might be needed in the future.
Now let us switch to describe the setup of the experiment. We first set step size η = 0.01, dimension d = 4, ρ = ‖w∗‖ = ‖A‖2 = 1, γ = 0.2 and gap = 5 × 10−3. Then we set A = diag([−γ;−γ + gap; a33; a44]), where the entries a33 and a44 are sampled uniformly random in [−γ + gap; ‖A‖2]. We draw w̃ = (A + ρ‖w∗‖Id)−ξθ, where θ ∼ N (0; Id) and log2 ξ is uniform on [−1, 1]. We set w∗ = ‖w∗‖ ‖w̃‖ w̃ and b = −(A + ρ‖w∗‖Id)w∗. The procedure makes w∗ the global minimizer of problem instance (A, b, ρ). Patterns shown on this figure exhibit for other random problem instances as well.

Review summaries:
The authors analysed the dynamics heavy ball momentum in non-convex optimization settings (namely the phase retrieval and the cubic-regularized problems). The results show that the algorithm leads to a faster convergence rate in comparison with gradient descent without inertial contributions. Although there are results in the context of convex optimization, the result in non-convex setting is novel and very interesting. 

This paper analyzes Polyak momentum in the deterministic case for two simple but important non-convex problems: phase retrieval and finding the cubic-regularized newton step. It is shown that in both cases the problems posses a “benign region” in which the objective “looks” a bit convex. In the phase retrieval case, assuming that the data is generated from a normal distribution, it is shown that having a non-zero momentum parameter in fact allows the momentum method to find the benign region faster, after which standard convex analysis can take over to find the optimal point.