Preserve streaming covariance precision near large offsets

algorithm_lab/online_covariance.py

@@ -4,7 +4,8 @@
 sample is a bool selecting n-1 rather than n as divisor. Require one pair
 for population or two for sample covariance. Return a float; unrepresentable
 intermediate state raises ValueError. O(n) time and O(1) state using a
-paired Welford recurrence; observations need not fit in memory.
+paired Welford recurrence after shifting by the first pair to preserve precision
+near large common offsets; observations need not fit in memory.
 
 >>> online_covariance([(1, 2), (2, 4), (3, 6)], sample=True)
 2.0
@@ -19,8 +20,12 @@ def online_covariance(observations, *, sample=False):
     if not isinstance(sample, bool):
         raise ValueError("sample must be a boolean")
     count, mean_x, mean_y, moment = 0, 0.0, 0.0, 0.0
+    origin_x = origin_y = 0.0
     for x, y in observations:
         x, y = finite(x), finite(y)
+        if count == 0:
+            origin_x, origin_y = x, y
+        x, y = x - origin_x, y - origin_y
         count += 1
         delta_x, delta_y = x - mean_x, y - mean_y
         mean_x += delta_x / count

tests/test_online_covariance.py

@@ -7,6 +7,16 @@
 
 
 class Tests(unittest.TestCase):
+    def test_exact_covariance_at_large_offsets(self):
+        # The centered observations are (-1, -1), (1, 1): sum products = 2.
+        for offset in [1e16, -1e16]:
+            values = [offset, offset + 2]
+            pairs = list(zip(values, values, strict=True))
+            self.assertEqual(online_covariance(iter(pairs), sample=True), 2)
+            self.assertEqual(online_covariance(iter(pairs)), 1)
+            reverse_pairs = zip(values, reversed(values), strict=True)
+            self.assertEqual(online_covariance(reverse_pairs, sample=True), -2)
+
     def test_standard_library_oracle(self):
         rng = random.Random(256)
         for n in range(2, 60):