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

SHOGUN Machine Learning Toolbox - Documentation