SHOGUN
v2.0.0
|
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 /******************************************************************************/