Polyak–Ruppert Averaged Iterates

• Polyak–Ruppert averaging is a stochastic approximation technique that averages iterates to achieve optimal asymptotic efficiency and improved finite-sample performance.
• It leverages optimal step-size strategies and sharp non-asymptotic MSE bounds to balance bias and variance, ensuring robust convergence.
• Extensions using the Kurdyka–Łojasiewicz inequality broaden its applicability to degenerate scenarios such as online logistic regression and recursive quantile estimation.

Polyak–Ruppert averaged iterates are a fundamental technique in stochastic approximation and online optimization, designed to improve the statistical efficiency and finite-sample behavior of stochastic gradient-type algorithms. The central idea is to construct an averaged estimator by taking a (typically uniform) average of the iterates produced by a stochastic approximation algorithm. This averaging scheme is known to achieve optimal asymptotic variance in many settings—often matching the Cramér–Rao lower bound—and recent developments have provided sharp non-asymptotic guarantees, expanded the framework to degenerate or non-convex cases, and adapted averaging to complex stochastic processes. The following sections synthesize core principles, sharp theoretical results, and implementation guidance for Polyak–Ruppert averaging, with special emphasis on modern non-asymptotic analysis and its broad applicability.

1. Definition, Classical Results, and Central Limit Theorem

Polyak–Ruppert averaging (sometimes called "stochastic averaged gradient descent," or SAGD) refers to the post-processing of an iterative stochastic algorithm—especially a stochastic gradient descent (SGD)—by averaging all past iterates: $$\hat{\theta}_n = \frac{1}{n} \sum_{k=1}^n \theta_k,$$ where $\theta_k$ are the iterates of SGD or a more general stochastic approximation algorithm.

The classical results, originating with Ruppert (1988) and Polyak & Juditsky (1992), showed that in strongly convex problems with decaying step-size (e.g., $\gamma_n = \gamma n^{-\beta}, ~\beta \in (1/2, 1)$), the averaged estimator $\hat{\theta}_n$ satisfies a central limit theorem: $$\sqrt{n}(\hat{\theta}_n - \theta^\star) \stackrel{d}{\to} \mathcal{N}\left(0, \Sigma^\star\right),$$ with $$\Sigma^\star = \left(D^2 f(\theta^\star)\right)^{-1} S^\star \left(D^2 f(\theta^\star)\right)^{-1}$$, where $D^2 f(\theta^\star)$ is the Hessian at the minimizer and $S^\star$ is the local covariance matrix of the stochastic gradients (Gadat et al., 2017). This is the same covariance as the Cramér–Rao lower bound for semi-parametric estimation, indicating asymptotic efficiency.

2. Non-Asymptotic Mean Square Error Bounds and Optimality

A major advance is the derivation of non-asymptotic, sharp mean square error (MSE) bounds. For general stochastic approximation under step-size $\gamma_n = \gamma n^{-\beta}$, $\beta \in (1/2,1)$, the following bound was established (Gadat et al., 2017): $$\mathbb{E} \left[\|\hat{\theta}_n - \theta^\star\|^2\right] \le \frac{\operatorname{Tr}(\Sigma^\star)}{n} + C n^{-r_\beta},$$ where $$r_\beta = \min\left\{\beta + 1/2, 2 - \beta\right\}$$ and $C$ is a problem-dependent constant.

• For $\beta=3/4$, the "optimal" trade-off yields $r_\beta=5/4$, so the dominant error vanishes at the optimal $O(1/n)$ rate, with the second-order term vanishing faster ($O(n^{-5/4})$). This matches the minimax optimality dictated by the local Cramér–Rao bound and is crucial for sharp finite-sample guarantees (Gadat et al., 2017).

3. Framework Extensions: Beyond Strong Convexity via Kurdyka–Łojasiewicz Inequality

Traditional sharp results depend on the objective $f$ being (uniformly) strongly convex. The framework was expanded by leveraging a global Kurdyka–Łojasiewicz (KL)–type inequality (Gadat et al., 2017), which only requires

$f(x)^{-r}\|\nabla f(x)\| \geq m > 0$

for large $f(x)$ and some $r \in [0, 1/2]$. This prevents the gradient from vanishing too rapidly and ensures sufficient reversion for convergence, even when $f$ is not strongly convex or is non-convex with benign geometry.

This analysis covers "pathological" cases, such as online logistic regression (convex but with degenerate curvature) and recursive quantile estimation (non-convex), where strong convexity fails but the KL inequality permits the recovery of the $O(1/n)$ leading term and asymptotic efficiency.

4. Structural Decomposition and High-Order Error Analysis

The error decomposition underpinning the non-asymptotic results combines spectral analysis of the iterative recursion with precise control of the local noise (variance/covariance structure). Averaged iterates satisfy a higher-order bias/error bound of order $O(n^{-r_\beta})$ with $r_\beta > 1$, which is significant for finite-sample performance. The averaging procedure crucially leverages the $(L^p, \sqrt{\gamma_n})$-consistency of the non-averaged iterates—a property that can be established under strong convexity or the more general KL-type condition with appropriate noise moment assumptions.

This higher-order analysis is also robust to non-convexity as long as the KL-type inequality is satisfied globally, allowing the technique to extend to a variety of stochastic approximation settings.

5. Practical Significance and Applications

• Statistical Efficiency and Finite-Sample Performance: The Polyak–Ruppert averaging not only makes the estimator asymptotically optimal (minimal variance) but—under sharp non-asymptotic control—ensures that the variance term (the only non-negligible term for large samples) dominates, with explicit control over the decay of higher-order terms, both in expectation and in probability.
• Non-Convex and Pathological Scenarios: The KL generalization enables the method to be used in settings such as online logistic regression (where curvature degenerates at infinity) and recursive quantile estimation (where the objective is non-convex), guaranteeing convergence with the same optimal rate as in the strongly convex case.
• Choice of Step-Size: The analysis yields quantitative guidance: for the fastest decay of the second-order term, $\beta = 3/4$ is optimal, balancing the trade-off between bias decay and variance.
• Statistical Inference: Because the method achieves exact asymptotic variance equal to the Cramér–Rao bound, it is an appropriate foundation for developing confidence intervals or sets for $\theta^\star$ in online and streaming learning settings.

6. Key Mathematical Summary

A concise statement of the result (Gadat et al., 2017):

Quantity	Formula/Condition
Averaged iterate	$\hat{\theta}_n = \frac{1}{n}\sum_{k=1}^{n}\theta_k$
Main MSE bound	$$\mathbb{E}\left[\|\hat{\theta}_n-\theta^\star\|^2\right] \leq \frac{\operatorname{Tr}(\Sigma^\star)}{n} + C n^{-r_\beta}$$
Rate exponent ($r_\beta$)	$$r_\beta = \min\left\{\beta + \frac{1}{2}, ~2-\beta\right\}$$, where $\gamma_n = \gamma n^{-\beta}$, $\beta \in (1/2,1)$
Optimal choice	$\beta = 3/4$ yields $r_\beta = 5/4$: $$\mathbb{E}[\|\hat{\theta}_n-\theta^\star\|^2] \leq \frac{\operatorname{Tr}(\Sigma^\star)}{n} + C n^{-5/4}$$
Covariance structure	$$\Sigma^\star = \left(D^2 f(\theta^\star)\right)^{-1} S^\star \left(D^2 f(\theta^\star)\right)^{-1}$$
Generalization (KL-type)	$$\liminf_{|x|\to\infty} f(x)^{-r}\|\nabla f(x)\| > 0$$, enables extension to non-strongly convex and weakly convex/nonconvex settings

7. Concluding Remarks

Polyak–Ruppert averaged iterates provide a general, statistically efficient, and robust approach for stochastic approximation, with precise non-asymptotic analysis now available even beyond the classically required strong convexity setting. The incorporation of KL-type conditions broadens the applicability to degenerate or non-convex landscapes, common in online learning and quantile estimation. These developments yield both optimal asymptotic rate and explicit, fast-decaying finite-sample control, making Polyak–Ruppert averaging a critical tool in modern stochastic optimization and statistical machine learning (Gadat et al., 2017).