qpssvmlib.c

matlab最新统计模式识别工具箱

    if( verb > 0 && (exitflag > 0 || (t % verb)==0 )) {
       mexPrintf("%d: UB=%.10f, LB=%.10f, UB-LB=%.10f, (UB-LB)/|UB|=%.10f \n",
        t, UB, LB, UB-LB, (UB!=0) ? (UB-LB)/ABS(UB) : 0);      
    }    

  }

  /* -- Find which stopping consition has been used -------- */
  if( UB-LB < tolabs ) exitflag = 1;
  else if(UB-LB < ABS(UB)*tolrel ) exitflag = 2;
  else exitflag = 0;

  /*----------------------------------------------------------   
    Set up outputs                                          
  ---------------------------------------------------------- */
  (*ptr_t) = t;
  (*ptr_History) = History;

  /*---------------------------------------------------------- 
    Clean up
  ---------------------------------------------------------- */
  mxFree( Hx );
  mxFree( d );
  
  return( exitflag ); 

}


/* --------------------------------------------------------------
 QPSSVM solver 

 Usage: exitflag = qpssvm_imdm( &get_col, diag_H, f, b, I, x, n, tmax, 
         tolabs, tolrel, &t, &History, verb );   
-------------------------------------------------------------- */
int qpssvm_imdm(const void* (*get_col)(long),
                  double *diag_H,
                  double *f,
                  double b,
                  uint16_T *I,
                  double *x,
                  long n,
                  long tmax,
                  double tolabs,
                  double tolrel,
                  long *ptr_t,
                  double **ptr_History,
                  long verb)
{
  double *x_nequ;
/*  double *Hx;*/
  double *d;
  double *History;
  double *col_u, *col_v;
  double *tmp_ptr;
  double LB;
  double UB;
  double tmp;
  double improv;
  double tmp_num;
  double tmp_den;
  double tau;
  double delta;
  double yu;
  long *inx;
  long *nk;
  long m;
  long t;
  long u;
  long v;
  long k;
  long i, j;
  long History_size;
  int exitflag;
  
  /* ------------------------------------------------------------ 
    Initialization                                               
  ------------------------------------------------------------ */

  /* count cumber of constraints */
  for( i=0, m=0; i < n; i++ ) m = MAX(m,I[i]);

  /* alloc and initialize x_nequ */
  x_nequ = mxCalloc(m, sizeof(double));
  if( x_nequ == NULL ) mexErrMsgTxt("Not enough memory.");

  /* alloc Inx */
  inx = mxCalloc(m*n, sizeof(double));
  if( inx == NULL ) mexErrMsgTxt("Not enough memory.");

  nk = mxCalloc(m, sizeof(double));
  if( nk == NULL ) mexErrMsgTxt("Not enough memory.");

  for( i=0; i < m; i++ ) x_nequ[i] = b;
  for( i=0; i < n; i++ ) {
     k = I[i]-1;
     x_nequ[k] -= x[i];
     inx[INDEX(nk[k],k,n)] = i;
     nk[k]++;
  }
    
  /* alloc History [2 x HISTORY_BUF] */
  History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
  History = mxCalloc(History_size*2,sizeof(double));
  if( History == NULL ) mexErrMsgTxt("Not enough memory.");

  /* alloc d [n x 1] */
  d = mxCalloc(n, sizeof(double));
  if( d == NULL ) mexErrMsgTxt("Not enough memory.");
 
  /* d = H*x + f; */
  for( i=0; i < n; i++ ) {
    if( x[i] > 0 ) {
      col_u = (double*)get_col(i);
      for( j=0; j < n; j++ ) {
          d[j] += col_u[j]*x[i];
      }
    }
  }
  for( i=0; i < n; i++ ) d[i] += f[i];
  
  /* UB = 0.5*x'*(f+d); */
  /* LB = 0.5*x'*(f-d); */
  for( i=0, UB = 0, LB=0; i < n; i++) {
    UB += x[i]*(f[i]+d[i]);
    LB += x[i]*(f[i]-d[i]);
  }
  UB = 0.5*UB;
  LB = 0.5*LB;

  /*
  for k=1:m,
    tmp = min(d(find(I==k)));
    if tmp < 0, LB = LB + b*tmp; end
  end
  */
  
  for( i=0; i < m; i++ ) {
    for( j=0, tmp = PLUS_INF; j < nk[i]; j++ ) {
      tmp = MIN(tmp, d[inx[INDEX(j,i,n)]]);
    }
    if( tmp < 0) LB += b*tmp;
  }
  
  /*
  for( i=0; i < m; i++ ) {
    for( j=0, tmp = PLUS_INF; j < n; j++ ) {
      if( I[j]-1 == i ) tmp = MIN(tmp, d[j]);
    }
    if( tmp < 0) LB += b*tmp;
  }*/

  exitflag = 0;
  t = 0;
  History[INDEX(0,0,2)] = LB;
  History[INDEX(1,0,2)] = UB;


  /* -- Main loop ---------------------------------------- */
  while( (exitflag == 0) && (t < tmax)) 
  {
    t++;

    exitflag = 1;
    for( k=0; k < m; k++ ) 
    {       
      /*
      inx = find(I==k);
      [tmp,u] = min(d(inx)); u = inx(u);
      */
        
     for( j=0,tmp = PLUS_INF, delta = 0; j < nk[k]; j++ ) {
        i = inx[INDEX(j,k,n)];
/*      for( i=0, tmp = PLUS_INF, delta = 0; i < n; i++ ) {
        if( I[i]-1 == k) {*/
        delta += x[i]*d[i];
        if( tmp > d[i] ) {
          tmp = d[i];
          u = i;
        }
      }

      /* if d(u) < 0, yu = b; else yu = 0; end  */
      if( d[u] < 0) yu = b; else yu = 0;
     
      /* delta = x(inx)'*d(inx) - yu*d(u); */
      delta -= yu*d[u];
            
      if( delta > tolabs/m && delta > tolrel*ABS(UB)/m) 
      {
         exitflag = 0;
         
         if( yu > 0 ) 
         {
           col_u = (double*)get_col(u);      

           improv = MINUS_INF;
           for( j=0; j < nk[k]; j++ ) {
             i = inx[INDEX(j,k,n)];
           
/*           for(i = 0; i < n; i++ ) {
             if( (I[i]-1 == k) && (i != u) && (x[i] > 0)) {              */
             if(x[i] > 0) {             
               
               tmp_num = x[i]*(d[i] - d[u]); 
               tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
               if( tmp_den > 0 ) {
                 if( tmp_num < tmp_den ) {
                    tmp = tmp_num*tmp_num / tmp_den;
                 } else {
                    tmp = tmp_num - 0.5 * tmp_den;
                 }
               }
               if( tmp > improv ) {
                 improv = tmp;
                 tau = MIN(1,tmp_num/tmp_den);
                 v = i;
               }
             }
           }

           tmp_num = -x_nequ[k]*d[u];
           if( tmp_num > 0 ) {
             tmp_den = x_nequ[k]*x_nequ[k]*diag_H[u];
             if( tmp_den > 0 ) {
               if( tmp_num < tmp_den ) {
                 tmp = tmp_num*tmp_num / tmp_den;
               } else {
                   tmp = tmp_num - 0.5 * tmp_den;
               }
             }
           } else {
             tmp = MINUS_INF; 
           }
           
           if( tmp > improv ) {
              tau = MIN(1,tmp_num/tmp_den);
              for( i = 0; i < n; i++ ) {             
                d[i] += x_nequ[k]*tau*col_u[i];
              }
             x[u] += tau*x_nequ[k];
             x_nequ[k] -= tau*x_nequ[k];
               
           } else {
            
             /* updating with the best line segment */
             col_v = (double*)get_col(v);
             for( i = 0; i < n; i++ ) {             
               d[i] += x[v]*tau*(col_u[i]-col_v[i]);
             }

             x[u] += tau*x[v];
             x[v] -= tau*x[v];
           }
         }
         else
         {
           improv = MINUS_INF;
           for( j=0; j < nk[k]; j++ ) {
             i = inx[INDEX(j,k,n)];
           
/*           for(i = 0; i < n; i++ ) {
             if( (I[i]-1 == k) && (x[i] > 0)) {*/
             if( x[i] > 0 && d[i] > 0) {
                
               tmp_num = x[i]*d[i]; 
               tmp_den = x[i]*x[i]*diag_H[i];
               if( tmp_den > 0 ) {
                 if( tmp_num < tmp_den ) {
                    tmp = tmp_num*tmp_num / tmp_den;
                 } else {
                    tmp = tmp_num - 0.5 * tmp_den;
                 }
               }
               if( tmp > improv ) {
                 improv = tmp;
                 tau = MIN(1,tmp_num/tmp_den);
                 v = i;
               }
             }    
           }

           /* updating with the best line segment */
           col_v = (double*)get_col(v);
           for( i = 0; i < n; i++ ) {             
             d[i] -= x[v]*tau*col_v[i];
           }

           x_nequ[k] += tau*x[v];
           x[v] -= tau*x[v];         
         }
                    
/*         for( i=0, UB = 0; i < n; i++) {
            UB += x[i]*(f[i]+d[i]);
         }
         UB = 0.5*UB;
 */
         UB = UB - improv;
      }
                   
      /* mexPrintf("t=%d,k=%d, u=%d, tau1=%f, den1=%f, num1=%f, delta=%f\n", 
             t,k,u,tau1,den1,num1,delta);*/

    }

    /* -- Computing LB --------------------------------------*/

    /*
    LB = 0.5*x'*(f-d);   
    for k=1:n,
      LB = LB + b*min(d(find(I==k)));
    end */
    
    for( i=0, UB = 0, LB=0; i < n; i++) {
       UB += x[i]*(f[i]+d[i]);
       LB += x[i]*(f[i]-d[i]);
    }
    UB = 0.5*UB;
    LB = 0.5*LB;

    for( k=0; k < m; k++ ) { 
      for( j=0,tmp = PLUS_INF; j < nk[k]; j++ ) {
        i = inx[INDEX(j,k,n)];

/*      for( j=0, tmp = PLUS_INF; j < n; j++ ) {
        if( I[j]-1 == i ) tmp = MIN(tmp, d[j]);*/
        tmp = MIN(tmp, d[i]);
      }
      if( tmp < 0) LB += b*tmp;
    }

    /* Store LB and UB */
    if( t < History_size ) {
      History[INDEX(0,t,2)] = LB;
      History[INDEX(1,t,2)] = UB;
    }
    else {
      tmp_ptr = mxCalloc((History_size+HISTORY_BUF)*2,sizeof(double));
      if( tmp_ptr == NULL ) mexErrMsgTxt("Not enough memory.");
      for( i = 0; i < History_size; i++ ) {
        tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
        tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
      }
      tmp_ptr[INDEX(0,t,2)] = LB;
      tmp_ptr[INDEX(1,t,2)] = UB;
  
      History_size += HISTORY_BUF;
      mxFree( History );
      History = tmp_ptr;
    }

    if( verb > 0 && (exitflag > 0 || (t % verb)==0 )) {
       mexPrintf("%d: UB=%.10f, LB=%.10f, UB-LB=%.10f, (UB-LB)/|UB|=%.10f \n",
        t, UB, LB, UB-LB, (UB!=0) ? (UB-LB)/ABS(UB) : 0);      
    }    

  }

  /* -- Find which stopping consition has been used -------- */
  if( UB-LB < tolabs ) exitflag = 1;
  else if(UB-LB < ABS(UB)*tolrel ) exitflag = 2;
  else exitflag = 0;

  /*----------------------------------------------------------   
    Set up outputs                                          
  ---------------------------------------------------------- */
  (*ptr_t) = t;
  (*ptr_History) = History;

  /*----------------------------------------------------------
    Clean up
  ---------------------------------------------------------- */
  mxFree( d );
  mxFree( inx );
  mxFree( nk );
  
  return( exitflag ); 

}