Changeset 323


Ignore:
Timestamp:
May 26, 2005, 10:54:30 PM (16 years ago)
Author:
Peter
Message:

major modifications in SVM::train() in an attempt to speed it up. Interface is changed, so from now if validation is needed it should be taken care of by the Kernel object.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/lib/svm/SVM.cc

    r307 r323  
    1010#include <c++_tools/utility/random_singleton.h>
    1111
     12#include <iostream>
     13#include <algorithm>
    1214#include <cmath>
    13 #include <limits.h>
     15#include <limits>
    1416#include <utility>
    1517#include <vector>
    … …  
    1921namespace svm { 
    2022
     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    utility::my_uniform_rng 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
    21139  SVM::SVM(const Kernel_MEV& kernel,
    22            const gslapi::vector& target,
    23            const std::vector<size_t>& train_set)
     140           const gslapi::vector& target)
     141           
    24142  : alpha_(target.size(),0),
    25143    bias_(0),
    26     c_(0),
     144    C_inverse_(0),
    27145    kernel_(kernel),
    28146    max_epochs_(10000000),
     147    output_(target.size(),0),
     148    sample_(target.size()),
    29149    target_(target),
    30150    trained_(false),
    31     train_set_(train_set),
    32     tolerance_(0.000001)
     151    tolerance_(0.00000001)
    33152       
    34153  {
    35154    if (max_epochs_>ULONG_MAX)
    36155      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;
     156  }
     157
     158
     159  bool SVM::train(void)
     160  {
    103161    // 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...
     162    assert(target_.size()==kernel_.size());
     163    assert(target_.size()==alpha_.size());
     164
     165    sample_.init(alpha_,tolerance_);
     166    gslapi::vector  E(target_.size());
     167    for (size_t i=0; i<E.size(); i++) {
     168      E(i)=0;
     169      for (size_t j=0; j<E.size(); j++)
     170        E(i) += kernel_mod(i,j)*target_(j)*alpha_(j);
     171      E(i)=E(i)-target_(i);
     172    }
     173    assert(target_.size()==E.size());
     174    assert(target_.size()==sample_.n());
     175
    125176    unsigned long int epochs = 0;
    126177    double alpha_new2;
    … …  
    128179    double u;
    129180    double v;
    130     bool stop_condition = false;
     181
    131182    // 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    
     183    while(choose(E)) {
     184      bounds(u,v);       
     185      double k = ( kernel_mod(sample_.value_first(), sample_.value_first()) +
     186                   kernel_mod(sample_.value_second(), sample_.value_second()) -
     187                   2*kernel_mod(sample_.value_first(), sample_.value_second()));
     188     
     189      double alpha_old1=alpha_(sample_.value_first());
     190      double alpha_old2=alpha_(sample_.value_second());
     191
     192      alpha_new2 = ( alpha_(sample_.value_second()) +
     193                     target_(sample_.value_second())*
     194                     ( E(sample_.value_first())-E(sample_.value_second()) )/k );
     195     
    167196      if (alpha_new2 > v)
    168197        alpha_new2 = v;
    169198      else if (alpha_new2<u)
    170199        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)) );
     200     
     201     
     202      // Updating the alphas
     203      // if alpha is 'zero' make the sample a non-support vector
     204      if (alpha_new2 < tolerance_){
     205        sample_.nsv_second();
     206      }
     207      else{
     208        sample_.sv_second();
     209      }
     210     
     211     
     212      alpha_new1 = (alpha_(sample_.value_first()) +
     213                    (target_(sample_.value_first()) *
     214                     target_(sample_.value_second()) *
     215                     (alpha_(sample_.value_second()) - alpha_new2) ));
     216           
     217      // if alpha is 'zero' make the sample a non-support vector
    179218      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    
     219        sample_.nsv_first();
     220      }
     221      else
     222        sample_.sv_first();
     223     
     224      alpha_(sample_.value_first()) = alpha_new1;
     225      alpha_(sample_.value_second()) = alpha_new2;
     226     
     227      // update E vector
     228      // Peter, perhaps one should only update SVs, but what happens in choose?
     229      for (size_t i=0; i<E.size(); i++) {
     230        E(i)+=( kernel_mod(i,sample_.value_first())*
     231                target_(sample_.value_first()) *
     232                (alpha_new1-alpha_old1) );
     233        E(i)+=( kernel_mod(i,sample_.value_second())*
     234                target_(sample_.value_second()) *
     235                (alpha_new2-alpha_old2) );
     236      }
     237           
    189238      epochs++; 
    190       if (epochs>max_epochs_){
     239      if (epochs>max_epochs_)
    191240        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 
     241    }
     242   
     243    trained_ = calculate_bias();
     244    return trained_;
     245  }
     246
     247
     248  bool SVM::choose(const theplu::gslapi::vector& E)
     249  {
     250    //std::cout << "e choose\n"  ;
     251    // First check for violation among SVs
     252    // E should be the same for all SVs
     253    // Choose that pair having largest violation/difference.
     254    sample_.update_second(0);
     255    sample_.update_first(0);
     256    if (sample_.nof_sv()>1){
     257      //std::cout << "there is SVs\n";
     258      double max = E(sample_(0));
     259      double min = max;
     260      for (size_t i=1; i<sample_.nof_sv(); i++){
     261        assert(alpha_(sample_(i))>tolerance_);
     262        if (E(sample_(i)) > max){
     263          max = E(sample_(i));
     264          sample_.update_second(i);
     265        }
     266        else if (E(sample_(i))<min){
     267          min = E(sample_(i));
     268          sample_.update_first(i);
     269        }
     270      }
     271      assert(alpha_(sample_.value_first())>tolerance_);
     272      assert(alpha_(sample_.value_second())>tolerance_);
     273
     274     
     275      if (E(sample_.value_second()) - E(sample_.value_first()) > 2*tolerance_){
     276        return true;
     277      }
     278     
     279      // If no violation check among non-support vectors
     280      //std::cout << "no violation SVs\n";
     281
     282      sample_.shuffle();
     283      //std::cout << "randomized\n";
     284     
     285      for (size_t i=sample_.nof_sv(); i<sample_.n();i++){
     286        //std::cout << "nr: " << i << std::endl;
     287        if (target_(sample_(i))==1){
     288          if(E(sample_(i)) < E(sample_.value_first()) - 2*tolerance_){
     289            sample_.update_second(i);
     290            return true;
     291          }
     292        }
     293        else{
     294          if(E(sample_(i)) > E(sample_.value_second()) + 2*tolerance_){
     295            sample_.update_first(i);
     296            return true;
     297          }
     298        }
     299      }
     300    }
     301
     302    // if no support vectors - special case
     303    else{
     304      for (size_t i=0; i<sample_.n(); i++) {
     305        if (target_(sample_(i))==1){
     306          for (size_t j=0; j<sample_.n(); j++) {
     307            if ( target_(sample_(j))==-1 &&
     308                 E(sample_(i)) < E(sample_(j))+2*tolerance_ ){
     309              sample_.update_first(i);
     310              sample_.update_second(j);
     311              return true;
     312            }
     313          }
     314        }
     315      }
     316    }
     317   
     318    //std::cout << "Done!" << std::endl;
     319    // If there is no violation then we should stop training
     320    return false;
     321
     322  }
    212323 
    213    
    214    
     324 
     325  void SVM::bounds( double& u, double& v) const
     326  {
     327    if (target_(sample_.value_first())!=target_(sample_.value_second())) {
     328      if (alpha_(sample_.value_second()) > alpha_(sample_.value_first())) {
     329        v = std::numeric_limits<double>::max();
     330        u = alpha_(sample_.value_second()) - alpha_(sample_.value_first());
     331      }
     332      else {
     333        v = (std::numeric_limits<double>::max() -
     334             alpha_(sample_.value_first()) +
     335             alpha_(sample_.value_second()));
     336        u = 0;
     337      }
     338    }
     339    else {       
     340      if (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) >
     341           std::numeric_limits<double>::max()) {
     342        u = (alpha_(sample_.value_second()) + alpha_(sample_.value_first()) -
     343              std::numeric_limits<double>::max());
     344        v =  std::numeric_limits<double>::max();   
     345      }
     346      else {
     347        u = 0;
     348        v = alpha_(sample_.value_first()) + alpha_(sample_.value_second());
     349      }
     350    }
     351  }
     352 
     353  bool SVM::calculate_bias(void)
     354  {
     355
     356    // calculating output without bias
     357    for (size_t i=0; i<output_.size(); i++) {
     358      output_(i)=0;
     359      for (size_t j=0; j<output_.size(); j++)
     360        output_(i)+=alpha_(j)*target_(j)*kernel_(i,j);
     361    }
     362
     363    if (!sample_.nof_sv()){
     364      std::cerr << "SVM::train() error: "
     365                << "Cannot calculate bias because there is no support vector"
     366                << std::endl;
     367      return false;
     368    }
     369
     370    // For samples with alpha>0, we have: target*output=1-alpha/C
     371    bias_=0;
     372    for (size_t i=0; i<sample_.nof_sv(); i++)
     373      bias_+= ( target_(sample_(i)) * (1-alpha_(sample_(i))*C_inverse_) -
     374                output_(sample_(i)) );
     375    bias_=bias_/sample_.nof_sv();
     376    for (size_t i=0; i<output_.size(); i++)
     377      output_(i) += bias_;
     378     
     379    return true;
     380  }
    215381
    216382}} // of namespace svm and namespace theplu

  • trunk/lib/svm/SVM.h

    r318 r323  
    1313namespace theplu {
    1414namespace svm { 
    15   ///
     15
     16  // @internal Class keeping track of which samples are support vectors and
     17  // not. The first nof_sv elements in the vector are indices of the
     18  // support vectors
     19  //
     20  class Index
     21  {
     22
     23  public:
     24    //Default Contructor
     25    Index();
     26
     27    //
     28    Index(const size_t);
     29
     30    //not implemented
     31    Index(const Index&);
     32
     33    // @return index_first
     34    inline size_t index_first(void) const
     35    { assert(index_first_<n()); return index_first_; }
     36
     37    // @return index_second
     38    inline size_t index_second(void) const { assert(index_second_<n());return index_second_; }
     39
     40    // synch the object against alpha
     41    void init(const gslapi::vector& alpha, const double);
     42
     43    // @return nof samples
     44    inline size_t n(void) const { return vec_.size(); }
     45
     46    // @return nof support vectors
     47    inline size_t nof_sv(void) const { return nof_sv_; }
     48
     49    // making first to an nsv. If already sv, nothing happens.
     50    void nsv_first(void);
     51
     52    // making first to an nsv. If already sv, nothing happens.
     53    void nsv_second(void);
     54
     55    // randomizes the nsv part of vector and sets index_first to
     56    // nof_sv_ (the first nsv)
     57    void shuffle(void);
     58
     59    // making first to a sv. If already sv, nothing happens.
     60    void sv_first(void);
     61
     62    // making first to a sv. If already sv, nothing happens.
     63    void sv_second(void);
     64
     65    //
     66    void update_first(const size_t);
     67
     68    //
     69    void update_second(const size_t);
     70
     71    // @return value_first
     72    inline size_t value_first(void) const { assert(value_first_<n()); return value_first_; }
     73
     74    // @return const ref value_second
     75    inline size_t value_second(void) const {  assert(value_first_<n()); return value_second_; }
     76
     77    inline size_t operator()(size_t i) const {
     78      assert(i<n()); assert(vec_[i]<n()); return vec_[i]; }
     79
     80  private:
     81    size_t index_first_;
     82    size_t index_second_;
     83    size_t nof_sv_;
     84    std::vector<size_t> vec_;
     85    size_t value_first_; // vec_[index_first_] exists for fast access
     86    size_t value_second_; // vec_[index_second_] exists for fast access
     87   
     88  };
     89
     90  ///
    1691  /// Class for SVM using Keerthi's second modification of Platt's SMO.
    17   /// @todo speed it up.
    1892  ///   
    1993  class SVM
    … …  
    2195 
    2296  public:
    23     ///
    24     /// Constructor taking the kernel matrix and the target vector as input
     97    ///
     98    /// Default constructor (not implemented)
     99    ///
     100    SVM();
     101
     102    ///
     103    /// Constructor taking the kernel and the target vector as input.
    25104    ///
    26105    SVM(const Kernel_MEV&,
    27         const gslapi::vector&,
    28         const std::vector<size_t>& = std::vector<size_t>());
     106        const gslapi::vector&);
     107
     108    ///
     109    /// @todo Copy constructor.
     110    ///
     111    SVM(const SVM&);
    29112         
    30113    ///
    31     /// Function returns \f$\alpha\f$
    32     ///
    33     inline gslapi::vector get_alpha(void) const { return alpha_; }
    34 
    35     ///
    36     /// Function returns the C-parameter
    37     ///
    38     inline double get_c(void) const { return c_; }
    39 
    40     ///
     114    /// The C-parameter is the balance term (see train()). A very
     115    /// large C means the training will be focused on getting samples
     116    /// correctly classified, with risk for overfitting and poor
     117    /// generalisation. A too small C will result in a training where
     118    /// misclassifications are not penalized. C is weighted with
     119    /// respect to the size, so \f$ n_+C_+ = n_-C_- \f$, meaning a
     120    /// misclassificaion of the smaller group is penalized
     121    /// harder. This balance is equivalent to the one occuring for
     122    /// regression with regularisation, or ANN-training with a
     123    /// weight-decay term. Default is C set to infinity.
     124    ///
     125    /// @returns mean of vector \f$ C_i \f$
     126    ///
     127    inline double C(void) const { return 1/C_inverse_; }
     128
     129    ///
     130    /// @return \f$\alpha\f$
     131    ///
     132    inline const gslapi::vector& alpha(void) const { return alpha_; }
     133
     134    ///
     135    /// Default is max_epochs set to 10,000,000.
     136    ///
    41137    /// @return number of maximal epochs
    42138    ///
    43139    inline long int max_epochs(void) const {return max_epochs_;}
    44140   
    45     ///
    46     /// Changing number of maximal epochs
    47     ///
    48     inline void max_epochs(const unsigned long int d) {max_epochs_=d;}
    49    
    50     ///
     141    ///
     142    /// The output is calculated as \f$ o_i = \sum \alpha_j t_j K_{ij}
     143    /// + bias \f$, where \f$ t \f$ is the target.
     144    ///
    51145    /// @return output
    52     /// @todo
    53     theplu::gslapi::vector output(void) const;
    54 
    55     ///
    56     /// Changing the C-parameter
    57     ///
    58     inline void set_c(const double c) {c_ = c;}
     146    ///
     147    inline const theplu::gslapi::vector& output(void) const { return output_; }
     148
     149    ///
     150    /// @todo Function sets \f$ \alpha=0 \f$ and makes SVM untrained.
     151    ///
     152    void reset(void);
    59153
    60154    ///
    61155    /// Training the SVM following Platt's SMO, with Keerti's
    62     /// modifacation. However the complete kernel is stored in
    63     /// memory. The reason for this is speed. When number of samples N
    64     /// is large this is not possible since the memory cost for the
    65     /// kernel scales N^2. In that case one should follow the SMO and
    66     /// calculate the kernel elements sequentially. Minimizing \f$
    67     /// \frac{1}{2}\sum
     156    /// modifacation. Minimizing \f$ \frac{1}{2}\sum
    68157    /// y_iy_j\alpha_i\alpha_j(K_{ij}+\frac{1}{C_i}\delta_{ij}) \f$,
    69     /// which corresponds to minimizing \f$ \sum w_i^2+\sum
    70     /// C_i\xi_i^2 \f$
    71     ///
    72 
    73     bool train(void);
    74    
     158    /// which corresponds to minimizing \f$ \sum w_i^2+\sum C_i\xi_i^2
     159    /// \f$. train_set defines which samples to use for the
     160    /// training. Default is to use all samples.
     161    ///
     162    bool train();
     163
     164       
    75165     
    76166  private:
    77167    gslapi::vector alpha_;
    78168    double bias_;
    79     double c_;
    80     const Kernel_MEV kernel_;
     169    double C_inverse_;
     170    const Kernel_MEV& kernel_;
    81171    unsigned long int max_epochs_;
    82     const gslapi::vector& target_;
     172    gslapi::vector output_;
     173    Index sample_;
     174    const gslapi::vector& target_;
    83175    bool trained_;
    84     std::vector<size_t> train_set_;
    85176    double tolerance_;
    86177   
     178
     179    ///
     180    /// Calculates bounds for alpha2
     181    ///
     182    void bounds(double&, double&) const;
     183
     184    ///
     185    /// @todo Function calculating bias
     186    ///
     187    /// @return true if successful
     188    ///
     189    bool calculate_bias(void);
     190
    87191    ///
    88192    ///   Private function choosing which two elements that should be
    89193    ///   updated. First checking for the biggest violation (output - target =
    90194    ///   0) among support vectors (alpha!=0). If no violation was found check
    91     ///   for sequentially among the other samples. If no violation there as
    92     ///   well, stop_condition is fullfilled.
    93     ///
    94     std::pair<size_t, size_t> choose(const theplu::gslapi::vector&,
    95                                      const theplu::gslapi::vector&,
    96                                      const theplu::gslapi::vector&,
    97                                      bool&);
    98 
    99        
     195    ///   sequentially among the other samples. If no violation there as
     196    ///   well training is completed
     197    ///
     198    ///  @return true if a pair of samples that violate the conditions
     199    ///  can be found
     200    ///
     201    bool choose(const theplu::gslapi::vector&);
     202
     203    ///
     204    /// @return kernel modified in with diagonal term (soft margin)
     205    ///
     206    inline double kernel_mod(const size_t i, const size_t j) const
     207    { return i!=j ? kernel_(i,j) : kernel_(i,j) + C_inverse_; }
    100208   
    101209  };
    102210
    103 
    104 
    105 
     211  ///
     212  /// @todo The output operator for the SVM class.
     213  ///
     214  std::ostream& operator<< (std::ostream& s, const SVM&);
     215 
     216 
    106217}} // of namespace svm and namespace theplu
    107218

  • trunk/test/svm_test.cc

    r307 r323  
    1414using namespace theplu;
    1515
    16 int main()
    17 
     16int main(const int argc,const char* argv[])
    1817{ 
    1918
    20 //  ifstream is("data/nm_kernel.txt");
    21 //  theplu::gslapi::matrix kernel_matrix(is);
    22 //  is.close();
    23  
    24 //  is.open("data/nm_target_bin.txt");
    25 //  theplu::gslapi::vector target(is);
    26 //  is.close();
     19  bool print = (argc>1 && argv[1]==std::string("-p"));
     20  bool ok = true;
    2721
    28 //  is.open("data/nm_alpha_linear_matlab.txt");
    29 //  theplu::gslapi::vector alpha_matlab(is);
    30 //  is.close();
     22  gslapi::matrix data2(2,3);
     23  data2(0,0)=0;
     24  data2(1,0)=0;
     25  data2(0,1)=0;
     26  data2(1,1)=1;
     27  data2(0,2)=1;
     28  data2(1,2)=0;
     29  gslapi::vector target2(3);
     30  target2(0)=-1;
     31  target2(1)=1;
     32  target2(2)=1;
     33  svm::KernelFunction* kf2 = new svm::PolynomialKernelFunction();
     34  svm::Kernel_MEV kernel2(data2,*kf2);
     35  assert(kernel2.size()==3);
     36  svm::SVM svm2(kernel2, target2);
     37  svm2.train();
    3138
    32 //  theplu::cpptools::SVM svm(kernel_matrix, target);
    33 //    svm.train();
     39  if (svm2.alpha()*target2){
     40    std::cerr << "condition not fullfilled" << std::endl;
     41    return -1;
     42  }
    3443
    35 //  theplu::gslapi::vector alpha = svm.get_alpha();
     44  if (svm2.alpha()(1)!=2 || svm2.alpha()(2)!=2){
     45    std::cerr << "wrong alpha" << std::endl;
     46    std::cerr << "alpha: " << svm2.alpha() <<  std::endl;
     47    std::cerr << "expected: 4 2 2" <<  std::endl;
     48
     49    return -1;
     50  }
     51
     52 
     53
     54  std::ifstream is("data/nm_data_centralized.txt");
     55  gslapi::matrix transposed_data(is);
     56  is.close();
     57  // Because how the kernel is treated is changed, data must be transposed.
     58  gslapi::matrix data=transposed_data;
     59
     60  svm::KernelFunction* kf = new svm::PolynomialKernelFunction();
     61  svm::Kernel_SEV kernel(data,*kf);
     62
     63
     64  is.open("data/nm_target_bin.txt");
     65  theplu::gslapi::vector target(is);
     66  is.close();
     67
     68  is.open("data/nm_alpha_linear_matlab.txt");
     69  theplu::gslapi::vector alpha_matlab(is);
     70  is.close();
     71
     72  theplu::svm::SVM svm(kernel, target);
     73  if (!svm.train()){
     74    ok=false;
     75    if (print)
     76      std::cerr << "Training failured" << std::endl;
     77  }
     78
     79  theplu::gslapi::vector alpha = svm.alpha();
    3680     
    37 //  // Comparing alpha to alpha_matlab
    38 //  theplu::gslapi::vector diff_alpha = alpha - alpha_matlab;
    39 //  if (diff_alpha*diff_alpha> pow(10.0,-10)){
    40 //    cerr << "Difference to matlab alphas too large/n";
    41 //    return -1;
    42 //  }
    43 //  // Comparing output to target
    44 //  theplu::gslapi::vector output = svm.get_output();
    45 //  double slack = 0;
    46 //  for (unsigned int i=0; i<target.size(); i++){
    47 //    if (output[i]*target[i] < 1){
    48 //      slack += 1 - output[i]*target[i];
    49 //    }
    50 //  }
    51 //  if (slack > pow(10.0,-10)){
    52 //    cerr << "Difference to matlab alphas too large/n";
    53 //    return -1;
    54 //  }
     81  // Comparing alpha to alpha_matlab
     82  theplu::gslapi::vector diff_alpha = alpha - alpha_matlab;
     83  if (diff_alpha*diff_alpha> 1e-10 ){
     84    if (print)
     85      std::cerr << "Difference to matlab alphas too large\n";
     86    ok=false;
     87  }
     88
     89  // Comparing output to target
     90  theplu::gslapi::vector output = svm.output();
     91  double slack = 0;
     92  for (unsigned int i=0; i<target.size(); i++){
     93    if (output[i]*target[i] < 1){
     94      slack += 1 - output[i]*target[i];
     95    }
     96  }
     97  if (slack > 1e-7){
     98    if (print){
     99      std::cerr << "Slack too large. Is the bias correct?\n";
     100      std::cerr << "slack: " << slack << std::endl;
     101    }
     102    ok = false;
     103  }
    55104 
    56   // testing on non solvable problem
    57   gslapi::matrix data(3,1);
    58   data(0,0)=-1;
    59   data(1,0)=10;
    60   data(2,0)=0;
    61    
    62   gslapi::vector target(3,1.0);
    63   target(2)=-1;
    64   svm::KernelFunction* kf =
    65     new svm::PolynomialKernelFunction();
    66   svm::Kernel_SEV kernel(data,*kf);
    67   //theplu::gslapi::matrix m=kernel.get();
    68   //std::cout << "kernel:\n" << m << "\n";
    69   svm::SVM svm(kernel,target);
    70   //  svm.set_c(1);
    71   //svm.train();
    72105 
    73   return 0;
     106  if(ok)
     107    return 0;
     108  return -1;
    74109 
    75110}
Note: See TracChangeset for help on using the changeset viewer.