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 namespace shogun
00081 {
00082 /* --------------------------------------------------------------
00083  QPSSVM solver 
00084 
00085  Usage: exitflag = qpssvm_solver( &get_col, diag_H, f, b, I, x, n, tmax, 
00086          tolabs, tolrel, &QP, &QD, verb );   
00087 -------------------------------------------------------------- */
00088 int8_t qpssvm_solver(const void* (*get_col)(uint32_t),
00089                   float64_t *diag_H,
00090                   float64_t *f,
00091                   float64_t b,
00092                   uint16_t *I,
00093                   float64_t *x,
00094                   uint32_t n,
00095                   uint32_t tmax,
00096                   float64_t tolabs,
00097                   float64_t tolrel,
00098                   float64_t *QP,
00099                   float64_t *QD,
00100                   uint32_t verb)
00101 {
00102   float64_t *x_nequ;
00103   float64_t *d;
00104   float64_t *col_u, *col_v;
00105   float64_t LB;
00106   float64_t UB;
00107   float64_t tmp;
00108   float64_t improv;
00109   float64_t tmp_num;
00110   float64_t tmp_den=0;
00111   float64_t tau=0;
00112   float64_t delta;
00113   float64_t yu;
00114   uint32_t *inx;
00115   uint32_t *nk;
00116   uint32_t m;
00117   uint32_t t;
00118   uint32_t u=0;
00119   uint32_t v=0;
00120   uint32_t k;
00121   uint32_t i, j;
00122   int8_t exitflag;
00123 
00124 
00125   /* ------------------------------------------------------------ 
00126     Initialization                                               
00127   ------------------------------------------------------------ */
00128 
00129   x_nequ=NULL;
00130   inx=NULL;
00131   nk=NULL;
00132   d=NULL;
00133 
00134   /* count cumber of constraints */
00135   for( i=0, m=0; i < n; i++ ) m = CMath::max(m,(uint32_t) I[i]);
00136 
00137   /* alloc and initialize x_nequ */
00138   x_nequ = (float64_t*) OCAS_CALLOC(m, sizeof(float64_t));
00139   if( x_nequ == NULL )
00140   {
00141       exitflag=-2;
00142       goto cleanup;
00143   }
00144 
00145   /* alloc Inx */
00146   inx = (uint32_t*) OCAS_CALLOC(m*n, sizeof(uint32_t));
00147   if( inx == NULL )
00148   {
00149       exitflag=-2;
00150       goto cleanup;
00151   }
00152 
00153   nk = (uint32_t*) OCAS_CALLOC(m, sizeof(uint32_t));
00154   if( nk == NULL )
00155   {
00156       exitflag=-2;
00157       goto cleanup;
00158   }
00159 
00160   for( i=0; i < m; i++ ) x_nequ[i] = b;
00161   for( i=0; i < n; i++ ) {
00162      k = I[i]-1;
00163      x_nequ[k] -= x[i];
00164      inx[INDEX2(nk[k],k,n)] = i;
00165      nk[k]++;
00166   }
00167     
00168   /* alloc d [n x 1] */
00169   d = (float64_t*) OCAS_CALLOC(n, sizeof(float64_t));
00170   if( d == NULL )
00171   {
00172       exitflag=-2;
00173       goto cleanup;
00174   }
00175  
00176   /* d = H*x + f; */
00177   for( i=0; i < n; i++ ) {
00178     if( x[i] > 0 ) {
00179       col_u = (float64_t*)get_col(i);
00180       for( j=0; j < n; j++ ) {
00181           d[j] += col_u[j]*x[i];
00182       }
00183     }
00184   }
00185   for( i=0; i < n; i++ ) d[i] += f[i];
00186   
00187   /* UB = 0.5*x'*(f+d); */
00188   /* LB = 0.5*x'*(f-d); */
00189   for( i=0, UB = 0, LB=0; i < n; i++) {
00190     UB += x[i]*(f[i]+d[i]);
00191     LB += x[i]*(f[i]-d[i]);
00192   }
00193   UB = 0.5*UB;
00194   LB = 0.5*LB;
00195 
00196   /*
00197   for k=1:m,
00198     tmp = min(d(find(I==k)));
00199     if tmp < 0, LB = LB + b*tmp; end
00200   end
00201   */
00202   
00203   for( i=0; i < m; i++ ) {
00204     for( j=0, tmp = OCAS_PLUS_INF; j < nk[i]; j++ ) {
00205       tmp = CMath::min(tmp, d[inx[INDEX2(j,i,n)]]);
00206     }
00207     if( tmp < 0) LB += b*tmp;
00208   }
00209   
00210   exitflag = 0;
00211   t = 0;
00212 
00213   /* -- Main loop ---------------------------------------- */
00214   while( (exitflag == 0) && (t < tmax)) 
00215   {
00216     t++;
00217 
00218     exitflag = 1;
00219     for( k=0; k < m; k++ ) 
00220     {       
00221       /*
00222       inx = find(I==k);
00223       [tmp,u] = min(d(inx)); u = inx(u);
00224       */
00225         
00226      for( j=0, tmp = OCAS_PLUS_INF, delta = 0; j < nk[k]; j++ ) {
00227         i = inx[INDEX2(j,k,n)];
00228         delta += x[i]*d[i];
00229         if( tmp > d[i] ) {
00230           tmp = d[i];
00231           u = i;
00232         }
00233       }
00234 
00235       /* if d(u) < 0, yu = b; else yu = 0; end  */
00236       if( d[u] < 0) yu = b; else yu = 0;
00237      
00238       /* delta = x(inx)'*d(inx) - yu*d(u); */
00239       delta -= yu*d[u];
00240             
00241       if( delta > tolabs/m && delta > tolrel*CMath::abs(UB)/m) 
00242       {
00243          exitflag = 0;
00244          
00245          if( yu > 0 ) 
00246          {
00247            col_u = (float64_t*)get_col(u);      
00248 
00249            improv = -OCAS_PLUS_INF;
00250            for( j=0; j < nk[k]; j++ ) {
00251              i = inx[INDEX2(j,k,n)];
00252            
00253 /*           for(i = 0; i < n; i++ ) {
00254              if( (I[i]-1 == k) && (i != u) && (x[i] > 0)) {              */
00255              if(x[i] > 0) {             
00256                
00257                tmp_num = x[i]*(d[i] - d[u]); 
00258                tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
00259                if( tmp_den > 0 ) {
00260                  if( tmp_num < tmp_den ) {
00261                     tmp = tmp_num*tmp_num / tmp_den;
00262                  } else {
00263                     tmp = tmp_num - 0.5 * tmp_den;
00264                  }
00265                }
00266                if( tmp > improv ) {
00267                  improv = tmp;
00268                  tau = CMath::min(1.0,tmp_num/tmp_den);
00269                  v = i;
00270                }
00271              }
00272            }
00273 
00274            tmp_num = -x_nequ[k]*d[u];
00275            if( tmp_num > 0 ) {
00276              tmp_den = x_nequ[k]*x_nequ[k]*diag_H[u];
00277              if( tmp_den > 0 ) {
00278                if( tmp_num < tmp_den ) {
00279                  tmp = tmp_num*tmp_num / tmp_den;
00280                } else {
00281                    tmp = tmp_num - 0.5 * tmp_den;
00282                }
00283              }
00284            } else {
00285              tmp = -OCAS_PLUS_INF; 
00286            }
00287            
00288            if( tmp > improv ) {
00289               tau = CMath::min(1.0,tmp_num/tmp_den);
00290               for( i = 0; i < n; i++ ) {             
00291                 d[i] += x_nequ[k]*tau*col_u[i];
00292               }
00293              x[u] += tau*x_nequ[k];
00294              x_nequ[k] -= tau*x_nequ[k];
00295                
00296            } else {
00297             
00298              /* updating with the best line segment */
00299              col_v = (float64_t*)get_col(v);
00300              for( i = 0; i < n; i++ ) {             
00301                d[i] += x[v]*tau*(col_u[i]-col_v[i]);
00302              }
00303 
00304              x[u] += tau*x[v];
00305              x[v] -= tau*x[v];
00306            }
00307          }
00308          else
00309          {
00310            improv = -OCAS_PLUS_INF;
00311            for( j=0; j < nk[k]; j++ ) {
00312              i = inx[INDEX2(j,k,n)];
00313            
00314 /*           for(i = 0; i < n; i++ ) {
00315              if( (I[i]-1 == k) && (x[i] > 0)) {*/
00316              if( x[i] > 0 && d[i] > 0) {
00317                 
00318                tmp_num = x[i]*d[i]; 
00319                tmp_den = x[i]*x[i]*diag_H[i];
00320                if( tmp_den > 0 ) {
00321                  if( tmp_num < tmp_den ) {
00322                     tmp = tmp_num*tmp_num / tmp_den;
00323                  } else {
00324                     tmp = tmp_num - 0.5 * tmp_den;
00325                  }
00326                }
00327                if( tmp > improv ) {
00328                  improv = tmp;
00329                  tau = CMath::min(1.0,tmp_num/tmp_den);
00330                  v = i;
00331                }
00332              }    
00333            }
00334 
00335            /* updating with the best line segment */
00336            col_v = (float64_t*)get_col(v);
00337            for( i = 0; i < n; i++ ) {             
00338              d[i] -= x[v]*tau*col_v[i];
00339            }
00340 
00341            x_nequ[k] += tau*x[v];
00342            x[v] -= tau*x[v];         
00343          }
00344 
00345          UB = UB - improv;
00346       }
00347                    
00348     }
00349 
00350     /* -- Computing LB --------------------------------------*/
00351 
00352     /*
00353     LB = 0.5*x'*(f-d);   
00354     for k=1:n,
00355       LB = LB + b*min(d(find(I==k)));
00356     end */
00357     
00358     for( i=0, UB = 0, LB=0; i < n; i++) {
00359        UB += x[i]*(f[i]+d[i]);
00360        LB += x[i]*(f[i]-d[i]);
00361     }
00362     UB = 0.5*UB;
00363     LB = 0.5*LB;
00364 
00365     for( k=0; k < m; k++ ) { 
00366       for( j=0,tmp = OCAS_PLUS_INF; j < nk[k]; j++ ) {
00367         i = inx[INDEX2(j,k,n)];
00368 
00369         tmp = CMath::min(tmp, d[i]);
00370       }
00371       if( tmp < 0) LB += b*tmp;
00372     }
00373 
00374     if( verb > 0 && (exitflag > 0 || (t % verb)==0 ))
00375     {
00376         float64_t gap=(UB!=0) ? (UB-LB)/CMath::abs(UB) : 0;
00377         SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(tolrel), 6);
00378     }
00379 
00380   }
00381 
00382   /* -- Find which stopping consition has been used -------- */
00383   if( UB-LB < tolabs ) exitflag = 1;
00384   else if(UB-LB < CMath::abs(UB)*tolrel ) exitflag = 2;
00385   else exitflag = 0;
00386 
00387   /*----------------------------------------------------------   
00388     Set up outputs                                          
00389   ---------------------------------------------------------- */
00390   *QP = UB;
00391   *QD = LB;
00392 
00393   /*----------------------------------------------------------
00394     Clean up
00395   ---------------------------------------------------------- */
00396 cleanup:
00397   OCAS_FREE( d );
00398   OCAS_FREE( inx );
00399   OCAS_FREE( nk );
00400   OCAS_FREE( x_nequ );  
00401   
00402   return( exitflag ); 
00403 
00404 }
00405 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation