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

Last change on this file since 559 was 559, checked in by Peter, 17 years ago

some changes in EB

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