Linear Stochastic Approximation

• Linear Stochastic Approximation (LSA) is a class of methods that use stochastic recursions to solve linear systems or fixed-point equations from noisy observations.
• LSA employs both ODE-based techniques and Lyapunov function analysis to ensure stability and derive finite-time error bounds under constant and diminishing step-size regimes.
• LSA underpins applications in reinforcement learning, distributed optimization, and variational Bayes by unifying bias correction, averaging methods, and noise analysis.

Linear stochastic approximation (LSA) is a class of stochastic recursions for solving a linear system or tracking the stable equilibrium of a linear mean field from noisy matrix–vector observations. In its canonical forms, the iterate is written either as

$$\theta_{n+1}=\theta_n-\alpha\bigl(\mathbf A_{n+1}\theta_n-\mathbf b_{n+1}\bigr)$$

or, with a different sign convention and explicit Markovian driver,

$$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$

The subject lies at the intersection of stochastic approximation, linear systems, Markov-process theory, and reinforcement learning, and modern work treats fixed- and diminishing-step regimes, i.i.d. and Markovian noise, iterate averaging, statistical inference, and distributed variants (Lakshminarayanan et al., 2017, Srikant et al., 2019, Durmus et al., 2021).

1. Core formulation and target equations

The population object behind LSA is a linear equation of the form

$\bar A\theta^\star=\bar b,$

or, in the mean-field convention used in several papers,

$\dot\theta=\bar A\theta+\bar b.$

After centering at the equilibrium, the limiting ODE is often written simply as

$\dot\theta=\bar A\theta,$

with $\theta^\star=0$ after recentering. The standard stability hypothesis is that $-\bar A$ is Hurwitz in the root-finding convention, or equivalently that $\bar A$ is Hurwitz in the mean-drift convention, so that the deterministic linear dynamics is globally asymptotically stable (Srikant et al., 2019, Durmus et al., 2021).

This formulation covers both static noisy linear systems and linear fixed-point equations. In the i.i.d. fixed-step literature, the random observations $(\mathbf A_n,\mathbf b_n)$ are unbiased estimators of $(A,b)$, and the recursion is analyzed as a stochastic solver of $$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$0 (Durmus et al., 2021). In the Markovian setting, the coefficients are functions of a Markov chain $$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$1 or $$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$2, and the mean system is defined with respect to the invariant distribution $$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$3, for example

$$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$4

For exponential-family variational Bayes, the same algebra appears through the covariance identity

$$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$5

so the stationary condition is an expected linear system $$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$6; this is the sense in which “stochastic linear regression” is LSA-like (Salimans et al., 2014).

A useful structural distinction runs through the literature. Some analyses target the last iterate of a fixed-step algorithm; others study Polyak–Ruppert averages

$$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$7

or tail averages over the second half of the trajectory. That distinction determines whether the leading error is a steady-state fluctuation of order $$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$8 or an averaged statistical fluctuation of order $$\Theta_{k+1}=\Theta_k+\epsilon\bigl(A(X_k)\Theta_k+b(X_k)\bigr).$$9 (Durmus et al., 2021, Durmus et al., 2022).

2. Stability mechanisms and analytical frameworks

The central stability device is a quadratic Lyapunov function built from the Lyapunov equation

$\bar A\theta^\star=\bar b,$0

or, with the opposite sign convention,

$\bar A\theta^\star=\bar b,$1

With $\bar A\theta^\star=\bar b,$2, one obtains a deterministic contraction for the mean ODE and, in the stochastic recursion, a drift inequality plus bias terms. In the Markovian constant-step analysis of linear stochastic approximation, this same quadratic $\bar A\theta^\star=\bar b,$3 is interpreted both as the standard Lyapunov function of the linear ODE and as the Stein function for steady-state performance bounds; in this linear setting, the Stein equation and the Lyapunov equation coincide (Srikant et al., 2019).

For i.i.d. fixed-step LSA, a complementary framework is random matrix-product analysis. Writing

$\bar A\theta^\star=\bar b,$4

the error decomposes into a transient term carried by $\bar A\theta^\star=\bar b,$5 and a fluctuation term that is a weighted sum of innovations. Tight control of moments of $\bar A\theta^\star=\bar b,$6 yields sharp non-asymptotic bounds even when the random matrices are non-symmetric and only mean-stable, not pathwise contractive (Durmus et al., 2021).

In diminishing-step stochastic approximation with Markovian noise, the dominant tool is the ODE method. The extension of the Borkar–Meyn stability theorem to Markovian noise studies

$\bar A\theta^\star=\bar b,$7

the scaled drift $\bar A\theta^\star=\bar b,$8, and the ODE at infinity

$\bar A\theta^\star=\bar b,$9

If $\dot\theta=\bar A\theta+\bar b.$0 has $\dot\theta=\bar A\theta+\bar b.$1 as a globally asymptotically stable equilibrium and the Markov noise satisfies either strong-law-type asymptotic-rate conditions or a $\dot\theta=\bar A\theta+\bar b.$2-drift/Poisson route, then the iterates are almost surely bounded and converge to a bounded invariant set of the mean ODE (Liu et al., 2024).

Markovian analyses additionally expose mixing time as an explicit parameter. In constant-step LSA with Markovian noise, one typically defines $\dot\theta=\bar A\theta+\bar b.$3 so that conditional expectations of $\dot\theta=\bar A\theta+\bar b.$4 and $\dot\theta=\bar A\theta+\bar b.$5 are within $\dot\theta=\bar A\theta+\bar b.$6 of stationarity after $\dot\theta=\bar A\theta+\bar b.$7, with geometric mixing giving $\dot\theta=\bar A\theta+\bar b.$8 (Srikant et al., 2019). This dependence survives in finite-time error floors, sample complexity, and, in some distributed settings, interaction with graph-dependent contraction factors (Lin et al., 2021).

3. Finite-time behavior, moments, and tail phenomena

The modern finite-time theory shows that fixed-step LSA has a sharp transient-versus-stationary decomposition. For Markovian LSA with constant step size $\dot\theta=\bar A\theta+\bar b.$9, one obtains a mean-square bound of the form

$\dot\theta=\bar A\theta,$0

so under geometric mixing the steady-state scale becomes $\dot\theta=\bar A\theta,$1 (Srikant et al., 2019). In the same analysis, lower-order moments of $\dot\theta=\bar A\theta,$2 satisfy Gaussian-style bounds up to order $\dot\theta=\bar A\theta,$3, but sufficiently high steady-state moments may fail to exist. The paper states explicitly that the stationary law is “not Gaussian, not sub-Gaussian, not even sub-exponential in general” (Srikant et al., 2019).

The i.i.d. fixed-step high-probability theory sharpens this picture. Under bounded random matrices, sub-Gaussian vector noise, and only the assumption that the mean matrix is Hurwitz, the last iterate admits a decomposition into a geometrically decaying transient, a leading fluctuation term with covariance given by the discrete Lyapunov equation

$\dot\theta=\bar A\theta,$4

and higher-order multiplicative-noise corrections. The dominant steady-state fluctuation is $\dot\theta=\bar A\theta,$5, matching the fixed-step CLT scale, but the confidence dependence is only polynomial under these weak assumptions; the paper proves that Gaussian or exponential concentration cannot hold in general (Durmus et al., 2021). This rules out a common misreading of diffusion limits: weak Gaussian limits as $\dot\theta=\bar A\theta,$6 do not imply Gaussian tails at a fixed nonzero step size.

For Polyak–Ruppert-averaged fixed-step LSA, the finite-time bounds are sharper still. With i.i.d. or uniformly geometrically ergodic Markov data, the leading term in the moment and high-probability bounds for the averaged iterate matches the local asymptotic minimax covariance, and the admissible fixed step size scales only with $\dot\theta=\bar A\theta,$7, not a polynomial in $\dot\theta=\bar A\theta,$8 (Durmus et al., 2022). In the Markov case, the remainder terms depend explicitly on the mixing time $\dot\theta=\bar A\theta,$9, and the optimized step size scales as $\theta^\star=0$0 in the averaged high-probability theory (Durmus et al., 2022).

4. Averaging, asymptotics, and statistical inference

Iterate averaging is the principal device that turns fixed-step LSA from a stationary-fluctuation procedure into a statistically efficient estimator. For i.i.d. LSA with constant step size and Polyak–Ruppert averaging, the averaged iterate satisfies an MSE bound of the form

$\theta^\star=0$1

so the dominant rate is $\theta^\star=0$2 for the averaged output. The same work also shows that a constant step size cannot, in general, be chosen uniformly over broad classes of LSA problems: instancewise admissible intervals exist, but not all bounded Hurwitz classes admit a universal constant step size (Lakshminarayanan et al., 2017).

Under Markovian noise, averaging remains asymptotically optimal in a finer sense. For constant-step Markovian LSA solving a linear fixed-point problem from a trajectory of an ergodic Markov chain, the averaged iterate has a non-asymptotic instance-dependent error bound whose leading term is

$\theta^\star=0$3

matching the local asymptotic minimax limit. The paper also proves a local minimax lower bound, establishing instance-optimality of the averaged estimator in the large-sample regime and showing that the worst-case scaling is $\theta^\star=0$4 up to logarithmic factors (Mou et al., 2021).

This asymptotic viewpoint has recently been extended from estimation to inference. For Polyak–Ruppert-averaged Markovian LSA with decreasing stepsizes, one-dimensional projections satisfy a non-asymptotic Berry–Esseen bound in Kolmogorov distance of order $\theta^\star=0$5 up to logarithms, and a multiplier subsample block bootstrap yields non-asymptotically valid confidence intervals, with rate $\theta^\star=0$6 up to logarithmic factors for the bootstrap approximation and $\theta^\star=0$7 up to logs for asymptotic variance estimation (Samsonov et al., 25 May 2025). In the i.i.d. decreasing-step setting, a later refinement replaces direct comparison with the Polyak–Juditsky limit by a two-stage comparison through the finite-$\theta^\star=0$8 covariance $\theta^\star=0$9, improving the multivariate Gaussian approximation rate in convex distance to $-\bar A$0 up to logarithms and yielding a multiplier bootstrap approximation rate up to $-\bar A$1 (Butyrin et al., 14 Oct 2025).

5. Constant stepsizes, Markovian bias, and extrapolation

A central distinction between i.i.d. and Markovian fixed-step LSA is the existence of a genuine stationary bias under dependence. Viewing the pair $-\bar A$2 as a time-homogeneous Markov chain, constant-step Markovian LSA converges geometrically in Wasserstein distance to a unique invariant distribution, and the stationary mean admits the expansion

$-\bar A$3

The leading term is linear in $-\bar A$4, and the paper emphasizes that this stands in contrast with the i.i.d. case, for which the bias vanishes. It also proves that Polyak–Ruppert tail averaging reduces variance but does not affect this bias (Huo et al., 2022).

The same line of work relates the bias magnitude to mixing. In the reversible-chain setting, the coefficients $-\bar A$5 are controlled by the absolute spectral gap, so the leading bias is roughly proportional to

$-\bar A$6

which the paper interprets as proportional to the Markov-chain mixing time up to logarithmic factors (Huo et al., 2022). This identifies the mechanism behind the bias: persistent state–iterate dependence induced by temporal correlation.

Richardson–Romberg extrapolation is then a natural bias-cancellation device. With multiple constant stepsizes $-\bar A$7 and weights $-\bar A$8 satisfying

$-\bar A$9

the first $\bar A$0 terms of the bias expansion cancel (Huo et al., 2022). For the two-stepsize case used later in the literature,

$\bar A$1

the $\bar A$2 term is removed. High-order analysis of constant-step Markovian LSA with PR averaging shows that this leading bias term cannot be eliminated by averaging alone, but RR does cancel it, and the leading fluctuation term of the RR iterate still aligns with the asymptotically optimal covariance matrix of vanilla averaged LSA (Levin et al., 7 Aug 2025).

This bias perspective also changes statistical inference. For averaged constant-step Markovian LSA, a CLT holds around the stationary mean $\bar A$3, not automatically around $\bar A$4. The resulting inference pipeline therefore uses either sufficiently small $\bar A$5, zero-bias model classes, or RR extrapolation to target $\bar A$6. The same work identifies important zero-bias settings, including independent multiplicative-noise models, linear regression with independent additive observation noise, and realizable linear TD learning (Huo et al., 2023). A common misconception is therefore incorrect: in Markovian fixed-step LSA, averaging alone does not generally debias the estimator.

6. Reinforcement learning, distributed variants, and other extensions

Reinforcement learning is one of the principal application domains of LSA. Constant-step TD(0) with linear function approximation can be written exactly in LSA form by taking

$\bar A$7

after centering around the TD fixed point; TD($\bar A$8) is handled analogously by augmenting the Markov state with the eligibility trace (Srikant et al., 2019). For diminishing-step algorithms with trace variables, the Markov-noise ODE method covers GTD($\bar A$9) and ETD($(\mathbf A_n,\mathbf b_n)$0) as linear stochastic approximations on enlarged state spaces, requiring only that the averaged drift matrix be Hurwitz rather than negative definite (Liu et al., 2024). In the non-asymptotic Markovian averaged theory, the full TD($(\mathbf A_n,\mathbf b_n)$1) family for $(\mathbf A_n,\mathbf b_n)$2 appears as a direct corollary, with the instance-dependent leading term making explicit how $(\mathbf A_n,\mathbf b_n)$3 changes the covariance structure (Mou et al., 2021).

Networked and multi-agent variants produce a different extension. In distributed LSA over time-varying directed graphs with merely row-stochastic interaction matrices, the correct aggregate is a time-varying weighted average defined by the absolute probability sequence $(\mathbf A_n,\mathbf b_n)$4, and the equilibrium becomes the convex combination

$(\mathbf A_n,\mathbf b_n)$5

rather than the straight average of local equilibria. Finite-time mean-square bounds are obtained by combining a single-agent Markovian LSA analysis for the weighted average with a network-disagreement analysis based on blockwise contraction and time-varying weighted quadratic comparison functions. When equal weighting is required on directed graphs, a push-sum-type algorithm restores the arithmetic average objective (Lin et al., 2021).

Two additional directions illustrate the breadth of the LSA viewpoint. In fixed-form variational Bayes, the exponential-family gradient can be rewritten as $(\mathbf A_n,\mathbf b_n)$6, with

$(\mathbf A_n,\mathbf b_n)$7

and the resulting “stochastic linear regression” interpretation motivates covariance-aware control variates and regression coefficients that, in the ideal matching exponential-family case, yield zero-variance gradient estimators (Salimans et al., 2014). In simulation-based approximate policy iteration with linear function approximation, the inner gradient-descent recursion that solves the least-squares policy-evaluation problem is itself a linear recursion of the form

$(\mathbf A_n,\mathbf b_n)$8

so the algorithm can be read as a nested or two-timescale extension of LSA rather than as classical TD learning (Winnicki et al., 2022).

Across these developments, LSA serves less as a single algorithm than as a unifying linear stochastic dynamical template. The common analytical themes are the Hurwitz stability of the mean drift, Lyapunov or ODE comparison, explicit separation of transient and stationary components, and the recognition that averaging, concentration, and bias correction play fundamentally different roles depending on whether the driving noise is i.i.d. or Markovian.