HOUDINI: Escaping from Moderately Constrained Saddles

1 Introduction

Achieving convergence of gradient descent to an (approximate) local minimum is a central question in non-convex optimization for machine learning. In recent years, breakthrough progress starting with the work of [Ge et al., 2015] has led to a flurry of results in this area (see e.g. [Jin et al., 2017; Du et al., 2017; Mokhtari et al., 2018; Jin et al., 2018; Carmon et al., 2017; Carmon and Duchi, 2018; Staib et al., 2019; Carmon and Duchi, 2020; Jin et al., 2021]), culminating in almost optimal bounds [Zhang and Li, 2021]. However, despite this success a key open question of [Ge et al., 2015] still remains unanswered - can gradient methods efficiently escape from saddle points in constrained non-convex optimization? In fact, even basic linear inequality constraints still remain an obstacle: "Dealing with inequality constraints is left as future work" [Ge et al., 2015] 1 . This is due to the NP-hardness of the related copositivity problem [Murty and Kabadi, 1987], which corresponds to the case when the number of constraints is linear in the dimension. In this paper, we make progress on this open question in the case when the number of constraints depends moderately on the dimension.

Consider a feasible set defined by k linear inequality constraints: S = {x ∈ Rd | Ax ≤ b}, where A ∈ Rk×d and b ∈ Rk . Let Bd(x, r) be a d-dimensional closed ball of radius r centred at x. We write B(x, r) when the dimension is clear from the context and drop x when x = 0. Our goal is to find minx∈S f(x) for the objective function f : Rd → R over S. We first introduce the standard smoothness assumption.

1 Using Lagrangian multipliers, equality constraints can be seen as reducing the dimension of the otherwise unconstrained problem.

Assumption 1 (Smoothness). The objective function f satisfies the following properties:

1. (First order) f has an L-Lipschitz gradient (f is Lsmooth): ∥∇f(x) -∇f(y)∥ ≤ L∥x-y∥, ∀x, y ∈ Rd .
2. (Second order) f has a ρ-Lipschitz Hessian: ∥∇2 f(x)- ∇2 f (y)∥ ≤ ρ∥x - y∥, ∀x, y ∈ Rd .

Definition 1 (Local minimum). For S ⊆ Rd , let f : Rd → R. A point x⋆ is a local minimum of f in S if and only if there exists r > 0 such that f(x) ≥ f(x⋆ ) for all x ∈ S∩B(x⋆ , r).

Since finding a local minimum is NP-hard even in the unconstrained case (see e.g. [Anandkumar and Ge, 2016] and the references within) the notion of a local minimum is typically relaxed as follows.

Definition 2 (Approximate local minimum). For S ⊆ Rd and f : Rd → R a point x⋆ is a (δ, r)-approximate local minimum if f(x) ≥ f(x⋆ ) - δ for all x ∈ S ∩ B(x⋆ , r).

For smooth functions, one can define stationary points in terms of the gradient and the eigenvalues instead:

Definition 3 ([Nesterov and Polyak, 2006; Jin et al., 2021]). A point x is an ε-second-order stationary point (ε-SOSP) if ∥∇f(x)∥ < ε and λmin(∇2 f(x)) > -√ρε, where λmin denotes the smallest eigenvalue.

When applying this definition to the constrained case, eigenvectors and eigenvalues are not well-defined since there might be no eigenvectors inside the feasible set, while an escaping direction might exist. Moreover, for f(x) = - 1 2 ∥x∥2 and any compact feasible set, the Hessian is -I at any point with λmin(-I) = -1. Hence an ε-SOSP doesn't exist according to the Definition 3, even though a local minimum exists. In fact, Definition 3 arises from the Taylor expansion, which justifies the choice of √ρϵ as the bound on the smallest eigenvalue. If the function has a ρ-Lipschitz Hessian:

To guarantee that the discrepancy between the function and its quadratic approximation is small relative to δ (from Definition 2), a natural choice of r is 3√δ/ρ, which bounds the discrepancy with Θ(δ). Therefore, based on the quadratic approximation, one can distinguish a (δ, r)-approximate local minimum from not a (cδ, r)-approximate local minimum for c < 1. By setting this r and selecting ε = 3√δ2ρ, we have √ρε = 3√δρ2 and for any h ∈ B(x, r) (see the full version):

Using the ball radius discussed above we arrive at the following version of Definition 2:

Definition 4 (Approximate SOSP). For S ⊆ Rd , let f : Rd → R be a twice-differentiable function with a ρ-Lipschitz Hessian. A point x⋆ is a δ-second-order stationary point (δSOSP) if for r = 3√δ/ρ:

Note that, while Definition 4 doesn't seemingly use second-order information, our choice of the ball radius guarantees that the function is close to its quadratic approximation. In particular, we can determine whether f is a cδ-SOSP for a constant c by only using second-order information.

1.1 Our Contribution

Our results hold for stochastic gradient descent (SGD):

Assumption 2. Access to a stochastic gradient oracle g(x):

1. (Unbiased expectation) E[g(x)] =∇f(x).

Our main result is the following theorem which quantifies the complexity of finding an approximate SOSP under a moderate number of linear inequality constraints, showing that this problem is solvable in polynomial time for k = O(log d). We refer to a function as (L, ρ)-smooth if it satisfies Assumption 1 and simply second-order smooth if both smoothness parameters are constant.

Theorem 5. Let S be a set defined by an intersection of k linear inequality constraints. Let f be a second-order smooth bounded function. Given access to a stochastic gradient oracle satisfying Assumption 2, there exists an algorithm which for any δ > 0 finds a δ-SOSP in ˜O( 1δ (d3 2k + d2 σ2 δ4/3 )) time using ˜O( dδ (1+ dσ2 δ4/3 )) stochastic gradient oracle calls. In the deterministic gradient case (σ = 0), the time complexity is ˜O( d3 2k δ ) and the number of gradient oracle calls is ˜O( dδ ).

The exponential dependence of time complexity on k in our results (not required in the oracle calls) is most likely unavoidable due to the following hardness result, which implies that when k = d then the complexity of this problem can't be polynomial in d under the standard hardness assumptions.

Remark 6 (Matrix copositivity [Murty and Kabadi, 1987]). For a quadratic function f(x) = xT Mx subject to constraints xi ≥ 0 for all i, it is NP-hard to decide whether there exists a solution with f(x) > 0.

Related results in convex optimization are covered in [Boyd and Vandenberghe, 2004; Bubeck and others, 2015]. Among related results in non-convex optimization, here we only focus on the algorithms using only gradient information.

1.2 RelatedWork

Unconstrained optimization. Recall that an ϵ-first-order stationary point (ϵ-FOSP) is defined so that ∥∇f(x)∥ ≤ ϵ. Analyses of convergence to an ε-FOSP are a cornerstone of non-convex optimization (see e.g. classic texts [Bertsekas, 1997; Nocedal and Wright, 1999]). Quantitative analysis of convergence to an ϵ-SOSP (Definition 3) started with the breakthrough work by [Ge et al., 2015] further refined in [Jin et al., 2017; Carmon and Duchi, 2018; Jin et al., 2018; Carmon and Duchi, 2020; Jin et al., 2021] and most recently in [Zhang and Li, 2021], who show an almost optimal bound. Due to the prevalence of SGD in deep learning, stochastic methods have attracted the most attention (see [AllenZhu, 2018; Allen-Zhu and Li, 2018; Fang et al., 2018; Tripuraneni et al., 2018; Xu et al., 2018; Zhou and Gu, 2020; Zhou et al., 2020] for the case of Lipschitz gradients and [Ge et al., 2015; Daneshmand et al., 2018] for non-Lipschitz gradients). For an in-depth summary of the previous work on unconstrained non-convex optimization we refer the reader to [Jin et al., 2021].

Constrained optimization. The case of equality constraints is typically reducible to the unconstrained case by using Lagrangian multipliers (see e.g. [Ge et al., 2015]). However, the general constrained case is substantially more challenging since even the definitions of stationarity require a substantial revision. For first-order convergence a rich literature exists, covering projected gradient, Frank-Wolfe, cubic regularization, etc (see e.g. [Mokhtari et al., 2018] and the references within). For second-order convergence, the landscape of existing work is substantially sparser due to NP-hardness (Remark 6, [Murty and Kabadi, 1987]). A large body of work focuses on achieving convergence using various forms of NP-oracles (see e.g. [Bian et al., 2015a; Cartis et al., 2018; Mokhtari et al., 2018; Haeser et al., 2019; Nouiehed and Razaviyayn, 2020]), while another approach is to define stationarity in terms of tight constraints only [Avdiukhin et al., 2019; Lu et al., 2020].

Relationship with other definitions of SOSP. As discussed in Remark 6, second-order constrained optimization is NP-hard due to the hardness of the matrix copositivity problem. Definitions of constrained SOSP in the previous work fall into two categories: 1) definitions of scaled stationary points, 1) definitions that only consider active constraints ("active constraints only" definitions), 2) definitions that preserve the NP-hardness of the problem and rely on NP-oracles to achieve polynomial-time convergence:

1. (Scaled) For the constraints x ≥ 0, [Bian et al., 2015b; O'Neill and Wright, 2020] consider the definition of scaled SOSP. The idea is to scale i-th coordinate by xi: i.e. instead of bounding ∇f(x) it bounds X∇f(x), where X = diag(x1, . . . , xd), and instead of eigenvalues ∇2 f(x) it considers eigenvalues of X∇2 f(x)X. For this definition, x ≥ 0 restricts possible applications.
2. ("Active constraints only") In [Avdiukhin et al., 2019; Lu et al., 2020] definitions analogous to Definition 3, and the second-order conditions are given with respect to the set of active (i.e. tight for the current iterate) con-

straints. This bypasses the NP-hardness since the point at which the hardness of the copositivity problem applies now becomes a stationary point by definition.

1. (NP-hard) In the results relying on NP-oracles (e.g. [Bian et al., 2015a; Mokhtari et al., 2018; Haeser et al., 2019; Nouiehed and Razaviyayn, 2020]) the complexity is shifted on solving black-box quadratic optimization problems of a certain type. A key advantage of these types of approaches is that they can handle an arbitrary number of constraints and hence promising in certain machine learning applications.

What is currently lacking in the state of the art is a quantitative analysis of the complexity of convergence to a secondorder stationary point, which shows full dependence on both the dimension and accuracy while defining stationarity with respect to the full set of constraints, instead of just active constraints only 2 . Our goal in Theorem 5 is to address this gap and give such an analysis.

1.3 Technical Overview

We address the NP-hardness of the copositivity problem by focusing on the case of a moderate number of constraints and arguing that it can be addressed using gradient-based methods. In order to streamline the presentation we first focus on the key challenge of escaping from a saddle point in a corner defined by the constraints when the underlying function is simply quadratic (Section 2). This is already enough to capture some of the key contributions, while more technical details and the full algorithm are given in Section 3.

Quadratic corner saddle point (Section 2). In this simplified scenario, the NP-hardness comes from the fact that the point we aim to find can lie in the intersection of an arbitrary subset of constraints. By doing an exhaustive search over this set of constraints (Algorithm 1) and enforcing them throughout the search process we are able to reduce to a setting similar to the equality-constrained case (Algorithm 2). We show different subsets of constraints enforced by the algorithm in an example in Figure 1. The key challenge is making this argument formal and arguing that this process converges to a constrained approximate SOSP as in Definition 4. This relies on performing a robust analysis of the properties of the

2 While unrefereed manuscripts ([Hsia and Sheu, 2018] and [Nouiehed et al., 2020], relying on [Hsia and Sheu, 2018] as a subroutine) do contain related results, our work differs in a number of important aspects. [Hsia and Sheu, 2018] and [Nouiehed et al., 2020] require access to the exact gradient and Hessian. We only assume access to the stochastic gradient oracle. Furthermore, [Hsia and Sheu, 2018] assumes access to matrix diagonalization. However, the diagonalization can only be found approximately, which compromises the stability of this approach, especially with respect to the linear term transformations. Our approach handles this issue by appropriate perturbation and using known results for matrix diagonalization. Finally, compared with [Hsia and Sheu, 2018], our main contribution is focused on solving a substantially different problem. While [Hsia and Sheu, 2018] find a global minimum of a quadratic problem, we find an approximate local minimum of an arbitrary smooth non-convex function.

-1

-2

-1

-2

-1

-2

-1

-2

points which lead to the hardness of copositivity. In particular, we show that after we guess the set of constraints correctly, the problem reduces to finding the smallest eigenvector (Lemma 9). Exact error analysis of the eigenvector process (Lemma 10, Lemma 11) is then required to complete the proof of the main theorem (Theorem 8).

General case (Section 3). In the algorithm for the general case, we address the three assumptions made in the quadratic corner saddle point case, while also handling the stochasticity in the gradient. The latter part is standard and is handled via variance reduction (see the full version). The full algorithm iterates the escape subroutine (Algorithm 3) until an escaping point is found. The escape subroutine first approximates the Hessian matrix using the gradient oracle and then performs an exhaustive search over the set of active constraints at the escaping point in a way similar to the quadratic corner case. After the correct subset of constraints is fixed the current iterate needs to be projected on this set of constraints, which also necessitates a recomputation of various related parameters. When this is done the problem is solved by a subroutine Algorithm 4, analogous to Algorithm 2 from the quadratic corner case.

However, since the function is no longer quadratic and the gradient can be large, several modifications are required. The algorithm tries to find an escaping direction within a ball of radius r, within which the function is well approximated by a quadratic by the smoothness assumption (as discussed above). First, the algorithm tries to escape using the gradient term. If that doesn't work then we consider two cases: 1) the solution is inside the ball of radius r = 3√δ/ρ (from Definition 4), 2) the solution is on the boundary of this ball. In the first case, the solution is a critical point and hence can be found using the Newton step. In the second case, we diagonalize the Hessian using an orthogonal transformation. This gives a quadratic function with a linear term whose critical points on the boundary of the ball can be found explicitly (up to the required precision) as roots of the corresponding polynomial. We note that the diagonalization performed in this case is the most computationally expensive step in the algorithm, resulting in polynomial dependence on the dimension.

Notation. For a set S let IntS be its interior and ∂S be its boundary. For x ∈ Rd and S ⊆ Rd , ProjS(x) = argminy∈S ∥x - y∥ is the projection of x on S. For a square matrix M with eigenvalues λ1 ≤ . . . ≤ λd, we denote λmin(M) = λ1 and λ|max|(M) = max(|λ1|, |λd|). For S = {x ∈ Rd | Ax ≤ b} and x ∈ S, we say that i-th constraint is active at x if A⊤i x = bi, where Ai is the i-th row of A. ˜O notation hides polylogarithmic dependence on all parameters, including error probability.

2 Quadratic Corner Saddle Point Case

We introduce the key ideas of the analysis in a simplified setting when: 1) the function f is quadratic, 2) the gradient is small, 3) the current iterate is located in a corner of the constraint space. Intuitively, this represents the key challenge of the constrained saddle escape problem since its NP-hardness comes from the hardness of the matrix copositivity problem in Remark 6 (i.e. the function is exactly quadratic and has no gradient at the current iterate which lies in the intersection of all constraints). We refer to this setting as the Quadratic Corner Saddle Point problem defined formally below. By shifting the coordinate system, w.l.o.g. we can assume that the saddle point is 0 and f(0) = 03 . If 0 is a (δ, r)-QCSP (as defined below), the objective can be decreased by δ within a ball of radius r.

Definition 7 (Quadratic Corner Saddle Point). Let S = {x | Ax ≤ 0}. For δ, r > 0 and function f(x) = 12 x⊤Mx, we say that a point 0 is a:

• (δ, r)-Quadratic Corner Saddle Point ((δ, r)-QCSP) if minx∈B(r)∩S f (x) < -δ.
• (δ, r)-boundary QCSP if minx∈B(r)∩∂S f(x) < -δ.

In this section, we show how to escape from a (δ, r)-QCSP (see the full version for the full proof):

Theorem 8 (Quadratic Corner Saddle Point Escape). Let δ, r > 0. Let f(x) = 1 2 x⊤Mx with λ|max|(M) ≤ L and

3 Algorithms 1 and 2 don't require saddle point x to be 0. All the statements are trivially adapted for the case when x is not 0

input: Starting point x, feasible set S = {y | A(y - x) ≤ 0 ∈ Rk } parameters: δ and r from definition of (δ, r)-QCSP for I ∈ 2[k] - every subset of constraints do 1 A← {y | A⊤i (y-x) = 0 for i ∈ I}, where Ai is the i-th row of A // Optimize in A 2 y ← FINDINSIDECORNER(x,A) 3 if y ∈ S and f(y) < f(x) - δ2 then 4 return y

5 return ⊥

Algorithm 2: FINDINSIDECORNER(x, A)

input: Corner x, affine subspace A with x ∈ A

parameters: δ and r from Definition 7, step size

η = 1L , number of iterations T = ˜O( Lr2 δ ) 1 Sample ξ ∼N(0, I), x0 ←ProjA(x+ ξ) 2 for t = 0, . . . , T -1 do 3 // Power method step 4 xt+1 ←ProjA(xt - η(∇f(xt) -∇f(x))) 5 e ← r xT-x ∥xT -x∥ 6 return x+e

let S = {x | Ax ≤ 0} be defined by k linear inequalities. If 0 is a (δ, r)-QCSP, then Algorithm 1 with probability at least 1 - ξ finds a point x ∈ S ∩ B(r) with f(x) < -Ω(δ) using O ( Lr2 k2k δ log 1ξ ) deterministic gradient oracle calls.

For the rest of the section, we assume that 0 is a (δ, r)- QCSP, i.e. minx∈S∩B(r) f(x) < -δ. We consider two cases depending on whether 0 is a (δ, r)-boundary QCSP.

Case 1: 0 is a (δ, r)-boundary QCSP. For a subset of inequality constraints I ⊆ [k] we define the subspace where these constraints are active: AI = {x | A⊤i x = 0 for all i ∈ I}. Let I be a maximal 4 subset of constraints such that minx∈AI∩B(r) f(x) < -δ. If P is a projection operator on AI, it suffices to optimize g(x) := f(Px) = 1 2 x⊤(PMP)x. Therefore, we reduced the original problem to minimizing a different quadratic form in the same feasible set. For any i ∈ I, A⊤i Px ≤ 0 holds trivially, since A⊤i y = 0 for any y ∈ A, and hence constraints from I can be ignored. If a constraint not from I is active in Px, then f(Px) ≥ -δ, since I is a maximal subset of constraints with minx∈AI∩B(r) < -δ. Therefore, this reduces Case 1 to the next case.

Case 2: 0 is not a (δ, r)-boundary QCSP. In this case, we show that any x ∈ B(r) with f(x) < -δ must lie in S, and for f(x) = 1 2 x⊤Mx it suffices to find the eigenvector corresponding to the smallest eigenvalue of M. We first show that there exists an eigenvector improving the objective.

4 As we don't know I, Algorithm 1 tries all subsets of constraints.

Lemma 9. Let f(x) = 12 x⊤Mx and S = {x | Ax ≤ 0}. If 0 is not a (δ, r)-boundary QCSP for δ, r > 0, then the following statements are equivalent:

1. 0 is (δ, r)-QCSP, i.e. minx∈S∩B(r) f(x) < -δ.
2. There exists an eigenvector e ofMsuch that e ∈ IntS∩ ∂B(r) and f(e) < -δ.

Finding an exact eigenvector might be impossible. Hence, we show that finding x ∈ B(r) with f(x) < -δ suffices, since either x or -x are in S.

Lemma 10. Let f(x) = 12 x⊤Mx and S = {x | Ax ≤ 0}. For δ, r > 0 and ˆx ∈ ∂B(r), if f(ˆx) < -δ and the following conditions hold, then either ˆx ∈ S or -ˆx ∈ S:

1. minx∈S∩B(r) f(x) < -δ, i.e. 0 is a (δ, r)-QCSP,
2. minx∈∂S∩B(r) f(x) ≥ -δ, i.e 0 is not a (δ, r)-boundary QCSP,

To show Lemma 10, we use the fact that, by Lemma 9, there exists an eigenvector e with e ∈ IntS ∩ ∂B(r) and f(e) < -δ. For the sake of contradiction, if both -ˆx and ˆx are not in S, at least one of them has a non-negative inner product with e. W.l.o.g. we assume ˆx⊤e ≥ 0, and we consider an arc on ∂B(r) connecting ˆx and e. Since e ∈ S and ˆx /∈ e, the arc intersects ∂S at some point. We show that any point x on the arc has f(x) < -δ, and hence this also holds for the point on the boundary, contradicting the assumption that 0 is not a (δ, r)-boundary QCSP, finishing the proof.

The main idea behind the algorithm is that Algorithm 2 emulates the power method on matrix I - ML . Hence, for any ε it allows us to find a vector x such that

Since λ|max|(M) ≤ L, all eigenvalues of I - ML are positive, and hence the power method approximates the eigenvector corresponding to the largest eigenvalue of I - ML , and hence to the smallest eigenvalue of M.

We know that f(x) < -δ, and we aim to find x ∈ B(r) with f(x) < -(1-ε)δ for a constant ε. Since there exists an eigenvector e ∈ ∂B(r) of M with 12 e⊤Me < -δ, we have λmin(M) < -2δ ∥e∥2 = -2δ r2 , and hence the largest eigenvalue of I - ML is at least 1 - λmin(M) L ≥ 1 + 2δ Lr2 . Finding x with 12 x⊤Mx < -(1-ε)δ is equivalent to finding x with x⊤(I - ML )x ≥ (1 + 2(1-ε)δ Lr2 )∥x∥2 , and the power method achieves this in O(log d + Lr2 εδ ) iterations (see proof of Lemma 11 in the full version).

Lemma 11. Let δ, r > 0, x ∈ Rd and ε ∈ (0, 1). Let f(x) = 1 2 x⊤Mx with λ|max|(M) ≤ L. Let A be a linear subspace of Rd . If minx∈A∩B(r) f(x) < -δ, then Algorithm 2 finds x ∈ A ∩ ∂B(r) with f(x) ≤ -(1 - ε)δ after T = ˜O ( Lr2 εδ ) iterations w.h.p.

Finally, we now prove Theorem 8. First, by exhaustive search, we guess a maximal subset of active constraints I such that subspace AI formed by these linear constraints has x ∈ B(r) ∩ S with f(x) < δ. Using Algorithm 2, we find y ∈ B(r) ∩ AI with f(y) < -(1 - ε)δ. Then y ∈ S by Lemma 11, since I is a maximal subset of constraints with an escape direction.

Algorithm 3: HOUDINIESCAPE(x, S, δ): Escaping

from a saddle point

4

5

8

9

input : Saddle point x, δ from definition of δ-SOSP, feasible set S = {y | Ay ≤ b ∈ Rk } output: either reports that x is a δ-SOSP or finds u ∈ S ∩ B(x, r) with f (u) < f (x) - Ω(δ) 1 Construct f ′ (x+h) - quadratic approximation of f(x + h) - using stochastic gradient oracle calls 2 for I - every subset of constraints do 3 A← {y | A⊤i x = bi, i ∈ I}, where Ai is the i-th row of A // Optimize in A Let p be the projection of x on A Let O ∈ Rd×dimA be an orthonormal basis of A 6 Define g(y) := f ′ (p+Oy) 7 Algorithm 4 tries to find an escape direction for g if Algorithm 4 finds direction y then return p + Oy

• 10 Didn't find escape direction: report that x is a δ-SOSP

Algorithm 4:

FINDINSIDE(x, δ, (M⊥, v⊥), (r⊥, S⊥))

input : δ from definition of δ-SOSP,

g(y) = 12 y⊤M⊥y + y⊤v⊥ - objective in

RdimA,

S⊥ ∩ B(r⊥) - feasible set in RdimA

output: Escaping direction, if exists

• 1 If any of the following candidates lies in the feasible set and decreases the objective by Ω(δ), return it:
• 2 Case 1. Large gradient: argminy∈S⊥∩B(r⊥) y⊤v⊥
• 3 Case 2. Solution in the interior: y ←-M-1 ⊥ v⊥
• 4 Case 3. Solution on the boundary:

5 Find orthogonal Q and diagonal

Λ = diag(λ1, . . . , λdimA) such that

• 6 Let ˜v ←Qv
• 7 We consider points with coordinates yi ← ˜ vi µj-λi for some µ
• 8 Find the values of µ for which the points have norm r⊥
• 9 The candidates are the points corresponding to these values of µ
• 10 If no candidate satisfies the condition, return ⊥

3 Main Result

Our main result is the following theorem that shows how to find a δ-SOSP.

Theorem 12. Let S = {x | Ax ≤ b} be a set defined by an intersection of k linear inequality constraints. Let f satisfy Assumptions 1 and 2 and let minx∈S f (x) = f ⋆ . Then there exists an algorithm which for δ > 0 finds a δ-SOSP in

˜O( f (x0)-f ⋆ δ d3 (2k + σ2 δ4/3 )) time using ˜O( f(x0)-f ⋆ δ (d+ d3 σ2 δ4/3 )) stochastic gradient oracle calls.

Our approach is outlined in Algorithms 3 and 4 5 :

• If x is not a δ-SOSP, Algorithm 3 finds an escape direction, i.e. a point y ∈ S ∩ B(x, r) which significantly decreases the objective value: f(y) < f(x) - Ω(δ). Therefore, if x0 is the initial point, our algorithm requires O( f (x0)-f ⋆ δ ) calls of Algorithm 3.
• Algorithm 3 enumerates all possible sets of active constraints

Recall that we aim to find a δ-SOSP, i.e. x ∈ S such that f(y) ≥ f(x) - δ, ∀y ∈ S ∩ B(x, r), where r = 3√δ/ρ.

Similarly to Section 2, we use a quadratic approximation and consider two cases depending on whether the minimizer y lies in the interior or on the boundary of S. The second case can be reduced to the first in a way similar to Section 2. To find the minimizer in the interior, we consider the following cases in Algorithm 3:

• Case 1. When the gradient is large, we ignore the quadratic term, and optimize the function based on the gradient. This covers the situation when the objective is sensitive to the change in the argument.
• Case 2. When the optimum lies in IntB(x, r), none of the constraints are active, and hence we can find the unique unconstrained critical point directly.
• Case 3. When the optimum lies in ∂B(x, r), the only active constraint is c(y) = ∥y - x∥2 - r2 = 0. By KKT conditions, for y to be the minimizer, there must exist µ such that ∇f(y) = µ∇c(y). We show that for each µ, there exists a unique y(µ) satisfying the condition above, and only for O(d) values of µ we have y(µ) ∈ ∂B(x, r), resulting in O(d) candidate solutions.

We now outline the proof of Theorem 12 (full proof is in the full version). As shown in Section 3.1, we can consider a quadratic approximation of the objective. By guessing which constraints are active at the minimizer x⋆ and enforcing these constraints, we restrict the function to some affine subspace A. By parameterizing A, we eliminate enforced constraints, and, since the rest of the constraints are not active at x⋆ , we need to optimize a quadratic function in the intersection of a ball and linear inequality constraints (see the full version).

3.1 Algorithm 3: Quadratic Approximation

The goal is to build a quadratic approximation of the objective with some additional properties (see below). To simplify the presentation, w.l.o.g. (by shifting the coordinate system) we assume that the current saddle point is 0 and consider the quadratic approximation of the function:

Since f is ρ-Lipschitz and r = 3√δ/ρ, in B(r) the quadratic approximation deviates from f by at most δ6 (see derivation before Definition 4). For a small value ξ (to be specified later), we instead analyze a noisy function

5 Algorithm 3 and 4 are simplified versions of rigorous algorithms in the full version

where:

1. v is a perturbed approximation of ∇f(0). The perturbation guarantees that w.h.p. all coordinates of v are sufficiently separated from 0 and linear systems of the form (M- µI)x = v with rank(M- µI) < d don't have any solutions, simplifying the analysis. To approximate the gradient, we use algorithm VRSG (see the full version), which w.h.p. estimates the gradient with precision ˜σ using ˜O( σ2 ˜σ2 ) stochastic gradient oracle calls.
2. M is a perturbed approximation of ∇2 f(0). The perturbation guarantees that M is non-degenerate with probability 1. Since the i-th column of ∇2 f (0) is by definition lim τ→0 ∇f (τ ei)-∇f (0) τ , using a sufficiently small τ and approximating ∇f (τ ei) and ∇f (0) using VRSG, we find good approximation of the Hessian.

Combining this with the derivation before Definition 4, we show the following:

Lemma 13. Let f satisfy Assumptions 1 and 2. Let f ′ (x) = f(0) + x⊤v + 1 2 x⊤Mx, where v and M are as in Algorithm 3. For δ > 0, r = 3√δ/ρ we have ∥f ′ (x) - f(x)∥ < δ2 for all x ∈ B(r) w.h.p.

Reducing Case x⋆ ∈ ∂S to Case x⋆ ∈ IntS. Similarly to Section 2, we reduce the case x⋆ ∈ ∂S to the case x⋆ ∈ IntS. If x⋆ ∈ ∂S, then there exist a non-empty set I of constraints active at x⋆ . Consider the iteration of Algorithm 3 where constraints from I are active. These active constraints define an affine subspace A, which e parameterize: if p = ProjA(x) and O ∈ Rd×dimA is an orthonormal basis of A, then any point in A can be represented as p +Oy for y ∈ RdimA. Defining g(y) = f ′ (p + Oy), minimizing f ′ in A ∩ S ∩ B(r) is equivalent to minimizing g in S⊥ ∩ B(r⊥), where:

1. S⊥ is a set of points y ∈ RdimA such that p+Oy ∈ S, namely S⊥ = {y | Ai(p+Oy) ≤ bi, i /∈ I}. Hence, S⊥ is defined by linear inequalities, similarly to S.
2. r⊥ is a radius such that condition y ∈ B(r⊥) is equivalent to p +Oy ∈ B(r). Since p is the projection of 0 on A, we have O⊤p = 0, and hence ∥p + Oy∥2 = ∥p∥2 + ∥Oy∥2 . Since O is an orthonormal basis of A, ∥Oy∥ = ∥y∥, and hence r⊥ = √r2 - ∥p∥2 .

For y⋆ such that x⋆ = p +Oy⋆ , no constraints from S⊥ are active, and hence x⋆ ∈ Int S⊥.

3.2 Algorithm 4: EscapingWhen y⋆ ∈ IntS⊥

In this section, we the find minimizer y⋆ of function g(y) = 1 2 y⊤M⊥y+y⊤v⊥ +C in S⊥ ∩B(r⊥), while assuming that y⋆ ∈ Int S⊥. Since the solutions we find can be approximate, we have to guarantee that the objective is not too sensitive to the change of its argument. It suffices to consider the case when ∥v⊥∥ is bounded, since for any y ∈ B(r⊥) and perturbation h there exists τ ∈ [0, 1] such that:

Figure 2: Cases of Algorithm 4

where we used that the objective is L-smooth and hence ∥∇g(y + τh)∥ ≤ ∥∇g(0)∥ + L∥y + τh∥. We consider the situation when ∥v⊥∥ is large as a separate case. Otherwise, for y⋆ , there are only two options: either y⋆ ∈ IntB(r⊥) or y⋆ ∈ ∂B(r⊥). Algorithm 4 handles these cases, as well as the case when ∥v∥ is large, separately.

Case 1: ∥v⊥∥ is large. If ∥v⊥∥ is large and we can find y with small y⊤v⊥, the linear term alone suffices to improve the objective. We show that, if such y doesn't exist, then g(y⋆ ) requires y⋆ ∈ ∂S⊥, which contradicts that y⋆ ∈ Int S⊥. Below we assume that ∥v⊥∥ is bounded.

Case 2: y ∈ IntB(r⊥). In this case, y⋆ is an unconstrained critical point of g, and hence it must satisfy ∇g(y) = 0, implying M⊥y + v⊥ = 0 which gives the unique solution y = -M-1 ⊥ v⊥. since M is a perturbed matrix, so is M⊥, and hence M⊥ is non-degenerate with probability 1. It remains to verify that y ∈ B(r⊥) ∩ S⊥ and y decreases the objective by Ω(δ).

Case 3: y ∈ ∂B(r⊥). Since the only active constraint at y⋆ is c(y) = 1 2 (∥y∥2 - r2⊥) = 0, by the KKT conditions, any critical point must satisfy ∇g(y) = µ∇c(y) for some µ ∈ R, which is equivalent to M⊥y + v⊥ = µy. Hence, for any fixed µ, y must be a solution of linear system (M⊥-µI)y = -v⊥. When M⊥ - µI is degenerate (i.e. µ is an eigenvalue of M⊥), due to perturbation of v⊥, the system doesn't have any solution with probability 1. It leaves us with the case whenM⊥-µI is non-degenerate, when there exists a unique solution y(µ) := -(M⊥ - µI)-1v⊥.

Diagonalization. Dependence of (M⊥ - µI)-1 on µ is non-trivial, but for a diagonal M⊥, the inverse can be found explicitly. Hence, we perform diagonalization of M⊥: we find orthogonal Q and diagonal Λ such that ∥M⊥ - Q⊤ΛQ∥ < ε in time O(d3 log 1/ε) [Pan and Chen, 1999]. Setting ε = O(δ/r2⊥), we guarantee that the function changes by at most O(δ) in B(r⊥). The function y(µ) = -(Q⊤ΛQµI)-1v⊥ can be written as Qy(µ) = -(Λ - µI)Qv⊥, and hence we work with rotated vectors ˜y(µ) := Qy(µ) and ˜v := Qv⊥.

Finding candidate µ. Since ˜y(µ) = -(M⊥-µI)-1 ˜v, for the i-th coordinate of ˜y(µ) we have ˜ yi(µ) = ˜ vi µ-λi . Since we are only interested in y(µ) ∈ ∂B(r⊥) and Q is an orthogonal matrix, we must have:

After multiplying the equation by ∏i (µ - λi)2 , we get an equation of the form p(µ) = 0, where p is a polynomial of degree 2d. We find roots µ1, . . . , µ2d of the polynomial in time O(d2 log d · log log 1/ε) [Pan, 1987], where ε is the required root precision. For each i, we compute y(µi) and verify whether it lies in B(r⊥) ∩ S⊥ and improves the objective by -Ω(δ).

Precision. We find the roots of the polynomial approximately, and when µ is close to λi for some i, even a small perturbation of µ can strongly affect yi(µ) = ˜ vi µ-λi . We solve this as follows: since ∥˜y(µ)∥ = r⊥, for all i we have |˜ yi(µ)| ≤ r⊥, implying |µ - λi| ≥ |˜ vi | r⊥ . Therefore, µ must be sufficiently far from any λi, where the lower bound on the distance depends on r⊥ ≤ 3√δ/ρ and on |˜ vi|. Noise added to v is preserved in ˜v, and each coordinate is sufficiently separated from 0 w.h.p. This reasoning is formalized in the full version.

4 Conclusion

In this paper, we have shown that it's possible to escape from a constrained second-order stationary point with the logarithmic number of constraints within polynomial time and using only a polynomial number of stochastic gradient oracle calls. We provide experimental results in the full version.

An open question is to determine the conditions that on one hand guarantee escaping from a saddle point in polynomial time even for the linear number of constraints, and on the other hand hold in practice. One such condition can be strict complementarity.

Another open question is handling non-linear constraints. We believe that it can be straightforwardly achieved using techniques from [Ge et al., 2015] by using assumptions on curvature and linear independence of the constraints. Finally, an interesting question would be a simpler algorithm for the general case, e.g. an algorithm resembling the approach from Section 2.