Talip7/scikit-learn-issues-multilabel huggingface dataset:
OPTICS.fit leaks memory when called under VS Code's built-in debugger

### Describe the bug

Running clustering algorithm with n_jobs parameter set to more than 1 thread causes memory leak each time algorithm is run.
This simple code causes additional memory leak at each loop cycle. The issue will not occur if i replace manifold reduction algorithm with precomputed features.

### Steps/Code to Reproduce

```python
import gc
import numpy as np
from sklearn.manifold import TSNE
from sklearn.cluster import OPTICS
import psutil
process = psutil.Process()


def main():
    data = np.random.random((100, 100))
    for _i in range(1, 50):
        points = TSNE().fit_transform(data)
        prediction = OPTICS(n_jobs=2).fit_predict(points)  # n_jobs!=1
        points = None
        prediction = None
        del prediction
        del points
        gc.collect()
        print(f"{process.memory_info().rss / 1e6:.1f} MB")


main()
```

### Expected Results

Program's memory usage nearly constant between loop cycles

### Actual Results

Program's memory usage increases infinitely

### Versions

```shell
System:
    python: 3.10.6 (tags/v3.10.6:9c7b4bd, Aug  1 2022, 21:53:49) [MSC v.1932 64 bit (AMD64)]
executable: .venv\Scripts\python.exe
   machine: Windows-10-10.0.26100-SP0

Python dependencies:
      sklearn: 1.6.0
          pip: 24.3.1
   setuptools: 63.2.0
        numpy: 1.25.2
        scipy: 1.14.1
       Cython: None
       pandas: 2.2.3
   matplotlib: 3.10.0
       joblib: 1.4.2
threadpoolctl: 3.5.0

Built with OpenMP: True

threadpoolctl info:
       user_api: openmp
   internal_api: openmp
    num_threads: 16
         prefix: vcomp
       filepath: .venv\Lib\site-packages\sklearn\.libs\vcomp140.dll
        version: None

       user_api: blas
   internal_api: openblas
    num_threads: 16
         prefix: libopenblas
       filepath: .venv\Lib\site-packages\numpy\.libs\libopenblas64__v0.3.23-246-g3d31191b-gcc_10_3_0.dll
        version: 0.3.23.dev
threading_layer: pthreads
   architecture: Cooperlake

       user_api: blas
   internal_api: openblas
    num_threads: 16
         prefix: libscipy_openblas
       filepath: .venv\Lib\site-packages\scipy.libs\libscipy_openblas-5b1ec8b915dfb81d11cebc0788069d2d.dll
        version: 0.3.27.dev
threading_layer: pthreads
   architecture: Cooperlake
```

COMMENT:
> There are many processes, which are not closing without pool resizing.

This is expected: the loky backend of joblib keeps the workers alive for a while to avoid the overhead of creating and terminating processes upon successive calls to `joblib.Parallel`. This would not be a problem is the memory usage per worker stays stable over time.

But note that the "leak" already seems to start before using the loky backend in the output. It's already there when using the threading backend. To make the problem easier to investigate, I would start with investigating the memory usage of a minimal reproducer that only uses thread-based parallelism, for instance using the scalene memory profiler.

Labels: Bug, Performance, Needs Investigation
A helpful warning when adding sparsity constraints to NMF

### Describe the issue linked to the documentation

