text:
\begin{document}

\maketitle
 \begin{abstract}
  We design an algorithm which finds an $\epsilon$-approximate stationary
  point (with $\|\nabla F(x)\|\le \epsilon$) using $O(\epsilon^{-3})$
  stochastic gradient and Hessian-vector products, matching guarantees
  that were previously available only under a stronger assumption of access
  to multiple queries with the same random seed.
 We prove a lower bound which establishes that this rate
is optimal and---surprisingly---that it cannot be improved using stochastic $p$th order methods
for any $p\ge 2$, even when the first $p$ derivatives of the objective
are Lipschitz. Together, these results characterize the complexity of
non-convex stochastic optimization with second-order methods and beyond.
Expanding our scope to the oracle complexity of finding
$(\epsilon,\gamma)$-approximate second-order stationary points, we
establish nearly matching upper and lower bounds for stochastic
second-order methods. Our lower bounds here are novel even in the noiseless
case.
\end{abstract}

\colt{
\begin{keywords}
  Stochastic optimization, non-convex optimization, second-order
  methods, variance reduction, Hessian-vector products.
\end{keywords}
}

\section{Introduction}\label{sec:intro}

Let $F:\R^d\to\R$ have Lipschitz continuous gradient and Hessian, and
consider the task of finding an $(\epsilon, \gamma)$-second-order
stationary point (SOSP), that is, $x\in\R^d$ such that
\begin{equation}\label{eq:sosp}
\norm{\grad F(x)} \le \eps
\quad\mbox{and}\quad
\hess F(x) \succeq - \gamma I.
\end{equation}
This task plays a central role in the study of non-convex optimization: for
functions satisfying a weak strict saddle condition \citep{ge2015escaping}, exact SOSPs
(with
$\epsilon=\gamma=0$) are
local minima, and therefore the condition \eqref{eq:sosp} serves
as a proxy for
approximate local optimality.\footnote{
	 However, it is NP-Hard to decide whether a SOSP
	 is a local minimum or a high-order saddle
	 point~\citep{murty1987some}.}
Moreover, for a growing set of non-convex optimization problems arising
in
machine learning, SOSPs are in fact \emph{global minima} \citep{ge2015escaping, ge2016matrix, sun2018geometric,
mc2019implicit}.
Consequently, there has been intense recent interest in the design of efficient
algorithms
for finding approximate SOSPs~\citep{jin2017how,allen2018how,carmon2018accelerated,fang2018spider,tripuraneni2018stochastic,xu2018first,fang2019sharp}.

In stochastic approximation tasks---particularly those motivated by
machine learning---access
to the objective function is often restricted to
stochastic estimates of its gradient; for each query point $x\in\R^d$ we observe
$\stocgrad{x,z}$, where $z\sim P_z$ is a random variable such that
\begin{equation}\label{eq:stoc-grad}
\E \bigl[\stocgrad{x,z}\bigr] = \grad F(x)
~~\mbox{and}~~
\E \, \norm{\stocgrad{x,z}  - \grad F(x)}^2 \le \sigma_1^2.
\end{equation}
This restriction typically arises due to computational considerations (when
$\stocgrad{\cdot{}, z}$ is much cheaper to compute than $\grad F(\cdot{})$, as in
empirical risk minimization or Monte Carlo simulation), or due to
fundamental online nature of the
 problem at hand (e.g., when $x$ represents a routing scheme
and $z$ represents traffic on a given day). However, for many problems with additional structure, we have
access to extra information. For example, we
often have access to stochastic second-order
information in the form of a Hessian estimator $\stochess{x,z}$ satisfying
\begin{equation}\label{eq:stoc-hess}
\E \bigl[\stochess{x,z}\bigr] = \hess F(x)
~~\mbox{and}~~
\E \, \opnorm{\stochess{x,z}  - \hess F(x)}^2 \le \sigma_2^2.
\end{equation}

In this paper, we characterize the extent to which the stochastic Hessian
information~\pref{eq:stoc-hess}, as well as higher-order information, contributes to the efficiency of
finding first- and second-order stationary points. We approach
this question from the
perspective of \emph{oracle complexity}~\citep{nemirovski1983problem},
which measures
efficiency by the number of queries to estimators of the
form~\pref{eq:stoc-grad}---and
possibly~\pref{eq:stoc-hess}---required to satisfy the
condition~\pref{eq:sosp}.

\subsection{Our Contributions}
\arxiv{\begin{figure}[t]
      \centering
      		\includegraphics[width=0.6\textwidth]{{figures/full_complexity.pdf}}
      \vspace{.85cm}
      		\caption{The \emph{elbow effect}: For stochastic
      			oracles, the optimal complexity sharply improves from $\veps^{-4}$
      			for
      			$p=1$ to $\veps^{-3}$ for $p=2$,
      			but there is no further improvement for
                        $p>2$. For noiseless oracles, the optimal
                        complexity begins at $\eps^{-2}$ for $p=1$ and smoothly
      			approaches
      			$\eps^{-1}$ as the derivative order
      			$p\to\infty$.}\label{fig:elbow}
\end{figure}}

\colt{
\captionsetup[figure]{font=scriptsize}
 \begin{figure}[t]
      \centering
      \begin{floatrow}
      	\ffigbox[0.43\textwidth]{
      		      		\centering
      		\resizebox{0.45\textwidth}{!}{
      		\includegraphics[width=0.5\textwidth]{{figures/full_complexity.pdf}}
              }
      \vspace{.85cm}
    }{

      		\caption{The ``elbow effect:'' for stochastic
      			oracles, optimal complexity sharply improves from $\veps^{-4}$
      			for
      			$p=1$ to $\veps^{-3}$ for $p=2$,
      			but has no further improvement for $p>2$. Noiseless
      			complexity begins at $\eps^{-2}$ for $p=1$ and smoothly
      			approaches
      			$\eps^{-1}$ as the derivative order
      			$p\to\infty$.}\label{fig:elbow}
      	}
        ~
      	\capbtabbox[0.5\textwidth]{
      		\resizebox{0.5\textwidth}{!}{
      		\begin{tabular}{p{3.8cm} p{1.4cm} p{1.6cm} p{2.5cm}}
      			\toprule
      			\makecell[lt]{Method} & \makecell[lt]{Requires \\
      				$\stochessNA$?} & \makecell[lt]{Complexity\\ bound} &
      			\makecell[lt]{Additional\\ assumptions} \\
      			\midrule
      			\makecell[lt]{SGD~\citep{ghadimi2013stochastic}\\~} & No
      			&
      			$O(\epsilon^{-4})$ &\\

      			\makecell[lt]{Restarted SGD~\citep{fang2019sharp}} & No &
      			$O(\eps^{-3.5})$ & \makecell[tl]{$\stocgradNA$ Lipschitz\\
      				almost surely}

      			\\
      			\makecell[lt]{Subsampled regularized\\
      				Newton~\citep{tripuraneni2018stochastic}} & Yes &
      				$O(\eps^{-3.5})$ & \\
      			\makecell[lt]{Recursive variance\\
      			reduction~\citep[e.g.,][]{fang2018spider}} & No &
      			$O(\epsilon^{-3})$
      			&\makecell[tl]{Mean-squared
      				\\smoothness,\\Sim.\ queries \\
      				(see \pref{app:oracle})}

      			\\
      			\Xhline{4\arrayrulewidth}
      			\addlinespace[3pt]
      			Hessian-vector recursive VR ({\pref{alg:epsilon_hvp}}) & Yes &
      			${O(\eps^{-3})}$& None\\
      			Subsampled Newton with VR ({\pref{alg:epsilon_cubic}}) & Yes &
      			${O(	\eps^{-3})}$& None\\
      			\bottomrule
      		\end{tabular}}
      	}{
      		 \caption{Comparison of guarantees for finding $\eps$-stationary
      			points
      			(i.e., $\E \norm{\grad F(x)}\le \eps$) for a function $F$ with
      			Lipschitz gradient and Hessian. See
      			\pref{table:summary_detailed} for detailed comparison.}
      		\label{table:summary}
      	}
    \end{floatrow}
\end{figure}
\captionsetup[figure]{font=normalsize}

}

We provide new upper and lower bounds on the
stochastic oracle complexity of finding $\eps$-stationary points and ($\epsilon,\gamma)$-SOSPs. In
brief, our main results are as follows.
\arxiv{
  \begin{itemize}
  }
  \colt{
  \begin{itemize}[leftmargin=\parindent]
    }
\item \textbf{Finding $\veps$-stationary points: The elbow effect.}
	We propose a new algorithm that finds an $\veps$-stationary point
	($\gamma=\infty$) with  $O(\veps^{-3})$ stochastic
	gradients and stochastic Hessian-vector products. We furthermore show
	that this guarantee is not improvable via a complementary
  $\Omega(\eps^{-3})$ lower bound. All previous algorithms achieving
  $O(\veps^{-3})$ complexity require ``multi-point'' queries, in
  which the algorithm can query stochastic gradients at multiple
  points for the same random seed. Moreover, we show that
  $\Omega(\eps^{-3})$ remains a lower bound for stochastic \pth-order
  methods for all $p\geq{}2$ and hence---in contrast to the deterministic
  setting---the optimal rates for higher-order methods exhibit an ``elbow
  effect''; see \pref{fig:elbow}.

