Stochastic Heavy Ball (SHB) Method

• Stochastic Heavy Ball (SHB) is a momentum-based optimization method that replaces exact gradients with stochastic estimates while maintaining an inertial term.
• It unifies several momentum schemes—including Polyak’s heavy-ball and Nesterov’s methods—across smooth, nonconvex, and strongly convex settings.
• Its performance and stability depend on regime-specific factors like noise structure, step-size adaptation, and batch size, influencing convergence and generalization.

Stochastic Heavy Ball (SHB) is the stochastic counterpart of Polyak’s heavy-ball method, obtained by replacing exact gradients with stochastic gradient estimates while retaining an inertial term. In its canonical discrete form, SHB updates

$$x_{k+1}=x_k-\alpha\,G(x_k;\xi_k)+\beta(x_k-x_{k-1}),$$

or, equivalently,

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$

with step size $\alpha>0$ and momentum parameter $\beta\in[0,1)$. The method is studied across smooth and nonsmooth optimization, quadratic regression, deep learning, inverse problems, consensus, and stochastic dynamical systems, with a recurring theme: inertial averaging can alter stability, bias decay, and landscape traversal, but its benefits are strongly regime-dependent rather than universal (Yan et al., 2018, Gadat et al., 2016).

1. Algorithmic forms and unified parametrizations

In the smooth empirical-risk setting, SHB is typically posed for

$f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$

with unbiased stochastic gradients $G(x_k;\xi_k)$ satisfying $\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$. The position form

$x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$

and the velocity form

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$

are treated as exactly equivalent. A widely used unification is the stochastic unified momentum (SUM) method

$$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$

where SHB is recovered by $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$0, stochastic Nesterov’s accelerated gradient (SNAG) by $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$1, and stochastic gradient (SG) by $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$2; in that parametrization, SG appears as SGD with effective step size $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$3 (Yan et al., 2018).

A second important representation is the Iterate Moving Average (IMA) form. There SHB is rewritten through auxiliary iterates

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$4

which is equivalent to the heavy-ball recursion when

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$5

This viewpoint is central in almost-sure last-iterate analyses and in Polyak-type adaptive step-size rules for momentum, because it makes the effective scaling by $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$6 explicit (Sebbouh et al., 2020, Oikonomou et al., 2024).

More general time-varying momentum schemes subsume SHB by allowing $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$7, $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$8, and related coupling coefficients to vary with time. In one such formulation,

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$9

and classical SHB is recovered from

$\alpha>0$0

by setting $\alpha>0$1 and $\alpha>0$2. This embedding is designed to cover both SHB and SNAG under a single stochastic-approximation analysis with time-varying parameters (Vidyasagar, 13 Jun 2025).

2. Smooth convergence theory across nonconvex, convex, and strongly convex regimes

For smooth nonconvex optimization, a central result is that SHB, SNAG, and SG all attain the same order of convergence in the training objective under standard assumptions: $\alpha>0$3-smoothness, bounded gradient norm, and bounded stochastic-gradient variance. In the SUM framework, if

$\alpha>0$4

then

$\alpha>0$5

The variant-dependent term enters only through $\alpha>0$6, giving a $\alpha>0$7 contribution for SHB and a $\alpha>0$8 contribution for SNAG, but the overall conclusion is explicit: SHB and SNAG have no advantage over SG in convergence rate for nonconvex training objectives under the stated assumptions (Yan et al., 2018).

Earlier stochastic-heavy-ball analyses place the method in the broader theory of stochastic approximations to second-order dissipative dynamics. In the formulation

$\alpha>0$9

with exponential memory $\beta\in[0,1)$0 or polynomial memory $\beta\in[0,1)$1, Lyapunov arguments yield almost sure boundedness, convergence of $\beta\in[0,1)$2 to zero, and critical-point convergence of $\beta\in[0,1)$3 under smooth coercive objectives. Under additional ellipticity or sub-Gaussian conditions on the noise, convergence to local minima is established almost surely. In strongly convex settings, the same framework gives non-asymptotic rates such as

$\beta\in[0,1)$4

