Repository: scikit-learn/scikit-learn
Pull request: DOC Add warning for using `"kmeans++"` init with high-dimensional data
<!--
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 #23566.

#### What does this implement/fix? Explain your changes.

As mentioned by @jeremiedbb in [this comment](https://github.com/scikit-learn/scikit-learn/issues/23566#issuecomment-1149962749), the `"kmeans++"` init picks centroids with proba inverse to square distance, which means it will pick points extremely isolated if there are any, which then stay their own centroids all along.

Instead of increasing the default value of `n_init`, we agreed on a small committee to warn the users to avoid using the `"kmeans++"` init with sparse high dimensional datasets.

#### Any other comments?

Maybe we should wait until #23528 is merged before merging the present PR.

Side effect: fixes small typo `k_init` -> `n_init` in `KMeans` docstring.

<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
`silhouette_score` fails whit `"k-means++"` init on very sparse data
### Describe the bug

Not sure this is a bug or something that just need to be better documented.

The `"k-means++"` init on very sparse data can select initial centroids than never get updated and this cause problems in the silhouette clustering evaluation with `sample_size=1000`.

The problem disappears when `sample_size` is increased or by setting particular random states.

### Steps/Code to Reproduce

```python
from sklearn.datasets import fetch_20newsgroups
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.cluster import KMeans
from sklearn import metrics

categories = [
    "alt.atheism",
    "talk.religion.misc",
    "comp.graphics",
    "sci.space",
]

dataset = fetch_20newsgroups(
    remove=("headers", "footers", "quotes"),
    subset="all",
    categories=categories,
    shuffle=True,
    random_state=42,
)

vectorizer = TfidfVectorizer(
    max_df=0.5,
    min_df=2,
    stop_words="english",
)

X = vectorizer.fit_transform(dataset.data)

km = KMeans(
    n_clusters=4,
    init="k-means++",
    max_iter=100,
    n_init=1,
    random_state=0,  #notice fixed random state
).fit(X)


metrics.silhouette_score(X, km.labels_, sample_size=1000, random_state=3) # works

metrics.silhouette_score(X, km.labels_, sample_size=1000, random_state=4) # fails
```

### Expected Results

Always return silhouette score.

### Actual Results

```python-traceback
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Untitled-1 in <module>
----> <a href='untitled:Untitled-1?line=37'>38</a> metrics.silhouette_score(X, km.labels_, sample_size=1000, random_state=43)

File ~/scikit-learn/sklearn/metrics/cluster/_unsupervised.py:117, in silhouette_score(X, labels, metric, sample_size, random_state, **kwds)
    [115](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=114)     else:
    [116](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=115)         X, labels = X[indices], labels[indices]
--> [117](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=116) return np.mean(silhouette_samples(X, labels, metric=metric, **kwds))

File ~/scikit-learn/sklearn/metrics/cluster/_unsupervised.py:231, in silhouette_samples(X, labels, metric, **kwds)
    [229](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=228) n_samples = len(labels)
    [230](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=229) label_freqs = np.bincount(labels)
--> [231](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=230) check_number_of_labels(len(le.classes_), n_samples)
    [233](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=232) kwds["metric"] = metric
    [234](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=233) reduce_func = functools.partial(
    [235](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=234)     _silhouette_reduce, labels=labels, label_freqs=label_freqs
    [236](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=235) )

File ~/scikit-learn/sklearn/metrics/cluster/_unsupervised.py:33, in check_number_of_labels(n_labels, n_samples)
     [22](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=21) """Check that number of labels are valid.
     [23](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=22) 
     [24](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=23) Parameters
   (...)
     [30](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=29)     Number of samples.
     [31](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=30) """
     [32](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=31) if not 1 < n_labels < n_samples:
---> [33](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=32)     raise ValueError(
     [34](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=33)         "Number of labels is %d. Valid values are 2 to n_samples - 1 (inclusive)"
     [35](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=34)         % n_labels
     [36](file:///home/arturoamor/scikit-learn/sklearn/metrics/cluster/_unsupervised.py?line=35)     )

ValueError: Number of labels is 1. Valid values are 2 to n_samples - 1 (inclusive)
```

### Versions

```shell
System:
    python: 3.9.9 | packaged by conda-forge | (main, Dec 20 2021, 02:40:17)  [GCC 9.4.0]
executable: /home/arturoamor/miniforge3/envs/dev-scikit-learn/bin/python
   machine: Linux-5.14.0-1038-oem-x86_64-with-glibc2.31

Python dependencies:
      sklearn: 1.2.dev0
          pip: 21.3.1
   setuptools: 60.5.0
        numpy: 1.22.1
        scipy: 1.7.3
       Cython: 0.29.26
       pandas: 1.3.5
   matplotlib: 3.5.1
       joblib: 1.1.0
threadpoolctl: 3.0.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /home/arturoamor/miniforge3/envs/dev-scikit-learn/lib/python3.9/site-packages/numpy.libs/libopenblas64_p-r0-2f7c42d4.3.18.so
        version: 0.3.18
threading_layer: pthreads
   architecture: SkylakeX
    num_threads: 8

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /home/arturoamor/miniforge3/envs/dev-scikit-learn/lib/python3.9/site-packages/scipy.libs/libopenblasp-r0-8b9e111f.3.17.so
        version: 0.3.17
threading_layer: pthreads
   architecture: SkylakeX
    num_threads: 8

       user_api: openmp
   internal_api: openmp
         prefix: libgomp
       filepath: /home/arturoamor/miniforge3/envs/dev-scikit-learn/lib/libgomp.so.1.0.0
        version: None
    num_threads: 8
```

Hints from issue discussion:
Since random_state is not set in KMeans in the reproducer, I can't reproduce it yet. We need to find a seed that trigger it

> Since random_state is not set in KMeans in the reproducer, I can't reproduce it yet. We need to find a seed that trigger it

I updated it.
Isn't it just bad luck ?

It just happens that it's selecting `sample_size` data points that all correspond to the same label ?

I can reproduce with this:
```py
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
from sklearn import metrics
 
X, _ = make_blobs(n_samples=10000, centers=2, n_features=2, random_state=0)
 
km = KMeans(n_clusters=2, init="random", n_init=1, random_state=0).fit(X)
 
metrics.silhouette_score(X, km.labels_, sample_size=5, random_state=48)
```

Here I set 2 clusters and only pick `5` samples to maximize the chance of picking all points corresponding to a same label.
In your examples, the clustering is highly non balanced: 
```py
np.unique(km.labels_, return_counts=True)
>>> (array([0, 1, 2, 3], dtype=int32), array([   1,    1, 3382,    3]))
```
Almost all points are in the same cluster, hence more chance to only pick points corresponding to a same label.
Looking at the results for different seeds of ``KMeans(n_clusters=4, init="k-means++", max_iter=100, n_init=1, random_state=seed).fit(X)``, this kind of unbalance happens sometimes. I guess it's because kmeans++ picks centroids with proba inverse to square distance, which means it will pick points extremely isolated if there are any, which then stay their own centroids all along.

Note that it doesn't mean it's a bad clustering. The inertia is close to the inertia obtained with a random init. It's just kind of  reflecting the structure of the dataset. 

It doesn't look like there's an issue with kmeans++. It's just that KMeans (euclidean distance in fact) is not really good for such high dimensional dataset.

Commit message: Merge branch 'main' of https://github.com/scikit-learn/scikit-learn into warn_kmeans

Patch to review:
diff --git a/sklearn/cluster/_kmeans.py b/sklearn/cluster/_kmeans.py
index 25cd8712c2f0c..1f77c8f8f7a67 100644
--- a/sklearn/cluster/_kmeans.py
+++ b/sklearn/cluster/_kmeans.py
@@ -1188,9 +1188,8 @@ class KMeans(_BaseKMeans):
             (n_clusters, n_features), default='k-means++'
         Method for initialization:
 
-        'k-means++' : selects initial cluster centers for k-mean
-        clustering in a smart way to speed up convergence. See section
-        Notes in k_init for more details.
+        'k-means++' : selects initial cluster centers for k-mean clustering in a
+        smart way to speed up convergence. See notes in `n_init` for more details.
 
         'random': choose `n_clusters` observations (rows) at random from data
         for the initial centroids.
@@ -1202,9 +1201,10 @@ class KMeans(_BaseKMeans):
         random state and return an initialization.
 
     n_init : 'auto' or int, default=10
-        Number of time the k-means algorithm will be run with different
-        centroid seeds. The final results will be the best output of
-        n_init consecutive runs in terms of inertia.
+        Number of times the k-means algorithm is run with different centroid
+        seeds. The final results is the best output of `n_init` consecutive runs
+        in terms of inertia. Several runs are recommended for sparse
+        high-dimensional problems (see :ref:`kmeans_sparse_high_dim`).
 
         When `n_init='auto'`, the number of runs will be 10 if using
         `init='random'`, and 1 if using `init='kmeans++'`.
@@ -1663,9 +1663,8 @@ class MiniBatchKMeans(_BaseKMeans):
             (n_clusters, n_features), default='k-means++'
         Method for initialization:
 
-        'k-means++' : selects initial cluster centers for k-mean
-        clustering in a smart way to speed up convergence. See section
-        Notes in k_init for more details.
+        'k-means++' : selects initial cluster centers for k-mean clustering in a
+        smart way to speed up convergence. See notes in `n_init` for more details.
 
         'random': choose `n_clusters` observations (rows) at random from data
         for the initial centroids.
@@ -1729,8 +1728,10 @@ class MiniBatchKMeans(_BaseKMeans):
 
     n_init : 'auto' or int, default=3
         Number of random initializations that are tried.
-        In contrast to KMeans, the algorithm is only run once, using the
-        best of the ``n_init`` initializations as measured by inertia.
+        In contrast to KMeans, the algorithm is only run once, using the best of
+        the `n_init` initializations as measured by inertia. Several runs are
+        recommended for sparse high-dimensional problems (see
+        :ref:`kmeans_sparse_high_dim`).
 
         When `n_init='auto'`, the number of runs will be 3 if using
         `init='random'`, and 1 if using `init='kmeans++'`.


Code review comments:
diff_hunk: @@ -1663,9 +1663,8 @@ class MiniBatchKMeans(_BaseKMeans):
             (n_clusters, n_features), default='k-means++'
         Method for initialization:
 
-        'k-means++' : selects initial cluster centers for k-mean
-        clustering in a smart way to speed up convergence. See section
-        Notes in k_init for more details.
+        'k-means++' : selects initial cluster centers for k-mean clustering in a
+        smart way to speed up convergence. See notes in `n_init` for more details.
path: sklearn/cluster/_kmeans.py
text: @user1:
Similarly:
```suggestion
        'k-means++' : selects initial cluster centroids using sampling 
        based on a empirical probably distribution of the points
        contribution to the overall inertia.
        This technique speeds up convergence, and is theoretically
        proven to be :math:`\mathcal{O}(\log k)`-optimal.
        See the description of `n_init` for more details.
```