Cubic-Regularized Newton Method

• Cubic-Regularized Newton (CRN) is a second-order optimization method that augments the Newton step with a cubic penalty to control the third-order Taylor remainder under Hessian Lipschitz continuity.
• It achieves second-order stationarity in nonconvex problems by escaping strict saddle points, with global guarantees ranging from O(1/k²) rates in convex settings to O(ε⁻³/²) for nonconvex objectives.
• CRN’s framework extends to various domains including manifold optimization, stochastic settings, and matrix problems through techniques like inexact subproblem solvers and Krylov subspace methods.

Cubic-Regularized Newton (CRN), often called cubic regularization (CR), is a second-order optimization method that replaces the raw Newton step by the global minimizer of a local quadratic model augmented with a cubic penalty. In its standard unconstrained form, for

$\min_{x\in\mathbb{R}^d} f(x),$

the method computes

$$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$

The cubic term is used because it stabilizes Newton steps when the Hessian is indefinite, noisy, or only approximately available, and under Hessian Lipschitz continuity it yields a model that controls the third-order Taylor remainder. For nonconvex problems, this makes CRN a canonical route to second-order stationary points rather than merely first-order critical points (Wang et al., 2018).

1. Core model and interpretation

The basic CRN construction starts from a twice continuously differentiable objective that is bounded below. The local model at $x_k$ is the second-order Taylor approximation plus a cubic regularizer, with $M>0$ acting as the regularization parameter. A central analytic fact is the Hessian-Lipschitz Taylor bound

$$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$

which explains why a cubic term with sufficiently large coefficient yields a global upper model of the objective. In this sense, the cubic penalty is not an ad hoc damping device: it is matched to the third-order remainder controlled by Hessian smoothness (Wang et al., 2018).

The method is designed to target second-order stationarity. In the exact scalar unconstrained setting, a second-order stationary point satisfies

$\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$

while the approximate notion used throughout the CRN literature typically requires

$$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$

This stronger criterion is the reason CRN is repeatedly studied for nonconvex optimization: it is explicitly intended to escape strict saddle points rather than merely reduce gradient norm (Wang et al., 2018).

The same modeling principle reappears in several normalizations and application domains. In convex optimization, the full-dimensional cubic model underlies the classical global rate $f(x_k)-f^\star = O(1/k^2)$, while in nonconvex optimization it supports the standard second-order stationarity guarantees. This breadth suggests that CRN is best viewed as a framework for curvature-aware model minimization, rather than as a single narrowly defined algorithm (Jiang et al., 2024).

2. Optimality conditions and global guarantees

At the minimizer of the cubic subproblem, the model satisfies a set of identities that are structurally central to CRN analysis. If

$$s=\arg\min_{u\in\mathbb{R}^d}\left\{\mathbf g^\top u+\frac12 u^\top \mathbf H u+\frac{M}{6}\|u\|^3\right\},$$

then

$\mathbf g+\mathbf H s+\frac{M}{2}\|s\|s=0,$

$$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$0

and

$$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$1

The first identity is the first-order optimality condition for the cubic model, the second is the second-order necessary condition, and the third gives a built-in model decrease proportional to $$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$2 (Wang et al., 2018).

These identities translate into the standard worst-case complexity bounds. In Euclidean nonconvex optimization, exact CRN and practical inexact CRN variants preserve the order-level rate $$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$3 for reaching second-order stationarity. One current-iterate formulation proves that after $$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$4 iterations the sequence contains a point $$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$5 such that

$$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$6

which is exactly the classical cubic-regularization scaling (Wang et al., 2018).

The same order is retained beyond Euclidean spaces. On Riemannian manifolds, a pullback-based CRN reaches a second-order $$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$7-stationary point in $$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$8 iterations under locally Lipschitz pullback Hessians, with local superlinear or quadratic behavior under stronger assumptions (Zhang et al., 2018). In matrix vector spaces, CRN applied to smooth unconstrained matrix objectives also attains $$s_{k+1}=\arg\min_{s\in\mathbb{R}^d}\left\{\nabla f(x_k)^\top s+\frac12 s^\top \nabla^2 f(x_k)s+\frac{M}{6}\|s\|^3\right\},\qquad x_{k+1}=x_k+s_{k+1}.$$9 global complexity and local quadratic convergence near non-degenerate minimizers (Garner et al., 2022).

For convex objectives, the guarantees sharpen. Full-dimensional CRN achieves $x_k$0 in the convex case, and subspace variants can recover that rate under favorable spectral conditions. Under strong convexity, several CRN analyses obtain linear convergence, while local quadratic convergence survives in nondegenerate regimes (Jiang et al., 2024).

