libqp_splx.cpp

Go to the documentation of this file.
00001 /*-----------------------------------------------------------------------
00002  * libqp_splx.c: solver for Quadratic Programming task with 
00003  * simplex constraints.
00004  *
00005  * DESCRIPTION
00006  *  The library provides function which solves the following instance of
00007  *  a convex Quadratic Programmin task:
00008  *  
00009  *   min QP(x):= 0.5*x'*H*x + f'*x  
00010  *    x
00011  *
00012  * subject to:   
00013  *   sum_{i in I_k} x[i] == b[k]  for all k such that S[k] == 0 
00014  *   sum_{i in I_k} x[i] <= b[k]  for all k such that S[k] == 1
00015  *                             x(i) >= 0 for all i=1:n
00016  *   
00017  *  where I_k = { i | I[i] == k}, k={1,...,m}.
00018  *
00019  * A precision of the found solution is controled by the input argumens
00020  * MaxIter, TolAbs, QP_TH and MaxIter which define the stopping conditions:
00021  * 
00022  *  nIter >= MaxIter     ->  exitflag = 0   Number of iterations
00023  *  QP-QD <= TolAbs      ->  exitflag = 1   Abs. tolerance (duality gap)
00024  *  QP-QD <= QP*TolRel   ->  exitflag = 2   Relative tolerance
00025  *  QP <= QP_TH          ->  exitflag = 3   Threshold on objective value
00026  *
00027  * where QP and QD are primal respectively dual objective values.
00028  *
00029  * INPUT ARGUMENTS
00030  *  get_col   function which returns pointer to the i-th column of H.
00031  *  diag_H [float64_t n x 1] vector containing values on the diagonal of H.
00032  *  f [float64_t n x 1] vector.
00033  *  b [float64_t n x 1] vector of positive numbers.
00034  *  I [uint16_T n x 1] vector containing numbers 1...m. 
00035  *  S [uint8_T n x 1] vector containing numbers 0 and 1.
00036  *  x [float64_t n x 1] solution vector; must be feasible.
00037  *  n [uint32_t 1 x 1] dimension of H.
00038  *  MaxIter [uint32_t 1 x 1] max number of iterations.
00039  *  TolAbs [float64_t 1 x 1] Absolute tolerance.
00040  *  TolRel [float64_t 1 x 1] Relative tolerance.
00041  *  QP_TH  [float64_t 1 x 1] Threshold on the primal value.
00042  *  print_state  print function; if == NULL it is not called.
00043  *
00044  * RETURN VALUE
00045  *  structure [libqp_state_T] 
00046  *  .QP [1 x 1] Primal objective value.
00047  *  .QD [1 x 1] Dual objective value.
00048  *  .nIter [1 x 1] Number of iterations.
00049  *  .exitflag [1 x 1] Indicates which stopping condition was used:
00050  *    -1  ... Not enough memory.
00051  *     0  ... Maximal number of iteations reached: nIter >= MaxIter.
00052  *     1  ... Relarive tolerance reached: QP-QD <= abs(QP)*TolRel
00053  *     2  ... Absolute tolerance reached: QP-QD <= TolAbs
00054  *     3  ... Objective value reached threshold: QP <= QP_TH.
00055  *
00056  * REFERENCE
00057  *  The algorithm is described in:
00058  *  V. Franc, V. Hlavac. A Novel Algorithm for Learning Support Vector Machines
00059  *   with Structured Output Spaces. Research Report K333 22/06, CTU-CMP-2006-04. 
00060  *   May, 2006. ftp://cmp.felk.cvut.cz/pub/cmp/articles/franc/Franc-TR-2006-04.ps
00061  *
00062  * Copyright (C) 2006-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
00063  * Center for Machine Perception, CTU FEL Prague
00064  *
00065  * This program is free software; you can redistribute it and/or
00066  * modify it under the terms of the GNU General Public 
00067  * License as published by the Free Software Foundation; 
00068  * Version 3, 29 June 2007
00069  *-------------------------------------------------------------------- */
00070 
00071 #include <math.h>
00072 #include <stdlib.h>
00073 #include <stdio.h>
00074 #include <string.h>
00075 #include <stdint.h>
00076 #include <limits.h>
00077 
00078 #include <shogun/lib/common.h>
00079 #include <shogun/classifier/svm/libqp.h>
00080 namespace shogun
00081 {
00082 
00083 libqp_state_T libqp_splx_solver(const float64_t* (*get_col)(uint32_t),
00084                   float64_t *diag_H,
00085                   float64_t *f,
00086                   float64_t *b,
00087                   uint32_t *I,
00088                   uint8_t *S,
00089                   float64_t *x,
00090                   uint32_t n,
00091                   uint32_t MaxIter,
00092                   float64_t TolAbs,
00093                   float64_t TolRel,
00094                   float64_t QP_TH,
00095                   void (*print_state)(libqp_state_T state))
00096 {
00097   float64_t *d;
00098   float64_t *col_u, *col_v;
00099   float64_t *x_neq;
00100   float64_t tmp;
00101   float64_t improv;
00102   float64_t tmp_num;
00103   float64_t tmp_den=0;
00104   float64_t tau=0;
00105   float64_t delta;
00106   uint32_t *inx;
00107   uint32_t *nk;
00108   uint32_t m;
00109   int32_t u=0;
00110   int32_t v=0;
00111   uint32_t k;
00112   uint32_t i, j;
00113   libqp_state_T state;
00114 
00115   
00116   /* ------------------------------------------------------------ 
00117     Initialization                                               
00118   ------------------------------------------------------------ */
00119   state.nIter = 0;
00120   state.QP = LIBQP_PLUS_INF;
00121   state.QD = -LIBQP_PLUS_INF;
00122   state.exitflag = 100;
00123 
00124   inx=NULL;
00125   nk=NULL;
00126   d=NULL;
00127   x_neq = NULL;
00128 
00129   /* count number of constraints */
00130   for( i=0, m=0; i < n; i++ ) 
00131     m = LIBQP_MAX(m,I[i]);
00132 
00133   /* auxciliary variables for tranforming equalities to inequalities */
00134   x_neq = (float64_t*) LIBQP_CALLOC(m, sizeof(float64_t));
00135   if( x_neq == NULL )
00136   {
00137       state.exitflag=-1;
00138       goto cleanup;
00139   }
00140 
00141   /* inx is translation table between variable index i and its contraint */
00142   inx = (uint32_t*) LIBQP_CALLOC(m*n, sizeof(uint32_t));
00143   if( inx == NULL )
00144   {
00145       state.exitflag=-1;
00146       goto cleanup;
00147   }
00148 
00149   /* nk is the number of variables coupled by i-th linear constraint */
00150   nk = (uint32_t*) LIBQP_CALLOC(m, sizeof(uint32_t));
00151   if( nk == NULL )
00152   {
00153       state.exitflag=-1;
00154       goto cleanup;
00155   }
00156 
00157   /* setup auxciliary variables */
00158   for( i=0; i < m; i++ ) 
00159     x_neq[i] = b[i];
00160 
00161 
00162   /* create inx and nk */
00163   for( i=0; i < n; i++ ) {
00164      k = I[i]-1;
00165      inx[LIBQP_INDEX(nk[k],k,n)] = i;
00166      nk[k]++;     
00167 
00168      if(S[k] != 0) 
00169        x_neq[k] -= x[i];
00170   }
00171     
00172   /* d = H*x + f is gradient*/
00173   d = (float64_t*) LIBQP_CALLOC(n, sizeof(float64_t));
00174   if( d == NULL )
00175   {
00176       state.exitflag=-1;
00177       goto cleanup;
00178   }
00179  
00180   /* compute gradient */
00181   for( i=0; i < n; i++ ) 
00182   {
00183     d[i] += f[i];
00184     if( x[i] > 0 ) {
00185       col_u = (float64_t*)get_col(i);      
00186       for( j=0; j < n; j++ ) {
00187           d[j] += col_u[j]*x[i];
00188       }
00189     }
00190   }
00191   
00192   /* compute state.QP = 0.5*x'*(f+d);
00193              state.QD = 0.5*x'*(f-d); */
00194   for( i=0, state.QP = 0, state.QD=0; i < n; i++) 
00195   {
00196     state.QP += x[i]*(f[i]+d[i]);
00197     state.QD += x[i]*(f[i]-d[i]);
00198   }
00199   state.QP = 0.5*state.QP;
00200   state.QD = 0.5*state.QD;
00201   
00202   for( i=0; i < m; i++ ) 
00203   {
00204     for( j=0, tmp = LIBQP_PLUS_INF; j < nk[i]; j++ ) 
00205       tmp = LIBQP_MIN(tmp, d[inx[LIBQP_INDEX(j,i,n)]]);
00206 
00207     if(S[i] == 0) 
00208       state.QD += b[i]*tmp;
00209     else
00210       state.QD += b[i]*LIBQP_MIN(tmp,0);
00211   }
00212   
00213   /* print initial state */
00214   if( print_state != NULL) 
00215     print_state( state );
00216 
00217   /* ------------------------------------------------------------ 
00218     Main optimization loop 
00219   ------------------------------------------------------------ */
00220   while( state.exitflag == 100 ) 
00221   {
00222     state.nIter ++;
00223 
00224     /* go over blocks of variables coupled by lin. constraint */
00225     for( k=0; k < m; k++ ) 
00226     {       
00227         
00228       /* compute u = argmin_{i in I_k} d[i] 
00229              delta =  sum_{i in I_k} x[i]*d[i] - b*min_{i in I_k} */
00230       for( j=0, tmp = LIBQP_PLUS_INF, delta = 0; j < nk[k]; j++ ) 
00231       {
00232         i = inx[LIBQP_INDEX(j,k,n)];
00233         delta += x[i]*d[i];
00234         if( tmp > d[i] ) {
00235           tmp = d[i];
00236           u = i;
00237         }
00238       }
00239 
00240       if(S[k] != 0 && d[u] > 0) 
00241         u = -1;
00242       else
00243         delta -= b[k]*d[u];
00244             
00245       /* if satisfied then k-th block of variables needs update */
00246       if( delta > TolAbs/m && delta > TolRel*LIBQP_ABS(state.QP)/m) 
00247       {         
00248         /* for fixed u select v = argmax_{i in I_k} Improvement(i) */
00249         if( u != -1 ) 
00250         {
00251           col_u = (float64_t*)get_col(u);
00252           improv = -LIBQP_PLUS_INF;
00253           for( j=0; j < nk[k]; j++ ) 
00254           {
00255             i = inx[LIBQP_INDEX(j,k,n)];
00256            
00257             if(x[i] > 0 && i != uint32_t(u)) 
00258             {
00259               tmp_num = x[i]*(d[i] - d[u]); 
00260               tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
00261               if( tmp_den > 0 ) 
00262               {
00263                 if( tmp_num < tmp_den ) 
00264                   tmp = tmp_num*tmp_num / tmp_den;
00265                 else 
00266                   tmp = tmp_num - 0.5 * tmp_den;
00267                  
00268                 if( tmp > improv ) 
00269                 { 
00270                   improv = tmp;
00271                   tau = LIBQP_MIN(1,tmp_num/tmp_den);
00272                   v = i;
00273                 } 
00274               }
00275             }
00276           }
00277 
00278           /* check if virtual variable can be for updated */
00279           if(x_neq[k] > 0 && S[k] != 0) 
00280           {
00281             tmp_num = -x_neq[k]*d[u]; 
00282             tmp_den = x_neq[k]*x_neq[k]*diag_H[u];
00283             if( tmp_den > 0 ) 
00284             {
00285               if( tmp_num < tmp_den ) 
00286                 tmp = tmp_num*tmp_num / tmp_den;
00287               else 
00288                 tmp = tmp_num - 0.5 * tmp_den;
00289                  
00290               if( tmp > improv ) 
00291               { 
00292                 improv = tmp;
00293                 tau = LIBQP_MIN(1,tmp_num/tmp_den);
00294                 v = -1;
00295               } 
00296             }
00297           }
00298 
00299           /* minimize objective w.r.t variable u and v */
00300           if(v != -1)
00301           {
00302             tmp = x[v]*tau;
00303             x[u] += tmp;
00304             x[v] -= tmp;
00305 
00306             /* update d = H*x + f */
00307             col_v = (float64_t*)get_col(v);
00308             for(i = 0; i < n; i++ )              
00309               d[i] += tmp*(col_u[i]-col_v[i]);
00310           }
00311           else
00312           {
00313             tmp = x_neq[k]*tau;
00314             x[u] += tmp;
00315             x_neq[k] -= tmp;
00316 
00317             /* update d = H*x + f */
00318             for(i = 0; i < n; i++ )              
00319               d[i] += tmp*col_u[i];
00320           }
00321         }
00322         else
00323         {
00324           improv = -LIBQP_PLUS_INF;
00325           for( j=0; j < nk[k]; j++ ) 
00326           {
00327             i = inx[LIBQP_INDEX(j,k,n)];
00328            
00329             if(x[i] > 0) 
00330             {
00331               tmp_num = x[i]*d[i]; 
00332               tmp_den = x[i]*x[i]*diag_H[i];
00333               if( tmp_den > 0 ) 
00334               {
00335                 if( tmp_num < tmp_den ) 
00336                   tmp = tmp_num*tmp_num / tmp_den;
00337                 else 
00338                   tmp = tmp_num - 0.5 * tmp_den;
00339                  
00340                 if( tmp > improv ) 
00341                 { 
00342                   improv = tmp;
00343                   tau = LIBQP_MIN(1,tmp_num/tmp_den);
00344                   v = i;
00345                 } 
00346               }
00347             }
00348           }
00349 
00350           tmp = x[v]*tau;
00351           x_neq[k] += tmp;
00352           x[v] -= tmp;
00353 
00354           /* update d = H*x + f */
00355           col_v = (float64_t*)get_col(v);
00356           for(i = 0; i < n; i++ )              
00357             d[i] -= tmp*col_v[i];
00358         }
00359 
00360         /* update objective value */
00361         state.QP = state.QP - improv;
00362       }
00363     }
00364     
00365     /* Compute primal and dual objectives */
00366     for( i=0, state.QP = 0, state.QD=0; i < n; i++) 
00367     {
00368        state.QP += x[i]*(f[i]+d[i]);
00369        state.QD += x[i]*(f[i]-d[i]);
00370     }
00371     state.QP = 0.5*state.QP;
00372     state.QD = 0.5*state.QD;
00373 
00374     for( k=0; k < m; k++ ) 
00375     { 
00376       for( j=0,tmp = LIBQP_PLUS_INF; j < nk[k]; j++ ) {
00377         i = inx[LIBQP_INDEX(j,k,n)];
00378         tmp = LIBQP_MIN(tmp, d[i]);
00379       }
00380       
00381       if(S[k] == 0) 
00382         state.QD += b[k]*tmp;
00383       else
00384         state.QD += b[k]*LIBQP_MIN(tmp,0);
00385     }
00386 
00387     /* print state */
00388     if( print_state != NULL) 
00389       print_state( state );
00390 
00391     /* check stopping conditions */
00392     if(state.QP-state.QD <= LIBQP_ABS(state.QP)*TolRel ) state.exitflag = 1;
00393     else if( state.QP-state.QD <= TolAbs ) state.exitflag = 2;
00394     else if( state.QP <= QP_TH ) state.exitflag = 3;
00395     else if( state.nIter >= MaxIter) state.exitflag = 0;
00396   }
00397 
00398   /*----------------------------------------------------------
00399     Clean up
00400   ---------------------------------------------------------- */
00401 cleanup:
00402   LIBQP_FREE( d );
00403   LIBQP_FREE( inx );
00404   LIBQP_FREE( nk );
00405   LIBQP_FREE( x_neq );
00406   
00407   return( state ); 
00408 }
00409 }
00410 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation