Online Covariance Estimation in Averaged SGD: Improved Batch-Mean Rates and Minimax Optimality via Trajectory Regression

Abstract: We study online covariance matrix estimation for Polyak--Ruppert averaged stochastic gradient descent (SGD). The online batch-means estimator of Zhu, Chen and Wu (2023) achieves an operator-norm convergence rate of $O(n^{-(1-α)/4})$, which yields $O(n^{-1/8})$ at the optimal learning-rate exponent $α\rightarrow 1/2^+$. A rigorous per-block bias analysis reveals that re-tuning the block-growth parameter improves the batch-means rate to $O(n^{-(1-α)/3})$, achieving $O(n^{-1/6})$. The modified estimator requires no Hessian access and preserves $O(d^2)$ memory. We provide a complete error decomposition into variance, stationarity bias, and nonlinearity bias components. A weighted-averaging variant that avoids hard truncation is also discussed. We establish the minimax rate $Θ(n^{-(1-α)/2})$ for Hessian-free covariance estimation from the SGD trajectory: a Le Cam lower bound gives $Ω(n^{-(1-α)/2})$, and a trajectory-regression estimator--which estimates the Hessian by regressing SGD increments on iterates--achieves $O(n^{-(1-α)/2})$, matching the lower bound. The construction reveals that the bottleneck is the sublinear accumulation of information about the Hessian from the SGD drift.

• The paper presents an improved batch-means estimator that, through optimal block-growth tuning and burn-in, achieves an operator-norm convergence rate of O(n^{-(1-α)/3}).
• It rigorously decomposes the estimation error into variance, stationarity bias, and nonlinearity bias, clarifying the bias-variance tradeoff in online covariance estimation.
• The proposed trajectory regression estimator attains minimax optimality, offering a fully online, Hessian-free approach for accurate uncertainty quantification in stochastic optimization.

Online Covariance Estimation in Polyak–Ruppert Averaged SGD: Improved Batch-Means Rates and Minimax Optimality

Problem Formulation and Motivation

Covariance estimation for Polyak–Ruppert averaged stochastic gradient descent (SGD) is fundamental for valid statistical inference in stochastic optimization, especially under large-scale or online settings where only first-order information is accessible. The asymptotic covariance $V = H^{-1} S H^{-1}$, depending on the unknown Hessian $H$ and gradient noise covariance $S$, governs the distributional behavior of estimators derived from averaged SGD. Classical methods that require access to the Hessian or explicit computation of all iterates may not be computationally feasible or compatible with online inference. This work investigates efficient, Hessian-free estimators of $V$ which can be computed online, and explores their statistical optimality.

Improved Batch-Means Estimators: Structural Advances

The traditional online batch-means estimator partitions iterates into blocks of increasing size and forms an average of block-wise covariance estimates. Zhu et al. formalized the variance-bias tradeoff for such estimators and established an operator-norm convergence rate of $O(n^{-(1-\alpha)/4})$ for learning-rate $\eta_t = \eta_0 t^{-\alpha}$, $1/2 < \alpha < 1$.

This work rigorously decomposes the estimation error into variance, stationarity bias, and nonlinearity bias. The primary advance is an analytic characterization of per-block bias, revealing it scales with the ratio of mixing time to block size. This leads to two principal results:

1. Optimal Block-Growth Tuning: By carefully tuning the block growth parameter $\beta$, the batch-means estimator achieves an improved operator-norm convergence rate of $O(n^{-(1-\alpha)/3})$ (i.e., $O(n^{-1/6})$ at $H$0), strictly faster than previous bounds.
2. Burn-In and Weighted Variants: Introducing a burn-in fraction parameter $H$1 (discarding a fraction of initial, nonstationary blocks) or adopting smooth block weightings achieves the same rate as hard truncation. While burn-in does not alter the asymptotic exponent, it significantly reduces finite-sample bias.

   Figure 1: Operator-norm error for different estimators and learning-rate exponents; improved rates with optimally tuned block-growth and burn-in are evident across all tested regimes.

The full error is given by the sum of variance $H$2, stationarity bias $H$3, and nonlinearity terms (subdominant at optimal tuning). Tuning $H$4 to balance variance and bias yields the stated rate, and the result holds consistently for both hard and soft block weighting schemes.

Minimax Lower Bounds and Trajectory Regression

A lower bound argument is constructed using Le Cam’s two-point method. By considering two quadratic objectives with identical noise but distinct Hessians and analyzing the KL divergence between the induced SGD trajectories, it is proven that no Hessian-free trajectory-based estimator can consistently outperform the rate $H$5, corresponding to $H$6 at $H$7. This serves as a fundamental limitation, strictly separating the best achievable rate from the i.i.d.\ optimal regime (which admits $H$8).

To match this lower bound, a trajectory regression estimator is proposed: the Hessian is estimated by regressing SGD increments $H$9 onto iterates $S$0, and the noise covariance is estimated from the regression residuals. The plug-in covariance then achieves the minimax rate. This estimator is fully online and Hessian-free, requiring only standard matrix inversions at termination.

Empirical Validation of Bias–Variance Structure

Experiments on synthetic quadratic objectives with varying $S$1 empirically substantiate three findings. First, the improved batch-means with optimal block growth and burn-in exhibits lower error and adheres to theoretical slopes. Second, burn-in significantly reduces bias attributable to early blocks, boosting performance for all sample sizes. Third, the variance-bias decomposition displays an early-block regime where bias dominates, justifying the exclusion of these blocks for practical inference.

Figure 2: Per-block bias and variance as functions of block index for $S$2 and $S$3; burn-in excludes the high-bias initial blocks from contributing to the final estimate.

Theoretical and Practical Implications

The improved analysis demonstrates that classical batch-means estimators, when carefully analyzed and optimally tuned, approach but do not reach the minimax rate, underscoring a practical tradeoff between robustness and optimality. The trajectory-regression approach achieves minimaxity but can be less robust to poor conditioning or model deviations. The explicit finite-sample constants quantified by the improved bias term suggest that moderate burn-in and block growth rates close to but above the mixing threshold are beneficial in operational deployments.

Practically, the results enable more accurate uncertainty quantification, including confidence intervals and hypothesis testing, for large-scale stochastic optimization problems where Hessian access is infeasible. They also lay the groundwork for future studies on high-dimensional ($S$4) behavior, Markovian data settings, and Berry–Esseen-type coverage guarantees.

Conclusion

This work establishes a rigorously optimized batch-means framework for online covariance estimation in SGD, demonstrating improved convergence rates via fine-grained bias decompositions and burn-in strategies. The analysis clarifies the minimax limits of Hessian-free estimation, positioning trajectory-regression as rate-optimal and batch-means as a highly practical, near-optimal method. Open questions remain concerning batch-means optimality in broader stochastic settings, extensions to non-Euclidean geometries, and the empirical impact of estimator choice on downstream inferential tasks.

---

References:

• "Online Covariance Estimation in Averaged SGD: Improved Batch-Mean Rates and Minimax Optimality via Trajectory Regression" (2604.10814)