Online Covariance Estimation for Averaged Stochastic Gradient Descent on Riemannian Manifolds

1 Introduction

• Existing online Euclidean batch-means estimators use increasing-size blocks and can operate without knowing the terminal sample size, but their constructions do not directly address manifold geometry.A Riemannian gradient lies in the tangent space at its current point, and the Euclidean update need not preserve the manifold constraint.Source: §1 · p. 2

1.1 Related work

• Prior work established asymptotic normality for averaged RSGD, while this paper focuses on estimating its limiting covariance as the tangent spaces evolve.The authors distinguish their focus from related inference for constrained and nonsmooth stochastic optimization, which they describe as evolving along other lines.Source: §1 · p. 3, §1.1 · p. 4

1 Introduction

• Changing averages also change tangent spaces, so the method recursively updates transported batch-average approximations rather than recomputing every inverse retraction.The paper separates Euclidean batch-means error from geometric error and shows the latter is higher order under standard conditions.Source: §1 · p. 3
• The paper constructs overlapping and non-overlapping online batch-means covariance estimators for averaged RSGD, requiring neither Riemannian Hessian computations nor a prespecified terminal sample size.Under the conditions of Theorem 3.1, an appropriate batch exponent yields a rate matching the Euclidean estimator of Zhu et al. (2023).Source: §1 · pp. 3–4

1.2 Paper organization

• The paper is organized from geometric preliminaries and estimator theory through inference, streaming PCA experiments, and future directions.Its notation section specifies conventions for inner products, norms, linear maps, outer products, and probability quantiles.Source: §1.2 · p. 4, §1.3 · p. 5

1.3 Notation

• Notation includes the outer product (u⊗v)w = ⟨v,w⟩u and uses z_{1−δ/2} and χ²_{d,1−δ} for normal and chi-squared quantiles.The paper also defines conditional expectation given the observations through iteration n.Source: §1.3 · p. 5

2.1 Riemannian geometry

• A retraction maps tangent vectors back to the manifold and has a local inverse that represents nearby points as tangent vectors; vector and parallel transport compare vectors across tangent spaces.Parallel transport preserves inner products and norms, while the vector transport considered here is defined using the retraction curve.Source: §2.1 · p. 6
• On the sphere, the normalization retraction follows a great circle through angle arctan(∥v∥), while the exponential map travels through angle ∥v∥.The sphere’s tangent space at x consists of vectors orthogonal to x.Source: §2.1 · p. 7

2.2 Averaged RSGD and asymptotic normality

• Averaged RSGD takes a retraction-based stochastic-gradient step and moves its running average toward each new iterate along a retraction curve with weight 1/n.In Euclidean space, these recursions reduce to the usual SGD update and recursive arithmetic average.Source: §2.2 · p. 7
• The assumptions keep iterates and averages in a compact region where inverse retractions are defined, and require local smoothness, positive-definite curvature, gradient-noise conditions, and moment bounds.The second-moment bound supports asymptotic normality; the fourth-moment bound is needed for covariance estimation.Source: §2.2 · pp. 8–10
• For RSGD iterates and recursive averages with γ_n=γn^-α (γ>0, α∈(1/2,1)), under Assumptions 2.1–2.3 and 2.4(a), √n eΔ_n converges in distribution to N(0,Σ), where Σ=A^-1SA^-1.Source: §2.2 · p. 10
• The limiting covariance depends on unknown operators, so practical inference requires consistent online covariance estimators that use RSGD iterates without evaluating the Riemannian Hessian.Combining those estimators with asymptotic normality yields confidence intervals and confidence regions.Source: §2.2 · p. 10

3.1 Online batch-means covariance estimators

• Batch means use retraction-based averaging, with deviations recursively recentered and transported; OL selects every batch, while NOL selects disjoint batches that partition the iterates.Source: §3 · p. 10, §3.1 · pp. 11–12

3 Online covariance estimation on Riemannian manifolds

• Because the target covariance acts on the unknown minimizer’s tangent space, the estimators are instead constructed on the tangent space at the current averaged iterate.Source: §3 · p. 10

