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

SHOGUN Machine Learning Toolbox - Documentation