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(DREAL* vector_y, CKernel* kernel, INT num_data, INT num_virt_data, INT num_classes, DREAL reg_const)
00079 : CSGObject()
00080 {
00081   m_num_classes=num_classes;
00082   m_num_virt_data=num_virt_data;
00083   m_reg_const = reg_const;
00084   m_num_data = num_data;
00085   m_vector_y = vector_y;
00086   m_kernel = kernel;
00087 
00088   Cache_Size = ((LONG) kernel->get_cache_size())*1024*1024/(sizeof(DREAL)*num_data);
00089   Cache_Size = CMath::min(Cache_Size, (LONG) num_data);
00090 
00091   SG_INFO("using %d kernel cache lines\n", Cache_Size);
00092   ASSERT(Cache_Size>=2);
00093 
00094   /* allocates memory for kernel cache */
00095   kernel_columns = new DREAL*[Cache_Size];
00096   cache_index = new DREAL[Cache_Size];
00097 
00098   for(INT i = 0; i < Cache_Size; i++ ) 
00099   {
00100     kernel_columns[i] = new DREAL[num_data];
00101     cache_index[i] = -2;
00102   }
00103   first_kernel_inx = 0;
00104 
00105 
00106 
00107   for(INT i = 0; i < 3; i++ )
00108   {
00109     virt_columns[i] = new DREAL[num_virt_data];
00110   }
00111   first_virt_inx = 0;
00112 
00113   diag_H = new DREAL[num_virt_data];
00114 
00115   for(INT i = 0; i < num_virt_data; i++ )
00116       diag_H[i] = kernel_fce(i,i);
00117 }
00118 
00119 CGMNPLib::~CGMNPLib()
00120 {
00121     for(INT i = 0; i < Cache_Size; i++ ) 
00122         delete[] kernel_columns[i];
00123 
00124     for(INT i = 0; i < 3; i++ ) 
00125         delete[] virt_columns[i];
00126 
00127     delete[] cache_index;
00128     delete[] kernel_columns;
00129 
00130     delete[] diag_H;
00131 }
00132 
00133 /* ------------------------------------------------------------
00134   Returns pointer at a-th column of the kernel matrix.
00135   This function maintains FIFO cache of kernel columns.
00136 ------------------------------------------------------------ */
00137 DREAL* CGMNPLib::get_kernel_col( INT a ) 
00138 {
00139   double *col_ptr;
00140   INT i;
00141   INT inx;
00142 
00143   inx = -1;
00144   for( i=0; i < Cache_Size; i++ ) {
00145     if( cache_index[i] == a ) { inx = i; break; }
00146   }
00147     
00148   if( inx != -1 ) {
00149     col_ptr = kernel_columns[inx];
00150     return( col_ptr );
00151   }
00152    
00153   col_ptr = kernel_columns[first_kernel_inx];
00154   cache_index[first_kernel_inx] = a;
00155 
00156   first_kernel_inx++;
00157   if( first_kernel_inx >= Cache_Size ) first_kernel_inx = 0;
00158 
00159   for( i=0; i < m_num_data; i++ ) {
00160     col_ptr[i] = m_kernel->kernel(i,a);
00161   }
00162 
00163   return( col_ptr );
00164 }
00165 
00166 /* ------------------------------------------------------------
00167   Computes index of input example and its class label from 
00168   index of virtual "single-class" example.
00169 ------------------------------------------------------------ */
00170 void CGMNPLib::get_indices2( INT *index, INT *c, INT i )
00171 {
00172    *index = i / (m_num_classes-1);
00173  
00174    *c= (i % (m_num_classes-1))+1;
00175    if( *c>= m_vector_y[ *index ]) (*c)++;
00176 
00177    return;
00178 }
00179 
00180 
00181 /* ------------------------------------------------------------
00182   Returns pointer at the a-th column of the virtual K matrix.
00183 
00184   (note: the b-th column must be preserved in the cache during 
00185    updating but b is from (a(t-2), a(t-1)) where a=a(t) and
00186    thus FIFO with three columns does not have to take care od b.)
00187 ------------------------------------------------------------ */
00188 DREAL* CGMNPLib::get_col( INT a, INT b )
00189 {
00190   INT i;
00191   double *col_ptr;
00192   double *ker_ptr;
00193   double value;
00194   INT i1,c1,i2,c2;
00195 
00196   col_ptr = virt_columns[first_virt_inx++];
00197   if( first_virt_inx >= 3 ) first_virt_inx = 0;
00198 
00199   get_indices2( &i1, &c1, a );
00200   ker_ptr = (double*) get_kernel_col( i1 );
00201 
00202   for( i=0; i < m_num_virt_data; i++ ) {
00203     get_indices2( &i2, &c2, i );
00204 
00205     if( KDELTA4(m_vector_y[i1],m_vector_y[i2],c1,c2) ) {
00206       value = (+KDELTA(m_vector_y[i1],m_vector_y[i2]) 
00207                -KDELTA(m_vector_y[i1],c2)
00208                -KDELTA(m_vector_y[i2],c1)
00209                +KDELTA(c1,c2)
00210               )*(ker_ptr[i2]+1);
00211     }
00212     else
00213     {
00214       value = 0;
00215     }
00216 
00217     if(a==i) value += m_reg_const; 
00218 
00219     col_ptr[i] = value;
00220   }
00221   
00222   return( col_ptr );
00223 }
00224 
00225 
00226 /* --------------------------------------------------------------
00227  GMNP solver based on improved MDM algorithm 1.
00228 
00229  Search strategy: u determined by common rule and v is 
00230  optimized.
00231 
00232  Usage: exitflag = gmnp_imdm( &get_col, diag_H, vector_c, dim,  
00233                   tmax, tolabs, tolrel, th, &alpha, &t, &History );
00234 -------------------------------------------------------------- */
00235 
00236 int CGMNPLib::gmnp_imdm(double *vector_c,
00237             INT dim, 
00238             INT tmax,
00239             double tolabs,
00240             double tolrel,
00241             double th,
00242             double *alpha,
00243             INT  *ptr_t,
00244             double **ptr_History,
00245             INT verb)
00246 {
00247   double LB;
00248   double UB;
00249   double aHa, ac;
00250   double tmp, tmp1;
00251   double Huu, Huv, Hvv;
00252   double min_beta, beta;
00253   double max_improv, improv;
00254   double lambda;
00255   double *History;
00256   double *Ha;
00257   double *tmp_ptr;
00258   double *col_u, *col_v;
00259   INT u=0;
00260   INT v=0;
00261   INT new_u=0;
00262   INT i;
00263   INT t;
00264   INT History_size;
00265   int exitflag;
00266 
00267   /* ------------------------------------------------------------ */
00268   /* Initialization                                               */
00269   /* ------------------------------------------------------------ */
00270 
00271   Ha = new double[dim];
00272   if( Ha == NULL ) SG_ERROR("Not enough memory.");
00273 
00274   History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
00275   History = new double[History_size*2];
00276   if( History == NULL ) SG_ERROR("Not enough memory.");
00277 
00278   /* inx = argmin(0.5*diag_H + vector_c ); */
00279   for( tmp1 =  PLUS_INF, i = 0; i < dim; i++ ) {
00280     tmp = 0.5*diag_H[i] + vector_c[i];
00281     if( tmp1 > tmp) {
00282       tmp1 = tmp;
00283       v = i;
00284     }
00285   }
00286 
00287   col_v = (double*)get_col(v,-1);
00288 
00289   for( min_beta = PLUS_INF, i = 0; i < dim; i++ ) 
00290   {
00291     alpha[i] = 0;
00292     Ha[i] = col_v[i];
00293 
00294     beta = Ha[i] + vector_c[i];
00295     if( beta < min_beta ) {
00296       min_beta = beta;
00297       u = i;
00298     }
00299   }
00300 
00301   alpha[v] = 1;
00302   aHa = diag_H[v];
00303   ac = vector_c[v];
00304 
00305   UB = 0.5*aHa + ac;
00306   LB = min_beta - 0.5*aHa;
00307 
00308   t = 0;
00309   History[INDEX(0,0,2)] = LB;
00310   History[INDEX(1,0,2)] = UB;
00311 
00312   if( verb ) {
00313     SG_PRINT("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00314       UB, LB, UB-LB,(UB-LB)/UB);
00315   }  
00316 
00317   /* Stopping conditions */
00318   if( UB-LB <= tolabs ) exitflag = 1;
00319   else if(UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00320   else if(LB > th ) exitflag = 3;
00321   else exitflag = -1;
00322 
00323   /* ------------------------------------------------------------ */
00324   /* Main optimization loop                                       */
00325   /* ------------------------------------------------------------ */
00326 
00327   col_u = (double*)get_col(u,-1);
00328   while( exitflag == -1 ) 
00329   {
00330     t++;     
00331 
00332     col_v = (double*)get_col(v,u);
00333 
00334     /* Adaptation rule and update */
00335     Huu = diag_H[u];
00336     Hvv = diag_H[v];
00337     Huv = col_u[v];
00338 
00339     lambda = (Ha[v]-Ha[u]+vector_c[v]-vector_c[u])/(alpha[v]*(Huu-2*Huv+Hvv));
00340     if( lambda < 0 ) lambda = 0; else if (lambda > 1) lambda = 1;
00341 
00342     aHa = aHa + 2*alpha[v]*lambda*(Ha[u]-Ha[v])+
00343                 lambda*lambda*alpha[v]*alpha[v]*(Huu-2*Huv+Hvv);
00344 
00345     ac = ac + lambda*alpha[v]*(vector_c[u]-vector_c[v]);
00346 
00347     tmp = alpha[v];
00348     alpha[u]=alpha[u]+lambda*alpha[v];
00349     alpha[v]=alpha[v]-lambda*alpha[v];
00350 
00351     UB = 0.5*aHa + ac;
00352     
00353 /*    max_beta = MINUS_INF;*/
00354     for( min_beta = PLUS_INF, i = 0; i < dim; i++ ) 
00355     {
00356        Ha[i] = Ha[i] + lambda*tmp*(col_u[i] - col_v[i]);
00357 
00358        beta = Ha[i]+ vector_c[i];
00359 
00360        if( beta < min_beta )
00361        { 
00362          new_u = i;
00363          min_beta = beta;
00364        }
00365     }    
00366 
00367     LB = min_beta - 0.5*aHa; 
00368     u = new_u;    
00369     col_u = (double*)get_col(u,-1);
00370 
00371     /* search for optimal v while u is fixed */
00372     for( max_improv =  MINUS_INF, i = 0; i < dim; i++ ) {
00373 
00374       if( alpha[i] != 0 ) {
00375         beta = Ha[i] + vector_c[i];
00376 
00377         if( beta >= min_beta ) {
00378 
00379           tmp = diag_H[u] - 2*col_u[i] + diag_H[i];
00380           if( tmp != 0 ) {
00381             improv = (0.5*(beta-min_beta)*(beta-min_beta))/tmp;
00382 
00383             if( improv > max_improv ) {
00384               max_improv = improv;
00385               v = i;
00386             }
00387           }
00388         }
00389       }
00390     }
00391 
00392     /* Stopping conditions */
00393     if( UB-LB <= tolabs ) exitflag = 1; 
00394     else if( UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
00395     else if(LB > th ) exitflag = 3;
00396     else if(t >= tmax) exitflag = 0; 
00397 
00398     /* print info */
00399     if(verb && (t % verb) == 0 ) {
00400       SG_PRINT("%d: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00401         t, UB, LB, UB-LB,(UB-LB)/UB);
00402     }  
00403 
00404     /* Store selected values */
00405     if( t < History_size ) {
00406       History[INDEX(0,t,2)] = LB;
00407       History[INDEX(1,t,2)] = UB;
00408     }
00409     else {
00410       tmp_ptr = new double[(History_size+HISTORY_BUF)*2];
00411       if( tmp_ptr == NULL ) SG_ERROR("Not enough memory.");
00412       for( i = 0; i < History_size; i++ ) {
00413         tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
00414         tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
00415       }
00416       tmp_ptr[INDEX(0,t,2)] = LB;
00417       tmp_ptr[INDEX(1,t,2)] = UB;
00418       
00419       History_size += HISTORY_BUF;
00420       delete[] History;
00421       History = tmp_ptr;
00422     }
00423   }
00424 
00425   /* print info about last iteration*/
00426   if(verb && (t % verb) ) {
00427     SG_PRINT("exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
00428       UB, LB, UB-LB,(UB-LB)/UB);
00429   }  
00430 
00431 
00432   /*------------------------------------------------------- */
00433   /* Set outputs                                            */
00434   /*------------------------------------------------------- */
00435   (*ptr_t) = t;
00436   (*ptr_History) = History;
00437 
00438   /* Free memory */
00439   delete[] Ha;
00440   
00441   return( exitflag ); 
00442 }
00443 
00444 /* ------------------------------------------------------------
00445   Retures (a,b)-th element of the virtual kernel matrix 
00446   of size [num_virt_data x num_virt_data]. 
00447 ------------------------------------------------------------ */
00448 double CGMNPLib::kernel_fce( INT a, INT b )
00449 {
00450   double value;
00451   INT i1,c1,i2,c2;
00452 
00453   get_indices2( &i1, &c1, a );
00454   get_indices2( &i2, &c2, b );
00455 
00456   if( KDELTA4(m_vector_y[i1],m_vector_y[i2],c1,c2) ) {
00457     value = (+KDELTA(m_vector_y[i1],m_vector_y[i2]) 
00458              -KDELTA(m_vector_y[i1],c2)
00459              -KDELTA(m_vector_y[i2],c1)
00460              +KDELTA(c1,c2)
00461             )*(m_kernel->kernel( i1, i2 )+1);
00462   }
00463   else
00464   {
00465     value = 0;
00466   }
00467 
00468   if(a==b) value += m_reg_const; 
00469 
00470   return( value );
00471 }

SHOGUN Machine Learning Toolbox - Documentation