KMeans.h
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief The k-means clustering algorithm.
6  *
7  *
8  *
9  * \author T. Glasmachers
10  * \date 2011
11  *
12  *
13  * \par Copyright 1995-2015 Shark Development Team
14  *
15  * <BR><HR>
16  * This file is part of Shark.
17  * <http://image.diku.dk/shark/>
18  *
19  * Shark is free software: you can redistribute it and/or modify
20  * it under the terms of the GNU Lesser General Public License as published
21  * by the Free Software Foundation, either version 3 of the License, or
22  * (at your option) any later version.
23  *
24  * Shark is distributed in the hope that it will be useful,
25  * but WITHOUT ANY WARRANTY; without even the implied warranty of
26  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27  * GNU Lesser General Public License for more details.
28  *
29  * You should have received a copy of the GNU Lesser General Public License
30  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
31  *
32  */
33 //===========================================================================
34 
35 
36 #ifndef SHARK_ALGORITHMS_KMEANS_H
37 #define SHARK_ALGORITHMS_KMEANS_H
38 
39 #include <shark/Core/DLLSupport.h>
40 #include <shark/Data/Dataset.h>
42 #include <shark/Models/RBFLayer.h>
45 
46 
47 namespace shark{
48 
49 
50 ///
51 /// \brief The k-means clustering algorithm.
52 ///
53 /// \par
54 /// The k-means algorithm takes vector-valued data
55 /// \f$ \{x_1, \dots, x_\ell\} \subset \mathbb R^d \f$
56 /// and splits it into k clusters, based on centroids
57 /// \f$ \{c_1, \dots, c_k\} \f$.
58 /// The result is stored in a Centroids object that can be used to
59 /// construct clustering models.
60 ///
61 /// \par
62 /// This implementation starts the search with the given centroids,
63 /// in case the provided centroids object (third parameter) contains
64 /// a set of k centroids. Otherwise the search starts from the first
65 /// k data points.
66 ///
67 /// \par
68 /// Note that the data set needs to include at least k data points
69 /// for k-means to work. This is because the current implementation
70 /// does not allow for empty clusters.
71 ///
72 /// \param data vector-valued data to be clustered
73 /// \param k number of clusters
74 /// \param centroids centroids input/output
75 /// \param maxIterations maximum number of k-means iterations; 0: unlimited
76 /// \return number of k-means iterations
77 ///
78 SHARK_EXPORT_SYMBOL std::size_t kMeans(Data<RealVector> const& data, std::size_t k, Centroids& centroids, std::size_t maxIterations = 0);
79 
80 ///
81 /// \brief The k-means clustering algorithm for initializing an RBF Layer
82 ///
83 /// \par
84 /// The k-means algorithm takes vector-valued data
85 /// \f$ \{x_1, \dots, x_\ell\} \subset \mathbb R^d \f$
86 /// and splits it into k clusters, based on centroids
87 /// \f$ \{c_1, \dots, c_k\} \f$.
88 /// The result is stored in a RBFLayer object that can be used to
89 /// construct clustering models.
90 ///
91 /// \par
92 /// This is just an alternative frontend to the version using Centroids. it creates a centroid object,
93 /// with as many clusters as are outputs in the RBFLayer and copis the result into the model.
94 ///
95 /// \par
96 /// Note that the data set needs to include at least k data points
97 /// for k-means to work. This is because the current implementation
98 /// does not allow for empty clusters.
99 ///
100 /// \param data vector-valued data to be clustered
101 /// \param model RBFLayer input/output
102 /// \param maxIterations maximum number of k-means iterations; 0: unlimited
103 /// \return number of k-means iterations
104 ///
105 SHARK_EXPORT_SYMBOL std::size_t kMeans(Data<RealVector> const& data, RBFLayer& model, std::size_t maxIterations = 0);
106 
107 template<class InputType>
108 KernelExpansion<InputType> kMeans(Data<InputType> const& dataset, std::size_t k, AbstractKernelFunction<InputType>& kernel, std::size_t maxIterations = 0){
109  if(!maxIterations)
110  maxIterations = std::numeric_limits<std::size_t>::max();
111 
112  std::size_t ell = dataset.numberOfElements();
113  RealMatrix kernelMatrix = calculateRegularizedKernelMatrix(kernel,dataset,0);
114  UIntVector clusterMembership(ell);
115  UIntVector clusterSizes(k,0);
116  RealVector ckck(k,0);
117 
118  //init cluster assignments
119  for(unsigned int i = 0; i != ell; ++i){
120  clusterMembership(i) = i % k;
121  }
122  DiscreteUniform<Rng::rng_type> uni(Rng::globalRng,0,k-1);
123  std::random_shuffle(clusterMembership.begin(),clusterMembership.end(),uni);
124  for(std::size_t i = 0; i != ell; ++i){
125  ++clusterSizes(clusterMembership(i));
126  }
127 
128  // k-means loop
129  std::size_t iter = 0;
130  bool equal = false;
131  for(; iter != maxIterations && !equal; ++iter) {
132  //the clustering model results in centers c_k= 1/n_k sum_i k(x_i,.) for all x_i points of cluster k
133  //we need to compute the squared distances between all centers and points that is
134  //d^2(c_k,x_i) = <c_k,c_k> -2 < c_k,x_i> + <x_i,x_i> for the i-th point.
135  //thus we precompute <c_k,c_k>= sum_ij k(x_i,x_j)/(n_k)^2 for all x_i,x_j points of cluster k
136  ckck.clear();
137  for(std::size_t i = 0; i != ell; ++i){
138  std::size_t c1 = clusterMembership(i);
139  for(std::size_t j = 0; j != ell; ++j){
140  std::size_t c2 = clusterMembership(j);
141  if(c1 != c2) continue;
142  ckck(c1) += kernelMatrix(i,j);
143  }
144  }
145  noalias(ckck) = safe_div(ckck,sqr(clusterSizes),0);
146 
147  UIntVector newClusterMembership(kernelMatrix.size1());
148  RealVector currentDistances(k);
149  for(std::size_t i = 0; i != ell; ++i){
150  //compute squared distances between the i-th point and the centers
151  //we skip <x_i,x_i> as it is always the same for all elements and we don't need it for comparison
152  noalias(currentDistances) = ckck;
153  for(std::size_t j = 0; j != ell; ++j){
154  std::size_t c = clusterMembership(j);
155  currentDistances(c) -= 2* kernelMatrix(i,j)/clusterSizes(c);
156  }
157  //choose the index with the smallest distance as new cluster
158  newClusterMembership(i) = (unsigned int) arg_min(currentDistances);
159  }
160  equal = boost::equal(
161  newClusterMembership,
162  clusterMembership
163  );
164  noalias(clusterMembership) = newClusterMembership;
165  //compute new sizes of clusters
166  clusterSizes.clear();
167  for(std::size_t i = 0; i != ell; ++i){
168  ++clusterSizes(clusterMembership(i));
169  }
170 
171  //if a cluster has size , assign a random point to it
172  for(unsigned int i = 0; i != k; ++i){
173  if(clusterSizes(i) == 0){
174  std::size_t elem = uni(ell-1);
175  --clusterSizes(clusterMembership(elem));
176  clusterMembership(elem)=i;
177  clusterSizes(i) = 1;
178  }
179  }
180  }
181 
182  //copy result in the expansion
183  KernelExpansion<InputType> expansion;
184  expansion.setStructure(&kernel,dataset,true,k);
185  expansion.offset() = -ckck;
186  expansion.alpha().clear();
187  for(std::size_t i = 0; i != ell; ++i){
188  std::size_t c = clusterMembership(i);
189  expansion.alpha()(i,c) = 2.0 / clusterSizes(c);
190  }
191 
192  return expansion;
193 }
194 }
195 #endif