Markovian Linear Stochastic Approximation

• Markovian linear stochastic approximation is a framework for analyzing linear recursions whose noise is generated by a dependent Markov chain, essential for applications like reinforcement learning.
• It employs tools such as Lyapunov and Poisson equations to derive stability conditions, finite-time error bounds, and to quantify stationary bias under constant stepsizes.
• The methodology supports algorithm design in areas including temporal-difference learning, distributed optimization, and two-timescale methods, ensuring robust convergence and instance-optimal performance.

Markovian linear stochastic approximation (MLSA) studies linear recursive schemes whose driving data are generated by a Markov chain rather than by independent sampling. A standard formulation is

$$\theta_{n+1}=\theta_n+\alpha_{n+1}\{-A(Z_{n+1})\theta_n+b(Z_{n+1})\},$$

where $(Z_n)$ is a Markov chain, $\alpha_n$ is a gain sequence, and the averaged coefficients under the stationary law $\pi$ are

$A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$

The canonical target is the solution $\theta^\star$ of $A\theta^\star=b$, or equivalently the stable equilibrium of the associated mean ODE when $-A$ is Hurwitz. The subject lies at the intersection of stochastic approximation, Markov-chain theory, and reinforcement learning, with core themes including stability on general state spaces, finite-time error bounds, stationary bias under constant stepsize, averaged-iterate optimality, and temporal-difference learning with linear function approximation (Durmus et al., 2021).

1. Formal model and canonical representations

The linear recursion appears in several equivalent notational forms. Constant-gain analyses often use

$$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$

with centered mean dynamics described by the ODE

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

after translating the equilibrium to the origin and assuming $(Z_n)$0 is Hurwitz (Srikant et al., 2019). In fixed-point formulations, one seeks $(Z_n)$1 satisfying $(Z_n)$2, and then studies the streaming recursion

$(Z_n)$3

along a single Markov trajectory (Mou et al., 2021).

A central structural object in MLSA is the random matrix product

$(Z_n)$4

which governs how initial error and future perturbations propagate. If $(Z_n)$5, then

$(Z_n)$6

with

$(Z_n)$7

This decomposition reduces finite-time analysis to control of a multiplicative term $(Z_n)$8 and an additive fluctuation term (Durmus et al., 2021).

Under Markovian sampling, the iterate process alone is generally not Markov. The natural time-homogeneous Markov object is the joint process $(Z_n)$9 or $\alpha_n$0, depending on the model. That viewpoint is fundamental in convergence-to-stationarity results, in bias analysis, and in distributed or multiscale extensions (Huo et al., 2022).

2. Stability theory and mean-field structure

The primary deterministic stability condition is Hurwitzness of the averaged linear drift. In the general-state-space setting, the assumption is

$\alpha_n$1

equivalently, there exists $\alpha_n$2 solving

$\alpha_n$3

In bounded finite-state settings, the same role is played by a matrix $\alpha_n$4 satisfying

$\alpha_n$5

These Lyapunov equations define weighted norms in which the mean dynamics contract (Durmus et al., 2021, Srikant et al., 2019).

A major development in MLSA is the relaxation of boundedness and compactness assumptions. One line of work assumes that the underlying Markov chain is irreducible and aperiodic and satisfies a super-Lyapunov drift condition, stronger than the usual Foster–Lyapunov drift, which implies $\alpha_n$6-uniform geometric ergodicity: $\alpha_n$7 for some $\alpha_n$8 and $\alpha_n$9. Under this framework, the entries of $\pi$0 may be unbounded, provided their growth is controlled through a polylogarithmic envelope involving $\pi$1. This is precisely the mechanism that permits general-state-space MLSA with unbounded features or rewards (Durmus et al., 2021).

A complementary stability framework extends the Borkar–Meyn theorem from martingale-difference noise to Markovian noise. For the general recursion

$\pi$2

with linear specialization $\pi$3, the averaged drift is

$\pi$4

