Stochastic Approximation (SA)

• Stochastic Approximation (SA) is an iterative framework that finds roots, fixed points, or equilibria using noisy feedback and adaptive step sizes.
• SA unifies methods like Robbins–Monro recursion, stochastic dual averaging, and retraction-based updates across Euclidean, manifold, and infinite-dimensional settings.
• SA theory addresses biased updates, Markovian and heavy-tailed noise, offering convergence guarantees, finite-time rates, and statistical inference tools.

Stochastic approximation (SA) is an iterative framework for finding the root of a mean field, a fixed point, or the solution of an operator equation or variational inequality when only noisy observations are available. In canonical form, the exact drift is replaced by a stochastic field and updated with either diminishing or constant step-size. Across the literature, SA appears as Robbins–Monro recursion, projected root-finding, stochastic dual averaging, retraction-based manifold iteration, and Markov-noise-driven schemes. Contemporary theory extends far beyond unbiased Euclidean stochastic gradients: it includes biased updates, state-dependent Markov chains, heavy-tailed and long-range dependent noise, constant-stepsize steady-state analysis, infinite-dimensional Banach spaces, stochastic Nash games, reinforcement learning, and statistical inference for SA outputs (Dieuleveut et al., 2023, Karimi et al., 2019, Wang et al., 15 Feb 2026, Karandikar et al., 2024).

1. Canonical recursion and problem classes

A standard SA objective is to solve

$h(w^*)=0$

when the mean field $h$ is not directly accessible. One canonical recursion is

$w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$

with step-size $\gamma_k>0$ and stochastic field $H(w,X)$ satisfying $\mathbb E[H(w,X)]\approx h(w)$ (Dieuleveut et al., 2023). A more general biased, Markov-noise-driven form is

$$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$

where $Y_{n+1}$ evolves as a state-dependent Markov chain with kernel $P_{x_n}$, $h(x)=\int H_x(y)\,\pi_x(dy)$, and $h$0 is an additional small bias (Karimi et al., 2019).

The same pattern covers several distinct mathematical targets. In stochastic Nash games, the equilibrium can be written as the solution of a stochastic variational inequality $h$1, with

$h$2

and SA becomes a projected stochastic operator iteration (Yousefian et al., 2013). In stochastic variational inequalities solved by stochastic dual averaging, the target is again an expected-value variational inequality with unique solution $h$3 (Yan et al., 2022). In reinforcement learning and signal processing, the mean field need not be of gradient type, which is one reason the “beyond-gradient” viewpoint is central in recent surveys (Dieuleveut et al., 2023).

Constraint handling is likewise part of the canonical picture. In Euclidean spaces one often uses projection, as in

$h$4

for projected SA with Markovian noise (Haque et al., 2024). On manifolds, the Euclidean increment is replaced by a retraction or exponential map, for example

$h$5

or

$h$6

which keeps the iterates on the constraint manifold by construction (Shah, 2017, Durmus et al., 2021).

Classical diminishing-step SA typically assumes

$h$7

a condition appearing repeatedly in Euclidean, manifold, and infinite-dimensional treatments (Martin et al., 2011, Shah, 2017, Karandikar et al., 2024). Constant-step SA forms a separate regime in which the iterates define a time-homogeneous Markov chain with a stationary law (Wang et al., 15 Feb 2026, Durmus et al., 2021).

2. Stability mechanisms: Lyapunov functions, ODE limits, and Poisson equations

The dominant analytical devices for SA are Lyapunov drift, ODE approximation, martingale arguments, and, for Markovian noise, Poisson-equation decompositions. In one widely used formulation, a continuously differentiable Lyapunov function $h$8 is required to satisfy inequalities such as

$h$9

and

$w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$0

together with $w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$1-smoothness of $w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$2 (Karimi et al., 2019). In the “beyond-gradient” framework, the key condition is a Lyapunov drift of the form

$w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$3

which yields a Robbins–Siegmund-type descent inequality and supports both non-asymptotic and almost-sure analyses (Dieuleveut et al., 2023).

The ODE method interprets SA as a noisy discretization of a deterministic flow. In Euclidean and manifold settings, if the linearly interpolated iterates track the ODE

$w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$4

then the limit set is almost surely an internally chain-transitive, invariant set of the flow (Shah, 2017). In the urn-model formulation, the analogous ODE is

$w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$5

and the set of limit points is almost surely a compact connected invariant set containing no strict internal attractor (Laruelle et al., 2013). For Newton’s recursive estimator of a mixing distribution, the ODE

$w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$6

is paired with a Kullback–Leibler Lyapunov function, giving global asymptotic stability of the target mixing distribution under identifiability (Martin et al., 2011).

Markovian dependence alters the proof structure because the noise is no longer a martingale difference. A standard remedy is the Poisson equation. In the biased SA framework, one assumes that

$w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$7

admits a bounded solution with suitable Lipschitz controls, allowing the error to be decomposed into martingale parts, telescoping boundary terms, and Lipschitz perturbations (Karimi et al., 2019). The same idea underlies finite-time analysis with unbounded Markovian noise, where a Poisson solution $w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$8 is used to rewrite $w_{k+1}=w_k+\gamma_{k+1}H(w_k,X_{k+1}),$9 and to control the resulting Lyapunov recursion without geometric-mixing assumptions (Haque et al., 2024). State-dependent performative prediction uses the same device to handle bias created by a controlled Markov chain whose transition kernel depends on the current parameter (Li et al., 2021).

These techniques show that SA is not restricted to unbiased i.i.d. noise. A persistent misconception is that SA theory fundamentally requires martingale-difference noise with bounded second moments. The current literature includes biased updates, controlled Markov chains, unbounded Markovian noise, and nonclassical temporal dependence, provided the Lyapunov and Poisson-equation machinery can be verified (Karimi et al., 2019, Haque et al., 2024, Chandak et al., 20 Mar 2026).

3. Convergence regimes and finite-time rates

Under strong monotonicity and Lipschitz continuity, classical almost-sure convergence results remain central. In the stochastic Nash-game setting, if $\gamma_k>0$0 is $\gamma_k>0$1-strongly monotone and $\gamma_k>0$2-Lipschitz, and if each player’s step-size sequence satisfies the standard summability conditions together with a coordination bound on relative step-sizes, then the iterates converge almost surely to the unique solution of $\gamma_k>0$3 (Yousefian et al., 2013). The same paper derives distributed adaptive step-size recursions by minimizing a one-step mean-square error upper bound, rather than prescribing a fixed harmonic rule.

Finite-time theory has expanded substantially. For biased SA with state-dependent Markov noise, randomized stopping with

$\gamma_k>0$4

yields

$\gamma_k>0$5

when $\gamma_k>0$6 (Karimi et al., 2019). In state-dependent performative prediction, the expected squared distance to the performative stable solution satisfies

$\gamma_k>0$7

under strong convexity in $\gamma_k>0$8, sensitivity control $\gamma_k>0$9, and Poisson-equation solvability for the controlled kernel (Li et al., 2021).

For unbounded Markovian noise, a general-purpose Lyapunov theorem gives an $H(w,X)$0 mean-square rate under $H(w,X)$1 and implies $H(w,X)$2 sample complexity to reach $H(w,X)$3 (Haque et al., 2024). In two-time-scale linear SA with finite-state Markovian noise, the mean squared error of the slow sequence decreases as

$H(w,X)$4

and the leading constant exactly matches the covariance in the corresponding CLT; the paper calls these tight finite-time bounds (Haque et al., 2023).

Nonclassical noise models now have explicit finite-time rates. For i.i.d. heavy-tailed noise with only a finite $H(w,X)$5-th moment, $H(w,X)$6,

$H(w,X)$7

while for long-range dependent noise with autocovariance tail parameter $H(w,X)$8,

$H(w,X)$9

These are presented as the first finite-time moment bounds in both settings (Chandak et al., 20 Mar 2026).

Multiple-sequence SA has also reached single-timescale rates. When all involved operators are strongly monotone, single-timescale MSSA achieves $\mathbb E[H(w,X)]\approx h(w)$0; when the secondary operators are strongly monotone but the main one is only of gradient form, the rate is $\mathbb E[H(w,X)]\approx h(w)$1 (Huang et al., 2024). For smooth and strongly convex stochastic optimization, risk-dependent SA can beat the usual $\mathbb E[H(w,X)]\approx h(w)$2 expectation bound when the minimal risk $\mathbb E[H(w,X)]\approx h(w)$3 is small, attaining

$\mathbb E[H(w,X)]\approx h(w)$4

for $\mathbb E[H(w,X)]\approx h(w)$5, and even

$\mathbb E[H(w,X)]\approx h(w)$6

with an epoch-based fixed-step variant (Zhang et al., 2019).

4. Constant-step SA, fluctuation theory, and statistical inference

Constant-stepsize SA does not converge to a point in the same sense as diminishing-step SA; instead, the relevant object is the stationary law of the induced Markov chain. In Euclidean settings, if

$\mathbb E[H(w,X)]\approx h(w)$7

with $\mathbb E[H(w,X)]\approx h(w)$8 having a unique root $\mathbb E[H(w,X)]\approx h(w)$9 and Hurwitz Jacobian $$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$0, then the stationary draw $$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$1 defines the centered-scaled variable

$$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$2

For both i.i.d. and Markovian noise, the law of $$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$3 is approximated in $$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$4 by a Gaussian $$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$5 with explicit non-asymptotic error

$$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$6

for small $$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$7 (Wang et al., 15 Feb 2026). The same analysis yields non-uniform Berry–Esseen-type tail bounds:

$$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$8

For merely convex objectives with flat minima, the correct scaling may be $$x_{n+1}=x_n-\gamma_{n+1}\bigl(H_{x_n}(Y_{n+1})+B_{n+1}\bigr),$$9 rather than $Y_{n+1}$0, and the limiting law becomes a non-Gaussian Gibbs law (Wang et al., 15 Feb 2026).

A Riemannian fixed-step theory gives an analogous stationary-distribution picture. Under Lyapunov conditions, the fixed-step chain on a manifold admits a unique invariant measure $Y_{n+1}$1, is geometrically ergodic, and the family $Y_{n+1}$2 converges weakly to the Dirac mass at the solution as $Y_{n+1}$3 (Durmus et al., 2021). The same work provides an asymptotic expansion of the bias,

$Y_{n+1}$4

and a CLT for the invariant law after rescaling in the tangent space.

Statistical inference for SA outputs has become an explicit topic in stochastic variational inequalities. For stochastic dual averaging with $Y_{n+1}$5, $Y_{n+1}$6, asymptotic normality is established both for the Polyak–Ruppert average and for the last iterate (Yan et al., 2022). The ergodic covariance is

$Y_{n+1}$7

while the non-ergodic covariance is described through a Lyapunov integral after decomposition of $Y_{n+1}$8. The same paper develops online plug-in and batch-means covariance estimators and constructs asymptotic confidence ellipsoids for the true SVI solution (Yan et al., 2022).

5. Non-Euclidean and infinite-dimensional formulations

SA on Riemannian manifolds replaces vector addition by a retraction. If $Y_{n+1}$9 satisfies $P_{x_n}$0 and $P_{x_n}$1, then the iteration

$P_{x_n}$2

inherits an ODE approximation theorem analogous to the Euclidean case (Shah, 2017). For embedded submanifolds one may use approximate retractions, such as SVD truncation for fixed-rank matrices or QR orthonormalization on the Stiefel manifold, and the limiting dynamics becomes the tangent projection of the ambient field (Shah, 2017). The same framework extends to non-differentiable constraint sets through the differential inclusion

$P_{x_n}$3

with tangent and normal cones defined in the usual variational sense (Shah, 2017).

Infinite-dimensional SA requires additional geometric care because, in a general Banach space, even the law of large numbers may fail. One formulation studies

$P_{x_n}$4

in a separable Banach space $P_{x_n}$5, under the same Robbins–Monro step-size conditions and a contractive condition around the unique root $P_{x_n}$6 (Karandikar et al., 2024). For i.i.d. Gaussian $P_{x_n}$7-valued noise, almost-sure convergence follows from Fernique’s theorem, the Ito–Nisio theorem, and a deterministic contraction argument. For $P_{x_n}$8 or Skorokhod space $P_{x_n}$9, independent $h(x)=\int H_x(y)\,\pi_x(dy)$0-martingale noise with uniformly bounded second moments is sufficient. For $h(x)=\int H_x(y)\,\pi_x(dy)$1-uniformly smooth Banach spaces, a Kolmogorov three-series type condition gives the more general result, and the paper also treats the $h(x)=\int H_x(y)\,\pi_x(dy)$2 case explicitly (Karandikar et al., 2024).

These developments show that Euclidean smoothness is not intrinsic to SA. Constraint geometry can be absorbed into projections, retractions, or differential inclusions, while probabilistic control can be transferred from finite-dimensional martingale theory to geometric conditions on the ambient space (Shah, 2017, Karandikar et al., 2024).

6. Applications, algorithmic specializations, and broader scope

SA is often identified with stochastic gradient descent, but the recent literature explicitly argues for a broader “beyond-gradient” viewpoint. Compressed stochastic gradient, stochastic expectation-maximization, temporal-difference learning, online EM, policy-gradient for average reward maximization, and variance-reduced SA-SPIDER all fit the generic SA mold once an appropriate mean field and Lyapunov function are specified (Dieuleveut et al., 2023, Karimi et al., 2019). In online EM for curved exponential families, the SA state is a sufficient statistic; in policy-gradient, the mean field approximates $h(x)=\int H_x(y)\,\pi_x(dy)$3 with bias $h(x)=\int H_x(y)\,\pi_x(dy)$4, leading to

$h(x)=\int H_x(y)\,\pi_x(dy)$5

under the Markov-noise theory (Karimi et al., 2019).

Game-theoretic and networked problems have produced several specialized SA schemes. In stochastic Nash games, the distributed adaptive SA method updates each player’s decision by projected stochastic gradients and updates each player’s own step-size through a quadratic recursion depending on problem parameters; the resulting sequences converge almost surely to equilibrium under standard SA noise assumptions and minimal coordination (Yousefian et al., 2013). For estimating the Price of Stability, the aR-IP-SeG method combines randomized block-coordinate stochastic extra-(sub)gradient with iterative penalization for a merely monotone SVI constraint, obtaining $h(x)=\int H_x(y)\,\pi_x(dy)$6 iteration and sample complexity for suitably defined infeasibility and suboptimality metrics (Jalilzadeh et al., 2022).

Statistical estimation has long been an SA domain. Newton’s recursive estimator for a finite mixing distribution can be written exactly in SA form, with the Kullback–Leibler divergence

$h(x)=\int H_x(y)\,\pi_x(dy)$7

serving as a Lyapunov function for the associated ODE; this gives almost-sure consistency, and the N+P modification estimates an additional unknown parameter jointly (Martin et al., 2011). Randomized urn models can also be reformulated as SA with remainder. In the two-color nonlinear urn, concave and convex skewing functions lead to different equilibrium structures, and non-convergence results for traps imply almost-sure convergence to attracting equilibria; the same SA viewpoint is used for adaptive asset allocation in finance (Laruelle et al., 2013).

Recent optimization variants further extend the algorithmic range. Multi-cut SA for stochastic convex composite optimization constructs cutting-plane models as maxima of selected one-cut models, controls the resulting bias by

$h(x)=\int H_x(y)\,\pi_x(dy)$8

and attains near-optimal $h(x)=\int H_x(y)\,\pi_x(dy)$9 convergence up to logarithms; computational experiments report that these methods generally outperform robust SA and stochastic dual averaging across all instances tested (Liang et al., 23 May 2025). In coupled multi-sequence settings, single-timescale MSSA has been applied to bilevel optimization and communication-efficient distributed learning with performance guarantees under Lipschitz continuity and strong monotonicity assumptions, without fixed-point smoothness (Huang et al., 2024).

Taken together, these strands indicate that SA is best understood not as a single algorithm but as a theory of noisy dynamical systems for root-finding, equilibrium computation, fixed-point approximation, and stochastic inference. A plausible implication is that the unifying objects are no longer the gradient and the diminishing step-size alone, but the triad of mean field, stability certificate, and noise model—a perspective made explicit across the current literature (Dieuleveut et al., 2023, Haque et al., 2024).