arXiv CSopen access

Finding Stationary Points by Comparisons

Finding Stationary Points by Comparisons Helin Wang1,3∗

arXiv:2606.27082v1 [cs.LG] 25 Jun 2026

1

Chenyi Zhang2,∗

Xiwen Tao1,3

Yexin Zhang3,4

Tongyang Li3,4,†

School of Electronics Engineering and Computer Science, Peking University 2 Computer Science Department, Stanford University 3 Center on Frontiers of Computing Studies, Peking University 4 School of Computer Science, Peking University

Abstract We study the problem of finding stationary points of non-convex functions when access to the objective is provided only through a comparison oracle that, given two points, outputs which has the larger function value. For a twice differentiable f : Rn → R with Lipschitz gradient e 2 /ϵ1.5 ) queries. and Hessian, we develop an algorithm that visits an ϵ-stationary point using O(n Our approach uses a subroutine that estimates the normalized Hessian to accuracy δ using e 2 log(1/δ)) queries. We further study this problem with a quantum comparison oracle model O(n where queries can be made in superpositions, and develop the first quantum algorithm that finds 1.5 e an ϵ-stationary point, which takes O(n/ϵ ) queries.

1

Introduction

We study the problem of finding a stationary point of a non-convex function f : Rn → R, which is a point x ∈ Rn satisfying ∥∇f (x)∥ ≤ ϵ, given an initial x(0) ∈ Rn with bounded initial function value gap, i.e., f (x(0) ) − inf x∈Rn f (x) ≤ ∆. Finding a stationary point is a natural objective in non-convex optimization, since computing the global optimum is NP hard in the worst case. Moreover, in problems including tensor decomposition [17], matrix completion [15], and regression with non-convex regularization [28], global optimality can be obtained by finding a second-order stationary point, which is a generalization of stationary points. As a foundational problem in optimization theory, finding a stationary point, also known as critical point computation [1, 2] or making the gradient small [3, 31], has been extensively studied. Given a gradient oracle, gradient descent finds an ϵ-stationary point of f using O(1/ϵ2 ) iterations, provided that f has Lipschitz gradient, i.e., ∥∇f (x) − ∇f (y)∥ ≤ L1 ∥x − y∥ for any x, y ∈ Rn [32]. More generally, if f has Lp -Lipschitz p-th order derivative, Birgin et al. [7] gives an algorithm using 1/p O(Lp ∆ϵ−(p+1)/p ) queries to a p-th order oracle, where a query at x returns all derivatives of f at x up to order p. Furthermore, Carmon et al. [9] establishes that these rates are optimal among dimension-independent algorithms. More recently, training of machine learning models solicits for even simpler information. For example, it is known that taking only signs of gradient descents still demonstrates good performance in training neural networks [6, 24, 26]. Moreover, in the breakthrough of large language models (LLMs), reinforcement learning from human feedback (RLHF) [11, 13, 14] played an important ∗ †

Equal contribution. Corresponding author. Email: [email protected]

1

role in training these LLMs, for instance on GPT-3 [36]. Compared to standard RL that applies function evaluation for rewards, RLHF is preference-based RL that only compares between options and determines which is better. There is emerging interest in preference-based RL, with a series of results [10, 35, 41, 43, 47, 49] establishing provable guarantees for learning a near-optimal policy from preference feedback. Furthermore, Wang et al. [45] proved that preference-based RL can be solved with small or no extra costs compared to those of standard reward-based RL for a wide range of models. More broadly, these developments reflect a longstanding interest in designing optimization methods using limited feedback. There has been a line of research on solving optimization problems using a comparison oracle; see the survey [23]. Formally, a function f : Rn → R is accessed through a comparison oracle Ofcomp : Rn × Rn → {−1, 1} that upon a pair of inputs (x, y) ∈ Rn × Rn , we get output: ( 1, f (x) ≥ f (y) Ofcomp (x, y) = (1) −1. f (x) ≤ f (y) with either output allowed when f (x) = f (y). Classical direct-search and pattern-search methods use such comparisons to accept or reject trial steps [4, 22], while the Nelder–Mead method [30] relies on comparisons among candidate points to drive the search. However, Nelder–Mead may converge to a nonstationary point even on continuously differentiable convex examples [29]. Based on that, Jamieson et al. [20] studied derivative-free optimization with Boolean comparison feedback and comparison-based line searches. GradientLess Descent [18] studies monotone-invariant zeroth-order optimization, and [43] studied ranking-based feedback, where the oracle provides ranking information over multiple candidate points. Beyond noiseless pairwise comparisons in Euclidean space, Saha et al. [38] studied batched and multiway preference feedback, and Ren et al. [37] considered comparisonbased optimization on Riemannian manifolds. Very recently, Ref. [42] proposed a function-free optimization framework in which the preference relation itself, rather than an underlying scalar objective, defines the optimization problem. For convex objectives, several works give provable guarantees under comparison or order-type feedback. Karabag et al. [21] proposed an ellipsoid-based method for smooth convex optimization using comparison oracles. Bergou et al. [5] and Gorbunov et al. [19] developed the stochastic threepoint (STP) methods using comparisons among randomly sampled trial points to select the next iterate. Dueling optimization [39, 40] studies convex optimization from noisy pairwise comparisons, and Lobanov et al. [27] developed accelerated methods in a similar setting. For nonconvex objectives, fewer comparison-query guarantees are known. The stochastic threepoint method of [5, 19] samples random search directions and chooses the next iterate using comparisons among trial points; for functions with Lipschitz gradients, this gives an O(n/ϵ2 ) comparison-query implementation for finding an ϵ-stationary point in expectation. These results leave open whether comparison access can exploit higher-order smoothness to improve the dependence on ϵ for finding stationary points. e 2 /ϵ1.5 ) queries to Ocomp and Our results. In this paper, we develop an algorithm that makes O(n f guarantees that one of the queried points is an ϵ-stationary point.1 Theorem 1 (Informal). Let f : Rn → R have L1 -Lipschitz gradient and L2 -Lipschitz Hessian. Given n x0 ∈ Rn satisfying f (x0 ) − inf x∈R with success probability at least 2/3, Algorithm 5 visits √ (x)2≤ ∆, 1.5 e an ϵ-stationary point using O(∆ L2 n /ϵ ) queries to Ofcomp . 1

e to hide poly-logarithmic factors in n, ϵ−1 , L1 , L2 , and ∆. We use O(·)

2

Up to poly-logarithmic factors, the ϵ-dependence of our rate matches the optimal rate of secondorder methods [9, 34]. Compared to [5, 19], our algorithm achieves an improved dependence on ϵ, at the cost of a worse dependence on the dimension n. It is unclear whether our bound, particularly the n2 dependence, is asymptotically optimal in any non-trivial regime of n and ϵ. Although [9] give relevant lower bounds in the dimension-independent setting, the dimension-dependent complexity of finding stationary points is still not well understood, even in the setting where the algorithm has access to gradient or Hessian information. We view developing matching lower bounds in the comparison setting as an interesting open problem. Note that here we only obtain the guarantee that our algorithm visits an ϵ-stationary point throughout the iterations instead of finding one, similar to prior works [5, 19]. This is due to the fact that in the comparison oracle model the algorithm only observes relative function values. Consequently, it is in general impossible to access the gradient norm or even to test whether a given point x is an ϵ-stationary point. For applications where a single point is needed, a natural heuristic is to return the iterate with the smallest function value, which can be identified by comparisons. While this does not certify stationarity in general, in many structured nonconvex problems near-optimality is closely tied to stationarity or second-order stationarity, such as tensor decomposition [16], matrix completion [15], and regression with nonconvex regularization [28]. Furthermore, we study the problem of finding stationary points in the quantum setting, where we can query a quantum comparison oracle comp Of,q |x⟩|y⟩|b⟩ = |x⟩|y⟩|b ⊕ 1{f (x) > f (y)}⟩

(2)

which performs comparison in quantum superpositions.2 Extending upon our classical algorithm, we obtain a quantum algorithm for finding stationary points using the quantum comparison oracle (2). Theorem 2 (Informal). Let f : Rn → R have L1 -Lipschitz gradient and L2 -Lipschitz Hessian. Given x0 ∈ Rn satisfying f (x√ 0 ) − inf x∈Rn (x) ≤ ∆, there exists a quantum algorithm that visits an comp e ϵ-stationary point using O(∆ L2 n/ϵ1.5 ) queries to Of,q .

Figure 1: The structure of our stationary point finding algorithm. In this figure the arrow represents the implementation relationship. The three purple frames represent the novel techniques that we propose, the yellow frame represents the known technique from [44], and the green frame represents our main theorem.

Hessian-vector product estimation (Section 3) Gradient direction estimation [44]

Hessian-gradient norm ratio estimation (Section 5)

Normalized Hessian estimation (Section 4)

Finding stationary points (Section 6) Techniques. Our approach for establishing Theorem 1 is summarized in Figure 1. Throughout, we use the comparison-based gradient direction estimation algorithm as a basic primitive. We first give an algorithm ComparisonHessVec (Algorithm 1) that, for any given x, y ∈ Rn , estimates the direction ∇2 f (x)·y. This is done by first estimating the directions of ∇f (x), ∇f (x+ry), 2

Quantum notations can be found in Section 2.

3

and ∇f (x − ry) using [44]. We then infer the direction of ∇f (x + ry) − ∇f (x − ry) using a geometric property: its intersection with ∇f (x) and ∇f (x + ry) as well as its intersection with ∇f (x) and ∇f (x − ry) give two segments of same length (see Figure 2). By choosing r sufficiently small, this direction provides a good approximation of the direction of ∇2 f (x) · y. Figure 2: The intuition of ComparisonHessVec (Algorithm 1) for computing Hessian-vector products using gradient directions.

∇f (x−ry) ∥∇f (x−ry)∥

∇f (x) ∥∇f (x)∥

∇f (x+ry) ∥∇f (x+ry)∥

direction of ∇2 f (x) · y

x

Then, we use ComparisonHessVec to implement ComparisonHE (Algorithm 2), which estimates the normalized Hessian ∇2 f (x)/∥∇2 f (x)∥. Because a column of ∇2 f (x) is given by multiplying a unit vector by ∇2 f (x), we can use ComparisonHessVec to estimate the normalized direction of each column. The remaining challenge is to determine the relative norms of different columns. For any two columns hi and hj of ∇2 f (x), this is accomplished by estimating the directions of hi , hj , and hi + hj , and then solving the resulting triangle. A corner case of the above argument arises when two columns hi and hj have nearly identical directions, in which case the error incurred by solving the triangle can be arbitrarily large. To address this issue, we introduce an additional perturbation along a direction sufficiently far from that of hi , which ensures that the resulting error remains controllable. Furthermore, we give an algorithm ComparisonRatio (Algorithm 4) that estimates the ratio between the norm of the Hessian and the norm of the gradient using comparisons. This is done by identifying a vector v of the approximate Hessian whose corresponding eigenvalue has the largest magnitude, estimating the directions of ∇f (x) and ∇f (x + µv), and then solving the resulting triangle. One source of error in this approximation is that v may differ from u, the eigenvector of the true Hessian ∇2 f (x) corresponding to the eigenvalue with the largest magnitude. We control this discrepancy by applying the Davis–Kahan theorem to bound ∥u − v∥, which yields a bound on the overall estimation error. Finally, we design a trust-region–type algorithm (Algorithm 5) such that, at any iteration xt which is not an ϵ-stationary point, it produces a new point xt+1 that decreases the function value by Ω(ϵ3/2 ). As shown in standard trust-region analysis [50], there exists a point x̂ ∈ Rn with √ ∥x̂ − xt ∥ = O( ϵ) such that f (x̂) − f (xt ) ≤ −Ω(ϵ3/2 ). Moreover, such a point can be found exactly given access to the normalized gradient, the normalized Hessian, and the ratio of their norms. We show that the approximations of these quantities produced by our previous algorithms are sufficient to recover a point x̂′ with the same asymptotic decrease in function value. Hence, after running this algorithm for Θ(1/ϵ1.5 ) iterations, it decreases the function value for at least Ω(∆) in expectation, thus we get a contradiction with a constant probability. This contradiction means that we have visited at least one stationary point during the iteration steps. The above argument can fail in certain corner cases, e.g., when the Hessian is rank-one, positive semidefinite, or negative semidefinite, or when the eigenvector v corresponding to the eigenvalue of largest magnitude is nearly aligned 4

with the gradient. In these cases, we additionally apply normalized gradient descent together with a binary line search to guarantee an Ω(ϵ−1.5 ) decrease in function value. Open questions. Our work leaves several open questions for future investigation. • First, is it possible to interpolate between our approach and the STP framework [5, 19] to achieve a sharper tradeoff between the dependencies on n and ϵ? • Second, it is natural to extend our results to stochastic comparison oracles, including the standard Bradley–Terry model [8]. • Third, it is worth investigating the practical applicability of our algorithm in practical settings. • Fourth, our approach may benefit from lazy Hessian-update techniques [12, 25]. In Appendix E, e 2 + n3/2 /ϵ3/2 ) queries to a zeroth-order we show that a lazy trust-region method using O(n oracle. Establishing an analogous result in the comparison-oracle model is a natural direction for future work.

2

Preliminaries

Basic notation. For any x ∈ Rn , we denote g and H as the gradient and Hessian of f at x, respectively, and omit the dependence on x when it is clear from context. For any vector v and n (x) := {y ∈ Rn : ∥y − x∥ ≤ R} and matrix A, let ∥v∥ and ∥A∥ be their ℓ2 norm. We denote BR n n SR (x) := {y ∈ R : ∥y − x∥ = R}. Quantum notation. We briefly introduce some basic notation of quantum computing that will be used in this paper. We use the Dirac notation |·⟩ to represent quantum states, which can be seen as column vectors. Qubit, which is the basic unit of quantum state, is represented as |ϕ⟩ = α|0⟩+β|1⟩ on the computational basis {|0⟩, |1⟩} with α, β ∈ C and |α|2 +|β|2 = 1. An n-qubit system is represented in a 2nP -dimensional vector space, with basis states {|0⟩, |1⟩}⊗n . An n-qubit state can be written as P |ψ⟩ = x∈{0,1}n αx |x⟩ with x |αx |2 = 1. A quantum algorithm can access a function via queries to a quantum oracle. For a classical function ϕ, the oracle Oϕ is defined as a unitary transformation that maps |x⟩|b⟩ to |x⟩|b ⊕ ϕ(x)⟩. This allows the oracle to be queried on superpositions of inputs, producing corresponding superpositions of outputs. Gradient direction estimation using comparisons. Recently, Tao et al. [44] developed classical and quantum algorithms for gradient direction estimation using comparisons. These algorithms serve as basic building blocks in our algorithms, stated below: Theorem 3 (Classical gradient direction estimation, Theorem 3.3 in [44]). Let f : Rn → R with L1 -Lipschitz gradient. Given query access to a comparison oracle Ofcomp (1) and any precision ϵ, there exists an algorithm ComparisonGE(x, ϵ, γ) such that, for any x ∈ Rn with ∥∇f (x)∥ ≥ γ, it ∇f (x) b satisfying g b − ∥∇f outputs a unit vector g (x)∥ ≤ ϵ with success probability at least 2/3, using O(n log(1/ϵ)) queries. Theorem 4 (Quantum gradient direction estimation, Theorem 4.2 in [44]). In the setting of Theorem 3, given any precision ϵ, there exists an algorithm ComparisonQGE(x, ϵ, γ) such that, for ∇f (x) b satisfying g b − ∥∇f any x ∈ Rn with ∥∇f (x)∥ ≥ γ, it outputs a unit vector g (x)∥ ≤ ϵ with success probability at least 8/15 − 2ϵ, using O(log(n/ϵ)) queries. 5

3

Hessian-Vector Product Estimation by Comparisons

In this section, we introduce ComparisonHessVec (Algorithm 1), an algorithm for approximating the direction of a Hessian–vector product using comparisons. Algorithm 1: ComparisonHessVec Input: f : Rn→ R, x, y ∈ Rn , precision  δ̂, γx , γy . √ √ γ δ̂ ϵ δ̂ γ γx , γx , √x , 20y√L . 1 Set r0 ← min 100L 1 100L2 20 L 2

2

L2 r2 ρr2 2 ĝ0 ← ComparisonGE(x, γ 0 , γx ), ĝ1 ← ComparisonGE(x + r0 y, γ 0 , γx /2), x x L r2 ĝ−1 ← ComparisonGE(x − r0 y, γ2x 0 , γx /2). p p 3 Set g = 1 − ⟨ĝ−1 , ĝ0 ⟩2 ĝ1 − 1 − ⟨ĝ1 , ĝ0 ⟩2 ĝ−1 . b (x, y, δ, γx , γy ) = g/∥g∥. 4 return u

Theorem 5. Let f : Rn → R with L1 -Lipschitz gradient and L2 -Lipschitz Hessian. For any γx , γy > 0 and x, y ∈ Rd satisfying p ∥∇f (x)∥ ≥ γx , ∥∇2 f (x)∥ ≥ L2 ϵ, ∥y∥ = 1,

|y1 | ≥ γy ,

where y1 := ⟨y, u1 ⟩ with u1 being a eigenvector of ∇2 f (x) with the largest magnitude eigenvalue, b satisfying Algorithm 1 outputs a vector u b− u

∇2 f (x) · y ≤ δ̂ ∥∇2 f (x) · y∥

using Õ n log 1/(γx γy2 δ̂ 2 ) queries. Intuitively, Algorithm 1 use Taylor expansion ∇f (x ± ry) ≈ ∇f (x) ± r∇2 f (x) · y and solving the triangle in Figure 3. Proof. Denote x+ := x + r0 y and x− := x − r0 y. Since f has L2 -Lipschitz Hessian, ∇f (x+ ) − ∇f (x) − r0 ∇2 f (x) · y ≤ L2 r02 /2, 2

∇f (x− ) − ∇f (x) + r0 ∇ f (x) · y

≤ L2 r02 /2.

(3) (4)

(x− ) Therefore, we have ∥∇f (x+ ) + ∇f (x− ) − 2∇f (x)∥ ≤ L2 r02 and ∇2 f (x) · y − ∇f (x+ )−∇f ≤ 2r0 L2 r0 /2. Furthermore, as r0 ≤ γx /(100L1 ) and f has L1 -Lipschitz gradient, both ∥∇f (x+ )∥ and ∥∇f (x− )∥ is at least γx −L1 γx /(100L1 ) = 0.99γx . We first understand how to approximate ∇2 f (x)·y ∇f (x+ ) ∇f (x− ) ∇f (x) by normalized vectors ∥∇f (x)∥ , ∥∇f (x+ )∥ , ∥∇f (x− )∥ , and then analyze the approximation error due to using ĝ0 , ĝ1 , ĝ−1 , respectively. By Lemma 15, we have

1 ∥∇f (x) − r0 ∇2 f (x) · y∥ r D E2 2∥∇f (x)∥ ∇f (x)+r0 ∇2 f (x)·y ∇f (x) 1 − ∥∇f , 2 (x)+r0 ∇ f (x)·y∥ ∥∇f (x)∥ =

1 ∥∇f (x) + r0 ∇2 f (x) · y∥ r D E2 . 2∥∇f (x)∥ ∇f (x)−r0 ∇2 f (x)·y ∇f (x) 1 − ∥∇f (x)−r0 ∇2 f (x)·y∥ , ∥∇f (x)∥ 6

We denote the value above as α. Because f is L2 -Hessian Lipschitz, ∥r0 ∇2 f (x) · y∥ ≤ r0 L2 . Since γx γx , ∥r0 ∇2 f (x) · y∥ ≤ 100 r0 ≤ 100L . Also by Lemma 14 we have 2   ∇f (x) + r0 ∇2 f (x) · y ∇f (x) ≥ 0.94, , ∥∇f (x) + r0 ∇2 f (x) · y∥ ∥∇f (x)∥   ∇f (x) − r0 ∇2 f (x) · y ∇f (x) ≥ 0.94, , ∥∇f (x) − r0 ∇2 f (x) · y∥ ∥∇f (x)∥ 0.99 which promises that α ≥ 2√1−0.94 ≥ 1. In arguments next, we say a vector u is ξ-close to a vector 2 v if ∥u − v∥ ≤ ξ for any ξ ≥ 0. We prove that the vector s   ∇f (x) ∇f (x− ) ∇f (x) 2 ∇f (x+ ) g̃1 := 1− +α· , ∥∇f (x)∥ ∥∇f (x− )∥ ∥∇f (x)∥ ∥∇f (x+ )∥ s !   ∇f (x+ ) ∇f (x) 2 ∇f (x− ) (5) − 1− , ∥∇f (x+ )∥ ∥∇f (x)∥ ∥∇f (x− )∥ 7L2 r02 γx -close to a vector proportional to ∇f (x+ ). This is because (3), (4), and Lemma 13 imply L2 r02 ∇f (x+ ) ∇f (x)+r0 ∇2 f (x)·y -close to each other, that ∥∇f are 0.99γ and ∥∇f (x+ )∥ (x)+r0 ∇2 f (x)·y∥ x

is

s 1−

∇f (x− ) ∇f (x) , ∥∇f (x− )∥ ∥∇f (x)∥

2

∇f (x+ ) ∥∇f (x+ )∥

(6)

is proportional to ∇f (x+ ), and the definition of α implies s   2 ∇f (x) − r0 ∇ f (x) · y ∇f (x) + r0 ∇2 f (x) · y ∇f (x) 2 ∇f (x) −α· · 1 − , ∥∇f (x)∥ ∥∇f (x) − r0 ∇2 f (x) · y∥ ∥∇f (x) + r0 ∇2 f (x) · y∥ ∥∇f (x)∥ (7)

= ∇f (x) + r0 ∇2 f (x) · y/(2∥∇f (x)∥). L r2

∇f (x+ ) The above vector is 4γ2 x0 -close to 2∥∇f (x)∥ by (3), and the error in above steps cumulates by at most 6L2 r02 6L2 r02 L2 r02 7L2 r02 0.99γx using Lemma 14. In total, the error is at most 0.99γx + 4γx ≤ γx . Furthermore, this vector L r2 proportional to ∇f (x+ ) that is 4γ2 x0 -close to (5) has norm at least (1 − 0.01)/2 = 0.495 because the γx coefficient in (6) is positive, while in the equality above we have ∥r0 ∇2 f (x) · y∥ ≤ 100 . Therefore, applying Lemma 13, the vector g̃1 in (5) satisfies

g̃1 /∥g̃1 ∥ − ∇f (x+ )/∥∇f (x+ )∥ ≤ 29L2 r02 /γx . Following the same logic, we can prove that the vector s   ∇f (x) ∇f (x− ) ∇f (x) 2 ∇f (x+ ) g̃−1 := −α 1− , ∥∇f (x)∥ ∥∇f (x− )∥ ∥∇f (x)∥ ∥∇f (x+ )∥ s   ∇f (x+ ) ∇f (x) 2 ∇f (x− ) +α 1− , ∥∇f (x+ )∥ ∥∇f (x)∥ ∥∇f (x− )∥

(8)

(9)

satisfies g̃−1 /∥g̃−1 ∥ − ∇f (x− )/∥∇f (x− )∥ ≤ 29L2 r02 /γx . 7

(10)

Furthermore, (7) implies that the vector g̃1 − g̃−1 is

14L2 r02 γx -close to

∇f (x) + r0 ∇2 f (x) · y ∇f (x) − r0 ∇2 f (x) · y − = r0 ∇2 f (x) · y/∥∇f (x)∥. (11) 2∥∇f (x)∥ 2∥∇f (x)∥ √ √ Given that ∥∇2 f (x)∥ ≥ L2 ϵ and |y1 | ≥ γy , ∥∇2 f (x) · y∥ ≥ L2 ϵγy . Therefore, the RHS of (11) √ has norm at least

r0 L2 ϵγy , and by Lemma 13 we have γx

√ √ g̃1 − g̃−1 ∇2 f (x) · y 14L2 r02 . r0 L2 ϵγy  14r0 L2 √ = . ≤ − ∥g̃1 − g̃−1 ∥ ∥∇2 f (x) · y∥ γx γx ϵγy ∇f (x) Finally, by Theorem 3, the error terms coming from ComparisonGE, ĝ0 − ∥∇f (x)∥ ,

(12) ∇f (x+ ) ĝ1 − ∥∇f (x+ )∥ ,

L2 r02 γx . Combined with (8) and (10), we know that 61L2 r02 −1 the vector g we obtained in Algorithm 1 is γx close to g̃1 −g̃ 2α . Since α ≥ 1, by Lemma 13 we

and

∇f (x− ) ĝ−1 − ∥∇f (x− )∥ , are all upper bounded by

have

g g̃1 − g̃−1 61L2 r02 − ≤ . ∥g∥ ∥g̃1 − g̃−1 ∥ γx

(13)

Combining (12) and (13), √ g ∇2 f (x) · y 61L2 r02 14r0 L2 √ − ≤ + . (14) ∥g∥ ∥∇2 f (x) · y∥ γx ϵγy √ √  γx γ δ̂ ϵ γx Our choice of r0 = min 100L1 , 100L2 , 20√γxLδ̂ , 20y√L guarantees that the RHS of (14) is at most δ̂. 2 2 In terms of query complexity, we make 3 calls to ComparisonGE. By Theorem 3, the total query

complexity is O n log nL2 L21 /γx γy2 ϵδ̂ 2 .

4

Robust Hessian Estimation by Comparisons

In this section, we show how to obtain an estimate of the normalized Hessian at a given point using comparisons. √ Theorem 6. Let x satisfy ∥H∥ ≥ L2 ϵ and
∥g∥ ≥ ϵ. There exist universal constants c0 , K > 0 b satisfying such that, Algorithm 2 makes Õ(n2 log 1/δ ) queries and outputs H √ b− H ≤ K n η 1/4 . H ∥H∥ with probability at least 2/3. We divide the analysis into a non-degenerate case described in Assumption 1 and two degenerate cases: either hi has a sufficiently small norm, or the directions of h1 and hi are too close. We handle the non-degenerate case in Section 4.1 and handle the two degenerate cases in Section 4.2 separately. Denote H = (h1 , . . . , hn ) with hi = Hei and define ri := ∥hDi ∥/∥h1 ∥ ∈ [0, 1]. loss of genE Without D E hi +hj hj h +h hi erality, we assume∥He1 ∥ = maxi ∥Hei ∥. Otherwise, by testing ∥hi ∥ , ∥hi +hj ∥ and ∥hj ∥ , ∥hii +hjj ∥ , we can roughly tell which one is larger and then relabel columns. The error only incurs a constant overhead in the complexity. 8

Algorithm 2: ComparisonHE(x, δ): Estimate Normalized Hessian at x with Precision δ by Comparisons Input: Target point x, target accuracy δ > 0 b Output: H √ 1 1 Set η ← c0 δ 4 /n2 , thresholds τα , τβ ← Θ( η), σ ← 8 η 2 , sample t ∼ Unif({1, . . . , n}). e ±σet ± 2 Define yi ← √i for all i ∈ [n]. 1+σ 2

e +e ±σe ± 3 Define zi ← 1√ i 2 t for all i ∈ {2, . . . , n}. 2+σ 4 for i = 1, . . . , n do √ √ b (x, ei , η/4, ϵ, η/2 n), gi+ ← u b (x, yi+ , η/4, ϵ, η/2 n), 5 Query gi0 ← u √ b (x, yi− , η/4, ϵ, η/2 n). gi− ← u + − 6 if ⟨gi , gi ⟩ ≤ 1 − τα then 7 Set gi ← 0. 8 9 10

else (v1 , v2 ) ← arg maxv1 ̸=v2 ∈{g0 ,g+ ,g− } ⟨v1 , v2 ⟩. +v2 Set gi ← ∥vv11 +v . 2∥

i

i

i

11 for i = 2, . . . , n do √ √ + 0 ←u b (x, e1 + ei , η/4, ϵ, η/2 n), g1i b (x, z+ 12 Query g1i ←u i , η/4, ϵ, η/2 n), √ − b (x, z− g1i ←u i , η/4, ϵ, η/2 n). + − 13 if ⟨g1i , g1i ⟩ ≤ 1 − τβ then 14 15 16 17

Set g1i ← 0. else (v1 , v2 ) ← arg maxv1 ̸=v2 ∈{g0 ,g+ ,g− } ⟨v1 , v2 ⟩. +v2 Set g1i ← ∥vv11 +v . 2∥

1i

1i

1i

18 Set r̂1 ← 1. 19 for i ← 2 to n do 20 21 22 23 24 25

α bi ← ⟨gi , g1 ⟩, βbi ← ⟨g1i , g1 ⟩. if βbi ≥ 1 − τβ then if α bi ≤ 1 − τα then r̂i ← 0. else r̂i ← PerturbAndSolve(i, η).

else Solve the following equation with (α, β) = (b αi , βbi ) and take the unique root r̂i ∈ [0, 1]. 2 2 2 2 2 28 (β − α )r̂i + 2α(β − 1)r̂i + (β − 1) = 0.
e ← 1 [r̂1 g1 , . . . , r̂n gn ] + [r̂1 g1 , . . . , r̂n gn ]⊤ . 29 H 2 b ← H/∥ e H∥. e 30 return H 26

27

9

Algorithm 3: PerturbAndSolve(i, η) 1 Set ρ ← η 1/3 . 2 for j ← 2 to n do 3 4 5 6 7 8 9 10 11 12

if j ̸= i and α bj ≤ 1 − τα and α b1j ≤ 1 − ρ2 then √ b (x, ei + ρej , η/4, ϵ, η/2 n). Query gi,j ← u α bi,j ← ⟨gi,j , g1 ⟩. √ b (x, e1 + ei + ρej , η/4, ϵ, η/2 n). Query g1,i,j ← u βbi,j ← ⟨g1,i,j , g1 ⟩. if ⟨gi,j , g1 ⟩ ≤ 1 − ρ then return r̂i ← 0. Solve the following equation with (α, β) = (b αi,j , βbi,j ) and take r̂i,j ∈ [0, 1 + ρ]. (β 2 − α2 )r̂i2 + 2α(β 2 − 1)r̂i + (β 2 − 1) = 0. return r̂i ← r̂i,j .

13 return r̂i ← 1.

Our estimator first estimates the directions hi /∥hi ∥ and (h1 + hi )/∥h1 + hi ∥ via Hessian-vector product queries at (perturbed) ei and e1 + ei , and then computes r̂i by applying Lemma 15 in parallelogram spanned by h1 , hi and h1 + hi in non-degenerate cases that the direction of h1 and hi are not close. ̸ 0 and The computation of r̂i can be roughly stated as follows: for each i ≥ 2 with hi = h1 + hi ̸= 0, define unit vectors u1 :=

h1 , ∥h1 ∥

ui :=

hi , ∥hi ∥

u1i :=

h1 + hi . ∥h1 + hi ∥

Let αi := ⟨ui , u1 ⟩ and βi := ⟨u1i , u1 ⟩. A direct computation yields 1 + r i αi βi = F (ri , αi ) := q . 1 + ri2 + 2ri αi

(15)

Squaring (15) gives the quadratic equation solved in Line 28 in Algorithm 2 (βi2 − αi2 )r̂i2 + 2αi (βi2 − 1)r̂i + (βi2 − 1) = 0.

(16)

b (see Line 29). By solving this equation, we can roughly estimate ri , thus we can form the target H The analysis has two basic cases: (i) When h1 and hi overlap relatively small in directions and their norms are not too small, i.e. αi and βi are bounded away from 1, we can solve (16) with bounded error on |r̂i − ri |. Otherwise, when αi and βi are close to 1, (16) becomes degenerate. There are two cases that may cause degeneracy: (ii) When ∥hi ∥ is relatively small, we can tell this case by testing if the angle between u1 and ui is relatively large, and we can directly treat ∥hi ∥ = 0 since the error of column vectors cause bounded error to Hessian estimation. (iii) When h1 and hi have close directions, to avoid uncontrollable error, we apply Algorithm 3 to enforce an angle separation 1 − α2 = Ω(τα ), in which we use perturbed ei + ρej to replace ei , and perturbed e1 + ei + ρej to replace e1 + ei to guarantee the non-degeneracy in the computation of the ratio r̂i . This makes the √ inverse map Lipschitz with coefficient O(1/ρ2 ), giving a noise term O(η/ τα ρ2 ). The perturbation √ √ changes the effective triangle by O(ρ/ τα ), producing a bias term O(ρ/ τα ). Hence 1 η |r̂i − ri | ≤ √ (C 2 + Bρ), τα ρ 10

optimizing it gives |r̂i − ri | = O(η 1/4 ). (See Lemma 4.) √ With constant probability (guaranteed by Lemma 18), our sampled t satisfies |⟨et , v⟩| ≥ 1/ n, where v is the eigenvector corresponding to minimum eigenvalue of H, here we use the isotropy of n √ coordinates. Once we can satisfy a 1/ n-overlap, our choice of gi and g1i satisfies ∥gi − ui ∥ ≤ η and ∥g1i − u1i ∥ ≤ η according to Lemma 20. Inner products satisfy |b αi − αi | ≤ O(η) and |βbi − βi | ≤ O(η). When αi is bounded away from 1 and ri is not tiny, the implicit map β = F (r, α) has slope |∂F/∂r| = Ω(1/r), so solving the quadratic (16) gives |r̂i − ri | = O(η/rmin τα ). (See Lemma 1.) e satisfying ∥H e − H ⋆ ∥ ≤ O(√nη 1/4 ), Finally, assembling columns r̂i gi yields a scale-free matrix H b − H . (See Lemma 16 and and normalization stability implies the same order bound for H ∥H∥ Lemma 17.)

4.1

The non-degenerate case

Here we give a deterministic bound for the case that directly solves (16) without perturbation. Assumption 1 (Non-degenerate case). Fix an index i ≥ 2 such that hi ̸= 0 and h1 + hi = ̸ 0. Assume: 1. (Not too parallel) αi ≤ 1 − γ for some γ ∈ (0, 1). 2. (Not too small) ri ≥ rmin for some rmin ∈ (0, 1]. Lemma 1. Under Assumption 1, let (α̂, β̂) satisfy |α̂ − αi | ≤ ε and |β̂ − βi | ≤ ε, with ε ≤ γ/10. Let r̂ be the (unique) solution in [0, 1] of β̂ = F (r, α̂) introduced in (15). Then |r̂ − ri | ≤

C0 rmin γ

ε,

for some universal constant C0 . Proof. By the mean value theorem, |F (r̂, α̂) − F (ri , α̂)| = |

∂F (ξ, α̂)| |r̂ − ri | ∂r

for some ξ between r̂ and ri . Therefore, |r̂ − ri | =

|β̂ − F (ri , α̂)| . | ∂F ∂r (ξ, α̂)|

We bound numerator and denominator respectively. By triangle inequality, |β̂ − F (ri , α̂)| ≤ |β̂ − βi | + |F (ri , αi ) − F (ri , α̂)| ≤ ϵ + Lα ϵ,

(17)

where Lα is a Lipschitz constant of F in α on the compact set r ∈ [0, 1], α ∈ [−1, 1]. We can bound Lα explicitly by computing ∂F/∂α and taking a supremum; a crude universal bound Lα ≤ 2 suffices (indeed, one can check |∂F/∂α| ≤ 2 on this domain). Thus the numerator is less than 3ϵ. On the other hand, We have ∂F ξ(1 − α̂2 ) (ξ, α̂) = . ∂r (1 + ξ 2 + 2ξ α̂)3/2

11

Since ξ lies between r̂ ∈ [0, 1] and ri ≥ rmin , we have ξ ≥ rmin /2 provided ϵ is small enough; we can ensure this by restricting to ϵ ≤ γ/10 and using continuity of the inverse map. Also, γ 9 α̂ ≤ αi + ϵ ≤ 1 − γ + 10 = 1 − 10 γ, so
9 9 γ · 1 = 10 γ. 1 − α̂2 = (1 − α̂)(1 + α̂) ≥ 10 Therefore

9 (rmin /2) · ( 10 γ) ∂F (ξ, α̂) ≥ ≥ c rmin γ ∂r 8

for a universal constant c > 0. Combining numerator and denominator gives |r̂ − ri | ≤

3ε C0 = ε c rmin γ rmin γ

with C0 := 3/c. Corollary 1. Under Assumption 1, the variables α bi , βbi in (16) of Algorithm 2 satisfies |b αi − αi | ≤ 3η b and |βi − βi | ≤ 3η, and thus C1 η |r̂i − ri | ≤ rmin τα for a constant C1 .

4.2

The degenerate cases

In this subsection, we discuss about how to deal with the two degenerate cases that lead to the failure of Lemma 1. In the case that hi has a sufficiently small norm, we introduce the following lemma. Lemma 2. If r̂i is set to 0, which means the algorithm enters this case only when the column is treated as negligible, i.e. ri ≤ rmin . We have ri ≤ rmin . Proof. If r̂i is set to 0, by Line 21 and Line 22 in Algorithm 2, we have α b i ≤ 1 − τα . Geometry tells us τα = cos(arctan rmin )

(18)

Plugging in α bi = ⟨gi , g1 ⟩ and (18) we get ri ≤ rmin .

In the case that the directions of h1 and hi are too close, we can find an index j such that the direction of hj and the direction of h1 (so to hi ) are sufficiently separated for the perturbation step. Note that if we cannot find the index j, it means that the rank of H is 1 up to a small error, which would be dealt with by line search step in Lemma 9, so we have the following assumption.

12

Assumption 2. Fix an index i such that hi ̸= 0 and h1 + hi ̸= 0. Assume ri ≥ rmin > 0. Algorithm 3 finds an index j such that   hi + ρhj h1 2 1− , ≥ c0 ρ2 (19) ∥hi + ρhj ∥ ∥h1 ∥ for some absolute c0 > 0. iv Define α := ⟨v, u1 ⟩ , w := ∥uu11 +r +ri v∥ , β := ⟨w, u1 ⟩, then β = F (ri , α). Algorithm 3 estimates v and w by oracle queries

b (ei + ρej ; η), v̂ := u

b (e1 + ei + ρej ; η), ŵ := u

and computes α̂ := ⟨v̂, g1 ⟩ ,

β̂ := ⟨ŵ, g1 ⟩ ,

then returns r̂ as the solution in [0, 1 + ρ] to β̂ = F (r̂, α̂). Lemma 3. Under Assumption 2, there exists a constant C such that for η sufficiently small, η |r̂ − ri | ≤ C . rmin ρ2 Proof. By Theorem 5, ∥v̂ − v∥ ≤ η, ∥ŵ − w∥ ≤ η, and also ∥g1 − u1 ∥ ≤ η. Apply Lemma 14 we have |α̂ − α| ≤ 3η, |β̂ − β| ≤ 3η. (20) and

rmin (1 − α2 ) ∂F ri (1 − α2 ) (ri , α) = ≥ . 2 ∂r 8 (1 + ri + 2ri α)3/2

Using the separation 1 − α2 ≥ c0 ρ2 , we obtain ∂F rmin c0 2 (ri , α) ≥ ρ ∂r 8

(21)

By the mean value theorem, |r̂ − ri | =

|β̂ − F (ri , α̂)| | ∂F ∂r (ξ)|

for some ξ between ri and r̂. For η sufficiently small, α̂ is close to α and thus 1 − α̂2 ≥ 12 (1 − α2 ) ≥ 1 2 2 c0 ρ , and ξ remains in a compact interval containing [rmin /2, 1 + ρ]. Therefore, the lower bound (21) holds for | ∂F ∂r (ξ)| as well ∂F rmin c0 2 (ξ, α̂) ≥ ρ . ∂r 16 Hence 16 |r̂ − ri | ≤ |β̂ − F (ri , α̂)|. rmin c0 ρ2 Similar to (17), we have |β̂ − F (ri , α̂)| ≤ |β̂ − β| + |F (ri , α) − F (ri , α̂)| ≤ 3η + 3Lα η. where Lα is a Lipschitz constant of F in α on the compact set r ∈ [0, 1], α ∈ [−1, 1]. Conclude that the numerator is ≤ C ′ η, and therefore η |r̂ − ri | ≤ C . rmin ρ2

13

Lemma 4. Under Assumption 2 and assuming ri ≥ rmin , the perturbation case estimator satisfies   η 1 C 2 + Bρ |r̂ − ri | ≤ rmin ρ for constants C, B depending only on c0 . Proof. The first term Cη/ρ2 follows from Lemma 3. The second term accounts for the fact that the perturbation uses v (direction of hi + ρhj ) rather than ui (direction of hi ). By Lemma 21, ∥v−ui ∥ ≤ O(ρ) (with explicit constant 4/rmin ). On the non-degenerate set enforced by 1−α2 ≥ c0 ρ2 , the mapping from input directions to the recovered ratio is Lipschitz; therefore this O(ρ) direction change induces an O(ρ) bias in r. Combining the two terms of error yields the stated bound.

4.3

Proof of Theorem 6

Here complete the proof of Theorem 6, the main theorem of this section. √ Proof. By Lemma 18, we have Pru∼Unif(S1 (0)) [|⟨u, v⟩| ≥ 1/ n] ≥ 2/5 where v is the eigenvalue of H with the largest absolute eigenvalue. Because of the isotropy of n coordinates, v ∼ Unif(S1 (0)). So we have √ √ Pr [|⟨u, v⟩| ≥ 1/ n] Pr [|⟨u, v⟩| ≥ 1/ n] = u∼Unif(S1 (0)),v∼Unif(S1 (0)) u∼Unif(S1 (0)) √ = Pr [|⟨u, v⟩| ≥ 1/ n] u∼Unif({e1 ,...,en }),v∼Unif(S1 (0)) √ = Pr [|⟨et , v⟩| ≥ 1/ n] t∼Unif({1,...,n}),v∼Unif(S1 (0))

≥ 2/5. Repeat for O(1) times and the success probability can be more than 2/3. Let H ⋆ be the scale-free e = [r̂1 g1 , . . . , r̂n gn ] and define E := H e − H ⋆. target matrix (66), so that H/∥H∥ = H ⋆ /∥H ⋆ ∥. Let H Fix i. If r̂i is produced by the non-degenerate solver, Corollary 1 gives |r̂i − ri | ≤

C η rmin τα

for some constant C. In the case that hi has a sufficiently small norm, by Lemma 2 we have |r̂i − ri | < rmin . In the case that the directions of h1 and hi are too close, (18) tells us rmin = tan (arccos(1 − τβ )). When η is sufficiently small, through Taylor expansion, f (x) = tan(arccos(1 − x)) ∼ √ i.e. rmin ∼ τβ . If r̂i is produced by perturbation, then Lemma 4 yields   1 η |r̂i − ri | ≤ C 2 + Bρ . rmin ρ √ Combine these three cases and set ρ = Θ(η 1/3 ), τα = τβ = Θ( η), we get |r̂i − ri | ≤ O(η 1/4 ), 14

√ 2x + O(x3/2 ),

and ∥Eei ∥ = ∥r̂i gi − ri ui ∥ ≤ ∥(r̂i − ri )gi ∥ + ∥ri (gi − ui )∥ ≤ |r̂i − ri | · ∥gi ∥ + ri η ≤ 2|r̂i − ri | + η, since ∥gi ∥ ≤ 1 + η ≤ 2 and ri ≤ 1. Thus ∥Eei ∥ ≤ K1 η 1/4 for a constant K1 . By Lemma 17, √ ∥E∥ ≤ K1 n η 1/4 . By Lemma 19, ∥H ⋆ ∥ ≥ 1. If η is small enough that ∥E∥ ≤ 1/2, then Lemma 16 gives e √ H H⋆ − ≤ 4∥E∥ ≤ 4K1 n η 1/4 . ⋆ e ∥H ∥ ∥H∥ Using H/∥H∥ = H ⋆ /∥H ⋆ ∥ then gives the stated bound with K := 4K1 . Corollary 2. Under the setting of Theorem 6, we can increase success probability to at least 1 − p using O(n2 log(1/δ) log(1/p)) queries for any p ∈ (0, 1) to get the same result. b (r) to be Proof. Consider the m = ⌈18 ln(1/p)⌉ independent base runs in the algorithm. Define H (r) (s) b b the output of the r-th iteration and set Nr = {s ∈ {1, . . . , m} : ∥H − H ∥ ≤ 2δ/3}. We finally b (r) with |Nr | ≥ m/2 and prove that it satisfies our requirements with success probability return any H b (r) − H∥ ≤ δ0 . Since at least 1 − p. Let Xr be the indicator that the r-th candidate satisfies ∥H Pr[Xr = 1] ≥ 2/3, Hoeffding’s inequality gives "m # X Pr Xr ≤ m/2 ≤ exp(−m/18) ≤ p. r=1

Conditioned on the complementary event, more than half of the candidates are good. Any two good candidates are within 2δ0 = 2δ/3 of each other, so every good candidate satisfies the selection condition |Nr | ≥ m/2. Thus the algorithm can select at least one candidate. Moreover, for any b (r) , the set Nr intersects the set of good candidates because both sets have size selected candidate H greater than m/2 or at least m/2 with a strict majority of good candidates. Therefore, for some b (s) , good H b (r) − H∥ ≤ ∥H b (r) − H b (s) ∥ + ∥H b (s) − H∥ ≤ 2δ/3 + δ/3 = δ. ∥H The total query complexity is m · O(n2 log 1δ ) = O(n2 log 1δ log p1 ).

5

Estimate the Ratio of the Hessian and the Gradient Norm 2

f (x)∥ In this section, we show that, under the following assumption, the ratio ∥∇ ∥∇f (x)∥ can be estimated using a single call to ComparisonHE (Algorithm 2) and two calls to ComparisonGE (Theorem 3).

Assumption 3.

1. The Hessian H is neither positive definite nor negative definite.

2. The Hessian H is not rank-one. 3. There exists a unit eigenvector u of H with eigenvalue λu such that   g arccos , u ≥ β and ∥Hu∥ = |λu | ≥ λ, ∥g∥ for some β ∈ (0, π/2) and λ > 0. 15

(22)

Algorithm 4: ComparisonRatio(x, δκ ): Estimate Gradient-Hessian Ratio at point x with precision δκ by Comparisons Input: Point x, lower bound of gradient ϵ precision δκ , degenerate parameters β, λ Output: Estimate κ̂ ≈ ∥H∥/∥g∥ 1 . 2. 2 µ ← tan2 β · ϵ2 δκ 2 /n)). 3 (ξ1 , ξ2 , ξ3 ) ← (O(ϵδκ ), O(ϵδκ ), O(ϵ2 δκ b ← ComparisonHE(x; ξ3 ). 4 H b associated with the eigenvalue λv with the largest 5 Compute a unit eigenvector v of H absolute value b1 ← ComparisonGE(x; ξ1 ). 6 g b2 ← ComparisonGE(x + µv; ξ2 ). 7 g q b2 ⟩2 1 − ⟨b g1 , g q 8 κ̂ ← . µ 1 − ⟨b g2 , v⟩2 9 return κ̂.

Theorem 7. Let f : Rn → R with L1 -Lipschitz gradient and L2 -Lipschitz Hessian. Assume that Assumption 3 holds. Given any precision δκ > 0, there exists an algorithm ComparisonRatio(x, δκ ) such that, for any x ∈ Rn with ∥g∥ ≥ ϵ, it outputs κ̂ satisfying κ̂ −

∥H∥ ≤ δκ . ∥g∥

(23)

using Õ(n2 log 1/δκ ) queries. A

∇f (x)

∇f (x + µv)

∇2 f (x) · µu + o(µ) B C 2 ∇ f (x) · µu E D Figure 3: The intuition for computing the ratio.

By Lemma 25, we have the second-order expansion for some small µ, ∇f (x + µv) = g + µHv + rµ ,

∥rµ ∥ ≤ 12 L2 µ2 .

(24)

b corresponding to the eigenvalue λv , which is the largest Here we choose v as an eigenvector ofDH E b magnitude eigenvalue satisfying arccos ∥gg11 ∥ , v ≥ β and ∥Hv∥ = |λv | ≥ λ. Otherwise, when finding the stationary point, we apply line search on direction v to handle this case. As shown in 16

Figure 3, we can use line segment BE to replace BC. By law of sines, we can get an estimation of the Hessian-gradient norm ratio ∥∇2 f (x) · µu∥ µλu ∥∇2 f (x)∥ = ≈µ . ∥∇f (x)∥ ∥∇f (x)∥ ∥∇f (x)∥ where u is an eigenvector of the real Hessian H appropriately chosen. A core lemma we use here is the Davis-Kahan Theorem (Lemma 22), which bounds the distance of two subspaces spanned by eigenvectors under perturbation. We use this tool to bound ∥u − v∥ and thus get an upper bound of the error |κ̂ − ∥H∥ ∥g∥ |. We will rigorously describe how to choose u in Section 5.1, here we directly claim that u has two properties:   3 (i) ∥u − v∥ ≤ O nξ ϵδκ ; (ii) ∥H∥ − ∥Hu∥ ≤ ϵδ2κ + 2ξ3 . ξ3 We can take ξ3 sufficiently small to apply Lemma 22, which gives ∥ sin(Θ(U, V ))∥ ≤ λ′ −λ ′ k

k+1

. With

these two properties we can get an upper bound of |κ̂ − ∥H∥ ∥g∥ |.

5.1

Analysis of errors when estimating the Hessian-gradient norm ratio

Here we give some lemmas to prove the properties of u, and first we describe how to choose u. b Assume that all n Denote v as the eigenvector with the largest magnitude eigenvalue of H. ′ ′ b eigenvalues of H are λ1 ≥ · · · ≥ λn (correspond to eigenvectors v1 , . . . , vn ), and due to Assumption 3, we can assume λ′1 > 0 and λ′n < 0. Without loss of generality, we suppose that λv = λ′1 here, another H∥ case is the same. By the Dirichlet Theorem, there exists λ′k − λ′k+1 ≥ ∥n−1 for some index k. We ′ ′ ′ b denote k as the first index that λk − λk+1 ≥ ϵδκ /2n. We call the set Λ = {λ1 , . . . , λ′k } a cluster, and denote V = span(v1 , . . . , vk ). Denote the eigenvalues of Hessian H as λ1 ≥ · · · ≥ λn (correspond to eigenvectors u1 , . . . , un ). By Weyl’s Theorem (Lemma 23) we know that |λi − λ′i | ≤ ξ3 . Similarly define the cluster Λ = {λ1 , . . . , λk′ } and space U = span(u1 , . . . , uk′ ). Denote the projection of v in U as u. b

Lemma 5. With u being chosen as above, we guarantee ∥Hv∥ − ∥Hu∥ ≤

ϵδκ + ξ3 . 2

As a result, we have ∥H∥ − ∥Hu∥ ≤

ϵδκ + 2ξ3 . 2

(25)

b − H∥ ≤ ξ3 , by Davis-Kahan Theorem (Lemma 22), we have k = k ′ and Proof. Since ∥H  
ξ3 nξ3 sin Θ(U, V ) ≤ =O , λk − λk+1 ϵδκ and hence

 ∥u − v∥ ≤ O

17

nξ3 ϵδκ

 .

(26)

Figure 4: The intuition for applying the Davis-Kahan Theorem (Lemma 22) to prove ∥H∥−∥Hu∥ ≤ ϵδ2κ +2ξ3 .

∥H∥

∥Hu∥ξ3

ξ3

eigenvalues of H, as a sequence λ′1 ≥ · · · ≥ λ′n

Projection b λ′k λv (∥H∥)

ϵδκ 2n

λ′k+1

b as a sequence λ1 ≥ · · · ≥ λn eigenvalues of H,

For our choice of u, We prove ∥H∥ − ∥Hu∥ ≤ ϵδ2κ + 2ξ3 . Our intuition is shown in Figure 4. We κ find the first pair of eigenvalue (λ′k , λ′k+1 ) that are far apart (distance longer than ϵδ 2n ), thus we can ′ ′ say λ1 , . . . , λk are close enough, and so λ1 , . . . , λk are close too. Then we can apply Davis-Kahan P Theorem to bound the distance between U and V . Since u ∈ U , we can write u = ki=1 si ui where Pk 2 i=1 si = 1. Therefore, k k X X Hu = si Hui = si λi ui , i=1

and

i=1

v v u k u k uX uX s2 λ2 ≥ λ t s2 = λ . ∥Hu∥ = t i

i

i=1

k

i

k

i=1

κ Since for any i < k we have λ′i − λ′i+1 ≤ ϵδ 2n , we can further derive that

λ′1 − λ′k ≤

k−1 X

λ′i − λ′i+1 ≤

i=1

(k − 1)ϵδκ ϵδκ ≤ . 2n 2

From |λi − λ′i | ≤ ξ3 , ∥H∥ − ∥Hu∥ = λ1 − ∥Hu∥ ≤ λ′1 + ξ3 − (λ′k − ξ3 ) ≤

ϵδκ + 2ξ3 . 2

(27)

Lemma 6. In the setting above and choose u as in Lemma 5, if√ two vectors ∇f (x + µv) + k2 v2 µ 3 and ∇f (x) + k1 v1 + η Ĥ · µv are colinear, then η = ∥H∥(1 + O( tan β ) + O( nξ ϵδκ )). Proof. Note that 1 ∇f (x + µv) = ∇f (x) + ∇2 f (x) · µv + O(µ2 ). 2 From the observation in Figure 3, if there exist η ′ such that two vectors ∇f (x + µv) + k2 v2 and ∇f (x) + k1 v1 + η ′ H · µu are colinear, then η′ =

|BD| |CD| sin ∠BCD sin ∠CED sin ∠BCD = · = · . |BE| |CE| sin ∠BCE sin ∠CDE sin ∠BCE

18

(28)

We compute two ratios in (28) separately √
sin ∠BCE + O( µ) sin ∠CED = sin ∠BCE sin ∠BCE √ √ sin ∠BCE cos O( µ) + cos ∠BCE sin O( µ) = sin √ ∠BCE µ) O( √ = cos O( µ) + . tan ∠BCE

(29)

We can then compute the angle (λu · µ)2 + (λu · µ + O(µ2 ))2 − O(µ4 ) 2λu · µ(λu · µ + O(µ2 ))   2  2(λu · µ)2 + O(µ3 ) µ = arccos 1−O 2 2(λu · µ) λu µ = arccos(1 − O(µ)) √ = O( µ).

⟨∇2 f (x) · µu, ∇2 f (x) · µu + O(µ2 )⟩ = arccos

√ π−O( µ)

(30)

√ π−O( µ)

So there exists an angle ζ such that ∠BEC = + ζ, ∠BCE = − ζ. From 2 2    √  √ π − O( µ) O( µ) x2 sin + ζ = cos −ζ and cos x = 1 − + O(x4 ), 2 2 2 we get √ 1 O( µ) √ sin ∠BEC = 1 − ( − ζ)2 + O(( µ + ζ)4 ), 2 2 √ 1 O( µ) √ sin ∠BCE = 1 − ( + ζ)2 + O(( µ + ζ)4 ). 2 2 Therefore, √ √ O( µ) 1 − 21 ( 2 − ζ)2 + O(( µ + ζ)4 ) sin ∠BEC √ = √ O( µ) sin ∠BCE 1 − 12 ( 2 + ζ)2 + O(( µ + ζ)4 ) 2  √ 2 !  √ O( µ) O( µ) 1 =1+ +ζ − −ζ + O(µ1.5 , ζ 3 ) 2 2 2 √ = 1 + O( µζ) + O(µ1.5 , ζ 3 ).

On the other hand, by the Law of Sines, sin ∠BEC ∥∇2 f (x) · µu + O(µ2 )∥ = = 1 + O(µ). sin ∠BCE ∥∇2 f (x) · µu∥ √ Comparing two equations with the same left-hand side, we know that: ζ = O( µ), i.e. ∠BCE = √ π 2 + O( µ). Naturally, we have √ tan ∠BCE = O(1/ µ). Plugging in (29), we obtain sin ∠CED √ = 1 + O( µ). sin ∠BCE 19

Denote ⟨∇f (x), ∇f (x + µv)⟩ = α. Applying the Law of Cosines, we can compute α as cos α =

∥∇f (x)∥2 + ∥∇f (x + µv)∥2 − ∥∇2 f (x) · µv + O(µ2 )∥2 . 2 · ∥∇f (x)∥ · ∥∇f (x + µv)∥

So we have

√ α = O( µ).

(31)

∠BCD Denote ⟨∇f (x), u⟩ = θ. The other ratio sin sin ∠CDE in (28) can be computed as

√ √ sin θ cos O( µ) + cos θ sin O( µ) sin ∠BCD = sin ∠CDE sin θ cos α − cos θ sin α √ 1 + O(µ) + cot θ · O( µ) = √ 1 + O(µ) − cot θ · O( µ) √ = 1 + cot θ · O( µ).

(32)

Combining (31) and (32), we have sin ∠BEC sin ∠BCD √ η = · = (1 + O(µ))(1 + cot θ · O( µ)) = 1 + O sin ∠BCE sin ∠CDE ′

 √

µ tan θ

 .

3 Since ∥u − v∥ ≤ O( nξ ϵδκ ), we can further derive that

  nξ  η √ 3 = (1 + O(µ))(1 + cot θ − O · O( µ)) ∥H∥ ϵδκ   √ µ   = 1+O 3 tan θ − O( nξ ) ϵδκ  √    µ nξ3 =1+O +O . tan β ϵδκ

5.2

Proof of Theorem 7

Here we give a complete proof of our result on computing the ratio of the Hessian and the gradient norm. Proof of Theorem 7. By Lemma 5 we have (25). Given Line 4 and Line 6 in Algorithm 4, by Theorem 6 and Theorem 3 we have Ĥ −

H ≤ ξ3 , ∥H∥

ĝ1 =

g + k 1 v1 , ∥g + k1 v1 ∥

∥k1 v1 ∥ ≤ O(ξ1 )∥g∥.

(Analogous notation holds for ĝ2 with error k2 v2 .) Observe that   H H H H · µu − Ĥ · µv ≤ µ ·u− ·v + · v − Ĥ · v ∥H∥ ∥H∥ ∥H∥ ∥H∥   H ≤ µ ∥u − v∥ + − Ĥ . ∥H∥ 20

By (26), we have ∥u − v∥ ≤ O(ξ3 /ϵδκ ). Therefore, we have   H µnξ3 . · µu − Ĥ · µv ≤ O ∥H∥ ϵδκ and

 H · µu − ∥H∥ · Ĥ · µv ≤ ∥H∥ · O

µnξ3 ϵδκ

(33) 

(34)

.

√ The geometric calculation in Lemma 6 gives η/∥H∥ = 1 + O( µ/ tan θ), where θ = ∠(g, u). Due to (34), replacing the exact increment direction Hµu by the estimated increment direction ∥H∥Ĥµv 3 perturbs the relevant triangle angles by at most O( nξ ϵδκ ), hence we can write      √  √ µ µ η nξ 3   = 1 + O +O , = 1 + O ∥H∥ tan β ϵδκ tan θ − O( nξ3 )

(35)

ϵδκ

3 using θ ≥ β and µnξ ϵδκ sufficiently small. Furthermore, by Lemma 15 and Lemma 26,

By Taylor expansion,

q 1 − ⟨ĝ1 , ĝ2 ⟩2 η · ∥Ĥ · µv∥ κ̂ = q = . µ∥g + k1 v1 ∥ 2 µ 1 − ⟨ĝ2 , v⟩

(36)

∇f (x + µu) − ∇f (x) = H · µu + O(µ2 ),

(37)

so by (27) we have 1 ∥∇f (x + µu) − ∇f (x)∥ ∥H · µu + O(µ2 )∥ = µ ∥g∥ µ∥g∥   ∥Hu∥ µ = +O ∥g∥ ∥g∥   µ ∥Hu∥ +O = ∥g∥ ∥g∥   ∥H∥ ϵδκ + ξ3 + µ = +O . ∥g∥ ∥g∥

(38)

Combine (36) and (38), κ̂ −

1 ∥∇f (x + µu) − ∇f (x)∥ η∥Ĥµv∥ ∥Hµu + O(µ2 )∥ = − µ ∥g∥ µ∥g + k1 v1 ∥ µ∥g∥ ≤

η∥Ĥµv∥ ∥H∥∥Ĥµv∥ ∥H∥∥Ĥµv∥ ∥H∥∥Ĥµv∥ − + − µ∥g + k1 v1 ∥ µ∥g + k1 v1 ∥ µ∥g + k1 v1 ∥ µ∥g∥ | {z } | {z } T1

+

T2

∥H∥∥Ĥµv∥ ∥Hµu∥ ∥Hµu∥ ∥Hµu + O(µ2 )∥ − + − . µ∥g∥ µ∥g∥ µ∥g∥ µ∥g∥ {z } | {z } | T3

T4

(39) 21

Using ∥Ĥµv∥ ≤ µ∥Ĥ∥ ≤ µ(1 + ξ3 ) ≤ 2µ and ∥g + k1 v1 ∥ ≥ ∥g∥ − ∥k1 v1 ∥ ≥ (1 − O(ξ1 ))∥g∥ ≥ 12 ∥g∥, |η − ∥H∥| · ∥Ĥµv∥ C ≤ |η − ∥H∥|. µ∥g + k1 v1 ∥ ∥g∥   √ µ 3 By (35), |η − ∥H∥| ≤ ∥H∥ O( tan β ) + O( nξ ) . Hence ϵδκ T1 =

∥H∥ T1 ≤ ∥g∥ Using the inequality

1 1 ∥g+k1 v1 ∥ − ∥g∥

T2 ≤

  √    µ nξ3 O +O . tan β ϵδκ

(40)

(41)

∥k1 v1 ∥ 1 v1 ∥| = |∥g∥−∥g+k ∥g∥∥g+k1 v1 ∥ ≤ ∥g∥∥g+k1 v1 ∥ , we obtain

∥H∥∥Ĥµv∥ ∥k1 v1 ∥ ∥H∥ · ≤ · O(ξ1 ). µ ∥g∥∥g + k1 v1 ∥ ∥g∥

(42)

Moreover, ∥H∥∥Ĥµv∥ − ∥Hµu∥ = ∥H∥µλv − µ∥Hu∥ ≤ ∥H∥ · O(µ(ϵδκ + ξ3 )). which gives T3 ≤

∥H∥ · O(µ(ϵδκ + ξ3 )) ∥H∥ = · O(ϵδκ + ξ3 ). µ∥g∥ ∥g∥

(43)

As for T4 , from (37) we have ∥O(µ2 )∥ ≤ Cµ2 and thus   µ ∥O(µ2 )∥ ≤O . T4 ≤ µ∥g∥ ∥g∥

(44)

Plugging (41)–(44) into (39) yields κ̂ −

1 ∥∇f (x + µu) − ∇f (x)∥ ∥H∥ ≤ µ ∥g∥ ∥g∥

  √     µ nξ3 µ O + O( ) + O(ξ1 ) + O(ϵδκ + ξ3 ) +O , tan β ϵδκ ∥g∥ (45)

Combine (38) and (45), we have   √     µ nξ3 µ ∥H∥ ∥H∥ κ̂ − ≤ O + O( ) + O(ξ1 ) + O(ϵδκ + ξ3 ) + O ∥g∥ ∥g∥ tan β ϵδκ ∥g∥  √ µ nξ3 ξ1 + ξ3 ≤O + 2 + + δκ . ϵ ϵ δκ ϵ

(46)

By our choice of parameters in Algorithm 4 and using the Lipschitzness, the right-hand side of (46) is at most O(δκ ), which completes our proof.

6

Finding Stationary Points

In this section, we combine the tools developed in the previous sections to construct Algorithm 5, which outputs a list of points that includes at least one ϵ-second-order stationary point, a notion that generalizes and is strictly stronger than an ϵ-stationary point. Definition 1. A point x ∈ Rn is an ϵ-second-order stationary point of f if p
∥∇f (x)∥ ≤ ϵ and λmin ∇2 f (x) ≥ − L2 ϵ. 22

Algorithm 5: ComparisonTR(ϵ, ξ1 , ξ2 , ξ3 , T ): Finding stationary point Input: Initial point x0 , precision ϵ, tolerances (ξ1 , ξ2 , ξ3 ) = (1/100, 1/100, 1/100n), number of iterations T Output: p A list of points containing at least one ϵ-second order stationary point. ϵ/L2 . 1 r ← q ϵ 2 δκ ← c0 L2 . 3 for t = 0, . . . , T − 1 do 4 5 6 7 8 9 10 11

ĝ ← ComparisonGE(xt ; ξ1 ). Ĥ ← ComparisonHE(xt ; ξ3 ). κ̂ ← ComparisonRatio(xt ; δκ ). Compute a unit eigenvector v of Ĥ with eigenvalue of largest magnitude. x(1) ← xt − rĝ. Find l⋆ ∈ [−r, r] approximately minimizing f (xt + lv). x(2) ← xt + l⋆ v. Define mx (p) = ⟨ĝ, p⟩ +

12 13

κ̂ ⊤ p Ĥp, 2

∥p∥ ≤ r.

Let p⋆ ∈ arg min∥p∥≤r mx (p). x(3) ← x + p⋆ . xt+1 = arg min{f (x(1) ), f (x(2) ), f (x(3) ), f (xt )}.

14 return {x0 , . . . , xT }.

Theorem 8. Given f : Rn → R has L1 -Lipschitz gradient and L2 -Lipschitz Hessian. Denote x⋆ = arg min f (x). Given x0 ∈ Rn satisfying f (x0 ) − f (x⋆ ) ≤ ∆, with success probability at least 2/3, Algorithm 5 outputs a list of O(1/ϵ1.5 ) points, in which at least one point is an ϵ-second order √ ⋆ ∆ L2 n2 0 )−f (x ) stationary point, using T = O( f (xC ) iterations. The overall query complexity is Õ( ). ′ ϵ1.5 ϵ1.5 √ If in some steps ∇f (xt ) ≤ ϵ and λmin (∇2 f (xt )) ≥ − L2 ϵ, we have already visited an ϵ-second order stationary point. Otherwise, at any iterate xt with ∥∇f (xt )∥ ≥ ϵ, Algorithm 5 finds the next iterate xt+1 that decreases the function value by Ω(ϵ3/2 ). If after limited steps we cannot visit a stationary point, the function value will decrease more than f (x0 ) − f (x⋆ ), which is a contradiction. The iteration rule constructs three candidate steps (gradient step, trust-region step on a normalized quadratic model, and rank-one line-search step), and then selects the candidate with minimal function value by comparison. If the Hessian is rank-one or the eigenvector with the largest absolute eigenvalue has significant overlap with the gradient, line search performs well. Else if the Hessian has too small norm, or positive/negative definite, a normalized gradient descent step decreases Ω(ϵ1.5 ). Else, in all the other cases, trust-region step guarantees the decrease. Lemma 7. For any iterate √ xt in Algorithm 5 that satisfies: (i) ∥∇f (xt )∥ ≥ ϵ, or (ii) ∥∇f (xt )∥ < ϵ 2 and λmin (∇ f (xt )) ≤ − L2 ϵ, the next iterate xt+1 satisfies f (xt+1 ) ≤ f (xt ) − Cϵ1.5 for some absolute constant C > 0. 23

(47)

6.1

Descent of each candidate

For any iterate xt with ∥∇f (xt )∥ ≥ ϵ, let H := ∇2 f (xt ) and g := ∇f (xt ). Algorithm 5 builds three candidates within the ball: (i) normalized gradient descent along −ĝ with step size in [0, r], (ii) a line search along an eigenvector v corresponding to the eigenvalue with the largest magnitude on [−r, r], (iii) a trust-region minimizer of the approximate quadratic mxt (p) = ⟨ĝ, p⟩+ κ̂2 p⊤ Ĥp with ∥p∥ ≤ r. Lemma 8 (Descent of Candidate 1). Given f : Rn → R has L1 -Lipschitz gradient and L2 -Lipschitz Hessian. For any iterate xt ∈ Rn with ∥g∥ ≥ ϵ, let ĝ be the output of ComparisonGE(x; ξ1 ) with ĝ −

g ≤ ξ1 , ∥g∥

ξ1 ≤

1 . 100

(48)

Assume that at least one of the following two conditions holds: √ (i) ( small Hessian norm) ∥H∥ ≤ L2 ϵ; √ (ii) ( H is definite) H ⪰ 0 or H ⪯ 0, and ∥H∥ ≥ L2 ϵ; Then there exists an absolute constant C > 0 such that f (x(1) ) ≤ f (xt ) − C ϵ3/2 .

(49)

Proof. We use the cubic Taylor remainder bound for any p ∈ Rn , 1 L2 f (xt + p) ≤ f (xt ) + g⊤ p + p⊤ Hp + ∥p∥3 . 2 6

(50)

Denote p1 = −rĝ. We will bound the linear and quadratic terms for p1 and then absorb the remainder term. Let ḡ := g/∥g∥. Since ∥ĝ − ḡ∥ ≤ ξ1 and both are unit vectors, 1 1 ⟨ḡ, ĝ⟩ = 1 − ∥ḡ − ĝ∥2 ≥ 1 − ξ12 ≥ 1 − ξ1 . 2 2 Therefore

g⊤ p1 = −r ⟨g, ĝ⟩ = −r ∥g∥ ⟨ḡ, ĝ⟩ ≤ −(1 − ξ1 ) r ∥g∥.

(51)

Since ∥p1 ∥ = r and ∥H∥ denotes the spectral norm, 1 1 ⊤ 1 p Hp1 ≤ ∥H∥ ∥p1 ∥2 = ∥H∥ r2 . 2 1 2 2

(52)

Plugging (51) and (52) into (50) with p = p1 yields 1 L2 3 f (x(1) ) − f (xt ) ≤ −(1 − ξ1 )r∥g∥ + ∥H∥r2 + r . 2 6 p Using ∥g∥ ≥ ϵ and r = ϵ/L2 gives ϵ3/2 r∥g∥ ≥ √ , L2

r2 =

ϵ , L2

24

L2 3 1 ϵ3/2 r = √ . 6 6 L2

(53)

Hence (53) becomes  1 ϵ 1  ϵ3/2 √ + ∥H∥ · . f (x(1) ) − f (xt ) ≤ − (1 − ξ1 ) − 6 L2 L2 2

(54)

We now argue that in each of the three cases (i)–(iii), the positive quadratic term is controlled so that the right hand side√of (54) is less than −C ϵ3/2 after choosing c as a sufficiently small constant. In Case (i), ∥H∥ ≤ L2 ϵ. Then the quadratic term in (54) is bounded by 1 1p ϵ 1 ϵ3/2 ϵ ≤ L2 ϵ · = √ . ∥H∥ · 2 L2 2 L2 2 L2 Therefore

  1 1 ϵ3/2 √ . f (x(1) ) − f (xt ) ≤ − (1 − ξ1 ) − − 6 2 L2

So we have

f (x(1) ) ≤ f (xt ) − C ϵ3/2

for an absolute constant C > 0. In Case (ii), H ⪯ 0 or H ⪰ 0. By L2 -Lipschitzness, for any p, 1 L2 f (xt + p) ≤ f (xt ) + g⊤ p + p⊤ Hp + ∥p∥3 . 2 6 Apply this with p = p1 = −rĝ to obtain f (x(1) ) − f (xt ) ≤ rg⊤ ĝ +

r2 ⊤ L2 3 ĝ H ĝ + |r| . 2 6

(55)

Let ḡ := g/∥g∥. Since ∥ĝ − ḡ∥ ≤ ξ1 and both are unit vectors, 1 1 ⟨ḡ, ĝ⟩ = 1 − ∥ḡ − ĝ∥2 ≥ 1 − ξ12 ≥ 1 − ξ1 . 2 2 Therefore

g⊤ p1 = −r ⟨g, ĝ⟩ = −r∥g∥ ⟨ḡ, ĝ⟩ ≤ −(1 − ξ1 ) r∥g∥.

If H ⪯ 0, we have so

1 ⊤ p Hp1 ≤ 0. 2 1   L2 3 1 (1) ⊤ f (x ) − f (xt ) ≤ g p1 + |r| ≤ −ϵ (1 − ξ1 ) − r = −Ω(ϵ3/2 ). 6 6

Else if H ⪰ 0, consider the 1-dimension function f⋆ on the direction ĝ. By the Lipschitzness of H, for any l ∈ [0, r], we have f⋆′′ (xt − lĝ) ≥ 0. It tells us that

f⋆′ (xt − lĝ) ≥ ϵ.

So the function value descent f (x(1) ) − f (xt ) ≤ −rϵ = −Ω(ϵ3/2 ).

25

Lemma 9 (Descent of Candidate 2). Let x(2) = xt + l⋆ v where l⋆ ∈ [−r, r] minimizes f (xt + tv) over [−r, r]. Then either ∥g∥ ≤ ϵ, or f (x(2) ) − f (xt ) ≤ −c2 ϵ3/2

(56)

for an absolute constant c2 > 0. Proof. Consider ϕ(t) = f (xt + lv) and expand around l = 0, we have 1 ϕ(l) = ϕ(0) + l ⟨g, v⟩ + l2 v⊤ Hv + O(L2 |l|3 ). 2 For | ⟨g, v⟩ | ≥ 21 ∥g∥, choosing l = −r sign(⟨g, v⟩) yields ϕ(l) − ϕ(0) ≤ −r · 21 ∥g∥ + O(L2 r3 ) ≤ −Ω(ϵ3/2 ) √ since ∥g∥ ≥ ϵ and L2 r3 = ϵ3/2 / L2 . If instead | ⟨g, v⟩ | < 12 ∥g∥, then g has a component orthogonal √

to v of magnitude at least 23 ∥g∥, and Candidate 1 already provides a decrease of order r∥g∥ up to quadratic/cubic corrections; therefore the best-of-all selection ensures the stated decrease. Lemma 10 (Descent of Candidate 3). Let p⋆ minimize the normalized quadratic model mxt (p) = ⟨g, p⟩ +

κ̂ ⊤ p Ĥ p 2

subject to ∥p∥ ≤ r,

Then either or

∥∇mxt (p⋆ )∥ ≤ ϵ

(57)

1 f (x(3) ) − f (xt ) = mxt (p⋆ ) − mxt (0) ≤ − r∥∇mxt (p⋆ )∥. 2

(58)

In particular, in the non-stationary case (58) implies mxt (p⋆ ) − mxt (0) ≤ −Ω(rϵ) = −Ω(ϵ3/2 ). Proof. This is a direct corollary of Lemma 28, taking q = mxt , b = g and A = κ̂Ĥ. √ Next, we prove that Algorithm 5 can go further to visit an (ϵ, ϵ)-second order stationary point. Our intuition is to construct a negative curvature candidate and prove that the trust region candidate is better than or equal to it. Lemma 11. For any iterate xt ∈ Rn , suppose that vmin is the unit eigenvector corresponding to the minimum eigenvalue of H, denoted√as λmin . If ⟨vmin , g⟩ ≥ 0, set vmin = −vmin . Let x(4) = xt + rvmin . If ∥g∥ ≤ ϵ and λmin ≤ − L2 ϵ, then f (x(4) ) − f (xt ) ≤ −c4 ϵ3/2 for an absolute constant c4 > 0. Proof. We use the cubic Taylor remainder bound for any p ∈ Rn , 1 L2 ∥p∥3 . f (xt + p) ≤ f (xt ) + g⊤ p + p⊤ Hp + 2 6

26

(59)

Plugging in p4 = rvmin we have 1 L2 3 ⊤ f (x(4) ) ≤ f (xt ) + rg⊤ vmin + r2 vmin Hvmin + r 2 6 r ϵ ⊤ 1 ϵ ⊤ L2 ϵ3/2 = f (xt ) + g vmin + vmin Hvmin + · L2 2 L2 6 L3/2 2 s   1 1 ϵ3 − + = f (xt ) − c4 ϵ3/2 . ≤ f (xt ) + L2 2 6

Lemma √ 12. For any iterate xt ∈ Rn and let g = ∇f (xt ) and H = ∇2 f (xt ). If ∥g∥ ≤ ϵ and λmin ≤ − L2 ϵ, then f (x(3) ) ≤ f (x(4) ) + c5 ϵ3/2 (60) for an absolute constant c5 > 0. Proof. Denote x(3) = xt + p3 , x(4) = xt + p4 . For our choice of x(3) , we have mxt (p3 ) ≤ mxt (p4 ), which is equivalent to Txt (xt + p3 ) ≤ Txt (xt + p4 ).

(61)

By Lemma 27 we have |f (xt + p3 ) − Txt (xt + p3 )| ≤

L2 L2 3 ∥p3 ∥3 = r , 6 6

(62)

|f (xt + p4 ) − Txt (xt + p4 )| ≤

L2 L2 3 ∥p4 ∥3 = r . 6 6

(63)

Combining (61), (62) and (63), we can conclude (60).

6.2

Proof of function value descent guaranty

Here, we prove the theoretical guarantee of our Algorithm 5 for a large descent when the current point is not a second-order stationary point. Proof of Lemma 7. Assume ∥g(xt )∥ ≥ ϵ. Let xt+1 = arg min{f (x(1) ), f (x(2) ), f (x(3) )} in Algorithm 5. By Lemma 27, for each candidate pi = x(i) − xt with ∥pi ∥ ≤ r,
L2 3 f (x(i) ) − f (xt ) ≤ Txt (xt + pi ) − Txt (xt ) + r . 3

(64)

We show that at least ) ≤ −c′ rϵ for a constant c′ > 0; √ one candidate has Txt (xt + pi ) − Txt (xt(i) L2 3 1 3/2 then since 3 r = 3 ϵ / L2 , for c sufficiently small, we obtain f (x ) − f (xt ) ≤ −Ω(ϵ3/2 ) and thus the same for xt+1 . If the Hessian is rank-one or ⟨g, v⟩ ≥ 21 ∥g∥, x(2) ensures the Ω(ϵ3/2 ) decrease by Lemma 9. Else, for cases in Lemma 8, Txt (xt + p1 ) − Txt (xt ) ≤ − 27

99 1 r∥g∥ + r2 ∥H∥. 100 2

Using the ratio ∥H∥/∥g∥ ≈ κ̂, we have 1 2 r ∥H∥ ≤ cr∥g∥ 2 at the scale r =

p ϵ/L2 , so the linear term dominates for small ϵ. Hence Txt (xt + p1 ) − Txt (xt ) ≤ −Ω(r∥g∥) ≤ −Ω(rϵ) = −Ω(ϵ3/2 ).

Then (64) yields

f (x(1) ) − f (xt ) ≤ −Ω(ϵ3/2 ). √ In the case that H is not rank-one, ∥H∥ ≥ L2 ϵ, and H is not positive or negative definite. Denote the eigenvector of H that corresponding to the largest magnitude eigenvalue as v and assume that ⟨v, g⟩ ≤ 21 ∥g∥. If (57) fails, then by Lemma 10, we have 1 mxt (p3 ) − mxt (0) ≤ − r∥∇mxt (p3 )∥. 2 Since the failure of (57) means ∥∇mxt (p3 )∥ > ϵ, we get 1 mxt (p3 ) − mxt (0) ≤ − rϵ = −Ω(ϵ3/2 ). 2 By the constant-factor accuracy of κ̂ and Ĥ, the model mxt is a constant-factor approximation of the quadratic model Txt when ∥p∥ ≤ r. Thus Txt (xt + p3 ) − Txt (xt ) ≤ −Ω(ϵ3/2 ), and (64) gives

f (x(3) − f (xt ) ≤ −Ω(ϵ3/2 ). √ Assume ∥g∥ ≤ ϵ and λmin (H) ≤ − L2 ϵ. By Lemma 12 we have f (x3 ) − f (xt ) ≤ −Ω(ϵ3/2 ).

In all cases, at least one candidate x(i) satisfies f (x(i) ) ≤ f (xt ) − Cϵ3/2 . Since xt+1 = arg min{f (x(1) ), f (x(2) ), f (x(3) )}, we have f (xt+1 ) ≤ f (x(i) ) for that index i, proving (47).

6.3

Proof of Theorem 8

Here we give a complete proof of our result main result on finding the stationary points. Proof of Theorem 8. We use contradiction to prove the theorem. If the statement is not true and any iterate xt is not an ϵ-stationary point, by Lemma 7 we have E[f (xt+1 ) − f (xt )] ≤ −

2C 1.5 ϵ , 3

where C is the constant in Lemma 7, since we have f (xt+1 ) − f (xt ) ≤ −Cϵ1.5 if Algorithm 2 succeeds, which happens with probability at least 2/3 by Theorem 6, and we have f (xt+1 ) − f (xt ) ≤ 0 otherwise. Sum from t = 0 to T − 1 we get E[f (xT ) − f (x0 )] =

T −1 X

E[f (xt+1 ) − f (xt )] = −

t=0

T −1 X t=0

28

C ′ ϵ1.5 = −T C ′ ϵ1.5 .

By Markov’s Inequality,   a − E[f (x0 ) − f (xT )] 2 1 ′ 1.5 ≥ = . Pr f (x0 ) − f (xT ) ≥ T C ϵ 1 ′ 1.5 2 3 a − 2TC ϵ Therefore, with probability at least 2/3, the total decrease of function value is larger than f (x0 ) − f (x⋆ ), which leads to a contradiction. Consequently, in T steps, we have visited at least an ϵ-stationary point. Note that we have used Comparison-GE in Theorem 3 as a subroutine to estimate the direction of the gradient. The query complexity of Comparison-GE is O(n log 1ϵ ), which incurs an O(n) overhead – intuitively, classical algorithms by comparisons are limited by the fact that we need Ω(n) comparisons to explore an n-dimensional space. An idea that can significantly reduce the query complexity of our algorithm is by replacing Comparison-GE with Comparison-QGE in Theorem 4, a quantum algorithm with log n dependence for gradient estimation. This implies the following corollary. Corollary 3. There exists a quantum algorithm that visits an ϵ-second order stationary point using √
L2 n O( ∆ ϵ1.5 log nLϵ1 L2 ) queries to a quantum comparison oracle (2).

Acknowledgements We thank the anonymous reviewers for their constructive feedback. HW, XT, YZ, and TL were supported by the National Natural Science Foundation of China (Grant Number 62372006).

References [1] Deeksha Adil, Brian Bullins, Aaron Sidford, and Chenyi Zhang, Balancing gradient and Hessian queries in non-convex optimization, The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025, arXiv:2510.20786 1 [2] Naman Agarwal, Zeyuan Allen-Zhu, Brian Bullins, Elad Hazan, and Tengyu Ma, Finding approximate local minima faster than gradient descent, Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, pp. 1195–1199, 2017, arXiv:1611.01146 1 [3] Zeyuan Allen-Zhu, How to make the gradients small stochastically: Even faster convex and nonconvex SGD, Advances in Neural Information Processing Systems, vol. 31, 2018, arXiv:1801.02982 1 [4] Charles Audet and John E. Dennis Jr, Mesh adaptive direct search algorithms for constrained optimization, SIAM Journal on Optimization 17 (2006), no. 1, 188–217. 2 [5] El Houcine Bergou, Eduard Gorbunov, and Peter Richtárik, Stochastic three points method for unconstrained smooth minimization, SIAM Journal on Optimization 30 (2020), no. 4, 2726–2749, arXiv:1902.03591 2, 3, 5 [6] Jeremy Bernstein, Yu-Xiang Wang, Kamyar Azizzadenesheli, and Animashree Anandkumar, signSGD: Compressed optimisation for non-convex problems, International Conference on Machine Learning, pp. 560–569, PMLR, 2018, arXiv:1802.04434 1

29

[7] Ernesto G. Birgin, J.L. Gardenghi, José Mario Martínez, Sandra Augusta Santos, and Ph.L. Toint, Worst-case evaluation complexity for unconstrained nonlinear optimization using high-order regularized models, Mathematical Programming 163 (2017), no. 1, 359–368, arXiv:1709.07180 1 [8] Ralph Allan Bradley and Milton E. Terry, Rank analysis of incomplete block designs: I. the method of paired comparisons, Biometrika 39 (1952), no. 3/4, 324–345. 5 [9] Yair Carmon, John C. Duchi, Oliver Hinder, and Aaron Sidford, Lower bounds for finding stationary points I, Mathematical Programming 184 (2020), no. 1, 71–120, arXiv:1710.11606 1, 3 [10] Xiaoyu Chen, Han Zhong, Zhuoran Yang, Zhaoran Wang, and Liwei Wang, Human-in-the-loop: Provably efficient preference-based reinforcement learning with general function approximation, International Conference on Machine Learning, pp. 3773–3793, PMLR, 2022, arXiv:2205.11140 2 [11] Paul F. Christiano, Jan Leike, Tom Brown, Miljan Martic, Shane Legg, and Dario Amodei, Deep reinforcement learning from human preferences, Advances in Neural Information Processing Systems, vol. 30, 2017, arXiv:1706.03741 1 [12] Nikita Doikov and Geovani Nunes Grapiglia, First and zeroth-order implementations of the regularized Newton method with lazy approximated Hessians, Journal of Scientific Computing 103 (2025), no. 1, 32, arXiv:2309.02412 5 [13] Flint Xiaofeng Fan, Cheston Tan, Yew-Soon Ong, Roger Wattenhofer, and Wei-Tsang Ooi, FedRLHF: A convergence-guaranteed federated framework for privacy-preserving and personalized RLHF, Proceedings of the 24th International Conference on Autonomous Agents and Multiagent Systems, pp. 713–721, 2025, arXiv:2412.15538 1 [14] Mudit Gaur, Amrit Singh Bedi, Raghu Pasupathy, and Vaneet Aggarwal, On the global convergence of online RLHF with neural parametrization, 2024, arXiv:2410.15610 1 [15] Rong Ge, Jason D. Lee, and Tengyu Ma, Matrix completion has no spurious local minimum, Advances in Neural Information Processing Systems, vol. 29, 2016, arXiv:1605.07272 1, 3 [16] Rong Ge and Tengyu Ma, Decomposing overcomplete 3rd order tensors using sum-of-squares algorithms, Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques (2015), 829, arXiv:1504.05287 3 [17] Rong Ge and Tengyu Ma, On the optimization landscape of tensor decompositions, Advances in Neural Information Processing Systems, vol. 30, 2017, arXiv:1706.05598 1 [18] Daniel Golovin, John Karro, Greg Kochanski, Chansoo Lee, Xingyou Song, and Qiuyi Zhang, Gradientless descent: High-dimensional zeroth-order optimization, International Conference on Learning Representations, 2020, arXiv:1911.06317 2 [19] Eduard Gorbunov, Adel Bibi, Ozan Sener, El Houcine Bergou, and Peter Richtarik, A stochastic derivative free optimization method with momentum, International Conference on Learning Representations, 2020, arXiv:1905.13278 2, 3, 5

30

[20] Kevin G. Jamieson, Robert Nowak, and Ben Recht, Query complexity of derivative-free optimization, Advances in Neural Information Processing Systems, vol. 25, 2012, arXiv:1209.2434 2 [21] Mustafa O. Karabag, Cyrus Neary, and Ufuk Topcu, Smooth convex optimization using subzeroth-order oracles, Proceedings of the AAAI Conference on Artificial Intelligence 35 (2021), no. 5, 3815–3822, arXiv:2103.00667 2 [22] Tamara G. Kolda, Robert Michael Lewis, and Virginia Torczon, Optimization by direct search: New perspectives on some classical and modern methods, SIAM Review 45 (2003), no. 3, 385–482. 2 [23] Jeffrey Larson, Matt Menickelly, and Stefan M. Wild, Derivative-free optimization methods, Acta Numerica 28 (2019), 287–404, arXiv:904.11585 2 [24] Xiuxian Li, Kuo-Yi Lin, Li Li, Yiguang Hong, and Jie Chen, On faster convergence of scaled sign gradient descent, IEEE Transactions on Industrial Informatics, IEEE, 2023, arXiv:2109.01806 1 [25] Chengchang Liu, Zongqi Wan, Jialin Zhang, Xiaoming Sun, and John C.S. Lui, Quantum speedups for minimax optimization and beyond, The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2026. 5 [26] Sijia Liu, Pin-Yu Chen, Xiangyi Chen, and Mingyi Hong, signSGD via zeroth-order oracle, International Conference on Learning Representations, 2019. 1 [27] Aleksandr Lobanov, Alexander Gasnikov, and Andrei Krasnov, Acceleration exists! Optimization problems when oracle can only compare objective function values, Advances in Neural Information Processing Systems, vol. 37, 2024, arXiv:2402.09014 2 [28] Po-Ling Loh and Martin J. Wainwright, Regularized M-estimators with nonconvexity: Statistical and algorithmic theory for local optima, The Journal of Machine Learning Research 16 (2015), no. 1, 559–616, arXiv:1305.2436 1, 3 [29] Ken I. M. McKinnon, Convergence of the Nelder–Mead simplex method to a nonstationary point, SIAM Journal on Optimization 9 (1998), no. 1, 148–158. 2 [30] John A. Nelder and Roger Mead, A simplex method for function minimization, The Computer Journal 7 (1965), no. 4, 308–313. 2 [31] Yurii Nesterov, How to make the gradients small, Optima. Mathematical Optimization Society Newsletter 88 (2012), 10–11. 1 [32] Yurii Nesterov, Introductory lectures on convex optimization: A basic course, vol. 87, Springer Science & Business Media, 2013. 1 [33] Yurii Nesterov and Boris Polyak, Cubic regularization of newton method and its global performance, Math. Program. 108 (2006), 177–205. 40 [34] Yurii Nesterov and Boris T. Polyak, Cubic regularization of Newton method and its global performance, Mathematical Programming 108 (2006), no. 1, 177–205. 3

31

[35] Ellen Novoseller, Yibing Wei, Yanan Sui, Yisong Yue, and Joel Burdick, Dueling posterior sampling for preference-based reinforcement learning, Conference on Uncertainty in Artificial Intelligence, pp. 1029–1038, PMLR, 2020, arXiv:1908.01289 2 [36] Long Ouyang, Jeffrey Wu, Xu Jiang, Diogo Almeida, Carroll Wainwright, Pamela Mishkin, Chong Zhang, Sandhini Agarwal, Katarina Slama, Alex Ray, et al., Training language models to follow instructions with human feedback, Advances in Neural Information Processing Systems 35 (2022), 27730–27744, arXiv:2203.02155 2 [37] Yuxuan Ren, Abhishek Roy, and Shiqian Ma, Riemannian dueling optimization, 2026, arXiv:2603.00023 2 [38] Aadirupa Saha, Vitaly Feldman, Yishay Mansour, and Tomer Koren, Faster convergence with multiway preferences, Proceedings of the 27th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, vol. 238, pp. 433–441, PMLR, 2024, arXiv:2312.11788 2 [39] Aadirupa Saha, Tomer Koren, and Yishay Mansour, Dueling convex optimization, International Conference on Machine Learning, pp. 9245–9254, PMLR, 2021, arXiv:2210.02562 2 [40] Aadirupa Saha, Tomer Koren, and Yishay Mansour, Dueling convex optimization with general preferences, International Conference on Machine Learning, 2025, arXiv:2210.02562 2 [41] Aadirupa Saha, Aldo Pacchiano, and Jonathan Lee, Dueling RL: Reinforcement learning with trajectory preferences, International Conference on Artificial Intelligence and Statistics, pp. 6263–6289, PMLR, 2023, arXiv:2111.04850 2 [42] Katya Scheinberg and Zikai Xiong, Function-free optimization via comparison oracles, 2026, arXiv:2604.26867 2 [43] Zhiwei Tang, Dmitry Rybin, and Tsung-Hui Chang, Zeroth-order optimization meets human feedback: Provable learning via ranking oracles, 2023, arXiv:2303.03751 2 [44] Xiwen Tao, Chenyi Zhang, Helin Wang, Yexin Zhang, and Tongyang Li, Optimal classical and quantum algorithms for gradient testing and estimation by comparisons, 2026, To appear in the Forty-Third International Conference on Machine Learning, arXiv:2405.11454v2 3, 4, 5, 36 [45] Yuanhao Wang, Qinghua Liu, and Chi Jin, Is RLHF more difficult than standard RL? a theoretical perspective, Thirty-seventh Conference on Neural Information Processing Systems, 2023, arXiv:2306.14111 2 [46] Hermann Von Weyl, Über beschränkte quadratische formen, deren differenz vollstetig ist, Rendiconti del Circolo Matematico di Palermo (1884-1940) 27 (1909), 373–392. 39 [47] Yichong Xu, Ruosong Wang, Lin Yang, Aarti Singh, and Artur Dubrawski, Preference-based reinforcement learning with finite-time guarantees, Advances in Neural Information Processing Systems 33 (2020), 18784–18794, arXiv:2006.08910 2 [48] Yi Yu, Tengyao Wang, and Richard J. Samworth, A useful variant of the Davis-Kahan theorem for statisticians, Biometrika 102 (2015), no. 2, 315–323, arXiv:1405.0680 38 [49] Banghua Zhu, Jiantao Jiao, and Michael Jordan, Principled reinforcement learning with human feedback from pairwise or k-wise comparisons, ICLR 2023 Workshop on Mathematical and Empirical Understanding of Foundation Models, 2023, arXiv:2301.11270 2 32

[50] Danny C. Sorensen, Newton’s method with a model trust region modification, SIAM Journal on Numerical Analysis 19 (1982), no. 2, 409–426. 4

33

A

Auxiliary Lemmas

In this appendix, we collect all auxiliary lemmas needed for our proofs.

A.1

Distance between normalized vectors

Lemma 13. If v, v′ ∈ Rn are two vectors such that ∥v∥ ≥ γ and ∥v − v′ ∥ ≤ τ , we have v′ 2τ v − ′ ≤ . ∥v∥ ∥v ∥ γ Proof. By the triangle inequality, we have v v′ v v′ v′ v′ + − ′ ≤ − − ′ ∥v∥ ∥v ∥ ∥v∥ ∥v∥ ∥v∥ ∥v ∥ ′ ∥v − v ∥ |∥v∥ − ∥v′ ∥|∥v′ ∥ = + ∥v∥ ∥v∥∥v′ ∥ τ τ 2τ ≤ + = . γ γ γ

Lemma 14. If v1 , v2 ∈ Rn are two vectors such that ∥v1 ∥, ∥v2 ∥ ≥ γ, and v1′ , v2′ ∈ Rn are another two vectors such that ∥v1 − v1′ ∥, ∥v2 − v2′ ∥ ≤ τ where 0 < τ < γ, we have    ′  v1 6τ v1 v2 v2′ − ≤ , , ′ . ′ ∥v1 ∥ ∥v2 ∥ ∥v1 ∥ ∥v2 ∥ γ Proof. By the triangle inequality, we have    ′     ′  v1 v2 v1 v2′ v1 v2 v1 v2′ , − , ≤ , − , ∥v1 ∥ ∥v2 ∥ ∥v1′ ∥ ∥v2′ ∥ ∥v1 ∥ ∥v2 ∥ ∥v1 ∥ ∥v2 ∥  ′   ′  v1 v′ v1 v2′ + , 2 − , . ∥v1 ∥ ∥v2 ∥ ∥v1′ ∥ ∥v2′ ∥ On the one hand, by the triangle inequality and the Cauchy-Schwarz inequality,   ′   v2 v1 v2′ 1 v1 , − , ≤ ( ⟨v1 , v2 ⟩ − ⟨v1 , v2′ ⟩ + ⟨v1 , v2′ ⟩ − ⟨v1′ , v2′ ⟩⟩ ) ∥v1 ∥ ∥v2 ∥ ∥v1 ∥ ∥v2 ∥ ∥v1 ∥∥v2 ∥ ∥v2 − v2′ ∥ ∥v1 − v1′ ∥∥v2′ ∥ ≤ + ∥v2 ∥ ∥v1 ∥∥v2 ∥ τ τ (γ + τ ) ≤ + . γ γ2 On the other hand, by the Cauchy-Schwarz inequality, |⟨v1′ , v2′ ⟩| ≤ ∥v1′ ∥∥v2′ ∥, and hence  ′   ′  1 v1 v2′ v1 v2′ 1 , − , ′ = |⟨v1′ , v2′ ⟩| − ′ ′ ∥v1 ∥ ∥v2 ∥ ∥v1 ∥ ∥v2 ∥ ∥v1 ∥∥v2 ∥ ∥v1 ∥∥v2′ ∥ ∥v1′ ∥∥v2′ ∥ ≤ −1 ∥v1 ∥∥v2 ∥   γ+τ 2 ≤ − 1. γ 34

In all, due to τ < γ,

v1 v2 , ∥v1 ∥ ∥v2 ∥

−

v1′ v2′ , ∥v1′ ∥ ∥v2′ ∥

τ τ (γ + τ ) ≤ + + γ γ2



γ+τ γ

2 −1=

2τ (2γ + τ ) 6τ ≤ . 2 γ γ

Lemma 15. For any nonzero vectors v, g ∈ Rn , v D E2 u u v+g v u 1 − ∥v+g∥ , ∥v∥ ∥v − g∥ u . D E2 = t ∥v + g∥ v−g v 1 − ∥v−g∥ , ∥v∥ Proof. We have v+g v 2 v 2 1 − ⟨ ∥v+g∥ , ∥v∥ ⟩ ∥v + g∥2 ∥v + g∥2 − ⟨v + g, ∥v∥ ⟩ = · v 2 v−g v 2 ∥v − g∥2 2 ∥v − g∥ − ⟨v − g, ∥v∥ ⟩ 1 − ⟨ ∥v−g∥ , ∥v∥ ⟩

=

2 ⟨v + g, v + g⟩ − (∥v∥ + ⟨v,g⟩ ∥v∥ ) 2 ⟨v − g, v − g⟩ − (∥v∥ − ⟨v,g⟩ ∥v∥ ) 2

=

A.2

) ∥v∥2 + ∥g∥2 + 2⟨v, g⟩ − (∥v∥2 + 2⟨v, g⟩ + ⟨v,g⟩ ∥v∥2 2

∥v∥2 + ∥g∥2 − 2⟨v, g⟩ − (∥v∥2 − 2⟨v, g⟩ + ⟨v,g⟩ ) ∥v∥2

= 1.

Distance between normalized matrices

Lemma 16. Let A, B ∈ Rn×n be nonzero and assume ∥A − B∥ ≤ ρ < 12 ∥B∥. Then A B 4ρ − ≤ . ∥A∥ ∥B∥ ∥B∥ Proof. Write

A B A−B − = +B ∥A∥ ∥B∥ ∥A∥



1 1 − ∥A∥ ∥B∥

 .

Since ∥A∥ ≥ ∥B∥ − ∥A − B∥ ≥ ∥B∥ − ρ ≥ 12 ∥B∥, we have ∥A∥−1 ≤ 2∥B∥−1 . Also, 1 1 |∥A∥ − ∥B∥| ρ 2ρ = − ≤ 1 = . ∥A∥ ∥B∥ ∥A∥∥B∥ ∥B∥2 ( 2 ∥B∥)∥B∥ Therefore, A B ∥A − B∥ 1 1 ρ 2ρ 4ρ ≤ ≤ 1 − + ∥B∥ − + ∥B∥ = . 2 ∥A∥ ∥B∥ ∥A∥ ∥A∥ ∥B∥ ∥B∥ ∥B∥ ∥B∥ 2

35

A.3

Column norm bound implies spectral norm bound

Lemma 17. Let E ∈ Rn×n . If ∥Eei ∥ ≤ εc for all i = 1, . . . , n, then √ ∥E∥ ≤ εc n. Proof. For any v ∈ Rn with ∥v∥ = 1, Ev =

n X

vi (Eei ),

i=1

so ∥Ev∥ ≤

n X

|vi |∥Eei ∥ ≤ εc

i=1

n X

n 1/2 √ √ X = εc n. vi2 |vi | ≤ εc n

i=1

√

i=1

Taking the supremum over ∥v∥ = 1 yields ∥E∥ ≤ εc n.

A.4

Inner Product Concentration for Random Vectors on Sphere

Here we give a lemma of the inner product concentration for random vectors on sphere proved in [44], stated below: Lemma 18 (Lemma 3 of [44]). Let n ≥ 5. For any x ∈ Rn , x ̸= 0, and any constant c > 0, there exists constant p1 and p2 which is independent of n, such that √ p1 ≤ Pr [|⟨y, x⟩| ≤ ∥x∥/(c n)] ≤ p2 . y∼Sn

where y is chosen from Sn uniformly at random. In particular, we have the inequality   24 Pr |⟨y, x⟩| ≤ ∥x∥ · √ ≥ 3/5. y∼Sn 25 n

B

(65)

Basic Lemmas for Robust Hessian Estimation

In this section, we provide some supplementary lemmas for the guaranties claimed in Section 4. Throughout, we assume that ∥He1 ∥ = maxi ∥Hei ∥ > 0. We can make this assumption since we can find a column with maximum norm after obtaining all the column norm ratios, and it suffices to denote this column as index 1.

B.1

Notation

Let H ∈ Rn×n be symmetric, H ̸= 0, with columns H = (h1 , . . . , hn ), hi := Hei . Define ri :=

∥hi ∥ ∈ [0, 1], ∥h1 ∥

r1 = 1.

When hi = ̸ 0, define the unit column direction ui := hi /∥hi ∥. When h1 + hi = ̸ 0, define the unit sum direction u1i := (h1 + hi )/∥h1 + hi ∥. Define the scale-free target matrix H ⋆ := [r1 u1 , . . . , rn un ] ∈ Rn×n . Then H = ∥h1 ∥H ⋆ and hence H/∥H∥ = H ⋆ /∥H ⋆ ∥. 36

(66)

Lemma 19. The matrix H ⋆ be in (66) satisfies ∥H ⋆ ∥ ≥ 1. Proof. Because r1 = 1 and u1 is a unit vector, the first column of H ⋆ equals H ⋆ e1 = u1 and thus ∥H ⋆ ∥ ≥ ∥H ⋆ e1 ∥ = ∥u1 ∥ = 1.

B.2

Error bound of column vector estimation

Here we give some lemmas for proving the fact that the small perturbation hi + ρhj can only affect the estimation of hi under a bounded error when ρ is bounded. Lemma 20. Our choice of gi and g1i in Algorithm 2, Line 7, Line 10, Line 14 and Line 17 satisfies ∥gi − ui ∥ ≤ η and ∥g1i − u1i ∥ ≤ η. Proof. We analyze gi here, g1i case is similar. Assume that v is an eigenvector of H corresponding to minimum eigenvalue, and |⟨et , v⟩| ≥ √1n . Fix an index i. Denote yi0 = ei , yi1 = ei + σet , yi2 = ei − σet . Observe that for any j ̸= k ∈ {0, 1, 2}, σ |⟨yij − yik , v⟩| ≥ √ . n This implies at least two of three queries satisfies σ |⟨yij , v⟩| ≥ √ . 2 n

(67)

σ For each index j satisfies (67), by taking γy = 2√ , γx = ϵ, accuracy δ̂ = η4 in Algorithm 1, we n

know at least two queries in the setting {yi0 , yi1 , yi2 } satisfies the condition |⟨yij , v⟩| ≥ γy . Therefore, according to Theorem 5, at least two of three inequalities hold: gi0 −

hi η ≤ , ∥hi ∥ 4

gi+ −

hi + σht η ≤ , ∥hi + σht ∥ 4

gi− −

hi − σht η ≤ . ∥hi − σht ∥ 4

In non-degenerate case, our choice (v1 , v2 ) ← arg maxv1 ̸=v2 ∈{g0 ,g+ ,g− } ⟨v1 , v2 ⟩ tells v1 and v2 are i

i

i

two success queries to Algorithm 1, take v1 = gi+ , v2 = gi− for example, other case is the same, from Lemma 13 we have v1 + v 2 hi + σht hi − σht hi + σht hi hi − σht − ≤ gi+ − + gi− − + − ∥v1 + v2 ∥ ∥hi ∥ ∥hi + σht ∥ ∥hi − σht ∥ ∥hi + σht ∥ ∥hi − σht ∥ η η 2η ≤ + + = η, 4 4 4 where the last inequality holds because ⟨gi+ , gi− ⟩ ≥ 1 − τα tell us ∥hi ∥/∥σht ∥ ≥ η4 . In the degenerate case ⟨gi+ , gi− ⟩ ≤ 1 − τα , we have ∥hi ∥/∥σht ∥ ≤ η4 , i.e., ∥hi ∥ ≤

η √ . 32 n

In this case, taking gi = 0 guarantees error on this column is O 37



√η n



.

The following lemma is to bound the error between the normalized vector and the perturbed one. Lemma 21. Let a, c ∈ Rn with a ̸= 0. For ρ ∈ (0, 1/2) define v :=

a + ρc , ∥a + ρc∥

u :=

a . ∥a∥

If ∥c∥ ≤ ∥h1 ∥ and ∥a∥ ≥ rmin ∥h1 ∥, then ∥v − u∥ ≤

4 rmin

ρ.

Proof. By triangle inequality, ∥a + ρc∥ ≥ ∥a∥ − ρ∥c∥ ≥ ∥a∥ − ρ∥h1 ∥. Using ∥a∥ ≥ rmin ∥h1 ∥ and ρ ≤ rmin /2, we get ∥a + ρc∥ ≥ 21 ∥a∥. Therefore,   a + ρc ρ a 1 1 v−u= + − =a − c. ∥a + ρc∥ ∥a∥ ∥a + ρc∥ ∥a∥ ∥a + ρc∥ Hence ∥v − u∥ ≤ ∥a∥ Also

1 ρ 1 − + ∥c∥. ∥a + ρc∥ ∥a∥ ∥a + ρc∥

1 1 ρ∥c∥ |∥a + ρc∥ − ∥a∥| − ≤ , = ∥a + ρc∥ ∥a∥ ∥a + ρc∥∥a∥ ∥a + ρc∥∥a∥

by which we can conclude that ∥v − u∥ ≤ ∥a∥

C

ρ∥c∥ ρ 2ρ∥c∥ 2ρ∥h1 ∥ 4ρ + ∥c∥ = ≤ 1 = . ∥a + ρc∥∥a∥ ∥a + ρc∥ ∥a + ρc∥ r ∥a∥ min 2

Basic Lemmas for Estimating the Hessian-to-Gradient Ratio

We provide some supplementary lemmas of the claims in Section 5 in this appendix. We first introduce the core lemma in our proof: e ∈ Rd×d be two symmetric Lemma 22 (Davis-Kahan Theorem, see e.g., Theorem 1 of [48]). Let A, A e ≤ ξ for some ξ > 0. For any a < b, denote S = {v1 , . . . , vk } and matrices satisfying ∥A − A∥ e e associated with eigenvalues S = {ṽ1 , . . . , ṽk } as the set of normalized eigenvectors of A and A contained in the interval [a, b] and [a − ξ, b + ξ] respectively, and denote V := span(S),

e Ve := span(S).

Then, if the remaining eigenvalues of A lie outside the interval [a − γ, b + γ], we have k = k̃ and   ξ sin Θ(V, Ve ) ≤ , γ where

  sin Θ(V, Ve ) := diag(sin θ1 (V, Ve ), . . . , sin θk (V, Ve ))⊤ .

38

Lemma 23 (Weyl’s Theorem, see [46]). Let A, B be Hermitian on inner product space V with dimension n, with spectrum ordered in descending order λ1 ≥ . . . ≥ λn , and let 1 ≤ i, j ≤ n be integers. If i + j ≤ n + 1, then λi+j−1 (A + B) ≤ λi (A) + λj (B). If n < i + j, then λi (A) + λj (B) ≤ λi+j−n (A + B). Corollary 4 (Spectral stability). In the setting of Lemma 23, let 1 ≤ k ≤ n be an integer, then |λk (A + B) − λk (A)| ≤ ∥B∥. We use the following conversion: if ẑ is a unit vector satisfying ∥ẑ − z/∥z∥∥ ≤ ξ, then one may write z+e ẑ = for some e with ∥e∥ ≤ ξ∥z∥. (68) ∥z + e∥ Lemma 24. If z ̸= 0 and ẑ is unit with ∥ẑ − z/∥z∥∥ ≤ ξ < 1, then there exists e with ∥e∥ ≤ 2ξ∥z∥ such that (68) holds. Proof. Let z̄ = z/∥z∥. Define e := ∥z∥(ẑ − z̄). Then ∥e∥ ≤ ξ∥z∥ and z + e = ∥z∥ẑ, hence ẑ = (z + e)/∥z + e∥. Applying Lemma 24 to the two gradient direction queries: ĝ1 =

g + k 1 v1 , ∥g + k1 v1 ∥

ĝ2 =

∇f (x + µv) + k2 v2 , ∥∇f (x + µv) + k2 v2 ∥

(69)

with ∥k1 v1 ∥ ≤ O(ξ1 )∥g∥ and ∥k2 v2 ∥ ≤ O(ξ2 )∥∇f (x + µv)∥. Lemma 25. Given Assumption 3, any unit vector v satisfies ∇f (x + µv) = ∇f (x) + µ∇2 f (x)v + rµ

∥rµ ∥ ≤ 12 L2 µ2 .

with

Proof. Standard integral remainder; included for completeness. Z µ Z µ 2 2 ∇f (x + µv) − ∇f (x) = ∇ f (x + tv) v dt = µ∇ f (x)v + (∇2 f (x + tv) − ∇2 f (x))v dt, 0

so ∥rµ ∥ ≤

0

Rµ

1 2 0 L2 t dt = 2 L2 µ .

Lemma 26. Assume the colinearity condition of Lemma 6. Let d1 := g + k1 v1 ,

d2 := ∇f (x + µv) + k2 v2 ,

b d := η µ∥H∥Hv.

If ∥dd22 ∥ = ∥dd11 +d +d∥ , then q 1 − ⟨ĝ1 , ĝ2 ⟩2 ∥d∥ = q . ∥d1 ∥ 1 − ⟨ĝ2 , u⟩2

(70)

Proof. Apply Lemma 15 to the pair (v, g) = (d1 , d): it expresses the norm ratio ∥d∥/∥d1 ∥ via the sines of the angles between d1 , d1 + d and d. Noting that ĝ1 = d1 /∥d1 ∥ and ĝ2 = d2 /∥d2 ∥ by (69), we obtain (70). 39

D

Basic Lemmas for Finding Stationary Points

We give the proof details of Section 6 in this appendix. First we refer to a lemma for Taylor approximation error analysis. Lemma 27 (Lemma 1 of [33]). Suppose that f : Rd → R has L2 -Lipschitz Hessian. Then for all x, p ∈ Rd , denote 1 Tx (x + p) = f (x) + ⟨∇f (x), p⟩ + p⊤ ∇2 f (x)p, 2 we have L2 |f (x + p) − Tx (x + p)| ≤ ∥p∥3 . (71) 6 We give a lemma for the descent guaranty of quadratic approximation. Lemma 28.√Given a quadratic q(p) = q(0) + b⊤ p + 12 p⊤ Ap with symmetric A, and assume that ∥A∥ ≥ L2 ϵ, and A is not positive or negative definite. Denote the eigenvector of A that corresponding to the largest magnitude eigenvalue as v and assume that ⟨v, b⟩ ≤ 12 ∥b∥. For r > 0 and let p⋆ ∈ arg min∥p∥≤r q(p). Then either: • (Interior) ∥p⋆ ∥ < r and ∇q(p⋆ ) = b + Ap⋆ = 0; • (Boundary) ∥p⋆ ∥ = r and q(p⋆ ) − q(0) ≤ − 12 r∥∇q(p⋆ )∥. Proof. If ∥p⋆ ∥ < r, first-order optimality gives b + Ap⋆ = 0. Otherwise ∥p⋆ ∥ = r and there exists λ ≥ 0 with b + Ap⋆ + λp⋆ = 0 and A + λI ⪰ 0. Then ∇q(p⋆ ) = b + Ap⋆ = −λp⋆ so ∥∇q(p⋆ )∥ = λr. Moreover, q(0) − q(p⋆ ) = −b⊤ p⋆ − 12 (p⋆ )⊤ Ap⋆ . From b + Ap⋆ + λp⋆ = 0 we get −b⊤ p⋆ = λr2 + (p⋆ )⊤ Ap⋆ . Using A + λI ⪰ 0 we have (p⋆ )⊤ Ap⋆ ≥ −λr2 , hence q(0) − q(p⋆ ) ≥ 21 λr2 = 21 r∥∇q(p⋆ )∥.

E

Lazy Trust Region Method with Evaluation Oracle

In this section, we discuss how to combine our approach and lazy-type methods preliminarily. Consider the following problem: Problem 1. Given an objective function f : Rn → R that has L1 -Lipschitz gradient and L2 -Lipschitz Hessian. The query oracle of f is zeroth-order, i.e., we are given the oracle Of : x 7→ f (x). The goal is to output a list containing an FOSP. In lazy-type methods, we update Hessian every m steps. We guarantee that in these m steps, the Hessian in the current point and delayed Hessian are close enough so that using the delayed Hessian introduces bearable error. According to this idea we obtain the following theorem: √ Theorem 9. Take m = n. Then in O( n/ϵ1.5 + n) iteration steps, with high probability, we can visit an ϵ-first order stationary point. p Proof of Theorem 9. Take stepsize r = c ϵ/n where c is an absolute constant that is small enough, and let xkm+t+1 = xkm+t + p. By Lemma 27, we have L2  ϵ  32 |f (xkm+t+1 ) − Txkm+t (p)| ≤ . 6 n 40

(72)

Algorithm 6: LazyTrustRegion(x0 , ϵ) Input: Starting point x0 , precision ϵ Output: an ϵ-first order stationary point xT 1 for k = 0, 1, . . . , T do e km = HessianEstimation(xkm ). 2 H 3 for t = 0, 1. . . . , m − 1 do 4 gkm+t ← GradientEstimation(xkm+t ). 5 xkm+t+1 ← o n e km (x − xkm+t ) min∥x−xkm+t ∥≤r ⟨∇f (x), x − xkm+t ⟩ + 12 (x − xkm+t )⊤ H ≜ min∥x−xkm+t ∥≤r {mkm+t (p)}. 6 return xT +1

Due to the L2 -Lipschitzness of Hessian, after t ≤ m steps, we have e km ∥ ≤ L2 ∥xkm+t − xkm ∥ ≤ L2 ∥Ĥkm+t − H

km+t−1 X

∥xi+1 − xi ∥ ≤ L2 · rm.

(73)

i=km

When ∥xt+1 − xt ∥ ≤ r, the error, denoted ∆err , can be estimated as 1 e km )(x − xkm+t ) ∆err := |Txkm+t (p) − mkm+t (p)| = (x − xkm+t )⊤ (Ĥkm+t − H 2 1 L2 3 ≤ r2 L2 · rm = r m = O(r3 · m). 2 2

(74)

From Lemma 28, we have r |Txkm+t (p) − f (xkm+t )| ≥ ∥∇f (xkm+t ) + Ĥkm+t · p∥ = O(rϵ). (75) 2 p Choosing c small enough for r = c ϵ/n, we obtain |∆km+t+1 − ∆km+t | ≥ 2∆err , thus we can get a √ descent at least Θ(ϵ1.5 / n). When ∥xt+1 − xt ∥ < r, we have e km (xkm+t+1 − xkm+t ) + (∇2 f (xkm+t ) − H e km )(xkm+t+1 − xkm+t ) ∇f (xkm+t+1 ) = ∇f (xkm+t ) + H e km )(xkm+t+1 − xkm+t ). = (∇2 f (xkm+t ) − H (76) As a result, we can obtain e km ∥ · ∥xkm+t+1 − xkm+t ∥ ≤ rd · r ≤ ϵ. ∥∇f (xkm+t+1 )∥ ≤ ∥∇2 f (xkm+t ) − H

(77)

√ In Algorithm 6, the total times that we update the gradient is Titer = O( n/ϵ1.5 ). Therefore, the total query complexity is      1.5  n 1 2 n · Titer + Titer · +1 n =O + n . (78) m ϵ1.5

41