# A randomized algorithm for nonconvex minimization with inexact evaluations
and complexity guarantees

Shuyao Li  Stephen Wright

###### Abstract

We consider minimization of a smooth nonconvex function with inexact oracle
access to gradient and Hessian (but not the function value) to achieve
$(\epsilon_{g},\epsilon_{H})$-approximate second-order optimality. A novel
feature of our method is that if an approximate direction of negative
curvature is chosen as the step, we choose its sense to be positive or
negative with equal probability. We also use relative inexactness measures on
gradient and Hessian and relax the coupling between the first- and second-
order tolerances $\epsilon_{g}$ and $\epsilon_{H}$. Our convergence analysis
includes both an expectation bound based on martingale analysis and a high-
probability bound based on concentration inequalities. We apply our algorithm
to empirical risk minimization problems and obtain gradient sample complexity.

##  1 Introduction

Given a smooth function $f:\mathbb{R}^{d}\to\mathbb{R}$, we seek a point $x$
that approximately satisfies optimality conditions to be a local minimizer of
$f$. Specifically, for small positive tolerances $\epsilon_{g}$ and
$\epsilon_{H}$, we say that $x$ is an
$(\epsilon_{g},\epsilon_{H})$-approximate second-order point when

$\|\nabla
f(x)\|\leq\epsilon_{g},\quad\lambda_{\min}\left(\nabla^{2}f(x)\right)\geq-\epsilon_{H}.$
(1)

We define an algorithm for finding such a point and analyze its complexity,
expressed in terms of $\epsilon_{g}$ and $\epsilon_{H}$, Lipschitz constants,
failure probablities, and algorithmic constants. Our approach is based on an
elementary method in which steps of only two types are taken: gradient descent
steps and negative curvature steps for the Hessian. It allows gradient and
Hessian evaluations to be inexact, where the amount of gradient inexactness
can be measured in a relative sense (that is, as a fraction of the size of the
quantity being approximated). Moreover, it does not require evaluation of the
value of the objective function $f$ 111Coined “objective-function-free
optimization (OFFO)” in, e.g., [GJT23], which is an essential element of line-
search and trust-region methods. In such methods, the function value
evaluation may need to be quite precise (and thus expensive and intolerant of
noise [CS18, CMS18, Bla+19, CSS19, BCS21]) when the amount of decrease in $f$
over a step is small. On the other hand, the absence of function information
makes adaptivity more difficult; most algorithms (but not all, an exception is
[GJT23]) tend to fall back on conservative choices of steplengths.

Apart from not requiring function values, a distinctive feature of our method
is that it “flips a coin” to decide on whether to move in a positive or
negative sense along a direction of negative curvature for an approximate
Hessian. (A more technical term for this strategy would be “Rademacher
randomness.”) Because of this feature, our step may not be a descent step for
$f$; this fact is accounted for in our analysis. (We can obtain reasonable
complexity results despite these non-descent steps because such steps are
taken only when the gradient is small, so the amount of any increase in $f$ is
limited.)

In later sections, we pay particular attention to the case in which $f$ has
“finite-sum” form

$f(x)=\frac{1}{N}\sum_{i=1}^{N}f_{i}(x),$ (2)

for a (typically large) value of $N$. We derive bounds on the number of
samples (evaluations of $\nabla f_{i}$ or $\nabla^{2}f_{i}$ for some $i$)
required to meet the required level of exactness for $\nabla f$ and possibly
$\nabla^{2}f$ at each iterate, and the total number of samples needed to find
the point satisfying (1), to high probability.

Finally, our approach to the analysis does not require a strong coupling
between the tolerances $\epsilon_{g}$ and $\epsilon_{H}$. We state several
corollaries to our main result to show how various choices of these
quantities, where both are expressed in terms of some underlying tolerance
$\epsilon$, give rise to various overall complexity expressions. In prior
works, a coupling of $\epsilon_{H}=O(\sqrt{\epsilon_{g}})$ is often assumed;
we find situations where choices for which
$\epsilon_{H}=o(\sqrt{\epsilon_{g}})$ produce better overall bounds.

###  1.1 Prior Works

#### Approximate Second-Order Points.

Smooth nonconvex optimization problems are generally intractable, and one
common approach is to seek approximate second-order stationary points as a
surrogate problem. This approach is advantageous in many practical settings
where local minima are in fact global minima, as in the case (under certain
assumptions) in tensor decomposition [Ge+15], dictionary learning [SQW17],
phase retrieval [SQW16], and low-rank matrix factorization problems [GJZ17]
such as sensing, completion, and robust PCA. See a discussion of this issue in
[Jin+17, Section 1]. A comprehensive technical review [ZQW20] discusses the
loss landscape of problems in which there are no spurious local minima.

There are first-order methods that make perturbations to gradient descent
[Jin+17, Jin+21] or stochastic gradient descent [Ge+15, Jin+21]) directions to
guarantee convergence to second-order stationary points. In some cases these
require careful hyperparameter tuning or high iteration complexity.

Second-order methods can make use of negative curvature directions of the
Hessian (instead of random noise) to escape saddle points more efficiently. We
focus on two types of computational complexity for this class of algorithms.
Iteration complexity describes the number of iterations needed to achieve the
target accuracy, where each iteration may involve solving a difficult trust
region problem or minimizing a cubic function. Operation complexity describes
the total number of gradient evaluations and Hessian-vector products needed to
attain target accuracy.

For trust-region methods, [CGT12] obtained iteration complexity
$O(\max\\{\epsilon_{H}^{-1}\epsilon_{g}^{-2},\epsilon_{H}^{-3}\\})$ and
[GYY16] obtained $O(\max\\{\epsilon_{g}^{-3},\epsilon_{H}^{-3}\\})$ iteration
complexity. For operation complexity, [Cur+21] designed a Newton-CG trust
region method that achieves $\tilde{O}\left({\epsilon_{g}^{7/4}}\right)$,
assuming $\epsilon_{H}=\epsilon_{g}^{1/2}$.

Cubic regularization [NP06] and its adaptive variants [CGT11] find approximate
second-order points with iteration complexity
$O(\max\\{\epsilon_{g}^{-3/2},\epsilon_{H}^{-3}\\})$ and
$O(\max\\{\epsilon_{g}^{-2},\epsilon_{H}^{-3}\\})$, respectively. With more
assumptions on the cubic subproblems, [CGT12] also obtained
$O(\max\\{\epsilon_{g}^{-3/2},\epsilon_{H}^{-3}\\})$ iteration complexity. A
recent algorithm described in [GJT23] matches this complexity while also
eliminating the need for function value evaluations, a favorable property that
our algorithm also enjoys. However, it lacks second-order guarantees when
derivatives are inexact. Unique to this algorithm is its adaptability (i.e.,
it does not require knowledge of Lipschitz parameters), traditionally a
feature necessitating function value evaluations. It achieves adaptivity by
choosing steps based on a steadily more conservative model function that
eventually may become more conservative than the usual Taylor-series-based
upper bound. For the best operation complexity to our knowledge, [Aga+17]
achieved $\tilde{O}\left({\epsilon_{g}^{-7/4}}\right)$, assuming
$\epsilon_{H}=\epsilon_{g}^{1/2}$.

Fixed-step and line-search methods are most relevant to our current paper.
[RW18, ROW20] use a line search method that takes a Newton-CG step when the
gradients are large and compute a negative curvature step when the gradients
are small. They obtain iteration complexity of
$O\left({\max\\{\epsilon_{g}^{-3/2},\epsilon_{H}^{-3}\\}}\right)$ and
$\tilde{O}\left({\max\\{\epsilon_{g}^{-3}\epsilon_{H}^{3},\epsilon_{H}^{-3}\\}}\right)$,
respectively, and their operation complexity is a factor of
$\epsilon_{H}^{-1/2}$ greater than their iteration complexity. [Car+18]
similarly combined accelerated gradient descent and minimum eigenvector to
achieve operation complexity $O\left({\epsilon_{g}^{-7/4}}\right)$, assuming
$\epsilon_{H}=\epsilon_{g}^{1/2}$.

#### Inexact derivatives.

There is increasing interest in finding approximate second-order stationary
points when exact derivatives are not available. Some papers directly model
inexactness in the gradient and Hessian oracles. Others pay special attention
to the sampling setting, where the objective function has a finite-sum
structure $f=\frac{1}{N}\sum_{i}^{N}f_{i}$ and the inexactness comes from
(batch) subsampling of $f_{i}$’s. Since the algorithms and their guarantees
are very different in the sampling and non-sampling settings, we review them
separately.

##### Sampling / stochastic settings.

Specific to stochastic settings, we use the term “total operation complexity”
to describe the total number of gradient or Hessian-vector product evaluations
of a single sample $f_{i}$. Perturbed stochastic gradient descent [Ge+15,
Jin+21] has total operation complexity
$\tilde{O}\left({1/\epsilon_{g}^{4}}\right)$ when
$\epsilon_{H}=\epsilon_{g}^{1/2}$.

[Arj+20] obtained approximate SOSPs with total operation complexity
$O\left({\epsilon_{g}^{-3}+\epsilon_{g}^{-2}\epsilon_{H}^{-2}+\epsilon_{H}^{-5}}\right)$.
With access to stochastic gradients and stochastic Hessian-vector products,
they improve on Natasha 2 [All18] and use stochastic gradient descent with a
variance reduction scheme when gradients are large, and Oja’s method to
compute a negative curvature direction when gradients are small. They also
utilize a step of coin-flipping similar to ours to decide the negative
curvature direction in their complex methods, but their inexactness is assumed
to come only from sampling and they only obtain a constant-probability bound.
If access to full stochastic Hessians is assumed, they improve on [Tri+17] and
consider cubic regularization methods to remove the dependence on the gradient
Lipschitz constant in their complexity bounds. (We do not focus on reducing
the dependence on gradient or Hessian Lipschitz constants.)

[Ber+21] considered batch subsampling settings, combining a Newton-type line
search method with a routine that uses minimum eigenvectors to escape saddle
points. They obtained $1/\epsilon_{g}^{3/2}$ iteration complexity (assuming
$\epsilon_{H}=\epsilon_{g}^{1/2}$), while each iteration might need
$\epsilon_{H}^{-1/2}$ Hessian-vector products. Only an expectation bound on
the iteration complexity was proved.

##### Non-sampling, deterministic inexact settings.

[XRM20] analyzed trust region methods and cubic regularization algorithms;
they considered the settings where Hessians are inexact but gradients and
function values are exact.

Although the results of [Tri+17] discussed above are stated only for sampling
settings, their analysis can be extended to the case where inexactness is for
reasons other than batch subsampling.

[Yao+22] considers a fixed-step approach that has guarantees for finding
approximate second-order stationary points, where the gradients and Hessians
are inexact. The settings are quite similar to ours, and we show that we can
tolerate more inexactness. We discuss the latter issue further in Section 3.

[PS20] considers the iteration complexity of finding an approximate first-
order stationary point for convex and nonconvex objectives. This line search
algorithm that requires access to the approximate function value. [JSX23]
analyzed the oracle complexity and also considered trust region methods, and
similarly obtained only first-order guarantees. We mention those papers
because they allow for gradient relative error as does our algorithm.

###  1.2 Notation

