Stochastic Adaptive Regularization Method with Cubics:
A High Probability Complexity Bound

We are interested in unconstrained optimization problems of the form

where $\phi$ is possibly nonconvex and satisfies the following condition:

We present and analyze a stochastic adaptive cubic regularization algorithm for computing a point $x$ such that $\left\lVert{\nabla\phi(x)}\right\rVert\leq\epsilon$, for some $\epsilon>0$, when $\phi(x)$ or its derivatives are not computable exactly. Specifically, we assume access to stochastic zeroth-, first- and second-order oracles, which are defined as follows.

Stochastic zeroth-order oracle ($\mathsf{SZO}(\epsilon_{f},\lambda,a)$). Given a point $x$, the oracle computes $f(x,\Xi(x))$, where $\Xi(x)$ is a random variable, whose distribution may depend on $x$, $\epsilon_{f},\lambda$ and $a$, that satisfies

for any $t>0$.

We view $x$ as the input to the oracle, $f(x,\Xi(x))$ as the output and the values $(\epsilon_{f},\lambda,a)$ as values intrinsic to the oracle. Thus $|f(x,\Xi(x))-\phi(x)|$ is a sub-exponential random variable with parameters $(\lambda,a)$, whose mean is bounded by some constant $\epsilon_{f}>0$.