\item \textbf{$(\veps,\gamma)$-stationary points: Improved algorithm and
  nearly matching lower bound.} We extend our algorithm to find $(\eps,\gamma)$-stationary points using $O(\veps^{-3}+
  \veps^{-2}\gamma^{-2} + \gamma^{-5})$ stochastic gradient and
  Hessian-vector products, and prove a nearly matching   $\Omega(\veps^{-3}+ \gamma^{-5})$ lower bound.
\end{itemize}

In the remainder of this section we overview our results in greater
detail. Unless otherwise stated, we assume $F$ has both Lipschitz
gradient and Hessian. To simplify the overview, we focus on dependence
on $\eps^{-1}$ and $\gamma^{-1}$ while keeping the other
parameters---namely the initial optimality
gap $F(x\ind{0})-\inf_{x\in\R^d} F(x)$, the Lipschitz constants of $\grad F$
and $\hess
F$,  and the variances of their estimators---held fixed. Our main theorems
give
explicit dependence on these parameters.

\subsubsection{First-order stationary points ($\gamma=\infty$)}
We first describe our developments for the task of finding
$\eps$-approximate first-order
stationary points (satisfying~\pref{eq:sosp} with $\gamma=\infty$), and
subsequently extend our results to general $\gamma$. The reader may also
refer to \pref{table:summary}  for a succinct comparison of upper bounds.

\paragraph{Variance reduction via Hessian-vector products: A new gradient estimator.}
Using stochastic gradients and stochastic Hessian-vector products as primitives, we
design a new variance-reduced gradient estimator. Plugging it into standard stochastic gradient descent (SGD), we obtain an algorithm that returns a
point
$\hx$ satisfying $\E \,\norm{\grad F(\hx)}\le \eps$
and requires
$O(\eps^{-3})$ stochastic gradient and HVP queries in expectation. In comparison, vanilla SGD requires
$O(\eps^{-4})$ queries~\citep{ghadimi2013stochastic}, and the previously
best known rate under our assumptions
was $O(\eps^{-3.5})$, by both cubic-regularized Newton's method and a
restarted variant of SGD \citep{tripuraneni2018stochastic,fang2019sharp}.

Our approach builds on a line of work by~\citet{fang2018spider,zhou2018stochastic,
	wang2018spiderboost,cutkosky2019momentum}
      that also develop algorithms with complexity $O(\eps^{-3})$, but
      require a  ``multi-point'' oracle
      in which algorithm can query the stochastic gradient at
      multiple points for the same random seed. Specifically, in the
      $n$-point variant of this model, the algorithm can query at the set of points $(x_1, \ldots, x_n)$ and receive
\begin{equation}
			\stocder{x_1, z}{}, \ldots, \stocder{x_n, z}{}, \quad \text{where}
			\quad z \overset{\iid}{\sim} P_z, \label{eq:n_point_query_model}
\end{equation}
and where the estimator $\stocder{x,z}{}$ is unbiased and has bounded
variance in the sense of \pref{eq:stoc-grad}. The aforementioned works achieve
$O(\eps^{-3})$ complexity using $n=2$ simultaneous queries, while our
new algorithm
achieves the same rate using $n=1$ (i.e., $z$ is drawn afresh at each
query), but using stochastic Hessian-vector products in addition to
stochastic gradients. However, we show in
\pref{app:oracle} that under the statistical assumptions made in these
works, the two-point stochastic gradient oracle model is
\emph{strictly} stronger than the single-point stochastic
gradient/Hessian-vector product oracle we consider here. On the other hand, unlike our
algorithm, these works do not require Lipschitz Hessian.