We take order notation $O\left({\cdot}\right)$ and $o(\cdot)$ to have their
usual meanings, and say that $f=\Theta(g)$ for two functions $f$ and $g$ if
both $f=O\left({g}\right)$ and $g=O\left({f}\right)$. We say that
$f=\tilde{O}\left({g}\right)$ if $f=O\left({g\log^{k}g}\right)$ for some
constant $k$.

For an event $A$, we use $\mathbf{1}_{A}$ to denote the indicator random
variable that equals 1 if the sample is in $A$ and equals 0 otherwise. We use
$\left\|{x}\right\|$ for the Euclidean norm of $x$. For $d\in\mathbb{Z}_{+}$,
we use $I_{d}$ to denote the $d\times d$ identity matrix. For matrix $A$, we
use $\left\|{A}\right\|$ and $\left\|{A}\right\|_{F}$ to denote the spectral
norm and Frobenius norm of $A$, respectively. We use $\lambda_{\max}(A)$ and
$\lambda_{\min}(A)$ to denote the maximum and minimum eigenvalue of $A$,
respectively, while $\operatorname{tr}(A)$ denotes the trace. For square
matrices $A$ and $B$ of the same dimension, We use $\langle
A,B\rangle=\operatorname{tr}(A^{\top}B)$ to denote the entry-wise inner
product.

######  Definition 1.1 (Lipschitz continuity).

A function $h:\mathcal{X}\rightarrow\mathcal{Y}$ is $\ell$-Lipschitz if

$\left\|{h(x_{1})-h(x_{2})}\right\|_{\mathcal{X}}\leq\ell\left\|{x_{1}-x_{2}}\right\|_{\mathcal{Y}}\quad\mbox{for
all $x_{1},x_{2}\in\mathcal{X}$}.$

For purposes of this definition, we endow Euclidean spaces with the Euclidean
norm. When $\mathcal{X}$ is a matrix space, we endow it with the Frobenius
norm. When $\mathcal{Y}$ is a matrix space, we endow it with the spectral
norm.

##  2 A randomized algorithm with expected descent

Here we define our algorithm (Algorithm 1) after stating our assumptions on
$f$ and on the approximate gradients and Hessians. We follow in Section 2.2
with a bound of the expected number of iterations required for the algorithm
to terminate. Section 2.3 describes a high-probability bound, that is, a bound
on the number of iterations of the algorithm required to stop at an
approximate second-order point with probability at least $1-\delta$, for a
specified (small, positive) value of $\delta$.

###  2.1 Algorithm and Assumptions

We begin with our assumptions on $f$ and on the approximate gradients and
Hessians used by the algorithm.

######  Assumption 2.1.

The function $f$ has a gradient $\nabla f$ that is $L$-Lipschitz continuous at
all points of interest and a Hessian that is $M$-Lipschitz continuous, and is
uniformly bounded below by a quantity $\bar{f}$.

Under this assumption, the following bounds on $f$ are well known:

######  Fact 2.2.

Assume $f$ satisfies Assumption 2.1. For any $x,d,p\in\mathbb{R}^{d}$ and
$\alpha\in\mathbb{R}$, we have

$\displaystyle f(x+\alpha d)$ $\displaystyle\leq f(x)+\alpha\nabla
f(x)^{T}d+\alpha^{2}\frac{L}{2}\|d\|^{2}$ (3) $\displaystyle f(x+p)$
$\displaystyle\leq f(x)+\nabla
f(x)^{T}p+\frac{1}{2}p^{T}\nabla^{2}f(x)p+\frac{1}{6}M\|p\|^{3}$ (4)

We suppose that when our algorithm requests a gradient evaluation at the
iterate $x_{k}$, it receives the approximate value $g_{k}$, and if it requests
a Hessian, it receives the inexact value $H_{k}$. We make the following
assumption on the inexactness in these quantities.

######  Assumption 2.3.

At all iterations $k$, the errors in the gradient and Hessian oracles satisfy

$\|g_{k}-\nabla{f}(x_{k})\|\leq\frac{1}{3}\max\\{\epsilon_{g},\|g_{k}\|\\},\quad\|H_{k}-\nabla^{2}f(x_{k})\|\leq\frac{2}{9}\epsilon_{H}.$

We could tolerate a greater degree of Hessian inexactness, specifically,
$\|H_{k}-\nabla^{2}f(x_{k})\|\leq\frac{2}{9}\max\\{\epsilon_{H},|\lambda_{\min}(H_{k})|\\}$.
However, since Hessian steps are taken infrequently in practice, this looser
tolerance does not have great practical consequence, so we opt for the simpler
analysis afforded by the bound of Assumption 2.3.

Our algorithm requires an approximate direction of minimal curvature for the
approximate Hessian $H_{k}$. Allowing for further error in the computation of
this unit direction and the associated eigenvalue, we make an additional
assumption.

######  Assumption 2.4.

Let $\hat{p}_{k}$ be an approximate unit eigenvector of $H_{k}$ corresponding
to the approximate minimum eigenvalue $\hat{\lambda}_{k}$, that is,
$\hat{\lambda}_{k}:=\hat{p}_{k}^{T}H_{k}\hat{p}_{k}$. We assume that

$|\hat{\lambda}_{k}-\lambda_{\min}(H_{k})|\leq\frac{1}{9}\epsilon_{H}.$

This assumption enables practical computation of minimum eigenvalue and its
corresponding eigenvector via only matrix-vector products. [ROW20, Appendix B]
used Lanczos and conjugate gradient to construct a minimum eigenvalue oracle
satisfying Assumption 2.4.

We are ready to define our algorithm, which adds Rademacher randomness in the
negative curvature steps to the second-order algorithm in [WR22, Section 3.6]
and allows gradient and Hessian inexactness.

1 Given target tolerance $\epsilon_{g}$, $\epsilon_{H}$, Lipschitzness
constants $L$, $M$, upper bound $\alpha\in[\epsilon_{H},L]$ on negative
curvature stepsizes, and initialization $x_{0}$

2 for  _$k=0,1,2,\dotsc$_ do

3 if  _$\|g_{k}\| >\epsilon_{g}$_ then

/* gradient step */

4 $x_{k+1}=x_{k}-\frac{1}{L}g_{k}$

5 else if  _$\hat{\lambda}_{k} <-\epsilon_{H}$_ then

/* negative curvature step */

6 Choose $\alpha_{k}\in[\epsilon_{H},\min(\alpha,|\hat{\lambda}_{k}|)]$

7 Draw $\sigma_{k}\leftarrow\pm 1$ with probability $\frac{1}{2}$

8 $x_{k+1}=x_{k}+\frac{2\alpha_{k}}{M}\sigma_{k}\hat{p}_{k}$

9 else

10 return $x_{k}$

11

Algorithm 1 Algorithm using approximate gradients and Hessians to find
$(\epsilon_{g},\epsilon_{H})$-approximate second-order points

We observe the following output guarantee in Algorithm 1:

######  Proposition 2.5.

If Algorithm 1 terminates and returns $x_{n}$, then $x_{n}$ is an
$(\frac{4}{3}\epsilon_{g},\frac{4}{3}\epsilon_{H})$ approximate second-order
stationary point.

###### Proof.

From the stopping criteria, we have $\left\|{g_{n}}\right\|\leq\epsilon_{g}$
and $\hat{\lambda}_{n}\geq-\epsilon_{H}$. Thus,

$\left\|{\nabla
f(x_{n})}\right\|\leq\left\|{g_{n}}\right\|+\frac{1}{3}\max\\{\epsilon_{g},\|g_{n}\|\\}\leq\epsilon_{g}+\frac{1}{3}\epsilon_{g}=\frac{4}{3}\epsilon_{g}.$

To bound $\lambda_{\min}(\nabla^{2}f(x_{n}))$, write
$U_{n}=H_{n}-\nabla^{2}f(x_{n})$. By Assumption 2.3,
$\lambda_{\min}(-U_{n})\geq-\left\|{-U_{n}}\right\|\geq-\frac{2}{9}\epsilon_{H}$.
By Assumption 2.4,
$\lambda_{\min}H_{n}\geq\hat{\lambda}_{n}-\frac{1}{9}\epsilon_{H}\geq-\epsilon_{H}-\frac{1}{9}\epsilon_{H}\geq-\frac{10}{9}\epsilon_{H}$.

We use Weyl’s theorem to conclude that

$\lambda_{\min}(\nabla^{2}f(x_{n}))\geq\lambda_{\min}(H_{n})+\lambda_{\min}(-U_{n})\geq-\frac{10}{9}\epsilon_{H}-\frac{2}{9}\epsilon_{H}=-\frac{4}{3}\epsilon_{H},$

as required. ∎

Note that this output guarantee only applies if Algorithm 1 terminates. We
proceed to derive the iteration complexity of this algorithm. The proof
technique is to find a lower bound on the amount of decrease per iteration. A
complexity result then follows from the fact that $f$ is bounded from below by
$\bar{f}$.

Write $v_{k}:=g_{k}-\nabla f(x_{k})$. When the algorithm takes a gradient
step, we have $\|g_{k}\|>\epsilon_{g}$. By Assumption 2.3, we have
$\|v_{k}\|\leq\frac{1}{3}\|g_{k}\|$, so that $\|\nabla
f(x_{k})\|\geq\frac{2}{3}\|g_{k}\|$. From (3), we have

$\displaystyle f(x_{k+1})$
$\displaystyle=f\left(x_{k}-\frac{1}{L}g_{k}\right)$ $\displaystyle\leq
f(x_{k})-\frac{1}{L}\nabla
f(x_{k})^{T}g_{k}+\frac{L}{2}\cdot\frac{1}{L^{2}}\|g_{k}\|^{2}$
$\displaystyle=f(x_{k})-\frac{1}{2L}\left(\|\nabla
f(x_{k})\|^{2}-\|v_{k}\|^{2}\right)$ $\displaystyle\leq
f(x_{k})-\frac{1}{2L}\left(\frac{4}{9}\|g_{k}\|^{2}-\frac{1}{9}\|g_{k}\|^{2}\right)$
$\displaystyle\leq f(x_{k})-\frac{1}{6L}\|g_{k}\|^{2}$ $\displaystyle\leq
f(x_{k})-\frac{1}{6L}\epsilon_{g}^{2},$ (5)

For the negative curvature step, recall that
$\epsilon_{H}\leq\alpha_{k}\leq\alpha\leq|\hat{\lambda}_{k}|\leq L$ and
$\|H_{k}-\nabla^{2}f(x_{k})\|\leq\frac{2}{9}\max\\{\epsilon_{H},|\hat{\lambda}_{k}|\\}$.
It follows from (4) that

$\displaystyle f(x_{k+1})$
$\displaystyle=f(x_{k}+\frac{2\alpha_{k}}{M}\sigma_{k}\hat{p}_{k})$
$\displaystyle\leq f(x_{k})+2\frac{\alpha_{k}}{M}\nabla
f(x_{k})^{T}\sigma_{k}\hat{p}_{k}+\frac{1}{2}\frac{4\alpha_{k}^{2}}{M^{2}}\hat{p}_{k}^{T}\nabla^{2}f(x_{k})\hat{p}_{k}+\frac{M}{6}\frac{8\alpha_{k}^{3}}{M^{3}}$
$\displaystyle=f(x_{k})+\frac{2\alpha_{k}^{2}}{M^{2}}\left(\hat{p}_{k}^{T}\nabla^{2}f(x_{k})\hat{p}_{k}+\frac{2}{3}\alpha_{k}\right)+2\frac{\alpha_{k}}{M}\nabla
f(x_{k})^{T}\sigma_{k}\hat{p}_{k}$

When the algorithm takes a negative curvature step, we have
$\hat{\lambda}_{k}<-\epsilon_{H}<0$, so by Assumption 2.3, we have
$\|H_{k}-\nabla^{2}f(x_{k})\|\leq\frac{2}{9}\epsilon_{H}\leq\frac{2}{9}|\hat{\lambda}_{k}|$.
It follows from the definition of operator norm and Cauchy-Schwarz that
$|\hat{p}_{k}^{T}H_{k}\hat{p}_{k}-\hat{p}_{k}^{T}\nabla^{2}f(x_{k})\hat{p}_{k}|\leq\|H_{k}\hat{p}_{k}-\nabla^{2}f(x_{k})\hat{p}_{k}\|\leq\frac{2|\hat{\lambda}_{k}|}{9}$,
so

$\hat{p}_{k}^{T}\nabla^{2}f(x_{k})\hat{p}_{k}\leq\frac{2}{9}|\hat{\lambda}_{k}|+{p}_{k}^{T}H_{k}\hat{p}_{k}=\frac{2}{9}|\hat{\lambda}_{k}|+(-|\hat{\lambda}_{k}|)=-\frac{7}{9}|\hat{\lambda}_{k}|.$

We thus have

$\displaystyle f(x_{k+1})$ $\displaystyle\leq
f(x_{k})+\frac{2\alpha_{k}^{2}}{M^{2}}\left(\hat{p}_{k}^{T}\nabla^{2}f(x_{k})\hat{p}_{k}+\frac{2}{3}\alpha_{k}\right)+2\frac{\alpha_{k}}{M}\nabla
f(x_{k})^{T}\sigma_{k}\hat{p}_{k}$ $\displaystyle\leq
f(x_{k})+\frac{2\alpha_{k}^{2}}{M^{2}}\left(-\frac{7}{9}|\hat{\lambda}_{k}|+\frac{2}{3}\alpha_{k}\right)+2\frac{\alpha_{k}}{M}\nabla
f(x_{k})^{T}\sigma_{k}\hat{p}_{k}$ $\displaystyle\leq
f(x_{k})-\frac{2\epsilon_{H}^{3}}{9M^{2}}+2\frac{\alpha_{k}}{M}\nabla
f(x_{k})^{T}\sigma_{k}\hat{p}_{k},$ (6)

where the final inequality follows from our assumption
$\alpha_{k}\in[\epsilon_{H},|\hat{\lambda}_{k}|]$.

###  2.2 Bound on expected stopping time

We have the following result for expected stopping time of Algorithm 1. Here
the expectation is taken with respect to the random variables $\sigma_{k}$
used at the negative curvature iterations. For purposes of this and later
results, we define

$C_{\epsilon}:=\min\left(\frac{\epsilon_{g}^{2}}{6L},\frac{2\epsilon_{H}^{3}}{9M^{2}}\right).$
(7)

######  Theorem 2.6.

Suppose the function $f$ satisfies Assumption 2.1, its inexact gradient and
Hessian computations satisfy Assumption 2.3, and minimum eigenvalue
computation satisfies Assumption 2.4. Let $T$ denote the iteration at which
Algorithm 1 terminates. Then $T<\infty$ almost surely and

$\mathbb{E}T\leq\frac{f(x_{0})-\bar{f}}{C_{\epsilon}},$ (8)

where $C_{\epsilon}$ is defined in (7).

The proof of this result is given below. It constructs a supermartingale 222A
supermartingale with respect to filtration
$\\{\mathcal{G}_{1},\mathcal{G}_{2},\dotsc\\}$ is a sequence of random
variables $\\{Y_{1},Y_{2},\dotsc\\}$ such that for all $k\in\mathbb{Z}_{+}$,
(i) $\mathbb{E}\left|{Y_{k}}\right|<\infty$, (ii) $Y_{k}$ is
$\mathcal{G}_{k}$-measurable, and (iii)
$\mathbb{E}(Y_{k+1}\,|\mathcal{G}_{k})\leq Y_{k}$. based on the function value
and uses a supermartingale convergence theorem and optional stopping theorem
to obtain the final result. A similar proof technique is used in [Ber+21] but
for a line-search algorithm. We collect several relevant facts about
supermartingales before proving the result.

First, we need to ensure the relevant supermartingale is well defined even
after the algorithm terminates, so that it is possible to let the index $k$ of
the supermartingale go to $\infty$.

######  Fact 2.7 ([Dur19, Theorem 4.2.9]).

If $T$ is a stopping time and $X_{k}$ is a supermartingale, then
$X_{\min(T,k)}$ is a supermartingale.

The following supermatingale convergence theorem will be used to ensure the
function value converges, so that the algorithm terminates with probability 1.

######  Fact 2.8 (Supermartingale convergence theorem, [Dur19, Theorem
4.2.12]).

If $X_{k}\geq 0$ is a supermartingale, then as $k\rightarrow\infty$, there
exists a random variable $X$ such that $X_{k}\rightarrow X\;a.s.$ and
$\mathbb{E}X\leq\mathbb{E}X_{0}$.

Finally, we will use optional stopping theorem to derive the expected
iteration complexity. Note that we use a version of the optional stopping
theorem specific to nonnegative supermartingales that does not require uniform
integrability.

######  Fact 2.9 (Optional stopping theorem, [Dur19, Theorem 4.8.4]).

If $X_{n}$ is a nonnegative supermartingale and $N\leq\infty$ is a stopping
time, then $\mathbb{E}X_{0}\geq\mathbb{E}X_{N}$

###### Proof of Theorem 2.6.

We first construct a supermartingale based on function values. Since
$\mathbb{E}\left[\sigma_{k}\right]=0$, linearity of expectation implies that
$\mathbb{E}\left[2\frac{\alpha_{k}}{M}\nabla
f(x_{k})^{T}\sigma_{k}\hat{p}_{k}\middle|x_{k}\right]=0$. We therefore have
from (6) that

$\displaystyle\mathbb{E}\left[f(x_{k+1})\middle|x_{k}\right]$
$\displaystyle\leq f(x_{k})-\frac{2\epsilon_{H}^{3}}{9M^{2}}.$

By combining with the (deterministic) first-order decrease estimate (5), we
have

$\mathbb{E}\left[f(x_{k+1})\middle|x_{k}\right]\leq
f(x_{k})-\min\left(\frac{\epsilon_{g}^{2}}{6L},\frac{2\epsilon_{H}^{3}}{9M^{2}}\right)=f(x_{k})-C_{\epsilon}.$

Consider the stochastic process $M_{k}:=f(x_{k})+kC_{\epsilon}$. We have

$\displaystyle\mathbb{E}\left[M_{k+1}\middle|x_{k}\right]$
$\displaystyle=\mathbb{E}\left[f(x_{k+1})+(k+1)C_{\epsilon}\middle|x_{k}\right]$
$\displaystyle\leq\mathbb{E}\left[f(x_{k})-C_{\epsilon}+(k+1)C_{\epsilon}\middle|x_{k}\right]$
$\displaystyle=\mathbb{E}\left[f(x_{k})+kC_{\epsilon}\middle|x_{k}\right]=M_{k}.$

We need to select a filtration to define the supermartingale $M_{k}$. We view
$x_{k}$ as random variables defined with respect to $\sigma_{i},i\leq k$.
Since $M_{k}$ is expressed as a function of $x_{k}$ only, we define the
filtration $\\{\mathcal{G}_{k}\\}$ to be the filtration generated by $x_{k}$,
and it naturally holds that $M_{k}$ is $\mathcal{G}_{k}$-measurable for all
$k$ and
$\mathbb{E}\left[M_{k+1}\middle|\mathcal{G}_{k}\right]=\mathbb{E}\left[M_{k+1}\middle|x_{k}\right]$.
Hence, $\\{M_{k}\\}$ is a supermartingale with respect to filtration
$\\{\mathcal{G}_{k}\\}$.

Let $T$ denote the iteration at which our algorithm stops. Since the decision
to stop at iteration $k$ depends only on $x_{k}$, we have
$\\{T=k\\}\in\mathcal{G}_{k}$, which implies $T$ is a stopping time.

We will use the supermartingale convergence theorem (Fact 2.8) to show that
$T<+\infty$ almost surely, since the function value cannot decrease
indefinitely as it is bounded by $\bar{f}$ from below. To apply Fact 2.8, we
need to let $k\rightarrow\infty$, so we need to transform $\\{M_{k}\\}$ to
obtain a supermartingale $\\{Y_{k}\\}$ that is well defined even after the
algorithm terminates.

It follows from Fact 2.7 that $Y_{k}:=M_{\min(k,T)}$ is also a
supermartingale. Since $Y_{k}\geq\bar{f}$, it follows from the supermartingale
convergence theorem (Fact 2.8) applied to $Y_{k}-\bar{f}$ that
$Y_{k}\rightarrow Y_{\infty}$ almost surely for some random variable
$Y_{\infty}$ with
$\mathbb{E}Y_{\infty}\leq\mathbb{E}Y_{0}=\mathbb{E}M_{0}=f(x_{0})<\infty$.
Hence $\mathbb{P}[Y_{\infty}=+\infty]=0$. On the other hand, as
$k\rightarrow\infty$, we have $kC_{\epsilon}\rightarrow\infty$, so
$T=+\infty\implies
Y_{k}=M_{k}\geq\bar{f}+kC_{\epsilon}\rightarrow\infty\implies
Y_{\infty}=+\infty$. Therefore we have $\mathbb{P}[T<+\infty]=1$.

We can then apply the optional stopping theorem (Fact 2.9) to $Y_{k}-\bar{f}$.
It follows that

$\bar{f}+\mathbb{E}T\cdot C_{\epsilon}\leq\mathbb{E}f(x_{T})+\mathbb{E}T\cdot
C_{\epsilon}=\mathbb{E}\left[M_{T}\right]=\mathbb{E}\left[Y_{T}\right]\leq\mathbb{E}[Y_{0}]=\mathbb{E}\left[M_{0}\right]=f(x_{0}),$

where the first equality uses $T<+\infty$ almost surely and the last
inequality is Fact 2.9. By reorganizing this bound, we obtain
$\mathbb{E}T\leq\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}$, as desired. ∎

######  Remark 2.10.

The algorithm in [WR22, Section 3.6] also finds an
$(\epsilon_{g},\epsilon_{H})$-approximate second-order stationary points
without allowing for gradient and Hessian inexactness. Consequently, the
algorithm there is deterministic; it does not need to introduce Rademacher
randomness in the negative curvature step. The iteration complexity for that
algorithm is the same as the bound on $\mathbb{E}T$ in Theorem 2.6 to within a
constant factor.

In summary, the expected iteration complexity of Algorithm 1 matches that of
the corresponding simple deterministic algorithm up to a constant factor, but
it has the advantage of allowing for gradient and Hessian inexactness.

###  2.3 High-probability bound on stopping time

In this section, we state and prove our main result concerning the number of
iterations required for Algorithm 1 to terminate, with probability $1-\delta$,
for some small positive $\delta$. We start by stating the main result, then
state several immediate corollaries concerning various choices of the
parameters $\epsilon_{g}$, $\epsilon_{H}$, and $\alpha_{k}$.

Recall that the gradient and negative curvature steps in Algorithm 1 have the
form

$\displaystyle x_{k+1}$ $\displaystyle=x_{k}-\frac{1}{L}g_{k},$ (9)
$\displaystyle x_{k+1}$
$\displaystyle=x_{k}+\frac{2\alpha_{k}}{M}\sigma_{k}\hat{p}_{k},$ (10)

respectively, where $\sigma_{k}=\pm 1$ with probability 1/2, and
$\alpha_{k}\in[\epsilon_{H},\min(\alpha,|\hat{\lambda}_{k}|)]$.

######  Theorem 2.11.

Suppose the function $f$ satisfies Assumption 2.1, its inexact gradient and
Hessian computations satisfy Assumption 2.3, and the minimum eigenvalue /
eigenvector computation satisfies Assumption 2.4. Then Algorithm 1 terminates
after $n=2\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)+BK$ iterations
with probability at least $1-\delta$, where $C_{\epsilon}$ is defined in (7)
and

$\displaystyle K$
$\displaystyle=\max\left\\{C\log\left(\frac{1}{\delta}\right),4\eta
C\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)^{\frac{1}{\eta}},4\eta^{2}C^{1+\frac{1}{\eta-1}}B^{\frac{1}{\eta-1}}\right\\}$
$\displaystyle B$ $\displaystyle=1+\frac{18\alpha L}{M\epsilon_{g}}$
$\displaystyle C$
$\displaystyle=2304\frac{M^{2}\alpha^{2}\epsilon_{g}^{2}}{\epsilon_{H}^{6}}$
$\displaystyle\eta$ $\displaystyle\geq 2\text{ is an integer we can choose
arbitrarily.}$

The proof below will reveal the motivation behind the choice of constants
$K,B$ and $C$. To write the result more compactly, we can say

$n=2\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)+BC\max\left\\{\log\left(\frac{1}{\delta}\right),4\eta\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)^{\frac{1}{\eta}},4\eta^{2}\left(BC\right)^{\frac{1}{\eta-1}}\right\\}.$
(11)

We consider several choices of $\epsilon_{g}$ and $\epsilon_{H}$ and find
their corresponding iteration complexity from this result. In the following
statements, the expressions $O\left({\cdot}\right)$ and
$\tilde{O}\left({\cdot}\right)$ hide dependence on the Lipschitz constants for
gradient ($L$) and Hessian ($M$). In each case we find conditions under which
the iteration complexity is approximately
$n=\tilde{O}(\frac{1}{\epsilon_{g}^{2}})$, which is within a logarithmic
factor of the complexity of a gradient descent algorithm for finding an point
satisfying approximate first-order optimality conditions.

######  Corollary 2.12.

When $\epsilon_{g}=\epsilon,\epsilon_{H}=\sqrt{\epsilon M}$, and
$\alpha_{k}=\epsilon_{H}$ for all $k$, we can define $\alpha=\epsilon_{H}$ in
Algorithm 1. Then $C=O(1)$ and $B=O(\frac{1}{\epsilon^{1/2}})$. Choosing
$\eta=2$ gives $K=\tilde{O}(\frac{1}{\epsilon})$ and
$n=\tilde{O}(\frac{1}{\epsilon^{2}})$.

######  Corollary 2.13.

When $\epsilon_{g}=\epsilon,\epsilon_{H}=\sqrt{\epsilon M}$, and
$\alpha_{k}=|\hat{\lambda}_{k}|$, we can take $\alpha=2L$. Then
$B=O(\frac{1}{\epsilon})$ and $C=O(\frac{1}{\epsilon})$. Choosing $\eta=11$
gives $K=\tilde{O}(\frac{1}{\epsilon^{1.2}})$ and
$n=\tilde{O}(\frac{1}{\epsilon^{2.2}})$. By choosing a larger value of $\eta$
we can move the final bound closer to $n=\tilde{O}(\frac{1}{\epsilon^{2}})$ at
the expense of a larger constant factor.

######  Corollary 2.14.

When $\epsilon_{g}$ and $\epsilon_{H}$ satisfy
$\frac{\epsilon_{g}^{2}}{6L}=\frac{2\epsilon_{H}^{3}}{9M^{2}}$, we can set
$\epsilon=\epsilon_{g}$ and take $\alpha_{k}=\alpha=\epsilon_{H}$ to obtain
$C=O(\frac{1}{\epsilon_{H}})=O\left(\frac{1}{\epsilon^{2/3}}\right)$ and
$B=O\left(\frac{1}{\epsilon^{1/3}}\right)$. Choosing $\eta=2$ gives
$n=\tilde{O}(\frac{1}{\epsilon^{2}})$.

######  Remark 2.15.

Theorem 2.11 and its analysis requires that Assumption 2.3 holds for all
iterations. If the probabilty that Assumption 2.3 fails on any given iteration
is at most $\xi$, we can use the union bound to ensure the probability that
Assumption 2.3 holds for all $n$ iterations to be at least $1-n\xi$. If the
failure probability tolerance for both types of failure (either that Algorithm
1 doesn’t terminate or that Assumption 2.3 fails at some iterations) is
$\Theta(\delta)$ in total, then we require $\xi=O(\delta/n)$. Since
$n=O\left(\frac{1}{\epsilon^{2}}\log\left(\frac{1}{\delta}\right)\right)=O\left(\frac{1}{\epsilon^{2}\delta}\right)$,
the condition $\xi=O(\delta^{2}\epsilon^{2})$ ensures a total failure
probability of $\Theta(\delta)$. We will elaborate on this point in Section 4
when we discuss using mini-batch sampling to achieve Assumption 2.3.

###### Proof of Theorem 2.11.

For the $n$ given in Theorem 2.11 and (11), define the event

$E:=\\{\mbox{Algorithm does not terminate before step }n\\}.$ (12)

We wish to show $\mathbb{P}(E)\leq\delta$.

In the discussions below, we use notation $\\{k\in NC\\}$ and $\\{k\in GD\\}$
to indicate the event that $k$-th step (from $x_{k}$ to $x_{k+1}$) is a
negative curvature step and a gradient descent step, respectively. Let $|NC|$
be the number of negative curvature steps the algorithm has taken before
termination, and let $|GD|$ be the number of total gradient descent steps.
Note that $|NC|$ and $|GD|$ are both random variables whose values depend on
the Rademacher randomness $\sigma_{k}$.

We combine the inequalities (5) and (6) to obtain the descent from both
gradient descent and negative curvature steps:

$\displaystyle f(x_{k+1})$ $\displaystyle\leq f(x_{k})+\mathbf{1}_{\\{k\in
NC\\}}\left(-\frac{2\epsilon_{H}^{3}}{9M^{2}}+\frac{2\alpha_{k}}{M}\sigma_{k}\nabla
f(x_{k})^{T}\hat{p}_{k}\right)+\mathbf{1}_{\\{k\in
GD\\}}\left(-\frac{\epsilon_{g}^{2}}{6L}\right).$

By telescoping the above inequality, we obtain

$\displaystyle\begin{split}\bar{f}-f(x_{0})&\leq f(x_{n})-f(x_{0})\\\
&\leq\sum_{k=0}^{n-1}\mathbf{1}_{\\{k\in
NC\\}}\left(-\frac{2\epsilon_{H}^{3}}{9M^{2}}+\frac{2\alpha_{k}}{M}\sigma_{k}\nabla
f(x_{k})^{T}\hat{p}_{k}\right)+\mathbf{1}_{\\{k\in
GD\\}}\left(-\frac{\epsilon_{g}^{2}}{6L}\right)\\\
&\leq|NC|\frac{2\alpha_{k}\epsilon_{g}}{{M}}\left(-\frac{\epsilon_{H}^{3}}{9M\alpha_{k}\epsilon_{g}}+\frac{1}{|NC|}\sum_{k\in
NC}\frac{\nabla
f(x_{k})^{T}\hat{p}_{k}}{\epsilon_{g}}\sigma_{k}\right)+|GD|\left(-\frac{\epsilon_{g}^{2}}{6L}\right)\end{split}$
(13)

We first prove two claims that help us understand the event $E$ defined in
(12) better. The first claim implies that $E\subset\\{|NC|\geq K\\}$.

######  Claim 2.16.

For all $K\in\mathbb{Z}_{+}$, if the algorithm does not terminate within
$n=2\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)+BK$ steps (that is,
event $E$ is realized), then we have $|NC|\geq K$ with probability 1.

###### Proof of Claim 2.16.

Assume to the contrary that $|NC|<K$. It follows from
$|GD|+|NC|=n=2\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)+BK$ that
$|GD|>2\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}+(B-1)K>\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}+\frac{18\alpha
L}{M\epsilon_{g}}K$. Then

$\displaystyle f(x_{0})-f(x_{n})$
$\displaystyle\geq\sum_{k=0}^{n-1}\mathbf{1}_{\\{k\in
NC\\}}\left(\frac{2\epsilon_{H}^{3}}{9M^{2}}-\frac{2\alpha_{k}}{M}\sigma_{k}\nabla
f(x_{k})^{T}\hat{p}_{k}\right)+\mathbf{1}_{\\{k\in
GD\\}}\left(\frac{\epsilon_{g}^{2}}{6L}\right)$
$\displaystyle\geq|NC|\left(\frac{2\epsilon_{H}^{3}}{9M^{2}}-\frac{2\alpha}{M}\|\nabla
f(x_{k})\|\right)+|GD|\frac{\epsilon_{g}^{2}}{6L}$
$\displaystyle\geq|NC|\left(-\frac{2\alpha}{M}\frac{4}{3}\epsilon_{g}\right)+|GD|\frac{\epsilon_{g}^{2}}{6L}$
$\displaystyle>K\left(-\frac{8\alpha\epsilon_{g}}{3M}\right)+\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}+\frac{18\alpha
L}{M\epsilon_{g}}K\right)\frac{\epsilon_{g}^{2}}{6L}$
$\displaystyle=K\left(\frac{\epsilon_{g}\alpha}{3M}\right)+\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\frac{\epsilon_{g}^{2}}{6L}$
$\displaystyle>\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\frac{\epsilon_{g}^{2}}{6L}\geq
f(x_{0})-\bar{f},$

where the second line is Cauchy-Schwarz, the third line uses the facts that
$\|\nabla f(x_{k})\|\leq\tfrac{4}{3}\epsilon_{g}$ for iterations $k\in NC$
(see Assumption 2.3 and Algorithm 1), the fourth line uses the assumption on
$|NC|$ and $|GD|$, and the final inequality is because
$C_{\epsilon}\leq\frac{\epsilon^{2}}{6L}$ by definition. This inequality
contradicts the lower bound of $f$, so we must have $|NC|\geq K$, as claimed.
∎

The next claim concerns the size of the “$NC$” term in (13) .

######  Claim 2.17.

Assume that the algorithm does not terminate within
$n=2\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)+BK$ steps. Then in the
coefficient of the “$NC$” term in (13), we have

$-\frac{\epsilon_{H}^{3}}{9M\alpha_{k}\epsilon_{g}}+\frac{1}{|NC|}\sum_{k\in
NC}\frac{\nabla
f(x_{k})^{T}\hat{p}_{k}}{\epsilon_{g}}\sigma_{k}>-\frac{\epsilon_{H}^{3}}{18M\alpha_{k}\epsilon_{g}}$

###### Proof of Claim 2.17.

Assume to the contrary that

$-\frac{\epsilon_{H}^{3}}{9M\alpha_{k}\epsilon_{g}}+\frac{1}{|NC|}\sum_{k\in
NC}\frac{\nabla
f(x_{k})^{T}\hat{p}_{k}}{\epsilon_{g}}\sigma_{k}\leq-\frac{\epsilon_{H}^{3}}{18M\alpha_{k}\epsilon_{g}}.$

Equation (13) becomes

$\displaystyle\bar{f}-f(x_{0})$ $\displaystyle\leq f(x_{n})-f(x_{0})$
$\displaystyle\leq|NC|\frac{2\alpha_{k}\epsilon_{g}}{{M}}\left(-\frac{\epsilon_{H}^{3}}{9M\alpha_{k}\epsilon_{g}}+\frac{1}{|NC|}\sum_{k\in
NC}\frac{\nabla
f(x_{k})^{T}\hat{p}_{k}}{\epsilon_{g}}\sigma_{k}\right)+|GD|\left(-\frac{\epsilon_{g}^{2}}{6L}\right)$
$\displaystyle\leq|NC|\frac{2\alpha_{k}\epsilon_{g}}{{M}}\left(-\frac{\epsilon_{H}^{3}}{18M\alpha_{k}\epsilon_{g}}\right)+|GD|\left(-\frac{\epsilon_{g}^{2}}{6L}\right)$
$\displaystyle=|NC|\left(-\frac{\epsilon_{H}^{3}}{9M^{2}}\right)+|GD|\left(-\frac{\epsilon_{g}^{2}}{6L}\right)$
$\displaystyle\leq|NC|\left(-\frac{\epsilon_{H}^{3}}{9M^{2}}\right)+|GD|\left(-\frac{\epsilon_{g}^{2}}{12L}\right)$
$\displaystyle\leq-n\left(\frac{1}{2}C_{\epsilon}\right).$

where the final inequality follows from the definition of $C_{\epsilon}$ in
(7) and the fact that $|GD|+|NC|=n$. We thus have $n\leq
2\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)$, which contradicts
$n=2\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)+BK$ since $B>1$ and
$K\geq 1$. ∎

For the event $E$ defined in (12), Claim 2.17 implies that

$E\subset\left\\{\frac{1}{|NC|}\sum_{k\in NC}\frac{\nabla
f(x_{k})^{T}\hat{p}_{k}}{\epsilon_{g}}\sigma_{k}>\frac{\epsilon_{H}^{3}}{18M\alpha_{k}\epsilon_{g}}\right\\}\subset\left\\{\frac{1}{|NC|}\sum_{k\in
NC}\frac{\nabla
f(x_{k})^{T}\hat{p}_{k}}{\epsilon_{g}}\sigma_{k}>\frac{\epsilon_{H}^{3}}{18M\alpha\epsilon_{g}}\right\\},$

where the second “$\subset$” is due to $\alpha_{k}\leq\alpha$. By combining
Claims 2.16 and 2.17, we obtain $E\subset F$ where

$F:=\left\\{\frac{1}{|NC|}\sum_{k\in NC}\frac{\nabla
f(x_{k})^{T}\hat{p}_{k}}{\epsilon_{g}}\sigma_{k}>\frac{\epsilon_{H}^{3}}{18M\alpha\epsilon_{g}}\right\\}\cap\\{|NC|\geq
K\\}.$

To show $\mathbb{P}(E)\leq\delta$, it suffices to show
$\mathbb{P}\left(F\right)\leq\delta$.

Let the subsequence $k_{i}$ enumerate $\\{k:k\in NC\\}$, that is, $\\{k:k\in
NC\\}=\\{k_{1},\dots,k_{|NC|}\\}$. (Formally, we can consider $k_{i}$’s as a
sequence of random variables that take non-negative integer values). We define

$\tau_{i}:=\frac{\nabla
f(x_{k_{i}})^{T}\hat{p}_{k_{i}}\sigma_{k_{i}}}{\epsilon_{g}},$ (14)

and extend the definition of $\tau_{i}$ for $i>|NC|$ by defining $\tau_{i}=0$
for all such $i$.

Observe that realizations of the random variable $|NC|$ on event $\\{|NC|\geq
K\\}$ can be $\\{K,K+1,\cdots,n\\}$. Therefore, we have

$\displaystyle F$ $\displaystyle\subset\left\\{\exists
s=K,K+1,\dots,n,\;\frac{1}{s}\sum_{i=1}^{s}\frac{\nabla
f(x_{k_{i}})^{T}\hat{p}_{k_{i}}}{\epsilon_{g}}\sigma_{k_{i}}>\frac{\epsilon_{H}^{3}}{18M\alpha\epsilon_{g}}\right\\}\cap\\{|NC|\geq
K\\}$
$\displaystyle=\cup_{s=K}^{n}\left\\{\frac{1}{s}\sum_{i=1}^{s}\tau_{{i}}>\frac{\epsilon_{H}^{3}}{18M\alpha\epsilon_{g}}\right\\}\cap\\{|NC|\geq
K\\}$
$\displaystyle\subset\cup_{s=K}^{n}\left\\{\frac{1}{s}\sum_{i=1}^{s}\tau_{{i}}>\frac{\epsilon_{H}^{3}}{18M\alpha\epsilon_{g}}\right\\}.$

By the union bound, to show that $\mathbb{P}(F)\leq\delta$, it suffices to
choose $K$ large enough that

$\mathbb{P}\left\\{\frac{1}{s}\sum_{i=1}^{s}\tau_{{i}}>\frac{\epsilon_{H}^{3}}{18M\alpha\epsilon_{g}}\right\\}\leq\frac{\delta}{n}\quad\mbox{for
all $s=K,K+1,\dots,n$.}$

We proceed to apply Azuma’s inequality to $X_{t}:=\sum_{i=1}^{t}\tau_{i}$ and
$X_{0}=0$. Define $\sigma$-algebra $\mathcal{F}_{t}$ to be generated by
Radamacher variables $\sigma_{k_{i}}$ for $i\leq t$ and let $\mathcal{F}_{0}$
be the trivial $\sigma$-algebra. Recall the negative curvature update is
$x_{k+1}=x_{k}+\frac{2\alpha_{k}}{M}\sigma_{k}\hat{p}_{k}$. For a given $i$,
we have that $x_{j},k_{i-1}<j\leq k_{i}$ are $\mathcal{F}_{i-1}$-measurable,
since the algorithm takes only gradient steps between $k_{i-1}$-th and
$k_{i}$-th iteration and no randomness is involved after drawing
$\sigma_{k_{i-1}}$ from the Rademacher distribution.

######  Fact 2.18 (Azuma’s inequality).

Let $\left\\{X_{0},X_{1},\cdots\right\\}$ be a martingale (or supermartingale)
with respect to filtration
$\left\\{\mathcal{F}_{0},\mathcal{F}_{1},\cdots\right\\}$. Assume there are
constants $0<c_{1},c_{2},\cdots<\infty$ such that

$-\frac{1}{2}c_{t}\leq X_{t}-X_{t-1}\leq\frac{1}{2}c_{t}$

almost surely. Then for all $\beta>0$, we have

$\mathbb{P}\left(X_{n}-X_{0}\geq\beta\right)\leq\exp\left(-\frac{2\beta^{2}}{\sum_{t=1}^{n}c_{t}^{2}}\right).$

Observe from (14) that

$\mathbb{E}[\tau_{i}|\mathcal{F}_{i-1}]=\mathbb{E}\left[\frac{\nabla
f(x_{k_{i}})\hat{p}_{k_{i}}\sigma_{k_{i}}}{\epsilon_{g}}\middle|\mathcal{F}_{i-1}\right]=\frac{\nabla
f(x_{k_{i}})\hat{p}_{k_{i}}}{\epsilon_{g}}\mathbb{E}\sigma_{k_{i}}=0,$

where the second equality is because $x_{k_{i}}$ is
$\mathcal{F}_{i-1}$-measurable. So
$\mathbb{E}[X_{t}|\mathcal{F}_{t-1}]=\mathbb{E}[X_{t-1}+\tau_{t}|\mathcal{F}_{t-1}]=X_{t-1}$,
which implies that $X_{t}$ is a martingale with respect to
$\\{\mathcal{F}_{t}\\}$.

From (14) and Assumption 2.3, we have

$\left|\tau_{i}\right|\leq\frac{\|\nabla
f(x_{k_{i}})\|}{\epsilon_{g}}\leq\frac{\|g_{k_{i}}\|+\max(\epsilon_{g},\|g_{k_{i}}\|)}{\epsilon_{g}}\leq\frac{4}{3},$

where the last inequality follows from the fact that iteration $k_{i}$ is a
negative curvature step, so that $\|g_{k_{i}}\|\leq\epsilon_{g}$. Thus Azuma’s
inequality implies that

$\displaystyle\mathbb{P}\left\\{\sum_{i=1}^{s}\tau_{i}>\frac{s\epsilon_{H}^{3}}{18M\alpha\epsilon_{g}}\right\\}\leq\exp\left(-\frac{2\left(\frac{s\epsilon_{H}^{3}}{18M\alpha\epsilon_{g}}\right)^{2}}{\left(\frac{8}{3}\right)^{2}s}\right)=\exp\left(-\frac{s\,\epsilon_{H}^{6}}{1152M^{2}\alpha^{2}\epsilon_{g}^{2}}\right)\leq\exp\left(-\frac{K\epsilon_{H}^{6}}{1152M^{2}\alpha^{2}\epsilon_{g}^{2}}\right)$

Recall $C:=2304\frac{M^{2}\alpha^{2}\epsilon_{g}^{2}}{\epsilon_{H}^{6}}$. To
make the above probability less than $\frac{\delta}{n}$, it suffices that

$K\geq\frac{C}{2}\log(n/\delta)=\frac{C}{2}(\log(A+BK)+\log(\frac{1}{\delta})),$

