Stochastic Cubic Regularized Newton (SCRN)

• SCRN is a stochastic optimization method that employs cubic regularization with approximate gradient and Hessian evaluations to reach near second-order stationarity in nonconvex settings.
• It forms a cubic model at each iteration by solving a subproblem that balances descent and stabilization in the presence of noisy second-order information.
• SCRN integrates mini-batch sampling, variance reduction, and momentum techniques to improve convergence rates and sample complexity over traditional methods.

Stochastic Cubic Regularized Newton (SCRN), also termed stochastic Cubic Newton or SCN in parts of the literature, is a class of stochastic second-order methods that extends cubic-regularized Newton to settings in which only noisy, subsampled, or otherwise approximate gradients and Hessians are available. The method is studied for both stochastic expectation problems $f(x)=\mathbb E_\xi[f(x;\xi)]$ and finite-sum objectives $f(x)=\frac1n\sum_{i=1}^n f_i(x)$, and in general nonconvex settings its target is typically an approximate local minimum or second-order stationary point rather than a global minimizer (Tripuraneni et al., 2017, Chayti et al., 2023).

1. Optimization model and stationarity notions

SCRN appears in two closely related optimization models. One is stochastic approximation, in which $f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$, and the other is finite-sum optimization, in which $f(x)=\frac1n\sum_{i=1}^n f_i(x)$. The first model is prominent in the original stochastic cubic regularization analysis, whereas the second underlies variance-reduced, memory-based, and helper-based variants (Tripuraneni et al., 2017, Zhou et al., 2018).

The standard target in nonconvex settings is an approximate second-order stationary point. One formulation requires

$$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$

which is the definition used in the stochastic cubic regularization analysis for smooth nonconvex objectives with $\rho$-Lipschitz Hessian (Tripuraneni et al., 2017). Finite-sum cubic-regularization papers often state the goal as an $(\epsilon,\sqrt{\epsilon})$-approximate local solution, meaning small gradient norm together with a Hessian lower bound of order $-\sqrt{\epsilon}$ (Zhang et al., 2018, Zhou et al., 2018).

Several later papers compress first- and second-order conditions into a single scalar stationarity measure. One representative definition is

$$\mu_M(x):=\max\Bigl(\|\nabla f(x)\|^{3/2},\; \frac{-\lambda_{\min}(\nabla^2 f(x))^3}{M^{3/2}}\Bigr),$$

or a closely related variant with $[-\lambda_{\min}(\nabla^2 f(x))]_+^3$. Small $f(x)=\frac1n\sum_{i=1}^n f_i(x)$0 means both small gradient and near-positive-semidefinite Hessian (Chayti et al., 2023, Chayti et al., 2024). This suggests that the cubic regularization literature treats second-order stationarity not as an auxiliary condition but as the primary termination criterion.

Across the literature, the regularity assumptions are centered on lower boundedness of $f(x)=\frac1n\sum_{i=1}^n f_i(x)$1, Lipschitz continuity of the gradient, and especially Lipschitz continuity of the Hessian. The Hessian-Lipschitz property is the condition that enables third-order Taylor control and justifies the cubic term in the local model (Tripuraneni et al., 2017, Kovalev et al., 2019).

2. Cubic model and step computation

At each iterate $f(x)=\frac1n\sum_{i=1}^n f_i(x)$2, SCRN forms a cubic model using approximate first- and second-order information. A representative update computes

$f(x)=\frac1n\sum_{i=1}^n f_i(x)$3

where $f(x)=\frac1n\sum_{i=1}^n f_i(x)$4 and $f(x)=\frac1n\sum_{i=1}^n f_i(x)$5 approximate $f(x)=\frac1n\sum_{i=1}^n f_i(x)$6 and $f(x)=\frac1n\sum_{i=1}^n f_i(x)$7, and $f(x)=\frac1n\sum_{i=1}^n f_i(x)$8 is the cubic regularization parameter (Chayti et al., 2023, Tripuraneni et al., 2017).

The cubic term is the defining difference from a stochastic Newton step. Under Hessian Lipschitzness, it stabilizes the model in directions of negative curvature and can make the model a valid upper bound. In the strongly convex finite-sum setting, the same rationale is expressed through Taylor remainder bounds controlled by the Hessian-Lipschitz constant $f(x)=\frac1n\sum_{i=1}^n f_i(x)$9, with $f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$0 ensuring a valid cubic upper model (Kovalev et al., 2019).

