Talip7/scikit-learn-issues-multilabel huggingface dataset:
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:
> 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).

Labels: New Feature
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:
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.

Labels: New Feature
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:
@lesteve
Good point, I didn't take care of that in the second setup. I updated the results now.

@ogrisel @lesteve 
The free threading works nicely but it is still in the development stage and I assume that it won't be compatible with plenty of libraries for quite some time. As @ogrisel pointed out, the oversubscription should be also taken care of. I still think that the threading should not be hardcoded as the joblib manual suggests but it should use the `prefer` keyword instead. Or at least, the documentation for `pairwise_kernels` and `pairwise_distances` should be updated and clearly state that it is just threading. I was surprised when I turned on what I thought was multiprocessing but it didn't help with my CPU-intensive calculation. The standard approach is that one can define the joblib backend from outside of the function using `with parallel_config(backend=...):` as I did in the examples above.

Labels: New Feature
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:
> I still think that the threading should not be hardcoded as the joblib manual suggests but it should use the prefer keyword instead.

+1 for using the `prefer` keyword. Feel free to open a PR.

Labels: New Feature
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:
Hi,

@ogrisel @lesteve 
the problem is that the function currently writes to the array in place which prohibits multiprocessing. We can use `require='sharedmem'` instead of `prefer='threads'` but that doesn't really help. So we would have to write to the array outside of the parallel environment as I do in the modified functions above to allow for multiprocessing. Btw, I found out that NOT enforcing the Fortran order (i.e. not using `order="F"`) makes the code slightly faster in this case so the difference between writing in place and not writing in place is not that big for threading as it looks from the figures above. Maybe the best way would be to recognize if multiprocessing has been requested from outside of the function and write in place only for the threading backend. What do you think? Also, what kind of slicing should I use for the PR? In my own code, I slice both X and Y as it is the best for multiprocessing and I use prefer='threads' for predefined metrics and prefer='processes' for user-defined metrics. 

Also, another thought: Will the python 3.13 support free-threading by default or will it need some special compilation? If it is the second case, then sklearn shouldn't just assume that the user compiles the python with some experimental feature.

BTW you can remove the Needs Benchmarks and Needs Reproducible Code tags.

Labels: New Feature
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:
> Also, another thought: Will the python 3.13 support free-threading by default or will it need some special compilation? If it is the second case, than sklearn shouldn't just assume that the user compiles the python with some experimental feature.

The official python 3.13 binaries provided by python.org will not be free-threaded by default but some conda channel (and maybe other distributions) will ship both free-threaded and gil-based python versions. scikit-learn will probably ship wheels for both ABIs and users will be free to install what they want. If conda-forge is ready in time we will also ship free-threaded scikit-learn packages on conda-forge.

I don't think that the full Python ecosystem will be 100% free-threading ready at the time of the CPython 3.13 release (in October) but I hope that most of the core scientific packages will be ready.

My point is that I would rather invest dev effort in making free-threading work as best as possible as nearly as possible (in 2024 or 2025) than invest efforts in working around GIL and process-based paralellism limitations in the medium term.

I am fine with short term pragmatic improvements for parallel based parallelism if they do not add to much complexity to our code bases (scikit-learn, joblib, loky). though.

Labels: New Feature
⚠️ CI failed on Wheel builder (last failure: Aug 01, 2024) ⚠️

