Online Covariance Estimation for Averaged Stochastic Gradient Descent on Riemannian Manifolds

Jingrong Guan; Liwei Jiang; Sen Na

stat.MLOct 24, 2026

Online Covariance Estimation for Averaged Stochastic Gradient Descent on Riemannian Manifolds

Authors: Jingrong Guan, Liwei Jiang, Sen Na

Organizations: Department of Industrial Engineering and Management Sciences, Northwestern University · Edwardson School of Industrial Engineering, Purdue University · School of Industrial and Systems Engineering, Georgia Institute of Technology

Abstract

We study online inference for averaged Riemannian stochastic gradient descent (RSGD) when the unknown parameter lies on a manifold. Although the averaged RSGD iterate is asymptotically normal, practical inference requires estimating the unknown covariance operator in its limiting distribution. Existing online batch-means methods cannot be applied directly because a general Riemannian manifold lacks an addition operation, and the tangent space containing the Riemannian gradient varies across iterations. Through vector transport, we construct overlapping and non-overlapping covariance estimators that adaptively combine information from successive batches while accounting for the changing tangent spaces. The proposed estimators are online, require neither Riemannian Hessian computations nor a prespecified terminal sample size, and attain the same convergence rates as their Euclidean counterparts (Zhu et al., 2023). Together with asymptotic normality, these estimators enable online statistical inference for RSGD on manifolds. We illustrate the proposed framework through a detailed study of streaming principal component analysis on Grassmann manifolds.

Figures & tables

Figure 1: Left: Assumption 2.1 (a). The inverse image

R_{x}^{-1}(\mathcal{X})

lies in the tangent ball

B_{x}(0_{x},\rho_{\mathrm{ret}})

for every

x\in\mathcal{X}

, with a common radius

\rho_{\mathrm{ret}}>\tan(\max_{x,y\in\mathcal{X}}d(x,y))

. Right: Assumption 2.1 (b). For any

x,y\in\mathcal{X}

, the normalization retraction maps the line segment

\{tR_{x}^{-1}(y):t\in[0,1]\}

\mathrm{T}_{x}\mathcal{M}

onto the short great circle arc from

x

y

, which remains in

\mathcal{X}

, illustrating its retraction convexity.

Figure 2: Construction of a confidence region on the manifold. The retraction

R_{\widetilde{x}_{n}}

maps the tangent ellipsoid

\mathcal{E}_{n,\diamond}(1-\delta)

to the shaded region

\mathcal{C}_{n,\diamond}(1-\delta)

, and its center

0_{\widetilde{x}_{n}}

\widetilde{x}_{n}

. The right panel separates the tangent space from the manifold to show corresponding points, connected by the light gray dashed lines. The estimator subscript

\diamond

is omitted in the figure.

Figure 3: Spherical representation of a confidence region on

G_{3,1}

. Left: the black and orange lines represent the estimate

[\widetilde{X}_{n}]

and a subspace

[z]

in the confidence region, respectively. Right: the polar retraction maps each tangent vector

v

to the line represented by

z=(\widetilde{X}_{n}+v)/\sqrt{1+\|v\|^{2}}

. These unit representatives form the shaded spherical cap depicting the confidence region. The subscript

\diamond

is omitted.

Figure 4: Gaussian streaming PCA: log mean covariance error (operator norm) against log sample size. Columns correspond to

G_{3,1}

G_{7,2}

, and

G_{10,3}

, and the top and bottom rows use eigengaps

1

and

0.25

, respectively. The colored curves compare OL and NOL with

c\in\{1,2,5\}

. Each mean is based on

200

independent runs. The black dashed lines indicate the reference slope

-1/8

Figure 5: Gaussian streaming PCA: empirical coverage (top) and mean length (bottom) of

95\%

confidence intervals for

h_{q}(x_{\star})

. Columns correspond to

G_{10,3}

and

G_{20,4}

, with

\Delta_{\lambda}=0.25

and

c=1

. The blue, orange, and purple curves show OL, NOL, and Oracle, respectively, with results based on

200

independent runs. The small plots inside the top panels show the later coverage values, and the red dashed lines mark the nominal coverage level

0.95

Empirical coverage
$G_{10,3}$ ( $d=21$ )	$n=50{,}000$	$n=100{,}000$	$n=200{,}000$	$n=300{,}000$
OL	0.835 (0.02625)	0.905 (0.02073)	0.910 (0.02024)	0.895 (0.02168)
NOL	0.860 (0.02454)	0.905 (0.02073)	0.930 (0.01804)	0.920 (0.01918)
Oracle	0.930 (0.01804)	0.965 (0.01300)	0.955 (0.01466)	0.945 (0.01612)
$G_{20,4}$ ( $d=64$ )	$n=50{,}000$	$n=100{,}000$	$n=200{,}000$	$n=300{,}000$
OL	0.860 (0.02454)	0.860 (0.02454)	0.890 (0.02212)	0.895 (0.02168)

Table 1: Empirical coverage of nominal

95\%

confidence intervals for

h_{q}(x_{\star})

in Gaussian streaming PCA, with

\Delta_{\lambda}=0.25

and

c=1

. Results are based on 200 independent runs, with Monte Carlo standard errors in parentheses. Here

d=k(p-k)

is the dimension of

G_{p,k}

