smoker_l1o.c

支持向量机作为统计学习理论的实现方法

/* --------------------------------------------------------------------
 [l1o_err, margin, error] = smoker_l1o( K, labels, C, eps, tol, verb )
 Leave One Out validation for SVM trained by SMO.

  mandatory inputs:
   K [N x N ] kernel matrix of N training patterns.
   labels [1 x N] pattern labels (1 for 1st class, 2 for 2nd class ).
   C [real] or [2 x real] one trade-off constant for both the classes
     or two constants for the first and the second class.

  optional inputs:
   eps [real] tolerance of KKT-conditions fulfilment (default 0.001).
   tol [real] minimal change of optimized Lagrangeians (default 0.001).
   verb [int] if 1 then the progress is displayed.

  mandatory outputs:
   l1o_err [real] leave one out error.

  optional outputs:
   margin [real] margin of the classifier of the whole problem.
   error  [real] training error of the whole problem.

 Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
 (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
 Written Vojtech Franc (diploma thesis) 02.11.1999

 Modifications:
 21-Oct-2001, V.Franc
 18-Oct-2001, V.Franc, adopted from svm_smocker.c.
 -------------------------------------------------------------------- */

#include "mex.h"
#include "matrix.h"
#include <math.h>
#include <stdlib.h>

/* if RANDOM is defined then a random element is used within
 optimization procedure as originally suggested. */
#define RANDOM

#define ZERO_LIM   1e-9     /* patterns with alpha <= ZERO_LIM are SV */

#define MAX(A,B)   (((A) > (B)) ? (A) : (B) )
#define MIN(A,B)   (((A) < (B)) ? (A) : (B) )

#define C(arg)   (target[arg] > 0 ? C1 : C2)
#define K(A,B)   (kernel_matrix[A+B*N])

/* --- Global variables ---------------------------------------------- */

long N = 0;                /* number of training patterns */
double C1=1;               /* trade-off constant for the 1st class */
double C2=1;               /* trade-off constant for the 2nd class */
double tolerance=0.001;    /* tolerance in KKT fulfilment  */
double eps=0.001;          /* minimal Lagrangeian change */
double *target;            /* pointer at labels */
double *error_cache;       /* error cache */

long out = 0;              /* the left out pattern */

double *kernel_matrix;     /* precomputed kernel matrix */

double *alpha;             /* Lagrange multipliers */
double *b;                 /* Bias (threshold) */


/* ==============================================================
 Implementation of Sequential Minimal Optimizer (SMO)
============================================================== */

/* --------------------------------------------------------------
 Computes value of the learned function for k-th pattern.
-------------------------------------------------------------- */
double learned_func( long k )
{
   double s = 0.;
   long i;
   for( i = 0; i < N; i++ ) {
     if( i != out) {
       if( alpha[i] > 0 )
           s += alpha[i]*target[i]*K(i,k);
     }
   }
  s -= *b;
  return( s );
}

/*---------------------------------*/
double margin()
{
   double s = 0.;
   long i,j;

   for( i = 0; i < N; i++ ) {
     for(j = 0; j < N; j++ ) {
       if( i != out & j != out & alpha[i] > 0 & alpha[j] > 0) {
           s += alpha[i]*alpha[j]*target[i]*target[j]*K(i,j);
       }
     }
   }
  
   if( s > 0)
      s = 1/sqrt(s);
  return( s );
}


/* --------------------------------------------------------------
 Optimizes objective function for i1-th and i2-th pattern.
-------------------------------------------------------------- */
long takeStep( long i1, long i2 ) {
   double y1, y2, s;
   long i;
   double alpha1, alpha2; 
   double a1, a2;
   double E1, E2, L, H, k11, k22, k12, eta, Lobj, Hobj;
   double gamma;
   double c1, c2;
   double t;
   double b1, b2, bnew;
   double delta_b;
   double t1, t2;

   if( i1 == i2 ) return( 0 );

   alpha1 = alpha[i1];
   y1 = target[i1];

   /*   E1 = error_cache[i1];  */
   if( alpha1 > 0 && alpha1 < C(i1) )
      E1 = error_cache[i1]; 
   else 
      E1 = learned_func(i1) - y1; 

   alpha2 = alpha[i2];
   y2 = target[i2];

   /*   E2 = error_cache[i2]; */
   if( alpha2 > 0 && alpha2 < C(i2) )
      E2 = error_cache[i2];  
   else  
      E2 = learned_func(i2) - y2;  

   s = y1 * y2;
   if(s < 0)
   {
      L = MAX(0, alpha2 - alpha1);
      H = MIN(C(i2), C(i1) + alpha2 - alpha1);
   }
   else
   {
     L = MAX(0, alpha2 + alpha1 - C(i1) );
     H = MIN(C(i2), alpha2 + alpha1);
   }

   if( L == H ) return( 0 );

   k11 = K(i1,i1);
   k12 = K(i1,i2);
   k22 = K(i2,i2);
   eta = 2 * k12 - k11 - k22;

   if( eta < 0 ) {
      a2 = alpha2 + y2 * (E2 - E1) / eta;
      if( a2 < L )
         a2 = L;
      else if( a2 > H )
         a2 = H;
   }
   else {
      c1 = eta/2;
      c2 = y2 * (E1-E2)- eta * alpha2;
      Lobj = c1 * L * L + c2 * L;
      Hobj = c1 * H * H + c2 * H;

      if( Lobj > Hobj+eps )
         a2 = L;
      else if( Lobj < Hobj-eps )
         a2 = H;
      else
         a2 = alpha2;
   }

   /*
   if( a2 < ZERO_LIM )
      a2 = 0;
   else if( a2 > C(i2)-ZERO_LIM )
      a2 = C(i2);
   */

   if( fabs(a2-alpha2) < eps*(a2+alpha2+eps )) return( 0 );

   a1 = alpha1 - s * (a2 - alpha2 );
   if( a1 < 0 ) {
      a2 += s * a1;
      a1 = 0;
   }
   else if( a1 > C(i1) ) {
      t = a1-C(i1);
      a2 += s * t;
      a1 = C(i1);
   }

   if( a1 > 0 && a1 < C(i1) )
      bnew = *b + E1 + y1 * (a1 - alpha1) * k11 + y2 * (a2 - alpha2) * k12;
   else {
      if( a2 > 0 && a2 < C(i2) )
         bnew = *b + E2 + y1 *(a1 - alpha1)*k12 + y2*(a2 - alpha2) * k22;
      else {
         b1 = *b + E1 + y1 * (a1 - alpha1) * k11 + y2 * (a2 - alpha2) * k12;
         b2 = *b + E2 + y1 * (a1 - alpha1) * k12 + y2 * (a2 - alpha2) * k22;
         bnew = (b1 + b2) / 2;
      }
   }

   delta_b = bnew - *b;
   *b = bnew;

   t1 = y1 * (a1-alpha1);
   t2 = y2 * (a2-alpha2);

   for( i = 0; i < N; i++ ) {
     if( i != out )
     {
       if (0 < alpha[i] && alpha[i] < C(i)) { 
          error_cache[i] +=  t1 * K(i1,i) + t2 * K(i2,i) - delta_b; 
       } 
     } 
   }

   error_cache[i1] = 0; 
   error_cache[i2] = 0; 
   
   /*
   for( i = 0; i < N; i++ ) { 
      error_cache[i] +=  t1 * K(i1,i) + t2 * K(i2,i) - delta_b; 
   } 
   */

   alpha[i1] = a1;  /* Store a1 in the alpha array.*/
   alpha[i2] = a2;  /* Store a2 in the alpha array.*/

   return( 1 );
}