where $A=2\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)$. To achieve
this, we make $\frac{K}{2}\geq\frac{C}{2}\log(\frac{1}{\delta})$ and
$\frac{K}{2}\geq\frac{C}{2}(\log(A+BK)$ simultaneously true.

The first condition gives $K\geq C\log(\frac{1}{\delta})$. The second
condition becomes $K\geq C(\log(A+BK))$.

######  Fact 2.19.

Given $\eta\geq 2$. For all $x>0$, we have

$\log(x)\leq\eta(x^{\frac{1}{\eta}}-1)<\eta x^{\frac{1}{\eta}}$

Therefore, for arbitrarily chosen $\eta\geq 2$, it suffices that $K\geq
C\eta(A+BK)^{\frac{1}{\eta}}$, which holds if both $\frac{K^{\eta}}{2}\geq
C^{\eta}\eta^{\eta}A$ and $\frac{K^{\eta}}{2}\geq C^{\eta}\eta^{\eta}BK$ are
true, which hold if

$\displaystyle K$ $\displaystyle\geq 2C\eta A^{\frac{1}{\eta}}$ $\displaystyle
K$ $\displaystyle\geq
4C^{1+\frac{1}{\eta-1}}\eta^{1+\frac{1}{\eta-1}}B^{\frac{1}{\eta-1}}$

We arrive at the following sufficient condition for $K$:

$K\geq\max\left\\{C\log(\frac{1}{\delta}),4C\eta\left(\frac{f(x_{0})-\bar{f}}{C_{\epsilon}}\right)^{\frac{1}{\eta}},4\eta^{2}C^{1+\frac{1}{\eta-1}}B^{\frac{1}{\eta-1}}\right\\}$

∎

###  2.4 Computation of Inexact Eigenvectors and Eigenvalues

We include a lemma demonstrating that Assumption 2.4 can be satisfied using
matrix-vector products. This bridges the gap between iteration complexity and
operation complexity, with the latter being defined as the cumulative count of
gradient evaluations and Hessian-vector products. The underpinning algorithm,
a variant of the Lanczos algorithm, is tailored to compute the minimum
eigenvalues of a real symmetric matrix $H$ without knowledge of an upper bound
on $\left\|{H}\right\|$.

######  Lemma 2.20 ([ROW20, Lemma 10]).

Let $H$ be a real symmetric matrix. Let $\epsilon_{H}>0$ and $\delta\in(0,1)$
be given. With probability at least $1-\delta$ over uniformly random
initialization on the unit sphere, there exists an algorithm that on input $H$
(but not $\left\|{H}\right\|$) outputs a unit vector $v$ such that

$v^{\top}Hv\leq\lambda_{\min}(H)+\frac{1}{9}\epsilon_{H}$ (15)

in at most

$\min\left\\{n,1+\max\left\\{\frac{1}{2}\ln\left(25n/\delta^{2}\right),\frac{3}{2}\ln\left(25n/\delta^{2}\right)\sqrt{\frac{\|H\|}{\epsilon}}\right\\}\right\\}$
(16)

matrix-vector multiplications by $H$

In other words, the total number of gradient evaluations and Hessian-vector
products required by our algorithm is $O\left({n\epsilon_{H}^{-1/2}}\right)$,
with $n$ as specified in (11):

######  Corollary 2.21.

Suppose the function $f$ satisfies Assumption 2.1 and its inexact gradient and
Hessian computations satisfy Assumption 2.3. Then, Algorithm 1 terminates
after $n\epsilon_{H}^{-1/2}$ gradient evaluations or Hessian-vector products
with probability at least $1-\delta$, where $n$ is defined in (11). When
$\epsilon_{g}=\epsilon$, $\epsilon_{H}=\sqrt{\epsilon M}$, and the algorithm
takes a short step $\alpha=\epsilon_{H}$, the operation complexity is
$\tilde{O}\left({\frac{1}{\epsilon^{9/4}}}\right)$, where the dependence on
Lipschitz constants is suppressed.

##  3 Comparison with previous results

While Section 1.1 reviewed prior works of finding approximate SOSPs with
inexact derivatives under varied settings and assumptions, [Yao+22] is
especially relevant to this work. It addresses a similar setting with inexact
gradients and Hessians. They simply selected negative curvature step direction
$\hat{p}_{k}$ so that $g_{k}^{T}\hat{p}_{k}\leq 0$ and focused on establishing
conditions for the step size $\alpha_{k}$ in the negative curvature step to
ensure sufficient decrease. Below is a simplified analysis highlighting the
differences in our results.

In our paper, we use $\frac{1}{3}\Delta_{g,k}$ to denote the upper bound for
gradient inexactness $\|\nabla f(x_{k})-g_{k}\|$ as per Assumption 2.3. Recall
$\Delta_{g,k}=\max\\{\epsilon_{g},\|g_{k}\|\\}$. We also define
$\frac{1}{3}\delta_{g,k}$ as the tolerance for gradient inexactness in their
work. We will first scatch their analysis and then argue that their work
required a stricter condition for $\delta_{g,k}$.

From Equation (6), we can deduce that for negative curvature steps, the
descent is

$\displaystyle f(x_{k+1})$ $\displaystyle\leq
f(x_{k})+\frac{2\alpha_{k}^{2}}{M^{2}}\left(-\frac{7}{9}|\hat{\lambda}_{k}|+\frac{2}{3}\alpha_{k}\right)+2\frac{\alpha_{k}}{M}\nabla
f(x_{k})^{T}\sigma_{k}\hat{p}_{k}$ (17)

In their approach, they did not employ Rademacher randomness $\sigma_{k}$ and
instead derived a condition for $\alpha_{k}$ to ensure sufficient decrease of
$f(x_{k+1})-f(x_{k})\leq-\frac{\alpha_{k}^{3}}{9M^{2}}\leq-\frac{\epsilon_{H}^{3}}{9M^{2}}$.
Notably, this sufficient decrease condition provides a similar convergence
rate (up to a constant factor) as our expected iteration complexity. .

To ensure sufficient descent:

$\frac{2\alpha_{k}^{2}}{M^{2}}\left(-\frac{7}{9}|\hat{\lambda}_{k}|+\frac{2}{3}\alpha_{k}\right)+2\frac{\alpha_{k}}{M}\nabla
f(x_{k})^{T}\hat{p}_{k}\leq-\frac{\alpha_{k}^{3}}{9M^{2}},$

using their choice of the sign of $\hat{p}_{k}$

$g_{k}^{T}\hat{p}_{k}\leq 0\implies\nabla
f(x_{k})^{T}\hat{p}_{k}=(g_{k}+\nabla
f(x_{k})-g_{k})^{T}\hat{p}_{k}\leq\|\nabla f(x_{k})-g_{k}\|=\delta_{g,k}/3,$

they had to enforce

$7\alpha_{k}^{2}-7|\hat{\lambda}_{k}|\alpha_{k}+3M\delta_{g,k}\leq 0.$

Minimizing the quadratic function on the left,
$\alpha_{k}=\frac{1}{2}|\hat{\lambda}_{k}|$ can be chosen, which leads to the
requirement $7|\hat{\lambda}_{k}|^{2}\geq 12M\delta_{g,k}$.

Hence, the gradient inexactness requirement for the negative curvature step in
their approach is given by:

$\delta_{g,k}\leq O\left(\frac{|\lambda_{k}|^{2}}{M}\right),$

where $\lambda_{k}$ represents the most negative eigenvalue of the Hessian at
the $k$-th step. Although this requirement appears adaptive, allowing the
algorithm to accommodate higher gradient inexactness when the saddle point has
a stronger negative curvature, it’s challenging to highlight its practicality
or identify a context where gradient inexactness varies based on the Hessian’s
negative eigenvalues. Notably, in [Yao+22] they settled for:

$\delta_{g,k}\leq O\left(\frac{\epsilon_{H}^{2}}{M}\right),$

which lacks the benefits of the relative inexactness measure we propose, given
by $\Delta_{g,k}=\max\\{\epsilon_{g},\|g_{k}\|\\}$.

However, this only addresses the gradient inexactness for the negative
curvature step. For the capped conjugate gradient step as presented in
[Yao+22], they also stipulated:

$\delta_{g,k}\leq\max\\{\epsilon_{g},\|g_{k}\|\\},$

the exact same as the requirement for gradient descent in our Algorithm 1.

Combining requirements for both negative curvature and gradient steps, they
deduced:

$\delta_{g,k}\leq\min\left\\{O\left(\frac{\epsilon_{H}^{2}}{M}\right),\max\\{\epsilon_{g},\|g_{k}\|\\}\right\\}.$
(18)

In contrast, our method requires the same gradient inexactness for both types
of steps: $\Delta_{g,k}\leq\max\\{\epsilon_{g},\|g_{k}\|\\}$.

When we note that [Yao+22] requires strong coupling between tolerances
$\epsilon_{g}$ and $\epsilon_{H}$: For the term
$O\left(\frac{\epsilon_{H}^{2}}{M}\right)$ not to dominate (18), they require
$\epsilon_{H}\geq\sqrt{\epsilon_{g}M}$. In practice, similar to other works
discussed in Section 1.1, [Yao+22] sets $\epsilon_{H}=\sqrt{\epsilon_{g}M}$.
As a consequence, (18) in [Yao+22] becomes
$\delta_{g,k}\leq\min\left\\{O\left(\frac{\epsilon_{H}^{2}}{M}\right),\max\\{\epsilon_{g},\|g_{k}\|\\}\right\\}=O(\epsilon_{g})$,
no longer tolerating _relative_ gradient inexactness.

There exist configurations of $\epsilon_{g}$ and $\epsilon_{H}$ where
$\epsilon_{H}\geq\sqrt{\epsilon_{g}M}$ does not hold. For instance, consider
$\frac{\epsilon_{g}^{2}}{6L}=\frac{2\epsilon_{H}^{3}}{9M^{2}}$. Ignoring
dependencies on $M,L$, and constants, we find
$\epsilon_{H}=O(\epsilon_{g}^{2/3})$, less than $O(\epsilon_{g}^{1/2})$. This
yields a stricter gradient inexactness condition as per Equation (18). In
contrast, our analysis remains resilient to varying couplings between
$\epsilon_{H}$ and $\epsilon_{g}$, consistently requiring gradient inexactness
to satisfy $\Delta_{g,k}=\max\\{\epsilon_{g},\|g_{k}\|\\}$.

A pertinent example arises in matrix recovery contexts, where
$\epsilon_{H}^{3}$ and $\epsilon_{g}^{2}$ are on the same order. The
aforementioned discussions thus gain significance because
$\epsilon_{H}=o(\sqrt{\epsilon_{g}M})$. We illustrate that the strong coupling
of $\epsilon_{H}$ and $\epsilon_{g}$ in [Yao+22], encapsulated in Equation
(18), makes obtaining approximate gradients a harder problem compared with our
Assumption 2.3.

We demonstrate an example that arises naturally in matrix recovery settings in
which $\epsilon_{H}^{3}$ and $\epsilon_{g}^{2}$ are on the same order, and the
above discussion therefore becomes meaningful because
$\epsilon_{H}=o(\sqrt{\epsilon_{g}M})$. We show that the strong coupling
between $\epsilon_{H}$ and $\epsilon_{g}$ required in [Yao+22] (summarized in
Equation (18)) makes obtaining approximate gradients a harder problem,
compared with our Assumption 2.3.

#### Matrix Factorization.

Let $M^{*}\in\mathbb{R}^{d\times d}$ be a real symmetric and positive
semidefinite matrix. Suppose the rank of $M^{*}$ is $r$, where $r<d$, and let
$U\in\mathbb{R}^{r\times d}$. The nonconvex objective function we consider is
$f(U)=\frac{1}{2}\left\|{UU^{\top}-M^{*}}\right\|_{F}^{2}$. The global minimum
$\min_{U\in\mathbb{R}^{d\times r}}f(U)$ solves the symmetric low-rank matrix
factorization problem. Let $\sigma^{\star}_{1}$ and ${\sigma^{\star}_{r}}$
denote the largest singular value and smallest nonzero singular value of
$M^{*}$ respectively. Let $\kappa=\sigma^{\star}_{1}/{\sigma^{\star}_{r}}$ be
the condition number.

The following strict saddle property suggests that approximate second-order
stationary points are practically useful.

######  Definition 3.1.

Function $f(\cdot)$ is $(\epsilon_{g},\epsilon_{H},\zeta)$-_strict saddle_ if,
for any $U$, at least one of following holds:

  * •

$\left\|{\nabla f(U)}\right\|\geq\epsilon_{g}$.

  * •

$\lambda_{\min}(\nabla^{2}f(U))\leq-\epsilon_{H}$.

  * •

$U$ is $\zeta$-close to $\mathcal{X}^{\star}$ — the set of local minima.

The strict saddle point property provides a guarantee that a point is
$\zeta$-close to the region of local optima if it’s an
$(\epsilon_{g},\epsilon_{H})$-approximate second-order stationary point. This
closeness is pivotal since, for example, a good initialization near local
minima is often essential for the convergence of a local algorithm.

The matrix recovery literature, for instance [Jin+17], provides the constants
for the strict saddle property of $f(\cdot)$ and its gradient Lipschitz and
Hessian Lipschitz constants.

######  Fact 3.2.

For $f(\cdot)$, all local minima are global minima. It satisfies
$(\frac{1}{24}({\sigma^{\star}_{r}})^{3/2},\frac{1}{3}{\sigma^{\star}_{r}},\frac{1}{3}({\sigma^{\star}_{r}})^{1/2})$-strict
saddle property. For any $\Gamma>\sigma^{\star}_{1}$, inside the region
$\\{U:\left\|{U}\right\|^{2}<\Gamma\\}$, $f(\cdot)$ has gradient Lipschitz
constant $L=16\Gamma$ and Hessian Lipschitz constant
$M=24\Gamma^{\frac{1}{2}}$.

The function $f$ is known to satisfy a regularity condition similar to strong
convexity inside $\frac{1}{3}{\sigma^{\star}_{r}}^{1/2}$-neighborhood of
$\mathcal{X}^{\star}$. Thus, simple gradient descent algorithms initialized at
an
$(\frac{1}{24}({\sigma^{\star}_{r}})^{3/2},\frac{1}{3}{\sigma^{\star}_{r}})$-approximate
SOSP can achieve local linear convergence.

While the analysis of such algorithms is beyond the scope of this paper, we
pinpoint and compare the target error $\epsilon_{g}$ and $\epsilon_{H}$
required in our work and in [Yao+22] so that an
$(\epsilon_{g},\epsilon_{H})$-approximate second-order stationary point is in
$\frac{1}{3}{\sigma^{\star}_{r}}^{1/2}$-neighborhood of $\mathcal{X}^{\star}$.
For Algorithm 1, there’s no coupling requirement between $\epsilon_{g}$ and
$\epsilon_{H}$, allowing us to set

$\epsilon_{g}=\frac{1}{24}{\sigma^{\star}_{r}}^{3/2},\quad\epsilon_{H}=\frac{1}{3}{\sigma^{\star}_{r}}.$
(19)

In contrast, in [Yao+22], the condition $\epsilon_{H}=\sqrt{\epsilon_{g}M}$
mandates:

$\epsilon_{g}\leq\frac{1}{24}{\sigma^{\star}_{r}}^{3/2},\quad\sqrt{\epsilon_{g}M}\leq\frac{1}{3}{\sigma^{\star}_{r}},$

leading to
$\epsilon_{g}=O\left({\frac{{\sigma^{\star}_{r}}^{2}}{\Gamma^{1/2}}}\right)=O\left({{\kappa^{-1/2}}}\right){{\sigma^{\star}_{r}}^{3/2}}$,
which is $\sqrt{\kappa}$-factor smaller than (19). Given the inexactness
assumption in [Yao+22] which demands
$\delta_{g,k}\leq\min\left\\{O\left(\epsilon_{g}\right),\max\\{\epsilon_{g},\|g_{k}\|\\}\right\\}=O(\epsilon_{g})$,
this $\sqrt{\kappa}$-factor directly makes their gradient inexactness
condition more challenging to satisfy, e.g., requiring more samples in a
stochastic or adversarial setting.

##  4 Sampling

We consider the finite sum problem (2) and investigate the number of samples
that is necessary to construct inexact gradient and Hessian oracles that
satisfy Assumption 2.3. We primarily emphasize the sample complexity needed
for our gradient inexactness condition, as it’s a notable advance from
previous works. Meanwhile, the Hessian inexactness condition resembles
previous studies, and works like [XRM20, Lemma 16] and [Tri+17, Lemma 4]
already cover the required sample complexity for Hessian inexactness.

To derive the gradient sub-sampling bound, we leverage Lemma 11 from [DKM06]:

######  Fact 4.1.

Suppose $A\in\mathbb{R}^{m\times N},B\in\mathbb{R}^{N\times
p},c\in\mathbb{Z}^{+}$ such that $1\leq c\leq N$. Construct
$C\in\mathbb{R}^{m\times c}$ and $R\in\mathbb{R}^{c\times p}$ by sampling
uniformly at random with replacement $c$ columns from $A$ and $c$
corresponding rows from $B$ respectively and rescaling by
$\sqrt{\frac{N}{c}}$. Let $\delta\in(0,1)$. Then with probability at least
$1-\delta$,

$\|AB-
CR\|_{F}\leq\sqrt{\frac{N}{c}}\left(\sum_{k=1}^{N}\left\|A_{:,k}\right\|_{2}^{2}\left\|B_{k,:}\right\|_{2}^{2}\right)^{1/2}+\frac{N}{\sqrt{c}}\sqrt{8\log(1/\delta)}\max_{\alpha}\left\|A_{:,\alpha}\right\|_{2}\left\|B_{\alpha,:}\right\|_{2}.$

We are ready to prove the following sample complexity result for gradient sub-
sampling:

######  Lemma 4.2 (Gradient sub-sampling).

For a given $x\in\mathbb{R}^{d}$, suppose there is an upper bound $G(x)$ such
that $\|\nabla f_{i}(x)\|_{2}\leq G(x)<\infty$ for all $i\in[N]$. For any
given $\xi\in(0,1)$, if $|S_{g}|\geq
16\left(1+\sqrt{8\log(1/\xi)}\right)^{2}\left(\frac{G(x)}{\max\\{\epsilon_{g},\left\|{\nabla
f(x)}\right\|\\}}\right)^{2}$ where $S_{g}$ is uniform with-replacement sub-
sampling indices, then for $g(x):=\frac{1}{|S_{g}|}\sum_{i\in S_{g}}\nabla
f_{i}(x)$, we have

$\mathbb{P}\left(\|\nabla
f(x)-g(x)\|_{2}\leq\frac{1}{3}\max\\{\epsilon_{g},\left\|{g(x)}\right\|\\}\right)\geq
1-\xi.$

###### Proof.

The proof technique was originally from [RM19, Lemma 3]. By the triangle
inequality, we have

$\|\nabla f(x)-g(x)\|_{2}\leq\frac{1}{4}\max\\{\epsilon_{g},\left\|{\nabla
f(x)}\right\|\\},$ (20)

which implies $\|\nabla
f(x)-g(x)\|_{2}\leq\frac{1}{3}\max\\{\epsilon_{g},\left\|{g(x)}\right\|\\}$.
It suffices to prove that (20) holds with high probability.

Write

$A:=[\nabla f_{1}(x)\;\nabla f_{2}(x)\cdots\nabla
f_{N}(x)]\in\mathbb{R}^{d\times N},\qquad
B:=\left[\frac{1}{N}\cdots\frac{1}{N}\right]^{T}\in\mathbb{R}^{N\times 1}.$

Then $\nabla f(x)=AB$ and $g(x)=CR$ using notations in Fact 4.1. With
probability at least $1-\xi$,

$\displaystyle\|\nabla f(x)-g(x)\|_{2}$ $\displaystyle=\|AB-CR\|_{F}$
$\displaystyle\leq\sqrt{\frac{N}{c}}\left(\sum_{k=1}^{N}\left\|A_{:,k}\right\|_{2}^{2}\left\|B_{k,:}\right\|_{2}^{2}\right)^{1/2}+\frac{N}{\sqrt{c}}\sqrt{8\log(1/\xi)}\max_{\alpha}\left\|A_{:,\alpha}\right\|_{2}\left\|B_{\alpha,:}\right\|_{2}$
$\displaystyle=\sqrt{\frac{1}{Nc}}\left(\sum_{k=1}^{N}\left\|\nabla
f_{k}(x)\right\|_{2}^{2}\right)^{1/2}+\frac{1}{\sqrt{c}}\sqrt{8\log(1/\xi)}\max_{\alpha}\left\|\nabla
f_{\alpha}(x)\right\|_{2}$
$\displaystyle\leq\sqrt{\frac{1}{c}}G(x)+\frac{1}{\sqrt{c}}\sqrt{8\log(1/\xi)}G(x)=\frac{G(x)}{\sqrt{c}}\left(1+\sqrt{8\log(1/\xi)}\right)$

Taking
$c=16\left(1+\sqrt{8\log(1/\xi)}\right)^{2}\left(\frac{G(x)}{\max\\{\epsilon_{g},\left\|{\nabla
f(x)}\right\|\\}}\right)^{2}$ gives the desired result. ∎

######  Remark 4.3.

In the context of Remark 2.15, we need $\xi=O(\delta^{2}\epsilon_{g}^{2})$.
This additional complexity only logarithmically impacts the sample size,
making it a mild concern.

######  Remark 4.4 (Discussions on different assumptions).

Using Fact 4.1, we obtained a sample complexity bound
$\tilde{O}\left({\left(\frac{G(x)}{\max\\{\epsilon_{g},\left\|{\nabla
f(x)}\right\|\\}}\right)^{2}}\right)$, assuming only an upper bound $G(x)$ of
the norm of sample gradients. If in addition we assume an upper bound
$V_{g}(x)$ for the variance of the sample gradients

$\mathbb{E}_{i\in[n]}\left\|{\nabla f_{i}(x)-\nabla f(x)}\right\|^{2}\leq
V_{g}(x),$

we can refine this bound using the matrix Berstein inequality to
$\tilde{O}\left({\max\left\\{\frac{V_{g}(x)}{\max\\{\epsilon_{g}^{2},\left\|{\nabla
f(x)}\right\|^{2}\\}},\frac{G(x)}{\max\\{\epsilon_{g},\left\|{\nabla
f(x)}\right\|\\}}\right\\}}\right)$ (similar to [Tri+17, Lemma 4]). If we only
assume the second moment upper bound $V_{g}(x)$ without $G(x)$, then we need
to employ Chebyshev’s inequality to obtain
$O\left({\frac{V_{g}(x)}{\xi^{2}\max\\{\epsilon_{g}^{2},\left\|{\nabla
f(x)}\right\|^{2}\\}}}\right)$, which has a polynomial dependence on $1/\xi$.

######  Remark 4.5.

With the same analysis method, [Yao+22, Lemma 2.13] derived the following
gradient subsampling bound

$|S_{g}|\geq
16\left(1+\sqrt{8\log(1/\xi)}\right)^{2}\left(\frac{G(x)}{\epsilon_{g}}\right)^{2}.$

Lemma 4.2 improves on this bound because we allow for relative gradient
inexactness.

We can also utilize sub-sampling to estimate Hessians. Our result does not
improve over [Yao+22] since our Hessian inexactness requirement in Assumption
2.3 is the same as theirs, but we still state the result for completeness.

######  Lemma 4.6 (Hessian sub-sampling [XRM20, Lemma 16]).

For a given $x\in\mathbb{R}^{d}$, suppose there is an upper bound $K(x)$ such
that $\|\nabla^{2}f_{i}(x)\|\leq K(x)<\infty$ for all $i\in[N]$. For any given
$\xi\in(0,1)$, if $|S_{H}|\geq
484\log\frac{2d}{\xi}\left(\frac{K(x)}{\epsilon_{H}}\right)^{2}$, where
$S_{H}$ is sub-sampling indices (with- or without-replacement), then for
$H(x):=\frac{1}{|S_{H}|}\sum_{i\in S_{H}}\nabla^{2}f_{i}(x)$, we have

$\mathbb{P}\left(\|\nabla^{2}f(x)-H(x)\|_{2}\leq\frac{2}{9}\max\\{\epsilon_{g},\left\|{g(x)}\right\|\\}\right)\geq
1-\xi.$

##  5 Conclusion

In this paper, we presented a simple randomized algorithm that finds an
$(\epsilon_{g},\epsilon_{H})$-approximate second-order stationary point of a
smooth nonconvex function with inexact oracle access to gradient and Hessian.
Compared to previous studies, our method stands out due to its enhanced
tolerance for inexactness. Specifically, the gradient inexactness measure can
be relative, and no coupling between $\epsilon_{g}$ and $\epsilon_{H}$ is
required. Both features make our inexactness criteria easier to satisfy, as
evidenced in the low-rank matrix factorization example and when we examine the
subsampling complexities of finite-sum problems.

For future directions, it is worth noting that our runtime analysis hinges on
the idea that first-order steps give sufficient descent in the function
values, with gradient descent being a typical example. Thus, the Radamacher
randomness (Lines 1 of Algorithm 1) could also be injected to the saddle point
escaping steps in other algorithms, like those in [Yao+22, Car+18, RW18,
ROW20]. These algorithms might offer better runtime performance because of
their refined descent methods but could face challenges in meeting inexactness
criteria due to their negative curvature analysis, which our paper greatly
improves.

## References

  * [Aga+17] Naman Agarwal, Zeyuan Allen-Zhu, Brian Bullins, Elad Hazan and Tengyu Ma  “Finding approximate local minima faster than gradient descent”  In _STOC’17—Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing_ ACM, New York, 2017, pp. 1195–1199  DOI: 10.1145/3055399.3055464
  * [All18] Zeyuan Allen-Zhu  “Natasha 2: Faster Non-Convex Optimization Than SGD”  In _Advances in Neural Information Processing Systems_ 31 Curran Associates, Inc., 2018
  * [Arj+20] Yossi Arjevani, Yair Carmon, John C. Duchi, Dylan J. Foster, Ayush Sekhari and Karthik Sridharan  “Second-Order Information in Non-Convex Stochastic Optimization: Power and Limitations”  In _Proceedings of Thirty Third Conference on Learning Theory_ 125, Proceedings of Machine Learning Research  PMLR, 2020, pp. 242–299  URL: https://proceedings.mlr.press/v125/arjevani20a.html
  * [BCS21] Albert S Berahas, Liyuan Cao and Katya Scheinberg  “Global convergence rate analysis of a generic line search algorithm with noise”  In _SIAM Journal on Optimization_ 31.2 SIAM, 2021, pp. 1489–1518
  * [Ber+21] El-houcine Bergou, Youssef Diouane, Vladimir Kunc, Vyacheslav Kungurtsev and Clément W. Royer  “A Subsampling Line-Search Method with Second-Order Results”, 2021  arXiv:1810.07211 [math.OC]
  * [Bla+19] Jose Blanchet, Coralia Cartis, Matt Menickelly and Katya Scheinberg  “Convergence rate analysis of a stochastic trust-region method via supermartingales”  In _INFORMS journal on optimization_ 1.2 INFORMS, 2019, pp. 92–119
  * [Car+18] Y. Carmon, J.. Duchi, O. Hinder and A. Sidford  “Accelerated methods for non-convex optimization”  In _siopt_ 28, 2018, pp. 1751–1772
  * [CGT11] Coralia Cartis, Nicholas I.. Gould and Philippe L. Toint  “Adaptive Cubic Regularisation Methods for Unconstrained Optimization. Part II: Worst-Case Function- and Derivative-Evaluation Complexity”  In _Math. Program._ 130.2 Berlin, Heidelberg: Springer-Verlag, 2011, pp. 295–319  DOI: 10.1007/s10107-009-0337-y
  * [CGT12] C. Cartis, N.I.M. Gould and Ph.L. Toint  “Complexity bounds for second-order optimality in unconstrained optimization”  In _Journal of Complexity_ 28, 2012, pp. 93–108  DOI: 10.1016/j.jco.2011.06.001
  * [CMS18] Ruobing Chen, Matt Menickelly and Katya Scheinberg  “Stochastic optimization using a trust-region method and random models”  In _Mathematical Programming_ 169 Springer, 2018, pp. 447–487
  * [CS18] Coralia Cartis and Katya Scheinberg  “Global convergence rate analysis of unconstrained optimization methods based on probabilistic models”  In _Mathematical Programming_ 169 Springer, 2018, pp. 337–375
  * [CSS19] Frank E Curtis, Katya Scheinberg and Rui Shi  “A stochastic trust region algorithm based on careful step normalization”  In _Informs Journal on Optimization_ 1.3 INFORMS, 2019, pp. 200–220
  * [Cur+21] F.. Curtis, D.. Robinson, C.. Royer and S.. Wright  “Trust-region Newton-CG with strong second-order complexity guarantees for nonconvex optimization”  In _SIAM Journal on Optimization_ 31, 2021, pp. 518–544
  * [DKM06] Petros Drineas, Ravi Kannan and Michael W. Mahoney  “Fast Monte Carlo algorithms for matrices. I. Approximating matrix multiplication”  In _SIAM J. Comput._ 36.1, 2006, pp. 132–157  DOI: 10.1137/S0097539704442684
  * [Dur19] Rick Durrett  “Probability—theory and examples” 49, Cambridge Series in Statistical and Probabilistic Mathematics  Cambridge University Press, Cambridge, 2019, pp. xii+419  DOI: 10.1017/9781108591034
  * [Ge+15] Rong Ge, Furong Huang, Chi Jin and Yang Yuan  “Escaping From Saddle Points — Online Stochastic Gradient for Tensor Decomposition”  In _Proceedings of The 28th Conference on Learning Theory_ 40, Proceedings of Machine Learning Research  Paris, France: PMLR, 2015, pp. 797–842  URL: https://proceedings.mlr.press/v40/Ge15.html
  * [GJT23] Serge Gratton, Sadok Jerad and Philippe L. Toint  “Convergence Properties of an Objective-Function-Free Optimization Regularization Algorithm, Including an $\mathcal{O}(\epsilon^{-3/2})$ Complexity Bound”  In _SIAM Journal on Optimization_ 33.3, 2023, pp. 1621–1646  DOI: 10.1137/22M1499522
  * [GJZ17] Rong Ge, Chi Jin and Yi Zheng  “No Spurious Local Minima in Nonconvex Low Rank Problems: A Unified Geometric Analysis”  In _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ICML’17  Sydney, NSW, Australia: JMLR.org, 2017, pp. 1233–1242
  * [GYY16] Geovani Nunes Grapiglia, J Yuan and Y Yuan  “On the worst-case complexity of nonlinear stepsize control algorithms for convex unconstrained optimization”  In _Optimization Methods and Software_ 31.3 Taylor & Francis, 2016, pp. 591–604
  * [Jin+17] Chi Jin, Rong Ge, Praneeth Netrapalli, Sham M. Kakade and Michael I. Jordan  “How to Escape Saddle Points Efficiently”  In _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ICML’17  Sydney, NSW, Australia: JMLR.org, 2017, pp. 1724–1732
  * [Jin+21] Chi Jin, Praneeth Netrapalli, Rong Ge, Sham M. Kakade and Michael I. Jordan  “On Nonconvex Optimization for Machine Learning: Gradients, Stochasticity, and Saddle Points”  In _J. ACM_ 68.2 New York, NY, USA: Association for Computing Machinery, 2021  DOI: 10.1145/3418526
  * [JSX23] Billy Jin, Katya Scheinberg and Miaolan Xie  “Sample Complexity Analysis for Adaptive Optimization Algorithms with Stochastic Oracles”, 2023  arXiv:2303.06838 [math.OC]
  * [NP06] Y. Nesterov and B.. Polyak  “Cubic regularization of Newton method and its global performance”  In _Mathematical Programming, Series A_ 108, 2006, pp. 177–205
  * [PS20] Courtney Paquette and Katya Scheinberg  “A stochastic line search method with expected complexity analysis”  In _SIAM J. Optim._ 30.1, 2020, pp. 349–376  DOI: 10.1137/18M1216250
  * [RM19] Farbod Roosta-Khorasani and Michael W. Mahoney  “Sub-sampled Newton methods”  In _Math. Program._ 174.1-2, Ser. B, 2019, pp. 293–326  DOI: 10.1007/s10107-018-1346-5
  * [ROW20] Clément W. Royer, Michael O’Neill and Stephen J. Wright  “A Newton-CG algorithm with complexity guarantees for smooth unconstrained optimization”  In _Math. Program._ 180.1-2, Ser. A, 2020, pp. 451–488  DOI: 10.1007/s10107-019-01362-7
  * [RW18] C.. Royer and S.. Wright  “Complexity analysis of second-order line-search algorithms for smooth nonconvex optimization”  In _SIAM Journal on Optimization_ 28, 2018, pp. 1448–1477
  * [SQW16] Ju Sun, Qing Qu and John Wright  “A geometric analysis of phase retrieval”  In _2016 IEEE International Symposium on Information Theory (ISIT)_ , 2016, pp. 2379–2383  DOI: 10.1109/ISIT.2016.7541725
  * [SQW17] Ju Sun, Qing Qu and John Wright  “Complete Dictionary Recovery Over the Sphere I: Overview and the Geometric Picture”  In _IEEE Trans. Inf. Theor._ 63.2 IEEE Press, 2017, pp. 853–884  DOI: 10.1109/TIT.2016.2632162
  * [Tri+17] Nilesh Tripuraneni, Mitchell Stern, Chi Jin, Jeffrey Regier and Michael I. Jordan  “Stochastic Cubic Regularization for Fast Nonconvex Optimization”  In _CoRR_ abs/1711.02838, 2017  arXiv: http://arxiv.org/abs/1711.02838
  * [WR22] Stephen J. Wright and Benjamin Recht  “Optimization for data analysis”  Cambridge University Press, Cambridge, 2022, pp. x+227  DOI: 10.1017/9781009004282
  * [XRM20] Peng Xu, Fred Roosta and Michael W. Mahoney  “Newton-Type Methods for Non-Convex Optimization under Inexact Hessian Information”  In _Math. Program._ 184.1-2 Berlin, Heidelberg: Springer-Verlag, 2020, pp. 35–70  DOI: 10.1007/s10107-019-01405-z
  * [Yao+22] Zhewei Yao, Peng Xu, Fred Roosta, Stephen J Wright and Michael W Mahoney  “Inexact Newton-CG algorithms with complexity guarantees”  In _IMA Journal of Numerical Analysis_ , 2022  DOI: 10.1093/imanum/drac043
  * [ZQW20] Yuqian Zhang, Qing Qu and John Wright  “From Symmetry to Geometry: Tractable Nonconvex Problems”  In _CoRR_ abs/2007.06753, 2020  arXiv: https://arxiv.org/abs/2007.06753