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

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

killed bug in SVM

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