Repository: scikit-learn/scikit-learn
Pull request: DOC clarify that n_jobs argument uses threading in pairwise_distances and pairwise_kernels
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/main/CONTRIBUTING.md
-->

#### Reference Issues/PRs
<!--
Example: Fixes #1234. See also #3456.
Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->
Fixes #29587 

#### What does this implement/fix? Explain your changes.

Improves the documentation of `pairwise_distances` and `pairwise_kernels` so that it is clear that they use threading when the `n_jobs` argument is used. One could expect the functions to be truly parallel using multiprocessing otherwise, which is not the case and that has a serial performance impact for user-defined metric functions (see #29587).

Problem statement:
Truly parallel execution of pairwise_kernels and pairwise_distances
### Describe the workflow you want to enable

Both `pairwise_kernels` and `pairwise_distances` functions call `_parallel_pairwise` function, which is (contrary to its name) not parallel as it enforces the threading backend. Therefore, these functions are terribly slow, especially for computationally expensive user-defined metrics. I understand that the reasons for the threading backend are possibly large memory demands and data communication overhead but I suggest a different approach. Also, the documentation for these functions talks about parallel execution and processes which is currently simply not true.

### Describe your proposed solution

The memory and data communication issues can be reduced by a smarter distribution of the input data to individual processes. Right now, only `Y` is sliced in the `_parallel_pairwise` function which is suboptimal for parallel processing. Both `X` and `Y` should be sliced to lower the demands for multiprocessing. For example for 100x100 `X` and `Y` distributed to 100 processes, we have to copy 100+1 inputs to every process when slicing only `Y` while only 10+10 when slicing both `X` and `Y`. As a result, multiprocessing can be allowed. Also, joblib does automatic memmapping in some cases.

Alternatively, at least the documentation for  `pairwise_kernels` and `pairwise_distances` should be corrected.

### Describe alternatives you've considered, if relevant

_No response_

### Additional context

_No response_

Hints from issue discussion:
> which is (contrary to its name) not parallel as it enforces the threading backend.
> ... Therefore, these functions are terribly slow, especially for computationally expensive user-defined metrics.

Some CPU intensive operations in numpy or pandas do release the GIL, so the threading backend might be a good choice.

Furthermore, CPython 3.13 will come with an optional build flag that should a allow to run without any GIL (also known as the "free threading" mode of Python). It's not yet operational (it can cause segfaults in some cases) but @ngoldbaum, @lesteve and others are making progress at identifying and fixing the remaining problems. See #28978 for a tracking issue on the scikit-learn side.

It could be interesting to see if that can help with your specific workload.

More generally, before deciding to embark on a specific refactoring strategy for `_parallel_pairwise`, it would be interesting to see working code for typical uses of `pairwise_kernels` and `pairwise_distances` where you identified performance problems, and maybe report some quick benchmark results on the impact of varying the `n_jobs` parameter.

> Also, joblib does automatic memmapping in some cases.

I am worried about relying on this too much. In retrospect I have the feeling that it's a bit of brittle black magic that is complex to debug and maintain and reason about its performance implications (both in terms of memory usage and computational overhead).

More generally, process based parallelism can lead to hard to debug situations (see the recently opened #29579 for instance), so I would rather rely more on thread-based parallelism than process-based parallelism for scikit-learn in the future.
Hi @ogrisel,
thanks for your reply.

> Some CPU intensive operations in numpy or pandas do release the GIL, so the threading backend might be a good choice.

True, that's why wrote that it is especially severe for user-defined CPU-bound metric functions. I will try to put together some benchmark for a common sklearn kernel to see the effect there.

> Furthermore, CPython 3.13 will come with an optional build flag that should a allow to run without any GIL (also known as the "free threading" mode of Python). It's not yet operational (it can cause segfaults in some cases) but @ngoldbaum, @lesteve and others are making progress at identifying and fixing the remaining problems. See #28978 for a tracking issue on the scikit-learn side.
>
> It could be interesting to see if that can help with your specific workload.

That would be super nice to have. The question is how long it might take to make it practically usable in sklearn.

> More generally, before deciding to embark on a specific refactoring strategy for `_parallel_pairwise`, it would be interesting to see working code for typical uses of `pairwise_kernels` and `pairwise_distances` where you identified performance problems, and maybe report some quick benchmark results on the impact of varying the `n_jobs` parameter.

I observed a huge impact for my user-defined CPU-intensive metric function but for most people, predefined sklearn metrics are more relevant so I will benchmark that first.

> > Also, joblib does automatic memmapping in some cases.
> 
> I am worried about relying on this too much. In retrospect I have the feeling that it's a bit of brittle black magic that is complex to debug and maintain and reason about its performance implications (both in terms of memory usage and computational overhead).

Totally agree, automatic memmapping was just a side note

> More generally, process based parallelism can lead to hard to debug situations (see the recently opened #29579 for instance), so I would rather rely more on thread-based parallelism than process-based parallelism for scikit-learn in the future.

Sure, threading is better if it provides the same speed-up, it is just very often not the case with python. That's why many sklearn functions use the default joblib backend based on multiprocessing. If the free threading proves to work well and will be available in sklearn, it might help a lot. But as of now, multiprocessing is the only way to get full parallelism.

I will soon provide some benchmark data.
>> Furthermore, CPython 3.13 will come with an optional build flag that should a allow to run without any GIL (also known as the "free threading" mode of Python)

> That would be super nice to have. The question is how long it might take to make it practically usable in sklearn.

This is already usable in the sense that we have a CI for CPython 3.13 free-threaded and all the scikit-learn tests pass. I you have time to try CPython 3.13 free-threaded on your particular use case, feed-back would be super useful and welcome :pray:.

As of today, you have a few different ways to install CPython 3.13 free-threaded locally (through conda, your Linux distribution package manager, python.org download, etc ...) see the [py-free-threading.github.io doc](https://py-free-threading.github.io/installing_cpython/). numpy, scipy, scikit-learn, and other projects have development wheels for CPython 3.13 free-threaded see [this doc](https://py-free-threading.github.io/tracking/).

With this kind of "ongoing work" things, there may be some caveats along the road, for example my understanding is that there is currently a hit on single-threaded performance in CPython 3.13 free-threaded, but this will be improved in the future. Not sure about the details, this may be quickly mentioned at one point in [Anthony Shaw's PyCon US 2024 talk video](https://www.youtube.com/watch?v=lCjsy7tb420&t=1825).

@ogrisel 
I benchmarked the `_parallel_pairwise` function using the vectorized sklearn `rbf_kernel` metric and also non-vectorized user-defined rbf metric (as a simple example). Therefore, the 1st case is IO-bound and the 2nd case is CPU-bound. I tested both using the original  `_parallel_pairwise` function and also modified versions allowing for multiprocessing for both Y slicing and X+Y slicing (not sure if my slicing is optimal though). Here are the results for the IO-bound case (linear and log scale):
<p float="left">
  <img src="https://github.com/user-attachments/assets/2f4bea78-e561-4dec-bb8f-a84d3deb5182" width="400" />
  <img src="https://github.com/user-attachments/assets/98b17986-8842-4cc8-a533-4f0f8e7e63f0" width="400" /> 
</p>

Multithreading works in general better there. However, for the CPU-bound case, multiprocessing is way better:
<p float="left">
  <img src="https://github.com/user-attachments/assets/b3c0f04d-5e13-4585-8d22-cb71f7532d7d" width="400" />
  <img src="https://github.com/user-attachments/assets/450288e6-2629-49ed-936a-04f84b8e4705" width="400" /> 
</p>

In fact, threading has no effect in the case of CPU-bound task. Also, slicing of both X and Y decreases the communication and memory overhead and maybe there is even better way of slicing than I used. You can find the code below. The problem I see is that the user does not have the option to choose between threading and multiprocessing. As the joblib manual states, one should not hardcode the backend and use `prefer` option instead. In such case, the user can change the backend from outside the function. The problem is that the current version of  `_parallel_pairwise` writes to the output array in place, which is the fastest way for threading, but also prevents using multiprocessing.

@lesteve 
Good to know. I hope I will find some time to test it. I am really thrilled about it.

```py
import numpy as np
from sklearn.metrics.pairwise import rbf_kernel, _parallel_pairwise, _return_float_dtype, euclidean_distances, _pairwise_callable
from sklearn.utils import gen_even_slices
from sklearn.utils.validation import _num_samples
import matplotlib.pyplot as plt
from joblib import parallel_config, effective_n_jobs, delayed, Parallel
from functools import partial
import time
import math

def rbf_kernel2(x, y, gamma=None):
    if gamma is None:
        gamma = 1.0 / len(x)
    k = 0.0
    for i in range(len(x)):
        k += (x[i]-y[i])**2
    k *= -gamma
    k = math.exp(k)
    return k

def _parallel_pairwise2(X, Y, func, n_jobs, **kwds):
    """Break the pairwise matrix in n_jobs even slices
    and compute them in parallel."""

    if Y is None:
        Y = X
    X, Y, dtype = _return_float_dtype(X, Y)

    if effective_n_jobs(n_jobs) == 1:
        return func(X, Y, **kwds)
    
    ret = np.empty((X.shape[0], Y.shape[0]), dtype=dtype, order="F")
    slices = list(gen_even_slices(_num_samples(Y), effective_n_jobs(n_jobs)))
    out = Parallel(n_jobs=n_jobs)(
        delayed(func)(X, Y[s], **kwds)
        for s in slices
    )
    i = 0
    for i, s in enumerate(slices):
        ret[:, s] = out[i]

    return ret

def _parallel_pairwise3(X, Y, func, n_jobs, **kwds):
    """Break the pairwise matrix in n_jobs even slices
    and compute them in parallel."""

    if Y is None:
        Y = X
    X, Y, dtype = _return_float_dtype(X, Y)

    if effective_n_jobs(n_jobs) == 1:
        return func(X, Y, **kwds)
    
    ret = np.empty((X.shape[0], Y.shape[0]), dtype=dtype, order="F")
    eff_jobs = effective_n_jobs(n_jobs)
    n_jobs1 = int(np.sqrt(eff_jobs))
    n_jobs2 = int(eff_jobs//n_jobs1)
    # print(n_jobs1, n_jobs2)
    slices1 = list(gen_even_slices(_num_samples(X), n_jobs1))
    slices2 = list(gen_even_slices(_num_samples(Y), n_jobs2))
    out = Parallel(n_jobs=n_jobs)(
        delayed(func)(X[s1], Y[s2], **kwds)
        for s1 in slices1 for s2 in slices2
    )
    i = 0
    for s1 in slices1:
        for s2 in slices2:
            ret[s1, s2] = out[i]
            i += 1

    return ret

for n_samples, metric, title in [(5000, rbf_kernel, 'IO-bound'), (100, partial(_pairwise_callable, metric=rbf_kernel2), 'CPU-bound')]:
    # n_samples = 1000
    n_jobs = np.arange(1,11)
    # metric = rbf_kernel
    # metric = partial(_pairwise_callable, metric=rbf_kernel2)
    X, Y = np.random.rand(n_samples, 500), np.random.rand(n_samples, 500)
    results = [[] for i in range(5)]
    for n in n_jobs:
        with parallel_config(backend='threading'):
            result = %timeit -o -n 1 -r 10 _parallel_pairwise(X, Y, metric, n_jobs=n)
            results[0].append(result)
            result = %timeit -o -n 1 -r 10 _parallel_pairwise2(X, Y, metric, n_jobs=n)
            results[1].append(result)
            result = %timeit -o -n 1 -r 10 _parallel_pairwise3(X, Y, metric, n_jobs=n)
            results[2].append(result)
        with parallel_config(backend='loky'):
            result = %timeit -o -n 1 -r 10 _parallel_pairwise2(X, Y, metric, n_jobs=n)
            results[3].append(result)
            result = %timeit -o -n 1 -r 10 _parallel_pairwise3(X, Y, metric, n_jobs=n)
            results[4].append(result)
        print()
    times = []
    for i in range(len(results)):
        times.append([np.mean(result.timings) for result in results[i]])
    for scale in ['linear', 'log']:
        plt.figure()
        plt.title('_parallel_pairwise, '+title+', '+str(scale)+' scale, '+str(n_samples)+' samples')
        plt.xlim(np.min(n_jobs),np.max(n_jobs))
        if scale=='linear':
            plt.ylim(0,np.max(times))
        plt.xlabel('# of CPUs')
        plt.ylabel('times [s]')
        plt.xscale(scale)
        plt.yscale(scale)
        plt.plot(n_jobs, times[0], label='threading: original')
        plt.plot(n_jobs, times[1], label='threading: mod. Y slicing')
        plt.plot(n_jobs, times[2], label='threading: mod. X+Y slicing')
        plt.plot(n_jobs, times[3], label='multiprocessing: mod. Y slicing')
        plt.plot(n_jobs, times[4], label='multiprocessing: mod. X+Y slicing')
        plt.legend()
        plt.show()
```

EDIT:
output of `$ python -c 'import sklearn; sklearn.show_versions()'`

```
System:
    python: 3.12.4 | packaged by Anaconda, Inc. | (main, Jun 18 2024, 15:12:24) [GCC 11.2.0]
executable: /user/srsen/bin/anaconda3/envs/cheml/bin/python
   machine: Linux-4.18.0-348.23.1.el8_5.x86_64-x86_64-with-glibc2.28

Python dependencies:
      sklearn: 1.4.2
          pip: 24.0
   setuptools: 69.5.1
        numpy: 1.26.4
        scipy: 1.13.1
       Cython: 3.0.10
       pandas: None
   matplotlib: 3.8.4
       joblib: 1.4.2
threadpoolctl: 3.5.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: mkl
    num_threads: 20
         prefix: libmkl_rt
       filepath: /user/srsen/bin/anaconda3/envs/cheml/lib/libmkl_rt.so.2
        version: 2023.1-Product
threading_layer: intel

       user_api: openmp
   internal_api: openmp
    num_threads: 20
         prefix: libiomp
       filepath: /user/srsen/bin/anaconda3/envs/cheml/lib/libiomp5.so
        version: None

       user_api: openmp
   internal_api: openmp
    num_threads: 20
         prefix: libgomp
       filepath: /user/srsen/bin/anaconda3/envs/cheml/lib/libgomp.so.1.0.0
        version: None
```

Thank you very much for your analysis. Looking forward to see the results with free-threading Python if someone has time to set it up.
@stepan-srsen could you add the output of `python -c 'import sklearn; sklearn.show_versions()'` to your earlier comment i.e. https://github.com/scikit-learn/scikit-learn/issues/29587#issuecomment-2260837369, this could be useful for further reference :pray:.

I am trying to run your snippet and full-disclosure for now I am getting weird results on vanilla Python, so I must be doing something wrong I need to have a a closer look ...

@lesteve @ogrisel 
Hi,
I found some time to test the free-threaded version and here are the results. For the IO-bound case without setting the `OPENBLAS_NUM_THREADS` environment variable (and therefore using also the internal OPENBLAS threading):

<p float="left">
  <img src="https://github.com/user-attachments/assets/62e6c94d-d64a-4ed4-9080-8e56f41b0a68" width="400" />
  <img src="https://github.com/user-attachments/assets/6ff39aed-ce51-4c88-9276-a97fcdc63fa9" width="400" /> 
</p>

Everything, even multiprocessing was tested with the GIL turned off. None of the approaches actually accelerates the calculation in this case/setup when allowing parallelism inside the `_parallel_pairwise` function. However, if I set `OPENBLAS_NUM_THREADS=1`, I get the following results:

<p float="left">
  <img src="https://github.com/user-attachments/assets/01ed338a-bd29-4c9f-810b-93bbee90ab3e" width="400" />
  <img src="https://github.com/user-attachments/assets/6852b534-2424-422b-84cf-fa74d8ad98fb" width="400" /> 
</p>

These results correspond to the results from the regular python setup above. Interestingly, the BLAS threading can speed up single-core execution but its interplay with the parallelism used within the `_parallel_pairwise` function can actually lead to worse results.

For the CPU-bound tasks, the picture is clear, free-threading performs the same as multiprocessing or even slightly better (for numbers of CPUs for which the slicing is suboptimal) and the best times are basically the same as in my previous setup. I set `OPENBLAS_NUM_THREADS=1` but it does not play a role in this use case.

<p float="left">
  <img src="https://github.com/user-attachments/assets/03af8554-4591-4f0a-8b11-b8a7d89146da" width="400" /> 
  <img src="https://github.com/user-attachments/assets/5815af75-588c-4e25-a8b1-08575676b4a5" width="400" />
</p>

Dependencies versions in the new setup:

```
System:
    python: 3.13.0b4 | packaged by Anaconda, Inc. | experimental free-threading build (main, Jul 24 2024, 20:51:55) [GCC 11.2.0]
executable: /user/srsen/bin/anaconda3/envs/nogil/bin/python
   machine: Linux-4.18.0-348.23.1.el8_5.x86_64-x86_64-with-glibc2.28

Python dependencies:
      sklearn: 1.6.dev0
          pip: 24.1.2
   setuptools: None
        numpy: 2.1.0.dev0
        scipy: 1.15.0.dev0
       Cython: None
       pandas: None
   matplotlib: 3.10.0.dev503+gf4f40ba46f
       joblib: 1.4.2
threadpoolctl: 3.5.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 20
         prefix: libscipy_openblas
       filepath: /user/srsen/bin/anaconda3/envs/nogil/lib/python3.13t/site-packages/numpy.libs/libscipy_openblas64_-99b71e71.so
        version: 0.3.27
threading_layer: pthreads
   architecture: Haswell

       user_api: blas
   internal_api: openblas
    num_threads: 20
         prefix: libscipy_openblas
       filepath: /user/srsen/bin/anaconda3/envs/nogil/lib/python3.13t/site-packages/scipy.libs/libscipy_openblas-c128ec02.so
        version: 0.3.27.dev
threading_layer: pthreads
   architecture: Haswell

       user_api: openmp
   internal_api: openmp
    num_threads: 20
         prefix: libgomp
       filepath: /user/srsen/bin/anaconda3/envs/nogil/lib/python3.13t/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0
        version: None
```
@stepan-srsen just to make sure you do set `MKL_NUM_THREADS=1` (vanilla Python run since you are using MKL in your environment) or `OPENBLAS_NUM_THREADS=1` (free-threaded Python since you are using OpenBLAS here)?

I think not setting `MKL_NUM_THREADS=1` was one of the reason I was seeing surprising results originally.

I was also wondering how the outer-loop parallelism performance (i.e. what you are showing) with the BLAS library in single-threaded mode, compares to no parallellism in the outer loop and let the BLAS library use all the cores, but I did not have time to look at this ...


> For the CPU-bound tasks, the picture is clearer, free-threading performs same as multiprocessing or even slightly better (for numbers of CPUs for which the slicing si suboptimal) and the best times are basically the same as in my previous setup.

That's already great news! Thanks very much for taking the time to run this.

For the "IO-bound" case it would be worth investigating more about the interaction with different implementations of multithreaded BLAS used internally by numpy and scikit-learn (via scipy).
Thanks for updating the benchmark results with `OPENBLAS_NUM_THREADS=1`. So all in all, free-threading seems to work as intended but might require some manual protection against oversubscription with BLAS (for instance using `threadpoolctl.threadpool_limit`) as we typically do for BLAS calls nested under Cython/prange/OpenMP in other parts of the scikit-learn code base.

Given those results I think we can keep scikit-learn `pairwise_kernels` unchanged and instead allocate more efforts in making scikit-learn and joblib work well with free-threading CPython and the threading backend by protecting against oversubscription.

Commit message: clarification of the documentation for n_jobs argument

Patch to review:
diff --git a/sklearn/metrics/pairwise.py b/sklearn/metrics/pairwise.py
index fd08850ce47b3..771da7e1a1d9a 100644
--- a/sklearn/metrics/pairwise.py
+++ b/sklearn/metrics/pairwise.py
@@ -1938,7 +1938,7 @@ def _dist_wrapper(dist_func, dist_matrix, slice_, *args, **kwargs):
 
 def _parallel_pairwise(X, Y, func, n_jobs, **kwds):
     """Break the pairwise matrix in n_jobs even slices
-    and compute them in parallel."""
+    and compute them using threading."""
 
     if Y is None:
         Y = X
@@ -2333,8 +2333,8 @@ def pairwise_distances(
 
     n_jobs : int, default=None
         The number of jobs to use for the computation. This works by breaking
-        down the pairwise matrix into n_jobs even slices and computing them in
-        parallel.
+        down the pairwise matrix into n_jobs even slices and computing them
+        using threading.
 
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
@@ -2597,8 +2597,8 @@ def pairwise_kernels(
 
     n_jobs : int, default=None
         The number of jobs to use for the computation. This works by breaking
-        down the pairwise matrix into n_jobs even slices and computing them in
-        parallel.
+        down the pairwise matrix into n_jobs even slices and computing them
+        using threading.
 
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors. See :term:`Glossary <n_jobs>`


Code review comments:
diff_hunk: @@ -2597,8 +2597,8 @@ def pairwise_kernels(
 
     n_jobs : int, default=None
         The number of jobs to use for the computation. This works by breaking
-        down the pairwise matrix into n_jobs even slices and computing them in
-        parallel.
+        down the pairwise matrix into n_jobs even slices and computing them
+        using threading.
path: sklearn/metrics/pairwise.py
text: @user1:
```suggestion
        using multithreading.
```

diff_hunk: @@ -2333,8 +2333,8 @@ def pairwise_distances(
 
     n_jobs : int, default=None
         The number of jobs to use for the computation. This works by breaking
-        down the pairwise matrix into n_jobs even slices and computing them in
-        parallel.
+        down the pairwise matrix into n_jobs even slices and computing them
+        using threading.
path: sklearn/metrics/pairwise.py
text: @user1:
```suggestion
        using multithreading.
```