Repository: scikit-learn/scikit-learn
Pull request: [MRG + 1] KMeans optimisation for array C/F contiguity


<!--
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/master/CONTRIBUTING.md#pull-request-checklist
-->

#### 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 #4582

#### What does this implement/fix? Explain your changes.

Use `order='C'` internally in `check_array` for KMeans optimization.

#### Any other comments?

I don't think this will require any tests. Requires benchmarking? apart from what was described on the issue itself?

<!--
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:
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 from issue discussion:
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)
```

Commit message: KMeans optimisation for array C/F contiguity

Fixes #4582

Patch to review:
diff --git a/sklearn/cluster/k_means_.py b/sklearn/cluster/k_means_.py
index 07e390701eaac..207affa84b618 100644
--- a/sklearn/cluster/k_means_.py
+++ b/sklearn/cluster/k_means_.py
@@ -170,6 +170,9 @@ def k_means(X, n_clusters, init='k-means++', precompute_distances='auto',
             algorithm="auto", return_n_iter=False):
     """K-means clustering algorithm.
 
+    Warning: The data will be converted to C ordering, which might cause a
+    memory copy if the given data is in fortran order.
+
     Read more in the :ref:`User Guide <k_means>`.
 
     Parameters
@@ -299,7 +302,7 @@ def k_means(X, n_clusters, init='k-means++', precompute_distances='auto',
 
     # Validate init array
     if hasattr(init, '__array__'):
-        init = check_array(init, dtype=X.dtype.type, copy=True)
+        init = check_array(init, dtype=X.dtype.type, copy=True, order='C')
         _validate_center_shape(X, n_clusters, init)
 
         if n_init != 1:
@@ -393,7 +396,7 @@ def _kmeans_single_elkan(X, n_clusters, max_iter=300, init='k-means++',
                          precompute_distances=True):
     if sp.issparse(X):
         raise ValueError("algorithm='elkan' not supported for sparse input X")
-    X = check_array(X, order="C")
+    X = check_array(X, order='C')
     random_state = check_random_state(random_state)
     if x_squared_norms is None:
         x_squared_norms = row_norms(X, squared=True)
@@ -716,6 +719,9 @@ class KMeans(BaseEstimator, ClusterMixin, TransformerMixin):
 
     Read more in the :ref:`User Guide <k_means>`.
 
+    Warning: The data will be converted to C ordering, which might cause a
+    memory copy if the given data is in fortran order.
+
     Parameters
     ----------
 
@@ -862,14 +868,16 @@ def __init__(self, n_clusters=8, init='k-means++', n_init=10,
 
     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])
+        X = check_array(X, accept_sparse='csr', dtype=[np.float64, np.float32],
+                        order='C')
         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)
+        X = check_array(X, accept_sparse='csr', dtype=FLOAT_DTYPES,
+                        order='C')
         n_samples, n_features = X.shape
         expected_n_features = self.cluster_centers_.shape[1]
         if not n_features == expected_n_features:
@@ -1225,6 +1233,9 @@ class MiniBatchKMeans(KMeans):
 
     Read more in the :ref:`User Guide <mini_batch_kmeans>`.
 
+    Warning: The data will be converted to C ordering, which might cause a
+    memory copy if the given data is in fortran order.
+
     Parameters
     ----------
 
@@ -1522,7 +1533,7 @@ def partial_fit(self, X, y=None):
 
         """
 
-        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)


Code review comments:
diff_hunk: @@ -393,7 +396,7 @@ def _kmeans_single_elkan(X, n_clusters, max_iter=300, init='k-means++',
                          precompute_distances=True):
     if sp.issparse(X):
         raise ValueError("algorithm='elkan' not supported for sparse input X")
-    X = check_array(X, order="C")
+    X = check_array(X, order='C')
path: sklearn/cluster/k_means_.py
text: @user1:
unnecessary change

diff_hunk: @@ -716,6 +719,9 @@ class KMeans(BaseEstimator, ClusterMixin, TransformerMixin):
 
     Read more in the :ref:`User Guide <k_means>`.
 
+    Warning: The data will be converted to C ordering, which might cause a
path: sklearn/cluster/k_means_.py
text: @user1:
I would avoid to issue a warning here. You can add a line of explanation under X, similarly to the some explanation regarding sparse matrix [CSR](http://scikit-learn.org/dev/modules/generated/sklearn.preprocessing.QuantileTransformer.html#sklearn.preprocessing.QuantileTransformer.fit)

diff_hunk: @@ -1225,6 +1233,9 @@ class MiniBatchKMeans(KMeans):
 
     Read more in the :ref:`User Guide <mini_batch_kmeans>`.
 
+    Warning: The data will be converted to C ordering, which might cause a
path: sklearn/cluster/k_means_.py
text: @user1:
I would avoid to issue a warning here. You can add a line of explanation under X, similarly to the some explanation regarding sparse matrix [CSR](http://scikit-learn.org/dev/modules/generated/sklearn.preprocessing.QuantileTransformer.html#sklearn.preprocessing.QuantileTransformer.fit)

diff_hunk: @@ -170,6 +170,9 @@ def k_means(X, n_clusters, init='k-means++', precompute_distances='auto',
             algorithm="auto", return_n_iter=False):
     """K-means clustering algorithm.
 
+    Warning: The data will be converted to C ordering, which might cause a
path: sklearn/cluster/k_means_.py
text: @user1:
I would avoid to issue a warning here. You can add a line of explanation under X, similarly to the some explanation regarding sparse matrix [CSR](http://scikit-learn.org/dev/modules/generated/sklearn.preprocessing.QuantileTransformer.html#sklearn.preprocessing.QuantileTransformer.fit)

Also a conversion will for sure make a memory copy so you should remove "might" :)