ENH: Parallelization support for pairwise correlation

meson.build

@@ -15,6 +15,18 @@ py = import('python').find_installation(pure: false)
 tempita = files('generate_pxi.py')
 versioneer = files('generate_version.py')
 
+# Find OpenMP dependency
+openmp_dep = dependency('openmp', required: false)
+if not openmp_dep.found()
+    # Try to find OpenMP using compiler flags
+    cc = meson.get_compiler('c')
+    if cc.has_argument('-fopenmp')
+        openmp_dep = declare_dependency(
+            compile_args: ['-fopenmp'],
+            link_args: ['-fopenmp'],
+        )
+    endif
+endif
 
 add_project_arguments('-DNPY_NO_DEPRECATED_API=0', language: 'c')
 add_project_arguments('-DNPY_NO_DEPRECATED_API=0', language: 'cpp')
@@ -28,6 +40,14 @@ add_project_arguments(
     language: 'cpp',
 )
 
+# Add OpenMP compile args if available
+if openmp_dep.found()
+    add_project_arguments('-DHAVE_OPENMP', language: 'c')
+    add_project_arguments('-DHAVE_OPENMP', language: 'cpp')
+    message('OpenMP support enabled')
+else
+    message('OpenMP not found - parallel code will run sequentially')
+endif
 
 if fs.exists('_version_meson.py')
     py.install_sources('_version_meson.py', subdir: 'pandas')

pandas/_libs/algos.pyx

@@ -1,5 +1,8 @@
 cimport cython
 from cython cimport Py_ssize_t
