$2B$ or Not $2B$: A Tale of Three Algorithms for Streaming: Covariance Estimation after Welford and Chan-Golub-LeVeque

\noindentMSC (2020): Primary 62H12; Secondary 62-08, 65F30, 65Y20.\\ Keywords: streaming covariance, online statistics, Welford algorithm, Chan--Golub--LeVeque, conformal prediction, numerical stability.

Introduction

Covariance estimation is a foundational task in statistics, machine learning, and signal processing. In the classical batch setting all $n$ observations reside in memory and the textbook formula

is applied once. Modern applications routinely violate the batch assumption: sensor streams, federated databases, and large language model training pipelines all generate data faster than can be accumulated. An online algorithm must update its estimate of $\bm{\Sigma}_t$ as each new observation $\mathbf{x}_t$ arrives, using $O(p^2)$ working memory and $O(p^2)$ work per update.

Three strategies address this problem, but a unified treatment covering algebraic equivalence, floating-point error analysis, and statistical uncertainty quantification has not previously appeared. We fill this gap.

\paragraph{The Gram (bariance) algorithm.} It is generally well known (e.g. reichel2025bariance) that the sample variance has an $O(n)$ computation via scalar sums, without explicit centering. Lifting this to the matrix case yields (ref), with the streaming update rule $\mathbf{G}_t \gets \mathbf{G}_{t-1} + \mathbf{x}_t\mathbf{x}_t^\top$ and $\mathbf{s}_t \gets \mathbf{s}_{t-1}+\mathbf{x}_t$. The on-demand formula $\widehat{\bm{\Sigma}}_t = (t\mathbf{G}_t - \mathbf{s}_t\mathbf{s}_t^\top)/[t(t-1)]$ executes in one level-3 BLAS call. Its weakness is catastrophic cancellation when the data mean is large relative to the variance: both $t\mathbf{G}_t$ and $\mathbf{s}_t\mathbf{s}_t^\top$ are $O(t^2\|\bar\mathbf{x}\|^2)$ while their difference is $O(t\|\bm{\Sigma}\|)$.

\paragraph{The Welford algorithm.} welford1962 proposed a shift-invariant one-pass recurrence for scalar variance, extended to covariance by chan1979. The algorithm maintains a running mean $\bm{\mu}_t$ and a matrix $\mathbf{M}_t$ of centred outer-product corrections, updating both in $O(p^2)$ per step. Because corrections are computed relative to the current mean, the method is immune to large constant shifts in the data.

\paragraph{The Chan--Golub--LeVeque (CGL) algorithm.} chan1979 introduced a binary merge formula for independently computed covariance summaries: given $(n_A,\bm{\mu}_A,\mathbf{M}_A)$ and $(n_B,\bm{\mu}_B,\mathbf{M}_B)$, the combined summary is computed exactly. This makes the algorithm tree-parallelisable and directly applicable to federated settings where nodes cannot share raw observations.

\paragraph{Conformal prediction for streaming covariance.} All three algorithms produce point estimates. To quantify uncertainty, we propose a split conformal prediction protocol vovk2005,angelopoulos2023 that constructs a finite-sample, distribution-free confidence interval for each entry $\Sigma_{kl}$ at every step $t \ge 2$ of the stream. The interval has guaranteed marginal coverage $1-\alpha$ for any $\alpha \in (0,1)$, with no distributional assumptions on the data. Experiments show that the coverage guarantee is tight across all three algorithms under well-conditioned data, but the Gram interval collapses (inflates catastrophically) under large data shifts, while Welford and CGL maintain valid coverage.

\paragraph{Contributions.}

Complete proofs for all results appear in Appendix (ref).

Notation and Setup

Let $\mathbf{x}_1,\mathbf{x}_2,\ldots$ be a stream of observations in $\mathbb{R}^p$, $p \ge 2$. After $t \ge 2$ steps the target is the unbiased sample covariance

We write $\mathbf{s}_t = \sum_{i=1}^t \mathbf{x}_i$, $\mathbf{G}_t = \sum_{i=1}^t \mathbf{x}_i\mathbf{x}_i^\top$, $\bm{\mu}_t = \mathbf{s}_t/t$, and $\mathbf{M}_t$ for the Welford/CGL correction matrix. The unit roundoff for IEEE 754 double precision is $\varepsilon_{\mathrm{mach}} = 2^{-53} \approx 1.11\times10^{-16}$. We write $\|\cdot\|_F$, $\|\cdot\|_2$, $\|\cdot\|_{\max}$ for the Frobenius, spectral, and entry-wise maximum norms. For a matrix $A$ the condition number is $\kappa(A) = \|A\|_2\|A^{-1}\|_2$.

Algorithms

Gram (Bariance) Algorithm

\noindentBLAS structure. Each update adds a rank-one correction to $\mathbf{G}$ (DSYR, level-2 BLAS) and a vector increment to $\mathbf{s}$ (DAXPY). Computing $\widehat{\bm{\Sigma}}_t$ on demand requires one outer product (DGER). In batch mode, replacing the loop with a single DSYRK call achieves the optimal level-3 BLAS structure.

Welford Algorithm

The critical property is that both $\bm{\Delta} = \mathbf{x}_t - \bm{\mu}_{t-1}$ and $\mathbf{x}_t - \bm{\mu}_t$ are computed as residuals from means of the same scale as the data, preventing cancellation regardless of the absolute value of $\bm{\mu}_t$.

Chan--Golub--LeVeque Algorithm

Setting block size $b=1$ (each block is a single observation) recovers Algorithm (ref) exactly; block size $b=n$ recovers the batch formula (ref). Because the merge operation is associative and commutative, the algorithm can be applied in any binary-tree order, enabling map-reduce computation.

Algebraic Equivalence

Floating-Point Error Analysis

All three algorithms are algebraically identical but differ in their numerical stability. We quantify two distinct phenomena: rounding accumulation over $t$ steps, and catastrophic cancellation induced by a non-zero data mean.

Rounding Accumulation

Write $\bar x = \|\bm{\mu}_t\|_\infty$ and $\sigma = \|\bm{\Sigma}_t\|_\infty^{1/2}$.

Proofs are given in Appendix (ref). The key contrast: Gram accumulates error at scale $X^2$ (the magnitude of raw observations), while Welford accumulates at scale $\sigma_k\sigma_l$ (the covariance of centred observations). The Welford bound is thus independent of the data mean, as formalised below.

Catastrophic Cancellation

For $c = 10^7$ and $\sigma=1$, the Gram bound gives $pc^2\varepsilon_{\mathrm{mach}} \approx p\cdot 10^{14}\cdot 10^{-16} = p\cdot 10^{-2}$, which becomes non-negligible relative to $\sigma^2=1$. At $c=10^{12}$, the error is $O(p\cdot 10^8)$, fully destroying the estimate (confirmed experimentally in Figure (ref)).

Summary of Bounds

Conformal Prediction for Streaming Covariance

Point estimates from any of the three algorithms carry no automatic uncertainty certificate. We now develop a finite-sample, distribution-free confidence interval for each covariance entry $\Sigma_{kl}$ at every step $t$ of the stream.

Background: Split Conformal Prediction

Split conformal prediction vovk2005,papadopoulos2002,angelopoulos2023 produces a $(1-\alpha)$-coverage prediction interval for a new observation using a held-out calibration set. Given calibration nonconformity scores $s_1,\ldots,s_m$, the conformal quantile is \[ \hat q_\alpha^+ = \mathrm{Quantile}\!\left(\{s_i\}_{i=1}^m;\, \frac{\lceil(m+1)(1-\alpha)\rceil}{m}\right), \] and the interval for a fresh test point has guaranteed marginal coverage $\Pr(\text{true value} \in \hat C_\alpha) \ge 1-\alpha$, with no distributional assumptions on $s_1,\ldots,s_m$ beyond exchangeability.

Protocol for Streaming Covariance

Fix target entry $(k,l)$, nominal coverage $1-\alpha$, and an algorithm $\mathcal{A}\in\{\text{Gram},\text{Welford},\text{CGL}\}$. Let $\Sigma_{kl}$ denote the true population covariance entry.

Computational Cost

The Gram algorithm is faster in batch mode than numpy.cov because it avoids forming and writing the centred matrix $\mathbf{X} - \mathbf{1}_n\bm{\mu}^\top$, saving $O(np)$ memory writes. In streaming mode, Welford computes two outer products per step (for $\bm{\Delta}$ and $\mathbf{x}_t - \bm{\mu}_t$) versus one for Gram, yielding a $\approx 2\times$ FLOP disadvantage; this gap is visible in practice when the outer product (DGER) is the bottleneck.

Experimental Protocol

\paragraph{Hardware and software.} All experiments ran on a macOS 13.0 ARM64 system with 10 CPU cores and OpenBLAS 0.3.31.188.0 via scipy-openblas (USE64BITINT, DYNAMIC_ARCH, NO_AFFINITY, neoversen1, MAX_THREADS=64) under Python 3.12.0 and NumPy 2.4.4 with a fixed random seed. All computations used IEEE 754 double precision.