and the scaled limiting drift is

$\pi$5

If $\pi$6 is a globally asymptotically stable equilibrium of the ODE $\pi$7, and if either strong-law-type averaging or a Poisson-equation/$\pi$8-drift framework holds for the Markov chain, then $\pi$9 almost surely. This places Markovian SA stability on essentially the same ODE-at-infinity footing as classical i.i.d. SA (Liu et al., 2024).

When intrinsic stability is difficult to verify, expanding projections provide an additional stabilization device. In that scheme, the update is projected back into an increasing family of feasible sets $A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$0, while the Markovian noise is controlled through a Poisson-equation decomposition. The approach also permits random stepsizes, which is useful when the family of Markov kernels is not smooth in the parameter (Andrieu et al., 2011).

3. Finite-time bounds, moment growth, and random matrix products

Finite-time MLSA theory has two recurrent forms: transient contraction inherited from the mean ODE, and a residual fluctuation term determined by gain size and mixing. In the constant-stepsize bounded setting with mixing time $A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$1, one obtains a mean-square bound of the form

$A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$2

More precisely, after a burn-in of order $A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$3, the drift of the quadratic Lyapunov function yields a contraction rate roughly $A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$4, and the recursion reaches its asymptotic neighborhood in

$A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$5

samples. The same analysis shows that lower-order moments up to order $A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$6 are controlled by Gaussian-like bounds, whereas sufficiently high moments may be infinite in steady state. This rules out a blanket sub-Gaussian interpretation of constant-gain MLSA, even when the mean ODE is stable (Srikant et al., 2019).

For general state spaces and possibly unbounded coefficients, the random-matrix-product viewpoint yields sharper nonasymptotic control. Under super-Lyapunov drift, controlled-growth conditions on $A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$7, and Hurwitzness of the averaged matrix, the $A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$8-th moment of the multiplicative term satisfies

$A=E_\pi[A(Z_0)],\qquad b=E_\pi[b(Z_0)].$9

where $\theta^\star$0. The resulting error bound for the recursion has a transient term and a fluctuation term: $\theta^\star$1 Moreover, the fluctuation admits the decomposition

$\theta^\star$2

where $\theta^\star$3 is the dominant $\theta^\star$4-order term, while $\theta^\star$5 is smaller. For $\theta^\star$6,

$\theta^\star$7

This separates the leading additive-functional fluctuation from the higher-order multiplicative remainder (Durmus et al., 2021).

With classical diminishing gains, explicit MSE asymptotics are also available from a Poisson-equation analysis. For

$\theta^\star$8

if every eigenvalue $\theta^\star$9 of $A\theta^\star=b$0 satisfies $A\theta^\star=b$1, equivalently $A\theta^\star=b$2 is Hurwitz, then

$A\theta^\star=b$3

so $A\theta^\star=b$4. If instead the least stable mode satisfies $A\theta^\star=b$5, then the decay rate can deteriorate to $A\theta^\star=b$6. This provides an explicit sense in which weak mean stability slows MLSA even under geometric Markov ergodicity (Chen et al., 2020).

4. Constant-stepsize stationary behavior, bias, and memory effects

A defining feature of constant-gain MLSA under Markovian sampling is that the stationary mean is generally biased. In the joint-chain formulation $A\theta^\star=b$7, one can prove convergence in Wasserstein distance to a unique stationary distribution under uniform ergodicity of the Markov chain, boundedness of $A\theta^\star=b$8 and $A\theta^\star=b$9, Hurwitz stability of the mean matrix, and a small-gain condition. The stationary bias then admits a full stepsize expansion: $-A$0 and, for sufficiently small $-A$1,

$-A$2

In particular,

$-A$3

This contrasts sharply with the i.i.d. case, where the stationary bias vanishes. In the reversible setting, the magnitude of the leading bias term is controlled by the absolute spectral gap, yielding the qualitative relation “faster mixing, smaller bias” (Huo et al., 2022).

