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

Last change on this file since 451 was 451, checked in by Peter, 16 years ago

directory svm -> classifier

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.0 KB
Line 
1// $Id: SVM.cc 451 2005-12-15 21:09:12Z peter $
2
3
4#include <c++_tools/classifier/SVM.h>
5
6#include <c++_tools/classifier/Kernel.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/random/random.h>
11
12#include <iostream>
13#include <algorithm>
14#include <cmath>
15#include <limits>
16#include <utility>
17#include <vector>
18
19
20namespace theplu {
21namespace classifier { 
22
23  Index::Index(void)
24    : nof_sv_(0), vec_(std::vector<size_t>(0))
25  {
26  }
27
28  Index::Index(const size_t n)
29    : nof_sv_(0), vec_(std::vector<size_t>(n))
30  {
31    for (size_t i=0; i<vec_.size(); i++)
32      vec_[i]=i;
33  }
34
35  void Index::init(const gslapi::vector& alpha, const double tol)
36  {
37    nof_sv_=0;
38    size_t nof_nsv=0;
39    for (size_t i=0; i<alpha.size(); i++) 
40      if (alpha(i)<tol){
41        nof_nsv++;
42        vec_[vec_.size()-nof_nsv]=i;
43      }
44      else{
45        vec_[nof_sv_]=i;
46        nof_sv_++;
47      }
48    assert(nof_sv_+nof_nsv==vec_.size());
49
50  }
51
52  void Index::sv_first(void)
53  {
54    // if already sv, do nothing
55    if (index_first_<nof_sv())
56      return;
57
58    // swap elements
59    if(index_second_==nof_sv_){
60      index_second_=index_first_;
61    }
62    vec_[index_first_]=vec_[nof_sv_];
63    vec_[nof_sv_]=value_first_;
64    index_first_ = nof_sv_;
65
66    nof_sv_++;
67
68  }
69
70  void Index::sv_second(void)
71  {
72    // if already sv, do nothing
73    if (index_second_<nof_sv())
74      return;
75
76    // swap elements
77    if(index_first_==nof_sv_){
78      index_first_=index_second_;
79    }
80
81    vec_[index_second_]=vec_[nof_sv_];
82    vec_[nof_sv_]=value_second_;
83    index_second_=nof_sv_;
84
85    nof_sv_++;
86  }
87
88  void Index::nsv_first(void)
89  {
90    // if already nsv, do nothing
91    if ( !(index_first_<nof_sv()) )
92      return;
93   
94    if(index_second_==nof_sv_-1)
95      index_second_=index_first_;
96    vec_[index_first_]=vec_[nof_sv_-1];
97    vec_[nof_sv_-1]=value_first_;
98    index_first_=nof_sv_-1;
99   
100    nof_sv_--;
101  }
102
103  void Index::nsv_second(void)
104  {
105    // if already nsv, do nothing
106    if ( !(index_second_<nof_sv()) )
107      return;
108
109    if(index_first_==nof_sv_-1)
110      index_first_=index_second_;
111    vec_[index_second_]=vec_[nof_sv_-1];
112    vec_[nof_sv_-1]=value_second_;
113    index_second_ = nof_sv_-1;
114   
115    nof_sv_--;
116  }
117
118
119  void Index::shuffle(void)
120  {
121    random::DiscreteUniform a;
122    random_shuffle(vec_.begin()+nof_sv_, vec_.end(), a); 
123  }
124
125  void Index::update_first(const size_t i)
126  {
127    assert(i<n());
128    index_first_=i;
129    value_first_=vec_[i];
130  }
131
132  void Index::update_second(const size_t i)
133  {
134    assert(i<n());
135    index_second_=i;
136    value_second_=vec_[i];
137  }
138
139  SVM::SVM(const Kernel& kernel, const gslapi::vector& target)
140  : alpha_(target.size(),0),
141    bias_(0),
142    C_inverse_(0),
143    kernel_(&kernel),
144    max_epochs_(10000000),
145    output_(target.size(),0),
146    sample_(target.size()),
147    target_(target),
148    trained_(false),
149    tolerance_(0.00000001)
150  {
151    if (max_epochs_>ULONG_MAX)
152      max_epochs_=ULONG_MAX;
153  }
154
155  SVM::SVM(const SVM& other)
156  : alpha_(other.alpha_),
157    bias_(other.bias_),
158    C_inverse_(other.C_inverse_),
159    kernel_(other.kernel_),
160    max_epochs_(other.max_epochs_),
161    output_(other.output_),
162    sample_(other.sample_),
163    target_(other.target_),
164    trained_(other.trained_),
165    tolerance_(other.tolerance_)
166  {
167  }
168
169  bool SVM::train(void) 
170  {
171    // initializing variables for optimization
172    assert(target_.size()==kernel_->size());
173    assert(target_.size()==alpha_.size());
174
175    sample_.init(alpha_,tolerance_);
176    gslapi::vector  E(target_.size());
177    for (size_t i=0; i<E.size(); i++) {
178      E(i)=0;
179      for (size_t j=0; j<E.size(); j++) 
180        E(i) += kernel_mod(i,j)*target_(j)*alpha_(j);
181      E(i)=E(i)-target_(i);
182    }
183    assert(target_.size()==E.size());
184    assert(target_.size()==sample_.n());
185
186    unsigned long int epochs = 0;
187    double alpha_new2;
188    double alpha_new1;
189    double u;
190    double v;
191
192    // Training loop
193    while(choose(E)) {
194      bounds(u,v);       
195      double k = ( kernel_mod(sample_.value_first(), sample_.value_first()) + 
196                   kernel_mod(sample_.value_second(), sample_.value_second()) - 
197                   2*kernel_mod(sample_.value_first(), sample_.value_second()));
198     
199      double alpha_old1=alpha_(sample_.value_first());
200      double alpha_old2=alpha_(sample_.value_second());
201
202      alpha_new2 = ( alpha_(sample_.value_second()) + 
203                     target_(sample_.value_second())*
204                     ( E(sample_.value_first())-E(sample_.value_second()) )/k );
205     
206      if (alpha_new2 > v)
207        alpha_new2 = v;
208      else if (alpha_new2<u)
209        alpha_new2 = u;
210     
211     
212      // Updating the alphas
213      // if alpha is 'zero' make the sample a non-support vector
214      if (alpha_new2 < tolerance_){
215        sample_.nsv_second();
216      }
217      else{
218        sample_.sv_second();
219      }
220     
221     
222      alpha_new1 = (alpha_(sample_.value_first()) + 
223                    (target_(sample_.value_first()) * 
224                     target_(sample_.value_second()) * 
225                     (alpha_(sample_.value_second()) - alpha_new2) ));
226           
227      // if alpha is 'zero' make the sample a non-support vector
228      if (alpha_new1 < tolerance_){
229        sample_.nsv_first();
230      }
231      else
232        sample_.sv_first();
233     
234      alpha_(sample_.value_first()) = alpha_new1;
235      alpha_(sample_.value_second()) = alpha_new2;
236     
237      // update E vector
238      // Peter, perhaps one should only update SVs, but what happens in choose?
239      for (size_t i=0; i<E.size(); i++) {
240        E(i)+=( kernel_mod(i,sample_.value_first())*
241                target_(sample_.value_first()) *
242                (alpha_new1-alpha_old1) );
243        E(i)+=( kernel_mod(i,sample_.value_second())*
244                target_(sample_.value_second()) *
245                (alpha_new2-alpha_old2) );
246      }
247           
248      epochs++; 
249      if (epochs>max_epochs_)
250        return false;
251    }
252   
253    trained_ = calculate_bias();
254    return trained_;
255  }
256
257
258  bool SVM::choose(const theplu::gslapi::vector& E)
259  {
260    //std::cout << "e choose\n"  ;
261    // First check for violation among SVs
262    // E should be the same for all SVs
263    // Choose that pair having largest violation/difference.
264    sample_.update_second(0);
265    sample_.update_first(0);
266    if (sample_.nof_sv()>1){
267      //std::cout << "there is SVs\n";
268      double max = E(sample_(0));
269      double min = max;
270      for (size_t i=1; i<sample_.nof_sv(); i++){ 
271        assert(alpha_(sample_(i))>tolerance_);
272        if (E(sample_(i)) > max){
273          max = E(sample_(i));
274          sample_.update_second(i);
275        }
276        else if (E(sample_(i))<min){
277          min = E(sample_(i));
278          sample_.update_first(i);
279        }
280      }
281      assert(alpha_(sample_.value_first())>tolerance_);
282      assert(alpha_(sample_.value_second())>tolerance_);
283
284     
285      if (E(sample_.value_second()) - E(sample_.value_first()) > 2*tolerance_){
286        return true;
287      }
288     
289      // If no violation check among non-support vectors
290      //std::cout << "no violation SVs\n";
291
292      sample_.shuffle();
293      //std::cout << "randomized\n";
294     
295      for (size_t i=sample_.nof_sv(); i<sample_.n();i++){
296        //std::cout << "nr: " << i << std::endl;
297        if (target_(sample_(i))==1){
298          if(E(sample_(i)) < E(sample_.value_first()) - 2*tolerance_){
299            sample_.update_second(i);
300            return true;
301          }
302        }
303        else{
304          if(E(sample_(i)) > E(sample_.value_second()) + 2*tolerance_){
305            sample_.update_first(i);
306            return true;
307          }
308        }
309      }
310    }
311
312    // if no support vectors - special case
313    else{
314      for (size_t i=0; i<sample_.n(); i++) {
315        if (target_(sample_(i))==1){
316          for (size_t j=0; j<sample_.n(); j++) {
317            if ( target_(sample_(j))==-1 && 
318                 E(sample_(i)) < E(sample_(j))+2*tolerance_ ){
319              sample_.update_first(i);
320              sample_.update_second(j);
321              return true;
322            }
323          }
324        }
325      }
326    }
327   
328    //std::cout << "Done!" << std::endl;
329    // If there is no violation then we should stop training
330    return false;
331
332  }
333 
334 
335  void SVM::bounds( double& u, double& v) const
336  {
337    if (target_(sample_.value_first())!=target_(sample_.value_second())) {
338      if (alpha_(sample_.value_second()) > alpha_(sample_.value_first())) {
339        v = std::numeric_limits<double>::max();
340        u = alpha_(sample_.value_second()) - alpha_(sample_.value_first());
341      }
342      else {
343        v = (std::numeric_limits<double>::max() - 
344             alpha_(sample_.value_first()) + 
345             alpha_(sample_.value_second()));
346        u = 0;
347      }
348    }
349    else {       
350      if (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) > 
351           std::numeric_limits<double>::max()) {
352        u = (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) - 
353              std::numeric_limits<double>::max());
354        v =  std::numeric_limits<double>::max();   
355      }
356      else {
357        u = 0;
358        v = alpha_(sample_.value_first()) + alpha_(sample_.value_second());
359      }
360    }
361  }
362 
363  bool SVM::calculate_bias(void)
364  {
365
366    // calculating output without bias
367    for (size_t i=0; i<output_.size(); i++) {
368      output_(i)=0;
369      for (size_t j=0; j<output_.size(); j++) 
370        output_(i)+=alpha_(j)*target_(j)*(*kernel_)(i,j);
371    }
372
373    if (!sample_.nof_sv()){
374      std::cerr << "SVM::train() error: " 
375                << "Cannot calculate bias because there is no support vector" 
376                << std::endl;
377      return false;
378    }
379
380    // For samples with alpha>0, we have: target*output=1-alpha/C
381    bias_=0;
382    for (size_t i=0; i<sample_.nof_sv(); i++) 
383      bias_+= ( target_(sample_(i)) * (1-alpha_(sample_(i))*C_inverse_) - 
384                output_(sample_(i)) );
385    bias_=bias_/sample_.nof_sv();
386    for (size_t i=0; i<output_.size(); i++) 
387      output_(i) += bias_;
388     
389    return true;
390  }
391
392}} // of namespace classifier and namespace theplu
Note: See TracBrowser for help on using the repository browser.