The algorithms that achieve complexity $O(\eps^{-3})$ using
two-point queries
work by estimating gradient differences of the form $\grad
F(x) - \grad F(x')$ using
$\stocgrad{x,z} - \stocgrad{x',z}$ and applying recursive variance
reduction~\citep{nguyen2017sarah}. Our primary algorithmic
contribution is a second-order stochastic estimator
for $\grad F(x) - \grad F(x')$ which avoids simultaneous queries while
maintaining comparable error guarantees. To derive our estimator, we note
that $\grad F(x) - \grad F(x')  = \int_0^1 \hess F(xt + x'(1-t))
(x-x') dt$, and
use $K$ queries to the stochastic
Hessian
estimator~\pref{eq:stoc-hess} to numerically approximate this
integral.\footnote{
	More precisely, our estimator~\pref{eq:key-estimator} only requires
	stochastic Hessian-vector
	products, whose computation is often roughly as expensive as that of a
	stochastic gradient~\citep{pearlmutter1994fast}.}
      Specifically, our estimator takes the form
\begin{equation}\label{eq:key-estimator}
\frac{1}{K}\sum_{k=0}^{K-1} \wh{\grad^{2}F}\bigl(x\cdot (1-\tfrac{k}{K}) + x' \cdot
\tfrac{k}{K},
z\ind{i} \bigr)
(x-x'),
\end{equation}
where $z\ind{i}\iidsim{}P_z$. Unlike the usual estimator $\stocgrad{x,z} - \stocgrad{x',z}$, the
estimator~\pref{eq:key-estimator} is biased.
Nevertheless, we show that
\emph{choosing $K$ dynamically} according to $K\propto \norm{x-x'}^2$
provides adequate control over both bias and variance while maintaining
the desired
query complexity. Combining the
integral estimator~\pref{eq:key-estimator} with recursive variance
reduction, we attain $O(\eps^{-3})$ complexity.

\paragraph{Demonstrating the power of second-order information.}
For functions with Lipschitz gradient and Hessian, we prove an
$\Omega(\eps^{-3.5})$ lower bound on the minimax oracle complexity
of algorithms
for finding stationary points using \emph{only} stochastic
gradients~\pref{eq:stoc-grad}.\footnote{We formally prove our results for the structured class
 of \emph{zero-respecting algorithms}~\citep{carmon2019lower_i}; the
 lower bounds
 extend to general randomized algorithms via similar arguments
 to~\citet
 {arjevani2019lower}.} This lower
bound is an extension of the results
of~\citet{arjevani2019lower}, who showed that for functions with Lipschitz
  gradient but \emph{not} Lipschitz Hessian, the optimal rate is
  $\Theta(\eps^{-4})$ using \textit{only} stochastic
  gradients~\pref{eq:stoc-grad}. Together with our new $O(\eps^{-3})$ upper bound, this lower
bound reveals
that stochastic Hessian-vector products offer an $\Omega(\eps^{-0.5})$
improvement in the oracle complexity for finding stationary points in
the single-point query model.
This contrasts the noiseless optimization setting, where
finite
gradient differences can approximate Hessian-vector products arbitrarily well,
meaning these oracle models are equivalent.
\arxiv{
\begin{table}
	\small
  \renewcommand{\arraystretch}{1.1}
      		\begin{tabular}{p{5.5cm} p{1.6cm} p{2cm} p{5.8cm}}
      			\toprule
      			\makecell[lt]{Method} & \makecell[lt]{Requires \\
      				$\stochessNA$?} & \makecell[lt]{Complexity\\ bound} &
      			\makecell[lt]{Additional\\ assumptions} \\
      			\midrule
      			\makecell[lt]{SGD \citep{ghadimi2013stochastic}} & No
      			&
      			$O(\epsilon^{-4})$ &\\[9pt]
      	      			\makecell[lt]{Restarted SGD \citep{fang2019sharp}} & No &
      			$O(\eps^{-3.5})$ & \makecell[tl]{$\stocgradNA$ Lipschitz\\
      				almost surely}
                  \\[9pt]
      			\makecell[lt]{Subsampled regularized
      				Newton\\\citep{tripuraneni2018stochastic}} & Yes$^{*}$ &
      				$O(\eps^{-3.5})$ & \\[9pt]
      			\makecell[lt]{Recursive variance
      			reduction\\\citep[e.g.,][]{fang2018spider}} & No &
      			$O(\epsilon^{-3})$
      			&\makecell[tl]{Mean-squared
      				smoothness, \\simultaneous queries
      				(see \pref{app:oracle})}
      			\\
      			\Xhline{4\arrayrulewidth}
      			\addlinespace[3pt]
                  \makecell[lt]{SGD with \mainalg\\({\pref{alg:epsilon_hvp}})} & Yes$^{*}$ &
      			${O(\eps^{-3})}$& None\\[9pt]
      			Subsampled Newton \\w/ \mainalg ({\pref{alg:epsilon_cubic}}) &
      			Yes &
      			${O(	\eps^{-3})}$& None\\
      			\bottomrule
      		\end{tabular}
      		 \caption{Comparison of guarantees for finding $\eps$-stationary
      			points
      			(i.e., $\E \norm{\grad F(x)}\le \eps$) for a function $F$ with
      			Lipschitz gradient and Hessian. See
      			\pref{table:summary_detailed} for explicit
                        dependence on problem parameters. Algorithms
                        marked with $^{*}$ require only stochastic Hessian-vector products.}
      		\label{table:summary}
		\end{table}
}

\paragraph{Demonstrating the limitations of higher-order
information ($p > 2$).}
For algorithms that can query both stochastic gradients and stochastic
Hessians, we prove a lower bound of $\Omega(\eps^{-3})$ on the oracle
complexity of finding an expected $\epsilon$-stationary point.
This proves that our $O(\eps^{-3})$ upper bound is optimal in the leading order term in
$\eps$, despite using only stochastic Hessian-vector products rather than
full stochastic Hessian queries.

Notably, our $\Omega(\eps^{-3})$ lower
bound extends to
settings where stochastic higher-order oracles are available, i.e,
when the first $p$ derivatives are Lipschitz and we have bounded-variance estimators $\{\stocder{\cdot,
\cdot}{q}\}_{q \le p}$. The lower bound holds for any finite $p$, and thus, as a function of the
oracle order $p$, the minimax complexity has an elbow
(\pref{fig:elbow}):
for $p=1$ the complexity is $\Theta(\eps^{-4})$~\citep{arjevani2019lower}
while for all $p\ge 2$ it is $\Theta(\eps^{-3})$. This means that
smoothness and stochastic
derivatives beyond the second-order cannot improve the leading term in rates of
convergence to stationarity, establishing a fundamental limitation of
stochastic high-order information. This highlights another contrast with the noiseless setting, where $p$th order methods enjoy
improved complexity for every $p$ \citep{carmon2019lower_i}.

As we discuss in
\pref{app:oracle}, for multi-point stochastic oracles
\pref{eq:n_point_query_model}, the
rate $O(\eps^{-3})$ is attainable even without stochastic Hessian
access. Moreover, our $\Omega(\eps^{-3})$ lower bound for stochastic
$p$th order oracles holds even when multi-point queries are
allowed. Consequently, when viewed through the lens of  worst-case oracle complexity, our lower bounds show that even stochastic Hessian information is not
helpful in the multi-point setting.

\subsubsection{Second-order stationary points}

\paragraph{Upper bounds for general $\gamma$.}
We incorporate our recursive variance-reduced Hessian-vector product-based gradient
estimator into an algorithm that combines SGD with negative curvature
search. Under the slightly stronger (relative to \pref{eq:stoc-hess})
assumption that the stochastic Hessians have almost surely bounded error,
we prove that---with constant probability---the algorithm returns an $(\epsilon, \gamma)$-SOSP after performing $O(\eps^{-3} +
\eps^{-2}\gamma^{-2} + \gamma^{-5})$  stochastic gradient and
Hessian-vector product queries.
\paragraph{A lower bound for finding second-order stationary points.}
We prove a minimax lower bound which establishes that the stochastic second-order oracle complexity of finding
$(\eps, \gamma)$-SOSPs is $\Omega(\eps^{-3} + \gamma^{-5})$.
Consequently, the algorithms we develop have optimal worst-case
complexity in the regimes $\gamma=O(\epsilon^{2/3})$ and $\gamma =
\Omega(\eps^{0.5})$. Compared to
our lower bounds for finding $\veps$-stationary points,
proving the $\Omega(\gamma^{-5})$ lower bound requires a more
substantial modification of the constructions
of~\cite{carmon2019lower_i} and \cite{arjevani2019lower}. In fact, our lower bound is new
even in the noiseless regime (i.e., $\sigma_1=\sigma_2=0$), where it
becomes $\Omega(\eps^{-1.5} + \gamma^{-3})$; this matches the
guarantee of the cubic-regularized Newton's
method~\citep{nesterov2006cubic} and consequently characterizes the
optimal rate for finding approximate SOSPs using noiseless second-order
methods.
\subsection{Further related work} We briefly survey additional upper and lower complexity bounds
related to our work and
place our results within their context.
The works of \citet{monteiro2013accelerated,arjevani2019oracle,agarwal2017lower} delineate the second-order
oracle complexity
of \emph{convex} optimization in the noiseless setting;
\cite{arjevani2017oracle} treat the finite-sum
setting.

For
functions with Lipschitz gradient and Hessian, oracle access to the Hessian
significantly accelerates convergence to $\varepsilon$-approximate global
minima, reducing the complexity from $\Theta(\varepsilon^{-0.5})$ to
$\Theta(\varepsilon^{-2/7})$.
However, since the hard instances for first-order
convex optimization are
quadratic~\citep{nemirovski1983problem,arjevani2016iteration,simchowitz2018randomized},
 assuming Lipschitz continuity of the Hessian does not improve the
complexity if one only has access to a first-order oracle. This contrasts the case for finding $\eps$-approximate stationary points of
\emph{non-convex}
functions with noiseless oracles. There, Lipschitz continuity of the
Hessian
improves the first-order oracle complexity from $\Theta(\eps^{-2})$ to
$O(\eps^{-1.75})$, with a lower bound of $\Omega(\eps^{-12/7})$ for
deterministic algorithms~\citep{carmon2017convex,carmon2019lower_ii}.
Additional access to full Hessian further improves this complexity to
$\Theta(\eps^{-1.5})$, and for \pth-order oracles with Lipschitz \pth
derivative, the complexity further improves to
$\Theta(\eps^{-(1+\frac{1}{p})})$ \citep{carmon2019lower_i}; see
\pref{fig:elbow}.

\subsection{Paper organization}
We formally introduce our notation
and oracle model in \pref{sec:prelim}.
\pref{sec:epsilon} contains our results concerning the complexity of finding
$\epsilon$-first-order stationary points: algorithmic upper bounds
(\pref{sec:epsilon_upper}) and algorithm-independent lower bounds
(\pref{sec:epsilon_lower}). Following a similar outline,
\pref{sec:epsilon_gamma} describes our upper and lower bounds for finding
$(\eps, \gamma)$-SOSPs. \arxiv{We conclude the paper in \pref{sec:conclusion}
with a discussion of directions for further research. }\colt{In
\pref{app:conclusion}, we discussion directions for future research. }Additional
technical comparison with related work is given in \pref{app:comparison}
and \ref{app:oracle}, and proofs are given in \pref{app:estimators}
through \pref{app:lower}.

\newcommand{\xmark}{\ding{55}}
\newcommand{\cmark}{\ding{51}}

\paragraph{Notation.}
We let $\cC^p$ denote the class of $p$-times
differentiable real-valued functions, and let $\grad^q F$
denote the $q$th derivative of a
given function $F\in \cC^p$ for $q\in\{1,\dots,p\}$.
Given a function $F\in\cC^1$, we let $
\grad_i F(x) \defeq \brk*{\grad
F(x)}_i = \frac{\partial}{\partial x_i} F(x)$.  When $F\in\cC^2$
is twice
differentiable, we
define, $\nabla^2_{ij} f(x)
\ldef \brk*{\grad^2 f(x)}_{ij} =
\frac{\partial^2}{\partial x_i \partial x_j} f(x)$,
and similarly
define $\brk*{\grad^p f(x)}_{i_1,i_2,\ldots,i_p}=\frac{\partial^p}{\partial
x_{i_1}\cdots \partial x_{i_p}} f(x)$ for
\pth-order
derivatives.
For a vector $x\in\bbR^{d}$, $\nrm*{x}$ denotes the Euclidean
norm and
$\nrm*{x}_\infty$ denotes the $\ls_\infty$ norm. For matrices
$A\in\bbR^{d\times{}d}$, $\nrm*{A}_{\op}$ denotes the operator
norm. More generally, for symmetric \pth
order tensors $T$, we define the operator norm via
$\nrm*{T}_{\op}=\sup_{\nrm*{v}=1}\abs*{\tri*{T,v^{\otimes{}p}}}$, and
we let $T[v\ind{1},\ldots,v\ind{p}]=\tri*{T,v\ind{1}\otimes\cdots\otimes{}v\ind{p}}$.
Note that for a vector $x\in\bbR^{d}$ the operator norm $\nrm*{x}_\op$
coincides with the Euclidean norm $\nrm{x}$. We let $\symd$ denote the
space of symmetric matrices in $\bbR^{d\times{}d}$. We let $\bbB_r(x)$ denote
the Euclidean ball of radius $r$ centered at $x\in\bbR^{d}$ (with
dimension clear from context). We adopt non-asymptotic
big-O notation, where $f=O(g)$ for $f,g:\cX\to\bbR_{+}$ if
$f(x)\leq{}Cg(x)$ for some constant $C>0$.

\section{Setup} \label{sec:prelim}

We
study the problem of finding $\eps$-stationary and
$(\eps,\gamma)$-second order stationary points in the standard
 oracle complexity framework \citep{nemirovski1983problem}, which we
briefly review here.

\paragraph{Function classes.}
We consider $p$-times differentiable functions satisfying
standard regularity conditions, and define
\begin{equation*}
  \Fclasstilp = \left\{ F:\R^{d}\to\R \left|
    \begin{array}{l}
      F \in \cC^p, \quad F(0) - \inf_x F(x) \leq \Delta, \\
      \nrm*{\grad^q{}F(x)-\grad{}^qF(y)}_\op \leq{}L_q\nrm*{x-y}
      ~\mbox{for~all~} x, y \in \R^d,~q\in[p]
    \end{array}
    \right. \!\! \right\},
\end{equation*}
so that $\Lip{1:p}\defeq(L_1,\dots,L_p)$ specifies the Lipschitz constants
of the $q$th order derivatives $\nabla^q F$ with respect to the operator norm.	We make no restriction on the ambient dimension $d$.
	\paragraph{Oracles.}
        For a given function $F\in\Fclasstilp$, we consider a class of stochastic
        \pth order oracles defined by a distribution $P_z$ over a measurable
        set $\cZ$ and an estimator
	\begin{equation} \label{eq:oracle_answer}
  		\estimator(x,z) \defeq \prn*{\wh{F}(x,z), \stocder{x,
  		z}{}, \stocder{x, z}{2}, \ldots, \stocder{x, z}{p}},
	\end{equation}
        where
	$\crl{\stocder{\cdot, z}{q}}_{q=0}^{p}$ are unbiased estimators of the respective
	derivatives. That is, for all $x$, {$\En_{z\sim{}P_z}[\wh{F}(x,z)]=F(x)$} and $\En_{z
	\sim P_z}[\stocder{x, z}{q}] = \grad^q
	F(x)$ for all $q\in[p]$.\footnote{For $p\ge2$ we assume without
	loss of generality that $\stocder{x,z}{p}$ is a symmetric tensor.}

	Given variance parameters $\sigma_{1:p} = (\sigma_1,
	\ldots, \sigma_{p})$, we define the \textit{oracle class}
	$\stocpOracle$ to be the set of all stochastic
	\pth-order oracles for which the variance of the derivative estimators
	satisfies
	\begin{equation} \label{eq:oracle_variance}
	\En_{z\sim P_z}\,{\nrm*{
	 \stocder{x, z}{q} - \nabla^q F(x) }_\op^2} \leq
	 \sigma_q^2,~~q\in [p].
       \end{equation}
       The upper bounds in this paper hold even when
       $\sigma_0^{2}\ldef{}\max_{x\in\bbR^{d}}\mbox{Var}(\wh{F}(x,z))$ is infinite,
       while our lower bounds hold when $\sigma_0=0$, so
       to reduce notation, we leave dependence on this parameter
       tacit.
	\paragraph{Optimization protocol.}
	We consider stochastic \pth-order optimization algorithms that access an
	unknown
	function $F\in\Fclasstilp$ through multiple rounds of queries to a
	stochastic \pth-order oracle $(\oracle_F^{p},P_z) \in \stocpOracle$.
	When queried at $x\ind{t}$ in round $t$, the oracle performs an
        independent draw of $z^{(t)} \sim P_z$ and answers with
        $\oracle_F^{p}(x\ind{t},z\ind{t})$. Algorithm queries depend on $F$
        only through the oracle answers; see e.g.\ \citet[Section
        2]{arjevani2019lower} for a more formal treatment.

\section{Complexity of finding first-order stationary points}
\label{sec:epsilon}

In this section we focus on the task of finding $\epsilon$-approximate
stationary points (satisfying $\norm{\grad F(x)}\le \epsilon$). As prior work
observes~\citep[cf.][]{carmon2017convex,allen2018how}, stationary point
search is a useful primitive for achieving the end goal of finding
second-order stationary points~\eqref{eq:sosp}. We begin with describing
algorithmic upper bounds on the complexity of finding stationary points
with stochastic second-order oracles, and then proceed to match their
leading terms with general $p$th order lower bounds.

\subsection{Upper bounds} \label{sec:epsilon_upper}

Our algorithms rely on \emph{recursive variance
reduction}~\citep{nguyen2017sarah}:
we sequentially estimate the
gradient at the points $\{x\ind{t}\}_{t\ge0}$ by accumulating cheap
estimators
 of $\grad F(x\ind{\tau})-\grad
F(x\ind{\tau-1})$ for $\tau=t_0+1,\ldots, t$, where at iteration $t_0$ we
reset the gradient estimator by computing a high-accuracy approximation
of $\grad F(x\ind{t_0})$ with many oracle
queries. Our implementation of recursive variance reduction,
\pref{alg:gradient_estimator}, differs from previous approaches~\citep{fang2018spider,zhou2018stochastic,wang2018spiderboost} in
three aspects.
  \colt{

 First, in \pref{line:hvp_update_step} we estimate differences of
	the form $\grad F(x\ind{\tau})-\grad F(x\ind{\tau-1})$ by averaging
	stochastic Hessian-vector products. This allows us to do away
        with multi-point queries and operate under
	weaker assumptions than prior work (see \pref{app:oracle}), but it also
	introduces bias to our estimator, which makes its analysis more involved.
	This is the key novelty in our algorithm.
         Second, rather than resetting the gradient estimator
	every fixed number of steps, we reset with a
	user-defined probability $\bias$ (\pref{line:initialization_est});
	this makes the estimator stateless and greatly simplifies its analysis,
	especially in our algorithms for finding SOSPs, where we use a
        varying value of $\bias$. Finally, we dynamically select
	the batch size $\numbabysteps$ for estimating gradient differences
	based on the distance between iterates (\pref{line:rvr-params}), while
	prior work uses a constant batch size. Our dynamic batch size scheme is
	crucial for controlling the bias in our  estimator,
        while still allowing for large step sizes as in~\citet{wang2018spiderboost}.
    }
\arxiv{
  \begin{enumerate}
	\item In \pref{line:hvp_update_step} we estimate differences of
	the form $\grad F(x\ind{\tau})-\grad F(x\ind{\tau-1})$ by averaging
	stochastic Hessian-vector products. This allows us to do away
        with multi-point queries and operate under
	weaker assumptions than prior work (see \pref{app:oracle}), but it also
	introduces bias to our estimator, which makes its analysis more involved.
	This is the key novelty in our algorithm.
	\item Rather than resetting the gradient estimator
	every fixed number of steps, we reset with a
	user-defined probability $\bias$ (\pref{line:initialization_est});
	this makes the estimator stateless and greatly simplifies its analysis,
	especially when we use a varying value of $\bias$ to find second-order
	stationary points.
	\item We dynamically select
	the batch size $\numbabysteps$ for estimating gradient differences
	based on the distance between iterates (\pref{line:rvr-params}), while
	prior work uses a constant batch size. Our dynamic batch size scheme is
	crucial for controlling the bias in our gradient estimator,
        while still allowing for large step sizes as in~\citet{wang2018spiderboost}.
      \end{enumerate}
      }

The core of our analysis is the following lemma, which bounds the
gradient estimation error and expected oracle complexity. To
state the lemma, we let $\{x\ind{t}\}_{t\ge0}$ be sequence of queries to
  \pref{alg:gradient_estimator}, and let ${g\ind{t} =
  \getgradientfunc_{\epsilon,b}(x\ind{t}, x\ind{t-1}, g\ind{t-1})}$
be the sequence of estimates it returns.
\begin{lemma}\label{lem:base-var-bound}
	For any oracle in  $\cO_{2}(F, \sigma_{1:2})$ and
  $F \in \cF_2\prn{\Delta, L_{1:2}}$, \pref{alg:gradient_estimator} guarantees that
	\begin{equation*}
	\E\, \norm{g\ind{t} - \grad F(x\ind{t})}^2 \le \epsilon^2
      \end{equation*}
       for all $t\geq{}1$.
Furthermore, conditional on $x\ind{t-1}$, $x\ind{t}$ and $g\ind{t-1}$, the
$t^\textup{th}$
 execution of \pref{alg:gradient_estimator}  with reset probability $\bias$
uses at most
\[
O\Big(1 + \bias \frac{\sigma_1^2}{\epsilon^2} + {\nrm[\big]{x\ind{t} -
x\ind{t-1}}^2 \cdot \frac{
{\sigma_2^2 + \epsilon L_2}}{\bias\epsilon^2}}\Big)
\]
stochastic gradient and Hessian-vector product queries in expectation.
\end{lemma}

We prove the lemma in \pref{app:estimators} by bounding the
per-step variance using the HVP oracle's variance
bound~\pref{eq:oracle_variance}, and by bounding the per-step bias
relative to $\grad
F(x\ind{t})-\grad F(x\ind{t-1})$ using the Lipschitz continuity of the
Hessian.

\colt{
\newlength{\textfloatsepsave}
\setlength{\textfloatsepsave}{\floatsep}
\setlength{\floatsep}{6pt}
\setlength{\textfloatsep}{0pt}
}

\begin{algorithm}[t]
	\caption{Recursive variance reduction with stochastic
		Hessian-vector products (\mainalg)}
	\label{alg:gradient_estimator}
	\begin{algorithmic}[1]
		\Statex{} \algcomment{Gradient estimator for $F \in \cF_2(\Delta, L_{1:2})$ given stochastic oracle in  $\cO_{2}(F, \sigma_{1:2})$.}
		\Function{\getgradientfunc$_{\newestparams,\bias}$}{$x$, $\xprev$,
			$\gprev$}:
		\State  Set $\numbabysteps  = \ceil*{\frac{5 \prn*{ \sigma_2^2 +
					L_2\epsilon}}{\bias \epsilon^2}\cdot \norm{x-\xprev}^2}$ and
					$n =
		\ceil*{\frac{5
				\sigma_1^2}{\epsilon^2}}$.\label{line:rvr-params}
		\State Sample $\coin \sim \bern(\bias)$.\label{line:coin-toss}
		\If {$\coin$ is  $1$ or $\gprev$ is $\bot$}{}
		\label{line:initialization_est}
		\State\label{line:fresh_start} Query the oracle $\outerbatchsize $
		times at $x$ and set $g \leftarrow \frac{1}{n}
		\sum_{j=1}^{\outerbatchsize} \stocder{x, z\ind{j}}{}$, where
		$z\ind{j}\overset{\mathrm{i.i.d.}}{\sim}P_z.$
		\Else
		\State Define $x\ind{k} \ldef{} { \frac{k}{K}  x +
			\prn*{1 - \frac{k}{K}} \xprev}$ for $k \in \crl{0,\ldots,K}$.
		\State \label{line:hvp_update_step} Query the oracle at the set of
		points $\prn*{x\ind{k}}_{k=0}^{K-1}$ to compute
		\item[]  {\begin{center}
				$g \leftarrow \gprev +  \sum_{k=1}^K \stocder{x\ind{k-1},
				z\ind{k}}{2} \prn*{x\ind{k} - x\ind{k-1}},\quad\text{where }
				z\ind{k}\overset{\mathrm{i.i.d.}}{\sim}P_z.$
		\end{center}}
		\EndIf
		\State \Return $g$.
		\EndFunction
	\end{algorithmic}
\end{algorithm}

\begin{algorithm}[t]
	\caption{Stochastic gradient descent with \mainalg}
	\label{alg:epsilon_hvp}
	\begin{algorithmic}[1]
		\Require
		Oracle $(\estimator[2], P_z) \in
		\stocpOracle[2]$ for $F \in \cF_2\prn*{\Delta, L_1, L_2}$.
		Precision parameter $\epsilon$.
		\State Set $\eta = \frac{1}{2 \sqrt{ L_1^2 + \sigma_2^2 + \epsilon
				L_2}}$, $T = \ceil*{\frac{2\Delta}{\eta \epsilon^2}}$, $\bias=\min \crl*{1,
				\frac{\eta \epsilon \sqrt{\sigma_2^2 + \epsilon
				L_2}}{\sigma_1}}$. \label{line:parameter_setting}
		\State Initialize $x\ind{0}, x\ind{1} \leftarrow 0$, $\gradest[0]
		\leftarrow \bot$.
		\For {$t = 1 ~\text{to}~ T$}
		\State $\gradest[t] \leftarrow$ $\getgradientfunc_{\newestparams,
		\bias}(x
		\ind{t}, x\ind{t-1}, \gradest[t-1])$.
		\label{line:gradient_estimate_epsilon_hvp}
		\State $x \ind{t+1} \leftarrow x \ind{t} - \eta \gradest[t]$.
		\label{line:hvp_update_rule}
		\EndFor
		\State \textbf{return} $\hx$ chosen uniformly at random from
		$\crl*{x\ind{t}}_{t=1}^{T}$.
	\end{algorithmic}
\end{algorithm}

Our first algorithm for finding $\eps$-stationary points, \pref{alg:epsilon_hvp}, is simply stochastic gradient descent using the
\mainalg gradient estimator (\pref{alg:gradient_estimator}); we bound its
complexity by
$O(\eps^{-3})$. Before stating the result formally, we briefly sketch the
analysis here\colt{ (see \pref{app:epsilon_hvp} for details)}. Standard
analysis of SGD with step size
$\eta\le \frac{1}{2L_1}$ shows that its iterates satisfy $\E \norm{\grad
F(x\ind{t})}^2
\le \frac{1}{\eta}\E [F(x\ind{t+1}) - F(x\ind{t})] + O(1) \cdot \E\,
\norm{g\ind{t}
- \grad F(x\ind{t})}^2$. Telescoping over $T$ steps, using
\pref{lem:base-var-bound} and substituting in the
initial
suboptimality bound $\Delta$, this implies that
\begin{equation} \label{eq:sgd-variance-bound}
\frac{1}{T} \sum_{t=0}^{T-1} \En{} \norm{\grad F(x\ind{t})}^2 \le
\frac{\Delta}{\eta T} + O(\eps^2).
\end{equation}
Taking  $T=\Omega(\frac{\Delta}{\eta \eps^2})$, we are guaranteed that
a uniformly selected iterate has expected norm $O(\veps)$.

\colt{\setlength{\textfloatsep}{\textfloatsepsave}}

To account for oracle complexity, we observe from
\pref{lem:base-var-bound} that $T$ calls to
\pref{alg:gradient_estimator} require at most $T(\frac{\sigma_1^2
\bias}{\epsilon^2} + 1) +
\sum_{t=1}^{T} \En {\norm{x\ind{t}-x\ind{t-1}}^2} \cdot
\big(\frac{\sigma_2^2 + L_2 \epsilon}{\bias\epsilon^2} \big)$ oracle
queries in expectation. Using $x\ind{t}-x\ind{t-1}=\eta g\ind{t-1}$,
\pref{lem:base-var-bound} and~\pref{eq:sgd-variance-bound} imply that
$\sum_{t=1}^{T} \En {\norm{x\ind{t}-x\ind{t-1}}^2} \le O(T\epsilon^2)$.
We
then choose
$\bias$ to out the  terms $T\big({\frac{\sigma_1^2
\bias}{\epsilon^2}}\big)$ and $T\big({\frac{\sigma_2^2 + L_2\epsilon}{b}}\big)$.
\arxiv{
This gives the
following
complexity guarantee,
which we prove in~\pref{app:epsilon_hvp}.
}
\colt{
  This gives the
following
 guarantee.
}

\begin{theorem}\label{thm:epsilon_hvp}
 For any function $F \in \cF_2\prn*{\Delta, L_1, L_2}$, stochastic
 second-order oracle in $\cO_{2}(F, \sigma_1, \sigma_2)$, and $\epsilon <
\min \crl*{\sigma_1, \sqrt{\Delta{}L_1}}$, with probability at least
$\frac{3}{4}$, \pref{alg:epsilon_hvp}  returns a point $\hx$ such that $
 \nrm*{\grad F(\hx)} \leq \epsilon$ and performs at most
 \begin{equation*}
 O\Big(
 \frac{\Delta \sigma_1 \sigma_2}{\epsilon^3} + \frac{\Delta
 L_2^{0.5}\sigma_1 }{\epsilon^{2.5}}
 + \frac{\Delta L_1}{\epsilon^2}
\Big)
 \end{equation*}
stochastic gradient and Hessian-vector product queries.
\end{theorem}

The oracle complexity of \pref{alg:epsilon_hvp} depends on the Lipschitz
parameters of $F$ only through lower-order terms in $\epsilon$, with the
leading term scaling only with the variance of the gradient and Hessian
estimators. In the low noise regime
where $\sigma_1 < \epsilon$ and $\sigma_2 < \max\{L_1,\sqrt{L_2\eps}\}$, the
complexity becomes $O(\Delta L_1 \epsilon^{-2} +
\Delta L_2^{0.5}\epsilon^{-1.5})$ which is simply the maximum of the noiseless guarantees for gradient descent and Newton's method. We remark,
however, that in the noiseless regime $\sigma_1=\sigma_2=0$, a slightly
better
guarantee $O(\Delta L_1^{0.5}L_2^{0.25} \epsilon^{-1.75} +
\Delta L_2^{0.5}\epsilon^{-1.5})$ is achievable~\citep{carmon2017convex}.

In the noiseless setting, any algorithm that
uses only first-order and Hessian-vector product queries must have
complexity scaling with $L_1$, but full Hessian access can remove this
dependence \citep{carmon2019lower_ii}. We show that
the same holds true in the stochastic setting:
\pref{alg:epsilon_cubic}, a subsampled cubic regularized trust-region
method using \pref{alg:gradient_estimator} for gradient estimation,
enjoys a complexity bound independent of $L_1$. \arxiv{We defer the
analysis to \pref{app:upper_bounds_e_L2} and state the guarantee as
follows.}\colt{We defer the
full description and analysis to \pref{app:upper_bounds_e_L2} and state the guarantee here.}
\begin{theorem}\label{thm:epsilon_cubic}
	For any function $F \in \cF_2\prn*{\Delta, \infty, L_2}$, stochastic
	second
	order oracle in $\cO_{2}(F, \sigma_1, \sigma_2)$, and $\epsilon <
	\sigma_1$, with probability at least $\frac{3}{4}$, \pref{alg:epsilon_cubic} returns a point $\hx$ such that
	$ \nrm*{\grad F(\hx)} \leq \epsilon$ and performs at most
	\begin{equation*}
	O\Big(
		\frac{\Delta \sigma_1 \sigma_2 }{\epsilon^3}\cdot{}\log^{0.5}d + \frac{\Delta
			L_2^{0.5}\sigma_1 }{\epsilon^{2.5}}
	\Big)
	\end{equation*}
	stochastic gradient and Hessian queries.
\end{theorem}

The guarantee of \pref{thm:epsilon_cubic} constitutes an improvement in
query complexity over
\pref{thm:epsilon_hvp} in the regime $L_1 \gtrsim
(1+\frac{\sigma_1}{\eps})(\sigma_2 + \sqrt{L_2 \epsilon})$. However,
depending on the problem, full stochastic Hessians can be up to $d$ times
more expensive to compute than stochastic Hessian-vector products.

\arxiv{
\begin{algorithm}[t]
   \caption{Subsampled cubic-regularized trust-region method with
     \mainalg}
	\label{alg:epsilon_cubic}
	\begin{algorithmic}[1]
		\Require
		 Oracle $(\estimator[2], P_z) \in
\stocpOracle[2]$ for $F \in \cF_2\prn*{\Delta, \infty, L_2}$.
		 Precision parameter $\epsilon$.
		\State Set $M = 5 \max\crl*{L_2, \frac{\epsilon\sigma_2^2\log (d)}{\sigma_1^2}}$, $\eta =  25 \sqrt{\frac{\epsilon}{M}}$, $T =
		\ceil*{\frac{5 \Delta}{3 \eta \epsilon}}$ and $n_H = \ceil*{\frac{22 \sigma_2^2 \eta^2 \log(d)}{\epsilon^2}}$.
		\State Set $\bias=\min\crl*{1, \frac{\eta \sqrt{\sigma_2^2 + \epsilon L_2}}{25 \sigma_1}}$.
		\State Initialize $x\ind{0}, x\ind{1} \leftarrow 0$,~ $\gradest[0] \leftarrow \bot$.
		\For {$t = 1 ~\text{to}~ T$}
		\State \label{line:Hessian_estimator}Query the oracle $n_H$ times at
		$x \ind{t}$ and compute
		\begin{equation*}
		\hessest[t] \leftarrow \frac{1}{n_H} \sum_{j=1}^{n_H}
		\stocder{x\ind{t}, z\ind{t, j}}{2},\quad\text{where}\quad{}z\ind{t,
		j}\overset{\mathrm{i.i.d.}}{\sim}P_z.
		\end{equation*}
		\State $\gradest[t] \leftarrow
		\getgradientfunc_{\newestparams, \bias}\prn*{x\ind{t}, x\ind{t-1}, g\ind{t-1}}$.\label{line:gradient_estimation_cubic}
		\State\label{line:constrained_cubic_step}Set the next point $x\ind{t+1}$ as
	\begin{equation*}
			x\ind{t+1} \gets \argmin_{y:\norm{y-x\ind{t}}\le\eta}
		\tri[\big]{\gradest[t], y - x\ind{t}} + \frac{1}{2}\tri[\big]{y-x\ind{t},
		\hessest[t](y - x\ind{t})} + \frac{M}{6}
		\nrm[\big]{y-x\ind{t}}^{3}.
	\end{equation*}
		\EndFor
		\State \textbf{return} $\hx$ chosen uniformly at random from
		$\crl*{x\ind{t}}_{t=2}^{T + 1}.$
\end{algorithmic}
\end{algorithm}

}

\subsection{Lower bounds} \label{sec:epsilon_lower}

	Having presented stochastic second-order methods with
	$O(\eps^{-3})$-complexity bound for finding $\eps$-stationary
	points, our we next show that this rates cannot be
	improved. In fact, we show that this rate is optimal even when one is
	given access to stochastic higher derivatives of \emph{any} order.
	We prove our lower bounds for the class of \emph{zero-respecting}
	algorithms, which subsumes the majority of existing
	optimization methods; see \pref{app:epsilon_lower} for a
        formal definition. We believe that existing
	techniques~\citep{carmon2019lower_i,arjevani2019lower} can
	strengthen our lower bounds to apply to general randomized
	algorithms; for brevity, we do not pursue it here.

	The lower bounds in this section closely follow a recent construction
	by~\citet[Section 3]{arjevani2019lower}, who prove lower bounds for
	stochastic first-order methods. To establish complexity bounds for
	$p$th-order methods, we extend the `probabilistic zero-chain' gradient
	estimator introduced in~\citet{arjevani2019lower} to high-order derivative
	estimators.The most technically demanding part of our proof is a careful scaling of
	the basic construction to simultaneously meet multiple Lipschitz continuity
	and variance constraints. Deferring the proof details
	to~\pref{app:epsilon_lower}, our
	lower bound is as follows.

	\begin{theorem}
		\label{thm:eps_so_zero_respecting}
			For all $p\in\N$, $\Delta, \Lip{1:p}, \sigma_{1:p}>0$ and $\eps \le
			O(\sigma_1)$, there exists $F \in \cF_p\prn*{\Delta, L_{1:p}}$
			and
			$(\estimator[p], P_z) \in \stocpOracle[p]$, such that for any
			\pth-order zero-respecting algorithm, the number of
			 queries required to obtain an
			$\eps$-stationary point with constant
                        probability is bounded from below by
			\begin{align}\label{eq:minimax_eps_p}
			\Omega(1)\cdot \frac{\Delta\sigma_1^2}{\epsilon^3}
			\min\crl*{
				\min_{q\in\{2,\ldots,p\}}\prn*{\frac{\sigma_q}{\sigma_1}}^{\frac{1}{q-1}},
				\min_{q'\in\{1,\ldots,p\}}\prn*{\frac{L_{q'}}{\eps }}^{1/q'}}.
			\end{align}
			A construction of dimension
			$\Theta\Bigl(\frac{\Delta
			}{\epsilon}
			\min\Bigl\{
			\min_{q\in\{2,\ldots,p\}}\prn*{\frac{\sigma_q}{\sigma_1}}^{\frac{1}{q-1}},
			\min_{q'\in\{1,\ldots,p\}}\prn*{\frac{L_{q'}}{\eps }}^{1/q'}\Bigr\}\Bigr)$ realizes
			this lower bound.
	\end{theorem}
			For second-order methods (with $p=2$),
		\pref{thm:eps_so_zero_respecting}
		specializes to the oracle complexity lower bound
		\begin{align}\label{eq:minimax_eps_2}
		\Omega(1)\cdot
		\min\crl*{
			\frac{\Delta\sigma_1\sigma_2}{\epsilon^3},
				\frac{\Delta L_2^{0.5}\sigma_1}{\epsilon^{3.5}},
				\frac{\Delta L_1 \sigma_1^2}{\epsilon^{4}}},
		\end{align}
		which is tight in that it matches (up to
		numerical constants) the convergence rate of~\pref{alg:epsilon_hvp}
		in the regime where
		$\Delta
		\sigma_1\sigma_2\eps^{-3}$ dominates both the upper bound
		in~\pref{thm:epsilon_hvp} and expression~\pref{eq:minimax_eps_2}. The lower
		bound~\pref{eq:minimax_eps_2} is also tight when the second-order
		information is not available or reliable
                ($\sigma_2$ is infinite or very
		large, respectively): Standard SGD matches the $\eps^{-4}$
		term~\citep{ghadimi2013stochastic}, while more
                sophisticated variants based on
                restarting~\citep{fang2019sharp} and normalized
                updates with momentum \citep{cutkosky2020momentum} match the
		$\eps^{-3.5}$ term
		(the former up to logarithmic factors)---neither of
		these algorithms requires stochastic second derivative estimation.

		\pref{thm:eps_so_zero_respecting} implies that while  higher-order
		methods (with $p>2$) might achieve better dependence on the
		variance parameters than the upper bounds
		for~\pref{alg:epsilon_hvp} or \pref{alg:epsilon_cubic}, they cannot improve the
		$\eps^{-3}$ scaling. This highlights a fundamental limitation for higher-order methods in
		stochastic non-convex optimization which does not exist in the
		noiseless case. Indeed, without noise the optimal rate for finding
		$\epsilon$-stationary point with a $p$th order method is
		$\Theta(\eps^{-1+\frac{1}{p}})$~\cite{carmon2019lower_i}; we
		illustrate this contrast in \pref{fig:elbow}.

       Altogether, the results presented in this section
                fully characterize (with respect to dependence on
                $\eps$) the complexity of finding $\eps$-stationary
                points with stochastic second-order methods and beyond in the
                single-point query model. We briefly remark that lower
                bound in \pref{eq:minimax_eps_p} immediately extends
                to multi-point queries, which shows that even
                second-order methods offer little benefit once two or
                more simultaneous queries are allowed.

\section{Complexity of finding second-order stationary points}
Having established rates of convergence for finding $\eps$-stationary
points, we now turn our attention to $(\eps,\gamma)$-second order
stationary points, which have the additional
requirement that $\lambda_{\min}(\hess F(x)) \ge -\gamma$, i.e.\ that $F$
is $\gamma$-weakly convex around $x$. This section follows the general
organization of the prequel: we first design and analyze an algorithm
with improved upper bounds, and then develop nearly-matching lower bounds that apply to a
broad class of algorithms.

\label{sec:epsilon_gamma}
\arxiv{\subsection{Upper bounds}}
\colt{\subsection{Upper bounds}}
\label{sec:epsilon_gamma_upper}

Our first contribution for this section is an algorithm that enjoys improved
complexity for finding $(\eps,\gamma)$-second-order stationary points,
and that achieves this using only stochastic gradient and Hessian-vector product queries.
To guarantee second-order stationarity, we follow the established technique
 of interleaving an algorithm for finding a first-order stationary point with
 negative curvature descent~\citep{carmon2017convex, allen2018how}.
 However, we employ a randomized variant of this approach. Specifically, at every iteration we
 flip a biased coin to determine whether to perform a stochastic
 gradient step or a stochastic negative curvature descent step.

Our algorithm estimates stochastic gradients using the \mainalg scheme
(\pref{alg:gradient_estimator}), where the value of the restart probability
$\bias$ depends on the type of the previous step (gradient or negative
curvature). To implement negative curvature
descent,
we apply Oja's
 method~\citep{oja1982simplified,allen2017follow} which detects
 directions of negative curvature using only stochastic Hessian-vector
 product queries. For technical reasons pertaining to the analysis of Oja's method, we require
 the stochastic
 Hessians to be bounded almost surely, i.e.,
 $\opnorm{\stochess{x,z}-\hess
 	F(x)}\le\bar{\sigma}_2$ a.s.; we let $\stocpOracleas$ denote
 	the  class of such bounded noise oracles. Under this assumption,
 \pref{alg:epsilon_gamma_hvp}---whose description is deferred to the \pref{app:upper_eg}---enjoys the following convergence
 guarantee.\footnote{The notation $\Otilde(\cdot)$ hides lower-order terms and logarithmic dependence on the dimension $d$. See the proof in \pref{app:upper_eg} for the complete description of the algorithm and the full
  complexity bound, including lower order terms.}

\newcommand{\sigt}{\tilde{\sigma}_2}
\begin{theorem}
\label{thm:epsilon_gamma_hvp} For any function $F \in \cF_2(\Delta, L_{1:2})$, stochastic Hessian-vector product oracle in $\stocpOracleas$,
$\epsilon \leq \min\crl*{\sigma_1, \sqrt{\Delta L_1}}$, and $\gamma \leq \min\crl*{\bar{\sigma}_2, L_1,
\sqrt{\epsilon{}L_2}}$, with probability at least $\frac{5}{8}$ \pref{alg:epsilon_gamma_hvp} returns a point $\xhat$ such that
\[
  \nrm*{\grad \func{\xhat}} \leq \epsilon \quad\text{
   and} \quad\lambda_{\min{}} \prn*{\grad^2 \func{\xhat}} \geq -\gamma,
\] and performs at most
 \begin{align*}
\widetilde{O} \prn*{\frac{  \Delta \sigma_1 \bsigma_2}{\epsilon^3}  +
\frac{\Delta L_2 \sigma_1\bsigma_2}{\gamma^2 \epsilon^2} +
\frac{\Delta L^2_2 \prn*{\bsigma_2 + L_1}^2}{\gamma^5} + \frac{\Delta
L_1}{\epsilon^2}}
\end{align*}
 stochastic gradient and Hessian-vector product queries.
\end{theorem}

Similar to the case for finding $\eps$-stationary points (see
discussion preceding \pref{thm:epsilon_cubic}), using full stochastic
Hessian information allows us to design an algorithm (\pref{alg:eg_cubic})
which removes the dependence on $L_1$ from the theorem above. Moreover, estimating negative curvature directly from empirical Hessian
estimates saves us the need to use Oja's method, which means that we do
not need the additional boundedness assumption on the stochastic
Hessian used by \pref{alg:epsilon_gamma_hvp}. \arxiv{We defer the complete description and
analysis for \pref{alg:eg_cubic} to
\pref{app:eg_cubic_upper}, and state its complexity guarantee
below.}\colt{We defer the complete description, complexity guarantee,
and for analysis for \pref{alg:eg_cubic} to
\pref{app:eg_cubic_upper}.}

\arxiv{
\begin{theorem}
\label{thm:eg_cubic_upper} For any function $F \in \cF_2(\Delta, \infty,
L_2)$, stochastic second order oracle in $\cO_{2}(F, \sigma_1, \sigma_2)$,
$\epsilon \leq \sigma_1$, and $\gamma \leq \min\crl[\big]{\sigma_2,
\sqrt{\epsilon{}L_2}, \Delta^{\frac{1}{3}} L_2^\frac{2}{3}}$, with probability at least $\frac{3}{5}$ \pref{alg:eg_cubic} returns a point $\xhat$ such that
   \[
   \lambda_{\min{}} \prn*{\grad^2 \func{\xhat}} \geq -\gamma, \quad\text{
   and }\quad\nrm*{\grad \func{\xhat}} \leq \epsilon,
\]
 and performs at most
 \begin{align*}
\widetilde{O} \prn*{\frac{  \Delta \sigma_1 \sigma_2}{\epsilon^3}  +
\frac{\Delta L_2 \sigma_1\sigma_2}{\gamma^2 \epsilon^2} +  \frac{\Delta
L^2_2 \sigma_2^2}{\gamma^5}}
\end{align*}
 stochastic gradient and Hessian queries.
\end{theorem}
}

\arxiv{\subsection{Lower bounds}}
\colt{\subsection{Lower bounds}}
\label{sec:epsilon_gamma_lower}

	\newcommand{\dsunscaled}{\Sigma_{2T}}
        We now develop lower complexity bounds for the
        task of finding $(\eps,\gamma)$-stationary points. To do so,
        we prove new lower bounds for the simpler sub-problem of finding a \emph{$\gamma$-weakly convex} point, i.e., a point $x$ such
	that $\lambdamin(\grad^2 F(x))\ge-\gamma$ (with no restriction on
	$\norm{\grad F(x)}$). Lower bounds for finding $(\eps,\gamma)$-SOSPs
	 follow as the maximum (or, equivalently, the sum) of lower
	bounds we develop here and the lower bounds for finding $\veps$-stationary
	points given in~\pref{thm:gamma_so_zero_respecting}. To see why this is
	so, let
	$F_\eps$ and $F_\gamma$ be hard instances for finding
	$\epsilon$-stationary and $\gamma$-weakly-convex points
	respectively, and consider the ``direct sum''
	$F_{\epsilon,\gamma}(x)\ldef\frac{1}{2}F_\epsilon(x_1, \ldots, x_d) +
	\frac{1}{2}F_\gamma(x_{d+1},\ldots,x_{2d})$; this is a hard instance for
	finding $(\epsilon,\gamma)$-SOSPs that inherits all the regularity
	properties of its constituent functions.

	The basic construction we use here is a modification of the zero-chain
	introduced in~\citet{carmon2019lower_i} (see
        \pref{eq:hard_function_eps_LB} in \pref{app:lower})
	in 	which large
	$\lambdamin(\grad^2 F(x))$
	is possible only when essentially none of
	the entries of $x$ is zero. Given $T > 0$, we define the hard function
	\begin{equation}
	\barF(x) \defeq \Psi(1)\Phibar(x_1) + \sum_{i=2}^T ~ \brk*{
	\Psi(-x_{i-1})
        \Phibar(-x_i) + \Psi(x_{i-1}) \Phibar(x_i)) },
      	\label{eq:hard_function_gamma_LB}
      \end{equation}
	where $\Psi(x)\ldef\exp(1-\frac{1}{(2x-1)^2})\indicator{x>\frac{1}{2}}$
	(as in~\citet{carmon2019lower_i}) and $\Phibar(x)\defeq  8\prn{
	e^{\frac{-x^2}{2}}-1}$.

	Our design for the function $\Phibar$ guarantees that any query whose last coordinate
	is zero has significant negative curvature, while maintaining the original
	chain structure which guarantees that zero-respecting algorithms require
	many queries before ``discovering'' the last coordinate. We complete the
	construction by specifying a collection of stochastic derivative estimators similar to
	those in~\pref{sec:epsilon_gamma_lower}, except for that we
        choose the stochastic gradient estimator $\wh{\grad
	\barFunscaled}$ to be exactly equal to $\grad \barFunscaled$, so
      that the lower bound holds even for $\sigma_1=0$;
	Appropriately scaling $\barFunscaled$ allows us to tune the
	Lipschitz constants of its derivatives and the variance of the
	estimators, thereby	establishing the following complexity bounds (see
	\pref{app:gamma_lower} for a full derivation).
		\begin{theorem}
			\label{thm:gamma_so_zero_respecting}
		Let $p\ge 2$ and $\Delta, \Lip{1:p}, \sigma_{1:p}>0$ be fixed.
		 If $\gamma \le
		O(\min\{\sigma_2,L_1\})$, then there exists ${F \in \cF_p\prn*{\Delta,
		L_{1:p}}}$
		and
		$(\estimator[p], P_z) \in \stocpOracle[p]$ such that for any
		stochastic \pth-order zero-respecting algorithm, the number of
		queries to $\estimator[p]$required to obtain a
		$\gamma$-weakly convex point with constant probability
                is at least
		\begin{align}\label{eq:minimax_gamma_p}
			\Omega(1)\cdot
			\begin{cases}
			\frac{\Delta \sigma_2^2 L_2^2}{\gamma^5}, & p=2, \\
			\frac{\Delta \sigma_2^2
			}{ \gamma^3 }
			\min\crl*{
				\min_{q\in\{3,\ldots,p\}}\prn*{\frac{\sigma_q}{\sigma_2}}^{\frac{2}{q-2}},
				\min_{q'\in\{2,\ldots,p\}}\prn*{\frac{\Lip{q'}}{\gamma}}^{\frac{2}{q'-1}}},
			& p>2.\\
			\end{cases}
		\end{align}
		\arxiv{A construction of dimension
		$\Theta\Bigl(\frac{\Delta}{\gamma}
				\min\Bigl\{
					\min_{q\in\{3,\ldots,p\}}\prn*{\frac{\sigma_q}{\sigma_2}}^{\frac{2}{q-2}},
					\min_{q'\in\{2,\ldots,p\}}\prn*{\frac{\Lip{q'}}{\gamma}}^{\frac{2}{q'-1}}\Bigr\}\Bigr)$
			 realizes the lower bound.}
	\end{theorem}
        \arxiv{
        \pref{thm:gamma_so_zero_respecting} is new even in the
        noiseless case  (in which $\sigma_{1} = \dots =
        \sigma_p=0$), where it specializes to
        \begin{equation}\label{eq:minimax_gamma_p_det}
          \Omega(1)\cdot\frac{\Delta
          }{ \gamma }
          \min_{q\in\{2,\ldots,p\}}\prn*{\frac{\Lip{q}}{\gamma}}^{\frac{2}{q-1}}.
        \end{equation}
For the class $\cF_p(\Delta,L_p)$, the lower bound
\pref{eq:minimax_gamma_p_det} further simplifies to
$\Delta{}L_p^{\frac{2}{p-1}}\gamma^{-\frac{p+1}{p-1}}$, which is attained by
the \pth-order regularization method given in \citet[Theorem
3.6]{cartis2017improved}.
}
\colt{
        \pref{thm:gamma_so_zero_respecting} is new even in the
        noiseless case  (in which $\sigma_{1} = \dots =
        \sigma_p=0$), where it specializes to
        $\frac{\Delta
          }{ \gamma }
          \min_{q\in\{2,\ldots,p\}}\prn*{\frac{\Lip{q}}{\gamma}}^{\frac{2}{q-1}}$. For
          the class $\cF_p(\Delta,L_p)$, this further simplifies to
$\Delta{}L_p^{\frac{2}{p-1}}\gamma^{-\frac{p+1}{p-1}}$, which is attained by
the \pth-order regularization method given in \citet[Theorem
3.6]{cartis2017improved}.
}
Together, these results characterize the
deterministic complexity of finding $\gamma$-weakly convex points with
noiseless \pth-order methods.

Returning to the stochastic setting, the bound in \pref{thm:gamma_so_zero_respecting}, when combined with \pref{thm:eps_so_zero_respecting}, implies the
following oracle complexity lower bound bound for finding $(\veps,\gamma)$-SOSPs with
	zero-respecting stochastic second-order methods ($p=2$):
	\begin{align}
		\Omega(1)\cdot\prn*{
	\min\crl*{
		\frac{\Delta\sigma_1\sigma_2}{\epsilon^3},
		\frac{\Delta L_2^{0.5}\sigma_1}{\epsilon^{3.5}},
		\frac{\Delta L_1 \sigma_1^2}{\epsilon^{4}}
	} + \frac{\Delta \sigma_2^2\Lip{2}^2	}{\gamma^5}}.
	\end{align}
	Our lower bound matches the $\epsilon^{-3} + \gamma^{-5}$ terms in
	the upper bound given by~\pref{thm:epsilon_gamma_hvp}, but does not
	match the mixed term $\eps^{-2}\gamma^{-2}$ appearing in the
        upper bound.\footnote{Young's
	inequality only gives $\epsilon^{-3}+\gamma^{-5} \ge
	\Omega(\epsilon^{-9/5}\gamma^{-2})$.} Overall, the rates match whenever
	$\gamma = \Omega(\epsilon^{0.5})$ or $\gamma =O(\epsilon^{2/3})$.

	\pref{thm:gamma_so_zero_respecting} is suggestive of another ``elbow'' phenomenon:
	In the stochastic regime, the rate does not improve beyond
        $\gamma^{-3}$ for $p\ge3$, while the optimal rate in the
        noiseless regime, $\gamma^{-\frac{p+1}{p-1}}$, continues
        improving for all $p$.\footnote{Indeed, when
	high-order noise moments are assumed finite, the term
	$\min_{q\in\{3,\ldots,p\}}\prn*{{\sigma_q}/{\sigma_2}}^{\frac{2}{q-2}}$
	can longer be disregarded. This, in turn, implies that for sufficiently
	small $\gamma$, one cannot improve over $\gamma^{-3}$-scaling, as seen
	by~\pref{eq:minimax_gamma_p}.} However, we are not yet aware of an
        algorithm using stochastic third-order information or higher
        that can achieve the $\gamma^{-3}$ complexity bound.

\arxiv{
\section{Conclusion} \label{sec:conclusion}

This paper provides a fairly complete picture of the worst-case oracle
complexity of finding stationary points with a stochastic second-order
oracle: for $\epsilon$-stationary points we characterize the leading term in
$\eps^{-1}$ exactly and for ($\epsilon,\gamma$)-SOSPs we characterize the
leading term in $\gamma^{-1}$ for a wide range of parameters.
Nevertheless, our results point to a number of open questions.

\paragraph{Benefits of higher-order information for $\gamma$-weakly convex
points.}

Our upper and lower bounds (in \pref{thm:eg_cubic_upper} and
\pref{thm:gamma_so_zero_respecting}) resolve the optimal rate to find an
$(\epsilon, \gamma)$-stationary point for $p=2$, i.e., when $F$ is
second-order smooth and the algorithm can query stochastic gradient and
Hessian information. Furthermore, \pref{thm:eps_so_zero_respecting}
shows that higher order information ($p \geq 3$) cannot improve the
dependence of the rate on the first-order stationarity parameter $\epsilon$.
However, our lower bound for dependence on $\gamma$ scales as
$\gamma^{-5}$ for $p= 2$, but scales as $\gamma^{-3}$ for $p \geq
3$. The weaker lower bound for $p \geq 3$ leaves
open the possibility of a stronger upper bound using third-order
information or higher.

\paragraph{Global methods.} For statistical learning and sample average
approximation problems, it is natural to consider problem instances of the
form $F(x)=\En\brk[\big]{\wh{F}(x,z)}$. For this setting, a more powerful
oracle
model is the \emph{global oracle}, in which samples
$z\ind{1},\ldots,z\ind{n}$ are
drawn i.i.d. and the learner observes the entire function
$\wh{F}(\cdot,z\ind{t})$
for
each $t\in[n]$. Global oracles are more powerful than stochastic \pth order
oracles for every $p$, and lead to improved rates in the
convex setting \citep{foster2019complexity}. Is it possible to beat the $\eps^{-3}$ elbow for such
oracles, or do our lower bounds extend to this setting?

\paragraph{Adaptivity and instance-dependent complexity.}
Our lower bounds show that stochastic higher-order methods cannot
improve the $\eps^{-3}$ oracle complexity attained with stochastic
gradients and Hessian-vector products. Furthermore, in the multi-point
query model, stochastic second-order information does not even lead to
improved rates over stochastic first-order information. However, these
conclusions could be artifacts of our worst-case point of view---are there
natural families of problem instances for which higher-order methods can
adapt to additional problem structure and obtain stronger
instance-dependent convergence guarantees? Developing a theory of
instance-dependent complexity that can distinguish adaptive algorithms
stands out as an exciting research prospect.

 }

\colt{
\subsection*{Discussion}
Due to space constraints, we defer conclusions and discussion to \pref{app:conclusion}.}

\newpage
\subsection*{Acknowledgements}
We thank Blake Woodworth and Nati Srebo for helpful discussions. YA
acknowledges partial support from the Sloan Foundation and Samsung
Research. JCD acknowledges support from the NSF CAREER award
CCF-1553086, ONR YIP N00014-19-2288, Sloan Foundation, NSF HDR
1934578 (Stanford Data Science Collaboratory), and the DAWN Consortium.
DF acknowledges the support of TRIPODS award 1740751.  KS
acknowledges support from NSF CAREER Award 1750575 and a Sloan
Research Fellowship.

\setlength{\bibsep}{6pt}


TITLE: prove that $\frac{a_1+a_2+\dots+a_m}{m}\ge \frac{n+1}{2}$
QUESTION [1 upvotes]: Let $m$ and $n$ be positive integers. Let $a_1, a_2,\dots, a_m$ be distinct elements of $\{1, 2, \ldots, n\}$ such that whenever $a_i + a_j \le n$ for some $i,j, 1 \le i \le j \le m$, there exists $k, 1 \le k \le m$, with $a_i + a_j = a_k$. Prove that
$$\frac{a_1+a_2+\dots+a_m}{m}\ge \frac{n+1}{2}$$
Trial: I am thinking using AM$\ge$GM or $1+2+\dots+n= \frac{n(n+1)}{2}$. But unable to fully solve the problem.

REPLY [1 votes]: Let $S=\{a_1,\ldots, a_m\}$ and $a=\min S$.
For $0\le k<a$, let $$S_k=\{\,x\in S\mid x\equiv k\pmod a\,\}.$$
If $S_k$ is not empty, then because $x\in S_k$ and $a\in S$ implies $x+a\in S$ unless $x+k>n$, its elements form an arithmetic progression of step-size $a$ from $\min S\ge a$ to $\max S\ge n-a+1$. It is well-know that this implies
$$ \sum_{x\in S_k}x=|S_k|\cdot\frac{\min S+\max S}2\ge |S_k|\cdot\frac{a+n-a+1}2=|S_k|\cdot\frac{n+1}2.$$
(This inequality also holds trivially in the case that $S_k=\emptyset$).
Therefore
$$a_1+\ldots+a_m=\sum_{x\in S} x=\sum_{k=0}^{a-1}\sum_{x\in S_k}x\ge\sum_{k=0}^{a-1}|S_k|\cdot \frac{n+1}2=|S|\cdot\frac{n+1}2$$
and after division by $|S|=m$
the original claim.