clusterer.cc

   1 /*
   2  *  clueless-kmeans is a variant of k-means 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-kmeans.
   9  *
  10  *  clueless-kmeans 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-kmeans 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   scalar_t dist = 0;
  39   for(int d = 0; d < _dim; d++) {
  40     dist += sq(_cluster_means[k][d] - x[d]) / (2 * _cluster_var[k][d]);
  41     dist += 0.5 * log(_cluster_var[k][d]);
  42     ASSERT(!isnan(dist) && !isinf(dist));
  43   }
  44   return dist;
  45 }
  46
  47 void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
  48   int *used = new int[nb_points];
  49   for(int k = 0; k < nb_points; k++) { used[k] = 0; }
  50   for(int k = 0; k < _nb_clusters; k++) {
  51     int l;
  52     do { l = int(drand48() * nb_points); } while(used[l]);
  53     for(int d = 0; d < _dim; d++) {
  54       _cluster_means[k][d] = points[l][d];
  55       _cluster_var[k][d] = 1.0;
  56     }
  57     used[l] = 1;
  58   }
  59   delete[] used;
  60 }
  61
  62 scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **points,
  63                                                  int nb_classes, int *labels,
  64                                                  scalar_t **gamma)  {
  65   int *associated_clusters = new int[nb_points];
  66   scalar_t total_dist = 0;
  67
  68   for(int n = 0; n < nb_points; n++) {
  69     scalar_t lowest_dist = 0;
  70     for(int k = 0; k < _nb_clusters; k++) {
  71       scalar_t dist = distance_to_centroid(points[n], k);
  72       if(k == 0 || dist <= lowest_dist) {
  73         lowest_dist = dist;
  74         associated_clusters[n] = k;
  75       }
  76     }
  77
  78     total_dist += lowest_dist;
  79   }
  80
  81   for(int n = 0; n < nb_points; n++) {
  82     for(int k = 0; k < _nb_clusters; k++) {
  83       gamma[n][k] = 0.0;
  84     }
  85     gamma[n][associated_clusters[n]] = 1.0;
  86   }
  87
  88   delete[] associated_clusters;
  89
  90   return total_dist;
  91 }
  92
  93 scalar_t Clusterer::baseline_lp_cluster_association(int nb_points, scalar_t **points,
  94                                                     int nb_classes, int *labels,
  95                                                     scalar_t **gamma)  {
  96   glp_prob *lp;
  97
  98   int *coeff_row = new int[nb_points * _nb_clusters + 1];
  99   int *coeff_col = new int[nb_points * _nb_clusters + 1];
 100   scalar_t *coeff_wgt = new scalar_t[nb_points * _nb_clusters + 1];
 101
 102   lp = glp_create_prob();
 103
 104   glp_set_prob_name(lp, "baseline_lp_cluster_association");
 105   glp_set_obj_dir(lp, GLP_MIN);
 106
 107   glp_add_rows(lp, nb_points);
 108
 109   for(int n = 1; n <= nb_points; n++) {
 110     glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
 111   }
 112
 113   glp_add_cols(lp, nb_points * _nb_clusters);
 114   for(int k = 1; k <= _nb_clusters; k++) {
 115     for(int n = 1; n <= nb_points; n++) {
 116       int i = n + nb_points * (k - 1);
 117       glp_set_obj_coef(lp, i, distance_to_centroid(points[n-1], k-1));
 118       glp_set_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
 119     }
 120   }
 121
 122   int n_coeff = 1;
 123
 124   for(int n = 1; n <= nb_points; n++) {
 125     for(int k = 1; k <= _nb_clusters; k++) {
 126       coeff_row[n_coeff] = n;
 127       coeff_col[n_coeff] = n + nb_points * (k - 1);
 128       coeff_wgt[n_coeff] = 1.0;
 129       n_coeff++;
 130     }
 131   }
 132
 133   glp_load_matrix(lp, nb_points * _nb_clusters, coeff_row, coeff_col, coeff_wgt);
 134
 135   glp_simplex(lp, NULL);
 136
 137   scalar_t total_dist = glp_get_obj_val(lp);
 138
 139   for(int k = 1; k <= _nb_clusters; k++) {
 140     for(int n = 1; n <= nb_points; n++) {
 141       int i = n + nb_points * (k - 1);
 142       gamma[n-1][k-1] = glp_get_col_prim(lp, i);
 143     }
 144   }
 145
 146   delete[] coeff_row;
 147   delete[] coeff_col;
 148   delete[] coeff_wgt;
 149
 150   glp_delete_prob(lp);
 151
 152   return total_dist;
 153 }
 154
 155 scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t **points,
 156                                                          int nb_classes, int *labels,
 157                                                          scalar_t **gamma,
 158                                                          int absolute_proportion)  {
 159   // N points
 160   // K clusters
 161   // dist(n,k) distance of samples n to cluster k
 162   //
 163   // We want to optimize the
 164   //
 165   // \gamma(n,k) is the association of point n to cluster k
 166   //
 167   // to minimize
 168   //
 169   // \sum_{n,k} \gamma(n,k) dist(n,k)
 170   //
 171   // under
 172   //
 173   // (A) \forall n, k, \gamma(n, k) >= 0
 174   // (B) \forall n, \sum_k \gamma(n, k) = 1
 175   // (C) \forall k, \sum_n \gamma(n, k) = N/K
 176
 177   glp_prob *lp;
 178
 179   // The coefficients for the constraints are passed to the glpk
 180   // functions with a sparse representation.
 181
 182   // ** GLPK USES INDEXES STARTING AT 1, NOT 0. **
 183
 184   int nb_coeffs;
 185
 186   if(absolute_proportion) {
 187     nb_coeffs = nb_points * _nb_clusters + nb_points * _nb_clusters;
 188   } else {
 189     nb_coeffs = nb_points * _nb_clusters + nb_points * nb_classes * _nb_clusters;
 190   }
 191
 192   int *coeff_row = new int[nb_coeffs + 1];
 193   int *coeff_col = new int[nb_coeffs + 1];
 194   scalar_t *coeff_wgt = new scalar_t[nb_coeffs + 1];
 195
 196   int n_coeff = 1;
 197
 198   scalar_t *nb_samples_per_class = new scalar_t[nb_classes];
 199   for(int c = 0; c < nb_classes; c++) {
 200     nb_samples_per_class[c] = 0.0;
 201   }
 202
 203   for(int n = 0; n < nb_points; n++) {
 204     nb_samples_per_class[labels[n]] += 1.0;
 205   }
 206
 207   lp = glp_create_prob();
 208
 209   glp_set_prob_name(lp, "uninformative_lp_cluster_association");
 210   glp_set_obj_dir(lp, GLP_MIN);
 211
 212   // We have one column per coefficient gamma
 213
 214   glp_add_cols(lp, nb_points * _nb_clusters);
 215
 216   // The column for gamma[n][k] point 1<=n<=nb_points and cluster
 217   // 1<=k<=_nb_clusters is nb_points * (k - 1) + n;
 218
 219   // The constraints (A) will be expressed by putting directly bounds
 220   // on the variables (i.e. one per column). So we need one row per
 221   // (B) constraint, and one per (C) constraint.
 222
 223   glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
 224
 225   // First, we set the weights for the objective function, and the (A)
 226   // constraints on the individual gammas
 227
 228   for(int k = 1; k <= _nb_clusters; k++) {
 229     for(int n = 1; n <= nb_points; n++) {
 230       int col = n + nb_points * (k - 1);
 231
 232       // The LP weight on the gammas for the global loss is the
 233       // normalized distance of that sample to the centroid of that
 234       // cluster
 235
 236       glp_set_obj_coef(lp, col, distance_to_centroid(points[n-1], k-1));
 237
 238       // The (A) constraints: Each column correspond to one of the
 239       // gamma, and it has to be in [0,1]
 240
 241       glp_set_col_bnds(lp, col, GLP_DB, 0.0, 1.0);
 242     }
 243   }
 244
 245   // The (B) constraints: for each point, the sum of its gamma is
 246   // equal to 1.0
 247
 248   for(int n = 1; n <= nb_points; n++) {
 249     int row = n;
 250     glp_set_row_bnds(lp, row, GLP_FX, 1.0, 1.0);
 251     for(int k = 1; k <= _nb_clusters; k++) {
 252       coeff_row[n_coeff] = row;
 253       coeff_col[n_coeff] = nb_points * (k - 1) + n;
 254       coeff_wgt[n_coeff] = 1.0;
 255       n_coeff++;
 256     }
 257   }
 258
 259   // The (C) constraints: For each pair cluster/class, the sum of the
 260   // gammas for this cluster and this class is equal to the number of
 261   // sample of that class, divided by the number of clusters
 262
 263   if(absolute_proportion) {
 264     for(int k = 1; k <= _nb_clusters; k++) {
 265       for(int c = 1; c <= nb_classes; c++) {
 266         int row = nb_points + (k - 1) * nb_classes + c;
 267         scalar_t tau = nb_samples_per_class[c-1] / scalar_t(_nb_clusters);
 268         glp_set_row_bnds(lp, row, GLP_FX, tau, tau);
 269         for(int n = 1; n <= nb_points; n++) {
 270           if(labels[n-1] == c - 1) {
 271             coeff_row[n_coeff] = row;
 272             coeff_col[n_coeff] = (k-1)  * nb_points + n;
 273             coeff_wgt[n_coeff] = 1.0;
 274             n_coeff++;
 275           }
 276         }
 277       }
 278     }
 279   } else {
 280     for(int k = 1; k <= _nb_clusters; k++) {
 281       for(int c = 1; c <= nb_classes; c++) {
 282         int row = nb_points + (k - 1) * nb_classes + c;
 283         glp_set_row_bnds(lp, row, GLP_FX, 0.0, 0.0);
 284         for(int n = 1; n <= nb_points; n++) {
 285           coeff_row[n_coeff] = row;
 286           coeff_col[n_coeff] = (k-1)  * nb_points + n;
 287           coeff_wgt[n_coeff] =
 288             (labels[n-1] == c - 1 ? 1.0 : 0.0)
 289             - scalar_t(nb_samples_per_class[c-1]) / scalar_t(nb_points);
 290           n_coeff++;
 291         }
 292       }
 293     }
 294   }
 295
 296   ASSERT(n_coeff == nb_coeffs + 1);
 297
 298   glp_load_matrix(lp, nb_coeffs, coeff_row, coeff_col, coeff_wgt);
 299
 300   // Now a miracle occurs
 301
 302   glp_simplex(lp, NULL);
 303
 304   // We retrieve the result
 305
 306   scalar_t total_dist = glp_get_obj_val(lp);
 307
 308   for(int k = 1; k <= _nb_clusters; k++) {
 309     for(int n = 1; n <= nb_points; n++) {
 310       int i = n + nb_points * (k - 1);
 311       gamma[n-1][k-1] = glp_get_col_prim(lp, i);
 312     }
 313   }
 314
 315   delete[] nb_samples_per_class;
 316   delete[] coeff_row;
 317   delete[] coeff_col;
 318   delete[] coeff_wgt;
 319   glp_delete_prob(lp);
 320
 321   return total_dist;
 322 }
 323
 324 void Clusterer::update_clusters(int nb_points, scalar_t **points, scalar_t **gamma)  {
 325   for(int k = 0; k < _nb_clusters; k++) {
 326
 327     for(int d = 0; d < _dim; d++) {
 328       _cluster_means[k][d] = 0.0;
 329       _cluster_var[k][d] = 0.0;
 330     }
 331
 332     scalar_t sum_gamma = 0;
 333     for(int n = 0; n < nb_points; n++) {
 334       sum_gamma += gamma[n][k];
 335       for(int d = 0; d < _dim; d++) {
 336         _cluster_means[k][d] += gamma[n][k] * points[n][d];
 337         _cluster_var[k][d] += gamma[n][k] * sq(points[n][d]);
 338       }
 339     }
 340
 341     ASSERT(sum_gamma >= 1);
 342
 343     for(int d = 0; d < _dim; d++) {
 344       if(sum_gamma >= 2) {
 345         _cluster_var[k][d] =
 346           (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
 347         _cluster_var[k][d] = max(scalar_t(min_cluster_variance), _cluster_var[k][d]);
 348       } else {
 349         _cluster_var[k][d] = 1;
 350       }
 351
 352       _cluster_means[k][d] /= sum_gamma;
 353     }
 354   }
 355 }
 356
 357 void Clusterer::train(int mode,
 358                       int nb_clusters, int dim,
 359                       int nb_points, scalar_t **points,
 360                       int nb_classes, int *labels,
 361                       int *cluster_associations) {
 362   deallocate_array<scalar_t>(_cluster_means);
 363   deallocate_array<scalar_t>(_cluster_var);
 364   _nb_clusters = nb_clusters;
 365   _dim = dim;
 366   _cluster_means = allocate_array<scalar_t>(_nb_clusters, _dim);
 367   _cluster_var = allocate_array<scalar_t>(_nb_clusters, _dim);
 368
 369   scalar_t **gammas = allocate_array<scalar_t>(nb_points, _nb_clusters);
 370
 371   if(nb_clusters > nb_points) abort();
 372
 373   initialize_clusters(nb_points, points);
 374
 375   scalar_t pred_total_distance, total_distance = FLT_MAX;
 376   int nb_rounds = 0;
 377
 378   do {
 379     pred_total_distance = total_distance;
 380
 381     switch(mode) {
 382
 383     case STANDARD_ASSOCIATION:
 384       total_distance =
 385         baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
 386       break;
 387
 388     case STANDARD_LP_ASSOCIATION:
 389       total_distance =
 390         baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
 391       break;
 392
 393     case UNINFORMATIVE_LP_ASSOCIATION:
 394       total_distance =
 395         uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas, 0);
 396       break;
 397
 398     case UNINFORMATIVE_LP_ASSOCIATION_ABSOLUTE:
 399       total_distance =
 400         uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas, 1);
 401       break;
 402
 403     default:
 404       cerr << "Unknown sample-cluster association mode." << endl;
 405       abort();
 406     }
 407
 408     cout << "TRAIN " << nb_rounds << " " << total_distance << endl;
 409     update_clusters(nb_points, points, gammas);
 410     nb_rounds++;
 411   } while(total_distance < min_iteration_improvement * pred_total_distance &&
 412           nb_rounds < max_nb_iterations);
 413
 414   if(cluster_associations) {
 415     for(int n = 0; n < nb_points; n++) {
 416       for(int k = 0; k < _nb_clusters; k++) {
 417         if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
 418           cluster_associations[n] = k;
 419         }
 420       }
 421     }
 422   }
 423
 424   deallocate_array<scalar_t>(gammas);
 425 }
 426
 427 int Clusterer::cluster(scalar_t *point) {
 428   scalar_t lowest_dist = 0;
 429   int associated_cluster = -1;
 430
 431   for(int k = 0; k < _nb_clusters; k++) {
 432     scalar_t dist = 0;
 433
 434     for(int d = 0; d < _dim; d++) {
 435       dist += sq(_cluster_means[k][d] - point[d]) / (2 * _cluster_var[k][d]);
 436       dist += 0.5 * log(_cluster_var[k][d]);
 437       ASSERT(!isnan(dist) && !isinf(dist));
 438     }
 439
 440     if(k == 0 || dist <= lowest_dist) {
 441       lowest_dist = dist;
 442       associated_cluster = k;
 443     }
 444   }
 445
 446   ASSERT(associated_cluster >= 0);
 447
 448   return associated_cluster;
 449 }