source: trunk/yat/classifier/SVM.cc @ 722

Last change on this file since 722 was 722, checked in by Markus Ringnér, 15 years ago

Closes #129

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date ID
File size: 10.6 KB
Line 
1// $Id$
2
3/*
4  Copyright (C) The authors contributing to this file.
5
6  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 2 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
21  02111-1307, USA.
22*/
23
24#include "SVM.h"
25#include "DataLookup2D.h"
26#include "yat/random/random.h"
27#include "yat/statistics/Averager.h"
28#include "yat/utility/matrix.h"
29#include "yat/utility/vector.h"
30
31#include <algorithm>
32#include <cassert>
33#include <cctype>
34#include <cmath>
35#include <limits>
36#include <utility>
37#include <vector>
38
39namespace theplu {
40namespace yat {
41namespace classifier { 
42
43  SVM::SVM(const KernelLookup& kernel, const Target& target)
44    : SupervisedClassifier(target),
45      alpha_(target.size(),0),
46      bias_(0),
47      C_inverse_(0),
48      kernel_(&kernel),
49      margin_(0),
50      max_epochs_(100000),
51      output_(target.size(),0),
52      owner_(false),
53      sample_(target.size()),
54      trained_(false),
55      tolerance_(0.00000001)
56  {
57#ifndef NDEBUG
58    assert(kernel.columns()==kernel.rows());
59    assert(kernel.columns()==alpha_.size());
60    for (size_t i=0; i<alpha_.size(); i++) 
61      for (size_t j=0; j<alpha_.size(); j++)
62        assert(kernel(i,j)==kernel(j,i));
63    for (size_t i=0; i<alpha_.size(); i++) 
64      for (size_t j=0; j<alpha_.size(); j++) 
65        assert((*kernel_)(i,j)==(*kernel_)(j,i));
66    for (size_t i = 0; i<kernel_->rows(); i++)
67      for (size_t j = 0; j<kernel_->columns(); j++)
68        if (std::isnan((*kernel_)(i,j)))
69          std::cerr << "SVM: Found nan in kernel: " << i << " " 
70                    << j << std::endl;
71#endif       
72  }
73
74  SVM::~SVM()
75  {
76    if (owner_)
77      delete kernel_;
78  }
79
80  const utility::vector& SVM::alpha(void) const
81  {
82    return alpha_;
83  }
84
85  double SVM::C(void) const
86  {
87    return 1.0/C_inverse_;
88  }
89
90  void SVM::calculate_margin(void)
91  {
92    margin_ = 0;
93    for(size_t i = 0; i<alpha_.size(); ++i){
94      margin_ += alpha_(i)*target(i)*kernel_mod(i,i)*alpha_(i)*target(i);
95      for(size_t j = i+1; j<alpha_.size(); ++j)
96        margin_ += 2*alpha_(i)*target(i)*kernel_mod(i,j)*alpha_(j)*target(j);
97    }
98  }
99
100  const DataLookup2D& SVM::data(void) const
101  {
102    return *kernel_;
103  }
104
105
106  double SVM::kernel_mod(const size_t i, const size_t j) const
107  {
108    return i!=j ? (*kernel_)(i,j) : (*kernel_)(i,j) + C_inverse_;
109  }
110
111  SupervisedClassifier* SVM::make_classifier(const DataLookup2D& data,
112                                             const Target& target) const
113  {
114    SVM* sc=0;
115    try {
116      const KernelLookup& kernel = dynamic_cast<const KernelLookup&>(data);
117      assert(data.rows()==data.columns());
118      assert(data.columns()==target.size());
119      sc = new SVM(kernel,target);
120
121    //Copy those variables possible to modify from outside
122    // Peter, in particular C
123    }
124    catch (std::bad_cast) {
125      std::cerr << "Warning: SVM::make_classifier only takes KernelLookup" 
126                << std::endl;
127    }
128    return sc;
129  }
130
131  long int SVM::max_epochs(void) const
132  {
133    return max_epochs_;
134  }
135
136  const theplu::yat::utility::vector& SVM::output(void) const
137  {
138    return output_;
139  }
140
141  void SVM::predict(const DataLookup2D& input, utility::matrix& prediction) const
142  {
143    // Peter, should check success of dynamic_cast
144    const KernelLookup& input_kernel = dynamic_cast<const KernelLookup&>(input);
145
146    assert(input.rows()==alpha_.size());
147    prediction = utility::matrix(2,input.columns(),0);
148    for (size_t i = 0; i<input.columns(); i++){
149      for (size_t j = 0; j<input.rows(); j++){
150        prediction(0,i) += target(j)*alpha_(j)*input_kernel(j,i);
151        assert(target(j));
152      }
153      prediction(0,i) = margin_ * (prediction(0,i) + bias_);
154    }
155
156    for (size_t i = 0; i<prediction.columns(); i++)
157      prediction(1,i) = -prediction(0,i);
158   
159    assert(prediction(0,0));
160  }
161
162  double SVM::predict(const DataLookup1D& x) const
163  {
164    double y=0;
165    for (size_t i=0; i<alpha_.size(); i++)
166      y += alpha_(i)*target_(i)*kernel_->element(x,i);
167
168    return margin_*(y+bias_);
169  }
170
171  double SVM::predict(const DataLookupWeighted1D& x) const
172  {
173    double y=0;
174    for (size_t i=0; i<alpha_.size(); i++)
175      y += alpha_(i)*target_(i)*kernel_->element(x,i);
176
177    return margin_*(y+bias_);
178  }
179
180  void SVM::reset(void)
181  {
182    trained_=false;
183    alpha_=utility::vector(target_.size(), 0);
184  }
185
186  int SVM::target(size_t i) const
187  {
188    return target_.binary(i) ? 1 : -1;
189  }
190
191  bool SVM::train(void) 
192  {
193    // initializing variables for optimization
194    assert(target_.size()==kernel_->rows());
195    assert(target_.size()==alpha_.size());
196
197    sample_.init(alpha_,tolerance_);
198    utility::vector   E(target_.size(),0);
199    for (size_t i=0; i<E.size(); i++) {
200      for (size_t j=0; j<E.size(); j++) 
201        E(i) += kernel_mod(i,j)*target(j)*alpha_(j);
202      E(i)-=target(i);
203    }
204    assert(target_.size()==E.size());
205    assert(target_.size()==sample_.size());
206
207    unsigned long int epochs = 0;
208    double alpha_new2;
209    double alpha_new1;
210    double u;
211    double v;
212
213    // Training loop
214    while(choose(E)) {
215      bounds(u,v);       
216      double k = ( kernel_mod(sample_.value_first(), sample_.value_first()) + 
217                   kernel_mod(sample_.value_second(), sample_.value_second()) - 
218                   2*kernel_mod(sample_.value_first(), sample_.value_second()));
219     
220      double alpha_old1=alpha_(sample_.value_first());
221      double alpha_old2=alpha_(sample_.value_second());
222      alpha_new2 = ( alpha_(sample_.value_second()) + 
223                     target(sample_.value_second())*
224                     ( E(sample_.value_first())-E(sample_.value_second()) )/k );
225     
226      if (alpha_new2 > v)
227        alpha_new2 = v;
228      else if (alpha_new2<u)
229        alpha_new2 = u;
230     
231      // Updating the alphas
232      // if alpha is 'zero' make the sample a non-support vector
233      if (alpha_new2 < tolerance_){
234        sample_.nsv_second();
235      }
236      else{
237        sample_.sv_second();
238      }
239     
240     
241      alpha_new1 = (alpha_(sample_.value_first()) + 
242                    (target(sample_.value_first()) * 
243                     target(sample_.value_second()) * 
244                     (alpha_(sample_.value_second()) - alpha_new2) ));
245           
246      // if alpha is 'zero' make the sample a non-support vector
247      if (alpha_new1 < tolerance_){
248        sample_.nsv_first();
249      }
250      else
251        sample_.sv_first();
252     
253      alpha_(sample_.value_first()) = alpha_new1;
254      alpha_(sample_.value_second()) = alpha_new2;
255     
256      // update E vector
257      // Peter, perhaps one should only update SVs, but what happens in choose?
258      for (size_t i=0; i<E.size(); i++) {
259        E(i)+=( kernel_mod(i,sample_.value_first())*
260                target(sample_.value_first()) *
261                (alpha_new1-alpha_old1) );
262        E(i)+=( kernel_mod(i,sample_.value_second())*
263                target(sample_.value_second()) *
264                (alpha_new2-alpha_old2) );
265      }
266           
267      epochs++; 
268      if (epochs>max_epochs_){
269        std::cerr << "WARNING: SVM: maximal number of epochs reached.\n";
270        calculate_bias();
271        calculate_margin();
272        return false;
273      }
274    }
275    calculate_margin();
276    trained_ = calculate_bias();
277    return trained_;
278  }
279
280
281  bool SVM::choose(const theplu::yat::utility::vector& E)
282  {
283    // First check for violation among SVs
284    // E should be the same for all SVs
285    // Choose that pair having largest violation/difference.
286    sample_.update_second(0);
287    sample_.update_first(0);
288    if (sample_.nof_sv()>1){
289
290      double max = E(sample_(0));
291      double min = max;
292      for (size_t i=1; i<sample_.nof_sv(); i++){ 
293        assert(alpha_(sample_(i))>tolerance_);
294        if (E(sample_(i)) > max){
295          max = E(sample_(i));
296          sample_.update_second(i);
297        }
298        else if (E(sample_(i))<min){
299          min = E(sample_(i));
300          sample_.update_first(i);
301        }
302      }
303      assert(alpha_(sample_.value_first())>tolerance_);
304      assert(alpha_(sample_.value_second())>tolerance_);
305
306      if (E(sample_.value_second()) - E(sample_.value_first()) > 2*tolerance_){
307        return true;
308      }
309     
310      // If no violation check among non-support vectors
311      sample_.shuffle();
312      for (size_t i=sample_.nof_sv(); i<sample_.size();i++){
313        if (target_.binary(sample_(i))){
314          if(E(sample_(i)) < E(sample_.value_first()) - 2*tolerance_){
315            sample_.update_second(i);
316            return true;
317          }
318        }
319        else{
320          if(E(sample_(i)) > E(sample_.value_second()) + 2*tolerance_){
321            sample_.update_first(i);
322            return true;
323          }
324        }
325      }
326    }
327
328    // if no support vectors - special case
329    else{
330      // to avoid getting stuck we shuffle
331      sample_.shuffle();
332      for (size_t i=0; i<sample_.size(); i++) {
333        if (target(sample_(i))==1){
334          for (size_t j=0; j<sample_.size(); j++) {
335            if ( target(sample_(j))==-1 && 
336                 E(sample_(i)) < E(sample_(j))+2*tolerance_ ){
337              sample_.update_first(i);
338              sample_.update_second(j);
339              return true;
340            }
341          }
342        }
343      }
344    }
345   
346    // If there is no violation then we should stop training
347    return false;
348
349  }
350 
351 
352  void SVM::bounds( double& u, double& v) const
353  {
354    if (target(sample_.value_first())!=target(sample_.value_second())) {
355      if (alpha_(sample_.value_second()) > alpha_(sample_.value_first())) {
356        v = std::numeric_limits<double>::max();
357        u = alpha_(sample_.value_second()) - alpha_(sample_.value_first());
358      }
359      else {
360        v = (std::numeric_limits<double>::max() - 
361             alpha_(sample_.value_first()) + 
362             alpha_(sample_.value_second()));
363        u = 0;
364      }
365    }
366    else {       
367      if (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) > 
368           std::numeric_limits<double>::max()) {
369        u = (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) - 
370              std::numeric_limits<double>::max());
371        v =  std::numeric_limits<double>::max();   
372      }
373      else {
374        u = 0;
375        v = alpha_(sample_.value_first()) + alpha_(sample_.value_second());
376      }
377    }
378  }
379 
380  bool SVM::calculate_bias(void)
381  {
382
383    // calculating output without bias
384    for (size_t i=0; i<output_.size(); i++) {
385      output_(i)=0;
386      for (size_t j=0; j<output_.size(); j++) 
387        output_(i)+=alpha_(j)*target(j) * (*kernel_)(i,j);
388    }
389
390    if (!sample_.nof_sv()){
391      std::cerr << "SVM::train() error: " 
392                << "Cannot calculate bias because there is no support vector" 
393                << std::endl;
394      return false;
395    }
396
397    // For samples with alpha>0, we have: target*output=1-alpha/C
398    bias_=0;
399    for (size_t i=0; i<sample_.nof_sv(); i++) 
400      bias_+= ( target(sample_(i)) * (1-alpha_(sample_(i))*C_inverse_) - 
401                output_(sample_(i)) );
402    bias_=bias_/sample_.nof_sv();
403    for (size_t i=0; i<output_.size(); i++) 
404      output_(i) += bias_;
405     
406    return true;
407  }
408
409  void SVM::set_C(const double C)
410  {
411    C_inverse_ = 1/C;
412  }
413
414
415}}} // of namespace classifier, yat, and theplu
Note: See TracBrowser for help on using the repository browser.