---
title: Polyak-Ruppert Averaged Estimator
url: https://www.emergentmind.com/topics/polyak-ruppert-averaged-estimator
type: topic
---

# Polyak-Ruppert Averaged Estimator

The Polyak–Ruppert averaged estimator is the averaged output of a stochastic approximation procedure, most commonly the arithmetic mean of successive iterates of Robbins–Monro or stochastic gradient descent. In its classical form, one seeks to minimize a smooth function \(f:\mathbb{R}^d\to\mathbb{R}\) with unique minimizer \(\theta^\star\), runs the stochastic recursion
\[
\theta_{n+1}=\theta_n-\gamma_{n+1}\nabla f(\theta_n;X_{n+1}),
\]
and then replaces the last iterate by
\[
\bar\theta_n=\frac1n\sum_{k=1}^n\theta_k.
\]
The central phenomenon is variance reduction: under the standard conditions of the Polyak–Juditsky theory, the averaged estimator attains the \(O(1/n)\) scale with asymptotically optimal covariance, and subsequent work has extended that conclusion to weakly convex, linear, Markovian, reinforcement-learning, multilevel, two-timescale, preconditioned, momentum, streaming, and zeroth-order settings [1709.03342][2109.07117][2112.14582].

## 1. Classical formulation and averaging schemes

In the basic stochastic-gradient setting, the Polyak–Ruppert estimator is the running average of raw stochastic approximation iterates. The canonical update is the Robbins–Monro recursion
\[
\theta_{n+1}=\theta_n-\gamma_{n+1}\nabla f(\theta_n)+\gamma_{n+1}\Delta M_{n+1},
\]
where \((\Delta M_{n+1})\) is a martingale increment with conditional covariance \(S(\theta_n)\), and the averaged estimator is \(\bar\theta_n=\tfrac1n\sum_{k=1}^n\theta_k\) [1709.03342]. In linear stochastic approximation, the same principle appears as averaging iterates of a linear recursion for solving \(A\theta^\star=b\), sometimes in full-average form and sometimes in tail-average form, for example
\[
\bar\theta_n=\frac1n\sum_{k=n}^{2n-1}\theta_k
\]
or, with burn-in, \(\thetā_n=(2/n)\sum_{k=n/2}^{n-1}\theta_k\) [2405.16644][2207.04475].

The averaging operator is not unique. In streaming stochastic approximation with time-varying mini-batches \(n_t\), the averaged estimator is weighted by the arriving batch sizes,
\[
\bar\theta_t=\frac{1}{N_t}\sum_{i=0}^{t-1}n_{i+1}\theta_i,\qquad N_t=\sum_{i=1}^t n_i,
\]
and this admits a recursive implementation [2109.07117]. In multilevel stochastic approximation, the Polyak–Ruppert average is written
\[
\bar\theta_n=\frac1{\bar b_n}\sum_{k=1}^n b_k\theta_k,\qquad \bar b_n=\sum_{k=1}^n b_k,
\]
with a nondecreasing weight sequence \(\{b_k\}\) [1912.08130]. In continuous-time quasi-stochastic approximation, the same idea becomes time averaging over a terminal window,
\[
\bar\theta_T=\frac{1}{T-T_0}\int_{T_0}^T\theta_t\,dt,
\]
with \(T_0<T\) a burn-in time [2206.00814].

These variants share the same structural purpose: they suppress the high-variance component of the raw recursion while preserving the first-order statistical information. This suggests that “Polyak–Ruppert averaging” is best understood as a family of averaging operators applied to stochastic approximation trajectories rather than as a single formula.

## 2. Structural assumptions and asymptotic efficiency

A standard analytical framework assumes \(f\in C^2(\mathbb{R}^d)\), \(L\)-Lipschitz gradient, positive-definite Hessian \(D^2f(\theta^\star)\) at the minimizer, and martingale-difference noise with conditional covariance \(S(\theta_n)\) converging to \(S^\star=S(\theta^\star)\) [1709.03342]. In the strongly convex case, one also imposes moment bounds such as
\[
E[|\Delta M_{n+1}|^{2p}\mid\mathcal F_n]\le \Sigma_p(1+f(\theta_n)^p),
\]
while a weaker alternative replaces global strong convexity by a Kurdyka–Łojasiewicz-type condition
\[
m\le f(x)^{-r}\|\nabla f(x)\|\le M,\qquad r\in[0,\tfrac12],
\]
with bounded \(D^2f\) [1709.03342].

