source: trunk/lib/classifier/SVM.cc @ 571

Last change on this file since 571 was 571, checked in by Peter, 15 years ago

added function to set C-Parameter in SVM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.5 KB
Line 
1// $Id: SVM.cc 571 2006-04-06 13:08:17Z peter $
2
3#include <c++_tools/classifier/SVM.h>
4
5#include <c++_tools/classifier/DataLookup2D.h>
6#include <c++_tools/classifier/InputRanker.h>
7#include <c++_tools/gslapi/matrix.h>
8#include <c++_tools/gslapi/vector.h>
9#include <c++_tools/statistics/Averager.h>
10#include <c++_tools/statistics/Score.h>
11#include <c++_tools/random/random.h>
12
13#include <algorithm>
14#include <cassert>
15#include <cctype>
16#include <cmath>
17#include <limits>
18#include <utility>
19#include <vector>
20
21
22namespace theplu {
23namespace classifier { 
24
25  SVM::SVM(const KernelLookup& kernel, const Target& target)
26    : SupervisedClassifier(target),
27      alpha_(target.size(),0),
28      bias_(0),
29      C_inverse_(0),
30      kernel_(&kernel),
31      margin_(0),
32      max_epochs_(100000),
33      output_(target.size(),0),
34      owner_(false),
35      sample_(target.size()),
36      trained_(false),
37      tolerance_(0.00000001)
38  {
39#ifndef NDEBUG
40    for (size_t i=0; i<alpha_.size(); i++) 
41      for (size_t j=0; j<alpha_.size(); j++) 
42        assert(kernel(i,j)==kernel(j,i));
43    for (size_t i=0; i<alpha_.size(); i++) 
44      for (size_t j=0; j<alpha_.size(); j++) 
45        assert((*kernel_)(i,j)==(*kernel_)(j,i));
46    for (size_t i = 0; i<kernel_->rows(); i++)
47      for (size_t j = 0; j<kernel_->columns(); j++)
48        if (std::isnan((*kernel_)(i,j)))
49          std::cerr << "SVM: Found nan in kernel: " << i << " " 
50                    << j << std::endl;
51#endif       
52
53  }
54
55  SVM::SVM(const KernelLookup& kernel, const Target& target, 
56           statistics::Score& score, const size_t nof_inputs)
57    : SupervisedClassifier(target, &score, nof_inputs),
58      alpha_(target.size(),0),
59      bias_(0),
60      C_inverse_(0),
61      margin_(0),
62      max_epochs_(10000000),
63      output_(target.size(),0),
64      sample_(target.size()),
65      trained_(false),
66      tolerance_(0.00000001)
67  {
68
69    if (kernel.weighted()){
70      const MatrixLookup* data = kernel.data();
71      const MatrixLookup* weights = kernel.weights();
72      ranker_ = new InputRanker(*data, target, score, *weights);
73      delete data;
74      delete weights;
75    }
76    else{
77      const MatrixLookup* data = kernel.data();
78      ranker_ = new InputRanker(*data, target, score);
79      delete data;
80    }
81
82    std::vector<size_t> index;
83    index.reserve(nof_inputs);
84    for (size_t i=0; i<nof_inputs; i++) 
85      index.push_back(ranker_->id(i));
86    kernel_ = kernel.selected(index);
87    assert(kernel_->rows()==kernel.rows());
88    assert(kernel_->columns()==kernel.columns());
89    owner_ = true;
90   
91  }
92
93 
94
95  SVM::~SVM()
96  {
97    if (owner_)
98      delete kernel_;
99  }
100
101
102  void SVM::calculate_margin(void)
103  {
104    margin_ = 0;
105    for(size_t i = 0; i<alpha_.size(); ++i){
106      margin_ += alpha_(i)*target(i)*kernel_mod(i,i)*alpha_(i)*target(i);
107      for(size_t j = i+1; j<alpha_.size(); ++j)
108        margin_ += 2*alpha_(i)*target(i)*kernel_mod(i,j)*alpha_(j)*target(j);
109    }
110  }
111
112
113  SupervisedClassifier* SVM::make_classifier(const DataLookup2D& data, 
114                                             const Target& target) const
115  {
116    assert(data.rows()==data.columns());
117    assert(data.columns()==target.size());
118    // Peter, should check success of dynamic_cast
119    const KernelLookup& tmp = dynamic_cast<const KernelLookup&>(data);
120    SVM* sc;
121    if (ranker_)
122      sc = new SVM(tmp,target,*score_,nof_inputs_);
123    else
124      sc = new SVM(tmp,target);
125
126
127    //Copy those variables possible to modify from outside
128    // Peter, in particular C
129    return sc;
130  }
131
132  void SVM::predict(const DataLookup2D& input, gslapi::matrix& prediction) const
133  {
134    // Peter, should check success of dynamic_cast
135    const KernelLookup input_kernel = dynamic_cast<const KernelLookup&>(input);
136
137    const KernelLookup* kernel_pointer;
138    if (ranker_) {// feature selection
139      std::vector<size_t> index;
140      index.reserve(nof_inputs_);
141      for (size_t i=0; i<nof_inputs_; i++) 
142        index.push_back(ranker_->id(i));
143      kernel_pointer = input_kernel.selected(index);
144    }
145    else 
146      kernel_pointer = &input_kernel;
147
148    assert(input.rows()==alpha_.size());
149    prediction = gslapi::matrix(2,input.columns(),0);
150    for (size_t i = 0; i<input.columns(); i++){
151      for (size_t j = 0; j<input.rows(); j++){
152        prediction(0,i) += target(j)*alpha_(j)*(*kernel_pointer)(j,i);
153        assert(target(j));
154      }
155      prediction(0,i) = margin_ * (prediction(0,i) + bias_);
156    }
157
158    for (size_t i = 0; i<prediction.columns(); i++)
159      prediction(1,i) = -prediction(0,i);
160   
161    if (ranker_)
162      delete kernel_pointer;
163    assert(prediction(0,0));
164  }
165
166  double SVM::predict(const DataLookup1D& x) const
167  {
168    double y=0;
169    for (size_t i=0; i<alpha_.size(); i++)
170      y += alpha_(i)*target_(i)*kernel_->element(x,i);
171
172    return margin_*(y+bias_);
173  }
174
175  double SVM::predict(const DataLookup1D& x, const DataLookup1D& w) const
176  {
177    double y=0;
178    for (size_t i=0; i<alpha_.size(); i++)
179      y += alpha_(i)*target_(i)*kernel_->element(x,w,i);
180
181    return margin_*(y+bias_);
182  }
183
184  bool SVM::train(void) 
185  {
186    // initializing variables for optimization
187    assert(target_.size()==kernel_->rows());
188    assert(target_.size()==alpha_.size());
189
190    sample_.init(alpha_,tolerance_);
191    gslapi::vector  E(target_.size(),0);
192    for (size_t i=0; i<E.size(); i++) {
193      for (size_t j=0; j<E.size(); j++) 
194        E(i) += kernel_mod(i,j)*target(j)*alpha_(j);
195      E(i)-=target(i);
196    }
197    assert(target_.size()==E.size());
198    assert(target_.size()==sample_.size());
199
200    unsigned long int epochs = 0;
201    double alpha_new2;
202    double alpha_new1;
203    double u;
204    double v;
205
206    // Training loop
207    while(choose(E)) {
208      bounds(u,v);       
209      double k = ( kernel_mod(sample_.value_first(), sample_.value_first()) + 
210                   kernel_mod(sample_.value_second(), sample_.value_second()) - 
211                   2*kernel_mod(sample_.value_first(), sample_.value_second()));
212     
213      double alpha_old1=alpha_(sample_.value_first());
214      double alpha_old2=alpha_(sample_.value_second());
215      alpha_new2 = ( alpha_(sample_.value_second()) + 
216                     target(sample_.value_second())*
217                     ( E(sample_.value_first())-E(sample_.value_second()) )/k );
218     
219      if (alpha_new2 > v)
220        alpha_new2 = v;
221      else if (alpha_new2<u)
222        alpha_new2 = u;
223     
224      // Updating the alphas
225      // if alpha is 'zero' make the sample a non-support vector
226      if (alpha_new2 < tolerance_){
227        sample_.nsv_second();
228      }
229      else{
230        sample_.sv_second();
231      }
232     
233     
234      alpha_new1 = (alpha_(sample_.value_first()) + 
235                    (target(sample_.value_first()) * 
236                     target(sample_.value_second()) * 
237                     (alpha_(sample_.value_second()) - alpha_new2) ));
238           
239      // if alpha is 'zero' make the sample a non-support vector
240      if (alpha_new1 < tolerance_){
241        sample_.nsv_first();
242      }
243      else
244        sample_.sv_first();
245     
246      alpha_(sample_.value_first()) = alpha_new1;
247      alpha_(sample_.value_second()) = alpha_new2;
248     
249      // update E vector
250      // Peter, perhaps one should only update SVs, but what happens in choose?
251      for (size_t i=0; i<E.size(); i++) {
252        E(i)+=( kernel_mod(i,sample_.value_first())*
253                target(sample_.value_first()) *
254                (alpha_new1-alpha_old1) );
255        E(i)+=( kernel_mod(i,sample_.value_second())*
256                target(sample_.value_second()) *
257                (alpha_new2-alpha_old2) );
258      }
259           
260      epochs++; 
261      if (epochs>max_epochs_){
262        std::cerr << "WARNING: SVM: maximal number of epochs reached.\n";
263        calculate_bias();
264        calculate_margin();
265        return false;
266      }
267    }
268    calculate_margin();
269    trained_ = calculate_bias();
270    return trained_;
271  }
272
273
274  bool SVM::choose(const theplu::gslapi::vector& E)
275  {
276    // First check for violation among SVs
277    // E should be the same for all SVs
278    // Choose that pair having largest violation/difference.
279    sample_.update_second(0);
280    sample_.update_first(0);
281    if (sample_.nof_sv()>1){
282
283      double max = E(sample_(0));
284      double min = max;
285      for (size_t i=1; i<sample_.nof_sv(); i++){ 
286        assert(alpha_(sample_(i))>tolerance_);
287        if (E(sample_(i)) > max){
288          max = E(sample_(i));
289          sample_.update_second(i);
290        }
291        else if (E(sample_(i))<min){
292          min = E(sample_(i));
293          sample_.update_first(i);
294        }
295      }
296      assert(alpha_(sample_.value_first())>tolerance_);
297      assert(alpha_(sample_.value_second())>tolerance_);
298
299      if (E(sample_.value_second()) - E(sample_.value_first()) > 2*tolerance_){
300        return true;
301      }
302     
303      // If no violation check among non-support vectors
304      sample_.shuffle();
305      for (size_t i=sample_.nof_sv(); i<sample_.size();i++){
306        if (target_.binary(sample_(i))){
307          if(E(sample_(i)) < E(sample_.value_first()) - 2*tolerance_){
308            sample_.update_second(i);
309            return true;
310          }
311        }
312        else{
313          if(E(sample_(i)) > E(sample_.value_second()) + 2*tolerance_){
314            sample_.update_first(i);
315            return true;
316          }
317        }
318      }
319    }
320
321    // if no support vectors - special case
322    else{
323      // to avoid getting stuck we shuffle
324      sample_.shuffle();
325      for (size_t i=0; i<sample_.size(); i++) {
326        if (target(sample_(i))==1){
327          for (size_t j=0; j<sample_.size(); j++) {
328            if ( target(sample_(j))==-1 && 
329                 E(sample_(i)) < E(sample_(j))+2*tolerance_ ){
330              sample_.update_first(i);
331              sample_.update_second(j);
332              return true;
333            }
334          }
335        }
336      }
337    }
338   
339    // If there is no violation then we should stop training
340    return false;
341
342  }
343 
344 
345  void SVM::bounds( double& u, double& v) const
346  {
347    if (target(sample_.value_first())!=target(sample_.value_second())) {
348      if (alpha_(sample_.value_second()) > alpha_(sample_.value_first())) {
349        v = std::numeric_limits<double>::max();
350        u = alpha_(sample_.value_second()) - alpha_(sample_.value_first());
351      }
352      else {
353        v = (std::numeric_limits<double>::max() - 
354             alpha_(sample_.value_first()) + 
355             alpha_(sample_.value_second()));
356        u = 0;
357      }
358    }
359    else {       
360      if (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) > 
361           std::numeric_limits<double>::max()) {
362        u = (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) - 
363              std::numeric_limits<double>::max());
364        v =  std::numeric_limits<double>::max();   
365      }
366      else {
367        u = 0;
368        v = alpha_(sample_.value_first()) + alpha_(sample_.value_second());
369      }
370    }
371  }
372 
373  bool SVM::calculate_bias(void)
374  {
375
376    // calculating output without bias
377    for (size_t i=0; i<output_.size(); i++) {
378      output_(i)=0;
379      for (size_t j=0; j<output_.size(); j++) 
380        output_(i)+=alpha_(j)*target(j) * (*kernel_)(i,j);
381    }
382
383    if (!sample_.nof_sv()){
384      std::cerr << "SVM::train() error: " 
385                << "Cannot calculate bias because there is no support vector" 
386                << std::endl;
387      return false;
388    }
389
390    // For samples with alpha>0, we have: target*output=1-alpha/C
391    bias_=0;
392    for (size_t i=0; i<sample_.nof_sv(); i++) 
393      bias_+= ( target(sample_(i)) * (1-alpha_(sample_(i))*C_inverse_) - 
394                output_(sample_(i)) );
395    bias_=bias_/sample_.nof_sv();
396    for (size_t i=0; i<output_.size(); i++) 
397      output_(i) += bias_;
398     
399    return true;
400  }
401
402  void SVM::set_C(const double C)
403  {
404    C_inverse_ = 1/C;
405  }
406
407  Index::Index(void)
408    : nof_sv_(0), vec_(std::vector<size_t>(0))
409  {
410  }
411
412  Index::Index(const size_t n)
413    : nof_sv_(0), vec_(std::vector<size_t>(n))
414  {
415    for (size_t i=0; i<vec_.size(); i++)
416      vec_[i]=i;
417  }
418
419  void Index::init(const gslapi::vector& alpha, const double tol)
420  {
421    nof_sv_=0;
422    size_t nof_nsv=0;
423    for (size_t i=0; i<alpha.size(); i++) 
424      if (alpha(i)<tol){
425        nof_nsv++;
426        vec_[vec_.size()-nof_nsv]=i;
427      }
428      else{
429        vec_[nof_sv_]=i;
430        nof_sv_++;
431      }
432    assert(nof_sv_+nof_nsv==vec_.size());
433
434  }
435
436  void Index::sv_first(void)
437  {
438    // if already sv, do nothing
439    if (index_first_<nof_sv())
440      return;
441
442    // swap elements
443    if(index_second_==nof_sv_){
444      index_second_=index_first_;
445    }
446    vec_[index_first_]=vec_[nof_sv_];
447    vec_[nof_sv_]=value_first_;
448    index_first_ = nof_sv_;
449
450    nof_sv_++;
451
452  }
453
454  void Index::sv_second(void)
455  {
456    // if already sv, do nothing
457    if (index_second_<nof_sv())
458      return;
459
460    // swap elements
461    if(index_first_==nof_sv_){
462      index_first_=index_second_;
463    }
464
465    vec_[index_second_]=vec_[nof_sv_];
466    vec_[nof_sv_]=value_second_;
467    index_second_=nof_sv_;
468
469    nof_sv_++;
470  }
471
472  void Index::nsv_first(void)
473  {
474    // if already nsv, do nothing
475    if ( !(index_first_<nof_sv()) )
476      return;
477   
478    if(index_second_==nof_sv_-1)
479      index_second_=index_first_;
480    vec_[index_first_]=vec_[nof_sv_-1];
481    vec_[nof_sv_-1]=value_first_;
482    index_first_=nof_sv_-1;
483   
484    nof_sv_--;
485  }
486
487  void Index::nsv_second(void)
488  {
489    // if already nsv, do nothing
490    if ( !(index_second_<nof_sv()) )
491      return;
492
493    if(index_first_==nof_sv_-1)
494      index_first_=index_second_;
495    vec_[index_second_]=vec_[nof_sv_-1];
496    vec_[nof_sv_-1]=value_second_;
497    index_second_ = nof_sv_-1;
498   
499    nof_sv_--;
500  }
501
502
503  void Index::shuffle(void)
504  {
505    random::DiscreteUniform a;
506    random_shuffle(vec_.begin()+nof_sv_, vec_.end(), a); 
507  }
508
509  void Index::update_first(const size_t i)
510  {
511    assert(i<size());
512    index_first_=i;
513    value_first_=vec_[i];
514  }
515
516  void Index::update_second(const size_t i)
517  {
518    assert(i<size());
519    index_second_=i;
520    value_second_=vec_[i];
521  }
522
523}} // of namespace classifier and namespace theplu
Note: See TracBrowser for help on using the repository browser.