\paragraph{Runtime measurement.} For each $( n, p )$ pair: (i) generate $\mathbf{X}\sim\mathcal{N}(0,I)$; (ii) run 3 warm-up calls; (iii) record 20--30 wall-clock times via time.perf_counter; (iv) remove outliers using the $1.5\times\mathrm{IQR}$ rule; (v) report the trimmed mean with a 95% bootstrap percentile interval (300 resamples).

\paragraph{Accuracy measurement.} For each data configuration the reference is numpy.cov. We report the entry-wise maximum error $\|\widehat{\bm{\Sigma}} - S_{\mathrm{ref}}\|_{\max}$ and the relative Frobenius error $\|\widehat{\bm{\Sigma}} - S_{\mathrm{ref}}\|_F/\|S_{\mathrm{ref}}\|_F$.

\paragraph{Data configurations.} We test four regimes: (a) Gaussian i.i.d., varying $n$ and $p$; (b) Student-$t_3$ (heavy-tailed), varying $n$; (c) prescribing a range of condition numbers $\kappa(\mathbf{X})\in[10,10^{14}]$; (d) near-singular data with smallest singular value $\sigma_{\min}\in[10^{-12},1]$.

\paragraph{Conformal experiments.} Calibration uses $m=600$ independent trajectories drawn from a $p=5$ multivariate distribution with Toeplitz covariance (entry $(i,j)$: $0.5^{|i-j|}$). Test coverage is estimated from 1200 fresh trajectories. The shift experiment varies the mean offset $c\in\{0,10^3,10^6,10^9,10^{12}\}$ on a $p=4$ identity-covariance distribution and evaluates coverage and interval width at $t=150$.

Results

Runtime

Figure (ref) shows wall-clock time versus $n$ for $p=10$ and $p=50$. Gram is the fastest batch method, consistently beating numpy.cov by a factor of $1.3$--$1.6\times$ due to the avoided $n\times p$ centering write. Welford is slowest by 50--80$\times$ in pure-Python form because the inner loop invokes numpy.outer once per row, incurring $O(n)$ Python-level calls; a compiled implementation would close this to $\approx2\times$ (the FLOP ratio from Table (ref)). CGL (recursive, block size 64) sits between the two.

Figure (ref) shows runtime versus dimension at fixed $n=4{,}000$. For large $p$ all batch methods are dominated by the $O(np^2)$ matrix multiply and their curves converge; Gram retains a small constant-factor advantage from avoiding the centred copy.

Figure (ref) shows the speed ratio of each method relative to Gram.

Numerical Accuracy: Gaussian Data

Figure (ref) shows accuracy on Gaussian i.i.d.\ data. All methods agree with numpy.cov to within $10^{-13}$ in both metrics, consistent with Propositions (ref)--(ref).

Heavy-Tailed and Ill-Conditioned Data

Figure (ref) tests Student-$t_3$ data. Errors remain at floating-point noise for all methods; heavy tails do not affect the relative comparison because all algorithms process the same finite-precision data.

Figure (ref) sweeps the condition number of $\mathbf{X}$. All methods degrade once $\kappa\gtrsim 10^7$, as predicted by standard floating-point theory; Welford and CGL are modestly more robust in the range $\kappa\in[10^6,10^{12}]$.

Catastrophic Cancellation Under Large Shift

Figure (ref) shifts all observations by a constant $c$ and measures error against the unshifted (true) covariance. Gram error grows as $c^2\varepsilon_{\mathrm{mach}}$ (Proposition (ref)), losing $\approx 9$--$10$ decimal digits by $c=10^{12}$. Welford maintains full double-precision accuracy throughout. CGL loses moderate accuracy, reflecting the inter-block cancellation in the merge formula.

Online Fidelity

Figure (ref) tracks $\|\widehat{\bm{\Sigma}}_t^{\mathcal{A}} - \widehat{\bm{\Sigma}}_t^{\mathrm{batch}}\|_{\max}$ as $t$ grows on a zero-mean stream. All methods converge to floating-point noise by $t\approx 100$; the early fluctuations are consistent with the $1/\sqrt{t}$ rate of the estimation error.

Conformal Prediction Results

\paragraph{Coverage and width under well-conditioned data.} Figure (ref) (left) shows empirical coverage versus $t$ for all three algorithms at nominal level $1-\alpha = 95\%$. All algorithms achieve the nominal level (dashed line) at every $t\ge 10$, confirming Theorem (ref). Coverage slightly exceeds 95% for small $t$ because the conformal quantile is conservative for small calibration sets. Figure (ref) (right) shows interval width: all three methods produce intervals of the same order, narrowing as $t$ grows. Gram produces marginally narrower intervals for well-conditioned data because its lower variance translates to tighter calibration scores.