3.2 Online implementation

• The online implementation maintains aggregate sums for retained batches and adjusts them as the overall average moves, avoiding revisiting individual batches.Source: §3.2 · pp. 13–14

3.3 Convergence rates

• Both OL and NOL converge to Σ at rates matching Euclidean counterparts; choosing β = 2/(1 − α) optimizes the rate, which approaches n^-1/8 as α approaches 1/2 from above.Source: §3.3 · p. 15

3.4 Proof strategy

• The proof separates geometry-induced errors from the Euclidean batch-means error; for α > 1/2, the geometry errors do not affect the convergence rate.Source: §3.4 · pp. 15–16

4.1 Confidence intervals and regions

• A scalar interval estimates h(x⋆) at exn using a gradient-based variance estimate and standard error scaled by 1/√n.The interval uses the (1−δ/2)-quantile of the standard normal distribution.Source: §4.1 · p. 17
• For joint inference, the covariance-shaped tangent ellipsoid is retracted from TexnM onto the manifold to form a confidence region.In Euclidean space, the interval and region reduce to the coordinatewise interval and ellipsoidal region of Zhu et al. (2023, Corollary 3.2).Source: §4.1 · pp. 17–18, Figure 2 · p. 17
• Under the stated conditions, the interval and region constructions have asymptotically correct coverage for smooth h with nonzero gradient at x⋆.If the covariance estimate is singular, the confidence region is set to {exn}.Source: §4.1 · p. 18

4.2 Streaming principal component analysis (PCA)

• In streaming PCA, the Grassmann manifold represents the target leading subspace, which is unique when λk > λk+1.Its dimension is k(p−k), and the section specifies the geometric maps and Riemannian gradient used to implement RSGD and its covariance estimators.Source: §4.2 · pp. 18–19
• For PCA, the polar retraction maps a tangent ellipsoid at the averaged estimate to subspaces, while hq measures the squared projection of q onto the target subspace.The scalar feature is independent of the orthonormal representative; its estimated variance sets the interval width.Source: §4.2 · pp. 19–21

5.1 Experimental setup

• Experiments generate i.i.d. Gaussian observations, set k = ⌊√p⌋, choose a positive eigengap Δλ, and use 200 independent runs per setting, each with N = 3×10^5 observations.The limiting covariance is available in closed form to quantify covariance estimation errors and compare confidence intervals with oracle intervals based on the true asymptotic variance.Source: §5.1 · p. 21, §5 · p. 21

5.2 Covariance estimation

• Both estimators’ covariance errors generally decrease with sample size, with log-log slopes near −1/8; smaller eigengaps have larger absolute errors.Because the limiting covariance also grows as the eigengap decreases, larger absolute error does not necessarily mean lower relative accuracy.Source: §5.2 · p. 22, Figure 4 · p. 23
• At n = 3×10^5, c = 1 gives the smallest mean covariance error for both estimators across all settings considered.The authors use c = 1 for the confidence-interval experiments.Source: §5.2 · p. 22

5.3 Statistical inference

• For inference, the study compares nominal 95% intervals for hq(x⋆) with oracle intervals across G10,3 and G20,4 at eigengap 0.25.The oracle intervals use the true scalar asymptotic variance; empirical coverage is the fraction of 200 intervals containing hq(x⋆).Source: §5.3 · pp. 22–23, Figure 5 · p. 24
• At n = 3×10^5, OL/NOL coverage is 0.895/0.920 on G10,3 and 0.895/0.895 on G20,4, below nominal 95% but about 90%.Mean interval lengths decrease with sample size and are similar for OL and NOL.Source: §5.3 · p. 24, Figure 5 · p. 24

6 Conclusion and future directions

