SVMOcas.cpp

Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Written (W) 2007-2008 Vojtech Franc
00008  * Written (W) 2007-2009 Soeren Sonnenburg
00009  * Copyright (C) 2007-2009 Fraunhofer Institute FIRST and Max-Planck-Society
00010  */
00011 
00012 #include "features/Labels.h"
00013 #include "lib/Mathematics.h"
00014 #include "lib/Time.h"
00015 #include "base/Parallel.h"
00016 #include "classifier/LinearClassifier.h"
00017 #include "classifier/svm/SVMOcas.h"
00018 #include "classifier/svm/libocas.h"
00019 #include "features/DotFeatures.h"
00020 #include "features/Labels.h"
00021 
00022 CSVMOcas::CSVMOcas(E_SVM_TYPE type)
00023 : CLinearClassifier(), use_bias(true), bufsize(3000), C1(1), C2(1),
00024     epsilon(1e-3), method(type)
00025 {
00026     w=NULL;
00027     old_w=NULL;
00028 }
00029 
00030 CSVMOcas::CSVMOcas(
00031     float64_t C, CDotFeatures* traindat, CLabels* trainlab)
00032 : CLinearClassifier(), use_bias(true), bufsize(3000), C1(C), C2(C),
00033     epsilon(1e-3)
00034 {
00035     w=NULL;
00036     old_w=NULL;
00037     method=SVM_OCAS;
00038     set_features(traindat);
00039     set_labels(trainlab);
00040 }
00041 
00042 
00043 CSVMOcas::~CSVMOcas()
00044 {
00045 }
00046 
00047 bool CSVMOcas::train()
00048 {
00049     SG_INFO("C=%f, epsilon=%f, bufsize=%d\n", get_C1(), get_epsilon(), bufsize);
00050     SG_DEBUG("use_bias = %i\n", get_bias_enabled()) ;
00051 
00052     ASSERT(labels);
00053     ASSERT(features);
00054     ASSERT(labels->is_two_class_labeling());
00055 
00056     int32_t num_train_labels=0;
00057     lab=labels->get_labels(num_train_labels);
00058     w_dim=features->get_dim_feature_space();
00059     int32_t num_vec=features->get_num_vectors();
00060 
00061     ASSERT(num_vec==num_train_labels);
00062     ASSERT(num_vec>0);
00063 
00064     delete[] w;
00065     w=new float64_t[w_dim];
00066     memset(w, 0, w_dim*sizeof(float64_t));
00067 
00068     delete[] old_w;
00069     old_w=new float64_t[w_dim];
00070     memset(old_w, 0, w_dim*sizeof(float64_t));
00071     bias=0;
00072     old_bias=0;
00073 
00074     tmp_a_buf=new float64_t[w_dim];
00075     cp_value=new float64_t*[bufsize];
00076     cp_index=new uint32_t*[bufsize];
00077     cp_nz_dims=new uint32_t[bufsize];
00078     cp_bias=new float64_t[bufsize];
00079     memset(cp_bias, 0, sizeof(float64_t)*bufsize);
00080 
00081     float64_t TolAbs=0;
00082     float64_t QPBound=0;
00083     int32_t Method=0;
00084     if (method == SVM_OCAS)
00085         Method = 1;
00086     ocas_return_value_T result = svm_ocas_solver( get_C1(), num_vec, get_epsilon(),
00087             TolAbs, QPBound, bufsize, Method, 
00088             &CSVMOcas::compute_W,
00089             &CSVMOcas::update_W, 
00090             &CSVMOcas::add_new_cut, 
00091             &CSVMOcas::compute_output,
00092             &CSVMOcas::sort,
00093             this);
00094 
00095     SG_INFO("Ocas Converged after %d iterations\n"
00096             "==================================\n"
00097             "timing statistics:\n"
00098             "output_time: %f s\n"
00099             "sort_time: %f s\n"
00100             "add_time: %f s\n"
00101             "w_time: %f s\n"
00102             "solver_time %f s\n"
00103             "ocas_time %f s\n\n", result.nIter, result.output_time, result.sort_time,
00104             result.add_time, result.w_time, result.solver_time, result.ocas_time);
00105 
00106     delete[] tmp_a_buf;
00107 
00108     uint32_t num_cut_planes = result.nCutPlanes;
00109 
00110     for (uint32_t i=0; i<num_cut_planes; i++)
00111     {
00112         delete[] cp_value[i];
00113         delete[] cp_index[i];
00114     }
00115 
00116     delete[] cp_value;
00117     cp_value=NULL;
00118     delete[] cp_index;
00119     cp_index=NULL;
00120     delete[] cp_nz_dims;
00121     cp_nz_dims=NULL;
00122     delete[] cp_bias;
00123     cp_bias=NULL;
00124 
00125     delete[] lab;
00126     lab=NULL;
00127 
00128     delete[] old_w;
00129     old_w=NULL;
00130 
00131     return true;
00132 }
00133 
00134 /*----------------------------------------------------------------------------------
00135   sq_norm_W = sparse_update_W( t ) does the following:
00136 
00137   W = oldW*(1-t) + t*W;
00138   sq_norm_W = W'*W;
00139 
00140   ---------------------------------------------------------------------------------*/
00141 float64_t CSVMOcas::update_W( float64_t t, void* ptr )
00142 {
00143   float64_t sq_norm_W = 0;         
00144   CSVMOcas* o = (CSVMOcas*) ptr;
00145   uint32_t nDim = (uint32_t) o->w_dim;
00146   float64_t* W=o->w;
00147   float64_t* oldW=o->old_w;
00148 
00149   for(uint32_t j=0; j <nDim; j++)
00150   {
00151       W[j] = oldW[j]*(1-t) + t*W[j];
00152       sq_norm_W += W[j]*W[j];
00153   }          
00154   o->bias=o->old_bias*(1-t) + t*o->bias;
00155   sq_norm_W += CMath::sq(o->bias);
00156 
00157   return( sq_norm_W );
00158 }
00159 
00160 /*----------------------------------------------------------------------------------
00161   sparse_add_new_cut( new_col_H, new_cut, cut_length, nSel ) does the following:
00162 
00163     new_a = sum(data_X(:,find(new_cut ~=0 )),2);
00164     new_col_H = [sparse_A(:,1:nSel)'*new_a ; new_a'*new_a];
00165     sparse_A(:,nSel+1) = new_a;
00166 
00167   ---------------------------------------------------------------------------------*/
00168 void CSVMOcas::add_new_cut(
00169     float64_t *new_col_H, uint32_t *new_cut, uint32_t cut_length,
00170     uint32_t nSel, void* ptr)
00171 {
00172     CSVMOcas* o = (CSVMOcas*) ptr;
00173     CDotFeatures* f = o->features;
00174     uint32_t nDim=(uint32_t) o->w_dim;
00175     float64_t* y = o->lab;
00176 
00177     float64_t** c_val = o->cp_value;
00178     uint32_t** c_idx = o->cp_index;
00179     uint32_t* c_nzd = o->cp_nz_dims;
00180     float64_t* c_bias = o->cp_bias;
00181 
00182     float64_t sq_norm_a;
00183     uint32_t i, j, nz_dims;
00184 
00185     /* temporary vector */
00186     float64_t* new_a = o->tmp_a_buf;
00187     memset(new_a, 0, sizeof(float64_t)*nDim);
00188 
00189     for(i=0; i < cut_length; i++) 
00190     {
00191         f->add_to_dense_vec(y[new_cut[i]], new_cut[i], new_a, nDim);
00192 
00193         if (o->use_bias)
00194             c_bias[nSel]+=y[new_cut[i]];
00195     }
00196 
00197     /* compute new_a'*new_a and count number of non-zerou dimensions */
00198     nz_dims = 0; 
00199     sq_norm_a = CMath::sq(c_bias[nSel]);
00200     for(j=0; j < nDim; j++ ) {
00201         if(new_a[j] != 0) {
00202             nz_dims++;
00203             sq_norm_a += new_a[j]*new_a[j];
00204         }
00205     }
00206 
00207     /* sparsify new_a and insert it to the last column of sparse_A */
00208     c_nzd[nSel] = nz_dims;
00209     if(nz_dims > 0)
00210     {
00211         c_idx[nSel]=new uint32_t[nz_dims];
00212         c_val[nSel]=new float64_t[nz_dims];
00213 
00214         uint32_t idx=0;
00215         for(j=0; j < nDim; j++ )
00216         {
00217             if(new_a[j] != 0)
00218             {
00219                 c_idx[nSel][idx] = j;
00220                 c_val[nSel][idx++] = new_a[j];
00221             }
00222         }
00223     }
00224 
00225     new_col_H[nSel] = sq_norm_a;
00226 
00227     for(i=0; i < nSel; i++)
00228     {
00229         float64_t tmp = c_bias[nSel]*c_bias[i];
00230         for(j=0; j < c_nzd[i]; j++)
00231             tmp += new_a[c_idx[i][j]]*c_val[i][j];
00232 
00233         new_col_H[i] = tmp;
00234     }
00235     //CMath::display_vector(new_col_H, nSel+1, "new_col_H");
00236     //CMath::display_vector((int32_t*) c_idx[nSel], (int32_t) nz_dims, "c_idx");
00237     //CMath::display_vector((float64_t*) c_val[nSel], nz_dims, "c_val");
00238 }
00239 
00240 void CSVMOcas::sort(float64_t* vals, uint32_t* idx, uint32_t size)
00241 {
00242     CMath::qsort_index(vals, idx, size);
00243 }
00244 
00245 /*----------------------------------------------------------------------
00246   sparse_compute_output( output ) does the follwing:
00247 
00248   output = data_X'*W;
00249   ----------------------------------------------------------------------*/
00250 void CSVMOcas::compute_output(float64_t *output, void* ptr)
00251 {
00252     CSVMOcas* o = (CSVMOcas*) ptr;
00253     CDotFeatures* f=o->features;
00254     int32_t nData=f->get_num_vectors();
00255 
00256     float64_t* y = o->lab;
00257 
00258     f->dense_dot_range(output, 0, nData, y, o->w, o->w_dim, 0.0);
00259 
00260     for (int32_t i=0; i<nData; i++)
00261         output[i]+=y[i]*o->bias;
00262     //CMath::display_vector(o->w, o->w_dim, "w");
00263     //CMath::display_vector(output, nData, "out");
00264 }
00265 
00266 /*----------------------------------------------------------------------
00267   sq_norm_W = compute_W( alpha, nSel ) does the following:
00268 
00269   oldW = W;
00270   W = sparse_A(:,1:nSel)'*alpha;
00271   sq_norm_W = W'*W;
00272   dp_WoldW = W'*oldW';
00273 
00274   ----------------------------------------------------------------------*/
00275 void CSVMOcas::compute_W(
00276     float64_t *sq_norm_W, float64_t *dp_WoldW, float64_t *alpha, uint32_t nSel,
00277     void* ptr )
00278 {
00279     CSVMOcas* o = (CSVMOcas*) ptr;
00280     uint32_t nDim= (uint32_t) o->w_dim;
00281     CMath::swap(o->w, o->old_w);
00282     float64_t* W=o->w;
00283     float64_t* oldW=o->old_w;
00284     memset(W, 0, sizeof(float64_t)*nDim);
00285     float64_t old_bias=o->bias;
00286     float64_t bias=0;
00287 
00288     float64_t** c_val = o->cp_value;
00289     uint32_t** c_idx = o->cp_index;
00290     uint32_t* c_nzd = o->cp_nz_dims;
00291     float64_t* c_bias = o->cp_bias;
00292 
00293     for(uint32_t i=0; i<nSel; i++)
00294     {
00295         uint32_t nz_dims = c_nzd[i];
00296 
00297         if(nz_dims > 0 && alpha[i] > 0)
00298         {
00299             for(uint32_t j=0; j < nz_dims; j++)
00300                 W[c_idx[i][j]] += alpha[i]*c_val[i][j];
00301         }
00302         bias += c_bias[i]*alpha[i];
00303     }
00304 
00305     *sq_norm_W = CMath::dot(W,W, nDim) + CMath::sq(bias);
00306     *dp_WoldW = CMath::dot(W,oldW, nDim) + bias*old_bias;
00307     //SG_PRINT("nSel=%d sq_norm_W=%f dp_WoldW=%f\n", nSel, *sq_norm_W, *dp_WoldW);
00308     
00309     o->bias = bias;
00310     o->old_bias = old_bias;
00311 }

SHOGUN Machine Learning Toolbox - Documentation