gmnplib.cpp

Go to the documentation of this file.
00001 /*-----------------------------------------------------------------------
00002  *
00003  * This program is free software; you can redistribute it and/or modify
00004  * it under the terms of the GNU General Public License as published by
00005  * the Free Software Foundation; either version 3 of the License, or
00006  * (at your option) any later version.
00007  *
00008  * Library of solvers for Generalized Nearest Point Problem (GNPP).
00009  *
00010  * Written (W) 1999-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
00011  * Copyright (C) 1999-2008 Center for Machine Perception, CTU FEL Prague 
00012  *
00013  *
00014 gmnplib.c: Library of solvers for Generalized Minimal Norm Problem (GMNP).
00015  
00016  Generalized Minimal Norm Problem to solve is
00017   
00018   min 0.5*alpha'*H*alpha + c'*alpha
00019 
00020   subject to  sum(alpha) = 1,  alpha(i) >= 0
00021   
00022  H [dim x dim] is symmetric positive definite matrix.
00023  c [dim x 1] is an arbitrary vector.
00024 
00025  The precision of the found solution is given by
00026  the parameters tmax, tolabs and tolrel which
00027  define the stopping conditions:
00028  
00029  UB-LB <= tolabs      ->  exit_flag = 1   Abs. tolerance.
00030  UB-LB <= UB*tolrel   ->  exit_flag = 2   Relative tolerance.
00031  LB > th              ->  exit_flag = 3   Threshold on lower bound.
00032  t >= tmax            ->  exit_flag = 0   Number of iterations.
00033 
00034  UB ... Upper bound on the optimal solution.
00035  LB ... Lower bound on the optimal solution.
00036  t  ... Number of iterations.
00037  History ... Value of LB and UB wrt. number of iterations.
00038 
00039 
00040  The following algorithms are implemented:
00041  ..............................................
00042 
00043  - GMNP solver based on improved MDM algorithm 1 (u fixed v optimized)
00044     exitflag = gmnp_imdm( &get_col, diag_H, vector_c, dim,  
00045                  tmax, tolabs, tolrel, th, &alpha, &t, &History, verb  );
00046 
00047   For more info refer to V.Franc: Optimization Algorithms for Kernel 
00048   Methods. Research report. CTU-CMP-2005-22. CTU FEL Prague. 2005.
00049   ftp://cmp.felk.cvut.cz/pub/cmp/articles/franc/Franc-PhD.pdf .
00050 
00051  Modifications:
00052  09-sep-2005, VF
00053  24-jan-2005, VF
00054  26-nov-2004, VF
00055  25-nov-2004, VF
00056  21-nov-2004, VF
00057  20-nov-2004, VF
00058  31-may-2004, VF
00059  23-Jan-2004, VF
00060 
00061 -------------------------------------------------------------------- */
00062 
00063 #include "classifier/svm/gmnplib.h"
00064 #include "lib/Mathematics.h"
00065 
00066 #include <string.h>
00067 #include <limits.h>
00068 
00069 #define HISTORY_BUF 1000000
00070 
00071 #define MINUS_INF INT_MIN
00072 #define PLUS_INF  INT_MAX
00073 
00074 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
00075 #define KDELTA(A,B) (A==B)
00076 #define KDELTA4(A1,A2,A3,A4) ((A1==A2)||(A1==A3)||(A1==A4)||(A2==A3)||(A2==A4)||(A3==A4))
00077 
00078 CGMNPLib::CGMNPLib(
00079     float64_t* vector_y, CKernel* kernel, int32_t num_data,
00080     int32_t num_virt_data, int32_t num_classes, float64_t reg_const)
00081 : CSGObject()
00082 {
00083   m_num_classes=num_classes;
00084   m_num_virt_data=num_virt_data;
00085   m_reg_const = reg_const;
00086   m_num_data = num_data;
00087   m_vector_y = vector_y;
00088   m_kernel = kernel;
00089 
00090   Cache_Size = ((int64_t) kernel->get_cache_size())*1024*1024/(sizeof(float64_t)*num_data);
00091   Cache_Size = CMath::min(Cache_Size, (int64_t) num_data);
00092 
00093   SG_INFO("using %d kernel cache lines\n", Cache_Size);
00094   ASSERT(Cache_Size>=2);
00095 
00096   /* allocates memory for kernel cache */
00097   kernel_columns = new float64_t*[Cache_Size];
00098   cache_index = new float64_t[Cache_Size];
00099 
00100   for(int32_t i = 0; i < Cache_Size; i++ )
00101   {
00102     kernel_columns[i] = new float64_t[num_data];
00103     cache_index[i] = -2;
00104   }
00105   first_kernel_inx = 0;
00106 
00107 
00108 
00109   for(int32_t i = 0; i < 3; i++ )
00110   {
00111     virt_columns[i] = new float64_t[num_virt_data];
00112   }
00113   first_virt_inx = 0;
00114 
00115   diag_H = new float64_t[num_virt_data];
00116 
00117   for(int32_t i = 0; i < num_virt_data; i++ )
00118       diag_H[i] = kernel_fce(i,i);
00119 }
00120 
00121 CGMNPLib::~CGMNPLib()
00122 {
00123     for(int32_t i = 0; i < Cache_Size; i++ ) 
00124         delete[] kernel_columns[i];
00125 
00126     for(int32_t i = 0; i < 3; i++ ) 
00127         delete[] virt_columns[i];
00128 
00129     delete[] cache_index;
00130     delete[] kernel_columns;
00131 
00132     delete[] diag_H;
00133 }
00134 
00135 /* ------------------------------------------------------------
00136   Returns pointer at a-th column of the kernel matrix.
00137   This function maintains FIFO cache of kernel columns.
00138 ------------------------------------------------------------ */
00139 float64_t* CGMNPLib::get_kernel_col( int32_t a )
00140 {
00141   float64_t *col_ptr;
00142   int32_t i;
00143   int32_t inx;
00144 
00145   inx = -1;
00146   for( i=0; i < Cache_Size; i++ ) {
00147     if( cache_index[i] == a ) { inx = i; break; }
00148   }
00149     
00150   if( inx != -1 ) {
00151     col_ptr = kernel_columns[inx];
00152     return( col_ptr );
00153   }
00154    
00155   col_ptr = kernel_columns[first_kernel_inx];
00156   cache_index[first_kernel_inx] = a;
00157 
00158   first_kernel_inx++;
00159   if( first_kernel_inx >= Cache_Size ) first_kernel_inx = 0;
00160 
00161   for( i=0; i < m_num_data; i++ ) {
00162     col_ptr[i] = m_kernel->kernel(i,a);
00163   }
00164 
00165   return( col_ptr );
00166 }
00167 
00168 /* ------------------------------------------------------------
00169   Computes index of input example and its class label from 
00170   index of virtual "single-class" example.
00171 ------------------------------------------------------------ */
00172 void CGMNPLib::get_indices2( int32_t *index, int32_t *c, int32_t i )
00173 {
00174    *index = i / (m_num_classes-1);
00175  
00176    *c= (i % (m_num_classes-1))+1;
00177    if( *c>= m_vector_y[ *index ]) (*c)++;
00178 
00179    return;
00180 }
00181 
00182 
00183 /* ------------------------------------------------------------
00184   Returns pointer at the a-th column of the virtual K matrix.
00185 
00186   (note: the b-th column must be preserved in the cache during 
00187    updating but b is from (a(t-2), a(t-1)) where a=a(t) and
00188    thus FIFO with three columns does not have to take care od b.)
00189 ------------------------------------------------------------ */
00190 float64_t* CGMNPLib::get_col( int32_t a, int32_t b )
00191 {
00192   int32_t i;
00193   float64_t *col_ptr;
00194   float64_t *ker_ptr;
00195   float64_t value;
00196   int32_t i1,c1,i2,c2;
00197 
00198   col_ptr = virt_columns[first_virt_inx++];
00199   if( first_virt_inx >= 3 ) first_virt_inx = 0;
00200 
00201   get_indices2( &i1, &c1, a );
00202   ker_ptr = (float64_t*) get_kernel_col( i1 );
00203 
00204   for( i=0; i < m_num_virt_data; i++ ) {
00205     get_indices2( &i2, &c2, i );
00206 
00207     if( KDELTA4(m_vector_y[i1],m_vector_y[i2],c1,c2) ) {
00208       value = (+KDELTA(m_vector_y[i1],m_vector_y[i2]) 
00209                -KDELTA(m_vector_y[i1],c2)
00210                -KDELTA(m_vector_y[i2],c1)
00211                +KDELTA(c1,c2)
00212               )*(ker_ptr[i2]+1);
00213     }
00214     else
00215     {
00216       value = 0;
00217     }
00218 
00219     if(a==i) value += m_reg_const; 
00220 
00221     col_ptr[i] = value;
00222   }
00223   
00224   return( col_ptr );
00225 }
00226 
00227 
00228 /* --------------------------------------------------------------
00229  GMNP solver based on improved MDM algorithm 1.
00230 
00231  Search strategy: u determined by common rule and v is 
00232  optimized.
00233 
00234  Usage: exitflag = gmnp_imdm( &get_col, diag_H, vector_c, dim,  
00235                   tmax, tolabs, tolrel, th, &alpha, &t, &History );
00236 -------------------------------------------------------------- */
00237 
00238 int8_t CGMNPLib::gmnp_imdm(float64_t *vector_c,
00239             int32_t dim, 
00240             int32_t tmax,
00241             float64_t tolabs,
00242             float64_t tolrel,
00243             float64_t th,
00244             float64_t *alpha,
00245             int32_t  *ptr_t,
00246             float64_t **ptr_History,
00247             int32_t verb)
00248 {
00249   float64_t LB;
00250   float64_t UB;
00251   float64_t aHa, ac;
00252   float64_t tmp, tmp1;
00253   float64_t Huu, Huv, Hvv;
00254   float64_t min_beta, beta;
00255   float64_t max_improv, improv;
00256   float64_t lambda;
00257   float64_t *History;
00258   float64_t *Ha;
00259   float64_t *tmp_ptr;
00260   float64_t *col_u, *col_v;
00261   int32_t u=0;
00262   int32_t v=0;
00263   int32_t new_u=0;
00264   int32_t i;
00265   int32_t t;
00266   int32_t History_size;
00267   int8_t exitflag;
00268 
00269   /* ------------------------------------------------------------ */
00270   /* Initialization                                               */
00271   /* ------------------------------------------------------------ */
00272 
00273   Ha = new float64_t[dim];
00274   if( Ha == NULL ) SG_ERROR("Not enough memory.");
00275 
00276   History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
00277   History = new float64_t[History_size*2];
00278   if( History == NULL ) SG_ERROR("Not enough memory.");
00279 
00280   /* inx = argmin(0.5*diag_H + vector_c ); */
00281   for( tmp1 =  PLUS_INF, i = 0; i < dim; i++ ) {
00282     tmp = 0.5*diag_H[i] + vector_c[i];
00283     if( tmp1 > tmp) {
00284       tmp1 = tmp;
00285       v = i;
00286     }
00287   }
00288 
00289   col_v = (float64_t*)get_col(v,-1);
00290 
00291   for( min_beta = PLUS_INF, i = 0; i < dim; i++ ) 
00292   {
00293     alpha[i] = 0;
00294     Ha[i] = col_v[i];
00295 
00296     beta = Ha[i] + vector_c[i];
00297     if( beta < min_beta ) {
00298       min_beta = beta;
00299       u = i;
00300     }
00301   }
00302 
00303   alpha[v] = 1;
00304   aHa = diag_H[v];
00305   ac = vector_c[v];
00306 
00307   UB = 0.5*aHa + ac;
00308   LB = min_beta - 0.5*aHa;
00309 
00310   t = 0;
00311   History[INDEX(0,0,2)] = LB;
00312   History[INDEX(1,0,2)] = UB;
00313 
00314   if( verb ) {
00315     SG_PRINT("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00316       UB, LB, UB-LB,(UB-LB)/UB);
00317   }  
00318 
00319   /* Stopping conditions */
00320   if( UB-LB <= tolabs ) exitflag = 1;
00321   else if(UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00322   else if(LB > th ) exitflag = 3;
00323   else exitflag = -1;
00324 
00325   /* ------------------------------------------------------------ */
00326   /* Main optimization loop                                       */
00327   /* ------------------------------------------------------------ */
00328 
00329   col_u = (float64_t*)get_col(u,-1);
00330   while( exitflag == -1 ) 
00331   {
00332     t++;     
00333 
00334     col_v = (float64_t*)get_col(v,u);
00335 
00336     /* Adaptation rule and update */
00337     Huu = diag_H[u];
00338     Hvv = diag_H[v];
00339     Huv = col_u[v];
00340 
00341     lambda = (Ha[v]-Ha[u]+vector_c[v]-vector_c[u])/(alpha[v]*(Huu-2*Huv+Hvv));
00342     if( lambda < 0 ) lambda = 0; else if (lambda > 1) lambda = 1;
00343 
00344     aHa = aHa + 2*alpha[v]*lambda*(Ha[u]-Ha[v])+
00345                 lambda*lambda*alpha[v]*alpha[v]*(Huu-2*Huv+Hvv);
00346 
00347     ac = ac + lambda*alpha[v]*(vector_c[u]-vector_c[v]);
00348 
00349     tmp = alpha[v];
00350     alpha[u]=alpha[u]+lambda*alpha[v];
00351     alpha[v]=alpha[v]-lambda*alpha[v];
00352 
00353     UB = 0.5*aHa + ac;
00354     
00355 /*    max_beta = MINUS_INF;*/
00356     for( min_beta = PLUS_INF, i = 0; i < dim; i++ ) 
00357     {
00358        Ha[i] = Ha[i] + lambda*tmp*(col_u[i] - col_v[i]);
00359 
00360        beta = Ha[i]+ vector_c[i];
00361 
00362        if( beta < min_beta )
00363        { 
00364          new_u = i;
00365          min_beta = beta;
00366        }
00367     }    
00368 
00369     LB = min_beta - 0.5*aHa; 
00370     u = new_u;    
00371     col_u = (float64_t*)get_col(u,-1);
00372 
00373     /* search for optimal v while u is fixed */
00374     for( max_improv =  MINUS_INF, i = 0; i < dim; i++ ) {
00375 
00376       if( alpha[i] != 0 ) {
00377         beta = Ha[i] + vector_c[i];
00378 
00379         if( beta >= min_beta ) {
00380 
00381           tmp = diag_H[u] - 2*col_u[i] + diag_H[i];
00382           if( tmp != 0 ) {
00383             improv = (0.5*(beta-min_beta)*(beta-min_beta))/tmp;
00384 
00385             if( improv > max_improv ) {
00386               max_improv = improv;
00387               v = i;
00388             }
00389           }
00390         }
00391       }
00392     }
00393 
00394     /* Stopping conditions */
00395     if( UB-LB <= tolabs ) exitflag = 1; 
00396     else if( UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00397     else if(LB > th ) exitflag = 3;
00398     else if(t >= tmax) exitflag = 0; 
00399 
00400     /* print info */
00401     if(verb && (t % verb) == 0 ) {
00402       SG_PRINT("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00403         t, UB, LB, UB-LB,(UB-LB)/UB);
00404     }  
00405 
00406     /* Store selected values */
00407     if( t < History_size ) {
00408       History[INDEX(0,t,2)] = LB;
00409       History[INDEX(1,t,2)] = UB;
00410     }
00411     else {
00412       tmp_ptr = new float64_t[(History_size+HISTORY_BUF)*2];
00413       if( tmp_ptr == NULL ) SG_ERROR("Not enough memory.");
00414       for( i = 0; i < History_size; i++ ) {
00415         tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
00416         tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
00417       }
00418       tmp_ptr[INDEX(0,t,2)] = LB;
00419       tmp_ptr[INDEX(1,t,2)] = UB;
00420       
00421       History_size += HISTORY_BUF;
00422       delete[] History;
00423       History = tmp_ptr;
00424     }
00425   }
00426 
00427   /* print info about last iteration*/
00428   if(verb && (t % verb) ) {
00429     SG_PRINT("exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00430       UB, LB, UB-LB,(UB-LB)/UB);
00431   }  
00432 
00433 
00434   /*------------------------------------------------------- */
00435   /* Set outputs                                            */
00436   /*------------------------------------------------------- */
00437   (*ptr_t) = t;
00438   (*ptr_History) = History;
00439 
00440   /* Free memory */
00441   delete[] Ha;
00442   
00443   return( exitflag ); 
00444 }
00445 
00446 /* ------------------------------------------------------------
00447   Retures (a,b)-th element of the virtual kernel matrix 
00448   of size [num_virt_data x num_virt_data]. 
00449 ------------------------------------------------------------ */
00450 float64_t CGMNPLib::kernel_fce( int32_t a, int32_t b )
00451 {
00452   float64_t value;
00453   int32_t i1,c1,i2,c2;
00454 
00455   get_indices2( &i1, &c1, a );
00456   get_indices2( &i2, &c2, b );
00457 
00458   if( KDELTA4(m_vector_y[i1],m_vector_y[i2],c1,c2) ) {
00459     value = (+KDELTA(m_vector_y[i1],m_vector_y[i2]) 
00460              -KDELTA(m_vector_y[i1],c2)
00461              -KDELTA(m_vector_y[i2],c1)
00462              +KDELTA(c1,c2)
00463             )*(m_kernel->kernel( i1, i2 )+1);
00464   }
00465   else
00466   {
00467     value = 0;
00468   }
00469 
00470   if(a==b) value += m_reg_const; 
00471 
00472   return( value );
00473 }

SHOGUN Machine Learning Toolbox - Documentation