• The overlapping and non-overlapping estimators are the first of their type for averaged RSGD and need neither Riemannian Hessians nor a prespecified terminal sample size.The geometric errors induced by the manifold structure are of smaller order than the leading batch-means error.Source: §6 · p. 24
• Under standard assumptions, an appropriate batch-size control sequence gives rate n^-(1-α)/4, matching the Euclidean estimator of Zhu et al. (2023).Combined with asymptotic normality, consistency yields asymptotically valid confidence intervals for smooth scalar functions and joint confidence regions for x⋆.Source: §6 · pp. 24–25
• Streaming PCA experiments on Grassmann manifolds support estimator consistency and demonstrate confidence-interval construction in higher-dimensional settings.Source: §6 · p. 25
• Future work includes faster online covariance estimators and lower bounds to characterize minimax rates for RSGD iterates on manifolds.The authors note that faster covariance estimation may improve finite-sample inference performance.Source: §6 · p. 25

Disclosure of AI usage

• OpenAI’s GPT 5.6 Sol and earlier models assisted with proofs of Lemmas A.6 and D.8, simulation code, and proofreading.Source: Disclosure of AI usage · p. 25

Supplemental Materials: Online Covariance Estimation for Averaged Stochastic Gradient Descent

• Appendices A–F cover technical lemmas, geometric maps, RSGD iterates, proofs of Theorem 3.1 and Corollary 3.2, inference results, and additional numerical experiments.Source: on Riemannian Manifolds · p. 29

A Technical lemmas

• For a symmetric positive-definite matrix and step sizes γ_k = γk^-α with α ∈ (1/2, 1), the matrix-product recursion yields a second-moment bound of order n^-α.The accompanying lemmas also control weighted products and partial sums of the recursion.Source: A Technical lemmas · pp. 29–31
• A quadratic comparison lemma bounds differences between finite weighted sums of random-vector outer products.It applies to square-integrable vectors with nonnegative deterministic weights.Source: A Technical lemmas · p. 32
• The OL and NOL batch constructions satisfy deterministic block bounds under the stated growth conditions, including for the polynomial schedule used in Theorem 3.1.These estimates cover batch lengths, normalizing factors, and weighted sums used in the estimator proofs.Source: A Technical lemmas · pp. 32–35

B Manifold lemmas

• Under the stated compactness and retraction assumptions, the inverse retraction map has a jointly C2 extension near X × X, with C3 extensions when the retraction has C3 regularity.The appendix also establishes uniform derivative bounds for these maps.Source: B Manifold lemmas · p. 35
• The manifold lemmas establish uniform comparison between the distance d(x, y) and the norm of the inverse retraction vector for x, y ∈ X.The bounds use compactness and the local behavior of the exponential map.Source: B Manifold lemmas · pp. 36–37
• Further retraction and transport estimates provide uniform derivative control and Lipschitz dependence under the specified assumptions.The lemmas use compactness, retraction properties, and regularity of the inverse extension.Source: B Manifold lemmas · pp. 37–38
• Transport along retraction curves and its inverse are C2 under Assumption 2.1(e).The transport is characterized as the solution of a linear differential equation with C2 coefficients.Source: B Manifold lemmas · p. 38

C Algorithmic estimates

• The transported stochastic gradient decomposes into the Hessian term AΔ_n−1, gradient noise at x⋆, and stochastic and deterministic remainders.The expansion defines the linearization around x⋆ used in the algorithmic estimates.Source: C Algorithmic estimates · pp. 40–41
• Recursive and arithmetic averages differ in first-moment norm by at most Cγ_s over an interval, with corresponding second-moment bounds.The bounds also cover global averages by taking s = 1 and j = n.Source: C Algorithmic estimates · pp. 45–46

D Proof of Theorem 3.1

• The batch definitions use iterates from t_j through j, with their recursive mean and its exact and approximate tangent-space representations.The operators are extended to either batch selection using averages defined earlier in the paper.Source: D Proof of Theorem 3.1 · pp. 46–47
• The four terms represent online vector updates, changes of base point, manifold averaging, and batch-means error in a fixed tangent space.The operators act on T_x⋆M, with parallel transport used to identify the relevant tangent spaces.Source: D Proof of Theorem 3.1 · p. 47
• For both OL and NOL, the proof bounds the fixed-tangent-space batch-means error separately before combining it with the other terms.The remaining three terms are each bounded by C n^(1/2−1/(2β)−α).Source: D Proof of Theorem 3.1 · pp. 46–47

