Escape saddle points by a simple gradient-descent based algorithm
Year: 2021

Abstract:
Escaping saddle points is a central research topic in nonconvex optimization. In this paper, we propose a simple gradient-based algorithm such that for a smooth function f : Rn ! R, it outputs an ✏-approximate second-order stationary point in Õ(log n/✏1.75) iterations. Compared to the previous state-of-the-art algorithms by Jin et al. with Õ(log n/✏2) or Õ(log n/✏1.75) iterations, our algorithm is polynomially better in terms of log n and matches their complexities in terms of 1/✏. For the stochastic setting, our algorithm outputs an ✏-approximate second-order stationary point in Õ(log n/✏4) iterations. Technically, our main contribution is an idea of implementing a robust Hessian power method using only gradients, which can find negative curvature near saddle points and achieve the polynomial speedup in log n compared to the perturbed gradient descent methods. Finally, we also perform numerical experiments that support our results.

heading: 1 Introduction
text: Nonconvex optimization is a central research area in optimization theory, since lots of modern machine learning problems can be formulated in models with nonconvex loss functions, including deep neural networks, principal component analysis, tensor decomposition, etc. In general, finding a global minimum of a nonconvex function is NP-hard in the worst case. Instead, many theoretical works focus on finding a local minimum instead of a global one, because recent works (both empirical and theoretical) suggested that local minima are nearly as good as global minima for a significant amount of well-studied machine learning problems; see e.g. [4, 12, 14, 15, 17, 18]. On the other hand, saddle points are major obstacles for solving these problems, not only because they are ubiquitous in high-dimensional settings where the directions for escaping may be few (see e.g. [5, 8, 11]), but also saddle points can correspond to highly suboptimal solutions (see e.g. [19, 28]).
Hence, one of the most important topics in nonconvex optimization is to escape saddle points. Specifically, we consider a twice-differentiable function f : Rn ! R such that
• f is `-smooth: krf(x1) rf(x2)k  `kx1 x2k 8x1,x2 2 Rn, • f is ⇢-Hessian Lipschitz: kH(x1) H(x2)k  ⇢kx1 x2k 8x1,x2 2 Rn;
here H is the Hessian of f . The goal is to find an ✏-approximate second-order stationary point x✏:1
krf(x✏)k  ✏, min(H(x✏)) p ⇢✏. (1)
⇤Corresponding author. Email: tongyangli@pku.edu.cn 1We can ask for an (✏1, ✏2)-approx. second-order stationary point s.t. krf(x)k  ✏1 and min(r2f(x))
✏2 in general. The scaling in (1) was adopted as a standard in literature [1, 6, 10, 20, 21, 22, 26, 29, 30, 31].
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In other words, at any ✏-approx. second-order stationary point x✏, the gradient is small with norm being at most ✏ and the Hessian is close to be positive semi-definite with all its eigenvalues p⇢✏.
Algorithms for escaping saddle points are mainly evaluated from two aspects. On the one hand, considering the enormous dimensions of machine learning models in practice, dimension-free or almost dimension-free (i.e., having poly(log n) dependence) algorithms are highly preferred. On the other hand, recent empirical discoveries in machine learning suggests that it is often feasible to tackle difficult real-world problems using simple algorithms, which can be implemented and maintained more easily in practice. On the contrary, algorithms with nested loops often suffer from significant overheads in large scales, or introduce concerns with the setting of hyperparameters and numerical stability (see e.g. [1, 7]), making them relatively hard to find practical implementations.
It is then natural to explore simple gradient-based algorithms for escaping from saddle points. The reason we do not assume access to Hessians is because its construction takes ⌦(n2) cost in general, which is computationally infeasible when the dimension is large. A seminal work along this line was by Ge et al. [12], which found an ✏-approximate second-order stationary point satisfying (1) using only gradients in O(poly(n, 1/✏)) iterations. This is later improved to be almost dimension-free Õ(log4 n/✏2) in the follow-up work [20],2 and the perturbed accelerated gradient descent algorithm [22] based on Nesterov’s accelerated gradient descent [27] takes Õ(log6 n/✏1.75) iterations. However, these results still suffer from a significant overhead in terms of log n. On the other direction, Refs. [3, 25, 30] demonstrate that an ✏-approximate second-order stationary point can be find using gradients in Õ(log n/✏1.75) iterations. Their results are based on previous works [1, 6] using Hessian-vector products and the observation that the Hessian-vector product can be approximated via the difference of two gradient queries. Hence, their implementations contain nested-loop structures with relatively large numbers of hyperparameters. It has been an open question whether it is possible to keep both the merits of using only first-order information as well as being close to dimension-free using a simple, gradient-based algorithm without a nested-loop structure [23]. This paper answers this question in the affirmative.
Contributions. Our main contribution is a simple, single-loop, and robust gradient-based algorithm that can find an ✏-approximate second-order stationary point of a smooth, Hessian Lipschitz function f : Rn ! R. Compared to previous works [3, 25, 30] exploiting the idea of gradient-based Hessian power method, our algorithm has a single-looped, simpler structure and better numerical stability. Compared to the previous state-of-the-art results with single-looped structures by [22] and [20, 21] using Õ(log6 n/✏1.75) or Õ(log4 n/✏2) iterations, our algorithm achieves a polynomial speedup in log n: Theorem 1 (informal). Our single-looped algorithm finds an ✏-approximate second-order stationary point in Õ(log n/✏1.75) iterations.
Technically, our work is inspired by the perturbed gradient descent (PGD) algorithm in [20, 21] and the perturbed accelerated gradient descent (PAGD) algorithm in [22]. Specifically, PGD applies gradient descents iteratively until it reaches a point with small gradient, which can be a potential saddle point. Then PGD generates a uniform perturbation in a small ball centered at that point and then continues the GD procedure. It is demonstrated that, with an appropriate choice of the perturbation radius, PGD can shake the point off from the neighborhood of the saddle point and converge to a second-order stationary point with high probability. The PAGD in [22] adopts a similar perturbation idea, but the GD is replaced by Nesterov’s AGD [27].
Our algorithm is built upon PGD and PAGD but with one main modification regarding the perturbation idea: it is more efficient to add a perturbation in the negative curvature direction nearby the saddle point, rather than the uniform perturbation in PGD and PAGD, which is a compromise since we generally cannot access the Hessian at the saddle due to its high computational cost. Our key observation lies in the fact that we do not have to compute the entire Hessian to detect the negative curvature. Instead, in a small neighborhood of a saddle point, gradients can be viewed as Hessianvector products plus some bounded deviation. In particular, GD near the saddle with learning rate 1/` is approximately the same as the power method of the matrix (I H/`). As a result, the most negative eigenvalues stand out in GD because they have leading exponents in the power method, and thus it approximately moves along the direction of the most negative curvature nearby the sad-
2The Õ notation omits poly-logarithmic terms, i.e., Õ(g) = O(g poly(log g)).
dle point. Following this approach, we can escape the saddle points more rapidly than previous algorithms: for a constant ✏, PGD and PAGD take O(log n) iterations to decrease the function value by ⌦(1/ log3 n) and ⌦(1/ log5 n) with high probability, respectively; on the contrary, we can first take O(log n) iterations to specify a negative curvature direction, and then add a larger perturbation in this direction to decrease the function value by ⌦(1). See Proposition 3 and Proposition 5. After escaping the saddle point, similar to PGD and PAGD, we switch back to GD and AGD iterations, which are efficient to decrease the function value when the gradient is large [20, 21, 22].
Our algorithm is also applicable to the stochastic setting where we can only access stochastic gradients, and the stochasticity is not under the control of our algorithm. We further assume that the stochastic gradients are Lipschitz (or equivalently, the underlying functions are gradient-Lipschitz, see Assumption 2), which is also adopted in most of the existing works; see e.g. [9, 20, 21, 35]. We demonstrate that a simple extended version of our algorithm takes O(log2 n) iterations to detect a negative curvature direction using only stochastic gradients, and then obtain an ⌦(1) function value decrease with high probability. On the contrary, the perturbed stochastic gradient descent (PSGD) algorithm in [20, 21], the stochastic version of PGD, takes O(log10 n) iterations to decrease the function value by ⌦(1/ log5 n) with high probability. Theorem 2 (informal). In the stochastic setting, our algorithm finds an ✏-approximate second-order stationary point using Õ(log2 n/✏4) iterations via stochastic gradients.
Our results are summarized in Table 1. Although the underlying dynamics in [3, 25, 30] and our algorithm have similarity, the main focus of our work is different. Specifically, Refs. [3, 25, 30] mainly aim at using novel techniques to reduce the iteration complexity for finding a second-order stationary point, whereas our work mainly focuses on reducing the number of loops and hyper-parameters of negative curvature finding methods while preserving their advantage in iteration complexity, since a much simpler structure accords with empirical observations and enables wider applications. Moreover, the choice of perturbation in [3] is based on the Chebyshev approximation theory, which may require additional nested-looped structures to boost the success probability. In the stochastic setting, there are also other results studying nonconvex optimization [16, 24, 32, 37, 13, 33, 36] from different perspectives than escaping saddle points, which are incomparable to our results.
It is worth highlighting that our gradient-descent based algorithm enjoys the following nice features:
• Simplicity: Some of the previous algorithms have nested-loop structures with the concern of practical impact when setting the hyperparameters. In contrast, our algorithm based on negative curvature finding only contains a single loop with two components: gradient descent (including AGD or SGD) and perturbation. As mentioned above, such simple structure is preferred in machine learning, which increases the possibility of our algorithm to find real-world applications.
• Numerical stability: Our algorithm contains an additional renormalization step at each iteration when escaping from saddle points. Although in theoretical aspect a renormalization step does not affect the output and the complexity of our algorithm, when finding negative curvature near
saddle points it enables us to sample gradients in a larger region, which makes our algorithm more numerically stable against floating point error and other errors. The introduction of renormalization step is enabled by the simple structure of our algorithm, which may not be feasible for more complicated algorithms [3, 25, 30].
• Robustness: Our algorithm is robust against adversarial attacks when evaluating the value of the gradient. Specifically, when analyzing the performance of our algorithm near saddle points, we essentially view the deflation from pure quadratic geometry as an external noise. Hence, the effectiveness of our algorithm is unaffected under external attacks as long as the adversary is bounded by deflations from quadratic landscape.
Finally, we perform numerical experiments that support our polynomial speedup in log n. We perform our negative curvature finding algorithms using GD or SGD in various landscapes and general classes of nonconvex functions, and use comparative studies to show that our Algorithm 1 and Algorithm 3 achieve a higher probability of escaping saddle points using much fewer iterations than PGD and PSGD (typically less than 1/3 times of the iteration number of PGD and 1/2 times of the iteration number of PSGD, respectively). Moreover, we perform numerical experiments benchmarking the solution quality and iteration complexity of our algorithm against accelerated methods. Compared to PAGD [22] and even advanced optimization algorithms such as NEON+ [30], Algorithm 2 possesses better solution quality and iteration complexity in various landscapes given by more general nonconvex functions. With fewer iterations compared to PAGD and NEON+ (typically less than 1/3 times of the iteration number of PAGD and 1/2 times of the iteration number of NEON+, respectively), our Algorithm 2 achieves a higher probability of escaping from saddle points.
Open questions. This work leaves a couple of natural open questions for future investigation:
• Can we achieve the polynomial speedup in log n for more advanced stochastic optimization algorithms with complexity Õ(poly(log n)/✏3.5) [2, 3, 10, 29, 31] or Õ(poly(log n)/✏3) [9, 34]?
• How is the performance of our algorithms for escaping saddle points in real-world applications, such as tensor decomposition [12, 17], matrix completion [14], etc.?
Broader impact. This work focuses on the theory of nonconvex optimization, and as far as we see, we do not anticipate its potential negative societal impact. Nevertheless, it might have a positive impact for researchers who are interested in understanding the theoretical underpinnings of (stochastic) gradient descent methods for machine learning applications.
Organization. In Section 2, we introduce our gradient-based Hessian power method algorithm for negative curvature finding, and present how our algorithms provide polynomial speedup in log n for both PGD and PAGD. In Section 3, we present the stochastic version of our negative curvature finding algorithm using stochastic gradients and demonstrate its polynomial speedup in log n for PSGD. Numerical experiments are presented in Section 4. We provide detailed proofs and additional numerical experiments in the supplementary material.

