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_

COMMENT:
@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
```

Labels: New Feature