loss_machine.cc

   1
   2 ///////////////////////////////////////////////////////////////////////////
   3 // This program is free software: you can redistribute it and/or modify  //
   4 // it under the terms of the version 3 of the GNU General Public License //
   5 // as published by the Free Software Foundation.                         //
   6 //                                                                       //
   7 // This program is distributed in the hope that it will be useful, but   //
   8 // WITHOUT ANY WARRANTY; without even the implied warranty of            //
   9 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU      //
  10 // General Public License for more details.                              //
  11 //                                                                       //
  12 // You should have received a copy of the GNU General Public License     //
  13 // along with this program. If not, see <http://www.gnu.org/licenses/>.  //
  14 //                                                                       //
  15 // Written by Francois Fleuret                                           //
  16 // (C) Idiap Research Institute                                          //
  17 //                                                                       //
  18 // Contact <francois.fleuret@idiap.ch> for comments & bug reports        //
  19 ///////////////////////////////////////////////////////////////////////////
  20
  21 #include "tools.h"
  22 #include "loss_machine.h"
  23
  24 LossMachine::LossMachine(int loss_type) {
  25   _loss_type = loss_type;
  26 }
  27
  28 void LossMachine::get_loss_derivatives(SampleSet *samples,
  29                                        scalar_t *responses,
  30                                        scalar_t *derivatives) {
  31
  32   switch(_loss_type) {
  33
  34   case LOSS_EXPONENTIAL:
  35     {
  36       for(int n = 0; n < samples->nb_samples(); n++) {
  37         derivatives[n] =
  38           - samples->label(n) * exp( - samples->label(n) * responses[n]);
  39       }
  40     }
  41     break;
  42
  43   case LOSS_EV_REGULARIZED:
  44     {
  45       scalar_t sum_pos = 0, sum_sq_pos = 0, nb_pos = 0, m_pos, v_pos;
  46       scalar_t sum_neg = 0, sum_sq_neg = 0, nb_neg = 0, m_neg, v_neg;
  47
  48       for(int n = 0; n < samples->nb_samples(); n++) {
  49         if(samples->label(n) > 0) {
  50           sum_pos += responses[n];
  51           sum_sq_pos += sq(responses[n]);
  52           nb_pos += 1.0;
  53         }
  54         else if(samples->label(n) < 0) {
  55           sum_neg += responses[n];
  56           sum_sq_neg += sq(responses[n]);
  57           nb_neg += 1.0;
  58         }
  59       }
  60
  61       m_pos = sum_pos / nb_pos;
  62       v_pos = sum_sq_pos/(nb_pos - 1) - sq(sum_pos)/(nb_pos * (nb_pos - 1));
  63
  64       scalar_t loss_pos = nb_pos * exp(v_pos/2 - m_pos);
  65
  66       m_neg = sum_neg / nb_neg;
  67       v_neg = sum_sq_neg/(nb_neg - 1) - sq(sum_neg)/(nb_neg * (nb_neg - 1));
  68
  69       scalar_t loss_neg = nb_neg * exp(v_neg/2 + m_neg);
  70
  71       for(int n = 0; n < samples->nb_samples(); n++) {
  72         if(samples->label(n) > 0) {
  73           derivatives[n] =
  74             ( - 1/nb_pos + (responses[n] - m_pos)/(nb_pos - 1)) * loss_pos;
  75         } else if(samples->label(n) < 0) {
  76           derivatives[n] =
  77             (   1/nb_neg + (responses[n] - m_neg)/(nb_neg - 1)) * loss_neg;
  78         }
  79       }
  80     }
  81
  82     break;
  83
  84   case LOSS_HINGE:
  85     {
  86       for(int n = 0; n < samples->nb_samples(); n++) {
  87         if(samples->label(n) != 0 && samples->label(n) * responses[n] < 1)
  88           derivatives[n] = 1;
  89         else
  90           derivatives[n] = 0;
  91       }
  92     }
  93     break;
  94
  95   case LOSS_LOGISTIC:
  96     {
  97       for(int n = 0; n < samples->nb_samples(); n++) {
  98         if(samples->label(n) == 0)
  99           derivatives[n] = 0.0;
 100         else
 101           derivatives[n] = samples->label(n) * 1/(1 + exp(samples->label(n) * responses[n]));
 102       }
 103     }
 104     break;
 105
 106   default:
 107     cerr << "Unknown loss type in BoostedClassifier::get_loss_derivatives."
 108          << endl;
 109     exit(1);
 110   }
 111
 112 }
 113
 114 scalar_t LossMachine::loss(SampleSet *samples, scalar_t *responses) {
 115   scalar_t l = 0;
 116
 117   switch(_loss_type) {
 118
 119   case LOSS_EXPONENTIAL:
 120     {
 121       for(int n = 0; n < samples->nb_samples(); n++) {
 122         l += exp( - samples->label(n) * responses[n]);
 123         ASSERT(!isinf(l));
 124       }
 125     }
 126     break;
 127
 128   case LOSS_EV_REGULARIZED:
 129     {
 130       scalar_t sum_pos = 0, sum_sq_pos = 0, nb_pos = 0, m_pos, v_pos;
 131       scalar_t sum_neg = 0, sum_sq_neg = 0, nb_neg = 0, m_neg, v_neg;
 132
 133       for(int n = 0; n < samples->nb_samples(); n++) {
 134         if(samples->label(n) > 0) {
 135           sum_pos += responses[n];
 136           sum_sq_pos += sq(responses[n]);
 137           nb_pos += 1.0;
 138         } else if(samples->label(n) < 0) {
 139           sum_neg += responses[n];
 140           sum_sq_neg += sq(responses[n]);
 141           nb_neg += 1.0;
 142         }
 143       }
 144
 145       l = 0;
 146
 147       if(nb_pos > 0) {
 148         m_pos = sum_pos / nb_pos;
 149         v_pos = sum_sq_pos/(nb_pos - 1) - sq(sum_pos)/(nb_pos * (nb_pos - 1));
 150         l += nb_pos * exp(v_pos/2 - m_pos);
 151       }
 152
 153       if(nb_neg > 0) {
 154         m_neg = sum_neg / nb_neg;
 155         v_neg = sum_sq_neg/(nb_neg - 1) - sq(sum_neg)/(nb_neg * (nb_neg - 1));
 156         l += nb_neg * exp(v_neg/2 + m_neg);
 157       }
 158
 159     }
 160     break;
 161
 162   case LOSS_HINGE:
 163     {
 164       for(int n = 0; n < samples->nb_samples(); n++) {
 165         if(samples->label(n) != 0) {
 166           if(samples->label(n) * responses[n] < 1)
 167             l += (1 - samples->label(n) * responses[n]);
 168         }
 169       }
 170     }
 171     break;
 172
 173   case LOSS_LOGISTIC:
 174     {
 175       for(int n = 0; n < samples->nb_samples(); n++) {
 176         if(samples->label(n) != 0) {
 177           scalar_t u = - samples->label(n) * responses[n];
 178           if(u > 20) {
 179             l += u;
 180           } if(u > -20) {
 181             l += log(1 + exp(u));
 182           }
 183         }
 184       }
 185     }
 186     break;
 187
 188   default:
 189     cerr << "Unknown loss type in LossMachine::loss." << endl;
 190     exit(1);
 191   }
 192
 193   return l;
 194 }
 195
 196 scalar_t LossMachine::optimal_weight(SampleSet *sample_set,
 197                                      scalar_t *weak_learner_responses,
 198                                      scalar_t *current_responses) {
 199
 200   switch(_loss_type) {
 201
 202   case LOSS_EXPONENTIAL:
 203     {
 204       scalar_t num = 0, den = 0, z;
 205       for(int n = 0; n < sample_set->nb_samples(); n++) {
 206         z = sample_set->label(n) * weak_learner_responses[n];
 207         if(z > 0) {
 208           num += exp( - sample_set->label(n) * current_responses[n]);
 209         } else if(z < 0) {
 210           den += exp( - sample_set->label(n) * current_responses[n]);
 211         }
 212       }
 213
 214       return 0.5 * log(num / den);
 215     }
 216     break;
 217
 218   case LOSS_EV_REGULARIZED:
 219     {
 220
 221       scalar_t u = 0, du = -0.1;
 222       scalar_t *responses = new scalar_t[sample_set->nb_samples()];
 223
 224       scalar_t l, prev_l = -1;
 225
 226       const scalar_t minimum_delta_for_optimization = 1e-5;
 227
 228       scalar_t shift = 0;
 229
 230       {
 231         scalar_t sum_pos = 0, sum_sq_pos = 0, nb_pos = 0, m_pos, v_pos;
 232         scalar_t sum_neg = 0, sum_sq_neg = 0, nb_neg = 0, m_neg, v_neg;
 233
 234         for(int n = 0; n < sample_set->nb_samples(); n++) {
 235           if(sample_set->label(n) > 0) {
 236             sum_pos += responses[n];
 237             sum_sq_pos += sq(responses[n]);
 238             nb_pos += 1.0;
 239           } else if(sample_set->label(n) < 0) {
 240             sum_neg += responses[n];
 241             sum_sq_neg += sq(responses[n]);
 242             nb_neg += 1.0;
 243           }
 244         }
 245
 246         if(nb_pos > 0) {
 247           m_pos = sum_pos / nb_pos;
 248           v_pos = sum_sq_pos/(nb_pos - 1) - sq(sum_pos)/(nb_pos * (nb_pos - 1));
 249           shift = max(shift, v_pos/2 - m_pos);
 250         }
 251
 252         if(nb_neg > 0) {
 253           m_neg = sum_neg / nb_neg;
 254           v_neg = sum_sq_neg/(nb_neg - 1) - sq(sum_neg)/(nb_neg * (nb_neg - 1));
 255           shift = max(shift, v_neg/2 + m_neg);
 256         }
 257
 258 //         (*global.log_stream) << "nb_pos = " << nb_pos << " nb_neg = " << nb_neg << endl;
 259
 260       }
 261
 262       int nb = 0;
 263
 264       while(nb < 100 && abs(du) > minimum_delta_for_optimization) {
 265         nb++;
 266
 267 //         (*global.log_stream) << "l = " << l << " u = " << u << " du = " << du << endl;
 268
 269         u += du;
 270         for(int s = 0; s < sample_set->nb_samples(); s++) {
 271           responses[s] = current_responses[s] + u * weak_learner_responses[s] ;
 272         }
 273
 274         {
 275           scalar_t sum_pos = 0, sum_sq_pos = 0, nb_pos = 0, m_pos, v_pos;
 276           scalar_t sum_neg = 0, sum_sq_neg = 0, nb_neg = 0, m_neg, v_neg;
 277
 278           for(int n = 0; n < sample_set->nb_samples(); n++) {
 279             if(sample_set->label(n) > 0) {
 280               sum_pos += responses[n];
 281               sum_sq_pos += sq(responses[n]);
 282               nb_pos += 1.0;
 283             } else if(sample_set->label(n) < 0) {
 284               sum_neg += responses[n];
 285               sum_sq_neg += sq(responses[n]);
 286               nb_neg += 1.0;
 287             }
 288           }
 289
 290           l = 0;
 291
 292           if(nb_pos > 0) {
 293             m_pos = sum_pos / nb_pos;
 294             v_pos = sum_sq_pos/(nb_pos - 1) - sq(sum_pos)/(nb_pos * (nb_pos - 1));
 295             l += nb_pos * exp(v_pos/2 - m_pos - shift);
 296           }
 297
 298           if(nb_neg > 0) {
 299             m_neg = sum_neg / nb_neg;
 300             v_neg = sum_sq_neg/(nb_neg - 1) - sq(sum_neg)/(nb_neg * (nb_neg - 1));
 301             l += nb_neg * exp(v_neg/2 + m_neg - shift);
 302           }
 303
 304         }
 305
 306         if(l > prev_l) du = du * -0.25;
 307         prev_l = l;
 308       }
 309
 310       delete[] responses;
 311
 312       return u;
 313     }
 314
 315   case LOSS_HINGE:
 316   case LOSS_LOGISTIC:
 317     {
 318
 319       scalar_t u = 0, du = -0.1;
 320       scalar_t *responses = new scalar_t[sample_set->nb_samples()];
 321
 322       scalar_t l, prev_l = -1;
 323
 324       const scalar_t minimum_delta_for_optimization = 1e-5;
 325
 326       int n = 0;
 327       while(n < 100 && abs(du) > minimum_delta_for_optimization) {
 328         n++;
 329         u += du;
 330         for(int s = 0; s < sample_set->nb_samples(); s++) {
 331           responses[s] = current_responses[s] + u * weak_learner_responses[s] ;
 332         }
 333         l = loss(sample_set, responses);
 334         if(l > prev_l) du = du * -0.25;
 335         prev_l = l;
 336       }
 337
 338       (*global.log_stream) << "END l = " << l << " du = " << du << endl;
 339
 340       delete[] responses;
 341
 342       return u;
 343     }
 344
 345   default:
 346     cerr << "Unknown loss type in LossMachine::optimal_weight." << endl;
 347     exit(1);
 348   }
 349
 350 }
 351
 352 void LossMachine::subsample(int nb, scalar_t *labels, scalar_t *responses,
 353                             int nb_to_sample, int *sample_nb_occurences, scalar_t *sample_responses,
 354                             int allow_duplicates) {
 355
 356   switch(_loss_type) {
 357
 358   case LOSS_EXPONENTIAL:
 359     {
 360       scalar_t *weights = new scalar_t[nb];
 361
 362       for(int n = 0; n < nb; n++) {
 363         if(labels[n] == 0) {
 364           weights[n] = 0;
 365         } else {
 366           weights[n] = exp( - labels[n] * responses[n]);
 367         }
 368         sample_nb_occurences[n] = 0;
 369         sample_responses[n] = 0.0;
 370       }
 371
 372       scalar_t total_weight;
 373       int nb_sampled = 0, sum_sample_nb_occurences = 0;
 374
 375       int *sampled_indexes = new int[nb_to_sample];
 376
 377       (*global.log_stream) << "Sampling " << nb_to_sample << " samples." << endl;
 378
 379       do {
 380         total_weight = robust_sampling(nb,
 381                                        weights,
 382                                        nb_to_sample,
 383                                        sampled_indexes);
 384
 385         for(int k = 0; nb_sampled < nb_to_sample && k < nb_to_sample; k++) {
 386           int i = sampled_indexes[k];
 387           if(allow_duplicates || sample_nb_occurences[i] == 0) nb_sampled++;
 388           sample_nb_occurences[i]++;
 389           sum_sample_nb_occurences++;
 390         }
 391       } while(nb_sampled < nb_to_sample);
 392
 393       (*global.log_stream) << "nb_sampled = " << nb_sampled << " nb_to_sample = " << nb_to_sample << endl;
 394
 395       (*global.log_stream) << "Done." << endl;
 396
 397       delete[] sampled_indexes;
 398
 399       scalar_t unit_weight = log(total_weight / scalar_t(sum_sample_nb_occurences));
 400
 401       for(int n = 0; n < nb; n++) {
 402         if(sample_nb_occurences[n] > 0) {
 403           if(allow_duplicates) {
 404             sample_responses[n] = - labels[n] * unit_weight;
 405           } else {
 406             sample_responses[n] = - labels[n] * (unit_weight + log(scalar_t(sample_nb_occurences[n])));
 407             sample_nb_occurences[n] = 1;
 408           }
 409         }
 410       }
 411
 412       delete[] weights;
 413
 414     }
 415     break;
 416
 417   default:
 418     cerr << "Unknown loss type in LossMachine::resample." << endl;
 419     exit(1);
 420   }
 421
 422
 423 }