The source of the bias is temporal dependence coupled with multiplicative noise. In a general constant-stepsize Markovian SA framework, the stationary mean satisfies

$-A$4

where $-A$5 is a stationary average of a disturbance-correction term $-A$6. For linear recursions, the bias representation simplifies, but it still need not vanish when $-A$7 depends on the Markov state. The same analysis shows that the joint parameter–disturbance process is geometrically ergodic in a topological sense, and that the asymptotic covariance of Polyak–Ruppert averages is

$-A$8

with the $-A$9 correction potentially amplified by ill-conditioning of the mean dynamics (Lauand et al., 2023).

This stationary-bias picture clarifies several misconceptions. Stable mean dynamics do not imply zero bias, and Polyak–Ruppert averaging does not remove the Markovian bias itself. Tail averaging reduces variance from order $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$0 to order $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$1, but its first-moment expansion retains the same $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$2 stationary component. By contrast, Richardson–Romberg extrapolation with $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$3 stepsizes cancels the first $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$4 terms in the bias expansion; in the two-stepsize case it reduces the leading bias from $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$5 to $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$6 (Huo et al., 2022).

In nonlinear constant-stepsize SA with Markovian data, the asymptotic bias decomposes into $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$7, where $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$8 is the Markovian term, $$\Theta_{k+1}=\Theta_k+\epsilon\big(A(X_k)\Theta_k+b(X_k)\big),$$9 is the nonlinearity term, and $\dot\theta=\bar A\theta,$0 is a compound interaction term. Specializing this framework to the linear case $\dot\theta=\bar A\theta,$1 yields $\dot\theta=\bar A\theta,$2, so $\dot\theta=\bar A\theta,$3 and only the pure Markovian contribution remains. This isolates linear MLSA as the regime in which memory effects appear without nonlinear interaction terms (Huo et al., 2024).

5. Averaging, optimality, diffusion limits, and statistical inference

Averaged MLSA under Markovian noise admits both sharp asymptotic characterizations and genuinely nonasymptotic inference guarantees. For a linear fixed-point equation observed along a single ergodic Markov trajectory, a constant-stepsize averaged SA estimator satisfies a bound whose leading term is

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

with a higher-order remainder of order $\dot\theta=\bar A\theta,$5 up to problem-dependent factors. The same work proves a matching local asymptotic minimax lower bound, so the averaged SA estimator is instance-optimal up to universal constants. The last iterate attains the correct worst-case scaling, but the averaged estimator captures the correct instance-dependent covariance structure (Mou et al., 2021).

For Polyak–Ruppert averages with decaying stepsizes $\dot\theta=\bar A\theta,$6, $\dot\theta=\bar A\theta,$7, Berry–Esseen theory is available under uniformly geometrically ergodic Markov sampling. In one-dimensional projections,

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

admits a Gaussian approximation with Kolmogorov error of order $\dot\theta=\bar A\theta,$9, up to logarithmic factors, and the optimized choice $(Z_n)$00 yields the rate

$(Z_n)$01

The same analysis establishes the non-asymptotic validity of a multiplier subsample bootstrap with block length $(Z_n)$02, giving approximation error of order $(Z_n)$03, and an overlapping-batch-means variance estimator with error

$(Z_n)$04

for suitable gain choices. These results turn asymptotic CLTs for MLSA into finite-sample confidence-interval statements under dependent sampling (Samsonov et al., 25 May 2025).

A diffusion perspective has recently been developed for constant-stepsize linear TD under Markovian noise. In that setting, the ODE

$(Z_n)$05

captures only the mean contraction, while the recursion is weakly approximated by the SDE

$(Z_n)$06

Here $(Z_n)$07 is the long-run Markov covariance defined through the Poisson equation

$(Z_n)$08

equivalently

$(Z_n)$09

The resulting stability bound has the form

$(Z_n)$10

so the familiar constant-stepsize error floor is interpreted as the balance between deterministic damping and Markovian long-run covariance. Near the fixed point, the scaled error converges to an Ornstein–Uhlenbeck approximation with stationary covariance solving

$(Z_n)$11

This makes the geometry of slowly contracting directions explicit (Forzo et al., 16 Jun 2026).

6. Reinforcement-learning instances and major extensions

Temporal-difference learning is the canonical MLSA application. For discounted policy evaluation with truncated eligibility traces,

$(Z_n)$12

the TD$(Z_n)$13 update

$(Z_n)$14

fits the linear SA template with augmented state $(Z_n)$15. Under a super-Lyapunov drift condition on the underlying Markov kernel and growth assumptions such as

$(Z_n)$16

the MLSA finite-time theory applies to TD(0), TD$(Z_n)$17 with truncated traces, and related linear value-function estimators on unbounded state spaces (Durmus et al., 2021). In the bounded finite-state setting, constant-stepsize finite-time guarantees for TD(0) and TD$(Z_n)$18 are available without assuming i.i.d. samples and without requiring projection (Srikant et al., 2019).

Two-timescale extensions are central for gradient-corrected RL methods. For linear two-timescale SA with Markovian noise, one line of analysis shows that Markov dependence does not alter convergence rates relative to martingale noise: it affects only the constants through Poisson-equation bounds. Under $(Z_n)$19 and $(Z_n)$20 Hurwitz and suitable stepsizes, the slow iterate satisfies

$(Z_n)$21

while the fast tracking error is $(Z_n)$22, and an asymptotic expansion yields a matching lower bound $(Z_n)$23 (2002.01268). A sharper finite-time theory for Markovian linear two-timescale SA identifies the exact CLT covariance in the leading term: $(Z_n)$24 with $(Z_n)$25 solving a coupled Lyapunov system. This produces tight $(Z_n)$26-order bounds for the slow variable and applies to TDC, GTD, GTD2, and Polyak averaging as a special case (Haque et al., 2023).

Almost-sure convergence without projection has also been established for two-timescale SA under Markovian noise. The key device is to control the fast iterate by the running maximum of the slow iterate rather than by the current slow iterate: $(Z_n)$27 Under standard two-timescale step-size separation and globally asymptotically stable limiting ODEs, this yields boundedness and convergence of the coupled iterates. As a concrete consequence, off-policy TDC$(Z_n)$28 with eligibility traces and linear function approximation converges almost surely under Markovian sampling, without projection and without requiring a compact noise space (Mahadevan et al., 29 May 2026).

MLSA also extends beyond single-agent policy evaluation. In distributed linear SA over uniformly strongly connected time-varying directed graphs, each agent follows a local SA recursion coupled through stochastic consensus weights. Finite-time mean-square error bounds are available for both constant and diminishing gains. When the interaction matrices are merely stochastic, the limiting ODE equilibrium corresponds to a $(Z_n)$29-weighted convex combination of local terms rather than the straight average; a push-sum-type variant restores the straight average under one-way communication and yields a separate finite-time guarantee (Lin et al., 2021).

Finally, MLSA has been used as a vehicle for estimating Markov-chain asymptotic variance through Poisson-equation reformulation. In the tabular setting and in linear function approximation, the recursive estimator is itself a linear SA with Markovian noise, and with diminishing stepsizes it achieves

$(Z_n)$30

The same construction leads to TD-like algorithms for average-reward RL when asymptotic variance is used as a risk measure (Agrawal et al., 2024).

Across these developments, three broad points recur. First, Markovian dependence is not a minor perturbation of i.i.d. noise: it can create $(Z_n)$31 stationary bias, alter long-run covariance, and invalidate light-tail heuristics. Second, those effects can nonetheless be analyzed sharply through Lyapunov equations, Poisson equations, and random matrix products. Third, the linear case is technically rich enough to support general-state-space stability theory, precise finite-time rates, inference procedures, and reinforcement-learning algorithms ranging from TD learning to distributed and two-timescale methods.