D.1 Batch-means error in the fixed tangent space

• The population covariance is Σ = A^-1SA^-1 = E[(A^-1ζ_1)⊗2], with A the Hessian at x⋆ and S the gradient covariance there.The noise vectors are defined as ζ_k := −grad f(x⋆, ξ_k) and are iid.Source: §D.1 · p. 47
• For the noise-based estimator, independence and centering make the uncentered estimator’s expectation equal to Σ; the proof then restores the centering term.The argument controls cross terms separately for OL and NOL using their batch and block structure.Source: §D.1 · pp. 48–50
• The section establishes bounds for noise-based, linearized, and RSGD covariance estimators under the stated assumptions for both OL and NOL.The proof compares estimators based on A^-1ζ_k, U_k, and Δ_k, then combines the component bounds.Source: §D.1 · pp. 51–52

D.2 Geometric errors

• The geometric-error analysis separately bounds recursive retraction averaging, base-point changes, and errors from online vector updates.The three comparisons are made against arithmetic averaging in T_x⋆M, transported batch quantities, and exact batch-vector updates, respectively.Source: §D.2 · pp. 53–54
• For either OL or NOL, the proof combines the four decomposition-term bounds; the final remainder is negligible relative to n^-1/(2β) when α > 1/2.Its ratio to n^-1/(2β) is n^(1/2−α), which tends to zero.Source: §D.2 · p. 56
• Choosing β = 2/(1−α) equates the theorem’s two decay exponents and gives the resulting exponent (1−α)/4.The proof states that this choice satisfies the required strict lower bound on β.Source: §D.2 · p. 56

E Proof of Proposition 4.1

• The transported covariance estimate converges in probability to Σ^-1, supporting the joint-region argument through Slutsky’s theorem.The proof also bounds the largest eigenvalue of the estimated covariance operator in probability.Source: E Proof of Proposition 4.1 · pp. 56–57
• The joint-region argument maps a tangent-space ellipsoid through the retraction and relates its containment of the transported x⋆ to coverage of x⋆.The proof uses a ball on which the retraction is one-to-one with probability tending to one.Source: E Proof of Proposition 4.1 · p. 57
• The scalar interval proof uses a first-order expansion of h composed with the retraction and concludes the scalar confidence-interval assertion.The argument assumes the directional variance ⟨grad h(x⋆), Σ grad h(x⋆)⟩ is positive.Source: E Proof of Proposition 4.1 · p. 57

F Additional numerical experiments

• The appendix uses Gaussian-model experiments with 200 independent runs per setting to examine covariance errors and statistical inference.Source: F Additional numerical experiments · p. 57

F.1 Covariance estimation in higher dimensions

• The experiments cover G15,3, G20,4, and G25,5, whose manifold dimensions are 36, 64, and 100, with eigengaps 1 and 0.25 and batch constants 1, 2, and 5.Source: §F.1 · p. 58
• At n = 3×10^5, c = 1 gives the lowest mean error for both OL and NOL across all six settings.The settings span three manifolds and two eigengaps.Source: Figure 6 · p. 58, §F.1 · p. 58

• The experiments construct nominal 95% intervals for hq(x⋆) using OL and NOL variance estimates with c = 1, and compare them with oracle intervals.Source: §F.2 · p. 58
• The experiments cover twelve settings at n = 3 × 10^5, and Figure 7 tracks coverage and mean interval length over sample size for G3,1 and G25,5.Source: §F.2 · p. 58, Figure 7 · p. 60
• At n = 3 × 10^5, OL and NOL coverage ranges from 0.880 to 0.925, while their mean interval lengths are similar within each setting.On G25,5 with eigengap 0.25, coverage is 0.885 for OL and 0.880 for NOL, below the nominal 0.95 level.Source: §F.2 · p. 58, Table 2 · p. 59
• Across the plotted sample sizes, OL and NOL have similar coverage and mean interval lengths, with coverage around 90% at larger sample sizes.Intervals become shorter on average as sample size increases.Source: §F.2 · p. 59, Figure 7 · p. 60