Under these assumptions, the asymptotic covariance of the averaged estimator takes the classical sandwich form. In the smooth optimization setting,
\[
\Sigma^\star=D^2f(\theta^\star)^{-1}S^\star D^2f(\theta^\star)^{-1},
\]
and the leading \(1/n\) term in the mean-squared error is exactly \(\mathrm{Tr}(\Sigma^\star)/n\), which matches the Cramér–Rao lower bound [1709.03342]. In streaming stochastic approximation, the corresponding constant is
\[
\Lambda=\mathrm{tr}\bigl(\nabla^2F(\theta^\star)^{-1}\Sigma\nabla^2F(\theta^\star)^{-1}\bigr),
\]
and the leading term \(\Lambda/N_t\) again attains the Cramér–Rao lower bound [2109.07117].

The same efficiency interpretation persists in reinforcement learning. For averaged tabular Q-learning, \(\bar Q_T\) is shown to be a regular asymptotically linear estimator for \(Q^\star\) with influence function
\[
\phi(X_t)=(-\gamma P^{\pi^\star})^{-1}Z_t,
\]
and no other regular estimator can attain smaller asymptotic variance than \(V=\mathrm{Var}(\phi(X_t))\) [2112.14582]. In this sense, Polyak–Ruppert averaging is not merely a variance-reduction heuristic; in the settings above it is an efficiency mechanism that recovers the optimal first-order covariance without prior knowledge of the strong-convexity parameter or Hessian [1709.03342].

## 3. Non-asymptotic theory and the extension beyond strong convexity

A central non-asymptotic result is the bound established for step sizes \(\gamma_n=\gamma n^{-\beta}\) with \(\beta\in(1/2,1)\). Assuming \(L^4\)-consistency of the raw iterates,
\[
E[\|\theta_n-\theta^\star\|^4]=O(\gamma_n^2),
\]
one obtains
\[
E[\|\bar\theta_n-\theta^\star\|^2]\le \frac{\mathrm{Tr}(\Sigma^\star)}{n}+C\,n^{-r_\beta},
\qquad
r_\beta=\min\{\beta+\tfrac12,\,2-\beta\}>1.
\]
Optimizing \(r_\beta\) yields \(\beta=\tfrac34\) and
\[
E[\|\bar\theta_n-\theta^\star\|^2]\le \frac{\mathrm{Tr}(\Sigma^\star)}{n}+O(n^{-5/4}),
\]
so the averaged estimator is non-asymptotically minimax-optimal up to a negligible second-order term [1709.03342].

A notable feature of that result is that it does not require uniform strong convexity. The same theorem applies when \(f\) satisfies the global Kurdyka–Łojasiewicz-type condition, including pathological examples such as online logistic regression and recursive quantile estimation. For logistic loss with bounded design, the paper verifies the global KL condition with \(r=0\) and derives
\[
E[\|\bar\theta_n-\theta^\star\|^2]\le \frac{\mathrm{Tr}(\Sigma^\star)}{n}+O(n^{-5/4}),
\]
while for recursive quantile estimation it gives
\[
E[(\bar\theta_n-q_\alpha)^2]\le \frac{\alpha(1-\alpha)}{p(q_\alpha)\,n}+O(n^{-5/4})
\]
with the Cramér–Rao-optimal leading term [1709.03342].

High-probability analogues are now available in more abstract stochastic approximation models. A general-purpose theorem shows that if one already has concentration bounds for the raw iterates, then the averaged iterate \(y_n=(1/(n+1))\sum_{i=0}^n x_i\) satisfies a sharp non-asymptotic bound with leading order \(O((\ln(1/\delta)+1)/n)\), and a matching lower-bound example shows that the resulting \(O(\sqrt{\ln(1/\delta)/n})\) scale cannot be improved in order, up to constants [2505.21796]. For linear stochastic approximation with i.i.d. data or uniformly geometrically ergodic Markov noise, finite-time high-probability bounds recover the Gaussian leading term with covariance \(Ā^{-1}\Sigma Ā^{-1}\), and all dependence on the dimension \(d\) enters only through \(\log d\) in the step size and in exponentially small transient terms [2207.04475].

## 4. Linear, streaming, multilevel, and two-timescale variants

Linear stochastic approximation has provided a particularly precise laboratory for Polyak–Ruppert analysis. In the constant-step Hurwitz case, the averaged iterate obeys a CLT with covariance
\[
\Sigma(\alpha)=\bar A^{-1}\bigl[\Sigma^\star+\Theta(\alpha)\bigr]\bar A^{-T},
\]
where \(\Theta(\alpha)\) is an \(O(\alpha)\) correction to the classical Polyak–Ruppert covariance; in the critical case where \(\bar A\) has non-negative real parts and is diagonalizable, the averaged procedure still achieves an \(O(1/T)\) mean-squared error rate [2004.04719]. This refined picture shows that PR averaging may preserve the optimal \(1/T\) scale even when the drift is only marginally stable.

