Repository: scikit-learn/scikit-learn
Issue category: Performance Issue

Issue:
KMeans optimisation for array C/F contiguity (was: Make sure that the output of PCA.fit_transform is C contiguous)
otherwise, I would rather use:

``` Python
pca.fit(X)
X_new = pca.transform(X)
```

Because of FORTRAN data for inner product is very slow. (for example, in the KMeans)

KMeans optimisation for array C/F contiguity (was: Make sure that the output of PCA.fit_transform is C contiguous)
otherwise, I would rather use:

``` Python
pca.fit(X)
X_new = pca.transform(X)
```

Because of FORTRAN data for inner product is very slow. (for example, in the KMeans)



Hints:
Can you give an example?
Maybe we should then rather change the behavior in KMeans.
I wouldn't change the output format, since we don't know what the user wants to do next.

They should be the same that the output and input format of PCA.fit_transform, isn't it?
Because PCA.transform is such.

In _k_means._assign_labels_array, the ddot will be very slow, when X is Fortran data.

``` Python
    for sample_idx in range(n_samples):
        min_dist = -1
        for center_idx in range(n_clusters):
            dist = 0.0
            # hardcoded: minimize euclidean distance to cluster center:
            # ||a - b||^2 = ||a||^2 + ||b||^2 -2 <a, b>
            dist += ddot(n_features, &X[sample_idx, 0], x_stride,
                         &centers[center_idx, 0], center_stride)
```

I have a large sample set, before the dimension reduction, each iteration only need a little more than a minute, but need to be 7 minute after dimension reduction.

2 important questions:
- How general is the requirement to be C contiguous? It speeds up for
  k-means, but will it not slow down for other algorithms?
- Can the change be made in k-means rather than in PCA? If you are not
  using copy_X=False, I don't see any argument for not making the change
  in k-means, and I would much prefer it.

- I just think it is a natural idea that the features of a sample are stored in a row.
- I don't know how to speeds up for k-means, in the case of using Fortran data. :(
- copy_X=False. This reduces the memory usage, but without speeds up for k-means. Anyway I like it too.

>   • I just think it is a natural idea that the features of a sample are stored in a row.

I think that you are assuming too much.

