source: trunk/c++_tools/classifier/SVM.cc @ 604

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

ref #96 Changed InputRanker? to return vector of index rather than element. Also added draft to FeatureSelection? class.

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