Repository: scikit-learn/scikit-learn
Pull request: Fix prevent infinite loop in KMeans
Fixes https://github.com/scikit-learn/scikit-learn/issues/28055

When there are more centers than non-duplicate samples, the relocation of empty cluster will trigger at every iteration preventing the algorithm to stop.

Problem statement:
Infinite Loop in K-means when relocating empty clusters
### Describe the bug

Relocating empty clusters in Kmeans is not working as expected in this edge case, where :
- There is duplicate entries.
- The number of clusters is equal to the number of entries.
- Very particular initial positions.

Kmeans is stuck in a infinite loop, and the only way to end it is the max number of iterations. 

I have suggestions to improve the relocation of empty clusters : 
- When selecting an entry to fill an empty cluster (the entry that has max_dist), we should not select an entry in a cluster that contains only 1 entry.
- After each relocation of an empty clusters, Centers needs to be updated.

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.cluster import KMeans

X = np.array([
    [100, 100], [100, 102], [106, 100], [107, 101], [107, 101], [108, 100]
])
init = np.array([
    [98, 104], [103, 100], [109, 102], [110, 100], [111, 100], [112, 100]
])
kmeans = KMeans(
    n_clusters=6,
    init=init,
    random_state=0,
    verbose=True,
    algorithm='lloyd',
    max_iter=10_000
).fit(X)
```

### Expected Results

We should only perform 1 iteration

```
Initialization complete
Iteration 0, inertia 40.0.
Iteration 1, inertia 0.0.
```

### Actual Results

Kmeans is stuck in an infinite loop

```
Initialization complete
Iteration 0, inertia 40.0.
Iteration 1, inertia 0.0.
Iteration 2, inertia 0.0.
Iteration 3, inertia 0.0.
Iteration 4, inertia 0.0.
Iteration 5, inertia 0.0.
Iteration 6, inertia 0.0.
Iteration 7, inertia 0.0.
Iteration 8, inertia 0.0.
Iteration 9, inertia 0.0.
Iteration 10, inertia 0.0.
Iteration 11, inertia 0.0.
Iteration 12, inertia 0.0.
Iteration 13, inertia 0.0.
Iteration 14, inertia 0.0.
Iteration 15, inertia 0.0.
Iteration 16, inertia 0.0.
Iteration 17, inertia 0.0.
Iteration 18, inertia 0.0.
Iteration 19, inertia 0.0.
```

### Versions

```shell
System:
   python: 3.9.0
   machine: Windows-10

Python dependencies:
      sklearn: 1.3.2
          pip: 20.2.3
   setuptools: 49.2.1
        numpy: 1.25.1
        scipy: 1.11.4
       Cython: None
       pandas: 2.0.3
   matplotlib: 3.7.2
       joblib: 1.3.2
threadpoolctl: 3.2.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 12
         prefix: libopenblas
        version: 0.3.23
threading_layer: pthreads

       user_api: openmp
   internal_api: openmp
    num_threads: 12
         prefix: vcomp
        version: None

       user_api: blas
   internal_api: openblas
    num_threads: 12
         prefix: libopenblas
        version: 0.3.21.dev
threading_layer: pthreads
```

Hints from issue discussion:
> Kmeans is stuck in a infinite loop, and the only way to end it is the max number of iterations.

Hi @brmel, I used the code example you provided, except I omitted the `max_iterations` parameter:

```
import numpy as np
from sklearn.cluster import KMeans

X = np.array([[100, 100], [100, 102], [106, 100], [107, 101], [107, 101], [108, 100]])
InitialPos = np.array([[98, 104], [103, 100], [109,102], [110,100], [111,100], [112,100]])
kmeans = KMeans(n_clusters=6, init=InitialPos, random_state=0, verbose=True, algorithm='lloyd').fit(X)
```

In this case, Kmeans isn't stuck in an infinite loop, but throws a ConvergenceWarning after 300 iterations:

```
/usr/local/lib/python3.10/dist-packages/sklearn/cluster/_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
/usr/local/lib/python3.10/dist-packages/sklearn/cluster/_kmeans.py:1416: RuntimeWarning: Explicit initial center position passed: performing only one init in KMeans instead of n_init=10.
  super()._check_params_vs_input(X, default_n_init=10)