3. Exact, inexact, and matrix-free solution of the cubic subproblem

Although the outer CRN framework is conceptually simple, the cubic subproblem is itself nontrivial. The canonical model

$x_k$1

may have multiple saddle points, a local maximizer, and poor local minima. A detailed subproblem analysis nevertheless shows that gradient descent can approximate the global minimum efficiently, with at most logarithmic dependence on the dimension after a small random perturbation removes the hard case. This result is important because it shows that CRN need not rely on dense exact inner solves; matrix-free first-order dynamics can be a principled inner solver for the cubic model (Carmon et al., 2016).

A separate issue is implementable inexactness. Earlier inexact CR analyses used a Hessian condition of the form

$x_k$2

which depends on the future step and is therefore not directly implementable. A later note replaces this by the current-iterate adaptive conditions

$x_k$3

and proves the same order-level convergence rate. The proof works by controlling total decrease over many iterations rather than per-iteration decrease, which makes the adaptive rule operational (Wang et al., 2018).

Several implementation lines push CRN farther from exact Hessian algebra. Hessian-free and zeroth-order variants use finite differences to approximate derivatives, combine this with an adaptive search over the regularization parameter and finite-difference step, and reuse a previously computed Hessian approximation for several iterations through lazy Hessian updates. In that setting, the reported global complexities are $x_k$4 function-and-gradient evaluations for the Hessian-free method and $x_k$5 function evaluations for the derivative-free method (Doikov et al., 2023).

At the quasi-Newton end of the spectrum, ARCLQN solves the standard cubic model exactly with a limited-memory quasi-Newton matrix $x_k$6 rather than the true Hessian. Using the compact form

$x_k$7

the subproblem is reduced to a scalar equation for $x_k$8, and the solver uses only norms inside the inner Newton loop. The resulting exact matrix-free subproblem solver exploits limited-memory structure rather than explicit Hessian construction (Forristal et al., 2022).

4. Geometry, constraints, and non-Euclidean extensions

CRN extends naturally to optimization on manifolds. For a Riemannian manifold $x_k$9 with retraction $M>0$0, the pullback

$M>0$1

turns the manifold problem into an unconstrained tangent-space problem, and CRN minimizes the cubic model

$M>0$2

Under compactness and local pullback-Hessian Lipschitz continuity, this retains the Euclidean $M>0$3 complexity; with a nondegenerate Riemannian Hessian, it becomes locally quadratic (Zhang et al., 2018).

A different extension arises in spectral-constrained matrix optimization. There, matrix functions such as $M>0$4, $M>0$5, and scalar spectral maps are used to rewrite constraints on eigenvalues as smooth unconstrained problems over symmetric matrices. The required derivative machinery is supplied by a second-order chain rule for compositions $M>0$6, after which the matrix CRN subproblem becomes

$M>0$7

This allows CRN to handle spectrally constrained objectives and, in the paper’s application, a fairness-oriented robust covariance estimation model related to Tyler’s M-estimator (Garner et al., 2022).

CRN has also been generalized to saddle-point models. For strongly-convex-strongly-concave $M>0$8, a cubic-regularized local saddle model adds a positive cubic term in the primal variable and a negative cubic term in the dual variable: $M>0$9 With a merit function $$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$0, this yields global linear and local quadratic convergence; for merely convex-concave problems, the same work introduces a homotopy continuation scheme (Huang et al., 2020).

Vector optimization introduces a further modification. For cone-ordered multiobjective problems, the cubic-regularized Newton direction is defined by minimizing a max-scalarized quadratic model plus $$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$1. In that setting the paper shows that the full step can be $$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$2-descent when $$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$3, even though the direction itself is not necessarily a descent direction. The resulting method does not use line search, has a global $$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$4 rate of convergence, and retains local q-quadratic convergence (Ghosh, 17 May 2025).

5. Scalable, stochastic, and distributed CRN

Large-scale and finite-sum settings have motivated stochastic and variance-reduced versions of CRN. A stochastic variance-reduced cubic-regularized Newton method uses a semi-stochastic gradient and a semi-stochastic Hessian tailored to cubic regularization, and proves convergence to an $$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$5-approximately local minimum within $$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$6 second-order oracle calls (Zhou et al., 2018). On embedded Riemannian manifolds, the analogous R-SVRC algorithm combines full gradient/Hessian refreshes at epoch boundaries with variance-reduced inner updates, obtains $$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$7 iteration complexity, and achieves second-order-oracle complexity

$$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$8

with exact or inexact cubic subproblem solutions (Zhang et al., 2020).

Subspace restriction is another major scaling device. Krylov CRN computes the cubic step inside the Krylov subspace

$$\left|f(y)-f(x)-\nabla f(x)^\top (y-x)-\frac12 (y-x)^\top \nabla^2 f(x)(y-x)\right| \le \frac{L}{6}\|y-x\|^3,$$9

using Lanczos to obtain an $\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$0-dimensional reduced problem whose dominant cost is $\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$1 Hessian-vector products. For convex optimization it proves the dimension-independent rate

$\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$2

and under favorable spectral structure it recovers the full $\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$3 CRN behavior (Jiang et al., 2024). In nonconvex coordinate-subspace form, SSCN applies cubic regularization on random coordinate blocks, interpolates between coordinate-descent-like and full-CRN behavior, and with adaptive sampling achieves the exact convergence rate $\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$4 to a second-order stationary point without sampling all coordinates (Zhao et al., 2024).

Momentum has been introduced in both deterministic and stochastic forms. CRm adds a momentum extrapolation and a monotone selection step to classical CRN, preserves the optimal $\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$5 worst-case rate to an $\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$6-second-order stationary point, and under a local error bound condition converges quadratically (Wang et al., 2018). In the stochastic setting, momentum-based Hessian estimators lead to Polyak-momentum and recursive-momentum SCRN methods with complexities

$\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$7

respectively, for the stochastic second-order stationarity notion used in that work (Yang et al., 17 Jul 2025).

CRN has also entered domain-specific large-scale settings. In federated learning, DP-FCRN combines local stochastic cubic models, Gaussian perturbation inside the local solver, and random-$\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$8 sparsification on the uplink; the paper’s qualitative conclusion is that sparsification reduces the needed noise by roughly a factor $\nabla f(x)=0,\qquad \nabla^2 f(x)\succeq 0,$9 and that second-order CRN updates outperform first-order DP Fed-SGD in both accuracy and convergence speed (Huo et al., 2024). In reinforcement learning, VR-CR-PN introduces Hessian-aided variance reduction for policy optimization and proves $$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$0 sample complexity to an $$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$1-second-order stationary point without importance sampling (Sun et al., 14 Jul 2025).

6. Asymptotic geometry, local acceleration, and comparative interpretation

Worst-case $$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$2 complexity is only the baseline behavior of CRN. Under the Kurdyka–Łojasiewicz (KL) property with exponent $$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$3, the asymptotic rate can be fully characterized. For the second-order stationarity measure

$$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$4

the regimes are finite termination when $$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$5, super-linear convergence when $$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$6, linear convergence when $$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$7, and sub-linear convergence when $$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$8. The same paper proves finite trajectory length and corresponding rates for the function gap, iterate distance, and distance to the second-order stationary set (Zhou et al., 2018).

A complementary local theory replaces nondegeneracy by a local error bound. Under

$$\|\nabla f(x)\|\le \epsilon,\qquad \lambda_{\min}(\nabla^2 f(x))\ge -\sqrt{\epsilon}.$$9

the CR iterates converge at least Q-quadratically to a second-order critical point, even when the solution set is degenerate or non-isolated. The same work shows that, without assuming convexity, this error-bound condition is equivalent to a quadratic growth condition under mild additional assumptions (Yue et al., 2018). This is one reason CRN remains relevant in nonconvex models with symmetry-induced degeneracy, such as phase retrieval and low-rank factorization.

For uniformly convex composite objectives whose smooth part has Hölder continuous Hessian, adaptive CRN admits a global linear rate governed by a degree-dependent second-order condition number

$f(x_k)-f^\star = O(1/k^2)$0

That analysis argues that CRN automatically achieves the best possible global complexity bound among the corresponding uniformly convex problem classes and clarifies why Newton-type conditioning can be substantially more favorable than gradient-method conditioning, especially when large quadratic components are present (Doikov et al., 2019).

Several common simplifications are therefore misleading. CRN is not merely Newton’s method with an added penalty term; across the literature, the cubic regularizer is the mechanism that makes the model globally well-behaved under Hessian Lipschitz continuity. CRN does not require exact dense Hessians in every implementation; implementable inexact conditions, finite-difference models, Krylov restriction, limited-memory quasi-Newton structure, variance reduction, and Hessian-vector-product inner solvers all preserve substantial parts of the theory. Nor is CRN confined to unconstrained Euclidean minimization: manifold optimization, spectral matrix problems, saddle-point systems, vector optimization, federated learning, and policy optimization all admit genuine CRN formulations. Taken together, these developments suggest that CRN is best understood as a general cubic-model methodology for second-order optimization rather than as a single fixed algorithmic template.