Use welford to calculate corr

Author: phofl

pandas/_libs/algos.pyx

@@ -268,7 +268,7 @@ def nancorr(const float64_t[:, :] mat, bint cov=False, minp=None):
         ndarray[float64_t, ndim=2] result
         ndarray[uint8_t, ndim=2] mask
         int64_t nobs = 0
-        float64_t vx, vy, sumx, sumy, sumxx, sumyy, meanx, meany, divisor
+        float64_t vx, vy, meanx, meany, divisor, prev_meany, prev_meanx, ssqdmx, ssqdmy, covxy
 
     N, K = (<object>mat).shape
 
@@ -283,37 +283,28 @@ def nancorr(const float64_t[:, :] mat, bint cov=False, minp=None):
     with nogil:
         for xi in range(K):
             for yi in range(xi + 1):
-                nobs = sumxx = sumyy = sumx = sumy = 0
+                nobs = ssqdmx = ssqdmy = covxy = meanx = meany = 0
                 for i in range(N):
                     if mask[i, xi] and mask[i, yi]:
                         vx = mat[i, xi]
                         vy = mat[i, yi]
                         nobs += 1
-                        sumx += vx
-                        sumy += vy
+                        prev_meanx = meanx
+                        prev_meany = meany
+                        meanx = meanx + 1 / nobs * (vx - meanx)
+                        meany = meany + 1 / nobs * (vy - meany)
+                        ssqdmx = ssqdmx + (vx - meanx) * (vx - prev_meanx)
+                        ssqdmy = ssqdmy + (vy - meany) * (vy - prev_meany)
+                        covxy = covxy + (vx - meanx) * (vy - prev_meany)
 
                 if nobs < minpv:
                     result[xi, yi] = result[yi, xi] = NaN
                 else:
-                    meanx = sumx / nobs
-                    meany = sumy / nobs
-
-                    # now the cov numerator
-                    sumx = 0
-                    for i in range(N):
-                        if mask[i, xi] and mask[i, yi]:
-                            vx = mat[i, xi] - meanx
-                            vy = mat[i, yi] - meany
-
-                            sumx += vx * vy
-                            sumxx += vx * vx
-                            sumyy += vy * vy
-
-                    divisor = (nobs - 1.0) if cov else sqrt(sumxx * sumyy)
+                    divisor = (nobs - 1.0) if cov else sqrt(ssqdmx * ssqdmy)
 
                     # numerical issues for constant columns
-                    if divisor > 1e-15:
-                        result[xi, yi] = result[yi, xi] = sumx / divisor
+                    if divisor != 0:
+                        result[xi, yi] = result[yi, xi] = covxy / divisor
                     else:
                         result[xi, yi] = result[yi, xi] = NaN
 
@@ -323,6 +314,7 @@ def nancorr(const float64_t[:, :] mat, bint cov=False, minp=None):
 # Pairwise Spearman correlation
 
 
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def nancorr_spearman(const float64_t[:, :] mat, Py_ssize_t minp=1) -> ndarray:

pandas/tests/frame/methods/test_cov_corr.py

@@ -208,7 +208,7 @@ def test_corr_item_cache(self):
         assert df["A"] is ser
         assert df.values[0, 0] == 99
 
-    @pytest.mark.parametrize("length", [2, 20, 200, 2000, 20000])
+    @pytest.mark.parametrize("length", [2, 20, 200, 2000])
     def test_corr_for_constant_columns(self, length):
         # GH: 37448
         df = DataFrame(length * [[0.4, 0.1]], columns=["A", "B"])
@@ -218,6 +218,15 @@ def test_corr_for_constant_columns(self, length):
         )
         tm.assert_frame_equal(result, expected)
 
+    def test_calc_corr_small_numbers(self):
+        # GH: 37452
+        df = DataFrame(
+            {"A": [1.0e-20, 2.0e-20, 3.0e-20], "B": [1.0e-20, 2.0e-20, 3.0e-20]}
+        )
+        result = df.corr()
+        expected = DataFrame({"A": [1.0, 1.0], "B": [1.0, 1.0]}, index=["A", "B"])
+        tm.assert_frame_equal(result, expected)
+
 
 class TestDataFrameCorrWith:
     def test_corrwith(self, datetime_frame):