add: parallelization support for pearson correlation

Author: gangula-karthik

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,
@@ -346,61 +349,101 @@ 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
         float64_t[:, ::1] result
         uint8_t[:, :] mask
-        int64_t nobs = 0
-        float64_t vx, vy, dx, dy, meanx, meany, divisor, ssqdmx, ssqdmy, covxy, val
+        int64_t nobs
+        float64_t vx, vy, dx, dy, meanx, meany
+        float64_t prev_meanx, prev_meany
+        float64_t ssqdmx, ssqdmy, covxy, divisor, val
 
     N, K = (<object>mat).shape
-    if minp is None:
-        minpv = 1
-    else:
-        minpv = <int64_t>minp
+    minpv = 1 if minp is None else <int64_t>minp
 
     result = np.empty((K, K), dtype=np.float64)
     mask = np.isfinite(mat).view(np.uint8)
 
-    with nogil:
-        for xi in range(K):
+    # Welford's method for the variance-calculation
+    # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
+
+    if use_parallel:
+        for xi in prange(K, schedule="dynamic", nogil=True):
             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
                 for i in range(N):
                     if mask[i, xi] and mask[i, yi]:
                         vx = mat[i, xi]
                         vy = mat[i, yi]
-                        nobs += 1
+                        nobs = 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
+                        meanx = meanx + dx / nobs
+                        meany = meany + dy / nobs
+                        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
+                    val = NaN
                 else:
-                    divisor = (nobs - 1.0) if cov else sqrt(ssqdmx * ssqdmy)
-
-                    # clip `covxy / divisor` to ensure coeff is within bounds
-                    if divisor != 0:
+                    if cov:
+                        divisor = nobs - 1.0
+                    else:
+                        divisor = sqrt(ssqdmx * ssqdmy)
+                    if divisor == 0.0:
+                        val = NaN
+                    else:
                         val = covxy / divisor
                         if not cov:
                             if val > 1.0:
                                 val = 1.0
-                            elif val < -1.0:
+                            if val < -1.0:
                                 val = -1.0
-                        result[xi, yi] = result[yi, xi] = val
-                    else:
+
+                result[xi, yi] = result[yi, xi] = val
+    else:
+        # Sequential version
+        with nogil:
+            for xi in range(K):
+                for yi in range(xi + 1):
+                    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
+                            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:
+                        divisor = (nobs - 1.0) if cov else sqrt(ssqdmx * ssqdmy)
+                        if divisor != 0:
+                            val = covxy / divisor
+                            if not cov:
+                                if val > 1.0:
+                                    val = 1.0
+                                if val < -1.0:
+                                    val = -1.0
+                            result[xi, yi] = result[yi, xi] = val
+                        else:
+                            result[xi, yi] = result[yi, xi] = NaN
 
-    return result.base
+    return result
 
 # ----------------------------------------------------------------------
 # Pairwise Spearman correlation

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,
+        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``.
 
+        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(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, parallel=parallel)
         elif method == "spearman":
             correl = libalgos.nancorr_spearman(mat, minp=min_periods)
         elif method == "kendall" or callable(method):