Stochastic first-order oracle ($\mathsf{SFO}(\kappa_{g})$. Given a point $x$ and constants $\mu_{1}>0$, $\delta_{1}\in[0,\frac{1}{2})$, the oracle computes $g(x,\Xi^{1}(x))$, such that

where $\Xi^{1}(x)$ is a random variable whose distribution may depend on $x$, $\mu_{1}$, $\delta_{1}$ and $\kappa_{g}$. We view $x$, $\mu_{1}$ and $\delta_{1}$ as inputs to the oracle, while $\kappa_{g}$ is intrinsic to the oracle.

Stochastic second-order oracle ($\mathsf{SSO}(\kappa_{H})$). Given a point $x$ and constants $\mu_{2}>0$, $\delta_{2}\in[0,\frac{1}{2})$, the oracle computes $H(x,\Xi^{2}(x))$, such that

where $\Xi^{2}(x)$ is a random variable whose distribution may depend on $x$, $\mu_{2}$, $\delta_{2}$ and $\kappa_{H}$. The norm on the matrix is the operator norm. $x$, $\mu_{2}$ and $\delta_{2}$ are inputs to the oracle, while $\kappa_{H}$ is intrinsic to the oracle.

Wherever possible, we will omit the dependence on $x$ and write $\Xi,\Xi^{1}$, and $\Xi^{2}$ instead of $\Xi(x),\Xi^{1}(x)$, and $\Xi^{2}(x)$, and we will use $f(x),g(x)$ and $H(x)$ to denote the outputs of the stochastic oracles for brevity.

Similarly, in [BSV14], [GRVZ18], a trust region method is analyzed under essentially $\mathsf{SFO}$ and $\mathsf{SSO}$, but with an exact zeroth-order oracle. Later in [CMS18, BCMS19] an expected complexity bound is derived for a trust region method with a stochastic zeroth-order oracle. In [CBS22] a high probability iteration complexity bound for first- and second-order stochastic trust region methods is derived under the same oracles we use. Recently, the same oracles were used within a stochastic quasi-Newton method in [MWX23], and a stochastic SQP-based method for nonlinear equality constrained problems in [BXZ23].

Adaptive regularization with cubics (ARC) methods are theoretically superior alternatives to line search and trust region methods, when applied to deterministic smooth functions, because of their optimal complexity of $O(\epsilon^{-3/2})$ vs $O(\epsilon^{-2})$ for finding $\epsilon$ stationary points [CGT11c, CGT11a]. There are many variants of adaptive cubic regularization methods under various assumptions and requirements on the function value, gradient, and Hessian estimates. Specifically, in [CGT11b, LLHT18, BGMT19, WZLL19, KL17, PJP20], the oracles are assumed to be deterministically bounded, with adaptive magnitude or errors. In [CS18, BG22, BGMT22, BGMT20], bounds on expected complexity are provided under exact or deterministically bounded zeroth-order oracle and the gradient and Hessian oracles similar to $\mathsf{SFO}$ and $\mathsf{SSO}$. A cubicly regularized method in a fully stochastic setting is analyzed in [TSJ⁺18]. The method is not adaptive, relying on the knowledge of the Lipschitz constants of $\nabla\phi(x)$, and therefore not requiring a zeroth-order oracle at all. However, the results in that paper are simply derived assuming that stochastic gradient and Hessian estimates are sufficiently accurate at each iteration. The final complexity bound only applies with probability that this holds true. No expected complexity bound can thus be derived.

The analysis presented here extends the stochastic settings and high probability results in [JSX21] and [CBS22] to the framework in [CS18]. However, this extension is far from trivial, as it requires careful modification of most of the elements of the existing analysis. We point out these modifications in the appropriate places in the paper.

The oracles used in this paper are essentially the same as in [JSX21] and [CBS22]. However, our assumption on the oracles is a bit stronger in this paper than in these two previous works. In particular, we assume that $\mathsf{SFO}$ and $\mathsf{SSO}$ are implementable for arbitrarily small values of $\mu_{1}$ and $\mu_{2}$, respectively. In contrast, the analysis in [JSX21] and [CBS22] allows for the case when these oracles cannot be implemented for arbitrarily small error bounds. We will discuss further in the paper, that even though SARC may impose small values of $\mu_{1}$ and $\mu_{2}$, this happens only with small probability.

We do not discuss the numerical performance of our method in this paper. Although deterministic implementations of ARC can be competitive with trust-region and line search methods when implemented with care, their efficiency in practice is highly dependent on the subproblem solver used. We expect similar behavior in the stochastic case and leave this study as a subject of future research.

The Stochastic Adaptive Regularization with Cubics (SARC) method is presented below as Algorithm 1. At each iteration $k$, given gradient estimate $g_{k}$, Hessian estimate $H_{k}$, and a regularization parameter $\sigma_{k}>0$, the following model is approximately minimized with respect to $s$ to obtain the trial step $s_{k}$:

The constant term $\phi(x_{k})$ is never computed and is used simply for presentation purposes. In the case of SARC, $g_{k}$ and $H_{k}$ are computed using $\mathsf{SFO}$ and $\mathsf{SSO}$ so as to satisfy certain accuracy requirements with sufficiently high probability, which will be specified in Section 3. We require the trial step $s_{k}$ to be an “approximate minimizer" of $m_{k}(x_{k}+s)$ in the sense that it has to satisfy:

and

where $\eta\in(0,1)$ is a user-chosen constant. The conditions are typical in the literature, e.g., in [CGT11c] and can be satisfied, for example, using algorithms in [CGT11b, CD19] to approximately minimize the model (4), as well as by any global minimizer of $m_{k}(x_{k}+s)$.

As in any variant of the ARC method, once $s_{k}$ is computed, the trial step is accepted (and $\sigma_{k}$ is decreased) if the estimated function value of $x_{k}^{+}=x_{k}+s_{k}$ is sufficiently smaller than that of $x_{k}$, when compared to the model value decrease. We call these iterations successful. Otherwise, the trial step is rejected (and $\sigma_{k}$ is increased). We call these iterations unsuccessful. In the case of SARC, however, function value estimates are obtained via $\mathsf{SZO}$ and the step acceptance criterion is modified by adding an "error correction" term of $2\epsilon_{f}^{\prime}$. This is because function value estimates have an irreducible error, so without this correction term, the algorithm may always reject improvement steps.

Algorithm 1 generates a stochastic process and we will analyze it in the next section. First, however, we state and prove several lemmas that establish the behavior of the algorithm for every realization.

A key concept that will be used in the analysis is the concept of a true iteration. Let $e_{k}=|f(x_{k})-\phi(x_{k})|$ and $e_{k}^{+}=|f(x_{k}^{+})-\phi(x_{k}^{+})|$.

We will now prove a sequence of results that hold for each realization of Algorithm 1, and are essential for the complexity analysis. The two key results are Corollary 1 and Lemma 5, where Corollary 1 shows that until an $\epsilon$-stationary point is reached, every true iteration with large enough $\sigma_{k}$ is successful, and Lemma 5 establishes the lower bound on function improvement on true and successful iterations. Lemmas 1, 2, 3 and 4 lay the building blocks for them: On every successful iteration, the function improvement is lower bounded in terms of the norm of the step (Lemma 1). There is a threshold value of $\sigma_{k}$ where any true iteration is either always successful or results in a very small step (Lemma 2). When an iteration is true and the step is not very small, the norm of the step is lower bounded in terms of $\|\nabla\phi(x_{k}^{+})\|$ (Lemma 3). Until an $\epsilon$-stationary point is reached, the step cannot be too small on true iterations (Lemma 4).

Note that for the above lemma to hold, $\bar{\sigma}$ does not need to include $L$ in the numerator. However, we will need another condition on $\bar{\sigma}$ later that will involve $L$; hence for simplicity of notation we introduced $\bar{\sigma}$ above to satisfy all necessary bounds.

We can now use these important properties of Algorithm 1 to show that stochastic process that the algorithm generates fits into the framework analyzed in [JSX21].

Algorithm 1 generates a stochastic process. Let $M_{k}$ denote the collection of random variables $\left\{\Xi_{k},\Xi_{k}^{+},\Xi_{k}^{1},\Xi_{k}^{2}\right\}$, whose realizations are $\left\{\xi_{k},\xi_{k}^{+},\xi_{k}^{1},\xi_{k}^{2}\right\}$. Let $\left\{{\cal F}_{k}:k\geq 0\right\}$ denote the filtration generated by $M_{0},M_{1},\ldots,M_{k}$. At iteration $k$, $X_{k}$ denotes the random iterate, $G_{k}$ is the random gradient approximation, $\mathsf{H}_{k}$ is the random Hessian approximation, $\Sigma_{k}$ is the random model regularization parameter. $S_{k}$ is the step computed for the random model, $f(X_{k},\Xi_{k})\text{ and }f(X_{k}^{+},\Xi_{k}^{+})$ are the random function estimates at the current point and the candidate point, respectively.

Conditioned on $X_{k}$ and $\Sigma_{k}$, the random quantities ${G}_{k}$ and $\mathsf{H}_{k}$ are determined by $\Xi_{k}^{1}$ and $\Xi_{k}^{2}$ respectively. The realization of $S_{k}$ depends on the realizations of ${G}_{k}$ and $\mathsf{H}_{k}$. The function estimates $f(X_{k},\Xi_{k})\text{ and }f(X_{k}^{+},\Xi_{k}^{+})$ are determined by $\Xi_{k},\Xi_{k}^{+}$, conditioned on $X_{k}$ and $X_{k}^{+}$. In summary, the stochastic process $\left\{\left({G}_{k},\mathsf{H}_{k},S_{k},f(X_{k},\Xi_{k}),f(X_{k}^{+},\Xi_{k}^{+}),X_{k},\Sigma_{k}\right)\right\}$ generated by the algorithm, with realization $\left\{\left(g_{k},H_{k},s_{k},f(x_{k},\xi_{k}),f(x_{k}^{+},\xi_{k}^{+}),x_{k},\sigma_{k}\right)\right\}$, is adapted to the filtration $\left\{{\cal F}_{k}:k\geq 0\right\}$.

We further define $E_{k}:=|f(X_{k},\Xi_{k})-\phi(X_{k})|$ and $E_{k}^{+}:=|f(X_{k}^{+},\Xi_{k}^{+})-\phi(X_{k}^{+})|$, with realizations $e_{k}$ and $e_{k}^{+}$. Let $\Theta_{k}:={\mathbbm{1}}\left\{\text{iteration $k$ is successful}\right\}$, and let $I_{k}:={\mathbbm{1}}\left\{\text{iteration $k$ is true}\right\}$. The indicator random variables $\Theta_{k}$ and $I_{k}$ are clearly measurable with respect to the filtration ${\cal F}_{k}$.

The next lemma shows that by construction of the algorithm, the stochastic model $m_{k}$ at iteration $k$ is “sufficiently accurate" with probability at least $1-\delta_{1}-\delta_{2}$.

Recall Definition 1 and that we denote the event of iteration $k$ being true by indicator random variable $I_{k}$. It is crucial for our analysis that $\mathbb{P}(I_{k}=1\mid\mathcal{F}_{k-1})\geq p>\frac{1}{2}$ for all $k$. We will later combine Lemma 6 with the properties of $\mathsf{SZO}$ for a bound on $\delta_{1}+\delta_{2}$ to ensure $p>\frac{1}{2}$.

The iteration complexity of our algorithm is defined as the following stopping time.

It is important to note that even if for some iteration $k$, $\left\lVert{\nabla\phi(X_{k}^{+})}\right\rVert\leq\epsilon$, this iteration may not be successful and thus $\left\lVert{\nabla\phi(X_{k+1})}\right\rVert$ may be greater than $\epsilon$. This is a consequence of the complexity analysis of cubic regularization methods that measure progress in terms of the gradient at the trial point and not at the current iterate, and thus is not specific to SARC. The stopping time is thus defined as the first time at which the algorithm computes a point at which the gradient of $\phi$ is less than $\epsilon$.

It is easy to see that $T_{\epsilon}$ is a stopping time of the stochastic process with respect to ${\cal F}_{k}$. Given a level of accuracy $\epsilon$, we aim to derive a bound on the iterations complexity $T_{\epsilon}$ with overwhelmingly high probability. In particular, we will show the number of iterations until the stopping time $T_{\epsilon}$ is a sub-exponential random variable, whose value (with high probability) scales as ${O}(\epsilon^{-3/2})$, similarly to the deterministic case. Towards that end, we define stochastic process $Z_{k}$ to measure the progress towards optimality.