In the streaming framework, Polyak–Ruppert averaging is combined with time-varying mini-batches \(n_t=\lceil C_\rho t^\rho\rceil\) and learning rates \(\gamma_t=C_\gamma n_t^\beta t^{-\alpha}\). Under smooth, convex, noise, Hessian-Lipschitz, and score-covariance assumptions, the averaged iterate achieves \(E\|\bar\theta_t-\theta^\star\|^2=O(N_t^{-1})\) in both the constant and time-varying mini-batch regimes, and the schedule \(\alpha=2/3,\beta=1/3\) fixes the effective exponent \(\phi=2/3\) for all \(\rho\in(-1,1)\), making the method robust to the arrival rate [2109.07117]. The paper emphasizes that increasing mini-batches lower the intercept of the error curve while averaging yields the flat slope \(-1\) on logarithmic scale [2109.07117].

Multilevel stochastic approximation combines PR averaging with a hierarchy of biased approximations \(F_k(\theta,U)\). Under mild step-size, weight, contractivity, variance, bias, and cost assumptions, the averaged estimator satisfies CLTs both in a slow regime with error \(N^{-\delta}\), \(\delta\in(0,1/2)\), and in a critical regime with error \(N^{-1/2}\sqrt{\log N}\), with sandwich covariance \(H^{-1}\Gamma(H^{-1})^\top\) up to explicit scalar prefactors [1912.08130]. The resulting error-versus-cost rates match classical multilevel Monte Carlo rates up to logarithms [1912.08130].

Two-timescale linear stochastic approximation introduces a further distinction between last-iterate and averaged behavior. For the PR average \(\bar w_N=(1/N)\sum_{k=1}^N w_k\) in a scheme with step sizes \(\beta_k\approx k^{-b}\), \(\gamma_k\approx k^{-a}\), \(1/2<a<b<1\), one obtains the convex-distance normal approximation bound
\[
d_C\bigl(\mathcal L(\sqrt N(\bar w_N-w^\star)),\,N(0,\Sigma_{PR})\bigr)
\le C\bigl(N^{-a/2}+N^{-(1-b)/2}\bigr).
\]
Optimizing \(b-a\to0\) gives rate \(O(N^{-1/4})\), and, unlike the last iterate, the PR-averaged rate is minimized when \(a\approx b\), so larger timescale separation degrades the averaged regime [2508.07928].

A further variant arises in quasi-stochastic approximation. In that deterministic-noise setting, PR averaging over \([T_0,T]\) yields
\[
\bar\theta_T=\theta^\star+a_T c(\kappa,\rho)Y^\star+O(T^{-2\rho}),
\]
and if the design ensures \(Y^\star=0\), then \(\|\bar\theta_T-\theta^\star\|=O(a_T^2)=O(T^{-2\rho})\). Taking \(\rho\) arbitrarily close to \(1\) gives \(\|\bar\theta_T-\theta^\star\|^2=O(T^{-4+\delta})\) for any \(\delta>0\) [2206.00814].

## 5. Gaussian approximation, bootstrap, and online inference

The classical asymptotic normality of PR averages has been strengthened into quantitative Gaussian approximation and bootstrap theory. For decreasing-step linear stochastic approximation, a Berry–Esseen bound in convex distance gives the fastest rate \(O(n^{-1/4})\) at \(\gamma=1/2\), and a multiplier bootstrap based on perturbed LSA recursions yields non-asymptotically valid confidence intervals with the same \(O(n^{-1/4})\) scale [2405.16644]. A later refinement improves the convex-distance approximation to
\[
O(n^{-1/3})
\]
when approximating \(\sqrt n(\bar\theta_n-\theta^\star)\) by the Gaussian law with Polyak–Juditsky covariance \(\Sigma_\infty\), and establishes bootstrap validity at rate \(O(n^{-1/2})\) up to logarithmic and smaller remainder terms [2510.12375].

For nonlinear strongly convex SGD, the averaged estimator admits a decomposition
\[
\sqrt n(\bar\theta_n-\theta^\star)=W_n+D_n,
\]
where \(W_n\) is a linearized sum of independent increments and \(D_n\) collects initial-condition, linearization, and Lipschitz-noise remainders. Using Gaussian approximation for nonlinear statistics, one obtains convex-distance CLT bounds of order up to \(1/\sqrt n\), together with a weighted-SGD multiplier bootstrap that avoids plugging in the limiting covariance \(\Sigma_\infty\) directly [2502.06719]. This places Polyak–Ruppert inference within a non-asymptotic framework rather than a purely qualitative CLT regime.

