2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming: Covariance Estimation after Welford and Chan–Golub–LeVeque1 Felix Reichel May 1, 2026
arXiv:2605.00247v1 [stat.CO] 30 Apr 2026
Abstract Three algorithms for computing the unbiased sample covariance matrix in a streaming or distributed setting are placed on a unified algebraic, numerical, and statistical foundation. The Gram algorithm, derived from the bariance reformulation of Reichel [8], maintains the running ∑𝑡 ∑𝑡 cross-product matrix 𝐆𝑡 = 𝑖=1 𝐱𝑖 𝐱⊤ and column-sum vector 𝐬𝑡 = 𝑖=1 𝐱𝑖 , yielding the unbi𝑖 ased covariance in 𝑂(𝑝2 ) per update. The Welford algorithm [10] propagates a running mean and outer-product corrections, achieving the same asymptotic cost with provably better numerical stability under large data shifts. The Chan–Golub–LeVeque (CGL) algorithm [2] supports block-parallel merging via an exact combination formula, making it the natural choice for distributed and map-reduce architectures. All three produce the same estimator in exact arithmetic; their finite-precision behaviour differs markedly. Beyond runtime and numerical comparisons, we introduce a conformal prediction framework for streaming covariance estimation that yields finite-sample, distribution-free confidence sets for each entry of the covariance matrix at any step 𝑡 of the data stream. Experiments confirm that the Gram algorithm is fastest for batch computation, Welford is uniquely robust to catastrophic cancellation under large mean shifts, CGL is optimal for distributed settings, and conformal intervals achieve the nominal coverage level across all three algorithms.
MSC (2020): Primary 62H12; Secondary 62-08, 65F30, 65Y20. Keywords: streaming covariance, online statistics, Welford algorithm, Chan–Golub–LeVeque, conformal prediction, numerical stability.
Contents 1 Introduction
2
2 Notation and Setup
4
3 Algorithms 3.1 Gram (Bariance) Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Welford Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Chan–Golub–LeVeque Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4 4 5 5
4 Algebraic Equivalence
6
5 Floating-Point Error Analysis 5.1 Rounding Accumulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Catastrophic Cancellation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6 6 7
1
Felix Reichel, Graduate Student, Department of Economics, Lund School of Economics and Management (LUSEM), Altenbergerstraße 69, 4040 Linz, Austria. Email: [email protected].
1
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
5.3
Summary of Bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
6 Conformal Prediction for Streaming Covariance 6.1 Background: Split Conformal Prediction . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Protocol for Streaming Covariance . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7 7 7
7 Computational Cost
8
8 Experimental Protocol
9
9 Results 9.1 Runtime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.2 Numerical Accuracy: Gaussian Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.3 Heavy-Tailed and Ill-Conditioned Data . . . . . . . . . . . . . . . . . . . . . . . . . . 9.4 Catastrophic Cancellation Under Large Shift . . . . . . . . . . . . . . . . . . . . . . 9.5 Online Fidelity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.6 Conformal Prediction Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9 9 9 10 12 12 13
10 Applications
14
11 Practitioner Decision Guide
15
12 Conclusion
16
13 Ommited Proofs 13.1 Proof of Theorem 3.1: Gram Identity . . . . . . . . . . . . . . . . . . . . . . . . . . . 13.2 Proof of Theorem 3.2: Welford Invariant . . . . . . . . . . . . . . . . . . . . . . . . . 13.3 Proof of Theorem 3.3: CGL Correctness . . . . . . . . . . . . . . . . . . . . . . . . . 13.4 Proofs of Floating-Point Bounds (Propositions 5.1–5.4) . . . . . . . . . . . . . . . . . 13.5 Proof of Proposition 5.4: Cancellation Bound . . . . . . . . . . . . . . . . . . . . . . 13.6 Proof of Theorem 6.1: Conformal Coverage . . . . . . . . . . . . . . . . . . . . . . . 13.7 Bariance Scalar Identity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17 17 17 18 18 19 20 21
1 Introduction Covariance estimation is a foundational task in statistics, machine learning, and signal processing. In the classical batch setting all 𝑛 observations reside in memory and the textbook formula Σ𝑛 =
) 1 ( ⊤ 𝐗 𝐗 − 𝑛−1 𝐬𝐬⊤ , 𝑛−1
𝐬 = 𝐗⊤ 𝟏𝑛 ,
(1)
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 Σ𝑡 as each new observation 𝐱𝑡 arrives, using 𝑂(𝑝2 ) working memory and 𝑂(𝑝2 ) work per update. Three strategies address this problem, but a unified treatment covering algebraic equivalence, floatingpoint error analysis, and statistical uncertainty quantification has not previously appeared. We fill
2
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
this gap. The Gram (bariance) algorithm. It is generally well known (e.g. Reichel [8]) that the sample variance has an 𝑂(𝑛) computation via scalar sums, without explicit centering. Lifting this to the matrix case yields (1), with the streaming update rule 𝐆𝑡 ← 𝐆𝑡−1 + 𝐱𝑡 𝐱𝑡⊤ and 𝐬𝑡 ← 𝐬𝑡−1 + 𝐱𝑡 . The on-demand formula Σ̂𝑡 = (𝑡𝐆𝑡 − 𝐬𝑡 𝐬⊤ 𝑡 )∕[𝑡(𝑡 − 1)] executes in one level-3 BLAS call. Its weakness is catastrophic cancellation when the data mean is large relative to the variance: both 𝑡𝐆𝑡 and 𝐬𝑡 𝐬⊤ 𝑡 are 2 2 𝑂(𝑡 ‖𝐱‖ ̄ ) while their difference is 𝑂(𝑡‖Σ‖). The Welford algorithm. Welford [10] proposed a shift-invariant one-pass recurrence for scalar variance, extended to covariance by Chan, Golub & LeVeque [2]. The algorithm maintains a running mean µ𝑡 and a matrix 𝐌𝑡 of centred outer-product corrections, updating both in 𝑂(𝑝2 ) per step. Because corrections are computed relative to the current mean, the method is immune to large constant shifts in the data. The Chan–Golub–LeVeque (CGL) algorithm. Chan, Golub & LeVeque [2] introduced a binary merge formula for independently computed covariance summaries: given (𝑛𝐴 , µ𝐴 , 𝐌𝐴 ) and (𝑛𝐵 , µ𝐵 , 𝐌𝐵 ), the combined summary is computed exactly. This makes the algorithm tree-parallelisable and directly applicable to federated settings where nodes cannot share raw observations. Conformal prediction for streaming covariance. All three algorithms produce point estimates. To quantify uncertainty, we propose a split conformal prediction protocol [1, 9] that constructs a finitesample, distribution-free confidence interval for each entry Σ𝑘𝑙 at every step 𝑡 ≥ 2 of the stream. The interval has guaranteed marginal coverage 1 − 𝛼 for any 𝛼 ∈ (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. Contributions. 1. Complete algorithmic descriptions with proofs of correctness and three-way algebraic equivalence (Sections 3–4). 2. A floating-point error analysis quantifying rounding accumulation and catastrophic cancellation for each algorithm (Section 5). 3. A split conformal prediction framework for streaming covariance entry estimation, with a finite-sample coverage guarantee (Section 6). 4. Benchmarks on x86-64 hardware with OpenBLAS covering runtime (varying 𝑛 and 𝑝), accuracy under Gaussian, heavy-tailed, ill-conditioned, and near-singular data, and conformal coverage under large shifts (Sections 8–9). 5. A practitioner decision guide (Section 11). Complete proofs for all results appear in Appendix 13.
3
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
2 Notation and Setup Let 𝐱1 , 𝐱2 , … be a stream of observations in ℝ𝑝 , 𝑝 ≥ 2. After 𝑡 ≥ 2 steps the target is the unbiased sample covariance 𝑡 𝑡 ∑ 1 ∑ ⊤ −1 Σ𝑡 = (𝐱 − 𝐱̄ 𝑡 )(𝐱𝑖 − 𝐱̄ 𝑡 ) , 𝐱̄ 𝑡 = 𝑡 𝐱𝑖 . (2) 𝑡 − 1 𝑖=1 𝑖 𝑖=1 ∑𝑡 ∑𝑡 We write 𝐬𝑡 = 𝑖=1 𝐱𝑖 , 𝐆𝑡 = 𝑖=1 𝐱𝑖 𝐱⊤ , µ𝑡 = 𝐬𝑡 ∕𝑡, and 𝐌𝑡 for the Welford/CGL correction matrix. 𝑖 The unit roundoff for IEEE 754 double precision is 𝜀mach = 2−53 ≈ 1.11 × 10−16 . We write ‖ ⋅ ‖𝐹 , ‖ ⋅ ‖2 , ‖ ⋅ ‖max for the Frobenius, spectral, and entry-wise maximum norms. For a matrix 𝐴 the condition number is 𝜅(𝐴) = ‖𝐴‖2 ‖𝐴−1 ‖2 .
3 Algorithms 3.1 Gram (Bariance) Algorithm Algorithm 1 Gram streaming covariance Require: Stream 𝐱1 , 𝐱2 , … ∈ ℝ𝑝 1: 𝐬 ← 𝟎𝑝 ; 𝐆 ← 𝟎𝑝×𝑝 ; 𝑡 ← 0 2: for each new observation 𝐱 do 3: 𝑡 ← 𝑡 + 1; 𝐬 ← 𝐬 + 𝐱; 𝐆 ← 𝐆 + 𝐱𝐱⊤ 4: if 𝑡 ≥ 2 then 5: return Σ̂Gram = (𝑡𝐆 − 𝐬𝐬⊤ )∕[𝑡(𝑡 − 1)] 𝑡 6: end if 7: end for Theorem 3.1 (Gram identity). For any 𝐱1 , … , 𝐱𝑡 ∈ ℝ𝑝 with 𝑡 ≥ 2, Σ𝑡 =
𝑡𝐆𝑡 − 𝐬𝑡 𝐬⊤ 𝑡 𝑡(𝑡 − 1)
.
Proof. See Appendix 13.1. BLAS structure. Each update adds a rank-one correction to 𝐆 (DSYR, level-2 BLAS) and a vector increment to 𝐬 (DAXPY). Computing Σ̂𝑡 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.
4
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
3.2 Welford Algorithm Algorithm 2 Welford streaming covariance Require: Stream 𝐱1 , 𝐱2 , … ∈ ℝ𝑝 1: µ ← 𝟎𝑝 ; 𝐌 ← 𝟎𝑝×𝑝 ; 𝑡 ← 0 2: for each new observation 𝐱 do 3: 𝑡 ←𝑡+1 4: ∆←𝐱−µ 5: µ ← µ + ∆∕𝑡 6: 𝐌 ← 𝐌 + ∆ (𝐱 − µ)⊤ 7: if 𝑡 ≥ 2 then return Σ̂Welf = 𝐌∕(𝑡 − 1) 𝑡 8: end if 9: end for
⊳ residual w.r.t. old mean ⊳ update mean in place ⊳ outer product
Theorem 3.2 (Welford invariant). After processing 𝐱1 , … , 𝐱𝑡 via Algorithm 2, 𝐌𝑡 =
𝑡 ∑
(𝐱𝑖 − µ𝑡 )(𝐱𝑖 − µ𝑡 )⊤ .
𝑖=1
Hence Σ̂Welf = Σ𝑡 . 𝑡 Proof. See Appendix 13.2. The critical property is that both ∆ = 𝐱𝑡 − µ𝑡−1 and 𝐱𝑡 − µ𝑡 are computed as residuals from means of the same scale as the data, preventing cancellation regardless of the absolute value of µ𝑡 .
3.3 Chan–Golub–LeVeque Algorithm Algorithm 3 CGL merge of two summaries Require: (𝑛𝐴 , µ𝐴 , 𝐌𝐴 ) and (𝑛𝐵 , µ𝐵 , 𝐌𝐵 ) 1: ∆ ← µ𝐵 − µ𝐴 ; 𝑛𝐴𝐵 ← 𝑛𝐴 + 𝑛𝐵 2: µ𝐴𝐵 ← (𝑛𝐴 µ𝐴 + 𝑛𝐵 µ𝐵 )∕𝑛𝐴𝐵 3: 𝐌𝐴𝐵 ← 𝐌𝐴 + 𝐌𝐵 + ∆∆⊤ ⋅ 𝑛𝐴 𝑛𝐵 ∕𝑛𝐴𝐵 4: return (𝑛𝐴𝐵 , µ𝐴𝐵 , 𝐌𝐴𝐵 )
Theorem 3.3 (CGL correctness). Let 𝐴, 𝐵 be disjoint index sets. Define ∆ = µ𝐵 − µ𝐴 . Then 𝐌𝐴∪𝐵 = 𝐌𝐴 + 𝐌𝐵 +
𝑛𝐴 𝑛𝐵 ∆∆⊤ . 𝑛𝐴 + 𝑛𝐵
Proof. See Appendix 13.3. Setting block size 𝑏 = 1 (each block is a single observation) recovers Algorithm 2 exactly; block size 𝑏 = 𝑛 recovers the batch formula (1). Because the merge operation is associative and commutative, the algorithm can be applied in any binary-tree order, enabling map-reduce computation.
5
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
4 Algebraic Equivalence Theorem 4.1 (Three-way equivalence in exact arithmetic). In exact arithmetic, Σ̂Gram = Σ̂Welf = 𝑡 𝑡 CGL = Σ𝑡 for all 𝑡 ≥ 2 and any block size 𝑏 ≥ 1 in the CGL algorithm. Σ̂𝑡 Proof. Gram equals Σ𝑡 by Theorem 3.1. Welford equals Σ𝑡 by Theorem 3.2. CGL with 𝑏 = 1 equals Welford by direct substitution: one observation 𝐱𝑡 is a block with 𝑛𝐵 = 1, 𝐌𝐵 = 𝟎, µ𝐵 = 𝐱𝑡 , and the merge formula reduces to the Welford update. CGL with 𝑏 > 1 follows by induction on the tree using Theorem 3.3. Remark 4.2 (Bariance connection). For 𝑝 = 1, Theorem 3.1 gives Σ̂𝑡 = (𝑡𝑆𝑥𝑥 −𝑆𝑥2 )∕[𝑡(𝑡 −1)], which is the bariance identity of Reichel [8]. The Gram algorithm is the natural matrix extension of the bariance.
5 Floating-Point Error Analysis All three algorithms are algebraically identical but differ in their numerical stability. We quantify two distinct phenomena: rounding accumulation over 𝑡 steps, and catastrophic cancellation induced by a non-zero data mean.
5.1 Rounding Accumulation 1∕2
Write 𝑥̄ = ‖µ𝑡 ‖∞ and 𝜎 = ‖Σ𝑡 ‖∞ . Proposition 5.1 (Gram rounding bound). Let |𝑥𝑖𝑘 | ≤ 𝑋 for all 𝑖, 𝑘. The entry-wise error of the Gram estimator satisfies 𝑥̄ 2 | || Gram 𝜀 . ||Σ̂𝑘𝑙,𝑡 − Σ𝑘𝑙,𝑡 ||| ≲ 𝑋 2 𝜀mach + 𝑡 − 1 mach Proposition 5.2 (Welford rounding bound). Under the same assumptions, | || Welf ||Σ̂𝑘𝑙,𝑡 − Σ𝑘𝑙,𝑡 ||| ≲ 𝜎𝑘 𝜎𝑙 𝜀mach , ̄ where 𝜎𝑘2 = Σ𝑘𝑘,𝑡 . independently of 𝑥, Proposition 5.3 (CGL rounding bound). For a balanced binary-tree merge of depth log2 𝑡, || CGL | ||Σ̂𝑘𝑙,𝑡 − Σ𝑘𝑙,𝑡 ||| ≲ 𝜎𝑘 𝜎𝑙 𝜀mach log2 𝑡. Proofs are given in Appendix 13.4. The key contrast: Gram accumulates error at scale 𝑋 2 (the magnitude of raw observations), while Welford accumulates at scale 𝜎𝑘 𝜎𝑙 (the covariance of centred observations). The Welford bound is thus independent of the data mean, as formalised below.
6
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
5.2 Catastrophic Cancellation Proposition 5.4 (Cancellation bound). Suppose 𝐱𝑖 = µ + 𝐳𝑖 with ‖µ‖2 = 𝑐 and 𝐳𝑖 ∼ (0, Σ), ‖Σ‖2 = 𝜎2 . Then ‖‖ Gram ‖ ‖‖Σ̂𝑡 − Σ𝑡 ‖‖‖ ≳ 𝑝𝑐2 𝜀mach , ‖ ‖𝐹 ‖‖ Welf ‖ ‖‖Σ̂𝑡 − Σ𝑡 ‖‖‖ = 𝑂(𝑝𝜎2 𝜀mach ), ‖ ‖𝐹 independently of 𝑐. Proof. See Appendix 13.5. For 𝑐 = 107 and 𝜎 = 1, the Gram bound gives 𝑝𝑐2 𝜀mach ≈ 𝑝 ⋅ 1014 ⋅ 10−16 = 𝑝 ⋅ 10−2 , which becomes non-negligible relative to 𝜎2 = 1. At 𝑐 = 1012 , the error is 𝑂(𝑝 ⋅ 108 ), fully destroying the estimate (confirmed experimentally in Figure 7).
5.3 Summary of Bounds Table 1: Floating-point error bounds. 𝑐 = ‖µ‖2 , 𝜎2 = ‖Σ‖2 , 𝜀mach = 2−53 . “Shift-invariant” means the bound does not grow with 𝑐.
Algorithm
Error bound
Scale
Shift-inv.
Parallel
Gram Welford CGL
𝑂(𝑐2 𝜀mach + 𝜎2 𝜀mach ) 𝑂(𝜎2 𝜀mach ) 𝑂(𝜎2 𝜀mach log 𝑡)
𝑋2 𝜎2 𝜎2 log 𝑡
✓ ✓
✓
6 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 Σ𝑘𝑙 at every step 𝑡 of the stream.
6.1 Background: Split Conformal Prediction Split conformal prediction [1, 7, 9] produces a (1 − 𝛼)-coverage prediction interval for a new observation using a held-out calibration set. Given calibration nonconformity scores 𝑠1 , … , 𝑠𝑚 , the conformal quantile is ⌈(𝑚 + 1)(1 − 𝛼)⌉ ), 𝑞̂𝛼+ = Quantile({𝑠𝑖 }𝑚 ; 𝑖=1 𝑚 and the interval for a fresh test point has guaranteed marginal coverage Pr(true value ∈ 𝐶̂ 𝛼 ) ≥ 1 − 𝛼, with no distributional assumptions on 𝑠1 , … , 𝑠𝑚 beyond exchangeability.
6.2 Protocol for Streaming Covariance Fix target entry (𝑘, 𝑙), nominal coverage 1 − 𝛼, and an algorithm 𝒜 ∈ {Gram, Welford, CGL}. Let Σ𝑘𝑙 denote the true population covariance entry. 7
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
(𝑗)
1. Calibration. Draw 𝑚 independent calibration trajectories {𝐱𝑖 }𝑡𝑖=1 , 𝑗 = 1, … , 𝑚, from the data-generating distribution. For each trajectory 𝑗 and each step 𝑡 of interest, compute (𝑗) ) | ( (𝑗) | 𝑠𝑗 (𝑡) = |||𝒜 𝐱1 , … , 𝐱𝑡 − Σ𝑘𝑙 |||. 𝑘𝑙
2. Quantile. Set 𝑞̂𝛼+ (𝑡) = Quantile({𝑠𝑗 (𝑡)}; ⌈(𝑚 + 1)(1 − 𝛼)⌉∕𝑚). 3. Interval. For a fresh test stream at step 𝑡, let 𝜎̂ 𝑘𝑙 (𝑡) be the algorithm’s estimate. The conformal interval is 𝐶̂ 𝛼 (𝑡) = [𝜎̂ 𝑘𝑙 (𝑡) − 𝑞̂𝛼+ (𝑡), 𝜎̂ 𝑘𝑙 (𝑡) + 𝑞̂𝛼+ (𝑡)]. Theorem 6.1 (Finite-sample coverage). Under the assumption that calibration and test trajectories are i.i.d. (exchangeable), the conformal interval 𝐶̂ 𝛼 (𝑡) satisfies ( ) Pr Σ𝑘𝑙 ∈ 𝐶̂ 𝛼 (𝑡) ≥ 1 − 𝛼. Proof. See Appendix 13.6. The proof applies the standard split-conformal coverage argument of Vovk, Gammerman & Shafer [9] with nonconformity score 𝑠 = |𝜎̂ 𝑘𝑙 (𝑡) − Σ𝑘𝑙 |. Remark 6.2 (Algorithm dependence of interval width). The coverage guarantee is algorithm-independent. However, the interval width 2𝑞̂𝛼+ (𝑡) reflects the variance of each estimator’s error distribution. For Gram under large shift 𝑐, the calibration scores 𝑠𝑗 (𝑡) are inflated by the cancellation term 𝑂(𝑐2 𝜀mach ), producing wide intervals or, if 𝑐 is not replicated in calibration, invalid coverage. Welford is immune to this effect. Remark 6.3 (Unknown Σ𝑘𝑙 ). In practice Σ𝑘𝑙 is unknown, so calibration scores are computed using the batch reference on each calibration trajectory (which uses the full 𝑚 observations as ground truth). For large 𝑚 this reference converges to Σ𝑘𝑙 at rate 𝑂(𝑚−1∕2 ). Alternatively, one can treat the batch reference as the target and the interval inherits the same finite-sample guarantee relative to that target.
7 Computational Cost Table 2: Per-observation cost in streaming mode. 𝑝 is the dimension. “Level” refers to the BLAS level of the dominant operation.
Algorithm
FLOPs per update
Memory (doubles)
BLAS level
Gram Welford CGL (𝑏 = 1) CGL (𝑏 > 1)
𝑝2 + 𝑝 (DSYR + DAXPY) 2𝑝2 + 2𝑝 (2×DGER) Same as Welford 2𝑝2 𝑛𝑏 ∕𝑏 + 𝑝2 per block
𝑝2 + 𝑝 𝑝2 + 𝑝 𝑝2 + 𝑝 + 𝑏𝑝 𝑝2 + 𝑝 + 𝑏𝑝
2 2 2 3 (DSYRK)
Gram batch numpy.cov
2𝑛𝑝2 (DSYRK) + 𝑝2 2𝑛𝑝2 + 2𝑛𝑝 (center + DSYRK)
𝑛𝑝 2𝑛𝑝
3 3
The Gram algorithm is faster in batch mode than numpy.cov because it avoids forming and writing the centred matrix 𝐗 − 𝟏𝑛 µ⊤ , saving 𝑂(𝑛𝑝) memory writes. In streaming mode, Welford computes two outer products per step (for ∆ and 𝐱𝑡 − µ𝑡 ) versus one for Gram, yielding a ≈ 2× FLOP disadvantage; this gap is visible in practice when the outer product (DGER) is the bottleneck. 8
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
8 Experimental Protocol 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. Runtime measurement. For each (𝑛, 𝑝) pair: (i) generate 𝐗 ∼ 𝒩(0, 𝐼); (ii) run 3 warm-up calls; (iii) record 20–30 wall-clock times via time.perf_counter; (iv) remove outliers using the 1.5 × IQR rule; (v) report the trimmed mean with a 95% bootstrap percentile interval (300 resamples). Accuracy measurement. For each data configuration the reference is numpy.cov. We report the entry-wise maximum error ‖Σ̂ − 𝑆ref ‖max and the relative Frobenius error ‖Σ̂ − 𝑆ref ‖𝐹 ∕‖𝑆ref ‖𝐹 . Data configurations. We test four regimes: (a) Gaussian i.i.d., varying 𝑛 and 𝑝; (b) Student-𝑡3 (heavy-tailed), varying 𝑛; (c) prescribing a range of condition numbers 𝜅(𝐗) ∈ [10, 1014 ]; (d) nearsingular data with smallest singular value 𝜎min ∈ [10−12 , 1]. Conformal experiments. Calibration uses 𝑚 = 600 independent trajectories drawn from a 𝑝 = 5 multivariate distribution with Toeplitz covariance (entry (𝑖, 𝑗): 0.5|𝑖−𝑗| ). Test coverage is estimated from 1200 fresh trajectories. The shift experiment varies the mean offset 𝑐 ∈ {0, 103 , 106 , 109 , 1012 } on a 𝑝 = 4 identity-covariance distribution and evaluates coverage and interval width at 𝑡 = 150.
9 Results 9.1 Runtime Figure 1 shows wall-clock time versus 𝑛 for 𝑝 = 10 and 𝑝 = 50. Gram is the fastest batch method, consistently beating numpy.cov by a factor of 1.3–1.6× due to the avoided 𝑛 × 𝑝 centering write. Welford is slowest by 50–80× in pure-Python form because the inner loop invokes numpy.outer once per row, incurring 𝑂(𝑛) Python-level calls; a compiled implementation would close this to ≈ 2× (the FLOP ratio from Table 2). CGL (recursive, block size 64) sits between the two. Figure 2 shows runtime versus dimension at fixed 𝑛 = 4,000. For large 𝑝 all batch methods are dominated by the 𝑂(𝑛𝑝2 ) matrix multiply and their curves converge; Gram retains a small constantfactor advantage from avoiding the centred copy. Figure 3 shows the speed ratio of each method relative to Gram.
9.2 Numerical Accuracy: Gaussian Data Figure 4 shows accuracy on Gaussian i.i.d. data. All methods agree with numpy.cov to within 10−13 in both metrics, consistent with Propositions 5.1–5.3.
9
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
A
Gram Welford Chan-Golub-LeVeque NumPy
20
15
10
5
0
Covariance construction, Dimension p = 50 Gram Welford Chan-Golub-LeVeque NumPy
40
Wall-clock time, milliseconds
25
Wall-clock time, milliseconds
B
Covariance construction, Dimension p = 10
30
20
10
0 0
2,000
4,000
6,000
8,000
10,000
0
2,000
4,000
Sample size n
6,000
8,000
10,000
Sample size n
Note: Lines report trimmed mean runtime. Shaded bands are bootstrap 95 percent confidence intervals.
Figure 1: Wall-clock time vs. sample size 𝑛. Each point is the trimmed mean over 20–30 repetitions after warm-up and 1.5 × IQR outlier removal; shaded bands are 95% bootstrap percentile intervals. Left: 𝑝 = 10. Right: 𝑝 = 50. Welford is omitted from the right panel because it is ≈ 60× slower and would dominate the scale.
Runtime as covariance dimension increases Gram Welford Chan-Golub-LeVeque NumPy
Wall-clock time, milliseconds
120
100
80
60
40
20
0 0
25
50
75
100
125
150
175
200
Dimension p, fixed sample size n = 4,000 Note: Lines report trimmed mean runtime. Shaded bands are bootstrap 95 percent confidence intervals.
Figure 2: Wall-clock time vs. dimension 𝑝 at 𝑛 = 4,000. Asymptotically, all batch methods are 𝑂(𝑛𝑝2 ) and track each other; the constant-factor gap reflects memory-bandwidth differences.
9.3 Heavy-Tailed and Ill-Conditioned Data Figure 5 tests Student-𝑡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. 10
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
Relative runtime benchmark
50
Welford divided by Gram Chan-G-L divided by Gram NumPy divided by Gram Parity with Gram
Runtime ratio
40
30
20
10
0 0
2,000
4,000
6,000
8,000
10,000
Sample size n, fixed dimension p = 50 Note: Values above one indicate that the comparison method is slower than Gram. Shaded bands use endpoint ratios from bootstrap confidence intervals.
Figure 3: Speed ratio vs. Gram at 𝑝 = 50. Values > 1 indicate the alternative is slower. numpy.cov is 1.4–1.6× slower; CGL (recursive, block 64) is 3–8× slower due to merge overhead. A
B
Gaussian accuracy, maximum absolute error
Gaussian accuracy, relative frobenius error Gram Welford Chan-Golub-LeVeque
Relative Frobenius error
Maximum absolute error
Gram Welford Chan-Golub-LeVeque
10−15
0
2,000
4,000
6,000
8,000
10−15
10,000
Sample size n, Gaussian data with p = 10
0
2,000
4,000
6,000
8,000
10,000
Sample size n, Gaussian data with p = 10
Note: Errors are computed relative to numpy.cov. Shaded bands are 95 percent Monte Carlo confidence intervals across replications.
Figure 4: Accuracy vs. 𝑛 on Gaussian i.i.d. data, 𝑝 = 10. Errors are at floatingpoint noise level (< 10−13 ) for all methods and all 𝑛 tested, confirming algebraic equivalence to numerical precision.
Figure 6 sweeps the condition number of 𝐗. All methods degrade once 𝜅 ≳ 107 , as predicted by standard floating-point theory; Welford and CGL are modestly more robust in the range 𝜅 ∈ [106 , 1012 ].
11
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
A
Heavy-tailed accuracy, relative frobenius error
Gram Welford Chan-Golub-LeVeque
−14
Gram Welford Chan-Golub-LeVeque
Relative Frobenius error
Maximum absolute error
10
B
Heavy-tailed accuracy, maximum absolute error
10−15
0
2,000
4,000
6,000
8,000
10−15
10,000
0
Sample size n, Student t data with 3 degrees of freedom and p = 8
2,000
4,000
6,000
8,000
10,000
Sample size n, Student t data with 3 degrees of freedom and p = 8
Note: Data are drawn from a Student t distribution with 3 degrees of freedom. Shaded bands are 95 percent Monte Carlo confidence intervals.
Figure 5: Accuracy on Student-𝑡3 data, 𝑝 = 8. Heavy tails increase absolute errors slightly (larger entries, more cancellation potential in Gram) but all methods remain numerically consistent at noise level. A
Maximum absolute error
106
B
Ill-conditioned data, maximum absolute error Gram Welford Chan-Golub-LeVeque
Ill-conditioned data, relative frobenius error Gram Welford Chan-Golub-LeVeque
10−15
Relative Frobenius error
10
10
102 10−2 10−6 10−10
6 × 10−16
4 × 10−16
3 × 10−16 10−14
102
104
106
108
1010
1012
1014
102
Condition number of the data matrix
104
106
108
1010
1012
1014
Condition number of the data matrix
Note: Synthetic matrices are constructed with prescribed singular values. Shaded bands are 95 percent Monte Carlo confidence intervals.
Figure 6: Accuracy vs. condition number 𝜅(𝐗), 𝑛 = 500, 𝑝 = 6. All methods degrade for 𝜅 ≳ 107 . Welford and CGL maintain lower error than Gram in the moderate ill-conditioning regime.
9.4 Catastrophic Cancellation Under Large Shift Figure 7 shifts all observations by a constant 𝑐 and measures error against the unshifted (true) covariance. Gram error grows as 𝑐2 𝜀mach (Proposition 5.4), losing ≈ 9–10 decimal digits by 𝑐 = 1012 . Welford maintains full double-precision accuracy throughout. CGL loses moderate accuracy, reflecting the inter-block cancellation in the merge formula.
9.5 Online Fidelity batch Figure 8 tracks ‖Σ̂𝒜 ‖max as 𝑡 grows on a zero-mean stream. All methods converge to floating𝑡 − Σ̂𝑡 √ point noise by 𝑡 ≈ 100; the early fluctuations are consistent with the 1∕ 𝑡 rate of the estimation error.
12
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
Maximum absolute error relative to centered reference
Cancellation under large additive shifts 10
11
Gram Welford Chan-Golub-LeVeque
107 103 10−1 10−5 10−9 10−13
101
103
105
107
109
1011
1013
Common additive shift Note: The same data matrix is shifted by a scalar constant before covariance is computed. Shaded bands are 95 percent Monte Carlo confidence intervals.
Figure 7: Catastrophic cancellation under large shift. Gram error grows as 𝑐2 𝜀mach , while Welford remains accurate to 𝑂(𝜎2 𝜀mach ) throughout. CGL sits between the two.
9.6 Conformal Prediction Results Coverage and width under well-conditioned data. Figure 9 (left) shows empirical coverage versus 𝑡 for all three algorithms at nominal level 1 − 𝛼 = 95%. All algorithms achieve the nominal level (dashed line) at every 𝑡 ≥ 10, confirming Theorem 6.1. Coverage slightly exceeds 95% for small 𝑡 because the conformal quantile is conservative for small calibration sets. Figure 9 (right) shows interval width: all three methods produce intervals of the same order, narrowing as 𝑡 grows. Gram produces marginally narrower intervals for well-conditioned data because its lower variance translates to tighter calibration scores. Coverage under large mean shift. Figure 10 tests conformal coverage when calibration and test data share the same shift 𝑐 (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 𝑐, reflecting the 𝑂(𝑐2 𝜀mach ) error term inflating the calibration scores. Welford and CGL maintain narrow intervals throughout.
13
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
Maximum absolute error relative to batch reference
Online covariance fidelity over the stream Gram Welford Chan-Golub-LeVeque
10−15
10−16 0
1,000
2,000
3,000
4,000
5,000
6,000
Observations processed Note: At each checkpoint, the streaming estimate is compared with batch numpy.cov on all observations seen so far. Shaded bands are 95 percent Monte Carlo confidence intervals.
Figure 8: Online covariance fidelity, 𝑝 = 8. All methods agree with the batch reference to floating-point noise by 𝑡 ≈ 100. A
B
Conformal coverage over the stream
Conformal interval width over the stream
1.000
Gram Welford Chan-Golub-LeVeque
1.4 0.975 1.2
Interval width
Empirical coverage
0.950 0.925 0.900 0.875
0.8
0.6
0.850
Gram Welford Chan-Golub-LeVeque 95 percent nominal
0.825 0.800
1.0
25
50
75
100
125
150
175
0.4
200
Stream length
25
50
75
100
125
150
175
200
Stream length
Note: Coverage bands are plus or minus 1.96 standard errors. The data-generating covariance is Toeplitz with correlation parameter 0.5.
Figure 9: Conformal coverage and width vs. stream length 𝑡. Left: all algorithms achieve the 95% nominal level (dashed). Shaded bands are ±1.96 standard errors of the empirical coverage. Right: interval width narrows as 𝑡 increases. Calibration: 𝑚 = 600 trajectories, Toeplitz Σ (𝑝 = 5).
10 Applications Federated / privacy-restricted learning. The Gram sufficient statistics (𝐆𝑡 , 𝐬𝑡 , 𝑡) 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 14
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
A
B
Conformal coverage under shifts 3.0
1.0
Gram Welford Chan-Golub-LeVeque
2.5
Interval width
0.8
Empirical coverage
Conformal interval width under shifts
1e9
0.6
0.4
2.0
1.5
1.0 95 percent nominal Gram Welford Chan-Golub-LeVeque
0.2
0.0
0
1e2
1e3
0.5
1e4
1e6
1e8
1e10
0.0
1e12
0
1e2
1e3
Common additive shift
1e4
1e6
1e8
1e10
1e12
Common additive shift
Note: Intervals are calibrated separately at each shift level using split conformal scores. Coverage bars include 95 percent binomial standard-error intervals.
Figure 10: Conformal coverage (left) and interval width (right) under large shift 𝑐, 𝑝 = 4, 𝑡 = 150. Coverage is maintained for all algorithms because calibration and test share the same shift (exchangeability holds). Interval width inflates exponentially for Gram as 𝑐 increases, reflecting Gram’s 𝑂(𝑐2 𝜀mach ) numerical error. Welford and CGL maintain narrow intervals.
the parallelism advantage and the numerical stability of Welford. Sandwich covariance for M-estimators. Let 𝐠𝑖 ∈ ℝ𝑝 be the score vector for observation 𝑖. The ̂ = 𝑛∕(𝑛 − 1) ⋅ (𝑛𝐆 − 𝐬𝐬⊤ )∕𝑛2 where heteroskedasticity-consistent sandwich covariance [6, 11] is Ω ∑ ∑ and 𝐬 = 𝑖 𝐠𝑖 . This is precisely the Gram formula applied to the score matrix, admitting 𝐆 = 𝑖 𝐠𝑖 𝐠⊤ 𝑖 a streaming update as each observation is processed. Panel / fixed-effects data. For panel unit 𝑖 observed at 𝑇 periods, the within-unit covariance is ̂ 𝑖 = (𝐗⊤ 𝐗𝑖 − 𝑇 −1 𝐬𝑖 𝐬⊤ )∕(𝑇 − 1). The streaming formulation processes units sequentially without Σ 𝑖 𝑖 forming the full 𝑁𝑇 × 𝑝 matrix. Conformal uncertainty in online inference. The conformal intervals of Section 6 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.
11 Practitioner Decision Guide 1∕2
Table 3: Algorithm selection. “Large shift”: ‖µ‖ ≫ ‖Σ‖2 .
Scenario Batch, BLAS available Streaming, zero/small mean Streaming, large shift Conformal interval width Distributed / federated Privacy: no raw-data access
Gram
Welford
CGL (𝑏=1)
CGL (𝑏 ≫ 1)
Best Best Avoid Avoid (large 𝑐) Good (agg.) Only
Slow Good Best Best Poor No
Medium Good Good Good Poor Partial
Good — — — Best Best
15
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
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 6 for principled uncertainty quantification.
12 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× 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 Higham [4] and Chan, Golub & LeVeque [2].
References [1] Angelopoulos, A. N. and Bates, S. (2023). Conformal prediction: A gentle introduction. Foundations and Trends in Machine Learning, 16(4):494–591. [2] Chan, T. F., Golub, G. H., and LeVeque, R. J. (1979). Updating formulae and a pairwise algorithm for computing sample variances. Technical Report STAN-CS-79-773, Stanford University. [3] Golub, G. H. and Van Loan, C. F. (2013). Matrix Computations, 4th ed. Johns Hopkins University Press. [4] Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms, 2nd ed. SIAM. [5] Lehmann, E. L. and Casella, G. (1998). Theory of Point Estimation, 2nd ed. Springer. [6] Newey, W. K. and West, K. D. (1987). A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica, 55(3):703–708. [7] Papadopoulos, H., Proedrou, K., Vovk, V., and Gammerman, A. (2002). Inductive confidence machines for regression. In Proc. ECML, pp. 345–356. [8] Reichel, F. (2025). High-Performance Variance-Covariance Matrix Construction Using an Uncentered Gram Formulation, Forthcoming in European Journal of Mathematics and Statistics (EJ-MATH).
16
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
[9] Vovk, V., Gammerman, A., and Shafer, G. (2005). Algorithmic Learning in a Random World. Springer. [10] Welford, B. P. (1962). Note on a method for calculating corrected sums of squares and products. Technometrics, 4(3):419–420. [11] White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4):817–838.
13 Ommited Proofs 13.1 Proof of Theorem 3.1: Gram Identity Proof. Expand the (𝑘, 𝑙) entry of the left side of (2) directly: 𝑡 ∑
(𝑥𝑖𝑘 − 𝑥̄ 𝑘 )(𝑥𝑖𝑙 − 𝑥̄ 𝑙 ) =
𝑖=1
𝑡 ∑
𝑥𝑖𝑘 𝑥𝑖𝑙 − 𝑥̄ 𝑘
𝑖=1
𝑡 ∑
𝑥𝑖𝑙 − 𝑥̄ 𝑙
𝑖=1
𝑡 ∑
𝑥𝑖𝑘 + 𝑡 𝑥̄ 𝑘 𝑥̄ 𝑙
𝑖=1
= 𝐺𝑘𝑙 − 𝑥̄ 𝑘 𝑠𝑙 − 𝑥̄ 𝑙 𝑠𝑘 + 𝑡𝑥̄ 𝑘 𝑥̄ 𝑙 𝑠𝑘 𝑠𝑙 𝑠𝑙 𝑠𝑘 𝑠𝑘 𝑠𝑙 = 𝐺𝑘𝑙 − − + 𝑡 𝑡 𝑡 𝑠𝑘 𝑠𝑙 = 𝐺𝑘𝑙 − , 𝑡 where we used 𝑥̄ 𝑘 = 𝑠𝑘 ∕𝑡. Dividing by 𝑡 − 1 and multiplying numerator and denominator by 𝑡 gives 𝐺𝑘𝑙 − 𝑠𝑘 𝑠𝑙 ∕𝑡 𝑡𝐺𝑘𝑙 − 𝑠𝑘 𝑠𝑙 = , 𝑡−1 𝑡(𝑡 − 1) which is the (𝑘, 𝑙) entry of (𝑡𝐆𝑡 − 𝐬𝑡 𝐬⊤ 𝑡 )∕[𝑡(𝑡 − 1)]. Since this holds for every (𝑘, 𝑙), the matrix identity follows. □
13.2 Proof of Theorem 3.2: Welford Invariant Proof. We proceed by induction on 𝑡. Base case (𝑡 = 1). At 𝑡 = 1: µ1 = 𝐱1 , ∆ = 𝐱1 − 𝟎 = 𝐱1 , and 𝐱1 − µ1 = 𝟎, so 𝐌1 = 𝟎. Thus, the ∑1 claimed sum 𝑖=1 (𝐱𝑖 − µ1 )(𝐱𝑖 − µ1 )⊤ = 𝟎. Inductive step. Assume 𝐌𝑡−1 = µ𝑡 = µ𝑡−1 + ∆𝑡 ∕𝑡.
∑𝑡−1 𝑖=1
(𝐱𝑖 − µ𝑡−1 )(𝐱𝑖 − µ𝑡−1 )⊤ . At step 𝑡, let ∆𝑡 = 𝐱𝑡 − µ𝑡−1 and
Then 𝐱𝑡 − µ𝑡 = ∆𝑡 − ∆𝑡 ∕𝑡 = (𝑡 − 1)∆𝑡 ∕𝑡. The update is 𝐌𝑡 = 𝐌𝑡−1 + ∆𝑡 (𝐱𝑡 − µ𝑡 )⊤ . We need to show 𝐌𝑡 =
17
∑𝑡 𝑖=1
(𝐱𝑖 − µ𝑡 )(𝐱𝑖 − µ𝑡 )⊤ .
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
Write 𝐱𝑖 − µ𝑡 = (𝐱𝑖 − µ𝑡−1 ) − (µ𝑡 − µ𝑡−1 ) = (𝐱𝑖 − µ𝑡−1 ) − ∆𝑡 ∕𝑡 for each 𝑖. Then 𝑡 ∑
𝑡−1
(𝐱𝑖 − µ𝑡 )(𝐱𝑖 − µ𝑡 )⊤ =
𝑖=1
∑
(𝐱𝑖 − µ𝑡 )(𝐱𝑖 − µ𝑡 )⊤ + (𝐱𝑡 − µ𝑡 )(𝐱𝑡 − µ𝑡 )⊤
𝑖=1 𝑡−1
=
∑[
(𝐱𝑖 − µ𝑡−1 ) −
∆𝑡 ][
∆ 𝑡 ]⊤
𝑡
𝑡
𝑖=1
Expanding the sum and using
∑𝑡−1 𝑖=1
∑
+
(𝑡−1)2 𝑡2
∆𝑡 ∆⊤ 𝑡 .
(𝐱𝑖 − µ𝑡−1 ) = 𝟎:
𝑡−1
=
(𝐱𝑖 − µ𝑡−1 ) −
(𝐱𝑖 − µ𝑡−1 )(𝐱𝑖 − µ𝑡−1 )⊤ +
𝑖=1
(𝑡 − 1)2 𝑡−1 ⊤ ∆ ∆ + ∆𝑡 ∆ ⊤ 𝑡 𝑡 𝑡 𝑡2 𝑡2
] (𝑡 − 1) [ 1 + (𝑡 − 1) ∆𝑡 ∆ ⊤ 𝑡 𝑡2 𝑡−1 = 𝐌𝑡−1 + ∆𝑡 ∆ ⊤ 𝑡 𝑡 ( 𝑡−1 )⊤ = 𝐌𝑡−1 + ∆𝑡 ∆𝑡 = 𝐌𝑡−1 + ∆𝑡 (𝐱𝑡 − µ𝑡 )⊤ = 𝐌𝑡 ,
= 𝐌𝑡−1 +
𝑡
completing the induction.
□
13.3 Proof of Theorem 3.3: CGL Correctness Proof. Write µ𝐴𝐵 = (𝑛𝐴 µ𝐴 + 𝑛𝐵 µ𝐵 )∕𝑛𝐴𝐵 . For 𝑖 ∈ 𝐴: 𝐱𝑖 − µ𝐴𝐵 = (𝐱𝑖 − µ𝐴 ) + (µ𝐴 − µ𝐴𝐵 ) = ∑ (𝐱𝑖 − µ𝐴 ) − 𝑛𝐵 ∆∕𝑛𝐴𝐵 . Summing outer products over 𝐴 and using 𝑖∈𝐴 (𝐱𝑖 − µ𝐴 ) = 𝟎: ∑
(𝐱𝑖 − µ𝐴𝐵 )(𝐱𝑖 − µ𝐴𝐵 )⊤ = 𝐌𝐴 +
𝑖∈𝐴
𝑛𝐴 𝑛𝐵2 2 𝑛𝐴𝐵
∆∆⊤ .
Similarly for 𝐵 (with µ𝐵 − µ𝐴𝐵 = 𝑛𝐴 ∆∕𝑛𝐴𝐵 ): ∑
(𝐱𝑖 − µ𝐴𝐵 )(𝐱𝑖 − µ𝐴𝐵
)⊤ = 𝐌
𝑖∈𝐵
𝐵+
2 𝑛𝐵 𝑛𝐴 2 𝑛𝐴𝐵
∆∆⊤ .
Adding: 𝐌𝐴∪𝐵 = 𝐌𝐴 + 𝐌𝐵 +
𝑛𝐴 𝑛𝐵 (𝑛𝐵 + 𝑛𝐴 ) 2 𝑛𝐴𝐵
∆∆⊤ = 𝐌𝐴 + 𝐌𝐵 +
𝑛𝐴 𝑛𝐵 ∆∆⊤ . 𝑛𝐴𝐵
□
13.4 Proofs of Floating-Point Bounds (Propositions 5.1–5.4) We use the standard model of floating-point arithmetic [4, Chapter 2]: for any operation ◦ ∈ {+, −, ×, ÷}, f l(𝑎◦𝑏) = (𝑎◦𝑏)(1 + 𝛿) with |𝛿| ≤ 𝜀mach . Sequential summation of 𝑛 terms bounded by 𝑀 satisfies 2 |𝑆̂𝑛 − 𝑆𝑛 | ≤ (𝑛 − 1)𝜀mach ⋅ 𝑛𝑀 + 𝑂(𝜀mach ) [4, Theorem 3.1]. ∑𝑡 Proof of Proposition 5.1. The (𝑘, 𝑙) entry of 𝐆̂ 𝑡 is computed as 𝐺̂ 𝑘𝑙 = 𝑖=1 𝑥𝑖𝑘 𝑥𝑖𝑙 in floating point. By the summation error bound each term 𝑥𝑖𝑘 𝑥𝑖𝑙 is bounded by 𝑋 2 , giving |𝐺̂ 𝑘𝑙 −𝐺𝑘𝑙 | ≤ (𝑡−1)𝜀mach ⋅𝑡𝑋 2 + 18
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
2 𝑂(𝜀mach ). The outer product 𝑠̂𝑘 𝑠̂𝑙 contributes an additional rounding at scale |𝑥̄ 𝑘 ||𝑥̄ 𝑙 |𝑡2 . Forming (𝑡𝐺̂ 𝑘𝑙 − 𝑠̂𝑘 𝑠̂𝑙 )∕[𝑡(𝑡 − 1)]:
|| ̂ Gram | 𝑡 |𝐺̂ 𝑘𝑙 − 𝐺𝑘𝑙 | + |𝑠̂𝑘 𝑠̂𝑙 − 𝑠𝑘 𝑠𝑙 | 2 + 𝑂(𝜀mach ) ||Σ𝑘𝑙 − Σ𝑘𝑙 ||| ≤ 𝑡(𝑡 − 1) (𝑡 − 1)𝜀mach 𝑡𝑋 2 + |𝑥̄ 𝑘 ||𝑥̄ 𝑙 |𝑡 2 𝜀mach 2 + 𝑂(𝜀mach ) 𝑡(𝑡 − 1) 𝑥̄ 𝑘 𝑥̄ 𝑙 ≲ 𝑋 2 𝜀mach + 𝜀 . □ 𝑡 − 1 mach ≤
Proof of Proposition 5.2. At each step 𝑡 the Welford algorithm computes ∆𝑡 = 𝐱𝑡 − µ̂ 𝑡−1 and 𝐱𝑡 − µ̂ 𝑡 ̂ 𝑡 (𝐱𝑡 − µ̂ 𝑡 )⊤ has entries of magnitude 𝑂(𝜎𝑘 𝜎𝑙 ) on average (in as residuals. Each outer product ∆ expectation over the trajectory). Sequential accumulation over 𝑡 − 1 steps gives 2 |Σ̂ Welf − Σ𝑘𝑙 | ≤ (𝑡 − 1)𝜀mach 𝜎𝑘 𝜎𝑙 + 𝑂(𝜀mach ), 𝑘𝑙
independent of 𝑥̄ 𝑘 or 𝑥̄ 𝑙 . See Chan, Golub & LeVeque [2], Theorem 3.1, for a detailed treatment. □ Proof of Proposition 5.3. A single merge call introduces rounding error |𝛿[𝐌𝐴𝐵 ]𝑘𝑙 | ≤ 3𝜀mach ⋅(𝑛𝐴 𝑛𝐵 ∕𝑛𝐴𝐵 )|∆𝑘 ||∆𝑙 |+ 2 𝑂(𝜀mach ) in the correction term. For a balanced binary tree of depth 𝑑 = log2 𝑡, there are 𝑡 − 1 total merges. At level 𝑗 (counting from leaves), there are 𝑡∕2𝑗 merges each with block sizes ≈ 2𝑗 , contributing 𝑂(𝜀mach 𝜎𝑘 𝜎𝑙 ) each. Summing over 𝑑 = log2 𝑡 levels: |Σ̂ CGL − Σ𝑘𝑙 | = 𝑂(𝜎𝑘 𝜎𝑙 𝜀mach log2 𝑡). 𝑘𝑙
□
Proof of Proposition 5.4. Suppose 𝐱𝑖 = µ + 𝐳𝑖 with 𝔼[𝐳𝑖 ] = 𝟎 and 𝔼[𝐳𝑖 𝐳𝑖⊤ ] = Σ. Then 𝐺𝑘𝑙 = 𝜇𝑘 𝜇𝑙 𝑡 + ∑ (𝑧 𝜇 + 𝑧𝑖𝑙 𝜇𝑘 + 𝑧𝑖𝑘 𝑧𝑖𝑙 ) and 𝑠𝑘 𝑠𝑙 ∕𝑡 = 𝜇𝑘 𝜇𝑙 𝑡 + (lower order). The dominant terms in 𝑡𝐺𝑘𝑙 and 𝑖 𝑖𝑘 𝑙 𝑠𝑘 𝑠𝑙 are both 𝑂(𝜇𝑘 𝜇𝑙 𝑡2 ), and their difference is 𝑂(𝑡Σ𝑘𝑙 ). In floating point, 𝐺̂ 𝑘𝑙 has rounding error 𝑂(𝑡𝑋 2 𝜀mach ) where 𝑋 ≈ |𝜇𝑘 | + 𝜎𝑘 . When |𝜇𝑘 | ≫ 𝜎𝑘 , 𝑋 ≈ |𝜇𝑘 |, and the error in Σ̂Gram is dominated 𝑘𝑙 2 2 2 2 by 𝑡|𝜇𝑘 | 𝜀mach ∕[𝑡(𝑡 − 1)] = |𝜇𝑘 | 𝜀mach ∕(𝑡 − 1) ≈ 𝜇𝑘 𝜀mach , which sums to 𝑂(𝑝𝑐 𝜀mach ) over the 𝑝 ‖ ‖ diagonal entries, giving ‖‖‖Σ̂Gram − Σ‖‖‖ ≳ 𝑝𝑐2 𝜀mach . ‖ ‖𝐹 For Welford, ∆𝑡 = 𝐳𝑡 + (µ − µ̂ 𝑡−1 ) is 𝑂(𝜎) once the running mean has converged (which happens geometrically fast). The outer product ∆𝑡 (𝐱𝑡 − µ̂ 𝑡 )⊤ is then 𝑂(𝜎2 ) regardless of 𝑐, giving ‖‖ Welf ‖ ‖‖Σ̂ − Σ‖‖‖ = 𝑂(𝑝𝜎2 𝜀mach ). □ ‖ ‖𝐹
13.5 Proof of Proposition 5.4: Cancellation Bound Proof. Write each observation as 𝐱𝑖 = µ + 𝐳𝑖 , where the centred component satisfies 𝔼[𝐳𝑖 ] = 𝟎 and 𝔼[𝐳𝑖 𝐳⊤ ] = Σ. For a fixed entry (𝑘, 𝑙), the Gram statistics are 𝑖 𝐺𝑘𝑙 =
𝑡 ∑ 𝑖=1
𝑥𝑖𝑘 𝑥𝑖𝑙 = 𝑡𝜇𝑘 𝜇𝑙 + 𝜇𝑘
𝑡 ∑ 𝑖=1
19
𝑧𝑖𝑙 + 𝜇𝑙
𝑡 ∑ 𝑖=1
𝑧𝑖𝑘 +
𝑡 ∑ 𝑖=1
𝑧𝑖𝑘 𝑧𝑖𝑙 ,
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
and 𝑠𝑘 𝑠𝑙 = (𝑡𝜇𝑘 +
𝑡 ∑
𝑧𝑖𝑘 ) (𝑡𝜇𝑙 +
𝑖=1
𝑡 ∑
𝑧𝑖𝑙 ) .
𝑖=1
Hence the two quantities entering the Gram numerator, 𝑡𝐺𝑘𝑙 and 𝑠𝑘 𝑠𝑙 , both contain the leading term 𝑡 2 𝜇𝑘 𝜇𝑙 . In exact arithmetic these leading terms cancel, leaving a centred quantity of order 𝑡Σ𝑘𝑙 . In floating-point arithmetic, however, the products and sums are formed at the raw data scale |𝜇𝑘 | + 𝜎𝑘 . Thus the rounding error in the Gram numerator contains a term of size ( ) 𝑂 𝑡2 |𝜇𝑘 𝜇𝑙 |𝜀mach , and division by 𝑡(𝑡 − 1) gives the entry-wise contribution || Gram | ||Σ̂𝑘𝑙,𝑡 − Σ𝑘𝑙,𝑡 ||| ≳ |𝜇𝑘 𝜇𝑙 |𝜀mach . Summing these contributions over the diagonal entries gives 𝑝
∑ 2 ‖‖ Gram ‖ ‖‖Σ̂𝑡 𝜇𝑘 𝜀mach . − Σ𝑡 ‖‖‖ ≳ ‖ ‖𝐹 𝑘=1
∑𝑝 Since 𝑘=1 𝜇𝑘2 = ‖µ‖22 = 𝑐2 , this is of order 𝑐2 𝜀mach ; in the common dense-shift case, where the mean contribution is spread across the 𝑝 coordinates at comparable scale, the Frobenius accumulation is written as ‖ ‖‖ Gram ‖‖Σ̂𝑡 − Σ𝑡 ‖‖‖ ≳ 𝑝𝑐2 𝜀mach . ‖𝐹 ‖ For Welford, the update is based on residuals ∆𝑡 = 𝐱𝑡 − µ𝑡−1 ,
𝐱𝑡 − µ 𝑡 .
The common shift µ cancels in these differences. Consequently, once the running mean has reached the scale of the sample mean, both residual factors are governed by the centred variables 𝐳𝑖 , not by the absolute location µ. Each outer-product correction therefore has entries of order 𝜎𝑘 𝜎𝑙 , and the floating-point error accumulates at covariance scale rather than raw-data scale: || Welf | ||Σ̂𝑘𝑙,𝑡 − Σ𝑘𝑙,𝑡 ||| = 𝑂(𝜎𝑘 𝜎𝑙 𝜀mach ). Taking the Frobenius norm over the 𝑝 × 𝑝 matrix yields ‖‖ Welf ‖ ‖‖Σ̂𝑡 − Σ𝑡 ‖‖‖ = 𝑂(𝑝𝜎2 𝜀mach ), ‖ ‖𝐹 which does not depend on 𝑐.
13.6 Proof of Theorem 6.1: Conformal Coverage (𝑗)
Proof. Let 𝑍1 , … , 𝑍𝑚 , 𝑍𝑚+1 be i.i.d. random variables, where 𝑍𝑗 = |Σ̂𝑡 𝑘𝑙 − Σ𝑘𝑙 | for 𝑗 ≤ 𝑚 (calibration) and 𝑍𝑚+1 is the test score. By exchangeability, for any fixed threshold 𝑞, Pr(𝑍𝑚+1 > 𝑞) = Pr(𝑍1 > 𝑞).
20
2𝐵 or Not 2𝐵: A Tale of Three Algorithms for Streaming
The conformal quantile 𝑞̂𝛼+ = Quantile({𝑍𝑗 }𝑚 ; ⌈(𝑚 + 1)(1 − 𝛼)⌉∕𝑚) is defined so that 𝑞̂𝛼+ ≥ 𝑍(𝑘∗ ) 𝑗=1 where 𝑘 ∗ = ⌈(𝑚+1)(1−𝛼)⌉ is the 𝑘 ∗ -th order statistic. By the standard conformal coverage argument [9, Theorem 1], Pr(𝑍𝑚+1 ≤ 𝑞̂𝛼+ ) ≥ 1 − 𝛼. (𝑚+1)
Since 𝑍𝑚+1 ≤ 𝑞̂𝛼+ is equivalent to |Σ̂𝑡 coverage guarantee follows. □
+ ̂ 𝑘𝑙 − Σ𝑘𝑙 | ≤ 𝑞̂𝛼 , which is the event that Σ𝑘𝑙 ∈ 𝐶𝛼 (𝑡), the
13.7 Bariance Scalar Identity Definition 13.1 (Bariance [8]). For 𝑛 ≥ 2 and 𝑥1 , … , 𝑥𝑛 ∈ ℝ, Bar(𝑥1 , … , 𝑥𝑛 ) =
∑ 1 (𝑥𝑖 − 𝑥𝑗 )2 . 2𝑛(𝑛 − 1) 𝑖≠𝑗
Proposition 13.2. Bar(𝑥) = (𝑛𝑆𝑥𝑥 − 𝑆𝑥2 )∕[𝑛(𝑛 − 1)], where 𝑆𝑥 = 1 ∑ ̄ 2. Bar(𝑥) = 𝜎̂ 2 = (𝑥 − 𝑥) 𝑖 𝑖
∑ 𝑖
𝑥𝑖 and 𝑆𝑥𝑥 =
∑ 𝑖
𝑥𝑖2 . Moreover,
𝑛−1
Proof. Expand (𝑥𝑖 − 𝑥𝑗 )2 = 𝑥𝑖2 − 2𝑥𝑖 𝑥𝑗 + 𝑥𝑗2 and sum over 𝑖 ≠ 𝑗. Using ∑ 𝑥 𝑥 = 𝑆𝑥2 − 𝑆𝑥𝑥 : 𝑖≠𝑗 𝑖 𝑗 ∑
∑ 𝑖≠𝑗
𝑥𝑖2 = (𝑛 − 1)𝑆𝑥𝑥 and
(𝑥𝑖 − 𝑥𝑗 )2 = 2(𝑛 − 1)𝑆𝑥𝑥 − 2(𝑆𝑥2 − 𝑆𝑥𝑥 ) = 2𝑛𝑆𝑥𝑥 − 2𝑆𝑥2 .
𝑖≠𝑗
Divide by 2𝑛(𝑛 − 1). The equality with 𝜎̂ 2 follows from algebra.
21
∑ 𝑖
̄ 2 = 𝑆𝑥𝑥 − 𝑆𝑥2 ∕𝑛 and standard (𝑥𝑖 − 𝑥)