clusterer.cc

   1 /*
   2  *  clueless-kmean is a variant of k-mean which enforces balanced
   3  *  distribution of classes in every cluster
   4  *
   5  *  Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
   6  *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
   7  *
   8  *  This file is part of clueless-kmean.
   9  *
  10  *  clueless-kmean is free software: you can redistribute it and/or
  11  *  modify it under the terms of the GNU General Public License
  12  *  version 3 as published by the Free Software Foundation.
  13  *
  14  *  clueless-kmean is distributed in the hope that it will be useful,
  15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  *  General Public License for more details.
  18  *
  19  *  You should have received a copy of the GNU General Public License
  20  *  along with selector.  If not, see <http://www.gnu.org/licenses/>.
  21  *
  22  */
  23
  24 #include "clusterer.h"
  25 #include <glpk.h>
  26
  27 Clusterer::Clusterer() {
  28   _cluster_means = 0;
  29   _cluster_var = 0;
  30 }
  31
  32 Clusterer::~Clusterer() {
  33   deallocate_array<scalar_t>(_cluster_means);
  34   deallocate_array<scalar_t>(_cluster_var);
  35 }
  36
  37 scalar_t Clusterer::distance_to_centroid(scalar_t *x, int k) {
  38   // We take the variance into account + the normalization term. This
  39   // is between k-mean and EM with a diagonal covariance
  40   scalar_t dist = 0;
  41   for(int d = 0; d < _dim; d++) {
  42     dist += sq(_cluster_means[k][d] - x[d]) / (2 * _cluster_var[k][d]);
  43     dist += 0.5 * log(_cluster_var[k][d]);
  44     ASSERT(!isnan(dist) && !isinf(dist));
  45   }
  46   return dist;
  47 }
  48
  49 void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
  50   int *used = new int[nb_points];
  51   for(int k = 0; k < nb_points; k++) { used[k] = 0; }
  52   for(int k = 0; k < _nb_clusters; k++) {
  53     int l;
  54     do { l = int(drand48() * nb_points); } while(used[l]);
  55     for(int d = 0; d < _dim; d++) {
  56       _cluster_means[k][d] = points[l][d];
  57       _cluster_var[k][d] = 1.0;
  58     }
  59     used[l] = 1;
  60   }
  61   delete[] used;
  62 }
  63
  64 scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **points,
  65                                                  int nb_classes, int *labels,
  66                                                  scalar_t **gamma)  {
  67   int *associated_clusters = new int[nb_points];
  68   scalar_t total_dist = 0;
  69
  70   for(int n = 0; n < nb_points; n++) {
  71     scalar_t lowest_dist = 0;
  72     for(int k = 0; k < _nb_clusters; k++) {
  73       scalar_t dist = distance_to_centroid(points[n], k);
  74       if(k == 0 || dist <= lowest_dist) {
  75         lowest_dist = dist;
  76         associated_clusters[n] = k;
  77       }
  78     }
  79
  80     total_dist += lowest_dist;
  81   }
  82
  83   for(int n = 0; n < nb_points; n++) {
  84     for(int k = 0; k < _nb_clusters; k++) {
  85       gamma[n][k] = 0.0;
  86     }
  87     gamma[n][associated_clusters[n]] = 1.0;
  88   }
  89
  90   delete[] associated_clusters;
  91
  92   return total_dist;
  93 }
  94
  95 scalar_t Clusterer::baseline_lp_cluster_association(int nb_points, scalar_t **points,
  96                                                     int nb_classes, int *labels,
  97                                                     scalar_t **gamma)  {
  98   glp_prob *lp;
  99
 100   int *ia = new int[nb_points * _nb_clusters + 1];
 101   int *ja = new int[nb_points * _nb_clusters + 1];
 102   scalar_t *ar = new scalar_t[nb_points * _nb_clusters + 1];
 103
 104   lp = glp_create_prob();
 105
 106   glp_set_prob_name(lp, "baseline_lp_cluster_association");
 107   glp_set_obj_dir(lp, GLP_MIN);
 108
 109   glp_add_rows(lp, nb_points);
 110
 111   for(int n = 1; n <= nb_points; n++) {
 112     glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
 113   }
 114
 115   glp_add_cols(lp, nb_points * _nb_clusters);
 116   for(int k = 1; k <= _nb_clusters; k++) {
 117     for(int n = 1; n <= nb_points; n++) {
 118       int i = n + nb_points * (k - 1);
 119       glp_set_obj_coef(lp, i, distance_to_centroid(points[n-1], k-1));
 120       glp_set_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
 121     }
 122   }
 123
 124   int l = 1;
 125
 126   for(int n = 1; n <= nb_points; n++) {
 127     for(int k = 1; k <= _nb_clusters; k++) {
 128       ia[l] = n;
 129       ja[l] = n + nb_points * (k - 1);
 130       ar[l] = 1.0;
 131       l++;
 132     }
 133   }
 134
 135   glp_load_matrix(lp, nb_points * _nb_clusters, ia, ja, ar);
 136
 137   glp_simplex(lp, NULL);
 138
 139   scalar_t total_dist = glp_get_obj_val(lp);
 140
 141   for(int k = 1; k <= _nb_clusters; k++) {
 142     for(int n = 1; n <= nb_points; n++) {
 143       int i = n + nb_points * (k - 1);
 144       gamma[n-1][k-1] = glp_get_col_prim(lp, i);
 145     }
 146   }
 147
 148   delete[] ia;
 149   delete[] ja;
 150   delete[] ar;
 151
 152   glp_delete_prob(lp);
 153
 154   return total_dist;
 155 }
 156
 157 scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t **points,
 158                                                          int nb_classes, int *labels,
 159                                                          scalar_t **gamma)  {
 160   glp_prob *lp;
 161
 162   int nb_coeffs = nb_points * _nb_clusters + nb_points * _nb_clusters;
 163
 164   int *ia = new int[nb_coeffs + 1];
 165   int *ja = new int[nb_coeffs + 1];
 166   scalar_t *ar = new scalar_t[nb_coeffs + 1];
 167
 168   scalar_t *nb_samples_per_class = new scalar_t[nb_classes];
 169   for(int c = 0; c < nb_classes; c++) {
 170     nb_samples_per_class[c] = 0.0;
 171   }
 172
 173   for(int n = 0; n < nb_points; n++) {
 174     nb_samples_per_class[labels[n]] += 1.0;
 175   }
 176
 177   lp = glp_create_prob();
 178
 179   glp_set_prob_name(lp, "uninformative_lp_cluster_association");
 180   glp_set_obj_dir(lp, GLP_MIN);
 181
 182   // We have one constraint per point and one per cluster/class
 183
 184   glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
 185
 186   // (A) For each point, the constraint is that the sum of its
 187   // association coefficients is equal to 1.0
 188
 189   for(int n = 1; n <= nb_points; n++) {
 190     int row = n;
 191     glp_set_row_bnds(lp, row, GLP_FX, 1.0, 1.0);
 192   }
 193
 194   // (B) For each pair cluster/class, the sum of the association
 195   // coefficient to this cluster for this class is equal to the number
 196   // of sample of that class, divided by the number of clusters
 197
 198   for(int k = 1; k <= _nb_clusters; k++) {
 199     for(int c = 1; c <= nb_classes; c++) {
 200       int row = nb_points + (k - 1) * nb_classes + c;
 201       scalar_t tau = nb_samples_per_class[c-1] / scalar_t(_nb_clusters);
 202       glp_set_row_bnds(lp, row, GLP_FX, tau, tau);
 203     }
 204   }
 205
 206   // Each one of the constraints above involve a linear combination of
 207   // all the association coefficients
 208
 209   glp_add_cols(lp, nb_points * _nb_clusters);
 210
 211   for(int k = 1; k <= _nb_clusters; k++) {
 212     for(int n = 1; n <= nb_points; n++) {
 213       int col = n + nb_points * (k - 1);
 214
 215       // The LP weight on this association coefficient for the global
 216       // loss is the normalized distance of that sample to the
 217       // centroid of that cluster
 218
 219       glp_set_obj_coef(lp, col, distance_to_centroid(points[n-1], k-1));
 220
 221       // And all the association coefficient is in [0,1]
 222
 223       glp_set_col_bnds(lp, col, GLP_DB, 0.0, 1.0);
 224     }
 225   }
 226
 227   int l = 1;
 228
 229   // We build the matrix of the LP coefficients
 230
 231   // The sums of the association coefficients per points for the
 232   // constraints (A) above.
 233
 234   for(int n = 1; n <= nb_points; n++) {
 235     for(int k = 1; k <= _nb_clusters; k++) {
 236       int row = n;
 237       ia[l] = row;
 238       ja[l] = nb_points * (k - 1) + n;
 239       ar[l] = 1.0;
 240       l++;
 241     }
 242   }
 243
 244   // And the sums of coefficients for each pair class/cluster for
 245   // constraint (B) above.
 246
 247   for(int k = 1; k <= _nb_clusters; k++) {
 248     for(int c = 1; c <= nb_classes; c++) {
 249       int row = nb_points + (k - 1) * nb_classes + c;
 250       for(int n = 1; n <= nb_points; n++) {
 251         if(labels[n-1] == c - 1) {
 252           ia[l] = row;
 253           ja[l] = (k-1)  * nb_points + n;
 254           ar[l] = 1.0;
 255           l++;
 256         }
 257       }
 258     }
 259   }
 260
 261   ASSERT(l == nb_coeffs + 1);
 262
 263   glp_load_matrix(lp, nb_coeffs, ia, ja, ar);
 264
 265   // Now a miracle occurs
 266
 267   glp_simplex(lp, NULL);
 268
 269   // We retrieve the result
 270
 271   scalar_t total_dist = glp_get_obj_val(lp);
 272
 273   for(int k = 1; k <= _nb_clusters; k++) {
 274     for(int n = 1; n <= nb_points; n++) {
 275       int i = n + nb_points * (k - 1);
 276       gamma[n-1][k-1] = glp_get_col_prim(lp, i);
 277     }
 278   }
 279
 280   delete[] nb_samples_per_class;
 281   delete[] ia;
 282   delete[] ja;
 283   delete[] ar;
 284   glp_delete_prob(lp);
 285
 286   return total_dist;
 287 }
 288
 289 void Clusterer::update_clusters(int nb_points, scalar_t **points, scalar_t **gamma)  {
 290   for(int k = 0; k < _nb_clusters; k++) {
 291
 292     for(int d = 0; d < _dim; d++) {
 293       _cluster_means[k][d] = 0.0;
 294       _cluster_var[k][d] = 0.0;
 295     }
 296
 297     scalar_t sum_gamma = 0;
 298     for(int n = 0; n < nb_points; n++) {
 299       sum_gamma += gamma[n][k];
 300       for(int d = 0; d < _dim; d++) {
 301         _cluster_means[k][d] += gamma[n][k] * points[n][d];
 302         _cluster_var[k][d] += gamma[n][k] * sq(points[n][d]);
 303       }
 304     }
 305
 306     ASSERT(sum_gamma >= 1);
 307
 308     for(int d = 0; d < _dim; d++) {
 309       if(sum_gamma >= 2) {
 310         _cluster_var[k][d] =
 311           (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
 312         _cluster_var[k][d] = max(scalar_t(min_cluster_variance), _cluster_var[k][d]);
 313       } else {
 314         _cluster_var[k][d] = 1;
 315       }
 316
 317       _cluster_means[k][d] /= sum_gamma;
 318     }
 319   }
 320 }
 321
 322 void Clusterer::train(int mode,
 323                       int nb_clusters, int dim,
 324                       int nb_points, scalar_t **points,
 325                       int nb_classes, int *labels,
 326                       int *cluster_associations) {
 327   deallocate_array<scalar_t>(_cluster_means);
 328   deallocate_array<scalar_t>(_cluster_var);
 329   _nb_clusters = nb_clusters;
 330   _dim = dim;
 331   _cluster_means = allocate_array<scalar_t>(_nb_clusters, _dim);
 332   _cluster_var = allocate_array<scalar_t>(_nb_clusters, _dim);
 333
 334   scalar_t **gammas = allocate_array<scalar_t>(nb_points, _nb_clusters);
 335
 336   if(nb_clusters > nb_points) abort();
 337
 338   initialize_clusters(nb_points, points);
 339
 340   scalar_t pred_total_distance, total_distance = FLT_MAX;
 341   int nb_rounds = 0;
 342
 343   do {
 344     pred_total_distance = total_distance;
 345
 346     switch(mode) {
 347
 348     case STANDARD_ASSOCIATION:
 349     total_distance =
 350       baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
 351       break;
 352
 353     case STANDARD_LP_ASSOCIATION:
 354     total_distance =
 355       baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
 356       break;
 357
 358     case UNINFORMATIVE_LP_ASSOCIATION:
 359     total_distance =
 360       uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
 361       break;
 362
 363     default:
 364       cerr << "Unknown sample-cluster association mode." << endl;
 365       abort();
 366     }
 367
 368     cout << "TRAIN " << nb_rounds << " " << total_distance << endl;
 369     update_clusters(nb_points, points, gammas);
 370     nb_rounds++;
 371   } while(total_distance < min_iteration_improvement * pred_total_distance &&
 372           nb_rounds < max_nb_iterations);
 373
 374   if(cluster_associations) {
 375     for(int n = 0; n < nb_points; n++) {
 376       for(int k = 0; k < _nb_clusters; k++) {
 377         if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
 378           cluster_associations[n] = k;
 379         }
 380       }
 381     }
 382   }
 383
 384   deallocate_array<scalar_t>(gammas);
 385 }
 386
 387 int Clusterer::cluster(scalar_t *point) {
 388   scalar_t lowest_dist = 0;
 389   int associated_cluster = -1;
 390
 391   for(int k = 0; k < _nb_clusters; k++) {
 392     scalar_t dist = 0;
 393
 394     for(int d = 0; d < _dim; d++) {
 395       dist += sq(_cluster_means[k][d] - point[d]) / (2 * _cluster_var[k][d]);
 396       dist += 0.5 * log(_cluster_var[k][d]);
 397       ASSERT(!isnan(dist) && !isinf(dist));
 398     }
 399
 400     if(k == 0 || dist <= lowest_dist) {
 401       lowest_dist = dist;
 402       associated_cluster = k;
 403     }
 404   }
 405
 406   ASSERT(associated_cluster >= 0);
 407
 408   return associated_cluster;
 409 }