Reinforcement learning provides a second major inference domain. In tabular synchronous Q-learning, the standardized partial-sum process
\[
S^T(r)=\frac1{\sqrt T}\sum_{t=1}^{\lfloor Tr\rfloor}(Q_t-Q^\star)
\]
converges weakly in \(D([0,1],\mathbb R^d)\) to \(V^{1/2}W(\cdot)\), and the induced FCLT implies a fully online inference method for reinforcement learning [2112.14582]. The same paper proves that \(\bar Q_T\) is semiparametric-efficient and gives a nonasymptotic \(\ell_\infty\)-error bound whose leading term matches the instance-dependent lower bound up to logarithms [2112.14582].

In zeroth-order optimization, the averaged stochastic zeroth-order gradient algorithm \(\bar x_n=(1/n)\sum_{k=0}^{n-1}x_k\) satisfies
\[
\sqrt n(\bar x_n-x_\nu^\star)\xrightarrow{d}N(0,\widetilde V)
\]
for the Gaussian-smoothed objective \(f_\nu\), with
\[
\widetilde V=(\nabla^2 f_\nu(x_\nu^\star))^{-1}S_\nu(\nabla^2 f_\nu(x_\nu^\star))^{-1}.
\]
The paper further constructs an online covariance estimator \(\widehat\Sigma_n\) via overlapping batches and derives asymptotically valid confidence intervals and ellipsoids [2102.05198].

## 6. Adaptive preconditioning, momentum, and limits of averaging

One recurrent misconception is that averaging automatically removes all first-order distortions introduced by algorithmic modifications. Recent work shows a more qualified picture. For time-varying preconditioning, an exact decomposition isolates the preconditioner sequence into a dynamic remainder \(R_n\), and the Polyak–Ruppert CLT
\[
\sqrt n(\bar x_n-x^\star)\xrightarrow{d}N(0,H^{-1}SH^{-1})
\]
holds when the stabilization rate \(\beta\) of \(M_t=(P_tH)^{-1}\) satisfies
\[
\beta>\frac{\alpha+1}{2}.
\]
Within the class of polynomial rate hypotheses used in the upper bound, this threshold is sharp; the paper also verifies \(\beta=1\) for SA-AdaGrad, SA-RMSProp, and SA-ONS, so these methods preserve the same first-order sandwich covariance, and under bounded inputs they also preserve the \(n^{-1/6}\) Wasserstein rate at \(\alpha^\star=2/3\) [2604.23498].

Momentum and non-convergent adaptive preconditioning require a different analysis. For SA-Adam with first-moment buffer \(m_t\), diagonal preconditioner \(P_t\), step sizes \(n_t=n_0t^{-\alpha}\), and momentum gains \(p_t=c_1t^{-\gamma}\) with \(\gamma\in(\alpha,1)\), the augmented-state analysis yields
\[
\sqrt n(I_n-x^\star)\xrightarrow{d}N(0,H^{-1}SH^{-1}),
\]
so the iterate-marginal covariance is exactly the plain SGD sandwich and the adaptivity is asymptotically invisible. The paper stresses that this holds for SA-Adam with sub-linearly vanishing momentum gain, not constant-\(\beta\) deployed Adam; with coupled \(L_2\) weight decay, the limit becomes the ridge-penalized sandwich [2606.17364].

A different limitation appears in constant-step LSA under Markovian noise. There, PR averaging reduces variance but does not remove the invariant-distribution bias: under stationarity,
\[
\lim_{n\to\infty}E[\theta_n^\alpha]-\theta^\star=\alpha\Delta+O(\alpha^{3/2}),
\]
and the paper states explicitly that the \(O(\alpha)\) bias cannot be eliminated by PR averaging [2508.05570]. Richardson–Romberg extrapolation,
\[
\bar\theta_N^{RR}=2\bar\theta_N^\alpha-\bar\theta_N^{2\alpha},
\]
cancels that linear bias term while preserving the \(N^{-1/2}\) variance term with asymptotically optimal covariance \(\Sigma_\infty\) [2508.05570].

Taken together, these results delineate the scope of the Polyak–Ruppert principle. Averaging reliably recovers the optimal first-order variance in a remarkably broad class of stochastic approximation algorithms, but its interaction with constant step sizes, Markov dependence, timescale separation, dynamic preconditioning, and momentum is structure-dependent rather than automatic.

Source: https://www.emergentmind.com/topics/polyak-ruppert-averaged-estimator