source: trunk/lib/svm/SVM.cc @ 307

Last change on this file since 307 was 307, checked in by Peter, 18 years ago

changed names of Kernels and made them work

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.4 KB
Line 
1// $Id: SVM.cc 307 2005-05-03 13:28:29Z peter $
2
3
4#include <c++_tools/svm/SVM.h>
5
6#include <c++_tools/svm/Kernel_MEV.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/utility/random_singleton.h>
11
12#include <cmath>
13#include <limits.h>
14#include <utility>
15#include <vector>
16
17
18namespace theplu {
19namespace svm { 
20
21  SVM::SVM(const Kernel_MEV& kernel, 
22           const gslapi::vector& target, 
23           const std::vector<size_t>& train_set) 
24  : alpha_(target.size(),0),
25    bias_(0),
26    c_(0),
27    kernel_(kernel),
28    max_epochs_(10000000),
29    target_(target),
30    trained_(false),
31    train_set_(train_set),
32    tolerance_(0.000001)
33       
34  {
35    if (max_epochs_>ULONG_MAX)
36      max_epochs_=ULONG_MAX;
37    if (!train_set_.size())
38      for (size_t i=0; i<target_.size(); i++)
39        train_set_.push_back(i); 
40  }
41
42
43
44  std::pair<size_t,size_t> 
45  SVM::choose(const theplu::gslapi::vector& E, 
46              const theplu::gslapi::vector& target_train,
47              const theplu::gslapi::vector& alpha_train,
48              bool& stop_condition)
49  {
50    std::pair<size_t, size_t> index(0,0);     
51   
52    // First check for violation among support vectors
53    double max = -10000000;
54    double min = 10000000;
55    // Peter, avoid looping over all alphas
56    for (size_t i=0; i<E.size(); i++){ 
57      if (alpha_train(i)){
58        if (E(i)>max){
59          max = E(i);
60          index.second = i;
61        }
62        if (E(i)<min){
63          min = E(i);
64          index.first = i;
65        }
66      }
67    }
68    //index = E.minmax_index(support_vectors);
69    if (E(index.second) - E(index.first) > 2*tolerance_)
70      return index;
71   
72    // If no violation check among non-support vectors
73    std::vector<size_t> tmp;
74    tmp.resize(E.size());
75    for (size_t i=0; i<E.size(); i++)
76      tmp[i]=i;
77    utility::my_uniform_rng a;
78    random_shuffle(tmp.begin(), tmp.end(), a);
79    for (size_t i=0; i<tmp.size(); i++){
80      if (target_train(tmp[i])==1){
81        if ( E(tmp[i]) < E(index.first) - 2*tolerance_){
82          index.second = tmp[i];
83          return index;
84        }
85      }
86      else if (target_train(tmp[i])==-1){
87        if ( E(tmp[i]) > E(index.second) + 2*tolerance_){
88          index.first = tmp[i];
89          return index;
90        }
91      }
92    }
93    // If there is no violation then we should stop
94    stop_condition = true;
95    return index;
96  }
97 
98 
99
100  bool SVM::train() 
101  {
102    using namespace std;
103    // initializing variables for optimization
104    gslapi::matrix kernel_train(train_set_.size(),train_set_.size());
105    gslapi::vector alpha_train(train_set_.size());
106    gslapi::vector target_train(train_set_.size());
107       
108    // copy submatrix to train on. Peter, might be unnecessary.
109    for (unsigned int i=0; i<train_set_.size(); i++) {
110      target_train(i) = target_(train_set_[i]);
111      for (unsigned int j=0; j<train_set_.size(); j++) 
112        kernel_train(i,j) = kernel_(train_set_[i],train_set_[j]); 
113    }
114    // Modify the diagonal of the kernel matrix (soft margin)
115    if (c_>0){
116      double d = 1/c_;
117      for (unsigned int i=0; i<kernel_train.columns(); i++)
118        kernel_train(i,i) += d; 
119    }
120 
121    gslapi::vector  E = ( kernel_train*target_train.mul(alpha_train) 
122                          - target_train ) ;
123
124    double upper_bound = pow(10., 32); //Peter, should be static...
125    unsigned long int epochs = 0;
126    double alpha_new2;
127    double alpha_new1;
128    double u;
129    double v;
130    bool stop_condition = false;
131    // Training loop
132    while(!stop_condition) {
133      pair<size_t, size_t> index = choose(E, target_train, 
134                                        alpha_train,stop_condition); 
135      // Defining upper and lower bounds for alpha
136      if (target_train[index.first]!=target_train[index.second]) {       
137        if (alpha_train[index.second] > alpha_train[index.first]) {
138          v = upper_bound;
139          u = alpha_train[index.second] - alpha_train[index.first];
140        }
141        else {
142          v = (upper_bound - alpha_train[index.first] + 
143               alpha_train[index.second]);
144          u = 0;
145        }
146      }
147      else {       
148        if (alpha_train[index.second] + alpha_train[index.first] > 
149            upper_bound) {
150          u = (alpha_train[index.second] + alpha_train[index.first] - 
151               upper_bound);
152          v = upper_bound;   
153        }
154        else {
155          u = 0;
156          v = alpha_train[index.first] + alpha_train[index.second];
157        }
158      }
159       
160      double k = ( kernel_train(index.first, index.first) + 
161                   kernel_train(index.second, index.second) - 
162                   2*kernel_train(index.first, index.second));
163 
164      alpha_new2 = (alpha_train[index.second] + target_train[index.second]*
165                    (E[index.first]-E[index.second])/k );
166   
167      if (alpha_new2 > v)
168        alpha_new2 = v;
169      else if (alpha_new2<u)
170        alpha_new2 = u;
171   
172      // Updating the alphas and some more
173      if (alpha_new2 < tolerance_)
174        alpha_new2=0;
175   
176      alpha_new1 = (alpha_train[index.first] + 
177                    (target_train[index.first] * target_train[index.second] * 
178                     (alpha_train[index.second] - alpha_new2)) );
179      if (alpha_new1 < tolerance_){
180        alpha_new1=0;
181      }
182
183      alpha_train[index.first] = alpha_new1;
184      alpha_train[index.second] = alpha_new2;
185   
186      E = kernel_train*target_train.mul(alpha_train) - target_train;
187      gslapi::vector ones(alpha_train.size(),1);
188   
189      epochs++; 
190      if (epochs>max_epochs_){
191        return false;
192      }
193    }
194
195    // Calculating the bias from support vectors. E = output - bias - target
196    statistics::Averager bias;
197    for (size_t i=0; i<E.size(); i++)
198      if (alpha_train(i))
199        bias.add(-E(i));
200    bias_ = bias.mean();
201 
202    trained_= true;
203    alpha_.set_all(0);
204    for (unsigned int i=0; i<train_set_.size(); i++) 
205      alpha_(train_set_[i]) = alpha_train(i);
206   
207    // training performed succesfully
208    return true;
209   
210  }
211
212 
213   
214   
215
216}} // of namespace svm and namespace theplu
Note: See TracBrowser for help on using the repository browser.