+from cython.parallel cimport (
+    prange,
+)
 from libc.math cimport (
     fabs,
     sqrt,
@@ -345,60 +348,125 @@ def kth_smallest(numeric_t[::1] arr, Py_ssize_t k) -> numeric_t:
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-@cython.cdivision(True)
-def nancorr(const float64_t[:, :] mat, bint cov=False, minp=None):
+def nancorr(
+    const float64_t[:, :] mat,
+    bint cov=False,
+    minp=None,
+    bint use_parallel=False
+):
     cdef:
-        Py_ssize_t i, xi, yi, N, K
-        int64_t minpv
+        Py_ssize_t i, j, xi, yi, N, K
+        bint minpv
         float64_t[:, ::1] result
-        uint8_t[:, :] mask
+        # Initialize to None since we only use in the no missing value case
+        float64_t[::1] means = None
+        float64_t[::1] ssqds = None
+        ndarray[uint8_t, ndim=2] mask
+        bint no_nans
         int64_t nobs = 0
-        float64_t vx, vy, dx, dy, meanx, meany, divisor, ssqdmx, ssqdmy, covxy, val
+        float64_t mean, ssqd, val
+        float64_t vx, vy, dx, dy, meanx, meany, divisor, ssqdmx, ssqdmy, covxy
 
     N, K = (<object>mat).shape
     if minp is None:
         minpv = 1
     else:
-        minpv = <int64_t>minp
-
+        minpv = <int>minp
     result = np.empty((K, K), dtype=np.float64)
     mask = np.isfinite(mat).view(np.uint8)
+    no_nans = mask.all()
 
-    with nogil:
-        for xi in range(K):
-            for yi in range(xi + 1):
-                # Welford's method for the variance-calculation
-                # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
-                nobs = ssqdmx = ssqdmy = covxy = meanx = meany = 0
+    # Computing the online means and variances is expensive - so if possible we can
+    # precompute these and avoid repeating the computations each time we handle
+    # an (xi, yi) pair
+    if no_nans:
+        means = np.empty(K, dtype=np.float64)
+        ssqds = np.empty(K, dtype=np.float64)
+        with nogil:
+            for j in range(K):
+                ssqd = mean = 0
                 for i in range(N):
-                    if mask[i, xi] and mask[i, yi]:
+                    val = mat[i, j]
+                    dx = val - mean
+                    mean += 1 / (i + 1) * dx
+                    ssqd += (val - mean) * dx
+                means[j] = mean
+                ssqds[j] = ssqd
+
+    # ONLY CHANGE: Add parallel option to the main correlation loop
+    if use_parallel:
+        for xi in prange(K, schedule="dynamic", nogil=True):
+            for yi in range(xi + 1):
+                covxy = 0
+                if no_nans:
+                    for i in range(N):
                         vx = mat[i, xi]
                         vy = mat[i, yi]
-                        nobs += 1
-                        dx = vx - meanx
-                        dy = vy - meany
-                        meanx += 1. / nobs * dx
-                        meany += 1. / nobs * dy
-                        ssqdmx += (vx - meanx) * dx
-                        ssqdmy += (vy - meany) * dy
-                        covxy += (vx - meanx) * dy
-
+                        covxy += (vx - means[xi]) * (vy - means[yi])
+                    ssqdmx = ssqds[xi]
+                    ssqdmy = ssqds[yi]
+                    nobs = N
+                else:
+                    nobs = ssqdmx = ssqdmy = covxy = meanx = meany = 0
+                    for i in range(N):
+                        # Welford's method for the variance-calculation
+                        # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
+                        if mask[i, xi] and mask[i, yi]:
+                            vx = mat[i, xi]
+                            vy = mat[i, yi]
+                            nobs = nobs + 1
+                            dx = vx - meanx
+                            dy = vy - meany
+                            meanx = meanx + (1 / nobs * dx)
+                            meany = meany + (1 / nobs * dy)
+                            ssqdmx = ssqdmx + ((vx - meanx) * dx)
+                            ssqdmy = ssqdmy + ((vy - meany) * dy)
+                            covxy = covxy + ((vx - meanx) * dy)
                 if nobs < minpv:
                     result[xi, yi] = result[yi, xi] = NaN
                 else:
                     divisor = (nobs - 1.0) if cov else sqrt(ssqdmx * ssqdmy)
-
-                    # clip `covxy / divisor` to ensure coeff is within bounds
                     if divisor != 0:
-                        val = covxy / divisor
-                        if not cov:
-                            if val > 1.0:
-                                val = 1.0
-                            elif val < -1.0:
-                                val = -1.0
-                        result[xi, yi] = result[yi, xi] = val
+                        result[xi, yi] = result[yi, xi] = covxy / divisor
+                    else:
+                        result[xi, yi] = result[yi, xi] = NaN
+    else:
+        with nogil:
+            for xi in range(K):
+                for yi in range(xi + 1):
+                    covxy = 0
+                    if no_nans:
+                        for i in range(N):
+                            vx = mat[i, xi]
+                            vy = mat[i, yi]
+                            covxy += (vx - means[xi]) * (vy - means[yi])
+                        ssqdmx = ssqds[xi]
+                        ssqdmy = ssqds[yi]
+                        nobs = N
                     else:
+                        nobs = ssqdmx = ssqdmy = covxy = meanx = meany = 0
+                        for i in range(N):
+                            # Welford's method for the variance-calculation
+                            # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
+                            if mask[i, xi] and mask[i, yi]:
+                                vx = mat[i, xi]
+                                vy = mat[i, yi]
+                                nobs += 1
+                                dx = vx - meanx
+                                dy = vy - meany
+                                meanx += 1 / nobs * dx
+                                meany += 1 / nobs * dy
+                                ssqdmx += (vx - meanx) * dx
+                                ssqdmy += (vy - meany) * dy
+                                covxy += (vx - meanx) * dy
+                    if nobs < minpv:
                         result[xi, yi] = result[yi, xi] = NaN
+                    else:
+                        divisor = (nobs - 1.0) if cov else sqrt(ssqdmx * ssqdmy)
+                        if divisor != 0:
+                            result[xi, yi] = result[yi, xi] = covxy / divisor
+                        else:
+                            result[xi, yi] = result[yi, xi] = NaN
 
     return result.base
 

pandas/_libs/meson.build

@@ -1,3 +1,15 @@
+# Add OpenMP dependency detection at the top of the file
+openmp_dep = dependency('openmp', required: false)
+if not openmp_dep.found()
+    cc = meson.get_compiler('c')
+    if cc.has_argument('-fopenmp')
+        openmp_dep = declare_dependency(
+            compile_args: ['-fopenmp'],
+            link_args: ['-fopenmp'],
+        )
+    endif
+endif
+
 _algos_take_helper = custom_target(
     'algos_take_helper_pxi',
     output: 'algos_take_helper.pxi',
@@ -70,7 +82,7 @@ libs_sources = {
     # numpy include dir is implicitly included
     'algos': {
         'sources': ['algos.pyx', _algos_common_helper, _algos_take_helper],
-        'deps': _khash_primitive_helper_dep,
+        'deps': [_khash_primitive_helper_dep, openmp_dep],
     },
     'arrays': {'sources': ['arrays.pyx']},
     'groupby': {'sources': ['groupby.pyx']},

pandas/core/frame.py

@@ -11264,6 +11264,7 @@ def corr(
         method: CorrelationMethod = "pearson",
         min_periods: int = 1,
         numeric_only: bool = False,
+        use_parallel: bool = False,
     ) -> DataFrame:
         """
         Compute pairwise correlation of columns, excluding NA/null values.
@@ -11292,6 +11293,11 @@ def corr(
             .. versionchanged:: 2.0.0
                 The default value of ``numeric_only`` is now ``False``.
 
+        use_parallel : bool, default False
+            Use parallel computation for Pearson correlation.
+            Only effective for large matrices where parallelization overhead
+            is justified by compute time savings.
+
         Returns
         -------
         DataFrame
@@ -11332,14 +11338,18 @@ def corr(
               dogs  cats
         dogs   1.0   NaN
         cats   NaN   1.0
+
+        >>> # Use parallel computation for large DataFrames
+        >>> large_df = pd.DataFrame(np.random.randn(10000, 100))
+        >>> corr_matrix = large_df.corr(use_parallel=True)
         """  # noqa: E501
         data = self._get_numeric_data() if numeric_only else self
         cols = data.columns
         idx = cols.copy()
         mat = data.to_numpy(dtype=float, na_value=np.nan, copy=False)
 
         if method == "pearson":
-            correl = libalgos.nancorr(mat, minp=min_periods)
+            correl = libalgos.nancorr(mat, minp=min_periods, use_parallel=use_parallel)
         elif method == "spearman":
             correl = libalgos.nancorr_spearman(mat, minp=min_periods)
         elif method == "kendall" or callable(method):

setup.py

@@ -438,6 +438,10 @@ def srcpath(name=None, suffix=".pyx", subdir="src"):
     "_libs.algos": {
         "pyxfile": "_libs/algos",
         "depends": _pxi_dep["algos"],
+        "extra_compile_args": extra_compile_args
+        + (["-fopenmp"] if not is_platform_windows() else ["/openmp"]),
+        "extra_link_args": extra_link_args
+        + (["-fopenmp"] if not is_platform_windows() else []),
     },
     "_libs.arrays": {"pyxfile": "_libs/arrays"},
     "_libs.groupby": {"pyxfile": "_libs/groupby"},