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