Initialization complete
Iteration 0, inertia 40.0.
Iteration 1, inertia 0.0.
Iteration 2, inertia 0.0.
...
Iteration 297, inertia 0.0.
Iteration 298, inertia 0.0.
Iteration 299, inertia 0.0.
/usr/local/lib/python3.10/dist-packages/sklearn/base.py:1152: ConvergenceWarning: Number of distinct clusters (5) found smaller than n_clusters (6). Possibly due to duplicate points in X.
  return fit_method(estimator, *args, **kwargs)
```
Hello @raj-pulapakura , 
Thank you for your response. 
In your case, Kmeans stopped at 300, because the default value for `max_iterations` is 300.
I believe Kmeans should stop in the second iteration because inertia is 0.0
I named the issue "infinite loop" because if someone want to reach the optimal solution and set `max_iterations `to a very big number (or inf), Kmeans will be stuck in an infinite loop. 
@brmel omds that makes me look so dumb 😅. 

Yeah, you make a valid point that the algorithm should stop when inertia becomes 0, it doesn't make sense to keep going.
Can we stop the algorithm when we get inertia 0? could this be a viable solution?
> Can we stop the algorithm when we get inertia 0? could this be a viable solution?

A simple `if` statement should do the trick, right?
We can try
Theoretically, the algorithm that does the relocation should stop automatically when inertia is equal to 0 and no need to add an if condition for that. 
The reason the algorithm is broken is mentioned in the description above (we fill an empty cluster by generating another empty cluster), and adding an if condition will only hide the bug.
Also, I don't know if we can reproduce another infinite loop even when inertia is not equal to 0. 
Indeed this is a bug. Thanks @brmel for the minimal reproducer.
The fact that you duplicated samples makes that we are in the ill-conditioned case where `n_clusters > n_samples`. Normally, we raise an error for this case so I'm not surprise that the algorithm does not work properly.

Then, I'm not sure if we can easily detect it without making a very expensive check to check for sample duplication.
@brmel I am entirely sure to understand what you are proposing here.

> When selecting an entry to fill an empty cluster (the entry that has max_dist), we should not select an entry in a cluster that contains only 1 entry.
> After each relocation of an empty clusters, Centers needs to be updated.

Debugging the `lloyd` loop, I can see that we have an "empty" centroid (the centroid for which we don't affect any sample) and this centroid will be changing label at each iteration which make the strict convergence fail. The non-strict convergence is also failing because between 2 iterations the change of label of the "empty" centroid implies a centroid shift.

I am not yet understanding what is the reason that the "empty" centroid is changing label (I thought that it was an `openmp` behaviour but the number of chunk is too low to even parallelize). @jeremiedbb will know for sure. I am pretty sure this is a corner case that we should be able to fix.

Commit message: Merge branch 'fix-infinite-loop-kmeans' of github.com:jeremiedbb/scikit-learn into fix-infinite-loop-kmeans

Patch to review:
diff --git a/doc/whats_new/v1.4.rst b/doc/whats_new/v1.4.rst
index c674a8619e076..fa5e6c094796c 100644
--- a/doc/whats_new/v1.4.rst
+++ b/doc/whats_new/v1.4.rst
@@ -20,6 +20,9 @@ Changelog
   :pr:`28121` by :user:`Pietro Peterlongo <pietroppeter>` and
   :user:`Yao Xiao <Charlie-XIAO>`.
 
+- |Fix| Avoid inifinte loop in :class:`cluster.KMeans` when the number of clusters is
+  larger than the number of non-duplicate samples.
+  :pr:`28165` by :user:`Jérémie du Boisberranger <jeremiedbb>`.
 
 .. _changes_1_4:
 
diff --git a/sklearn/cluster/_k_means_common.pyx b/sklearn/cluster/_k_means_common.pyx
index 151af55076b7b..7c9c1bb54eaae 100644
--- a/sklearn/cluster/_k_means_common.pyx
+++ b/sklearn/cluster/_k_means_common.pyx
@@ -192,6 +192,11 @@ cpdef void _relocate_empty_clusters_dense(
         int new_cluster_id, old_cluster_id, far_idx, idx, k
         floating weight
 
+    if np.max(distances) == 0:
+        # Happens when there are more clusters than non-duplicate samples. Relocating
+        # is pointless in this case.
+        return
+
     for idx in range(n_empty):
 
         new_cluster_id = empty_clusters[idx]
@@ -241,6 +246,11 @@ cpdef void _relocate_empty_clusters_sparse(
             X_indices[X_indptr[i]: X_indptr[i + 1]],
             centers_old[j], centers_squared_norms[j], True)
 
+    if np.max(distances) == 0:
+        # Happens when there are more clusters than non-duplicate samples. Relocating
+        # is pointless in this case.
+        return
+
     cdef:
         int[::1] far_from_centers = np.argpartition(distances, -n_empty)[:-n_empty-1:-1].astype(np.int32)
 
@@ -274,12 +284,18 @@ cdef void _average_centers(
         int n_features = centers.shape[1]
         int j, k
         floating alpha
+        int argmax_weight = np.argmax(weight_in_clusters)
 
     for j in range(n_clusters):
         if weight_in_clusters[j] > 0:
             alpha = 1.0 / weight_in_clusters[j]
             for k in range(n_features):
                 centers[j, k] *= alpha
+        else:
+            # For convenience, we avoid setting empty clusters at the origin but place
+            # them at the location of the biggest cluster.
+            for k in range(n_features):
+                centers[j, k] = centers[argmax_weight, k]
 
 
 cdef void _center_shift(
diff --git a/sklearn/cluster/tests/test_k_means.py b/sklearn/cluster/tests/test_k_means.py
index 5b0c7ab9aace8..43e48e3aff594 100644
--- a/sklearn/cluster/tests/test_k_means.py
+++ b/sklearn/cluster/tests/test_k_means.py
@@ -1352,3 +1352,20 @@ def test_sample_weight_zero(init, global_random_seed):
     # (i.e. be at a distance=0 from it)
     d = euclidean_distances(X[::2], clusters_weighted)
     assert not np.any(np.isclose(d, 0))
+
+
+@pytest.mark.parametrize("array_constr", data_containers, ids=data_containers_ids)
+@pytest.mark.parametrize("algorithm", ["lloyd", "elkan"])
+def test_relocating_with_duplicates(algorithm, array_constr):
+    """Check that kmeans stops when there are more centers than non-duplicate samples
+
+    non-regression test for issue #28055
+    """
+    X = np.array([[0, 0], [1, 1], [1, 1], [1, 0], [0, 1]])
+    km = KMeans(n_clusters=5, init=X, algorithm=algorithm)
+
+    msg = r"Number of distinct clusters \(4\) found smaller than n_clusters \(5\)"
+    with pytest.warns(ConvergenceWarning, match=msg):
+        km.fit(array_constr(X))
+
+    assert km.n_iter_ == 1


Code review comments:
diff_hunk: @@ -20,6 +20,9 @@ Changelog
   :pr:`28121` by :user:`Pietro Peterlongo <pietroppeter>` and
   :user:`Yao Xiao <Charlie-XIAO>`.
 
+- |Fix| Avoid inifinte loop in :class:`cluster.KMeans` when the number of clusters is
path: doc/whats_new/v1.4.rst
text: @user1:

```suggestion
- |Fix| Avoid infinite loop in :class:`cluster.KMeans` when the number of clusters is
```