qpssvmlib.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  * Library of solvers for QP task required in StructSVM learning.
00008  *
00009  * Written (W) 1999-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
00010  * Copyright (C) 1999-2008 Center for Machine Perception, CTU FEL Prague
00011 --------------------------------------------------------------------
00012 Synopsis:
00013 
00014   exitflag = qpssvm_solver( &get_col, diag_H, f, b, I, x, n, tmax,
00015              tolabs, tolrel, &t, &History, verb );
00016 
00017   exitflag = qpssvm_solver( &get_col, diag_H, f, b, I, x, n, tmax,
00018              tolabs, tolrel, &QP, &QD, verb );
00019 Description:
00020 
00021  It solves the following QP task:
00022 
00023    min 0.5*x'*H*x + f'*x
00024     x
00025 
00026  subject to 
00027 
00028    sum(x(find(I==k))) <= b   for all k=1:max(I)
00029    x >= 0
00030 
00031  where I is a set of positive indices from (1 to max(I)).
00032 
00033  A precision of the found solution is given by the parameters tmax,
00034  tolabs and tolrel which define the stopping conditions:
00035 
00036  UB-LB <= tolabs      ->  exitflag = 1   Abs. tolerance.
00037  UB-LB <= UB*tolrel   ->  exitflag = 2   Relative tolerance.
00038  t >= tmax            ->  exitflag = 0   Number of iterations.
00039 
00040  UB ... Upper bound on the optimal solution, i.e., Q_P.
00041  LB ... Lower bound on the optimal solution, i.e., Q_D.
00042  t  ... Number of iterations.
00043 
00044 
00045 Inputs/Outputs:
00046 
00047  const void* (*get_col)(uint32_t) retunr pointer to i-th column of H
00048  diag_H [float64_t n x n] diagonal of H.
00049  f [float64_t n x 1] is an arbitrary vector.
00050  b [float64_t 1 x 1] scalar
00051  I [uint16_T n x 1] Indices (1..max(I)); max(I) <= n
00052  x [float64_t n x 1] solution vector (inital solution).
00053  n [uint32_t 1 x 1] dimension of H.
00054  tmax [uint32_t 1 x 1] Max number of steps.
00055  tolrel [float64_t 1 x 1] Relative tolerance.
00056  tolabs [float64_t 1 x 1] Absolute tolerance.
00057  t [uint32_t 1 x 1] Number of iterations.
00058  History [float64_t 2 x t] Value of LB and UB wrt. number of iterations.
00059  verb [int 1 x 1] if > 0 then prints info every verb-th iteation.
00060 
00061  For more info refer to TBA
00062 
00063  Modifications:
00064  01-Oct-2007, VF
00065  20-Feb-2006, VF
00066  18-feb-2006, VF
00067 
00068 -------------------------------------------------------------------- */
00069 
00070 #include <math.h>
00071 #include <stdlib.h>
00072 #include <stdio.h>
00073 #include <string.h>
00074 #include <stdint.h>
00075 #include <limits.h>
00076 
00077 #include "classifier/svm/libocas_common.h"
00078 #include "classifier/svm/qpssvmlib.h"
00079 
00080 /* --------------------------------------------------------------
00081  QPSSVM solver 
00082 
00083  Usage: exitflag = qpssvm_solver( &get_col, diag_H, f, b, I, x, n, tmax, 
00084          tolabs, tolrel, &QP, &QD, verb );   
00085 -------------------------------------------------------------- */
00086 int8_t qpssvm_solver(const void* (*get_col)(uint32_t),
00087                   float64_t *diag_H,
00088                   float64_t *f,
00089                   float64_t b,
00090                   uint16_t *I,
00091                   float64_t *x,
00092                   uint32_t n,
00093                   uint32_t tmax,
00094                   float64_t tolabs,
00095                   float64_t tolrel,
00096                   float64_t *QP,
00097                   float64_t *QD,
00098                   uint32_t verb)
00099 {
00100   float64_t *x_nequ;
00101   float64_t *d;
00102   float64_t *col_u, *col_v;
00103   float64_t LB;
00104   float64_t UB;
00105   float64_t tmp;
00106   float64_t improv;
00107   float64_t tmp_num;
00108   float64_t tmp_den=0;
00109   float64_t tau=0;
00110   float64_t delta;
00111   float64_t yu;
00112   uint32_t *inx;
00113   uint32_t *nk;
00114   uint32_t m;
00115   uint32_t t;
00116   uint32_t u=0;
00117   uint32_t v=0;
00118   uint32_t k;
00119   uint32_t i, j;
00120   int8_t exitflag;
00121 
00122 
00123   /* ------------------------------------------------------------ 
00124     Initialization                                               
00125   ------------------------------------------------------------ */
00126 
00127   x_nequ=NULL;
00128   inx=NULL;
00129   nk=NULL;
00130   d=NULL;
00131 
00132   /* count cumber of constraints */
00133   for( i=0, m=0; i < n; i++ ) m = MAX(m,I[i]);
00134 
00135   /* alloc and initialize x_nequ */
00136   x_nequ = (float64_t*) OCAS_CALLOC(m, sizeof(float64_t));
00137   if( x_nequ == NULL )
00138   {
00139       exitflag=-2;
00140       goto cleanup;
00141   }
00142 
00143   /* alloc Inx */
00144   inx = (uint32_t*) OCAS_CALLOC(m*n, sizeof(uint32_t));
00145   if( inx == NULL )
00146   {
00147       exitflag=-2;
00148       goto cleanup;
00149   }
00150 
00151   nk = (uint32_t*) OCAS_CALLOC(m, sizeof(uint32_t));
00152   if( nk == NULL )
00153   {
00154       exitflag=-2;
00155       goto cleanup;
00156   }
00157 
00158   for( i=0; i < m; i++ ) x_nequ[i] = b;
00159   for( i=0; i < n; i++ ) {
00160      k = I[i]-1;
00161      x_nequ[k] -= x[i];
00162      inx[INDEX2(nk[k],k,n)] = i;
00163      nk[k]++;
00164   }
00165     
00166   /* alloc d [n x 1] */
00167   d = (float64_t*) OCAS_CALLOC(n, sizeof(float64_t));
00168   if( d == NULL )
00169   {
00170       exitflag=-2;
00171       goto cleanup;
00172   }
00173  
00174   /* d = H*x + f; */
00175   for( i=0; i < n; i++ ) {
00176     if( x[i] > 0 ) {
00177       col_u = (float64_t*)get_col(i);
00178       for( j=0; j < n; j++ ) {
00179           d[j] += col_u[j]*x[i];
00180       }
00181     }
00182   }
00183   for( i=0; i < n; i++ ) d[i] += f[i];
00184   
00185   /* UB = 0.5*x'*(f+d); */
00186   /* LB = 0.5*x'*(f-d); */
00187   for( i=0, UB = 0, LB=0; i < n; i++) {
00188     UB += x[i]*(f[i]+d[i]);
00189     LB += x[i]*(f[i]-d[i]);
00190   }
00191   UB = 0.5*UB;
00192   LB = 0.5*LB;
00193 
00194   /*
00195   for k=1:m,
00196     tmp = min(d(find(I==k)));
00197     if tmp < 0, LB = LB + b*tmp; end
00198   end
00199   */
00200   
00201   for( i=0; i < m; i++ ) {
00202     for( j=0, tmp = OCAS_PLUS_INF; j < nk[i]; j++ ) {
00203       tmp = MIN(tmp, d[inx[INDEX2(j,i,n)]]);
00204     }
00205     if( tmp < 0) LB += b*tmp;
00206   }
00207   
00208   exitflag = 0;
00209   t = 0;
00210 
00211   /* -- Main loop ---------------------------------------- */
00212   while( (exitflag == 0) && (t < tmax)) 
00213   {
00214     t++;
00215 
00216     exitflag = 1;
00217     for( k=0; k < m; k++ ) 
00218     {       
00219       /*
00220       inx = find(I==k);
00221       [tmp,u] = min(d(inx)); u = inx(u);
00222       */
00223         
00224      for( j=0, tmp = OCAS_PLUS_INF, delta = 0; j < nk[k]; j++ ) {
00225         i = inx[INDEX2(j,k,n)];
00226         delta += x[i]*d[i];
00227         if( tmp > d[i] ) {
00228           tmp = d[i];
00229           u = i;
00230         }
00231       }
00232 
00233       /* if d(u) < 0, yu = b; else yu = 0; end  */
00234       if( d[u] < 0) yu = b; else yu = 0;
00235      
00236       /* delta = x(inx)'*d(inx) - yu*d(u); */
00237       delta -= yu*d[u];
00238             
00239       if( delta > tolabs/m && delta > tolrel*ABS(UB)/m) 
00240       {
00241          exitflag = 0;
00242          
00243          if( yu > 0 ) 
00244          {
00245            col_u = (float64_t*)get_col(u);      
00246 
00247            improv = -OCAS_PLUS_INF;
00248            for( j=0; j < nk[k]; j++ ) {
00249              i = inx[INDEX2(j,k,n)];
00250            
00251 /*           for(i = 0; i < n; i++ ) {
00252              if( (I[i]-1 == k) && (i != u) && (x[i] > 0)) {              */
00253              if(x[i] > 0) {             
00254                
00255                tmp_num = x[i]*(d[i] - d[u]); 
00256                tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
00257                if( tmp_den > 0 ) {
00258                  if( tmp_num < tmp_den ) {
00259                     tmp = tmp_num*tmp_num / tmp_den;
00260                  } else {
00261                     tmp = tmp_num - 0.5 * tmp_den;
00262                  }
00263                }
00264                if( tmp > improv ) {
00265                  improv = tmp;
00266                  tau = MIN(1,tmp_num/tmp_den);
00267                  v = i;
00268                }
00269              }
00270            }
00271 
00272            tmp_num = -x_nequ[k]*d[u];
00273            if( tmp_num > 0 ) {
00274              tmp_den = x_nequ[k]*x_nequ[k]*diag_H[u];
00275              if( tmp_den > 0 ) {
00276                if( tmp_num < tmp_den ) {
00277                  tmp = tmp_num*tmp_num / tmp_den;
00278                } else {
00279                    tmp = tmp_num - 0.5 * tmp_den;
00280                }
00281              }
00282            } else {
00283              tmp = -OCAS_PLUS_INF; 
00284            }
00285            
00286            if( tmp > improv ) {
00287               tau = MIN(1,tmp_num/tmp_den);
00288               for( i = 0; i < n; i++ ) {             
00289                 d[i] += x_nequ[k]*tau*col_u[i];
00290               }
00291              x[u] += tau*x_nequ[k];
00292              x_nequ[k] -= tau*x_nequ[k];
00293                
00294            } else {
00295             
00296              /* updating with the best line segment */
00297              col_v = (float64_t*)get_col(v);
00298              for( i = 0; i < n; i++ ) {             
00299                d[i] += x[v]*tau*(col_u[i]-col_v[i]);
00300              }
00301 
00302              x[u] += tau*x[v];
00303              x[v] -= tau*x[v];
00304            }
00305          }
00306          else
00307          {
00308            improv = -OCAS_PLUS_INF;
00309            for( j=0; j < nk[k]; j++ ) {
00310              i = inx[INDEX2(j,k,n)];
00311            
00312 /*           for(i = 0; i < n; i++ ) {
00313              if( (I[i]-1 == k) && (x[i] > 0)) {*/
00314              if( x[i] > 0 && d[i] > 0) {
00315                 
00316                tmp_num = x[i]*d[i]; 
00317                tmp_den = x[i]*x[i]*diag_H[i];
00318                if( tmp_den > 0 ) {
00319                  if( tmp_num < tmp_den ) {
00320                     tmp = tmp_num*tmp_num / tmp_den;
00321                  } else {
00322                     tmp = tmp_num - 0.5 * tmp_den;
00323                  }
00324                }
00325                if( tmp > improv ) {
00326                  improv = tmp;
00327                  tau = MIN(1,tmp_num/tmp_den);
00328                  v = i;
00329                }
00330              }    
00331            }
00332 
00333            /* updating with the best line segment */
00334            col_v = (float64_t*)get_col(v);
00335            for( i = 0; i < n; i++ ) {             
00336              d[i] -= x[v]*tau*col_v[i];
00337            }
00338 
00339            x_nequ[k] += tau*x[v];
00340            x[v] -= tau*x[v];         
00341          }
00342 
00343          UB = UB - improv;
00344       }
00345                    
00346     }
00347 
00348     /* -- Computing LB --------------------------------------*/
00349 
00350     /*
00351     LB = 0.5*x'*(f-d);   
00352     for k=1:n,
00353       LB = LB + b*min(d(find(I==k)));
00354     end */
00355     
00356     for( i=0, UB = 0, LB=0; i < n; i++) {
00357        UB += x[i]*(f[i]+d[i]);
00358        LB += x[i]*(f[i]-d[i]);
00359     }
00360     UB = 0.5*UB;
00361     LB = 0.5*LB;
00362 
00363     for( k=0; k < m; k++ ) { 
00364       for( j=0,tmp = OCAS_PLUS_INF; j < nk[k]; j++ ) {
00365         i = inx[INDEX2(j,k,n)];
00366 
00367         tmp = MIN(tmp, d[i]);
00368       }
00369       if( tmp < 0) LB += b*tmp;
00370     }
00371 
00372     if( verb > 0 && (exitflag > 0 || (t % verb)==0 ))
00373     {
00374         float64_t gap=(UB!=0) ? (UB-LB)/ABS(UB) : 0;
00375         SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(tolrel), 6);
00376     }
00377 
00378   }
00379 
00380   /* -- Find which stopping consition has been used -------- */
00381   if( UB-LB < tolabs ) exitflag = 1;
00382   else if(UB-LB < ABS(UB)*tolrel ) exitflag = 2;
00383   else exitflag = 0;
00384 
00385   /*----------------------------------------------------------   
00386     Set up outputs                                          
00387   ---------------------------------------------------------- */
00388   *QP = UB;
00389   *QD = LB;
00390 
00391   /*----------------------------------------------------------
00392     Clean up
00393   ---------------------------------------------------------- */
00394 cleanup:
00395   OCAS_FREE( d );
00396   OCAS_FREE( inx );
00397   OCAS_FREE( nk );
00398   OCAS_FREE( x_nequ );  
00399   
00400   return( exitflag ); 
00401 
00402 }
00403 

SHOGUN Machine Learning Toolbox - Documentation