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

Last change on this file since 959 was 959, checked in by Peter, 14 years ago

Fixed so NBC and SVM are throwing when unexpected DataLookup2D is
apssed to make_classifier or predict.

Speeding up NBC::predict by separating weighted code from
unweighted. Also fixed some bugs in NBC.

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