heading: 2 A simple algorithm for negative curvature finding
text: We show how to find negative curvature near a saddle point using a gradient-based Hessian power method algorithm, and extend it to a version with faster convergence rate by replacing gradient descents by accelerated gradient descents. The intuition works as follows: in a small enough region nearby a saddle point, the gradient can be approximately expressed as a Hessian-vector product formula, and the approximation error can be efficiently upper bounded, see Eq. (6). Hence, using only gradients information, we can implement an accurate enough Hessian power method to find negative eigenvectors of the Hessian matrix, and further find the negative curvature nearby the saddle.

heading: 2.1 Negative curvature finding based on gradient descents
text: We first present an algorithm for negative curvature finding based on gradient descents. Specifically, for any x̃ 2 Rn with min(H(x̃))  p ⇢✏, it finds a unit vector ê such that êTH(x̃)ê  p⇢✏/4.
Algorithm 1: Negative Curvature Finding(x̃, r,T ). 1 y0 Uniform(Bx̃(r)) where Bx̃(r) is the `2-norm ball centered at x̃ with radius r; 2 for t = 1, ...,T do 3 yt yt 1 kyt 1k `r rf(x̃+ ryt 1/kyt 1k) rf(x̃) ;
4 Output yT /r.
Proposition 3. Suppose the function f : Rn ! R is `-smooth and ⇢-Hessian Lipschitz. For any 0 < 0  1, we specify our choice of parameters and constants we use as follows:
T = 8` p ⇢✏
· log ⇣ `
0
r n
⇡⇢✏
⌘ , r = ✏
8`
r ⇡
n 0. (2)
Suppose that x̃ satisfies min(r2f(x̃))  p ⇢✏. Then with probability at least 1 0, Algorithm 1 outputs a unit vector ê satisfying
êTH(x)ê  p ⇢✏/4, (3)
using O(T ) = Õ ⇣
lognp ⇢✏
⌘ iterations, where H stands for the Hessian matrix of function f .
Proof. Without loss of generality we assume x̃ = 0 by shifting Rn such that x̃ is mapped to 0. Define a new n-dimensional function
hf (x) := f(x) hrf(0),xi , (4)
for the ease of our analysis. Since hrf(0),xi is a linear function with Hessian being 0, the Hessian of hf equals to the Hessian of f , and hf (x) is also `-smooth and ⇢-Hessian Lipschitz. In addition, note that rhf (0) = rf(0) rf(0) = 0. Then for all x 2 Rn,
rhf (x) =
Z 1
⇠=0 H(⇠x) · x d⇠ = H(0)x+
Z 1
⇠=0 (H(⇠x) H(0)) · x d⇠. (5)
Furthermore, due to the ⇢-Hessian Lipschitz condition of both f and hf , for any ⇠ 2 [0, 1] we have kH(⇠x) H(0)k  ⇢kxk, which leads to
krhf (x) H(0)xk  ⇢kxk 2 . (6)
Observe that the Hessian matrix H(0) admits the following eigen-decomposition:
H(0) = nX
i=1
iuiu T i , (7)
where the set {ui}ni=1 forms an orthonormal basis of Rn. Without loss of generality, we assume the eigenvalues 1, 2, . . . , n corresponding to u1,u2, . . . ,un satisfy
1  2  · · ·  n, (8)
in which 1  p ⇢✏. If n  p ⇢✏/2, Proposition 3 holds directly. Hence, we only need to prove the case where n > p ⇢✏/2, in which there exists some p, p0 with
p  p ⇢✏ < p+1, p0  p ⇢✏/2 < p0+1. (9)
We use Sk, S? to separately denote the subspace of Rn spanned by {u1,u2, . . . ,up}, {up+1,up+2, . . . ,un}, and use S0k, S 0 ? to denote the subspace of Rn spanned by {u1,u2, . . . ,up0}, {up0+1,up+2, . . . ,un}. Furthermore, we define yt,k := Pp
i=1 hui,ytiui, yt,? := Pn i=p hui,ytiui, yt,k0 := Pp0 i=1 hui,ytiui, yt,?0 := Pn i=p0 hui,ytiui respectively to denote the component of yt in Line 3 in the subspaces Sk, S?, S0k, S 0 ?, and let ↵t := kyt,kk/kytk. Observe that
Pr n ↵0 0 p ⇡/n o Pr n |y0,1|/r 0 p ⇡/n o , (10)
where y0,1 := hu1,y0i denotes the component of y0 along u1. Consider the case where ↵0 0 p ⇡/n, which can be achieved with probability
Pr ⇢ ↵0 r ⇡
n 0
1
r ⇡
n 0 · Vol(Bn 10 (1)) Vol(Bn0 (1)) 1
r ⇡
n 0 ·
r n
⇡ = 1 0. (11)
We prove that there exists some t0 with 1  t0  T such that kyt0,?0k/kyt0k  p ⇢✏/(8`). (12) Assume the contrary, for any 1  t  T , we all have kyt,?0k/kytk > p ⇢✏/(8`). Then ky0t,?k satisfies the following recurrence formula: kyt+1,?0k  (1 + p ⇢✏/(2`))kyt,?0k+ k ?0k  (1 + p ⇢✏/(2`) + k k/kyt,?0k)kyt,?0k, (13)
where := kytkr` (rhf (ryt/kytk) H(0) · (ryt/kytk)) and k k/kytk  ⇢r/` due to (6). Hence,
kyt+1,?0k  ⇣ 1 + p ⇢✏
2` + k k
kyt,?0k
⌘ kyt,?0k  ⇣ 1 + p ⇢✏
2` + · 8⇢r p ⇢✏
⌘ kyt+1,?0k, (14)
which leads to kyt,?0k  ky0,?0k(1 + p ⇢✏/(2`) + 8⇢r/ p ⇢✏)t  ky0,?0k(1 + 5 p ⇢✏/(8`))t, 8t 2 [T ]. (15)
Similarly, we can have the recurrence formula for kyt,kk:
kyt+1,kk (1 + p ⇢✏/(2`))kyt,kk k kk (1 + p ⇢✏/(2`) k k/(↵tkytk))kyt,kk. (16)
Considering that k k/kytk  ⇢r/` due to (6), we can further have kyt+1,kk (1 + p ⇢✏/(2`) ⇢r/(↵t`))kyt,kk. (17)
Intuitively, kyt,kk should have a faster increasing rate than kyt,?k in this gradient-based Hessian power method, ignoring the deviation from quadratic approximation. As a result, the value the value ↵t = kyt,kk/kytk should be non-decreasing. It is demonstrated in Lemma 17 in Appendix B that, even if we count this deviation in, ↵t can still be lower bounded by some constant ↵min:
↵t ↵min = 0
4
r ⇡
n , 81  t  T . (18)
by which we can further deduce that kyt,kk ky0,kk(1 + p ⇢✏/` ⇢r/(↵min`)) t ky0,kk(1 + 7 p ⇢✏/(8`))t, 81  t  T . (19)
Observe that kyT ,?0k
kyT ,kk  ky0,?0k ky0,kk · ⇣1 + 5p⇢✏/(8`) 1 + 7 p ⇢✏/(8`) ⌘T  1 0
r n
⇡ ⇣1 + 5p⇢✏/(8`) 1 + 7 p ⇢✏/(8`) ⌘T  p ⇢✏ 8` . (20)
Since kyT ,kk  kyT k, we have kyT ,?0k/kyT k  p ⇢✏/(8`), contradiction. Hence, there here exists some t0 with 1  t0  T such that kyt0,?0k/kyt0k  p ⇢✏/(8`). Consider the normalized vector ê = yt0/r, we use ê?0 and êk0 to separately denote the component of ê in S0? and S 0 k. Then, kê?0k  p ⇢✏/(8`) whereas kêk0k 1 ⇢✏/(8`)2. Then,
êTH(0)ê = (ê?0 + êk0) T H(0)(ê?0 + êk0) = ê T ?0H(0)ê?0 + ê T k0H(0)êk0 (21)
since H(0)ê?0 2 S0? and H(0)êk0 2 S 0 k. Due to the `-smoothness of the function, all eigenvalue of the Hessian matrix has its absolute value upper bounded by `. Hence,
êT?0H(0)ê?0  `kê T ?0k 2 2 = ⇢✏/(64` 2). (22) Further according to the definition of Sk, we have
êTk0H(0)êk0  p ⇢✏kêk0k 2 /2. (23)
Combining these two inequalities together, we can obtain
êTH(0)ê = êT?H(0)ê?0 + ê T k0H(0)êk0  p ⇢✏kêk0k 2 /2 + ⇢✏/(64`2)  p ⇢✏/4. (24)
Remark 4. In practice, the value of kytk can become large during the execution of Algorithm 1. To fix this, we can renormalize yt to have `2-norm r at the ends of such iterations, and this does not influence the performance of the algorithm.

heading: 2.2 Faster negative curvature finding based on accelerated gradient descents
text: In this subsection, we replace the GD part in Algorithm 1 by AGD to obtain an accelerated negative curvature finding subroutine with similar effect and faster convergence rate, based on which we further implement our Accelerated Gradient Descent with Negative Curvature Finding (Algorithm 2). Near any saddle point x̃ 2 Rn with min(H(x̃))  p ⇢✏, Algorithm 2 finds a unit vector ê such that êTH(x̃)ê  p⇢✏/4.
Algorithm 2: Perturbed Accelerated Gradient Descent with Accelerated Negative Curvature Finding(x0, ⌘, ✓, , s,T 0, r0)
1 tperturb 0, z0 x0, x̃ x0, ⇣ 0; 2 for t = 0, 1, 2, ..., T do 3 if krf(xt)k  ✏ and t tperturb > T then 4 x̃ = xt; 5 xt Uniform(Bx̃(r0)) where Uniform(Bx̃(r0)) is the `2-norm ball centered at x̃ with radius r0, zt xt, ⇣ rf(x̃), tperturb t; 6 if t tperturb = T 0 then 7 ê := xt x̃kxt x̃k ; 8 xt x̃ f 0ê(x̃)
4|f 0ê(x̃)|
q ✏ ⇢ · ê, zt xt, ⇣ = 0;
9 xt+1 zt ⌘(rf(zt) ⇣); 10 vt+1 xt+1 xt; 11 zt+1 xt+1 + (1 ✓)vt+1; 12 if tperturb 6= 0 and t tperturb < T 0 then 13 zt+1 x̃+ r0 ·
zt+1 x̃ kzt+1 x̃k , xt+1 x̃+ r 0 · xt+1 x̃ kzt+1 x̃k ;
14 else
15 if f(xt+1)  f(zt+1) + hrf(zt+1),xt+1 zt+1i 2 kzt+1 xt+1k
2 then
16 (xt+1,vt+1) NegativeCurvatureExploitation(xt+1,vt+1, s)3; 17 zt+1 xt+1 + (1 ✓)vt+1;
The following proposition exhibits the effectiveness of Algorithm 2 for finding negative curvatures near saddle points: Proposition 5. Suppose the function f : Rn ! R is `-smooth and ⇢-Hessian Lipschitz. For any 0 < 0  1, we specify our choice of parameters and constants we use as follows:
⌘ := 1
4` ✓ :=
(⇢✏)1/4
4 p `
T 0 := 32 p `
(⇢✏)1/4 log
⇣ `
0
r n
⇢✏
⌘
:= ✓ 2
⌘ s := 4⇢ r 0 :=
0✏
32
r ⇡
⇢n (25)
then for a point x̃ satisfying min(r2f(x̃))  p ⇢✏, if running Algorithm 2 with the uniform perturbation in Line 5 being added at t = 0, the unit vector ê in Line 7 obtained after T 0 iterations satisfies:
P ⇣ êTH(x)ê  p ⇢✏/4 ⌘ 1 0. (26)
The proof of Proposition 5 is similar to the proof of Proposition 3, and is deferred to Appendix B.2.

heading: 2.3 Escaping saddle points using negative curvature finding
text: In this subsection, we demonstrate that our Algorithm 1 and Algorithm 2 with the ability to find negative curvature near saddle points can further escape saddle points of nonconvex functions. The
3 This NegativeCurvatureExploitation (NCE) subroutine was originally introduced in [22, Algorithm 3] and is called when we detect that the current momentum vt coincides with a negative curvature direction of zt. In this case, we reset the momentum vt and decide whether to exploit this direction based on the value of kvtk.
intuition works as follows: we start with gradient descents or accelerated gradient descents until the gradient becomes small. At this position, we compute the negative curvature direction, described by a unit vector ê, via Algorithm 1 or the negative curvature finding subroutine of Algorithm 2. Then, we add a perturbation along this direction of negative curvature and go back to gradient descents or accelerated gradient descents with an additional NegativeCurvatureExploitation subroutine (see Footnote 3). It has the following guarantee: Lemma 6. Suppose the function f : Rn ! R is `-smooth and ⇢-Hessian Lipschitz. Then for any point x0 2 Rn, if there exists a unit vector ê satisfying êTH(x0)ê  p ⇢✏ 4 where H stands for the Hessian matrix of function f , the following inequality holds:
f ⇣ x0
f 0 ê(x0)
4|f 0ê(x0)|
r ✏ ⇢ · ê ⌘  f(x0)
1
384
s ✏3
⇢ , (27)
where f 0ê stands for the gradient component of f along the direction of ê.
Proof. Without loss of generality, we assume x0 = 0. We can also assume hrf(0), êi  0; if this is not the case we can pick ê instead, which still satisfies ( ê)TH(x0)( ê)  p ⇢✏ 4 . In practice, to figure out whether we should use ê or ê, we apply both of them in (27) and choose the one with smaller function value. Then, for any x = xêê with some xê > 0, we have @
2f @x2ê (x)  p ⇢✏ 4 + ⇢xê
due to the ⇢-Hessian Lipschitz condition of f . Hence, @f
@xê (x)  f 0ê(0)
p ⇢✏
4 xê + ⇢x
2 ê, (28)
by which we can further derive that
f(xêê) f(0)  f 0 ê(0)xê
p ⇢✏
8 x 2 ê +
⇢ 3 x 3 ê 
p ⇢✏
8 x 2 ê +
⇢ 3 x 3 ê. (29)
Settings xê = 14 q ✏ ⇢ gives (27).
We give the full algorithm details based on Algorithm 1 in Appendix C.1. Along this approach, we achieve the following: Theorem 7 (informal, full version deferred to Appendix C.3). For any ✏ > 0 and a constant 0 <  1, Algorithm 2 satisfies that at least one of the iterations xt will be an ✏-approximate secondorder stationary point in
Õ ⇣ (f(x0) f⇤) ✏1.75 · log n ⌘
(30)
iterations, with probability at least 1 , where f⇤ is the global minimum of f .
Intuitively, the proof of Theorem 7 has two parts. The first part is similar to the proof of [22, Theorem 3], which shows that PAGD uses Õ(log6 n/✏1.75) iterations to escape saddle points. We show that there can be at most Õ( f/✏1.75) iterations with the norm of gradient larger than ✏ using almost the same techniques, but with slightly different parameter choices. The second part is based on the negative curvature part of Algorithm 2, our accelerated negative curvature finding algorithm. Specifically, at each saddle point we encounter, we can take Õ(log n/✏1/4) iterations to find its negative curvature (Proposition 5), and add a perturbation in this direction to decrease the function value by O(✏1.5) (Lemma 6). Hence, the iterations introduced by Algorithm 4 can be at most Õ logn
✏1.5 · 1 ✏0.25 = Õ(log n/✏1.75), which is simply an upper bound on the overall iteration number. The detailed proof is deferred to Appendix C.3. Remark 8. Although Theorem 7 only demonstrates that our algorithms will visit some ✏- approximate second-order stationary point during their execution with high probability, it is straightforward to identify one of them if we add a termination condition: once Negative Curvature Finding (Algorithm 1 or Algorithm 2) is applied, we record the position xt0 and the function value decrease due to the following perturbation. If the function value decrease is larger than 1384 q ✏3
⇢ as per Lemma 6, then the algorithms make progress. Otherwise, xt0 is an ✏-approximate second-order stationary point with high probability.

heading: 3 Stochastic setting
text: In this section, we present a stochastic version of Algorithm 1 using stochastic gradients, and demonstrate that it can also be used to escape saddle points and obtain a polynomial speedup in log n compared to the perturbed stochastic gradient (PSGD) algorithm in [21].

heading: 3.1 Stochastic negative curvature finding
text: In the stochastic gradient descent setting, the exact gradients oracle rf of function f cannot be accessed. Instead, we only have unbiased stochastic gradients g(x; ✓) such that
rf(x) = E✓⇠D[g(x; ✓)] 8x 2 Rn, (31) where D stands for the probability distribution followed by the random variable ✓. Define
⇣(x; ✓) := g(x; ✓) rf(x) (32) to be the error term of the stochastic gradient. Then, the expected value of vector ⇣(x; ✓) at any x 2 Rn equals to 0. Further, we assume the stochastic gradient g(x, ✓) also satisfies the following assumptions, which were also adopted in previous literatures; see e.g. [9, 20, 21, 35]. Assumption 1. For any x 2 Rn, the stochastic gradient g(x; ✓) with ✓ ⇠ D satisfies:
Pr[(kg(x; ✓) rf(x)k t)]  2 exp t 2 /(2 2) , 8t 2 R. (33)
Assumption 2. For any ✓ 2 supp(D), g(x; ✓) is ˜̀-Lipschitz for some constant ˜̀:
kg(x1; ✓) g(x2; ✓)k  ˜̀kx1 x2k, 8x1,x2 2 Rn. (34)
Assumption 2 emerges from the fact that the stochastic gradient g is often obtained as an exact gradient of some smooth function,
g(x; ✓) = rf(x; ✓). (35) In this case, Assumption 2 guarantees that for any ✓ ⇠ D, the spectral norm of the Hessian of f(x; ✓) is upper bounded by ˜̀. Under these two assumptions, we can construct the stochastic version of Algorithm 1, as shown in Algorithm 3.
Algorithm 3: Stochastic Negative Curvature Finding(x0, rs,Ts,m). 1 y0 0, L0 rs; 2 for t = 1, ...,Ts do 3 Sample ✓ (1) , ✓ (2) , · · · , ✓ (m) ⇠ D;
4 g(yt 1) 1 m Pm j=1 g(x0 + yt 1; ✓(j)) g(x0; ✓(j)) ;
5 yt yt 1 1 ` (g(yt 1) + ⇠t/Lt 1), ⇠t ⇠ N
⇣ 0, r
2 s d I
⌘ ;
6 Lt kytk rs Lt 1 and yt yt · rskytk ;
7 Output yT /rs.
Similar to the non-stochastic setting, Algorithm 3 can be used to escape from saddle points and obtain a polynomial speedup in log n compared to PSGD algorithm in [21]. This is quantitatively shown in the following theorem: Theorem 9 (informal, full version deferred to Appendix D.2). For any ✏ > 0 and a constant 0 <  1, our algorithm4 based on Algorithm 3 using only stochastic gradient descent satisfies that at least one of the iterations xt will be an ✏-approximate second-order stationary point in
Õ ⇣ (f(x0) f⇤) ✏4 · log2 n ⌘
(36)
iterations, with probability at least 1 , where f⇤ is the global minimum of f . 4Our algorithm based on Algorithm 3 has similarities to the Neon2online algorithm in [3]. Both algorithms find a second-order stationary point for stochastic optimization in Õ(log2 n/✏4) iterations, and we both apply directed perturbations based on the results of negative curvature finding. Nevertheless, our algorithm enjoys simplicity by only having a single loop, whereas Neon2online has a nested loop for boosting their confidence.

heading: 4 Numerical experiments
text: In this section, we provide numerical results that exhibit the power of our negative curvature finding algorithm for escaping saddle points. More experimental results can be found in Appendix E. All the experiments are performed on MATLAB R2019b on a computer with Six-Core Intel Core i7 processor and 16GB memory, and their codes are given in the supplementary material.
Comparison between Algorithm 1 and PGD. We compare the performance of our Algorithm 1 with the perturbed gradient descent (PGD) algorithm in [21] on a test function f(x1, x2) = 116x 4 1 1 2x 2 1 + 9 8x 2 2 with a saddle point at (0, 0). The advantage of Algorithm 1 is illustrated in Figure 1.
Comparison between Algorithm 3 and PSGD. We compare the performance of our Algorithm 3 with the perturbed stochastic gradient descent (PSGD) algorithm in [21] on a test function f(x1, x2) = (x31 x 3 2)/2 3x1x2 + (x 2 1 + x 2 2) 2 /2. Compared to the landscape of the previous experiment, this function is more deflated from a quadratic field due to the cubic terms. Nevertheless, Algorithm 3 still possesses a notable advantage compared to PSGD as demonstrated in Figure 2.

heading: Acknowledgements
text: We thank Jiaqi Leng for valuable suggestions and generous help on the design of numerical experiments. TL was supported by the NSF grant PHY-1818914 and a Samsung Advanced Institute of Technology Global Research Partnership.

Review summaries:
This paper provides an improved method for computing stationary points of nonconvex functions. The authors give an algorithm which, given an $L$-smooth $\rho$-Hessian Lipschitz nonconvex function $f$, computes an $(O(\epsilon), O(\sqrt{\rho \epsilon}))$-second order stationary point of $f$ using $O(\eps^{-7/4} \log n)$ gradient computations of $f$. This matches the $\epsilon$ dependence  achieved by several previous works in the area (for example [6,7,20]) but  improves the dependence on $\log n$, where $n$ is the dimension, by polynomial factors. The algorithm is also of similar simplicity to the previous state-of-the-art method given in [20], and combines the perturbed AGD framework given there with some clever observations to achieve the result.

In this work, the authors study the escaping saddle points problem for nonconvex optimization. Previously Jin et al. proposed a perturbed accelerated gradient descent (AGD) method which finds $\epsilon$-second-order stationary points in $O(\log^6 n/\epsilon^{1.75})$ number of iterations, where $n$ is the dimension. By utilizing a refined perturbed negative curvature estimator, the authors successfully improves the dependence of dimension from $\log^6 n$ to $\log n$. Experiment results suggest that their algorithm performs well in practice.