Add random kmeans init class.

Author: trivialfis

src/gpu/kmeans/kmeans_init.cu

@@ -5,7 +5,6 @@
 
 #include <thrust/device_vector.h>
 
-#include <cub/device/device_select.cuh>
 #include <cub/device/device_histogram.cuh>
 
 #include <random>
@@ -250,10 +249,28 @@ KmMatrix<T> GreedyRecluster<T>::recluster(KmMatrix<T>& _centroids, size_t _k) {
 
 }  // namespace detail
 
+
+/* ============ KmeansRandomInit Class member functions ============ */
+
+template <typename T>
+KmMatrix<T> KmeansRandomInit<T>::operator()(KmMatrix<T>& _data, size_t _k) {
+
+  KmMatrix<T> dist = generator_impl_->generate(_k);
+  MulOp<T>().mul(dist, dist, (T)_data.rows() - 1);
+
+  dist = dist.reshape(dist.cols(), dist.rows());
+
+  KmMatrix<T> centroids = _data.rows(dist);
+
+  return centroids;
+}
 
 
 /* ============== KmeansLlInit Class member functions ============== */
 
+// Although the paper suggested calculating the probability independently,
+// but due to zero distance between a point and itself, the already selected
+// points will have very low probability.
 template <typename T, template <class> class ReclusterPolicy >
 KmMatrix<T> KmeansLlInit<T, ReclusterPolicy>::probability(
     KmMatrix<T>& _data, KmMatrix<T>& _centroids) {
@@ -360,17 +377,18 @@ KmeansLlInit<T, ReclusterPolicy>::operator()(KmMatrix<T>& _data, size_t _k) {
 
   T cost = SumOp<T>().sum(prob);
 
-  size_t max_iter = std::max(T(MAX_ITER), std::log(cost));
+  size_t max_iter = std::max((size_t)(MAX_ITER),
+                             (size_t)std::ceil(std::log(cost)));
   for (size_t i = 0; i < max_iter; ++i) {
     prob = probability(_data, centroids);
     KmMatrix<T> new_centroids = sample_centroids(_data, prob);
     centroids = stack(centroids, new_centroids, KmMatrixDim::ROW);
   }
 
   if (centroids.rows() < _k) {
-    // FIXME: When n_centroids < k
-    // Get random selection in?
-    M_ERROR("Not implemented.");
+    KmMatrix<T> new_centroids = KmeansRandomInit<T>(generator_)(_data,
+        _k - centroids.rows());
+    centroids = stack(centroids, new_centroids, KmMatrixDim::ROW);
   }
 
   centroids = ReclusterPolicy<T>::recluster(centroids, k_);
@@ -385,9 +403,12 @@ KmeansLlInit<T, ReclusterPolicy>::operator()(KmMatrix<T>& _data, size_t _k) {
       KmMatrix<T>& data, KmMatrix<T>& centroids);                       \
   template KmMatrix<T> KmeansLlInit<T>::sample_centroids(               \
       KmMatrix<T>& data, KmMatrix<T>& centroids);                       \
+  template KmMatrix<T> KmeansRandomInit<T>::operator()(                 \
+      KmMatrix<T>& _data, size_t _k);
 
 INSTANTIATE(float)
 INSTANTIATE(double)
+INSTANTIATE(int)
 
 #undef INSTANTIATE
 
@@ -403,15 +424,14 @@ namespace detail {
       KmMatrix<T>& _centroids_dot,                              \
       KmMatrix<T>& _distance_pairs);                            \
   template KmMatrix<T> PairWiseDistanceOp<T>::operator()(       \
-      KmMatrix<T>& _data,                                       \
-      KmMatrix<T>& _centroids);                                 \
+      KmMatrix<T>& _data, KmMatrix<T>& _centroids);             \
 
 INSTANTIATE(float)
 INSTANTIATE(double)
+INSTANTIATE(int)
 
 #undef INSTANTIATE
 }
-// FIXME: int is not supported due to random kernel
 
 }  // namespace Kmeans
 }  // namespace H2O4GPU

src/gpu/kmeans/kmeans_init.cuh

@@ -73,6 +73,38 @@ class KmeansInitBase {
   virtual KmMatrix<T> operator()(KmMatrix<T>& data, size_t k) = 0;
 };
 
+/*
+ * Random initialization.
+ * @tparam Numeric data type.
+ */
+template <typename T>
+class KmeansRandomInit : public KmeansInitBase<T> {
+ private:
+  int seed_;
+  std::unique_ptr<GeneratorBase<T>> generator_impl_;
+
+ public:
+  /*
+   * @param seed Random seed for generating centroids.
+   */
+  KmeansRandomInit(size_t _seed) :
+      seed_(_seed), generator_impl_ (new UniformGenerator<T>) {}
+
+  /*
+   * @param gen Unique pointer to Random generator for generating centroids.
+   */
+  KmeansRandomInit(std::unique_ptr<GeneratorBase<T>>& _gen) :
+      generator_impl_(std::move(_gen)) {}
+
+  virtual ~KmeansRandomInit() override {}
+
+  /*
+   * @param data Data points stored in row major matrix.
+   * @param k Number of centroids.
+   */
+  virtual KmMatrix<T> operator()(KmMatrix<T>& data, size_t k) override;
+};
+
 /*
  * Each instance of KmeansLlInit corresponds to one dataset, if a new data set
  * is used, users need to create a new instance.
@@ -82,7 +114,7 @@ class KmeansInitBase {
  * Scalable K-Means++
  * </a>
  *
- * @tparam Data type, supported types are float and double.
+ * @tparam Numeric data type.
  */
 template <
   typename T,
@@ -107,8 +139,8 @@ struct KmeansLlInit : public KmeansInitBase<T> {
 
   KmMatrix<T> probability(KmMatrix<T>& data, KmMatrix<T>& centroids);
  public:
-  // sample_centroids/recluster should not be part of the interface, but
-  // following error is generated when put in private section:
+  // sample_centroids should not be part of the interface, but following error
+  // is generated when put in private section:
   // The enclosing parent function ("sample_centroids") for an extended
   // __device__ lambda cannot have private or protected access within its class
   KmMatrix<T> sample_centroids(KmMatrix<T>& data, KmMatrix<T>& centroids);

tests/cpp/gpu/kmeans/test_kmeans_init.cu

@@ -36,6 +36,26 @@ struct GeneratorMock : GeneratorBase<T> {
   }
 };
 
+TEST(KmeansRandom, Init) {
+  thrust::host_vector<float> h_data (20);
+  for (size_t i = 0; i < 20; ++i) {
+    h_data[i] = i * 2;
+  }
+  KmMatrix<float> data (h_data, 4, 5);
+  std::unique_ptr<GeneratorBase<float>> gen (new GeneratorMock<float>());
+  KmeansRandomInit<float> init (gen);
+
+  auto res = init(data, 2);
+
+  std::vector<float> h_sol =
+      {
+        30, 32, 34, 36, 38,
+        0,  2,  4,  6,  8
+      };
+  KmMatrix<float> sol (h_sol, 2, 5);
+  ASSERT_TRUE(sol == res);
+}
+
 // r --gtest_filter=KmeansLL.PairWiseDistance
 TEST(KmeansLL, PairWiseDistance) {