For an exact solution of the cubic subproblem, the standard optimality conditions are

$f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$1

These identities recur throughout the theory because they connect the step norm, local curvature, and stationarity of the next iterate (Zhang et al., 2018, Roy et al., 2020). In the interpolation-based SCRN analysis, they yield the model decrease bound

$f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$2

which is the basic cubic decrease mechanism (Roy et al., 2020).

A recurring analytical pattern is that if the iterate is not yet approximately second-order stationary, then the cubic subproblem admits a step of nontrivial norm, and the model decrease transfers to true objective decrease. Conversely, if the cubic step is sufficiently small, then the iterate is already near second-order stationarity. This “large step implies descent, small step implies stationarity” dichotomy underlies both the original stochastic analysis and later inexact, variance-reduced, and adaptive variants (Tripuraneni et al., 2017, Scheinberg et al., 2023).

3. Estimator design and algorithmic families

The simplest SCRN instantiation uses mini-batch stochastic estimators. In the original stochastic cubic regularization method, the algorithm samples independent mini-batches $f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$3 and $f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$4, forms

$f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$5

and solves the cubic model using $f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$6 and the Hessian operator $f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$7. The method is explicitly designed to rely only on stochastic gradients and stochastic Hessian-vector products, with the latter used in place of explicitly forming a dense Hessian matrix (Tripuraneni et al., 2017).

A major finite-sum direction is variance reduction. The paper "Stochastic Variance-Reduced Cubic Regularized Newton Method" introduced a semi-stochastic gradient together with a semi-stochastic Hessian specifically designed for cubic regularization, and established a complexity improvement over earlier cubic regularization schemes (Zhou et al., 2018). A more explicit stagewise design appears in adaptive SVRC, where a reference point $f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$8 is used to build estimators

$f(x)=\mathbb E_{\xi\sim\mathcal D}[f(x;\xi)]$9

$f(x)=\frac1n\sum_{i=1}^n f_i(x)$0

with adaptive sample sizes based on the previous step norm $f(x)=\frac1n\sum_{i=1}^n f_i(x)$1 rather than the unknown current step (Zhang et al., 2018).

A distinct finite-sum formulation is the memory-based stochastic cubic Newton method for strongly convex objectives. It stores points $f(x)=\frac1n\sum_{i=1}^n f_i(x)$2 for each component, builds second-order Taylor models around those stored points, and updates only a random subset $f(x)=\frac1n\sum_{i=1}^n f_i(x)$3 of the memories after each iteration. In its simplest form, each iteration computes the gradient and Hessian of a single randomly selected function only. The method is biased rather than unbiased, and its analysis is based on a tailored Lyapunov function rather than estimator unbiasedness (Kovalev et al., 2019).

The helper framework generalizes these constructions further. It represents $f(x)=\frac1n\sum_{i=1}^n f_i(x)$4 as $f(x)=\frac1n\sum_{i=1}^n f_i(x)$5, models the cheaper helper $f(x)=\frac1n\sum_{i=1}^n f_i(x)$6 locally, and handles noisy and possibly biased gradient and Hessian estimates, arbitrary batch sizes, variance reduction, and lazy Hessian updates in one analysis. This framework is presented as a unifying theory for stochastic and variance-reduced cubic Newton methods, and also applies to auxiliary learning (Chayti et al., 2023).

Recent work adds momentum to the estimator design. One line uses transported gradient momentum together with heavy-ball Hessian momentum and proves global convergence for arbitrary batch size, including one stochastic sample per iteration (Chayti et al., 2024). Another introduces Polyak momentum and recursive momentum for Hessian estimation, combined with a small-error stochastic gradient estimator and a potential function of the form

$f(x)=\frac1n\sum_{i=1}^n f_i(x)$7

to analyze stochastic second-order stationarity in moment form (Yang et al., 17 Jul 2025).

4. Complexity landscape

Reported complexity guarantees depend strongly on the oracle model: some papers count stochastic gradients and Hessian-vector products, some count higher-order oracle calls, some count Hessian samples, and some separate arithmetic complexity from oracle complexity. The principal claims therefore require direct comparison only within matching oracle models.

Setting	Reported guarantee	Paper
General smooth nonconvex SCRN with stochastic gradients and Hessian-vector products	$f(x)=\frac1n\sum_{i=1}^n f_i(x)$8 stochastic gradient and stochastic Hessian-vector product evaluations to find an $f(x)=\frac1n\sum_{i=1}^n f_i(x)$9-approximate local minimum	(Tripuraneni et al., 2017)
Finite-sum stochastic variance-reduced cubic regularized Newton	$$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$0 second-order oracle calls to an $$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$1-approximately local minimum	(Zhou et al., 2018)
SCRN under interpolation-like conditions	$$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$2 higher-order oracle calls to an $$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$3-local-minimizer	(Roy et al., 2020)
Adaptive variance-reduced subsampled Newton with cubic regularization	$$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$4 expected Hessian sample complexity	(Zhang et al., 2018)
Finite-sum EMA/SARAH cubic Newton	$$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$5 total stochastic oracle calls to an $$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$6-SOSP	(Pasechnyuk-Vilensky et al., 9 Oct 2025)

The 2018 variance-reduced method explicitly states that $$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$7 second-order oracle complexity outperforms the state-of-the-art cubic regularization algorithms including subsampled cubic regularization (Zhou et al., 2018). The interpolation-based analysis is positioned between generic stochastic and deterministic cubic Newton: without interpolation-like assumptions, the cited stochastic complexity is $$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$8; with the Strong Growth Condition it improves to $$\|\nabla f(x)\|\le \epsilon \quad\text{and}\quad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\rho\epsilon},$$9; the deterministic cubic-regularized Newton benchmark remains $\rho$0 (Roy et al., 2020).

Momentum-based analyses introduce a different stationarity regime. With stochastic second-order stationarity defined through moments,

$\rho$1

SCRN-PM attains complexity

$\rho$2

while SCRN-RM improves this to

$\rho$3

under an additional mean-cubed smoothness condition on the stochastic Hessian (Yang et al., 17 Jul 2025).

Another recent development is the proof of global convergence for nonconvex stochastic cubic Newton with arbitrary batch size. The momentum-based analysis in this line does not present the result as a single classical $\rho$4-oracle complexity bound, but it does establish convergence to second-order stationary points even when only one stochastic data sample is used per iteration (Chayti et al., 2024).

5. Structural regimes, local behavior, and improved rates

SCRN is sensitive to additional structure in the objective class. In sufficiently smooth and strongly convex finite-sum problems, the memory-based stochastic Newton and stochastic cubically regularized Newton methods achieve local linear-quadratic behavior. The SCN analysis uses the Lyapunov quantity

$\rho$5

and proves a local contraction

$\rho$6

under a suitable local initialization condition. When $\rho$7, the scheme reduces to the deterministic cubic Newton setting and the rate becomes locally superlinear or quadratic-like (Kovalev et al., 2019).

In over-parameterized or interpolation-like regimes, the Strong Growth Condition

$\rho$8

causes stochastic gradient noise to vanish near stationary points. Under this assumption, SCRN improves from the cited generic stochastic complexity $\rho$9 to $(\epsilon,\sqrt{\epsilon})$0, but still does not match the deterministic $(\epsilon,\sqrt{\epsilon})$1 rate. The same work explicitly argues that gradient-level interpolation is not sufficient to close the gap and suggests that Hessian interpolation-like assumptions may be necessary (Roy et al., 2020). This is one of the clearest articulated limitations in the SCRN literature.

Another important structural regime is gradient dominance. For functions satisfying

$(\epsilon,\sqrt{\epsilon})$2

SCRN improves the best-known sample complexity of stochastic gradient descent. The stated total sample complexity is $(\epsilon,\sqrt{\epsilon})$3 for $(\epsilon,\sqrt{\epsilon})$4 and $(\epsilon,\sqrt{\epsilon})$5 for $(\epsilon,\sqrt{\epsilon})$6 (Masiha et al., 2022). For $(\epsilon,\sqrt{\epsilon})$7, a variance-reduced SCRN with time-varying batch sizes reduces the average sample complexity to $(\epsilon,\sqrt{\epsilon})$8 for stochastic gradients and $(\epsilon,\sqrt{\epsilon})$9 for stochastic Hessians (Masiha et al., 2022).

The same paper extends this perspective to policy-based reinforcement learning under a weak gradient dominance condition of the form

$-\sqrt{\epsilon}$0

and reports that SCRN improves over stochastic policy gradient by a factor of $-\sqrt{\epsilon}$1 in the $-\sqrt{\epsilon}$2 case (Masiha et al., 2022). This suggests that cubic regularization can exploit global error-bound structure, not merely local negative-curvature information.

6. Extensions, implementations, and unresolved technical issues

The SCRN paradigm has been extended well beyond the Euclidean, single-machine baseline. On embedded Riemannian submanifolds, the Riemannian stochastic variance-reduced cubic regularized Newton method transports gradient and Hessian information between tangent spaces, solves a cubic model in the tangent space, and achieves $-\sqrt{\epsilon}$3 iteration complexity for an $-\sqrt{\epsilon}$4-second-order stationary point. With the choice $-\sqrt{\epsilon}$5, its total second-order oracle complexity is

$-\sqrt{\epsilon}$6

and the same rate is preserved under an inexact cubic subproblem solve (Zhang et al., 2020).

Distributed and nested formulations have also been studied. In Distributed Stochastic Cubic-Regularized Newton (DiSCRN), the outer model is

$-\sqrt{\epsilon}$7

with a consensus constraint $-\sqrt{\epsilon}$8. A local stopping criterion for the inner problem guarantees $-\sqrt{\epsilon}$9, and the outer objective is shown to decrease monotonically with high probability under the stated assumptions (Anderson et al., 2020).

A different extension replaces the full ambient space by random low-dimensional subspaces. The Stochastic Subspace Cubic Newton method is convex rather than nonconvex, but it is explicitly described as both a stochastic extension of cubically regularized Newton and a second-order enhancement of stochastic subspace descent. As the minibatch or subspace dimension varies, its global rate interpolates between stochastic coordinate descent and cubic regularized Newton (Hanzely et al., 2020).

High-probability adaptive regularization with cubics provides a complementary route. In this line, the method uses noisy function-value, gradient, and Hessian oracles together with an acceptance ratio

$$\mu_M(x):=\max\Bigl(\|\nabla f(x)\|^{3/2},\; \frac{-\lambda_{\min}(\nabla^2 f(x))^3}{M^{3/2}}\Bigr),$$0

and obtains a high-probability iteration bound of order $$\mu_M(x):=\max\Bigl(\|\nabla f(x)\|^{3/2},\; \frac{-\lambda_{\min}(\nabla^2 f(x))^3}{M^{3/2}}\Bigr),$$1, matching deterministic ARC/SCRN at the iteration level (Scheinberg et al., 2023).

Implementation remains a central technical issue. The original stochastic cubic regularization method emphasizes that stochastic Hessian-vector products can be computed as efficiently as stochastic gradients, which avoids explicit Hessian formation (Tripuraneni et al., 2017). The helper framework introduces lazy Hessian reuse to reduce arithmetic complexity in large dimension (Chayti et al., 2023). Re$$\mu_M(x):=\max\Bigl(\|\nabla f(x)\|^{3/2},\; \frac{-\lambda_{\min}(\nabla^2 f(x))^3}{M^{3/2}}\Bigr),$$2MCN develops a matrix-free Hutchinson variant for Hessian estimation and a fast inner solver based on a secular equation and conjugate gradients, again targeting large-scale finite-sum problems where full Hessians are impractical (Pasechnyuk-Vilensky et al., 9 Oct 2025).

A persistent unresolved issue is the gap between the best deterministic and stochastic second-order rates. Interpolation-like gradient assumptions improve SCRN substantially but do not recover deterministic cubic Newton complexity, and the literature explicitly points to missing Hessian-level structure as a plausible reason (Roy et al., 2020). The current trajectory of the field therefore combines three themes: more aggressive control of stochastic Hessian error, cheaper matrix-free or lazy implementations, and analytical frameworks that can accommodate biased, momentum-based, or auxiliary-information-driven estimators without losing second-order guarantees (Chayti et al., 2023, Chayti et al., 2024).