gpm.cpp

Go to the documentation of this file.
00001 /******************************************************************************
00002  ***        GPDT - Gradient Projection Decomposition Technique              ***
00003  ******************************************************************************
00004  ***                                                                        ***
00005  *** GPDT is a C++ software designed to train large-scale Support Vector    ***
00006  *** Machines for binary classification in both scalar and distributed      ***
00007  *** memory parallel environments. It uses the Joachims' problem            ***
00008  *** decomposition technique to split the whole quadratic programming (QP)  ***
00009  *** problem into a sequence of smaller QP subproblems, each one being      ***
00010  *** solved by a suitable gradient projection method (GPM). The presently   ***
00011  *** implemented GPMs are the Generalized Variable Projection Method        ***
00012  *** GVPM (T. Serafini, G. Zanghirati, L. Zanni, "Gradient Projection       ***
00013  *** Methods for Quadratic Programs and Applications in Training Support    ***
00014  *** Vector Machines"; Optim. Meth. Soft. 20, 2005, 353-378) and the        ***
00015  *** Dai-Fletcher Method DFGPM (Y. Dai and R. Fletcher,"New Algorithms for  ***
00016  *** Singly Linear Constrained Quadratic Programs Subject to Lower and      ***
00017  *** Upper Bounds"; Math. Prog. to appear).                                 ***
00018  ***                                                                        ***
00019  *** Authors:                                                               ***
00020  ***  Thomas Serafini, Luca Zanni                                           ***
00021  ***   Dept. of Mathematics, University of Modena and Reggio Emilia - ITALY ***
00022  ***   serafini.thomas@unimo.it, zanni.luca@unimo.it                        ***
00023  ***  Gaetano Zanghirati                                                    ***
00024  ***   Dept. of Mathematics, University of Ferrara - ITALY                  ***
00025  ***   g.zanghirati@unife.it                                                ***
00026  ***                                                                        ***
00027  *** Software homepage: http://dm.unife.it/gpdt                             ***
00028  ***                                                                        ***
00029  *** This work is supported by the Italian FIRB Projects                    ***
00030  ***      'Statistical Learning: Theory, Algorithms and Applications'       ***
00031  ***      (grant RBAU01877P), http://slipguru.disi.unige.it/ASTA            ***
00032  *** and                                                                    ***
00033  ***      'Parallel Algorithms and Numerical Nonlinear Optimization'        ***
00034  ***      (grant RBAU01JYPN), http://dm.unife.it/pn2o                       ***
00035  ***                                                                        ***
00036  *** Copyright (C) 2004-2008 by T. Serafini, G. Zanghirati, L. Zanni.       ***
00037  ***                                                                        ***
00038  ***                     COPYRIGHT NOTIFICATION                             ***
00039  ***                                                                        ***
00040  *** Permission to copy and modify this software and its documentation      ***
00041  *** for internal research use is granted, provided that this notice is     ***
00042  *** retained thereon and on all copies or modifications. The authors and   ***
00043  *** their respective Universities makes no representations as to the       ***
00044  *** suitability and operability of this software for any purpose. It is    ***
00045  *** provided "as is" without express or implied warranty.                  ***
00046  *** Use of this software for commercial purposes is expressly prohibited   ***
00047  *** without contacting the authors.                                        ***
00048  ***                                                                        ***
00049  *** This program is free software; you can redistribute it and/or modify   ***
00050  *** it under the terms of the GNU General Public License as published by   ***
00051  *** the Free Software Foundation; either version 3 of the License, or      ***
00052  *** (at your option) any later version.                                    ***
00053  ***                                                                        ***
00054  *** This program is distributed in the hope that it will be useful,        ***
00055  *** but WITHOUT ANY WARRANTY; without even the implied warranty of         ***
00056  *** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          ***
00057  *** GNU General Public License for more details.                           ***
00058  ***                                                                        ***
00059  *** You should have received a copy of the GNU General Public License      ***
00060  *** along with this program; if not, write to the Free Software            ***
00061  *** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.              ***
00062  ***                                                                        ***
00063  *** File:     gpm.cpp                                                      ***
00064  *** Type:     scalar                                                       ***
00065  *** Version:  1.0                                                          ***
00066  *** Date:     October, 2005                                                ***
00067  *** Revision: 1                                                            ***
00068  ***                                                                        ***
00069  *** SHOGUN adaptions  Written (W) 2006-2008 Soeren Sonnenburg              ***
00070  ******************************************************************************/
00071 #include <stdio.h>
00072 #include <stdlib.h>
00073 #include <string.h>
00074 #include <math.h>
00075 #include "classifier/svm/gpdt.h"
00076 #include "lib/io.h"
00077 
00078 #define maxvpm           30000  /* max number of method iterations allowed  */
00079 #define maxprojections   200
00080 #define in               8000   /* max size of the QP problem to solve      */
00081 #define alpha_max        1e10
00082 #define alpha_min        1e-10
00083 
00084 extern unsigned int Randnext;
00085 #define ThRand    (Randnext = Randnext * 1103515245L + 12345L)
00086 #define ThRandPos ((Randnext = Randnext * 1103515245L + 12345L) & 0x7fffffff)
00087 
00088 int InnerProjector(int method, int n, int *iy, double e, double *qk,
00089                    double l, double u, double *x, double &lambda);
00090 
00091 /* Uncomment if you want to allocate vectors on the stack  *
00092  * instead of the heap. On some architectures this helps   *
00093  * improving speed, but may generate a stack overflow      */
00094 // #define VARIABLES_ON_STACK
00095 
00096 /* Uncomment if you want to use the adaptive steplength    *
00097    in the GVPM solver                                      */
00098 #define VPM_ADA
00099 
00100 
00101 /******************************************************************************
00102  *** Generalized Variable Projection Method (T. Serafini, G. Zanghirati,    ***
00103  *** L. Zanni, "Gradient Projection Methods for Quadratic Programs and      ***
00104  *** Applications in Training Support Vector Machines";                     ***
00105  *** Optim. Meth. Soft. 20, 2005, 353-378)                                  ***
00106  ******************************************************************************/
00107 int gvpm(int Projector, int n, float *vecA, double *b, double c, double e, 
00108          int *iy, double *x, double tol, int *ls, int *proj)
00109 {
00110   int     i, j, iter, it, it2, luv, info;
00111   double  gd, max, normd, dAd, lam, lamnew, alpha, kktlam, ak, bk;
00112 
00113   int     lscount = 0, projcount = 0;
00114   double  eps     = 1.0e-16;
00115   double  DELTAsv, ProdDELTAsv;
00116   double  lam_ext;
00117 
00118   /* solver-specific settings */
00119 #ifdef VPM_ADA
00120   int     nc = 1, ia1 = -1;
00121   double  alpha1, alpha2;
00122 #endif
00123 
00124   /* allocation-dependant settings */
00125 #ifdef VARIABLES_ON_STACK
00126 
00127   int     ipt[in], ipt2[in], uv[in];
00128   double  g[in], y[in], tempv[in], d[in], Ad[in], t[in];
00129 
00130 #else
00131 
00132   int     *ipt, *ipt2, *uv;
00133   double  *g, *y, *tempv, *d, *Ad, *t;
00134 
00135   /*** array allocations ***/
00136   ipt   = (int    *)malloc(n * sizeof(int   ));
00137   ipt2  = (int    *)malloc(n * sizeof(int   ));
00138   uv    = (int    *)malloc(n * sizeof(int   ));
00139   g     = (double *)malloc(n * sizeof(double));
00140   y     = (double *)malloc(n * sizeof(double));
00141   d     = (double *)malloc(n * sizeof(double));
00142   Ad    = (double *)malloc(n * sizeof(double));
00143   t     = (double *)malloc(n * sizeof(double));
00144   tempv = (double *)malloc(n * sizeof(double));
00145 
00146 #endif
00147 
00148   DELTAsv = EPS_SV * c;
00149   if (tol <= 1.0e-5 || n <= 20)
00150       ProdDELTAsv = 0.0F;
00151   else
00152       ProdDELTAsv = EPS_SV * c;
00153 
00154   for (i = 0; i < n; i++)
00155       tempv[i] = -x[i];
00156   lam_ext = 0.0;
00157   projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, x, lam_ext);
00158 
00159   /* compute g = A*x + b in sparse form          *
00160    * (inline computation for better perfomrance) */
00161   {
00162     float *tempA;
00163 
00164     it = 0;
00165     for (i = 0; i < n; i++)
00166         if (fabs(x[i]) > ProdDELTAsv*1e-2)
00167             ipt[it++] = i;
00168 
00169     memset(t, 0, n*sizeof(double));
00170     for (i = 0; i < it; i++)
00171     {
00172         tempA = vecA + ipt[i]*n;
00173         for (j = 0; j < n; j++)
00174             t[j] += (tempA[j] * x[ipt[i]]);
00175     }
00176   }
00177 
00178   for (i = 0; i < n; i++)
00179   {
00180     g[i] = t[i] + b[i],
00181     y[i] = g[i] - x[i];
00182   }
00183 
00184   projcount += InnerProjector(Projector, n, iy, e, y, 0, c, tempv, lam_ext);
00185 
00186   max = alpha_min;
00187   for (i = 0; i < n; i++)
00188   {
00189       y[i] = tempv[i] - x[i];
00190       if (fabs(y[i]) > max)
00191           max = fabs(y[i]);
00192   }
00193 
00194   if (max < c*tol*1e-3)
00195   {
00196       lscount = 0;
00197       iter    = 0;
00198       goto Clean;
00199   }
00200 
00201   alpha = 1.0 / max;
00202 
00203   for (iter = 1; iter <= maxvpm; iter++)
00204   {
00205       for (i = 0; i < n; i++)
00206           tempv[i] = alpha*g[i] - x[i];
00207 
00208       projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, y, lam_ext);
00209 
00210       gd = 0.0;
00211       for (i = 0; i < n; i++)
00212       {
00213           d[i] = y[i] - x[i];
00214           gd  += d[i] * g[i];
00215       }
00216 
00217       /* compute Ad = A*d  or  Ad = Ay-t depending on their sparsity  *
00218        * (inline computation for better perfomrance)                  */
00219       {
00220          float *tempA;
00221 
00222          it = it2 = 0;
00223          for (i = 0; i < n; i++)
00224              if (fabs(d[i]) > (ProdDELTAsv*1.0e-2))
00225                  ipt[it++] = i;
00226          for (i = 0; i < n; i++)
00227              if (fabs(y[i]) > ProdDELTAsv)
00228                  ipt2[it2++] = i;
00229 
00230          memset(Ad, 0, n*sizeof(double));
00231          if (it < it2) // Ad = A*d
00232          {
00233              for (i = 0; i < it; i++)
00234              {
00235                  tempA = vecA + ipt[i]*n;
00236                  for (j = 0; j < n; j++)
00237                      Ad[j] += (tempA[j] * d[ipt[i]]);
00238              }
00239          }
00240          else          // Ad = A*y - t
00241          {
00242             for (i = 0; i < it2; i++)
00243             {
00244                 tempA = vecA + ipt2[i]*n;
00245                 for (j = 0; j < n; j++)
00246                     Ad[j] += (tempA[j] * y[ipt2[i]]);
00247             }
00248             for (j = 0; j < n; j++)
00249                 Ad[j] -= t[j];
00250          }
00251       }
00252 
00253       normd = 0.0;
00254       for (i = 0; i < n; i++)
00255           normd += d[i] * d[i];
00256 
00257       dAd = 0.0;
00258       for (i = 0; i < n; i++)
00259           dAd += d[i]*Ad[i];
00260 
00261       if (dAd > eps*normd && gd < 0.0)
00262       {
00263           lam = lamnew = -gd/dAd;
00264           if (lam > 1.0 || lam < 0.0)
00265               lam = 1.0;
00266           else
00267               lscount++;
00268 
00269 #ifdef VPM_ADA
00270 
00271           /*** use the adaptive switching rule for steplength selection ***/
00272 
00273           // compute alpha1 = (d'* (d.*diaga)) / (d'*Ad);
00274           alpha1 = normd / dAd;
00275 
00276           // alpha2 = d'*Ad / (Ad' * (Ad./diaga));
00277           alpha2 = 0.0;
00278           for (i = 0; i < n; i++)
00279                alpha2 += Ad[i] * Ad[i];
00280           alpha2 = dAd / alpha2;
00281 
00282           if ( (nc > 2
00283                 && (
00284                      (ia1 == 1
00285                       && (
00286                            lamnew < 0.1 || (alpha1 > alpha && alpha2 < alpha)
00287                          )
00288                      )
00289                      ||
00290                      (ia1 == -1
00291                       && (
00292                            lamnew > 5.0 || (alpha1 > alpha && alpha2 < alpha)
00293                          )
00294                      )
00295                    )
00296                )
00297                || nc > 9 )
00298           {
00299               ia1 = -ia1;
00300               nc  = 0;
00301           }
00302 
00303           if (ia1 == 1)
00304               alpha = alpha1;
00305           else
00306               alpha = alpha2;
00307 
00308           if (alpha < alpha_min)
00309               alpha = alpha_min;
00310           else if (alpha > alpha_max)
00311               alpha = alpha_max;
00312 
00313           nc++;
00314 
00315 #else
00316 
00317           /*** use the fixed switching rule for steplength selection ***/
00318 
00319           if ((iter % 6) < 3) // alpha = d'*Ad / (Ad' * (Ad./diaga));
00320           {
00321               alpha = 0.0;
00322               for (i = 0; i < n; i++)
00323                   alpha += Ad[i] * Ad[i];
00324               alpha = dAd / alpha;
00325           }
00326           else                // alpha = (d'* (d.*diaga)) / (d'*Ad);
00327           {
00328               alpha = 0.0;
00329               for (i = 0; i < n; i++)
00330                   alpha += d[i] * d[i];
00331               alpha = alpha / dAd;
00332           }
00333 
00334 #endif
00335 
00336       }
00337       else
00338       {
00339           lam   = 1.0;
00340           alpha = alpha_max;
00341       }
00342 
00343       for (i = 0; i < n; i++)
00344       {
00345           x[i] = x[i] + lam*d[i];
00346           t[i] = t[i] + lam*Ad[i];
00347           g[i] = b[i] + t[i];
00348       }
00349 
00350       /*** stopping criterion based on KKT conditions ***/
00351       bk = 0.0;
00352       ak = 0.0;
00353       for (i = 0; i < n; i++)
00354       {
00355           bk +=  x[i] * x[i];
00356           ak +=  d[i] * d[i];
00357       }
00358 
00359       if (lam*sqrt(ak) < tol*10 * sqrt(bk))
00360       {
00361           it     = 0;
00362           luv    = 0;
00363           kktlam = 0.0;
00364           for (i = 0; i < n; i++)
00365           {
00366               if (x[i] > DELTAsv && x[i] < c-DELTAsv)
00367               {
00368                   ipt[it++] = i;
00369                   kktlam    = kktlam - iy[i]*g[i];
00370               }
00371               else
00372                   uv[luv++] = i;
00373           }
00374 
00375           if (it == 0)
00376           {
00377               if (lam*sqrt(ak) < tol*0.5 * sqrt(bk))
00378               goto Clean;
00379           }
00380           else
00381           {
00382               kktlam = kktlam/it;
00383               info   = 1;
00384               for (i = 0; i < it; i++)
00385                   if (fabs(iy[ipt[i]]*g[ipt[i]]+kktlam) > tol)
00386                   {
00387                       info = 0;
00388                       break;
00389                   }
00390 
00391               if (info == 1)
00392                   for (i = 0; i < luv; i++)
00393                       if (x[uv[i]] <= DELTAsv)
00394                       {
00395                           if (g[uv[i]] + kktlam*iy[uv[i]] < -tol)
00396                           {
00397                               info = 0;
00398                               break;
00399                           }
00400                       }
00401                       else
00402                       {
00403                           if (g[uv[i]] + kktlam*iy[uv[i]] > tol)
00404                           {
00405                               info = 0;
00406                               break;
00407                           }
00408                       }
00409 
00410               if (info == 1)
00411                   goto Clean;
00412           }
00413       } // stopping rule based on the norm of d_k
00414   }
00415 
00416   SG_SWARNING("GVPM exits after maxvpm = %d iterations.\n", maxvpm);
00417 
00418 Clean:
00419 
00420   /*** allocation-dependant freeing ***/
00421 #ifndef VARIABLES_ON_STACK
00422   free(t);
00423   free(uv);
00424   free(ipt2);
00425   free(ipt);
00426   free(g);
00427   free(y);
00428   free(tempv);
00429   free(d);
00430   free(Ad);
00431 #endif
00432 
00433   if (ls != NULL)   *ls   = lscount;
00434   if (proj != NULL) *proj = projcount;
00435   return(iter);
00436 }
00437 
00438 /******************************************************************************
00439  *** Dai-Fletcher QP solver (Y. Dai and R. Fletcher,"New Algorithms for     ***
00440  *** Singly Linear Constrained Quadratic Programs Subject to Lower and      *** 
00441  *** Upper Bounds"; Math. Prog. to appear)                                  ***
00442  ******************************************************************************/
00443 int FletcherAlg2A(int Projector, int n, float *vecA, double *b, double c, 
00444                   double e, int *iy, double *x, double tol, int *ls, int *proj)
00445 {
00446   int    i, j, iter, it, it2, luv, info, lscount = 0, projcount = 0;
00447   double gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam, lam_ext;
00448   double eps     = 1.0e-16;
00449   double DELTAsv, ProdDELTAsv;
00450 
00451   /*** variables for the adaptive nonmonotone linesearch ***/
00452   int    L, llast;
00453   double fr, fbest, fv, fc, fv0;
00454 
00455   /*** allocation-dependant settings ***/
00456 #ifdef VARIABLES_ON_STACK
00457 
00458   int    ipt[in], ipt2[in], uv[in];
00459   double g[in], y[in], tempv[in], d[in], Ad[in], t[in],
00460          xplus[in], tplus[in], sk[in], yk[in];
00461 #else
00462 
00463   int    *ipt, *ipt2, *uv;
00464   double *g, *y, *tempv, *d, *Ad, *t, *xplus, *tplus, *sk, *yk;
00465 
00466   /*** arrays allocation ***/
00467   ipt   = (int    *)malloc(n * sizeof(int   ));
00468   ipt2  = (int    *)malloc(n * sizeof(int   ));
00469   uv    = (int    *)malloc(n * sizeof(int   ));
00470   g     = (double *)malloc(n * sizeof(double));
00471   y     = (double *)malloc(n * sizeof(double));
00472   tempv = (double *)malloc(n * sizeof(double));
00473   d     = (double *)malloc(n * sizeof(double));
00474   Ad    = (double *)malloc(n * sizeof(double));
00475   t     = (double *)malloc(n * sizeof(double));
00476   xplus = (double *)malloc(n * sizeof(double));
00477   tplus = (double *)malloc(n * sizeof(double));
00478   sk    = (double *)malloc(n * sizeof(double));
00479   yk    = (double *)malloc(n * sizeof(double));
00480 
00481 #endif
00482 
00483   DELTAsv = EPS_SV * c;
00484   if (tol <= 1.0e-5 || n <= 20)
00485       ProdDELTAsv = 0.0F;
00486   else
00487       ProdDELTAsv = EPS_SV * c;
00488 
00489   for (i = 0; i < n; i++)
00490       tempv[i] = -x[i];
00491 
00492   lam_ext = 0.0;
00493   
00494   projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, x, lam_ext);
00495 
00496   // g = A*x + b;
00497   // SparseProd(n, t, A, x, ipt);
00498   {
00499     float *tempA;
00500 
00501     it = 0;
00502     for (i = 0; i < n; i++)
00503         if (fabs(x[i]) > ProdDELTAsv)
00504             ipt[it++] = i;
00505 
00506     memset(t, 0, n*sizeof(double));
00507     for (i = 0; i < it; i++)
00508     {
00509          tempA = vecA + ipt[i] * n;
00510          for (j = 0; j < n; j++)
00511              t[j] += (tempA[j]*x[ipt[i]]);
00512     }
00513   }
00514 
00515   for (i = 0; i < n; i++)
00516   {
00517     g[i] = t[i] + b[i],
00518     y[i] = g[i] - x[i];
00519   }
00520 
00521   projcount += InnerProjector(Projector, n, iy, e, y, 0, c, tempv, lam_ext);
00522 
00523   max = alpha_min;
00524   for (i = 0; i < n; i++)
00525   {
00526       y[i] = tempv[i] - x[i];
00527       if (fabs(y[i]) > max)
00528           max = fabs(y[i]);
00529   }
00530 
00531   if (max < c*tol*1e-3)
00532   {
00533       lscount = 0;
00534       iter    = 0;
00535       goto Clean;
00536   }
00537 
00538   alpha = 1.0 / max;
00539 
00540   fv0   = 0.0;
00541   for (i = 0; i < n; i++)
00542       fv0 += x[i] * (0.5*t[i] + b[i]);
00543 
00544   /*** adaptive nonmonotone linesearch ***/
00545   L     = 2;
00546   fr    = alpha_max;
00547   fbest = fv0;
00548   fc    = fv0;
00549   llast = 0;
00550   akold = bkold = 0.0;
00551 
00552   for (iter = 1; iter <= maxvpm; iter++)
00553   {
00554       for (i = 0; i < n; i++)
00555           tempv[i] = alpha*g[i] - x[i];
00556 
00557       projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, y, lam_ext);
00558 
00559       gd = 0.0;
00560       for (i = 0; i < n; i++)
00561       {
00562           d[i] = y[i] - x[i];
00563           gd  += d[i] * g[i];
00564       }
00565 
00566       /* compute Ad = A*d  or  Ad = A*y - t depending on their sparsity */
00567       {
00568          float *tempA;
00569 
00570          it = it2 = 0;
00571          for (i = 0; i < n; i++)
00572              if (fabs(d[i]) > (ProdDELTAsv*1.0e-2))
00573                  ipt[it++]   = i;
00574          for (i = 0; i < n; i++)
00575              if (fabs(y[i]) > ProdDELTAsv)
00576                  ipt2[it2++] = i;
00577 
00578          memset(Ad, 0, n*sizeof(double));
00579          if (it < it2) // compute Ad = A*d
00580          {
00581             for (i = 0; i < it; i++)
00582             {
00583                 tempA = vecA + ipt[i]*n;
00584                 for (j = 0; j < n; j++)
00585                     Ad[j] += (tempA[j] * d[ipt[i]]);
00586             }
00587          }
00588          else          // compute Ad = A*y-t
00589          {
00590             for (i = 0; i < it2; i++)
00591             {
00592                 tempA = vecA + ipt2[i]*n;
00593                 for (j = 0; j < n; j++)
00594                     Ad[j] += (tempA[j] * y[ipt2[i]]);
00595             }
00596             for (j = 0; j < n; j++)
00597                 Ad[j] -= t[j];
00598          }
00599       }
00600 
00601       ak = 0.0;
00602       for (i = 0; i < n; i++)
00603           ak += d[i] * d[i];
00604 
00605       bk = 0.0;
00606       for (i = 0; i < n; i++)
00607           bk += d[i]*Ad[i];
00608 
00609       if (bk > eps*ak && gd < 0.0)    // ak is normd
00610           lamnew = -gd/bk;
00611       else
00612           lamnew = 1.0;
00613 
00614       fv = 0.0;
00615       for (i = 0; i < n; i++)
00616       {
00617           xplus[i] = x[i] + d[i];
00618           tplus[i] = t[i] + Ad[i];
00619           fv      += xplus[i] * (0.5*tplus[i] + b[i]);
00620       }
00621 
00622       if ((iter == 1 && fv >= fv0) || (iter > 1 && fv >= fr))
00623       {
00624           lscount++;
00625           fv = 0.0;
00626           for (i = 0; i < n; i++)
00627           {
00628               xplus[i] = x[i] + lamnew*d[i];
00629               tplus[i] = t[i] + lamnew*Ad[i];
00630               fv      += xplus[i] * (0.5*tplus[i] + b[i]);
00631           }
00632       }
00633 
00634       for (i = 0; i < n; i++)
00635       {
00636           sk[i] = xplus[i] - x[i];
00637           yk[i] = tplus[i] - t[i];
00638           x[i]  = xplus[i];
00639           t[i]  = tplus[i];
00640           g[i]  = t[i] + b[i];
00641       }
00642 
00643       // update the line search control parameters
00644 
00645       if (fv < fbest)
00646       {
00647           fbest = fv;
00648           fc    = fv;
00649           llast = 0;
00650       }
00651       else
00652       {
00653           fc = (fc > fv ? fc : fv);
00654           llast++;
00655           if (llast == L)
00656           {
00657               fr    = fc;
00658               fc    = fv;
00659               llast = 0;
00660           }
00661       }
00662 
00663       ak = bk = 0.0;
00664       for (i = 0; i < n; i++)
00665       {
00666           ak += sk[i] * sk[i];
00667           bk += sk[i] * yk[i];
00668       }
00669 
00670       if (bk < eps*ak)
00671           alpha = alpha_max;
00672       else
00673       {
00674           if (bkold < eps*akold)
00675               alpha = ak/bk;
00676           else
00677               alpha = (akold+ak)/(bkold+bk);
00678 
00679           if (alpha > alpha_max)
00680               alpha = alpha_max;
00681           else if (alpha < alpha_min)
00682               alpha = alpha_min;
00683       }
00684 
00685       akold = ak;
00686       bkold = bk;
00687 
00688       /*** stopping criterion based on KKT conditions ***/
00689 
00690       bk = 0.0;
00691       for (i = 0; i < n; i++)
00692           bk +=  x[i] * x[i];
00693 
00694       if (sqrt(ak) < tol*10 * sqrt(bk))
00695       {
00696           it     = 0;
00697           luv    = 0;
00698           kktlam = 0.0;
00699           for (i = 0; i < n; i++)
00700           {
00701               if ((x[i] > DELTAsv) && (x[i] < c-DELTAsv))
00702               {
00703                   ipt[it++] = i;
00704                   kktlam    = kktlam - iy[i]*g[i];
00705               }
00706               else
00707                   uv[luv++] = i;
00708           }
00709 
00710           if (it == 0)
00711           {
00712               if (sqrt(ak) < tol*0.5 * sqrt(bk))
00713                   goto Clean;
00714           }
00715           else
00716           {
00717 
00718               kktlam = kktlam/it;
00719               info   = 1;
00720               for (i = 0; i < it; i++)
00721                   if ( fabs(iy[ipt[i]] * g[ipt[i]] + kktlam) > tol )
00722                   {
00723                       info = 0;
00724                       break;
00725                   }
00726 
00727               if (info == 1)
00728                   for (i = 0; i < luv; i++)
00729                       if (x[uv[i]] <= DELTAsv)
00730                       {
00731                           if (g[uv[i]] + kktlam*iy[uv[i]] < -tol)
00732                           {
00733                               info = 0;
00734                               break;
00735                           }
00736                       }
00737                       else
00738                       {
00739                           if (g[uv[i]] + kktlam*iy[uv[i]] > tol)
00740                           {
00741                               info = 0;
00742                               break;
00743                           }
00744                       }
00745 
00746               if (info == 1)
00747                   goto Clean;
00748           }
00749       }
00750   }
00751 
00752   SG_SWARNING("Dai-Fletcher method exits after maxvpm = %d iterations.\n", maxvpm);
00753 
00754 Clean:
00755 
00756 #ifndef VARIABLES_ON_STACK
00757   free(sk);
00758   free(yk);
00759   free(tplus);
00760   free(xplus);
00761   free(t);
00762   free(uv);
00763   free(ipt2);
00764   free(ipt);
00765   free(g);
00766   free(y);
00767   free(tempv);
00768   free(d);
00769   free(Ad);
00770 #endif
00771 
00772   if (ls != NULL)   *ls   = lscount;
00773   if (proj != NULL) *proj = projcount;
00774   return(iter);
00775 
00776 }
00777 
00778 /******************************************************************************/
00779 /*** Encapsulating method to call the chosen Gradient Projection Method     ***/
00780 /******************************************************************************/
00781 int gpm_solver(int Solver, int Projector, int n, float *A, double *b, double c, 
00782                double e, int *iy, double *x, double tol, int *ls, int *proj)
00783 {
00784   /*** Uncomment the following if you need to scale the QP Hessian matrix
00785    *** before calling the chosen solver
00786   int    i, j;
00787   float  *ptrA;
00788   double max, s;
00789 
00790   max = fabs(A[0][0]);
00791   for (i = 1; i < n; i++)
00792       if (fabs(A[i][i]) > max)
00793           max = fabs(A[i][i]);
00794 
00795   s    = 1.0 / max;
00796   ptrA = vecA;
00797   for (i = 0; i < n; i++)
00798       for (j = 0;j < n; j++)
00799           *ptrA++ = (float)(A[i][j]*s);
00800 
00801   if (Solver == SOLVER_FLETCHER)
00802       j = FletcherAlg2A(n, vecA, b, c/s, e/s, iy, x, tol, ls);
00803   else
00804       j = gvpm(n, vecA, b, c/s, e/s, iy, x, tol, ls);
00805 
00806   for (i = 0; i < n; i++)
00807       x[i] *= s;
00808   ***/
00809 
00810   /*** calling the chosen solver with unscaled data ***/
00811   if (Solver == SOLVER_FLETCHER)
00812     return FletcherAlg2A(Projector, n, A, b, c, e, iy, x, tol, ls, proj);
00813   else
00814     return gvpm(Projector, n, A, b, c, e, iy, x, tol, ls, proj);
00815 }
00816 
00817 /******************************************************************************
00818  *** Piecewise linear monotone target function for the Dai-Fletcher         *** 
00819  *** projector (Y. Dai and R. Fletcher, "New Algorithms for Singly Linear   ***
00820  *** Constrained Quadratic Programs Subject to Lower and Upper Bounds";     ***
00821  *** Math. Prog. to appear)                                                 ***
00822  ******************************************************************************/
00823 double ProjectR(double *x, int n, double lambda, int *a, double b, double *c, 
00824                 double l, double u)
00825 {
00826   int    i;
00827   double r = 0.0;
00828 
00829   for (i = 0; i < n; i++)
00830   {
00831       x[i] = -c[i] + lambda*(double)a[i];
00832       if (x[i] >= u) x[i] = u;
00833       else if (x[i] < l) x[i] = l;
00834       r += (double)a[i]*x[i];
00835   }
00836 
00837   return (r - b);
00838 }
00839 
00840 /******************************************************************************
00841  *** Dai-Fletcher QP projector (Y. Dai and R. Fletcher, "New Algorithms for ***
00842  *** Singly Linear Constrained Quadratic Programs Subject to Lower and      ***
00843  *** Upper Bounds"; Math. Prog. to appear)                                  ***
00844  ******************************************************************************/
00845 /***                                                                        ***
00846  *** Solves the problem        min  x'*x/2 + c'*x                           ***
00847  ***                       subj to  a'*x - b = 0                            ***
00848  ***                                l <= x <= u                             ***
00849  ******************************************************************************/
00850 int ProjectDai(int n, int *a, double b, double *c, double l, double u, 
00851                double *x, double &lam_ext)
00852 {
00853   double lambda, lambdal, lambdau, dlambda, lambda_new, tol_lam;
00854   double r, rl, ru, s, tol_r;
00855   int    iter;
00856 
00857   tol_lam = 1.0e-11;
00858   tol_r   = 1.0e-10 * sqrt((u-l)*(double)n);
00859   lambda  = lam_ext;
00860   dlambda = 0.5;
00861   iter    = 1;
00862   b       = -b;
00863 
00864   // Bracketing Phase
00865   r = ProjectR(x, n, lambda, a, b, c, l, u);
00866   if (fabs(r) < tol_r)
00867       return 0;
00868 
00869   if (r < 0.0)
00870   {
00871       lambdal = lambda;
00872       rl      = r;
00873       lambda  = lambda + dlambda;
00874       r       = ProjectR(x, n, lambda, a, b, c, l, u);
00875       while (r < 0.0)
00876       {
00877          lambdal = lambda;
00878          s       = rl/r - 1.0;
00879          if (s < 0.1) s = 0.1;
00880          dlambda = dlambda + dlambda/s;
00881          lambda  = lambda + dlambda;
00882          rl      = r;
00883          r       = ProjectR(x, n, lambda, a, b, c, l, u);
00884       }
00885       lambdau = lambda;
00886       ru      = r;
00887   }
00888   else
00889   {
00890       lambdau = lambda;
00891       ru      = r;
00892       lambda  = lambda - dlambda;
00893       r       = ProjectR(x, n, lambda, a, b, c, l, u);
00894       while (r > 0.0)
00895       {
00896          lambdau = lambda;
00897          s       = ru/r - 1.0;
00898          if (s < 0.1) s = 0.1;
00899          dlambda = dlambda + dlambda/s;
00900          lambda  = lambda - dlambda;
00901          ru      = r;
00902          r       = ProjectR(x, n, lambda, a, b, c, l, u);
00903       }
00904     lambdal = lambda;
00905     rl      = r;
00906   }
00907 
00908 
00909   // Secant Phase
00910   s       = 1.0 - rl/ru;
00911   dlambda = dlambda/s;
00912   lambda  = lambdau - dlambda;
00913   r       = ProjectR(x, n, lambda, a, b, c, l, u);
00914 
00915   while (   fabs(r) > tol_r 
00916          && dlambda > tol_lam * (1.0 + fabs(lambda)) 
00917          && iter    < maxprojections                )
00918   {
00919      iter++;
00920      if (r > 0.0)
00921      {
00922          if (s <= 2.0)
00923          {
00924              lambdau = lambda;
00925              ru      = r;
00926              s       = 1.0 - rl/ru;
00927              dlambda = (lambdau - lambdal) / s;
00928              lambda  = lambdau - dlambda;
00929          }
00930          else
00931          {
00932              s          = ru/r-1.0;
00933              if (s < 0.1) s = 0.1;
00934              dlambda    = (lambdau - lambda) / s;
00935              lambda_new = 0.75*lambdal + 0.25*lambda;
00936              if (lambda_new < (lambda - dlambda))
00937                  lambda_new = lambda - dlambda;
00938              lambdau    = lambda;
00939              ru         = r;
00940              lambda     = lambda_new;
00941              s          = (lambdau - lambdal) / (lambdau - lambda);
00942          }
00943      }
00944      else
00945      {
00946          if (s >= 2.0)
00947          {
00948              lambdal = lambda;
00949              rl      = r;
00950              s       = 1.0 - rl/ru;
00951              dlambda = (lambdau - lambdal) / s;
00952              lambda  = lambdau - dlambda;
00953          }
00954          else
00955          {
00956              s          = rl/r - 1.0;
00957              if (s < 0.1) s = 0.1;
00958              dlambda    = (lambda-lambdal) / s;
00959              lambda_new = 0.75*lambdau + 0.25*lambda;
00960              if (lambda_new > (lambda + dlambda))
00961                  lambda_new = lambda + dlambda;
00962              lambdal    = lambda;
00963              rl         = r;
00964              lambda     = lambda_new;
00965              s          = (lambdau - lambdal) / (lambdau-lambda);
00966          }
00967      }
00968      r = ProjectR(x, n, lambda, a, b, c, l, u);
00969   }
00970 
00971   lam_ext = lambda;
00972   if (iter >= maxprojections)
00973       SG_SERROR("Projector exits after max iterations: %d\n", iter);
00974 
00975   return (iter);
00976 }
00977 
00978 #define SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
00979 
00980 /*** Median computation using Quick Select ***/
00981 double quick_select(double *arr, int n)
00982 {
00983   int low, high ;
00984   int median;
00985   int middle, l, h;
00986 
00987   low    = 0; 
00988   high   = n-1;
00989   median = (low + high) / 2;
00990   
00991   for (;;)
00992   {
00993     if (high <= low)
00994         return arr[median];
00995 
00996     if (high == low + 1)
00997     {
00998         if (arr[low] > arr[high])
00999             SWAP(arr[low], arr[high]);
01000         return arr[median];
01001     }
01002 
01003     middle = (low + high) / 2;
01004     if (arr[middle] > arr[high]) SWAP(arr[middle], arr[high]);
01005     if (arr[low]    > arr[high]) SWAP(arr[low],    arr[high]);
01006     if (arr[middle] > arr[low])  SWAP(arr[middle], arr[low]);
01007 
01008     SWAP(arr[middle], arr[low+1]);
01009 
01010     l = low + 1;
01011     h = high;
01012     for (;;)
01013     {
01014       do l++; while (arr[low] > arr[l]);
01015       do h--; while (arr[h]   > arr[low]);
01016       if (h < l)
01017           break;
01018       SWAP(arr[l], arr[h]);
01019     }
01020 
01021     SWAP(arr[low], arr[h]);
01022     if (h <= median)
01023         low = l;
01024     if (h >= median)
01025         high = h - 1;
01026   }
01027 }
01028 
01029 /******************************************************************************
01030  *** Pardalos-Kovoor projector (P.M. Pardalos and N. Kovoor, "An Algorithm  ***
01031  *** for a Singly Constrained Class of Quadratic Programs Subject to Upper  ***
01032  *** and Lower Bounds"; Math. Prog. 46, 1990, 321-328).                     ***
01033  ******************************************************************************
01034  *** Solves the problem                                                     ***
01035  ***                       min  x'*x/2 + qk'*x                              ***
01036  ***                   subj to  iy'*x + e = 0                               ***
01037  ***                            l <= x <= u                                 ***
01038  ***                            iy(i) ~= 0                                  ***
01039  ******************************************************************************/
01040 
01041 int Pardalos(int n, int *iy, double e,   double *qk,
01042                              double low, double up,  double *x)
01043 {
01044   int    i, l, iter; /* conters    */
01045   int    luv, lxint; /* dimensions */
01046   double d, xmin, xmax, xmold, xmid, xx, ts, sw, s, s1, testsum;
01047 
01048   /*** allocation-dependant settings ***/
01049 #ifdef VARIABLES_ON_STACK
01050   int    uv[in], uvt[in];
01051   double xint[2*in+2], xint2[2*in+2], a[in], b[in], at[in], bt[in];
01052   double newdia[in], newdt[in];
01053 #else
01054 
01055   int    *uv, *uvt;
01056   double *xint, *xint2, *a, *b, *at, *bt, *newdia, *newdt;
01057 
01058   /*** arrays allocation ***/
01059   uv     = (int    *)malloc(n * sizeof(int           ));
01060   uvt    = (int    *)malloc(n * sizeof(int           ));
01061   a      = (double *)malloc(n * sizeof(double        ));
01062   b      = (double *)malloc(n * sizeof(double        ));
01063   at     = (double *)malloc(n * sizeof(double        ));
01064   bt     = (double *)malloc(n * sizeof(double        ));
01065   newdia = (double *)malloc(n * sizeof(double        ));
01066   newdt  = (double *)malloc(n * sizeof(double        ));
01067   xint   = (double *)malloc((2*n + 2) * sizeof(double));
01068   xint2  = (double *)malloc((2*n + 2) * sizeof(double));
01069 
01070 #endif
01071 
01072   d = 0.0;
01073   for (i = 0; i < n; i++)
01074       d += iy[i] * qk[i];
01075   d = 0.5 * (d-e);
01076 
01077   for (i = 0; i < n; i++)
01078   {
01079       /* The following computations should divide by iy[i] instead           *
01080        * of multiply by iy[i], but this is correct for binary classification *
01081        * with labels -1 and 1.                                               */
01082       if (iy[i] > 0)
01083       {
01084           a[i] = ((qk[i] + low) * iy[i]) * 0.5;
01085           b[i] = ((up + qk[i]) * iy[i]) * 0.5;
01086       }
01087       else
01088       {
01089           b[i] = ((qk[i] + low) * iy[i]) * 0.5;
01090           a[i] = ((up + qk[i]) * iy[i]) * 0.5;
01091       }
01092       newdia[i] = (iy[i]*iy[i]);
01093   }
01094 
01095   xmin = -1e33;
01096   xmax = 1e33;
01097 
01098   /* arrays initialization */
01099   for (i = 0; i < n; i++)
01100   {
01101       uv[i]     = i;    /* contains the unset variables */
01102       xint[i]   = a[i];
01103       xint[n+i] = b[i];
01104   }
01105 
01106   xmid        = xmin;
01107   xint[2*n]   = xmin;
01108   xint[2*n+1] = xmax;
01109   ts          = 0.0;
01110   sw          = 0.0;
01111   luv         = n;
01112   lxint       = 2*n+2;
01113 
01114   iter = 0;
01115   do {
01116      for (i = 0; i < luv; i++)
01117      {
01118          at[i]    = a[uv[i]];
01119          bt[i]    = b[uv[i]];
01120          newdt[i] = newdia[uv[i]];
01121      }
01122 
01123      xmold = xmid;
01124      xmid = quick_select(xint, lxint);
01125      if (xmold == xmid)
01126          xmid = xint[(int)(ThRandPos % lxint)];
01127 
01128      s  = ts;
01129      s1 = sw;
01130      for (i = 0; i < luv; i++)
01131      {
01132          if (xmid > bt[i])
01133              s  += newdt[i]*bt[i];
01134          else if (xmid < at[i])
01135              s  += newdt[i]*at[i];
01136          else
01137              s1 += newdt[i];
01138      }
01139 
01140      testsum = s + s1*xmid;
01141      if (testsum <= (d+(1e-15)))
01142          xmin = xmid;
01143      if (testsum >= (d-(1e-15)))
01144          xmax = xmid;
01145 
01146      l = 0;
01147      for (i = 0; i < lxint; i++)
01148          if((xint[i] >= xmin) && (xint[i] <= xmax))
01149             xint2[l++] = xint[i];
01150      lxint = l;
01151      memcpy(xint, xint2, lxint*sizeof(double));
01152 
01153      l = 0;
01154      for (i = 0; i < luv; i++)
01155      {
01156          if (xmin >= bt[i])
01157              ts += newdt[i]*bt[i];
01158          else if (xmax <= at[i])
01159              ts += newdt[i]*at[i];
01160          else if ((at[i] <= xmin) && (bt[i] >= xmax))
01161              sw += newdt[i];
01162          else
01163              uvt[l++] = uv[i];
01164     }
01165     luv = l;
01166     memcpy(uv, uvt, luv*sizeof(int));
01167     iter++;
01168   } while(luv != 0 && iter < maxprojections);
01169 
01170   if (sw == 0)
01171       xx = xmin;
01172   else
01173       xx = (d-ts) / sw;
01174 
01175   for (i = 0; i < n; i++)
01176   {
01177       if (b[i] <= xmin)
01178           x[i] = b[i];
01179       else if (a[i] >= xmax)
01180           x[i] = a[i];
01181       else if ((a[i]<=xmin) && (xmax<=b[i]))
01182           x[i] = xx;
01183       else
01184           SG_SWARNING("Inner solver troubles...\n");
01185   }
01186 
01187   for (i = 0; i < n; i++)
01188       x[i] = (2.0*x[i]*iy[i]-qk[i]);
01189 
01190 #ifndef VARIABLES_ON_STACK
01191   free(newdt);
01192   free(newdia);
01193   free(a);
01194   free(b);
01195   free(uvt);
01196   free(uv);
01197   free(bt);
01198   free(at);
01199   free(xint2);
01200   free(xint);
01201 #endif
01202 
01203   return(iter);
01204 }
01205 
01206 /******************************************************************************/
01207 /*** Wrapper method to call the selected inner projector                    ***/
01208 /******************************************************************************/
01209 int InnerProjector(int method, int n, int *iy, double e, double *qk,
01210                    double l, double u, double *x, double &lambda)
01211 {
01212   if (method == 0)
01213       return Pardalos(n, iy, e, qk, l, u, x);
01214   else
01215       return ProjectDai(n, iy, e, qk, l, u, x, lambda);
01216 }
01217 
01218 /******************************************************************************/
01219 /*** End of gpm.cpp file                                                    ***/
01220 /******************************************************************************/

SHOGUN Machine Learning Toolbox - Documentation