>   • I don't know how to speeds up for k-means, in the case of using Fortran
>     data. :(

It's just a question of calling "np.ascontiguousarray" on the array when
making the copy.

>   • copy_X=False. This reduces the memory usage,

OK. We would have to do a copy here.

I am opposed to changing the behavior of PCA without a more exhausive
review of what it implies speedwise in the codebase: what algorithms
perform best with C or Fortran ordered data.

Ok, I agree with you, and thank you for your explanation.

@zhaipro I think it would be appreciated if you could benchmark k-means with different memory layouts and see how it performs. We could change the memory layout there, which would solve your problem. Btw, you might be interested in this PR #2008 which will make k-means quite a bit faster.

benchmark:
C order time: 1.370000 inertia: 422652.759578
F order time: 6.536000 inertia: 422652.759578

I have a large sample set, so that the `precompute_distances = False`, 
by the following code in k_means:

``` Python
if precompute_distances == 'auto':
    n_samples = X.shape[0]
    precompute_distances = (n_clusters * n_samples) < 12e6
```

Here, To show my problem, I set directly `precompute_distances = False`.

``` Python
import numpy as np
from sklearn.cluster import KMeans
from time import time

def bench_kmeans(name, data):
    start = time()
    km = KMeans(n_clusters=200, init='random', n_init=1, max_iter=1,
                copy_x=False, random_state=42, precompute_distances=False).fit(data)
    print("%s time: %f inertia: %f" % (name, time() - start, km.inertia_))


np.random.seed(0)
data = np.random.random(3000*1000).reshape((3000, 1000))
# for C order
bench_kmeans(name='C order', data=data)
# for F order
data = np.asfortranarray(data)
bench_kmeans(name='F order', data=data)
```

Is that true for slim and fat data and different number of clusters? Also, have you tried the elkan branch?

Cool, for elkan alg.
@amueller However, it requires more memory?

```
('shape:', (3000, 1000))
n_clusters  alg name    time    inertia:
50  lloyd   C order 0.426000    453724.729456
50  lloyd   F order 1.782000    453724.729456
50  elkan   C order 0.588000    244227.752545
50  elkan   F order 0.652000    244227.752545
100 lloyd   C order 0.744000    442351.514140
100 lloyd   F order 3.488000    442351.514140
100 elkan   C order 0.892000    239220.149899
100 elkan   F order 0.942000    239220.149899
200 lloyd   C order 1.401000    422652.759578
200 lloyd   F order 6.694000    422652.759578
200 elkan   C order 1.608000    229660.875075
200 elkan   F order 1.632000    229660.875075
('shape:', (1000, 3000))
n_clusters  alg name    time    inertia:
50  lloyd   C order 0.510000    453021.642152
50  lloyd   F order 1.416000    453021.642152
50  elkan   C order 0.630000    236247.135296
50  elkan   F order 0.692000    236247.135296
100 lloyd   C order 1.042000    427321.624518
100 lloyd   F order 2.580000    427321.624518
100 elkan   C order 1.152000    222766.273428
100 elkan   F order 1.212000    222766.273428
200 lloyd   C order 1.908000    378199.959299
200 lloyd   F order 5.046000    378199.959299
200 elkan   C order 1.964000    196655.300444
200 elkan   F order 2.069000    196655.300444
```

``` Python
import numpy as np
from sklearn.cluster import KMeans
from time import time
from itertools import product

def bench_kmeans(name, data, alg, n_clusters):
    start = time()
    km = KMeans(algorithm=alg, n_clusters=n_clusters, init='random', n_init=1, max_iter=1,
                copy_x=False, random_state=42, precompute_distances=False).fit(data)
    print("%d\t%s\t%s\t%f\t%f" % (n_clusters, alg, name, time() - start, km.inertia_))


def test_kmeans(data):
    c_data = data
    f_data = np.asfortranarray(data)
    print('n_clusters\talg\tname\ttime\tinertia:')
    for n_clusters, alg, (name, data) in product((50, 100, 200),
                                    ('lloyd', 'elkan'),
                                     zip(('C order', 'F order'), (c_data, f_data))):
        bench_kmeans(name=name, data=data, alg=alg, n_clusters=n_clusters)

np.random.seed(0)
data = np.random.random(3000*1000).reshape((3000, 1000))
print('shape:', data.shape)
test_kmeans(data)
data = data.reshape((1000, 3000))
print('shape:', data.shape)
test_kmeans(data)
```

Thanks a lot for the benchmark.
I'm quite surprised at the difference in inertia between the two implementations... maybe I messed up the stopping criterion again?
The memory requirement for elkan shouldn't be substantial.

It looks to me like for lloyd (the current implementation) we should `copy("C")` by default (or rather do that in `check_array`
I'm not sure if we should make this optional or add a parameter. I think I'd actually not add a parameter but just add to the docs "the data will be converted to C ordering, which might cause a memory copy if the data is given in fortran order" something like that.

> I'm quite surprised at the difference in inertia between the two implementations...

@amueller I think maybe the reason is that I have only one iteration.

makes sense.

Can someone give me a little explanation of what needs to be done for this issue? Since it seems like the issue description (of title was later updated), perhaps the issue description was not updated.

In particular, there is term ['C contiguous'](https://www.ibm.com/support/knowledgecenter/SSAT4T_14.1.0/com.ibm.xlf141.linux.doc/language_ref/contiguous.html) (probably right link) what does that refer, is that fortran term? machine learning term? or [numpy.ascountiguous](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.ascontiguousarray.html).

Also if a summary of issue discussion can be given it had be great (since following the issue discussion is a little confusing for me). I'm assuming that this issue needs work, even though the PR https://github.com/scikit-learn/scikit-learn/pull/5414 was merged, which seems to be for accelerating K means algorithm.
Summary: KMeans, when using check_array, should declare order='C'.

See https://docs.scipy.org/doc/numpy-1.12.0/glossary.html#term-row-major
Can you give an example?
Maybe we should then rather change the behavior in KMeans.
I wouldn't change the output format, since we don't know what the user wants to do next.

They should be the same that the output and input format of PCA.fit_transform, isn't it?
Because PCA.transform is such.

In _k_means._assign_labels_array, the ddot will be very slow, when X is Fortran data.

``` Python
    for sample_idx in range(n_samples):
        min_dist = -1
        for center_idx in range(n_clusters):
            dist = 0.0
            # hardcoded: minimize euclidean distance to cluster center:
            # ||a - b||^2 = ||a||^2 + ||b||^2 -2 <a, b>
            dist += ddot(n_features, &X[sample_idx, 0], x_stride,
                         &centers[center_idx, 0], center_stride)
```

I have a large sample set, before the dimension reduction, each iteration only need a little more than a minute, but need to be 7 minute after dimension reduction.

2 important questions:
- How general is the requirement to be C contiguous? It speeds up for
  k-means, but will it not slow down for other algorithms?
- Can the change be made in k-means rather than in PCA? If you are not
  using copy_X=False, I don't see any argument for not making the change
  in k-means, and I would much prefer it.

- I just think it is a natural idea that the features of a sample are stored in a row.
- I don't know how to speeds up for k-means, in the case of using Fortran data. :(
- copy_X=False. This reduces the memory usage, but without speeds up for k-means. Anyway I like it too.

>   • I just think it is a natural idea that the features of a sample are stored in a row.

I think that you are assuming too much.

>   • I don't know how to speeds up for k-means, in the case of using Fortran
>     data. :(

It's just a question of calling "np.ascontiguousarray" on the array when
making the copy.

>   • copy_X=False. This reduces the memory usage,

OK. We would have to do a copy here.

I am opposed to changing the behavior of PCA without a more exhausive
review of what it implies speedwise in the codebase: what algorithms
perform best with C or Fortran ordered data.

Ok, I agree with you, and thank you for your explanation.

@zhaipro I think it would be appreciated if you could benchmark k-means with different memory layouts and see how it performs. We could change the memory layout there, which would solve your problem. Btw, you might be interested in this PR #2008 which will make k-means quite a bit faster.

benchmark:
C order time: 1.370000 inertia: 422652.759578
F order time: 6.536000 inertia: 422652.759578

I have a large sample set, so that the `precompute_distances = False`, 
by the following code in k_means:

``` Python
if precompute_distances == 'auto':
    n_samples = X.shape[0]
    precompute_distances = (n_clusters * n_samples) < 12e6
```

Here, To show my problem, I set directly `precompute_distances = False`.

``` Python
import numpy as np
from sklearn.cluster import KMeans
from time import time

def bench_kmeans(name, data):
    start = time()
    km = KMeans(n_clusters=200, init='random', n_init=1, max_iter=1,
                copy_x=False, random_state=42, precompute_distances=False).fit(data)
    print("%s time: %f inertia: %f" % (name, time() - start, km.inertia_))


np.random.seed(0)
data = np.random.random(3000*1000).reshape((3000, 1000))
# for C order
bench_kmeans(name='C order', data=data)
# for F order
data = np.asfortranarray(data)
bench_kmeans(name='F order', data=data)
```

Is that true for slim and fat data and different number of clusters? Also, have you tried the elkan branch?

Cool, for elkan alg.
@amueller However, it requires more memory?

```
('shape:', (3000, 1000))
n_clusters  alg name    time    inertia:
50  lloyd   C order 0.426000    453724.729456
50  lloyd   F order 1.782000    453724.729456
50  elkan   C order 0.588000    244227.752545
50  elkan   F order 0.652000    244227.752545
100 lloyd   C order 0.744000    442351.514140
100 lloyd   F order 3.488000    442351.514140
100 elkan   C order 0.892000    239220.149899
100 elkan   F order 0.942000    239220.149899
200 lloyd   C order 1.401000    422652.759578
200 lloyd   F order 6.694000    422652.759578
200 elkan   C order 1.608000    229660.875075
200 elkan   F order 1.632000    229660.875075
('shape:', (1000, 3000))
n_clusters  alg name    time    inertia:
50  lloyd   C order 0.510000    453021.642152
50  lloyd   F order 1.416000    453021.642152
50  elkan   C order 0.630000    236247.135296
50  elkan   F order 0.692000    236247.135296
100 lloyd   C order 1.042000    427321.624518
100 lloyd   F order 2.580000    427321.624518
100 elkan   C order 1.152000    222766.273428
100 elkan   F order 1.212000    222766.273428
200 lloyd   C order 1.908000    378199.959299
200 lloyd   F order 5.046000    378199.959299
200 elkan   C order 1.964000    196655.300444
200 elkan   F order 2.069000    196655.300444
```

``` Python
import numpy as np
from sklearn.cluster import KMeans
from time import time
from itertools import product

def bench_kmeans(name, data, alg, n_clusters):
    start = time()
    km = KMeans(algorithm=alg, n_clusters=n_clusters, init='random', n_init=1, max_iter=1,
                copy_x=False, random_state=42, precompute_distances=False).fit(data)
    print("%d\t%s\t%s\t%f\t%f" % (n_clusters, alg, name, time() - start, km.inertia_))


def test_kmeans(data):
    c_data = data
    f_data = np.asfortranarray(data)
    print('n_clusters\talg\tname\ttime\tinertia:')
    for n_clusters, alg, (name, data) in product((50, 100, 200),
                                    ('lloyd', 'elkan'),
                                     zip(('C order', 'F order'), (c_data, f_data))):
        bench_kmeans(name=name, data=data, alg=alg, n_clusters=n_clusters)

np.random.seed(0)
data = np.random.random(3000*1000).reshape((3000, 1000))
print('shape:', data.shape)
test_kmeans(data)
data = data.reshape((1000, 3000))
print('shape:', data.shape)
test_kmeans(data)
```

Thanks a lot for the benchmark.
I'm quite surprised at the difference in inertia between the two implementations... maybe I messed up the stopping criterion again?
The memory requirement for elkan shouldn't be substantial.

It looks to me like for lloyd (the current implementation) we should `copy("C")` by default (or rather do that in `check_array`
I'm not sure if we should make this optional or add a parameter. I think I'd actually not add a parameter but just add to the docs "the data will be converted to C ordering, which might cause a memory copy if the data is given in fortran order" something like that.

> I'm quite surprised at the difference in inertia between the two implementations...

@amueller I think maybe the reason is that I have only one iteration.

makes sense.

Can someone give me a little explanation of what needs to be done for this issue? Since it seems like the issue description (of title was later updated), perhaps the issue description was not updated.

In particular, there is term ['C contiguous'](https://www.ibm.com/support/knowledgecenter/SSAT4T_14.1.0/com.ibm.xlf141.linux.doc/language_ref/contiguous.html) (probably right link) what does that refer, is that fortran term? machine learning term? or [numpy.ascountiguous](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.ascontiguousarray.html).

Also if a summary of issue discussion can be given it had be great (since following the issue discussion is a little confusing for me). I'm assuming that this issue needs work, even though the PR https://github.com/scikit-learn/scikit-learn/pull/5414 was merged, which seems to be for accelerating K means algorithm.
Summary: KMeans, when using check_array, should declare order='C'.

See https://docs.scipy.org/doc/numpy-1.12.0/glossary.html#term-row-major

Patch:
diff --git a/doc/whats_new/v0.20.rst b/doc/whats_new/v0.20.rst
index a7ae9f7415243..b4a8046d02956 100644
--- a/doc/whats_new/v0.20.rst
+++ b/doc/whats_new/v0.20.rst
@@ -125,6 +125,13 @@ Classifiers and regressors
   only require X to be an object with finite length or shape.
   :issue:`9832` by :user:`Vrishank Bhardwaj <vrishank97>`.
 
+Cluster
+
+- :class:`cluster.KMeans`, :class:`cluster.MiniBatchKMeans` and
+  :func:`cluster.k_means` passed with ``algorithm='full'`` now enforces
+  row-major ordering, improving runtime.
+  :issue:`10471` by :user:`Gaurav Dhingra <gxyd>`.
+
 Preprocessing
 
 - :class:`preprocessing.PolynomialFeatures` now supports sparse input.
diff --git a/sklearn/cluster/k_means_.py b/sklearn/cluster/k_means_.py
index 07e390701eaac..5b6866b61d748 100644
--- a/sklearn/cluster/k_means_.py
+++ b/sklearn/cluster/k_means_.py
@@ -22,6 +22,7 @@
 from ..utils.extmath import row_norms, squared_norm, stable_cumsum
 from ..utils.sparsefuncs_fast import assign_rows_csr
 from ..utils.sparsefuncs import mean_variance_axis
+from ..utils.validation import _num_samples
 from ..utils import check_array
 from ..utils import check_random_state
 from ..utils import as_float_array
@@ -175,7 +176,9 @@ def k_means(X, n_clusters, init='k-means++', precompute_distances='auto',
     Parameters
     ----------
     X : array-like or sparse matrix, shape (n_samples, n_features)
-        The observations to cluster.
+        The observations to cluster. It must be noted that the data
+        will be converted to C ordering, which will cause a memory copy
+        if the given data is not C-contiguous.
 
     n_clusters : int
         The number of clusters to form as well as the number of
@@ -230,10 +233,12 @@ def k_means(X, n_clusters, init='k-means++', precompute_distances='auto',
 
     copy_x : boolean, optional
         When pre-computing distances it is more numerically accurate to center
-        the data first.  If copy_x is True, then the original data is not
-        modified.  If False, the original data is modified, and put back before
-        the function returns, but small numerical differences may be introduced
-        by subtracting and then adding the data mean.
+        the data first.  If copy_x is True (default), then the original data is
+        not modified, ensuring X is C-contiguous.  If False, the original data
+        is modified, and put back before the function returns, but small
+        numerical differences may be introduced by subtracting and then adding
+        the data mean, in this case it will also not ensure that data is
+        C-contiguous which may cause a significant slowdown.
 
     n_jobs : int
         The number of jobs to use for the computation. This works by computing
@@ -280,7 +285,14 @@ def k_means(X, n_clusters, init='k-means++', precompute_distances='auto',
         raise ValueError('Number of iterations should be a positive number,'
                          ' got %d instead' % max_iter)
 
-    X = as_float_array(X, copy=copy_x)
+    # avoid forcing order when copy_x=False
+    order = "C" if copy_x else None
+    X = check_array(X, accept_sparse='csr', dtype=[np.float64, np.float32],
+                    order=order, copy=copy_x)
+    # verify that the number of samples given is larger than k
+    if _num_samples(X) < n_clusters:
+        raise ValueError("n_samples=%d should be >= n_clusters=%d" % (
+            _num_samples(X), n_clusters))
     tol = _tolerance(X, tol)
 
     # If the distances are precomputed every job will create a matrix of shape
@@ -392,8 +404,7 @@ def _kmeans_single_elkan(X, n_clusters, max_iter=300, init='k-means++',
                          random_state=None, tol=1e-4,
                          precompute_distances=True):
     if sp.issparse(X):
-        raise ValueError("algorithm='elkan' not supported for sparse input X")
-    X = check_array(X, order="C")
+        raise TypeError("algorithm='elkan' not supported for sparse input X")
     random_state = check_random_state(random_state)
     if x_squared_norms is None:
         x_squared_norms = row_norms(X, squared=True)
@@ -768,12 +779,14 @@ class KMeans(BaseEstimator, ClusterMixin, TransformerMixin):
         If None, the random number generator is the RandomState instance used
         by `np.random`.
 
-    copy_x : boolean, default True
+    copy_x : boolean, optional
         When pre-computing distances it is more numerically accurate to center
-        the data first.  If copy_x is True, then the original data is not
-        modified.  If False, the original data is modified, and put back before
-        the function returns, but small numerical differences may be introduced
-        by subtracting and then adding the data mean.
+        the data first.  If copy_x is True (default), then the original data is
+        not modified, ensuring X is C-contiguous.  If False, the original data
+        is modified, and put back before the function returns, but small
+        numerical differences may be introduced by subtracting and then adding
+        the data mean, in this case it will also not ensure that data is
+        C-contiguous which may cause a significant slowdown.
 
     n_jobs : int
         The number of jobs to use for the computation. This works by computing
@@ -860,14 +873,6 @@ def __init__(self, n_clusters=8, init='k-means++', n_init=10,
         self.n_jobs = n_jobs
         self.algorithm = algorithm
 
-    def _check_fit_data(self, X):
-        """Verify that the number of samples given is larger than k"""
-        X = check_array(X, accept_sparse='csr', dtype=[np.float64, np.float32])
-        if X.shape[0] < self.n_clusters:
-            raise ValueError("n_samples=%d should be >= n_clusters=%d" % (
-                X.shape[0], self.n_clusters))
-        return X
-
     def _check_test_data(self, X):
         X = check_array(X, accept_sparse='csr', dtype=FLOAT_DTYPES)
         n_samples, n_features = X.shape
@@ -885,13 +890,14 @@ def fit(self, X, y=None):
         Parameters
         ----------
         X : array-like or sparse matrix, shape=(n_samples, n_features)
-            Training instances to cluster.
+            Training instances to cluster. It must be noted that the data
+            will be converted to C ordering, which will cause a memory
+            copy if the given data is not C-contiguous.
 
         y : Ignored
 
         """
         random_state = check_random_state(self.random_state)
-        X = self._check_fit_data(X)
 
         self.cluster_centers_, self.labels_, self.inertia_, self.n_iter_ = \
             k_means(
@@ -944,7 +950,6 @@ def fit_transform(self, X, y=None):
         # np.array or CSR format already.
         # XXX This skips _check_test_data, which may change the dtype;
         # we should refactor the input validation.
-        X = self._check_fit_data(X)
         return self.fit(X)._transform(X)
 
     def transform(self, X):
@@ -1351,7 +1356,9 @@ def fit(self, X, y=None):
         Parameters
         ----------
         X : array-like or sparse matrix, shape=(n_samples, n_features)
-            Training instances to cluster.
+            Training instances to cluster. It must be noted that the data
+            will be converted to C ordering, which will cause a memory copy
+            if the given data is not C-contiguous.
 
         y : Ignored
 
@@ -1361,8 +1368,8 @@ def fit(self, X, y=None):
                         dtype=[np.float64, np.float32])
         n_samples, n_features = X.shape
         if n_samples < self.n_clusters:
-            raise ValueError("Number of samples smaller than number "
-                             "of clusters.")
+            raise ValueError("n_samples=%d should be >= n_clusters=%d"
+                             % (n_samples, self.n_clusters))
 
         n_init = self.n_init
         if hasattr(self.init, '__array__'):
@@ -1516,13 +1523,14 @@ def partial_fit(self, X, y=None):
         Parameters
         ----------
         X : array-like, shape = [n_samples, n_features]
-            Coordinates of the data points to cluster.
+            Coordinates of the data points to cluster. It must be noted that
+            X will be copied if it is not C-contiguous.
 
         y : Ignored
 
         """
 
-        X = check_array(X, accept_sparse="csr")
+        X = check_array(X, accept_sparse="csr", order="C")
         n_samples, n_features = X.shape
         if hasattr(self.init, '__array__'):
             self.init = np.ascontiguousarray(self.init, dtype=X.dtype)


Test patch:
diff --git a/sklearn/cluster/tests/test_k_means.py b/sklearn/cluster/tests/test_k_means.py
index f0f2b56bd591c..de8772d761e22 100644
--- a/sklearn/cluster/tests/test_k_means.py
+++ b/sklearn/cluster/tests/test_k_means.py
@@ -169,7 +169,8 @@ def _check_fitted_model(km):
     assert_greater(km.inertia_, 0.0)
 
     # check error on dataset being too small
-    assert_raises(ValueError, km.fit, [[0., 1.]])
+    assert_raise_message(ValueError, "n_samples=1 should be >= n_clusters=%d"
+                         % km.n_clusters, km.fit, [[0., 1.]])
 
 
 def test_k_means_plus_plus_init():
@@ -750,6 +751,11 @@ def test_k_means_function():
     # to many clusters desired
     assert_raises(ValueError, k_means, X, n_clusters=X.shape[0] + 1)
 
+    # kmeans for algorithm='elkan' raises TypeError on sparse matrix
+    assert_raise_message(TypeError, "algorithm='elkan' not supported for "
+                         "sparse input X", k_means, X=X_csr, n_clusters=2,
+                         algorithm="elkan")
+
 
 def test_x_squared_norms_init_centroids():
     """Test that x_squared_norms can be None in _init_centroids"""