Currently the documentation of [NMF](https://scikit-learn.org/dev/modules/generated/sklearn.decomposition.NMF.html), as well as extensions like the [MiniBatchNMF](https://scikit-learn.org/dev/modules/generated/sklearn.decomposition.MiniBatchNMF.html) provide useful comments and warnings for beginners. For example, what type of initialization is suited based on whether sparsity is desired etc.

One thing that is not however addressed is the scale ambiguity of solving NMF. Specifically, if one desires sparsity on one of the matrices, say W, one should make sure the norm of the other matrix, in this case H, is controlled. Otherwise, a trivial solution would be a rescaled version of W and H, where the norm of W is decreased and the norm of H is increased. This would give the same exact (dot-product) output, while reducing both the L1 and L2 norm of the W matrix. If the user doesn't manually inspect the norms of H later, they may be mislead on what is actually happening. They may think that they have a more sparse factorization, whereas for the most part, they have arrived at a similar solution, just that the matrices have been rescaled. This would really hinder the actual sparsity of the factorization. You can also find this issue discussed in the last paragraph on the first page of [this paper](https://arxiv.org/pdf/2207.06316).

### Suggest a potential alternative/fix

I checked the literature and people often choose between doing Projected Gradient Descent (i.e to project the other matrix to a specified norm so that the model doesn't cheat) or doing a norm regularization on the other matrix as well. Since adding PGD would be too much of a change, I think letting the user know and maybe encouraging them to also add a sparsity constraint on the other matrix is the way to go.

I think there should be a simple warning when either one of `alpha_w` or `alpha_h` is enabled while the other is zero. It would simply warn the user about the scale ambiguity of the NMF. Otherwise, I think having it as part of the documentation would also work. But as of now, it can be very misleading, especially for people new to the NMF.

If you agree, I would love to make a pull request on this, as it would also be my first PR to a large project :)

COMMENT:
Thanks for the report. What you propose makes sense to me, although I would appreciate it if you could illustrate the problem on a small synthetic case to make sure that would describe happens in practice. I tried quickly, but I got convergence warnings when setting `alpha_h` to 0. But maybe it depends on the particular data.

Labels: Documentation
A helpful warning when adding sparsity constraints to NMF

### Describe the issue linked to the documentation

Currently the documentation of [NMF](https://scikit-learn.org/dev/modules/generated/sklearn.decomposition.NMF.html), as well as extensions like the [MiniBatchNMF](https://scikit-learn.org/dev/modules/generated/sklearn.decomposition.MiniBatchNMF.html) provide useful comments and warnings for beginners. For example, what type of initialization is suited based on whether sparsity is desired etc.

One thing that is not however addressed is the scale ambiguity of solving NMF. Specifically, if one desires sparsity on one of the matrices, say W, one should make sure the norm of the other matrix, in this case H, is controlled. Otherwise, a trivial solution would be a rescaled version of W and H, where the norm of W is decreased and the norm of H is increased. This would give the same exact (dot-product) output, while reducing both the L1 and L2 norm of the W matrix. If the user doesn't manually inspect the norms of H later, they may be mislead on what is actually happening. They may think that they have a more sparse factorization, whereas for the most part, they have arrived at a similar solution, just that the matrices have been rescaled. This would really hinder the actual sparsity of the factorization. You can also find this issue discussed in the last paragraph on the first page of [this paper](https://arxiv.org/pdf/2207.06316).

### Suggest a potential alternative/fix

I checked the literature and people often choose between doing Projected Gradient Descent (i.e to project the other matrix to a specified norm so that the model doesn't cheat) or doing a norm regularization on the other matrix as well. Since adding PGD would be too much of a change, I think letting the user know and maybe encouraging them to also add a sparsity constraint on the other matrix is the way to go.

I think there should be a simple warning when either one of `alpha_w` or `alpha_h` is enabled while the other is zero. It would simply warn the user about the scale ambiguity of the NMF. Otherwise, I think having it as part of the documentation would also work. But as of now, it can be very misleading, especially for people new to the NMF.

If you agree, I would love to make a pull request on this, as it would also be my first PR to a large project :)

COMMENT:
Thanks for the reply! I will post a synthetic code over this weekend.

Labels: Documentation
A helpful warning when adding sparsity constraints to NMF

### Describe the issue linked to the documentation

Currently the documentation of [NMF](https://scikit-learn.org/dev/modules/generated/sklearn.decomposition.NMF.html), as well as extensions like the [MiniBatchNMF](https://scikit-learn.org/dev/modules/generated/sklearn.decomposition.MiniBatchNMF.html) provide useful comments and warnings for beginners. For example, what type of initialization is suited based on whether sparsity is desired etc.

One thing that is not however addressed is the scale ambiguity of solving NMF. Specifically, if one desires sparsity on one of the matrices, say W, one should make sure the norm of the other matrix, in this case H, is controlled. Otherwise, a trivial solution would be a rescaled version of W and H, where the norm of W is decreased and the norm of H is increased. This would give the same exact (dot-product) output, while reducing both the L1 and L2 norm of the W matrix. If the user doesn't manually inspect the norms of H later, they may be mislead on what is actually happening. They may think that they have a more sparse factorization, whereas for the most part, they have arrived at a similar solution, just that the matrices have been rescaled. This would really hinder the actual sparsity of the factorization. You can also find this issue discussed in the last paragraph on the first page of [this paper](https://arxiv.org/pdf/2207.06316).

### Suggest a potential alternative/fix

I checked the literature and people often choose between doing Projected Gradient Descent (i.e to project the other matrix to a specified norm so that the model doesn't cheat) or doing a norm regularization on the other matrix as well. Since adding PGD would be too much of a change, I think letting the user know and maybe encouraging them to also add a sparsity constraint on the other matrix is the way to go.

I think there should be a simple warning when either one of `alpha_w` or `alpha_h` is enabled while the other is zero. It would simply warn the user about the scale ambiguity of the NMF. Otherwise, I think having it as part of the documentation would also work. But as of now, it can be very misleading, especially for people new to the NMF.

If you agree, I would love to make a pull request on this, as it would also be my first PR to a large project :)

COMMENT:
@m-parchami Hi, Could you illustrate the problem with code?

Labels: Documentation
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
I can reproduce the issue with numpy-dev (no issue with numpy 2.2 latest release)

```
pytest sklearn/metrics/tests/test_pairwise.py -k extreme_values
```

```
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=0
E       
E       Mismatched elements: 1 / 1 (100%)
E       Max absolute difference among violations: 2.74553895e-05
E       Max relative difference among violations: 0.00027451
E        ACTUAL: array([[0.100044]], dtype=float32)
E        DESIRED: array([[0.100017]])
```

A git bisect seems to indicate that https://github.com/numpy/numpy/pull/27883 is the culprit ... I commented in https://github.com/numpy/numpy/pull/27883#issuecomment-2554267160.

We are comparing `scipy.cdist` against our implementation of pairwise_distances:
https://github.com/scikit-learn/scikit-learn/blob/4ad187a7401f939c4d9cd27090c4258b5d810650/sklearn/metrics/tests/test_pairwise.py#L967-L989

This only fails with `dtype=float32` and `dim=1000000`. The comparison used to pass with `1e-5` and now the relative diff is `2.7e-4` so something looks kind of fishy.

`scipy.cdist` does not change but our implemention does somehow with https://github.com/numpy/numpy/pull/27883.

- when failure:
```
distances=array([[0.10004405]], dtype=float32)
expected=array([[0.10001659]])
```

- when success
```
distances=array([[0.10001709]], dtype=float32)
expected=array([[0.10001659]])
```

Labels: Bug
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
Moving here...
> Actually, complete guess, but maybe in row_norms which does the equivalent of np.sqrt((X * X).sum(axis=1)) with a np.einsum?
> https://github.com/scikit-learn/scikit-learn/blob/4ad187a7401f939c4d9cd27090c4258b5d810650/sklearn/utils/extmath.py#L47-L82

So, I am a bit surprised by two things:
* As far as I can tell, at least after calculating the differences, we should be using float64, and I am not sure if float64 should lead to this big a difference.
* I expected `einsum` to always use naive summation and in that case there should be no difference between the branches.  So if this is einsum and the difference is coming from there that may really need a closer look.

Labels: Bug
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
Thanks for your feed-back! In case there is no further insights from @jjerphan or @jeremiedbb, I'll try to debug numpy vs numpy-dev to track down were the difference comes from so that it would be less about guessing.

Labels: Bug
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
I do not have time to look at it, unfortunately.

Labels: Bug
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
May try to step through the code here to see the likely change.  I am still very confident this is just a tiny precision change (very likely to the better even) due to different chunk sizes getting summed up.
It might be nice to just find out where this is happening for curiosity.

Labels: Bug
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
OK, it is indeed the einsum, and unlike normal summation which has a better precision, this indeed has a worse one.

The reason, is that the einsum code works the following way:
```
def inner(in1, in2, out):
    accum = 0
    for i1, i2 in zip(in1, in2):
        accum += i1 * i2

    out[0] += accum
```
Now, before the change, that for loop would go in chunks of 8196 and then add the result, which means that `accum` is bounded by ~8196.  But after the change, the for loop will go to the full 100_000.
That is good for `np.sum` because it uses pairwise summation, but worse here, because any chunking gives a bit of the pairwise summation advantage.

While slightly worse for `np.einsum("i,i->i")` we could just force that chunking in einsum for the "precision benefit".

Labels: Bug
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
Still not sure if the precision ask here is just a bit high (due to catastrophic cancellation that the docs warn about) considering that this already uses 64bit floats internally, but the behavior is "reverted" in NumPy now.
(In fact, if it was using 32bit floats internally, the NumPy change would have improved things a lot, because einsum uses a 64bit temporary accumulation.)

The precision here also potentially depends on the memory layout if X and Y are 2-dimensional...

Labels: Bug
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
## CI is no longer failing! ✅

[Successful run](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73041&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf) on Dec 23, 2024

Labels: Bug
⚠️ CI failed on Linux_Nightly.pylatest_pip_scipy_dev (last failure: Dec 22, 2024) ⚠️

**CI is still failing on [Linux_Nightly.pylatest_pip_scipy_dev](https://dev.azure.com/scikit-learn/scikit-learn/_build/results?buildId=73034&view=logs&j=dfe99b15-50db-5d7b-b1e9-4105c42527cf)** (Dec 22, 2024)
- test_euclidean_distances_extreme_values[1000000-float32-0.0001-1e-05]

COMMENT:
Thanks @seberg. Since the CI is green again after the merge of numpy/numpy#28043, I think we can close.

> Still not sure if the precision ask here is just a bit high (...)
> The precision here also potentially depends on the memory layout if X and Y are 2-dimensional...

Let's keep that in mind in case future changes cause this test to fail again.

Labels: Bug
Sensitivity Analysis with Random Forest Moel

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/19112

<div type='discussions-op-text'>

<sup>Originally posted by **lesteve** January  5, 2021</sup>
## 👋 Welcome!
  
We’re using Discussions as a place to connect with other members of our community. We hope that you:
  * Ask questions you’re wondering about.
  * Share ideas.
  * Engage with other community members.
  * Welcome others and are open-minded. Remember that this is a community we
  build together 💪.

To get started, comment below with an introduction of yourself and tell us about what you do with this community.

Note: we enabled the Github Discussions feature experimentally. We will be monitoring it and evaluating how it goes.</div>

COMMENT:
My input dataset includes total precipitation, snow water equivalence (SWE), soil moisture, soil temperature, air temperature, and stream water level. The GWL and input features appear to show some seasonal trends

Labels: Needs Triage
Sensitivity Analysis with Random Forest Moel

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/19112

<div type='discussions-op-text'>

<sup>Originally posted by **lesteve** January  5, 2021</sup>
## 👋 Welcome!
  
We’re using Discussions as a place to connect with other members of our community. We hope that you:
  * Ask questions you’re wondering about.
  * Share ideas.
  * Engage with other community members.
  * Welcome others and are open-minded. Remember that this is a community we
  build together 💪.

To get started, comment below with an introduction of yourself and tell us about what you do with this community.

Note: we enabled the Github Discussions feature experimentally. We will be monitoring it and evaluating how it goes.</div>

COMMENT:
Hi everyone, I am a graduate student in Geophysics and my research is on investigating factors driving groundwater depletion. I am currently working on a sensitivity analysis to identify the main drivers of groundwater depletion using Random Forest. My input dataset includes total precipitation, snow water equivalence (SWE), soil moisture, soil temperature, air temperature, and stream water level. The GWL and input features appear to show some seasonal trends I have the following concerns:1. Can the RF regression model deal with trended inputs (X) and outputs (Y)?  Can RF model non-stationary data, e.g., Y=GWL of bedrock well is trending downward?  If not, is de-trending of Y needed as a preprocessing step before input to the RF model?  Or, since your RF model is fit to existing data, without making a forecast (a new inputs X --> a new Y), do we need to worry about the non-stationary Y?  2) X exhibit different magnitudes:  preprocessing by scaling first before input to the RF model?  Should each input feature be scaled to 0-1 before training/testing? What splitting technique is recommended for this kind of analysis?

Labels: Needs Triage
Sensitivity Analysis with Random Forest Moel

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/19112

<div type='discussions-op-text'>

<sup>Originally posted by **lesteve** January  5, 2021</sup>
## 👋 Welcome!
  
We’re using Discussions as a place to connect with other members of our community. We hope that you:
  * Ask questions you’re wondering about.
  * Share ideas.
  * Engage with other community members.
  * Welcome others and are open-minded. Remember that this is a community we
  build together 💪.

To get started, comment below with an introduction of yourself and tell us about what you do with this community.

Note: we enabled the Github Discussions feature experimentally. We will be monitoring it and evaluating how it goes.</div>

COMMENT:
It seems like you are asking a generic machine learning question. The issue tracker is for scikit-learn bugs and new features. Maybe try [Cross-Validated](https://stats.stackexchange.com/) or a more appropriate StackExchange forum?

Unfortunately, we need to close this issue as this issue tracker is a communication tool used for the development of scikit-learn.

Labels: Needs Triage
Sensitivity Analysis with Random Forest Moel

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/19112

<div type='discussions-op-text'>

<sup>Originally posted by **lesteve** January  5, 2021</sup>
## 👋 Welcome!
  
We’re using Discussions as a place to connect with other members of our community. We hope that you:
  * Ask questions you’re wondering about.
  * Share ideas.
  * Engage with other community members.
  * Welcome others and are open-minded. Remember that this is a community we
  build together 💪.

To get started, comment below with an introduction of yourself and tell us about what you do with this community.

Note: we enabled the Github Discussions feature experimentally. We will be monitoring it and evaluating how it goes.</div>

COMMENT:
For completeness this kind of questions is fine for [Github Discussions](https://github.com/scikit-learn/scikit-learn/discussions) as well, maybe this is what you tried to do but you ended up creating an issue instead somehow ...

Labels: Needs Triage