\paragraph{Coverage under large mean shift.} Figure (ref) tests conformal coverage when calibration and test data share the same shift $c$ (left bars) for entry $(1,1)$ (the variance). All three algorithms maintain coverage because the shift is the same in calibration and test sets, so the nonconformity scores are exchangeable. However, Gram's interval widths (right panel) grow dramatically with $c$, reflecting the $O(c^2\varepsilon_{\mathrm{mach}})$ error term inflating the calibration scores. Welford and CGL maintain narrow intervals throughout.

Applications

\paragraph{Federated / privacy-restricted learning.} The Gram sufficient statistics $(\mathbf{G}_t,\mathbf{s}_t,t)$ can be computed locally at each node and aggregated without sharing raw observations. For privacy-sensitive settings where inter-node differences are large, the CGL protocol with block merging provides both the parallelism advantage and the numerical stability of Welford.

\paragraph{Sandwich covariance for M-estimators.} Let $\mathbf{g}_i\in\mathbb{R}^p$ be the score vector for observation $i$. The heteroskedasticity-consistent sandwich covariance white1980,newey1987 is $\hat\Omega = n/(n-1)\cdot (n\mathbf{G} - \mathbf{s}\mathbf{s}^\top)/n^2$ where $\mathbf{G} = \sum_i\mathbf{g}_i\mathbf{g}_i^\top$ and $\mathbf{s} = \sum_i\mathbf{g}_i$. This is precisely the Gram formula applied to the score matrix, admitting a streaming update as each observation is processed.

\paragraph{Panel / fixed-effects data.} For panel unit $i$ observed at $T$ periods, the within-unit covariance is $\hat\bm{\Sigma}_i = (\mathbf{X}_i^\top\mathbf{X}_i - T^{-1}\mathbf{s}_i\mathbf{s}_i^\top)/(T-1)$. The streaming formulation processes units sequentially without forming the full $NT\times p$ matrix.

\paragraph{Conformal uncertainty in online inference.} The conformal intervals of Section (ref) provide rigorous uncertainty quantification for covariance-based online learning algorithms (e.g., online PCA, online Mahalanobis distance) at every step of the stream, not just asymptotically.

Practitioner Decision Guide

Use the Gram algorithm when (a) data are zero-mean or pre-centred, (b) BLAS is well-tuned, or (c) only sufficient statistics can be stored. Use Welford when numerical stability under unknown shifts matters or when simplicity of implementation is paramount. Use CGL with large block size when data are distributed across nodes or when parallel computation is available. Combine any algorithm with the conformal protocol of Section (ref) for principled uncertainty quantification.

Conclusion

We have given a unified treatment of three streaming covariance algorithms---Gram, Welford, and Chan--Golub--LeVeque---covering algebraic equivalence, floating-point error analysis, computational cost, and a new conformal prediction framework. The key empirical findings are: (i) Gram is fastest for batch computation with BLAS, at a factor of $1.4$--$1.6\times$ over numpy.cov; (ii) Welford is the uniquely stable choice under large data shifts, maintaining double-precision accuracy where Gram loses up to 9 decimal digits; (iii) conformal intervals achieve the nominal 95% coverage for all three algorithms on well-conditioned data, but Gram's interval widths inflate catastrophically under large shifts while Welford's remain tight. Together, these results offer practitioners a clear basis for choosing among the three algorithms.

Additional Acknowledgements

The author thanks the numerical analysis literature on which this work builds, particularly higham2002 and chan1979.

Ommited Proofs

Proof of Theorem (ref): Gram Identity

Proof of Theorem (ref): Welford Invariant

Proof of Theorem (ref): CGL Correctness

Proofs of Floating-Point Bounds (Propositions (ref)--(ref))

We use the standard model of floating-point arithmetic higham2002: for any operation $\circ\in\{+,-,\times,\div\}$, $\mathrm{fl}(a\circ b) = (a\circ b)(1+\delta)$ with $|\delta|\le\varepsilon_{\mathrm{mach}}$. Sequential summation of $n$ terms bounded by $M$ satisfies $|\hat S_n - S_n| \le (n-1)\varepsilon_{\mathrm{mach}}\cdot nM + O(\varepsilon_{\mathrm{mach}}^2)$ higham2002.

Proof of Proposition (ref): Cancellation Bound

Proof of Theorem (ref): Conformal Coverage

Bariance Scalar Identity