Appendix figures & tables3 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 6: Gaussian streaming PCA: log mean covariance error (operator norm) against log sample size. Columns correspond to

G_{15,3}

G_{20,4}

, and

G_{25,5}

, and the top and bottom rows use eigengaps

1

and

0.25

, respectively. As in Figure 4 , darker and lighter curves represent OL and NOL, with blue, orange, and green indicating

c=1,2,5

, respectively. Each mean is based on

200

independent runs. The black dashed lines indicate the reference slope

-1/8

		Coverage			Mean length $\times 10^{3}$
$G_{p,k}$	$\Delta_{\lambda}$	OL	NOL	Oracle	OL	NOL	Oracle
$G_{3,1}$	1	0.920 (0.019)	0.920 (0.019)	0.980 (0.010)	4.425	4.447	4.953
$G_{3,1}$	0.25	0.920 (0.019)	0.925 (0.019)	0.980 (0.010)	11.305	11.722	13.424
$G_{7,2}$	1	0.915 (0.020)	0.910 (0.020)	0.960 (0.014)	3.267	3.252	3.488
$G_{7,2}$	0.25	0.915 (0.020)	0.895 (0.022)	0.955 (0.015)	4.581	4.628	5.012
$G_{10,3}$	1	0.915 (0.020)	0.925 (0.019)	0.955 (0.015)	2.224	2.273	2.415

Appendix

Table 2: Empirical coverage and mean length of nominal

95\%

confidence intervals for

h_{q}(x_{\star})

in Gaussian streaming PCA, with

n=3\times 10^{5}

and

c=1

. Results are based on

200

independent runs, with Monte Carlo standard errors in parentheses. Mean lengths are multiplied by

10^{3}

Figure 7: Gaussian streaming PCA: empirical coverage and mean length of

95\%

confidence intervals for

h_{q}(x_{\star})

, with

c=1

. Columns correspond to

G_{3,1}

and

G_{25,5}

; the upper and lower pairs of rows use eigengaps

1

and

0.25

, respectively, with coverage above mean length in each pair. The blue, orange, and purple curves show OL, NOL, and Oracle, respectively, with results based on

200

independent runs. The small plots inside the coverage panels show the later coverage values, and the red dashed lines mark the nominal coverage level

0.95

Covariance Estimation

Covariance estimation develops methods for learning dependence and uncertainty structure from noisy, incomplete, or high-dimensional data.

See all 61 papers

Apr 12, 2026

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

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.

Yijin Ni, Xiaoming Huo

H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology

Polyak-Ruppert Averaging Covariance Estimation Minimax Estimation Stochastic Gradient Descent

Ask PDF

Apr 23, 2026

Refining Covariance Matrix Estimation in Stochastic Gradient Descent Through Bias Reduction

We study online inference and asymptotic covariance estimation for the stochastic gradient descent (SGD) algorithm. While classical methods (such as plug-in and batch-means estimators) are available, they either require inaccessible second-order (Hessian) information or suffer from slow convergence. To address these challenges, we propose a novel, fully online de-biased covariance estimator that eliminates the need for second-order derivatives while significantly improving estimation accuracy. Our method employs a bias-reduction technique to achieve a convergence rate of

n^{(α-1)/2} \sqrt{\log n}

, outperforming existing Hessian-free alternatives.

Ziyang Wei, Wanrong Zhu, Jingyang Lyu +1

Department of Statistics University of Chicago Chicago, IL 60637, USA · Department of Statistics University of California, Irvine Irvine, CA 92617, USA · Department of Statistics University of Wisconsin–Madison Madison, WI 53706, USA

Covariance Estimation Stochastic Gradient Descent

Ask PDF

Apr 3, 2026

Inversion-Free Natural Gradient Descent on Riemannian Manifolds

The natural gradient method is a central tool for statistical optimisation, but its broader application is hindered by the assumption of a Euclidean parameter space, the repeated estimation of the Fisher information matrix (FIM), and the computational cost of its subsequent inversion. This paper proposes an intrinsic, inversion-free natural gradient method for statistical models whose parameters lie on general Riemannian manifolds. Formulating statistical optimisation in this non-Euclidean setting allows for the natural enforcement of parameter constraints, the elimination of non-identifiable parameters, and the exploitation of geodesic convexity. Our algorithm is based on a moving approximation of the inverse FIM, which is maintained directly on the manifold. This approximation is efficiently updated with new score vectors using low-rank matrix identities. We prove almost-sure convergence rates of

O(\log s / s^α)

for the sequence of iterates, and a similar rate for the approximate FIM. A limited-memory variant with sub-quadratic storage complexity is further proposed for large-scale applications. We demonstrate the efficacy of our method on variational Bayes within the Bures-Wasserstein manifold, normalising flows on the Stiefel manifold, and reduced-rank logistic regression.

Dario Draca, Takuo Matsubara, Minh-Ngoc Tran

School of Mathematics and Statistics, University of Sydney, Australia · University of Sydney Business School, Australia

Natural Gradient Methods Riemannian Optimization

Ask PDF

Abstract

Covariance Estimation

Explore similar work

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

Refining Covariance Matrix Estimation in Stochastic Gradient Descent Through Bias Reduction

Inversion-Free Natural Gradient Descent on Riemannian Manifolds