for exponential memory and $\beta\in[0,1)$5, together with sharper $\beta\in[0,1)$6-type behavior when $\beta\in[0,1)$7 and an explicit spectral condition is met (Gadat et al., 2016).

Almost-sure last-iterate theory was later developed directly for SHB. In the convex, smooth case, with IMA parameters

$\beta\in[0,1)$8

the last iterate satisfies

$\beta\in[0,1)$9

and $f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$0 almost surely. This yields

$f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$1

for $f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$2, and

$f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$3

in the overparametrized case with constant $f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$4. A related line of work gives last-iterate almost-sure and high-probability rates for SHB under globally convex or nonconvex objectives with $f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$5-Hölder gradients, including

$f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$6

in the nonconvex case and

$f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$7

with probability at least $f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$8 in the convex Lipschitz-gradient case (Sebbouh et al., 2020, Hudiani, 9 Jul 2025).

The time-varying unified theory extends these guarantees to biased stochastic gradients and conditional variances that may grow unboundedly with time. There the sufficient conditions are natural generalizations of Robbins–Monro and Kiefer–Wolfowitz–Blum conditions:

$f(x)=\frac{1}{n}\sum_{i=1}^n \ell(x,q_i),$9

with $G(x_k;\xi_k)$0 giving

$G(x_k;\xi_k)$1

Under KL′-type assumptions one further gets $G(x_k;\xi_k)$2 and $G(x_k;\xi_k)$3 almost surely, and under PL one obtains rates of the form

$G(x_k;\xi_k)$4

for every $G(x_k;\xi_k)$5 below an explicit threshold determined by the bias and variance decay exponents (Vidyasagar, 13 Jun 2025).

3. Acceleration, sample complexity, and the setting-dependence of momentum gains

A major strand of the SHB literature argues that, on stochastic quadratics, momentum does not generically improve sample complexity. One lower-bound analysis for quadratic optimization shows that SHB and ASG require

$G(x_k;\xi_k)$6

iterations to achieve $G(x_k;\xi_k)$7 under persistent martingale-difference noise, and that vanilla SGD matches the same order. A separate random-matrix analysis of large random least-squares problems gives an exact deterministic characterization of the loss trajectory and proves that small-batch SHB with fixed momentum parameter $G(x_k;\xi_k)$8 is equivalent, at the level of average-case loss dynamics, to SGD with effective step size $G(x_k;\xi_k)$9. In that model, fixed-momentum SHB provides no actual performance improvement over SGD when step sizes are adjusted correctly (Ganesh et al., 2021, Paquette et al., 2021).

Positive acceleration results appear only under more specific structural assumptions. For strongly convex quadratic objectives under anisotropic gradient noise,

$\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$0

and a multistage geometric step-decay schedule

$\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$1

SHB with

$\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$2

achieves last-iterate excess risk

$\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$3

yielding $\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$4 acceleration of the bias term while preserving near-minimax variance scaling. A related analysis shows that for strongly convex quadratics, SHB with the deterministic heavy-ball parameters

$\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$5

accelerates when the mini-batch size exceeds an explicit $\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$6-dependent threshold

$\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$7

Under $\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$8, one gets accelerated neighborhood convergence

$\mathbb{E}[G(x_k;\xi_k)]=\nabla f(x_k)$9

and a multistage schedule yields

$x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$0

in expectation for the $x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$1 distance to the minimizer on quadratics, alongside

$x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$2

for general smooth strongly convex objectives in squared norm (Pan et al., 2023, Dang et al., 2024).

Taken together, these results suggest that claims about SHB “acceleration” are inseparable from the noise model, batch regime, stepsize schedule, and objective class. In the stochastic nonconvex deep-learning analysis of SUM, acceleration over SG is absent for training convergence (Yan et al., 2018); in quadratic sample-complexity analyses, fixed-momentum gains disappear after effective-step normalization (Ganesh et al., 2021, Paquette et al., 2021); in contrast, anisotropic noise, large batches, or carefully designed stagewise schedules can recover accelerated bias decay (Pan et al., 2023, Dang et al., 2024).

4. Stability, generalization, and landscape properties

Although SHB does not improve nonconvex training rates over SG in the unified deep-learning analysis, it does improve uniform stability. For two runs on neighboring datasets, the stability recursion depends on cumulative gradient weights

$x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$3

Plugging in the special cases gives

$x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$4

hence

$x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$5

so SHB is the most stable of the three and SG the least stable. Under identical iteration counts, adding momentum in SHB or SNAG reduces the uniform-stability coefficient and can improve generalization even though the training convergence rate remains $x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$6. In the reported CIFAR-10 experiments, the best initial step size was $x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$7 for SHB and $x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$8 for SNAG and SG; with the same initial $x_{k+1}=x_k-\alpha G(x_k;\xi_k)+\beta(x_k-x_{k-1})$9, training-error curves were very close across methods, but the absolute train–test gap was smaller for SHB and SNAG than for SG (Yan et al., 2018).

A mean-field treatment on wide two- and three-layer neural networks adds a distinct notion of stability. In that setting SHB is modeled as a second-order mean-field dynamics, and the finite-width trajectory is compared to a kinetic mean-field limit with dimension-free non-asymptotic bounds. This yields dropout-stability bounds of order

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$0

for two-layer networks, with analogous three-layer bounds, and also establishes low-loss connectivity between SHB solutions. Under universal approximation, full-support initialization, and a stated convergence mode, the three-layer mean-field heavy-ball dynamics converges to the global minimum of the unregularized population risk (Wu et al., 2022).

These analyses separate two frequently conflated questions. One concerns optimization speed, where SHB may or may not improve rates depending on regime. The other concerns the geometry and robustness of the attained solution, where the literature records smaller train–test gaps, stronger uniform stability, dropout-stability, and low-loss connectivity under SHB dynamics (Yan et al., 2018, Wu et al., 2022).

5. Nonsmooth analysis and continuous-time stochastic dynamics

SHB has also been analyzed beyond smooth Euclidean gradients. For nonsmooth nonconvex objectives

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$1

one formulation uses

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$2

or equivalently

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$3

under diminishing stepsizes, bounded iterates, and an exponential-memory regime $$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$4. The analysis is built on semialgebraic or definable assumptions, conservative gradients, and the differential inclusion

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$5

with Lyapunov function

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$6

Essential accumulation points and minimal-energy accumulation points satisfy

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$7

almost surely, and under definable distributional assumptions every accumulation point has this property while $$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$8 converges. The same framework identifies “artificial critical points” arising from conservative calculus and proves that randomized initializations avoid them almost surely, so convergence is to Clarke critical points on a full-measure set of initial conditions (Le, 2023).

Continuous-time stochastic SHB has been studied through diffusion approximations of small-step heavy-ball recursions. In one formulation,

$$v_{k+1}=\beta v_k-\alpha G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1}$$9

and after rescaling time,

$$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$0

For Morse functions with the strict saddle and strong saddle properties, the expected exit time from a saddle neighborhood scales as

$$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$1

where

$$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$2

and a chain of $$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$3 saddles yields

$$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$4

This suggests fast saddle escape in the diffusion model (Hu et al., 2017).

A separate Hilbert-space theory studies stochastic Heavy Ball systems directly at the trajectory level. For convex minimization with Lipschitz gradient,

$$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$5

square-integrable diffusion and

$$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$6

imply almost sure weak convergence of $$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$7 to a minimizer, $$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$8, and

$$y_{k+1}=x_k-\alpha G(x_k),\qquad y^s_{k+1}=x_k-s\alpha G(x_k),\qquad x_{k+1}=y_{k+1}+\beta(y^s_{k+1}-y^s_k),$$9

With time rescaling, this system becomes a stochastic Su–Boyd–Candès-type vanishing-damping dynamics and yields

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$00

for $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$01. The same paper extends the analysis to monotone equations via a corrected SHB system and shows equivalence to a stochastic Fast OGDA-type second-order dynamics with residual and gap rates

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$02

(Bot et al., 3 Oct 2025).

6. Linear algebra, PCA, and distributed consensus

In linear systems, SHB admits a sketch-and-project interpretation. For the stochastic reformulation of a consistent system $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$03,

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$04

the adaptive SHBM update

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$05

chooses $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$06 and $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$07 by projecting the minimum-norm solution onto the affine span generated by the current stochastic gradient and the previous step. This yields linear convergence in expectation with a contraction factor strictly tighter than the corresponding basic stochastic Polyak method whenever the geometric gain $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$08. In the deterministic case, the same recursion becomes a conjugate-gradient-type method; in fact, for $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$09 it reduces to CGNE, and the stochastic generalization induces a stochastic conjugate-gradient framework (Zeng et al., 2023).

For PCA, the heavy-ball idea appears in variance-reduced power iteration. Given covariance

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$10

the variance-reduced heavy-ball power update uses

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$11

and

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$12

with normalization after each step. The resulting method enjoys global linear convergence in expectation to the leading eigenvector, with a per-epoch contraction factor that splits into a deterministic accelerated term and a stochastic variance term. The analysis emphasizes that large mini-batches reduce the covariance variance parameter $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$13, allowing larger $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$14 and making acceleration especially visible when the eigengap is small (Kim et al., 2019).

SHB also acts as a randomized gossip method for average consensus. If $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$15 is an incidence matrix encoding the equalities $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$16 on graph edges, then the sketch-and-project SHB update

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$17

becomes a distributed momentum gossip protocol. In the pairwise case, when edge $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$18 is selected,

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$19

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$20

while nonselected nodes perform only the momentum update. In the block case, every node in a sampled connected component moves toward the component average plus the same heavy-ball term. The protocol preserves the network average when $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$21, and the expected error contracts linearly under explicit spectral conditions on the mean sketch matrix $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$22 (Loizou et al., 2018).

7. Inverse problems, adaptive step-sizes, and regularization viewpoints

For linear ill-posed inverse problems, one stochastic heavy-ball scheme uses

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$23

with

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$24

Its iterate-moving-average representation

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$25

enables regularization analysis in expectation. If $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$26 and the stopping index satisfies $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$27 with $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$28, then

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$29

Under the benchmark source condition

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$30

one obtains

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$31

and the a priori choice $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$32 gives

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$33

(Jin et al., 2024).

A more general adaptive SHB for nonlinear ill-posed systems incorporates a strongly convex penalty $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$34, dual variables $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$35, an adaptive step size, an adaptive momentum coefficient, and a discrepancy-inspired stopping rule. The update reads

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$36

with

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$37

and $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$38 chosen from a residual-based minimum formula involving $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$39. The momentum parameter is computed by minimizing a quadratic upper bound,

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$40

subject to stability thresholds. The resulting Bregman-descent inequality

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$41

yields almost sure finite termination of the discrepancy-style stopping rule, almost sure convergence for exact data, and regularization in both almost-sure and expectation senses for noisy data (Gu et al., 13 May 2026).

Adaptive step-size design has also been transplanted from stochastic Polyak step-sizes to SHB. Using the IMA viewpoint and the identity $$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$42, three momentum-aware rules are proposed:

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$43

for MomSPS$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$44,

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$45

for MomDecSPS, and

$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$46

for MomAdaSPS. MomSPS$$v_{k+1}=\beta v_k-\alpha\,G(x_k;\xi_k), \qquad x_{k+1}=x_k+v_{k+1},$$47 converges to a neighborhood for convex smooth problems without interpolation, and to the true solution at a fast rate matching deterministic HB under interpolation. MomDecSPS and MomAdaSPS are described as the first adaptive step-sizes for SHB that guarantee convergence to the exact minimizer without a priori problem parameters and without interpolation (Oikonomou et al., 2024).

Across these application domains, SHB functions not only as a momentum heuristic but as a regularization mechanism, a projection-based accelerator, and a continuous-time inertial model. The common algebraic skeleton is simple, but the provable behavior depends sharply on how stochasticity enters, how parameters are scheduled, and whether the ambient problem is a smooth objective, a nonsmooth definable landscape, a spectral estimation task, or an inverse problem with explicit regularization (Jin et al., 2024, Gu et al., 13 May 2026, Oikonomou et al., 2024).