/*
- * clueless-kmean is a variant of k-mean which enforces balanced
+ * clueless-kmeans is a variant of k-means which enforces balanced
* distribution of classes in every cluster
*
* Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
* Written by Francois Fleuret <francois.fleuret@idiap.ch>
*
- * This file is part of clueless-kmean.
+ * This file is part of clueless-kmeans.
*
- * clueless-kmean is free software: you can redistribute it and/or
+ * clueless-kmeans is free software: you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* version 3 as published by the Free Software Foundation.
*
- * clueless-kmean is distributed in the hope that it will be useful,
+ * clueless-kmeans is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
deallocate_array<scalar_t>(_cluster_var);
}
+scalar_t Clusterer::distance_to_centroid(scalar_t *x, int k) {
+ scalar_t dist = 0;
+ for(int d = 0; d < _dim; d++) {
+ dist += sq(_cluster_means[k][d] - x[d]) / (2 * _cluster_var[k][d]);
+ dist += 0.5 * log(_cluster_var[k][d]);
+ ASSERT(!isnan(dist) && !isinf(dist));
+ }
+ return dist;
+}
+
+void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
+ int *used = new int[nb_points];
+ for(int k = 0; k < nb_points; k++) { used[k] = 0; }
+ for(int k = 0; k < _nb_clusters; k++) {
+ int l;
+ do { l = int(drand48() * nb_points); } while(used[l]);
+ for(int d = 0; d < _dim; d++) {
+ _cluster_means[k][d] = points[l][d];
+ _cluster_var[k][d] = 1.0;
+ }
+ used[l] = 1;
+ }
+ delete[] used;
+}
+
scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **points,
int nb_classes, int *labels,
scalar_t **gamma) {
for(int n = 0; n < nb_points; n++) {
scalar_t lowest_dist = 0;
for(int k = 0; k < _nb_clusters; k++) {
- scalar_t dist = 0;
-
- for(int d = 0; d < _dim; d++) {
- dist += sq(_cluster_means[k][d] - points[n][d]) / (2 * _cluster_var[k][d]);
- dist += 0.5 * log(_cluster_var[k][d]);
- ASSERT(!isnan(dist) && !isinf(dist));
- }
-
+ scalar_t dist = distance_to_centroid(points[n], k);
if(k == 0 || dist <= lowest_dist) {
lowest_dist = dist;
associated_clusters[n] = k;
scalar_t **gamma) {
glp_prob *lp;
- int *ia = new int[nb_points * _nb_clusters + 1];
- int *ja = new int[nb_points * _nb_clusters + 1];
- scalar_t *ar = new scalar_t[nb_points * _nb_clusters + 1];
+ int *coeff_row = new int[nb_points * _nb_clusters + 1];
+ int *coeff_col = new int[nb_points * _nb_clusters + 1];
+ scalar_t *coeff_wgt = new scalar_t[nb_points * _nb_clusters + 1];
lp = glp_create_prob();
for(int k = 1; k <= _nb_clusters; k++) {
for(int n = 1; n <= nb_points; n++) {
int i = n + nb_points * (k - 1);
-
- scalar_t dist = 0;
-
- for(int d = 0; d < _dim; d++) {
- dist += sq(_cluster_means[k-1][d] - points[n-1][d]) / (2 * _cluster_var[k-1][d]);
- dist += 0.5 * log(_cluster_var[k-1][d]);
- }
-
- glp_set_obj_coef(lp, i, dist);
+ glp_set_obj_coef(lp, i, distance_to_centroid(points[n-1], k-1));
glp_set_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
}
}
- int l = 1;
+ int n_coeff = 1;
for(int n = 1; n <= nb_points; n++) {
for(int k = 1; k <= _nb_clusters; k++) {
- ia[l] = n;
- ja[l] = n + nb_points * (k - 1);
- ar[l] = 1.0;
- l++;
+ coeff_row[n_coeff] = n;
+ coeff_col[n_coeff] = n + nb_points * (k - 1);
+ coeff_wgt[n_coeff] = 1.0;
+ n_coeff++;
}
}
- glp_load_matrix(lp, nb_points * _nb_clusters, ia, ja, ar);
+ glp_load_matrix(lp, nb_points * _nb_clusters, coeff_row, coeff_col, coeff_wgt);
glp_simplex(lp, NULL);
}
}
- delete[] ia;
- delete[] ja;
- delete[] ar;
+ delete[] coeff_row;
+ delete[] coeff_col;
+ delete[] coeff_wgt;
glp_delete_prob(lp);
scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t **points,
int nb_classes, int *labels,
- scalar_t **gamma) {
+ scalar_t **gamma,
+ int absolute_proportion) {
+ // N points
+ // K clusters
+ // dist(n,k) distance of samples n to cluster k
+ //
+ // We want to optimize the
+ //
+ // \gamma(n,k) is the association of point n to cluster k
+ //
+ // to minimize
+ //
+ // \sum_{n,k} \gamma(n,k) dist(n,k)
+ //
+ // under
+ //
+ // (A) \forall n, k, \gamma(n, k) >= 0
+ // (B) \forall n, \sum_k \gamma(n, k) = 1
+ // (C) \forall k, \sum_n \gamma(n, k) = N/K
+
glp_prob *lp;
- int nb_coeffs = nb_points * _nb_clusters + nb_points * _nb_clusters;
+ // The coefficients for the constraints are passed to the glpk
+ // functions with a sparse representation.
+
+ // ** GLPK USES INDEXES STARTING AT 1, NOT 0. **
+
+ int nb_coeffs;
- int *ia = new int[nb_coeffs + 1];
- int *ja = new int[nb_coeffs + 1];
- scalar_t *ar = new scalar_t[nb_coeffs + 1];
+ if(absolute_proportion) {
+ nb_coeffs = nb_points * _nb_clusters + nb_points * _nb_clusters;
+ } else {
+ nb_coeffs = nb_points * _nb_clusters + nb_points * nb_classes * _nb_clusters;
+ }
+
+ int *coeff_row = new int[nb_coeffs + 1];
+ int *coeff_col = new int[nb_coeffs + 1];
+ scalar_t *coeff_wgt = new scalar_t[nb_coeffs + 1];
+
+ int n_coeff = 1;
scalar_t *nb_samples_per_class = new scalar_t[nb_classes];
for(int c = 0; c < nb_classes; c++) {
glp_set_prob_name(lp, "uninformative_lp_cluster_association");
glp_set_obj_dir(lp, GLP_MIN);
- glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
+ // We have one column per coefficient gamma
- for(int n = 1; n <= nb_points; n++) {
- glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
- }
+ glp_add_cols(lp, nb_points * _nb_clusters);
- for(int k = 1; k <= _nb_clusters; k++) {
- for(int c = 1; c <= nb_classes; c++) {
- int row = nb_points + (k - 1) * nb_classes + c;
- scalar_t tau = nb_samples_per_class[c-1] / scalar_t(_nb_clusters);
- // cout << "tau " << k << " " << c << " " << tau << endl;
- glp_set_row_bnds(lp, row, GLP_FX, tau, tau);
- }
- }
+ // The column for gamma[n][k] point 1<=n<=nb_points and cluster
+ // 1<=k<=_nb_clusters is nb_points * (k - 1) + n;
- glp_add_cols(lp, nb_points * _nb_clusters);
+ // The constraints (A) will be expressed by putting directly bounds
+ // on the variables (i.e. one per column). So we need one row per
+ // (B) constraint, and one per (C) constraint.
+
+ glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
+
+ // First, we set the weights for the objective function, and the (A)
+ // constraints on the individual gammas
for(int k = 1; k <= _nb_clusters; k++) {
for(int n = 1; n <= nb_points; n++) {
- int r = n + nb_points * (k - 1);
+ int col = n + nb_points * (k - 1);
- scalar_t dist = 0;
+ // The LP weight on the gammas for the global loss is the
+ // normalized distance of that sample to the centroid of that
+ // cluster
- for(int d = 0; d < _dim; d++) {
- dist += sq(_cluster_means[k-1][d] - points[n-1][d]) / (2 * _cluster_var[k-1][d]);
- dist += 0.5 * log(_cluster_var[k-1][d]);
- }
+ glp_set_obj_coef(lp, col, distance_to_centroid(points[n-1], k-1));
+
+ // The (A) constraints: Each column correspond to one of the
+ // gamma, and it has to be in [0,1]
- glp_set_obj_coef(lp, r, dist);
- glp_set_col_bnds(lp, r, GLP_DB, 0.0, 1.0);
+ glp_set_col_bnds(lp, col, GLP_DB, 0.0, 1.0);
}
}
- int l = 1;
+ // The (B) constraints: for each point, the sum of its gamma is
+ // equal to 1.0
for(int n = 1; n <= nb_points; n++) {
+ int row = n;
+ glp_set_row_bnds(lp, row, GLP_FX, 1.0, 1.0);
for(int k = 1; k <= _nb_clusters; k++) {
- int row = n;
- ia[l] = row;
- ja[l] = nb_points * (k - 1) + n;
- ar[l] = 1.0;
- l++;
+ coeff_row[n_coeff] = row;
+ coeff_col[n_coeff] = nb_points * (k - 1) + n;
+ coeff_wgt[n_coeff] = 1.0;
+ n_coeff++;
}
}
- for(int k = 1; k <= _nb_clusters; k++) {
- for(int c = 1; c <= nb_classes; c++) {
- int row = nb_points + (k - 1) * nb_classes + c;
- for(int n = 1; n <= nb_points; n++) {
- if(labels[n-1] == c - 1) {
- ia[l] = row;
- ja[l] = (k-1) * nb_points + n;
- ar[l] = 1.0;
- l++;
+ // The (C) constraints: For each pair cluster/class, the sum of the
+ // gammas for this cluster and this class is equal to the number of
+ // sample of that class, divided by the number of clusters
+
+ if(absolute_proportion) {
+ for(int k = 1; k <= _nb_clusters; k++) {
+ for(int c = 1; c <= nb_classes; c++) {
+ int row = nb_points + (k - 1) * nb_classes + c;
+ scalar_t tau = nb_samples_per_class[c-1] / scalar_t(_nb_clusters);
+ glp_set_row_bnds(lp, row, GLP_FX, tau, tau);
+ for(int n = 1; n <= nb_points; n++) {
+ if(labels[n-1] == c - 1) {
+ coeff_row[n_coeff] = row;
+ coeff_col[n_coeff] = (k-1) * nb_points + n;
+ coeff_wgt[n_coeff] = 1.0;
+ n_coeff++;
+ }
+ }
+ }
+ }
+ } else {
+ for(int k = 1; k <= _nb_clusters; k++) {
+ for(int c = 1; c <= nb_classes; c++) {
+ int row = nb_points + (k - 1) * nb_classes + c;
+ glp_set_row_bnds(lp, row, GLP_FX, 0.0, 0.0);
+ for(int n = 1; n <= nb_points; n++) {
+ coeff_row[n_coeff] = row;
+ coeff_col[n_coeff] = (k-1) * nb_points + n;
+ coeff_wgt[n_coeff] =
+ (labels[n-1] == c - 1 ? 1.0 : 0.0)
+ - scalar_t(nb_samples_per_class[c-1]) / scalar_t(nb_points);
+ n_coeff++;
}
}
}
}
- ASSERT(l == nb_coeffs + 1);
+ ASSERT(n_coeff == nb_coeffs + 1);
+
+ glp_load_matrix(lp, nb_coeffs, coeff_row, coeff_col, coeff_wgt);
- glp_load_matrix(lp, nb_coeffs, ia, ja, ar);
+ // Now a miracle occurs
glp_simplex(lp, NULL);
+ // We retrieve the result
+
scalar_t total_dist = glp_get_obj_val(lp);
for(int k = 1; k <= _nb_clusters; k++) {
}
}
- // { // ******************************* START ***************************
-// #warning Test code added on 2013 Feb 07 20:32:05
- // // for(int n = 0; n < nb_points; n++) {
- // // scalar_t sum = 0;
- // // for(int k = 0; k < _nb_clusters; k++) {
- // // ASSERT(gamma[n][k] >= 0 && gamma[n][k] <= 1);
- // // sum += gamma[n][k];
- // // }
- // // cout << sum << endl;
- // // }
-
- // scalar_t *sum_gamma = new scalar_t[nb_classes];
-
- // for(int k = 0; k < _nb_clusters; k++) {
- // for(int c = 0; c < nb_classes; c++) { sum_gamma[c] = 0.0; }
- // for(int n = 0; n < nb_points; n++) {
- // sum_gamma[labels[n]] += gamma[n][k];
- // }
- // cout << "CLUSTER" << k;
- // for(int c = 0; c < nb_classes; c++) {
- // cout << " " << sum_gamma[c];
- // }
- // cout << endl;
- // }
-
- // delete sum_gamma;
-
- // } // ******************************** END ****************************
-
delete[] nb_samples_per_class;
- delete[] ia;
- delete[] ja;
- delete[] ar;
+ delete[] coeff_row;
+ delete[] coeff_col;
+ delete[] coeff_wgt;
glp_delete_prob(lp);
return total_dist;
}
-void Clusterer::baseline_update_clusters(int nb_points, scalar_t **points, scalar_t **gamma) {
+void Clusterer::update_clusters(int nb_points, scalar_t **points, scalar_t **gamma) {
for(int k = 0; k < _nb_clusters; k++) {
for(int d = 0; d < _dim; d++) {
for(int d = 0; d < _dim; d++) {
if(sum_gamma >= 2) {
- _cluster_var[k][d] = (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
+ _cluster_var[k][d] =
+ (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
+ _cluster_var[k][d] = max(scalar_t(min_cluster_variance), _cluster_var[k][d]);
} else {
_cluster_var[k][d] = 1;
}
- _cluster_var[k][d] = max(0.01, _cluster_var[k][d]);
- _cluster_means[k][d] /= sum_gamma;
- }
- }
-}
-void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
- int *used = new int[nb_points];
- for(int k = 0; k < nb_points; k++) { used[k] = 0; }
- for(int k = 0; k < _nb_clusters; k++) {
- int l;
- do { l = int(drand48() * nb_points); } while(used[l]);
- for(int d = 0; d < _dim; d++) {
- _cluster_means[k][d] = points[l][d];
- _cluster_var[k][d] = 1.0;
+ _cluster_means[k][d] /= sum_gamma;
}
- used[l] = 1;
}
- delete[] used;
}
-void Clusterer::train(int nb_clusters, int dim,
+void Clusterer::train(int mode,
+ int nb_clusters, int dim,
int nb_points, scalar_t **points,
int nb_classes, int *labels,
int *cluster_associations) {
do {
pred_total_distance = total_distance;
- total_distance =
- // baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
- // baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
- uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
+
+ switch(mode) {
+
+ case STANDARD_ASSOCIATION:
+ total_distance =
+ baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
+ break;
+
+ case STANDARD_LP_ASSOCIATION:
+ total_distance =
+ baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
+ break;
+
+ case UNINFORMATIVE_LP_ASSOCIATION:
+ total_distance =
+ uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas, 0);
+ break;
+
+ case UNINFORMATIVE_LP_ASSOCIATION_ABSOLUTE:
+ total_distance =
+ uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas, 1);
+ break;
+
+ default:
+ cerr << "Unknown sample-cluster association mode." << endl;
+ abort();
+ }
+
cout << "TRAIN " << nb_rounds << " " << total_distance << endl;
- baseline_update_clusters(nb_points, points, gammas);
+ update_clusters(nb_points, points, gammas);
nb_rounds++;
} while(total_distance < min_iteration_improvement * pred_total_distance &&
nb_rounds < max_nb_iterations);
- {
- cout << "TOTAL_NB_SAMPLES";
- for(int c = 0; c < nb_classes; c++) {
- int nb_samples = 0;
- for(int n = 0; n < nb_points; n++) {
- if(labels[n] == c) {
- nb_samples++;
- }
- }
- cout << " " << nb_samples;
- }
- cout << endl;
-
- for(int k = 0; k < _nb_clusters; k++) {
- cout << "CLUSTER_GAMMA_SUM " << k << " :";
- for(int c = 0; c < nb_classes; c++) {
- scalar_t sum = 0.0;
- for(int n = 0; n < nb_points; n++) {
- if(labels[n] == c) {
- sum += gammas[n][k];
- }
+ if(cluster_associations) {
+ for(int n = 0; n < nb_points; n++) {
+ for(int k = 0; k < _nb_clusters; k++) {
+ if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
+ cluster_associations[n] = k;
}
- cout << " " << sum;
- }
- cout << endl;
- }
- }
-
- for(int n = 0; n < nb_points; n++) {
- for(int k = 0; k < _nb_clusters; k++) {
- if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
- cluster_associations[n] = k;
}
}
}
projects / clueless-kmeans.git / blobdiff