/* --------------------------------------------------------------
 Finds the second Lagrange multiplayer to be optimize.
-------------------------------------------------------------- */
long examineExample( long i1 )
{
   double y1, alpha1, E1, r1;
   double tmax;
   double E2, temp;
   long k, i2;
   long k0;

   y1 = target[i1];
   alpha1 = alpha[i1];

   /*   E1 = error_cache[i1]; */

   if( alpha1 > 0 && alpha1 < C(i1) )
      E1 = error_cache[i1];  
   else  
      E1 = learned_func(i1) - y1; 

   r1 = y1 * E1;
   if(( r1 < -tolerance && alpha1 < C(i1) )
      || (r1 > tolerance && alpha1 > 0)) {
    /* Try i2 by three ways; if successful, then immediately return 1; */

      for( i2 = (-1), tmax = 0, k = 0; k < N; k++ ) {
        if( k != out ) 
        {
           if( alpha[k] > 0 && alpha[k] < C(k) ) {
              E2 = error_cache[k];
              temp = fabs(E1 - E2);
              if( temp > tmax ) {
                 tmax = temp;
                 i2 = k;
              }
           }
        }
      }
      if( i2 >= 0 ) {
          if( takeStep(i1,i2) )
             return( 1 );
      }

#ifdef RANDOM
      for( k0 = rand(), k = k0; k < N + k0; k++ ) {
         i2 = k % N;
#else
      for( k = 0; k < N; k++) {
         i2 = k;
#endif
         if( i2 != out)
         {
            if( alpha[i2] > 0 && alpha[i2] < C(i2) ) {
               if( takeStep(i1,i2) )
                  return( 1 );
            }
         }
      }

#ifdef RANDOM
      for( k0 = rand(), k = k0; k < N + k0; k++ ) {
         i2 = k % N;
#else
      for( k = 0; k < N; k++) {
         i2 = k;
#endif
         if( i2 != out)
         {
            if( takeStep(i1,i2) )
              return( 1 );
         }
      }

   } /* if( ... ) */

   return( 0 );
}

/* --------------------------------------------------------------
 Main SMO optimization cycle.
-------------------------------------------------------------- */
void runSMO( void )
{
   long numChanged = 0;
   long examineAll = 1;
   long k;

   while( numChanged > 0 || examineAll ) {
      numChanged = 0;

      if( examineAll ) {
         for( k = 0; k < N; k++ ) {
           if( k != out)
           {
              numChanged += examineExample( k );
           }
         }
      }
      else {
         for( k = 0; k < N; k++ ) {
           if( k != out)
           {
              if( alpha[k] != 0 && alpha[k] != C(k) )
                 numChanged += examineExample( k );
           }
         }
      }

      if( examineAll == 1 )
         examineAll = 0;
      else if( numChanged == 0 )
         examineAll = 1;
   }
}

/* ==============================================================
 Main MEX function - interface to Matlab.
============================================================== */
void mexFunction( int nlhs, mxArray *plhs[],
		  int nrhs, const mxArray*prhs[] )
{
   long i, j;
   double *labels12, *initAlpha, *nsv, *nerr, *tmp;
   double *init_alpha, *mar, *tot_err;
   int verb = 0;              /* verbosity */

   /* ---- check number of input arguments  ------------- */

   if(nrhs < 3)
      mexErrMsgTxt("Not enough input arguments.");
   if(nlhs < 1)
      mexErrMsgTxt("Not enough output arguments.");


   /* ---- check input arguments  ----------------------- */
   /* [Alpha,bias,nsv,nerr] = smocker( K, I, C, eps, tol, Alpha ) */

   /* data matrix [N x N ] */
   if( !mxIsNumeric(prhs[0]) || !mxIsDouble(prhs[0]) ||
       mxIsEmpty(prhs[0])    || mxIsComplex(prhs[0]) )
      mexErrMsgTxt("Input K must be a real matrix.");

   /* vector of labels (1,2) */
   if( !mxIsNumeric(prhs[1]) || !mxIsDouble(prhs[1]) ||
       mxIsEmpty(prhs[1])    || mxIsComplex(prhs[1]) ||
       (mxGetN(prhs[1]) != 1 && mxGetM(prhs[1]) != 1))
      mexErrMsgTxt("Input I must be a real vector.");

   /*  one or two real trade-off constant(s)  */
   if( !mxIsNumeric(prhs[2]) || !mxIsDouble(prhs[2]) ||
       mxIsEmpty(prhs[2])    || mxIsComplex(prhs[2]) ||
       (mxGetN(prhs[2]) != 1  && mxGetM(prhs[2]) != 1 ))
      mexErrMsgTxt("Input C must be one or two real scalar(s).");
   else {
      if( MAX( mxGetN(prhs[2]), mxGetM(prhs[2])) == 1 ) {
         C1 = mxGetScalar(prhs[2]);
         C2 = C1;
      }
      else
      {
         tmp = mxGetPr(prhs[2]);
         C1 = tmp[0];
         C2 = tmp[1];
      }
   }

   /* real parameter eps */
   if( nrhs >= 4 ) {
      if( !mxIsNumeric(prhs[3]) || !mxIsDouble(prhs[3]) ||
         mxIsEmpty(prhs[3])    || mxIsComplex(prhs[3]) ||
         mxGetN(prhs[3]) != 1  || mxGetM(prhs[3]) != 1 )
         mexErrMsgTxt("Input eps must be a real scalar.");
      else
         eps = mxGetScalar(prhs[3]);   /* take eps argument */
   }

   /* real parameter tol */
   if(nrhs >= 5) {
      if( !mxIsNumeric(prhs[4]) || !mxIsDouble(prhs[4]) ||
         mxIsEmpty(prhs[4])    || mxIsComplex(prhs[4]) ||
         mxGetN(prhs[4]) != 1  || mxGetM(prhs[4]) != 1 )
         mexErrMsgTxt("Input tol must be a real scalar.");
      else
         tolerance = mxGetScalar(prhs[4]);  /* take tolerance argument */
   }

   /* verb */
   if( nrhs >= 6 ) {
      if( !mxIsNumeric(prhs[5]) || !mxIsDouble(prhs[5]) ||
         mxIsEmpty(prhs[5])    || mxIsComplex(prhs[5]) ||
         mxGetN(prhs[5]) != 1  || mxGetM(prhs[5]) != 1 )
         mexErrMsgTxt("Input verb must be an integer scalar.");
      else
         verb = mxGetScalar(prhs[5]);   /* take eps argument */
   }

   /* ---- init variables ------------------------------- */

   kernel_matrix = mxGetPr(prhs[0]);  /* pointer at kernel matrix */
   N = mxGetN(prhs[0]);               /* number of data */
   labels12 = mxGetPr(prhs[1]);       /* labels (1,2) */

   /* allocate memory for targets (labels) (1,-1) */
   if( (target = mxCalloc(N, sizeof(double) )) == NULL) {
      mexErrMsgTxt("Not enough memory.");
   }

   /* transform labels12 (1,2) from to targets (1,-1) */
   for( i = 0; i < N; i++ ) {
      target[i] = - labels12[i]*2 + 3;
   }

   /* create output variable for bias */
   if( (b = mxCalloc(1, sizeof(double) )) == NULL) {
      mexErrMsgTxt("Not enough memory.");
   }
   *b = 0;

   /* allocate memory for error_cache */
   if( (error_cache = mxCalloc(N, sizeof(double) )) == NULL) {
      mexErrMsgTxt("Not enough memory for error cache.");
   }

   /* create vector for Lagrangeians */
   if( (alpha = mxCalloc(N, sizeof(double) )) == NULL) {
      mexErrMsgTxt("Not enough memory.");
   }

   for( i = 0; i < N; i++ ) {
        alpha[i] = 0; 
   }

   /*---------------------------------------------------*/
   /* finds solution for all patterns and store it */

   if( (init_alpha = mxCalloc(N, sizeof(double) )) == NULL) {
      mexErrMsgTxt("Not enough memory.");
   } 
   out = -1;

   runSMO();

   /* computes margin of the solution of the whole problem */
   if( nlhs >= 2) {

     plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
     mar = mxGetPr(plhs[1]);
     *mar = margin();
   }

   /* computes total error on the whole problem */ 
   if( nlhs >= 3) {

     plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL);
     tot_err = mxGetPr(plhs[2]);
     
     *tot_err = 0;
     for(i=1; i < N; i++) {
         if( alpha[i] >= C(i)-tolerance ) (*tot_err)++;
     }
     *tot_err= (*tot_err) / N;
   
   }


   for(i = 0; i < N; i++ ) {
     init_alpha[i] = alpha[i];
    }
  
   plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
   nerr = mxGetPr(plhs[0]);

   *nerr=0;

   for(i = 0; i < N; i++ ) 
   {
     /* rebuilds rule for SVs */
     if( init_alpha[i] != 0 )
     {
     
       out = i;

       if( verb == 1) 
         mexPrintf( ".");
       else if( verb ==2)
           mexPrintf("%d: ",i);
 
       /* set init solution */   
       *b = 0;
       for( j=0; j < N; j++) {
           alpha[j]=0;
           error_cache[j] = 0;
       }
    
       runSMO();

       

       /* classifies the left out pattern */
       if( target[out]*learned_func( out ) <0)
         {
            (*nerr)++;
            if(verb == 2) mexPrintf("err\n");
         }
       else
         {
            if(verb == 2) mexPrintf("ok:\n");
         }
     }
    }

   (*nerr) = (*nerr)/N;

     if(verb == 1) mexPrintf("\n");

   /* ----- free memory --------------------------------------- */
   mxFree( error_cache );
   mxFree( target );
}