**CI is still failing on [Wheel builder](https://github.com/scikit-learn/scikit-learn/actions/runs/10191626508)** (Aug 01, 2024)

COMMENT:
Issue when getting pandas development wheel for Python 3.13 free-threaded, which is not available so it is trying to build from source and fails. Note the `tar.gz` (rather than the `.whl` in the [build log](https://github.com/scikit-learn/scikit-learn/actions/runs/10155774547/job/28083048523):
```
  Collecting pandas
    Downloading https://pypi.anaconda.org/scientific-python-nightly-wheels/simple/pandas/3.0.0.dev0%2B1268.g9c8c685f48/pandas-3.0.0.dev0%2B1268.g9c8c685f48.tar.gz (4.3 MB)
```

Last time I looked there was not really a way to tell `pip` to only consider wheels not `.tar.gz` so that at least the error would be less cryptic but maybe I have missed it. `--no-binary :none:` is actually not enough.
 
Don't know why pandas cp313-313t is not available, looks like there was some issue with Python 3.13 (not only free-threaded but also vanilla 3.13) https://github.com/pandas-dev/pandas/actions/runs/10155668934/job/28082792259 ...

You could argue that we could have used the previous nightly wheel but somehow looks like this pandas is cleaning up his anaconda.org bucket and the previous wheel is not there anymore?

Labels: Build / CI
⚠️ CI failed on Wheel builder (last failure: Aug 01, 2024) ⚠️

**CI is still failing on [Wheel builder](https://github.com/scikit-learn/scikit-learn/actions/runs/10191626508)** (Aug 01, 2024)

COMMENT:
Not sure what the short-term option is, either:
- revert https://github.com/scikit-learn/scikit-learn/pull/29572
- try to install pandas and if that fails still run the tests without pandas installed, this is still better than nothing ... but this may be confusing that the test can run sometimes with pandas sometimes without (only for the Python 3.13 free-threaded build to be explicit)

Labels: Build / CI
⚠️ CI failed on Wheel builder (last failure: Aug 01, 2024) ⚠️

**CI is still failing on [Wheel builder](https://github.com/scikit-learn/scikit-learn/actions/runs/10191626508)** (Aug 01, 2024)

COMMENT:
Maybe we can wait 48h to see if pandas dev can fix their Python 3.13 nightly builds. If it keeps failing, then we might consider reverting #29572.

Labels: Build / CI
⚠️ CI failed on Wheel builder (last failure: Aug 01, 2024) ⚠️

**CI is still failing on [Wheel builder](https://github.com/scikit-learn/scikit-learn/actions/runs/10191626508)** (Aug 01, 2024)

COMMENT:
I opened an issue about this https://github.com/pandas-dev/pandas/issues/59372, let's see what they say!

Labels: Build / CI
⚠️ CI failed on Wheel builder (last failure: Aug 01, 2024) ⚠️

**CI is still failing on [Wheel builder](https://github.com/scikit-learn/scikit-learn/actions/runs/10191626508)** (Aug 01, 2024)

COMMENT:
Actually I opened https://github.com/scikit-learn/scikit-learn/pull/29598 that should fix it. It turns out that pandas has wheels available for cp313-313t but that they are older than the sdist so `pip` choses the sdist ... `--only-binary :all:` is the way to tell pip to only use wheels.

Labels: Build / CI
Try running examples in parallel during doc build

We are already using sphinx-gallery 0.17 which has added the feature to run examples in parallel see https://github.com/sphinx-gallery/sphinx-gallery/pull/877. See [sphinx-gallery doc](https://sphinx-gallery.github.io/stable/configuration.html#parallel) for how to configure it.

matplotlib is currently trying it and it seems to show interesting improvements in their CI see https://github.com/matplotlib/matplotlib/pull/28617#issuecomment-2252108112.

I expect that for scikit-learn the speed-up may be a little bit less than for matplotlib since some examples are already using multiple cores (e.g. with `n_jobs=2`). I had a quick look during the sphinx-gallery PR and it was making the doc a bit quicker locally: https://github.com/sphinx-gallery/sphinx-gallery/pull/877#issuecomment-2184534205.

General directions:
- configure sphinx-gallery to use 2 cores in `doc/conf.py`
```py
sphinx_gallery_conf = {
    ...
    'parallel': 2,
}
```
- open a PR with `[doc build]` commit to do a full build
- also generate the doc locally e.g. with `spin docs clean` + `spin docs html` and see how much sphinx-gallery parallel settings make a difference

COMMENT:
> With parallel=12, this seems to hang (probably need to setup OPENBLAS_NUM_THREADS and other environment variables).

From your comment, I suspect there is oversubscription with parallelizing sphinx gallery. With `n_jobs=2`, it'll set `OPENBLAS_NUM_THREADS` to `joblib.cpu_count() // 2`.

Since [CicleCI](https://circleci.com/docs/configuration-reference/#docker-execution-environment) has 2 cores, I suspect we could set `OPENBLAS_NUM_THREADS=1` and use 2 cores for sphinx-gallery. This should work for a majority of examples. The only issue would be examples that set `n_jobs=-1` or `n_jobs>=2`, which would oversubscribe.

Labels: Build / CI
Try running examples in parallel during doc build

We are already using sphinx-gallery 0.17 which has added the feature to run examples in parallel see https://github.com/sphinx-gallery/sphinx-gallery/pull/877. See [sphinx-gallery doc](https://sphinx-gallery.github.io/stable/configuration.html#parallel) for how to configure it.

matplotlib is currently trying it and it seems to show interesting improvements in their CI see https://github.com/matplotlib/matplotlib/pull/28617#issuecomment-2252108112.

I expect that for scikit-learn the speed-up may be a little bit less than for matplotlib since some examples are already using multiple cores (e.g. with `n_jobs=2`). I had a quick look during the sphinx-gallery PR and it was making the doc a bit quicker locally: https://github.com/sphinx-gallery/sphinx-gallery/pull/877#issuecomment-2184534205.

General directions:
- configure sphinx-gallery to use 2 cores in `doc/conf.py`
```py
sphinx_gallery_conf = {
    ...
    'parallel': 2,
}
```
- open a PR with `[doc build]` commit to do a full build
- also generate the doc locally e.g. with `spin docs clean` + `spin docs html` and see how much sphinx-gallery parallel settings make a difference

COMMENT:
Yep this is a good point to bear in mind, while we are at it setting `OMP_NUM_THREADS=1` is probably a good idea as well.

Side-comment: matplotlib use a CircleCI large runner for its doc build (4 cores 8GB RAM, rather than 2 cores 4GB RAM) , so I guess this may be something to look at for scikit-learn but there may be some caveats.

Labels: Build / CI
Small discrepancy with not implemented MPS backend function

The discrepancy can be reproduced with the following snippet:

```python
# https://docs.scipy.org/doc/scipy/dev/api-dev/array_api.html
import os
os.environ["SCIPY_ARRAY_API"] = "1"
# os.environ["PYTORCH_ENABLE_MPS_FALLBACK"] = "1"

import torch

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn import config_context
from sklearn.pipeline import make_pipeline
from sklearn.decomposition import PCA
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_validate

X, y = load_diabetes(return_X_y=True)
X, y = X.astype(np.float32), y.astype(np.float32)

assert torch.backends.mps.is_available()
mps_device = torch.device("mps")
X = torch.from_numpy(X).to(device=mps_device)
y = torch.from_numpy(y).to(device=mps_device)

model = make_pipeline(PCA(svd_solver="full"), Ridge(solver="svd"))
with config_context(array_api_dispatch=True):
    cv_results = cross_validate(model, X, y)
cv_results
```

When using the following model:

```
model = make_pipeline(PCA(svd_solver="full"), Ridge(solver="svd"))
```

PyTorch will call `svd` solver and this is not available for MPS. In this case, it will fallback **automatically** on CPU and raise a warning.

However, for the following model:

```
model = make_pipeline(PCA(svd_solver="covariance_eigh"), Ridge(solver="svd"))
```

PyTorch will call `eigh` but this is also not supported for the MPS backend. However, in this case, scikit-learn is raising an error and request to explicitly set the environment variable `PYTORCH_ENABLE_MPS_FALLBACK`.

This is slightly weird that we are to consistent on how to treat those similar cases. We should check what is the reason and decide if we prefer to raise a warning or an error.

COMMENT:
In the `eigh` case this is the end of the traceback I get
```
  File "/Users/xxx/git/scikit-learn/sklearn/decomposition/_pca.py", line 612, in _fit_full
    eigenvals, eigenvecs = xp.linalg.eigh(C)
                           ^^^^^^^^^^^^^^^^^
NotImplementedError: The operator 'aten::_linalg_eigh.eigenvalues' is not currently implemented for the MPS device. If you want this op to be added in priority during the prototype phase of this feature, please comment on https://github.com/pytorch/pytorch/issues/77764. As a temporary fix, you can set the environment variable `PYTORCH_ENABLE_MPS_FALLBACK=1` to use the CPU as a fallback for this op. WARNING: this will be slower than running natively on MPS.
```
I think this means the exception is raised in PyTorch.

But yes, it is weird that sometimes it automatically falls back to CPU and sometimes it doesn't. I added a vote for it to be implemented in the issue the error links https://github.com/pytorch/pytorch/issues/77764#issuecomment-2252455129

Labels: Array API