SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gpm.cpp
Go to the documentation of this file.
1 /******************************************************************************
2  *** GPDT - Gradient Projection Decomposition Technique ***
3  ******************************************************************************
4  *** ***
5  *** GPDT is a C++ software designed to train large-scale Support Vector ***
6  *** Machines for binary classification in both scalar and distributed ***
7  *** memory parallel environments. It uses the Joachims' problem ***
8  *** decomposition technique to split the whole quadratic programming (QP) ***
9  *** problem into a sequence of smaller QP subproblems, each one being ***
10  *** solved by a suitable gradient projection method (GPM). The presently ***
11  *** implemented GPMs are the Generalized Variable Projection Method ***
12  *** GVPM (T. Serafini, G. Zanghirati, L. Zanni, "Gradient Projection ***
13  *** Methods for Quadratic Programs and Applications in Training Support ***
14  *** Vector Machines"; Optim. Meth. Soft. 20, 2005, 353-378) and the ***
15  *** Dai-Fletcher Method DFGPM (Y. Dai and R. Fletcher,"New Algorithms for ***
16  *** Singly Linear Constrained Quadratic Programs Subject to Lower and ***
17  *** Upper Bounds"; Math. Prog. to appear). ***
18  *** ***
19  *** Authors: ***
20  *** Thomas Serafini, Luca Zanni ***
21  *** Dept. of Mathematics, University of Modena and Reggio Emilia - ITALY ***
22  *** serafini.thomas@unimo.it, zanni.luca@unimo.it ***
23  *** Gaetano Zanghirati ***
24  *** Dept. of Mathematics, University of Ferrara - ITALY ***
25  *** g.zanghirati@unife.it ***
26  *** ***
27  *** Software homepage: http://dm.unife.it/gpdt ***
28  *** ***
29  *** This work is supported by the Italian FIRB Projects ***
30  *** 'Statistical Learning: Theory, Algorithms and Applications' ***
31  *** (grant RBAU01877P), http://slipguru.disi.unige.it/ASTA ***
32  *** and ***
33  *** 'Parallel Algorithms and Numerical Nonlinear Optimization' ***
34  *** (grant RBAU01JYPN), http://dm.unife.it/pn2o ***
35  *** ***
36  *** Copyright (C) 2004-2008 by T. Serafini, G. Zanghirati, L. Zanni. ***
37  *** ***
38  *** COPYRIGHT NOTIFICATION ***
39  *** ***
40  *** Permission to copy and modify this software and its documentation ***
41  *** for internal research use is granted, provided that this notice is ***
42  *** retained thereon and on all copies or modifications. The authors and ***
43  *** their respective Universities makes no representations as to the ***
44  *** suitability and operability of this software for any purpose. It is ***
45  *** provided "as is" without express or implied warranty. ***
46  *** Use of this software for commercial purposes is expressly prohibited ***
47  *** without contacting the authors. ***
48  *** ***
49  *** This program is free software; you can redistribute it and/or modify ***
50  *** it under the terms of the GNU General Public License as published by ***
51  *** the Free Software Foundation; either version 3 of the License, or ***
52  *** (at your option) any later version. ***
53  *** ***
54  *** This program is distributed in the hope that it will be useful, ***
55  *** but WITHOUT ANY WARRANTY; without even the implied warranty of ***
56  *** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ***
57  *** GNU General Public License for more details. ***
58  *** ***
59  *** You should have received a copy of the GNU General Public License ***
60  *** along with this program; if not, write to the Free Software ***
61  *** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. ***
62  *** ***
63  *** File: gpm.cpp ***
64  *** Type: scalar ***
65  *** Version: 1.0 ***
66  *** Date: October, 2005 ***
67  *** Revision: 1 ***
68  *** ***
69  *** SHOGUN adaptions Written (W) 2006-2008 Soeren Sonnenburg ***
70  ******************************************************************************/
71 #include <stdio.h>
72 #include <stdlib.h>
73 #include <string.h>
74 #include <math.h>
75 #include <shogun/lib/external/gpdt.h>
76 #include <shogun/io/SGIO.h>
77 
78 namespace shogun
79 {
80 #define maxvpm 30000 /* max number of method iterations allowed */
81 #define maxprojections 200
82 #define in 8000 /* max size of the QP problem to solve */
83 #define alpha_max 1e10
84 #define alpha_min 1e-10
85 
86 extern uint32_t Randnext;
87 #define ThRand (Randnext = Randnext * 1103515245L + 12345L)
88 #define ThRandPos ((Randnext = Randnext * 1103515245L + 12345L) & 0x7fffffff)
89 
90 int32_t InnerProjector(
91  int32_t method, int32_t n, int32_t *iy, float64_t e, float64_t *qk, float64_t l,
92  float64_t u, float64_t *x, float64_t &lambda);
93 
94 /* Uncomment if you want to allocate vectors on the stack *
95  * instead of the heap. On some architectures this helps *
96  * improving speed, but may generate a stack overflow */
97 // #define VARIABLES_ON_STACK
98 
99 /* Uncomment if you want to use the adaptive steplength *
100  in the GVPM solver */
101 #define VPM_ADA
102 
103 
104 /******************************************************************************
105  *** Generalized Variable Projection Method (T. Serafini, G. Zanghirati, ***
106  *** L. Zanni, "Gradient Projection Methods for Quadratic Programs and ***
107  *** Applications in Training Support Vector Machines"; ***
108  *** Optim. Meth. Soft. 20, 2005, 353-378) ***
109  ******************************************************************************/
110 int32_t gvpm(
111  int32_t Projector, int32_t n, float32_t *vecA, float64_t *b, float64_t c,
112  float64_t e, int32_t *iy, float64_t *x, float64_t tol, int32_t *ls,
113  int32_t *proj)
114 {
115  int32_t i, j, iter, it, it2, luv, info;
116  float64_t gd, max, normd, dAd, lam, lamnew, alpha, kktlam, ak, bk;
117 
118  int32_t lscount = 0, projcount = 0;
119  float64_t eps = 1.0e-16;
120  float64_t DELTAsv, ProdDELTAsv;
121  float64_t lam_ext;
122 
123  /* solver-specific settings */
124 #ifdef VPM_ADA
125  int32_t nc = 1, ia1 = -1;
126  float64_t alpha1, alpha2;
127 #endif
128 
129  /* allocation-dependant settings */
130 #ifdef VARIABLES_ON_STACK
131 
132  int32_t ipt[in], ipt2[in], uv[in];
133  float64_t g[in], y[in], tempv[in], d[in], Ad[in], t[in];
134 
135 #else
136 
137  int32_t *ipt, *ipt2, *uv;
138  float64_t *g, *y, *tempv, *d, *Ad, *t;
139 
140  /*** array allocations ***/
141  ipt = SG_MALLOC(int32_t, n);
142  ipt2 = SG_MALLOC(int32_t, n);
143  uv = SG_MALLOC(int32_t, n);
144  g = SG_MALLOC(float64_t, n);
145  y = SG_MALLOC(float64_t, n);
146  d = SG_MALLOC(float64_t, n);
147  Ad = SG_MALLOC(float64_t, n);
148  t = SG_MALLOC(float64_t, n);
149  tempv = SG_MALLOC(float64_t, n);
150 
151 #endif
152 
153  DELTAsv = EPS_SV * c;
154  if (tol <= 1.0e-5 || n <= 20)
155  ProdDELTAsv = 0.0F;
156  else
157  ProdDELTAsv = EPS_SV * c;
158 
159  for (i = 0; i < n; i++)
160  tempv[i] = -x[i];
161  lam_ext = 0.0;
162  projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, x, lam_ext);
163 
164  /* compute g = A*x + b in sparse form *
165  * (inline computation for better perfomrance) */
166  {
167  float32_t *tempA;
168 
169  it = 0;
170  for (i = 0; i < n; i++)
171  if (fabs(x[i]) > ProdDELTAsv*1e-2)
172  ipt[it++] = i;
173 
174  memset(t, 0, n*sizeof(float64_t));
175  for (i = 0; i < it; i++)
176  {
177  tempA = vecA + ipt[i]*n;
178  for (j = 0; j < n; j++)
179  t[j] += (tempA[j] * x[ipt[i]]);
180  }
181  }
182 
183  for (i = 0; i < n; i++)
184  {
185  g[i] = t[i] + b[i],
186  y[i] = g[i] - x[i];
187  }
188 
189  projcount += InnerProjector(Projector, n, iy, e, y, 0, c, tempv, lam_ext);
190 
191  max = alpha_min;
192  for (i = 0; i < n; i++)
193  {
194  y[i] = tempv[i] - x[i];
195  if (fabs(y[i]) > max)
196  max = fabs(y[i]);
197  }
198 
199  if (max < c*tol*1e-3)
200  {
201  lscount = 0;
202  iter = 0;
203  goto Clean;
204  }
205 
206  alpha = 1.0 / max;
207 
208  for (iter = 1; iter <= maxvpm; iter++)
209  {
210  for (i = 0; i < n; i++)
211  tempv[i] = alpha*g[i] - x[i];
212 
213  projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, y, lam_ext);
214 
215  gd = 0.0;
216  for (i = 0; i < n; i++)
217  {
218  d[i] = y[i] - x[i];
219  gd += d[i] * g[i];
220  }
221 
222  /* compute Ad = A*d or Ad = Ay-t depending on their sparsity *
223  * (inline computation for better perfomrance) */
224  {
225  float32_t *tempA;
226 
227  it = it2 = 0;
228  for (i = 0; i < n; i++)
229  if (fabs(d[i]) > (ProdDELTAsv*1.0e-2))
230  ipt[it++] = i;
231  for (i = 0; i < n; i++)
232  if (fabs(y[i]) > ProdDELTAsv)
233  ipt2[it2++] = i;
234 
235  memset(Ad, 0, n*sizeof(float64_t));
236  if (it < it2) // Ad = A*d
237  {
238  for (i = 0; i < it; i++)
239  {
240  tempA = vecA + ipt[i]*n;
241  for (j = 0; j < n; j++)
242  Ad[j] += (tempA[j] * d[ipt[i]]);
243  }
244  }
245  else // Ad = A*y - t
246  {
247  for (i = 0; i < it2; i++)
248  {
249  tempA = vecA + ipt2[i]*n;
250  for (j = 0; j < n; j++)
251  Ad[j] += (tempA[j] * y[ipt2[i]]);
252  }
253  for (j = 0; j < n; j++)
254  Ad[j] -= t[j];
255  }
256  }
257 
258  normd = 0.0;
259  for (i = 0; i < n; i++)
260  normd += d[i] * d[i];
261 
262  dAd = 0.0;
263  for (i = 0; i < n; i++)
264  dAd += d[i]*Ad[i];
265 
266  if (dAd > eps*normd && gd < 0.0)
267  {
268  lam = lamnew = -gd/dAd;
269  if (lam > 1.0 || lam < 0.0)
270  lam = 1.0;
271  else
272  lscount++;
273 
274 #ifdef VPM_ADA
275 
276  /*** use the adaptive switching rule for steplength selection ***/
277 
278  // compute alpha1 = (d'* (d.*diaga)) / (d'*Ad);
279  alpha1 = normd / dAd;
280 
281  // alpha2 = d'*Ad / (Ad' * (Ad./diaga));
282  alpha2 = 0.0;
283  for (i = 0; i < n; i++)
284  alpha2 += Ad[i] * Ad[i];
285  alpha2 = dAd / alpha2;
286 
287  if ( (nc > 2
288  && (
289  (ia1 == 1
290  && (
291  lamnew < 0.1 || (alpha1 > alpha && alpha2 < alpha)
292  )
293  )
294  ||
295  (ia1 == -1
296  && (
297  lamnew > 5.0 || (alpha1 > alpha && alpha2 < alpha)
298  )
299  )
300  )
301  )
302  || nc > 9 )
303  {
304  ia1 = -ia1;
305  nc = 0;
306  }
307 
308  if (ia1 == 1)
309  alpha = alpha1;
310  else
311  alpha = alpha2;
312 
313  if (alpha < alpha_min)
314  alpha = alpha_min;
315  else if (alpha > alpha_max)
316  alpha = alpha_max;
317 
318  nc++;
319 
320 #else
321 
322  /*** use the fixed switching rule for steplength selection ***/
323 
324  if ((iter % 6) < 3) // alpha = d'*Ad / (Ad' * (Ad./diaga));
325  {
326  alpha = 0.0;
327  for (i = 0; i < n; i++)
328  alpha += Ad[i] * Ad[i];
329  alpha = dAd / alpha;
330  }
331  else // alpha = (d'* (d.*diaga)) / (d'*Ad);
332  {
333  alpha = 0.0;
334  for (i = 0; i < n; i++)
335  alpha += d[i] * d[i];
336  alpha = alpha / dAd;
337  }
338 
339 #endif
340 
341  }
342  else
343  {
344  lam = 1.0;
345  alpha = alpha_max;
346  }
347 
348  for (i = 0; i < n; i++)
349  {
350  x[i] = x[i] + lam*d[i];
351  t[i] = t[i] + lam*Ad[i];
352  g[i] = b[i] + t[i];
353  }
354 
355  /*** stopping criterion based on KKT conditions ***/
356  bk = 0.0;
357  ak = 0.0;
358  for (i = 0; i < n; i++)
359  {
360  bk += x[i] * x[i];
361  ak += d[i] * d[i];
362  }
363 
364  if (lam*sqrt(ak) < tol*10 * sqrt(bk))
365  {
366  it = 0;
367  luv = 0;
368  kktlam = 0.0;
369  for (i = 0; i < n; i++)
370  {
371  if (x[i] > DELTAsv && x[i] < c-DELTAsv)
372  {
373  ipt[it++] = i;
374  kktlam = kktlam - iy[i]*g[i];
375  }
376  else
377  uv[luv++] = i;
378  }
379 
380  if (it == 0)
381  {
382  if (lam*sqrt(ak) < tol*0.5 * sqrt(bk))
383  goto Clean;
384  }
385  else
386  {
387  kktlam = kktlam/it;
388  info = 1;
389  for (i = 0; i < it; i++)
390  if (fabs(iy[ipt[i]]*g[ipt[i]]+kktlam) > tol)
391  {
392  info = 0;
393  break;
394  }
395 
396  if (info == 1)
397  for (i = 0; i < luv; i++)
398  if (x[uv[i]] <= DELTAsv)
399  {
400  if (g[uv[i]] + kktlam*iy[uv[i]] < -tol)
401  {
402  info = 0;
403  break;
404  }
405  }
406  else
407  {
408  if (g[uv[i]] + kktlam*iy[uv[i]] > tol)
409  {
410  info = 0;
411  break;
412  }
413  }
414 
415  if (info == 1)
416  goto Clean;
417  }
418  } // stopping rule based on the norm of d_k
419  }
420 
421  SG_SWARNING("GVPM exits after maxvpm = %d iterations.\n", maxvpm);
422 
423 Clean:
424 
425  /*** allocation-dependant freeing ***/
426 #ifndef VARIABLES_ON_STACK
427  SG_FREE(t);
428  SG_FREE(uv);
429  SG_FREE(ipt2);
430  SG_FREE(ipt);
431  SG_FREE(g);
432  SG_FREE(y);
433  SG_FREE(tempv);
434  SG_FREE(d);
435  SG_FREE(Ad);
436 #endif
437 
438  if (ls != NULL) *ls = lscount;
439  if (proj != NULL) *proj = projcount;
440  return(iter);
441 }
442 
443 /******************************************************************************
444  *** Dai-Fletcher QP solver (Y. Dai and R. Fletcher,"New Algorithms for ***
445  *** Singly Linear Constrained Quadratic Programs Subject to Lower and ***
446  *** Upper Bounds"; Math. Prog. to appear) ***
447  ******************************************************************************/
448 int32_t FletcherAlg2A(
449  int32_t Projector, int32_t n, float32_t *vecA, float64_t *b, float64_t c,
450  float64_t e, int32_t *iy, float64_t *x, float64_t tol, int32_t *ls,
451  int32_t *proj)
452 {
453  int32_t i, j, iter, it, it2, luv, info, lscount = 0, projcount = 0;
454  float64_t gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam, lam_ext;
455  float64_t eps = 1.0e-16;
456  float64_t DELTAsv, ProdDELTAsv;
457 
458  /*** variables for the adaptive nonmonotone linesearch ***/
459  int32_t L, llast;
460  float64_t fr, fbest, fv, fc, fv0;
461 
462  /*** allocation-dependant settings ***/
463 #ifdef VARIABLES_ON_STACK
464 
465  int32_t ipt[in], ipt2[in], uv[in];
466  float64_t g[in], y[in], tempv[in], d[in], Ad[in], t[in],
467  xplus[in], tplus[in], sk[in], yk[in];
468 #else
469 
470  int32_t *ipt, *ipt2, *uv;
471  float64_t *g, *y, *tempv, *d, *Ad, *t, *xplus, *tplus, *sk, *yk;
472 
473  /*** arrays allocation ***/
474  ipt = SG_MALLOC(int32_t, n);
475  ipt2 = SG_MALLOC(int32_t, n);
476  uv = SG_MALLOC(int32_t, n);
477  g = SG_MALLOC(float64_t, n);
478  y = SG_MALLOC(float64_t, n);
479  tempv = SG_MALLOC(float64_t, n);
480  d = SG_MALLOC(float64_t, n);
481  Ad = SG_MALLOC(float64_t, n);
482  t = SG_MALLOC(float64_t, n);
483  xplus = SG_MALLOC(float64_t, n);
484  tplus = SG_MALLOC(float64_t, n);
485  sk = SG_MALLOC(float64_t, n);
486  yk = SG_MALLOC(float64_t, n);
487 
488 #endif
489 
490  DELTAsv = EPS_SV * c;
491  if (tol <= 1.0e-5 || n <= 20)
492  ProdDELTAsv = 0.0F;
493  else
494  ProdDELTAsv = EPS_SV * c;
495 
496  for (i = 0; i < n; i++)
497  tempv[i] = -x[i];
498 
499  lam_ext = 0.0;
500 
501  projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, x, lam_ext);
502 
503  // g = A*x + b;
504  // SparseProd(n, t, A, x, ipt);
505  {
506  float32_t *tempA;
507 
508  it = 0;
509  for (i = 0; i < n; i++)
510  if (fabs(x[i]) > ProdDELTAsv)
511  ipt[it++] = i;
512 
513  memset(t, 0, n*sizeof(float64_t));
514  for (i = 0; i < it; i++)
515  {
516  tempA = vecA + ipt[i] * n;
517  for (j = 0; j < n; j++)
518  t[j] += (tempA[j]*x[ipt[i]]);
519  }
520  }
521 
522  for (i = 0; i < n; i++)
523  {
524  g[i] = t[i] + b[i],
525  y[i] = g[i] - x[i];
526  }
527 
528  projcount += InnerProjector(Projector, n, iy, e, y, 0, c, tempv, lam_ext);
529 
530  max = alpha_min;
531  for (i = 0; i < n; i++)
532  {
533  y[i] = tempv[i] - x[i];
534  if (fabs(y[i]) > max)
535  max = fabs(y[i]);
536  }
537 
538  if (max < c*tol*1e-3)
539  {
540  lscount = 0;
541  iter = 0;
542  goto Clean;
543  }
544 
545  alpha = 1.0 / max;
546 
547  fv0 = 0.0;
548  for (i = 0; i < n; i++)
549  fv0 += x[i] * (0.5*t[i] + b[i]);
550 
551  /*** adaptive nonmonotone linesearch ***/
552  L = 2;
553  fr = alpha_max;
554  fbest = fv0;
555  fc = fv0;
556  llast = 0;
557  akold = bkold = 0.0;
558 
559  for (iter = 1; iter <= maxvpm; iter++)
560  {
561  for (i = 0; i < n; i++)
562  tempv[i] = alpha*g[i] - x[i];
563 
564  projcount += InnerProjector(Projector, n, iy, e, tempv, 0, c, y, lam_ext);
565 
566  gd = 0.0;
567  for (i = 0; i < n; i++)
568  {
569  d[i] = y[i] - x[i];
570  gd += d[i] * g[i];
571  }
572 
573  /* compute Ad = A*d or Ad = A*y - t depending on their sparsity */
574  {
575  float32_t *tempA;
576 
577  it = it2 = 0;
578  for (i = 0; i < n; i++)
579  if (fabs(d[i]) > (ProdDELTAsv*1.0e-2))
580  ipt[it++] = i;
581  for (i = 0; i < n; i++)
582  if (fabs(y[i]) > ProdDELTAsv)
583  ipt2[it2++] = i;
584 
585  memset(Ad, 0, n*sizeof(float64_t));
586  if (it < it2) // compute Ad = A*d
587  {
588  for (i = 0; i < it; i++)
589  {
590  tempA = vecA + ipt[i]*n;
591  for (j = 0; j < n; j++)
592  Ad[j] += (tempA[j] * d[ipt[i]]);
593  }
594  }
595  else // compute Ad = A*y-t
596  {
597  for (i = 0; i < it2; i++)
598  {
599  tempA = vecA + ipt2[i]*n;
600  for (j = 0; j < n; j++)
601  Ad[j] += (tempA[j] * y[ipt2[i]]);
602  }
603  for (j = 0; j < n; j++)
604  Ad[j] -= t[j];
605  }
606  }
607 
608  ak = 0.0;
609  for (i = 0; i < n; i++)
610  ak += d[i] * d[i];
611 
612  bk = 0.0;
613  for (i = 0; i < n; i++)
614  bk += d[i]*Ad[i];
615 
616  if (bk > eps*ak && gd < 0.0) // ak is normd
617  lamnew = -gd/bk;
618  else
619  lamnew = 1.0;
620 
621  fv = 0.0;
622  for (i = 0; i < n; i++)
623  {
624  xplus[i] = x[i] + d[i];
625  tplus[i] = t[i] + Ad[i];
626  fv += xplus[i] * (0.5*tplus[i] + b[i]);
627  }
628 
629  if ((iter == 1 && fv >= fv0) || (iter > 1 && fv >= fr))
630  {
631  lscount++;
632  fv = 0.0;
633  for (i = 0; i < n; i++)
634  {
635  xplus[i] = x[i] + lamnew*d[i];
636  tplus[i] = t[i] + lamnew*Ad[i];
637  fv += xplus[i] * (0.5*tplus[i] + b[i]);
638  }
639  }
640 
641  for (i = 0; i < n; i++)
642  {
643  sk[i] = xplus[i] - x[i];
644  yk[i] = tplus[i] - t[i];
645  x[i] = xplus[i];
646  t[i] = tplus[i];
647  g[i] = t[i] + b[i];
648  }
649 
650  // update the line search control parameters
651 
652  if (fv < fbest)
653  {
654  fbest = fv;
655  fc = fv;
656  llast = 0;
657  }
658  else
659  {
660  fc = (fc > fv ? fc : fv);
661  llast++;
662  if (llast == L)
663  {
664  fr = fc;
665  fc = fv;
666  llast = 0;
667  }
668  }
669 
670  ak = bk = 0.0;
671  for (i = 0; i < n; i++)
672  {
673  ak += sk[i] * sk[i];
674  bk += sk[i] * yk[i];
675  }
676 
677  if (bk < eps*ak)
678  alpha = alpha_max;
679  else
680  {
681  if (bkold < eps*akold)
682  alpha = ak/bk;
683  else
684  alpha = (akold+ak)/(bkold+bk);
685 
686  if (alpha > alpha_max)
687  alpha = alpha_max;
688  else if (alpha < alpha_min)
689  alpha = alpha_min;
690  }
691 
692  akold = ak;
693  bkold = bk;
694 
695  /*** stopping criterion based on KKT conditions ***/
696 
697  bk = 0.0;
698  for (i = 0; i < n; i++)
699  bk += x[i] * x[i];
700 
701  if (sqrt(ak) < tol*10 * sqrt(bk))
702  {
703  it = 0;
704  luv = 0;
705  kktlam = 0.0;
706  for (i = 0; i < n; i++)
707  {
708  if ((x[i] > DELTAsv) && (x[i] < c-DELTAsv))
709  {
710  ipt[it++] = i;
711  kktlam = kktlam - iy[i]*g[i];
712  }
713  else
714  uv[luv++] = i;
715  }
716 
717  if (it == 0)
718  {
719  if (sqrt(ak) < tol*0.5 * sqrt(bk))
720  goto Clean;
721  }
722  else
723  {
724 
725  kktlam = kktlam/it;
726  info = 1;
727  for (i = 0; i < it; i++)
728  if ( fabs(iy[ipt[i]] * g[ipt[i]] + kktlam) > tol )
729  {
730  info = 0;
731  break;
732  }
733 
734  if (info == 1)
735  for (i = 0; i < luv; i++)
736  if (x[uv[i]] <= DELTAsv)
737  {
738  if (g[uv[i]] + kktlam*iy[uv[i]] < -tol)
739  {
740  info = 0;
741  break;
742  }
743  }
744  else
745  {
746  if (g[uv[i]] + kktlam*iy[uv[i]] > tol)
747  {
748  info = 0;
749  break;
750  }
751  }
752 
753  if (info == 1)
754  goto Clean;
755  }
756  }
757  }
758 
759  SG_SWARNING("Dai-Fletcher method exits after maxvpm = %d iterations.\n", maxvpm);
760 
761 Clean:
762 
763 #ifndef VARIABLES_ON_STACK
764  SG_FREE(sk);
765  SG_FREE(yk);
766  SG_FREE(tplus);
767  SG_FREE(xplus);
768  SG_FREE(t);
769  SG_FREE(uv);
770  SG_FREE(ipt2);
771  SG_FREE(ipt);
772  SG_FREE(g);
773  SG_FREE(y);
774  SG_FREE(tempv);
775  SG_FREE(d);
776  SG_FREE(Ad);
777 #endif
778 
779  if (ls != NULL) *ls = lscount;
780  if (proj != NULL) *proj = projcount;
781  return(iter);
782 
783 }
784 
785 /******************************************************************************/
786 /*** Encapsulating method to call the chosen Gradient Projection Method ***/
787 /******************************************************************************/
788 int32_t gpm_solver(
789  int32_t Solver, int32_t Projector, int32_t n, float32_t *A, float64_t *b,
790  float64_t c, float64_t e, int32_t *iy, float64_t *x, float64_t tol,
791  int32_t *ls, int32_t *proj)
792 {
793  /*** Uncomment the following if you need to scale the QP Hessian matrix
794  *** before calling the chosen solver
795  int32_t i, j;
796  float32_t *ptrA;
797  float64_t max, s;
798 
799  max = fabs(A[0][0]);
800  for (i = 1; i < n; i++)
801  if (fabs(A[i][i]) > max)
802  max = fabs(A[i][i]);
803 
804  s = 1.0 / max;
805  ptrA = vecA;
806  for (i = 0; i < n; i++)
807  for (j = 0;j < n; j++)
808  *ptrA++ = (float32_t)(A[i][j]*s);
809 
810  if (Solver == SOLVER_FLETCHER)
811  j = FletcherAlg2A(n, vecA, b, c/s, e/s, iy, x, tol, ls);
812  else
813  j = gvpm(n, vecA, b, c/s, e/s, iy, x, tol, ls);
814 
815  for (i = 0; i < n; i++)
816  x[i] *= s;
817  ***/
818 
819  /*** calling the chosen solver with unscaled data ***/
820  if (Solver == SOLVER_FLETCHER)
821  return FletcherAlg2A(Projector, n, A, b, c, e, iy, x, tol, ls, proj);
822  else
823  return gvpm(Projector, n, A, b, c, e, iy, x, tol, ls, proj);
824 }
825 
826 /******************************************************************************
827  *** Piecewise linear monotone target function for the Dai-Fletcher ***
828  *** projector (Y. Dai and R. Fletcher, "New Algorithms for Singly Linear ***
829  *** Constrained Quadratic Programs Subject to Lower and Upper Bounds"; ***
830  *** Math. Prog. to appear) ***
831  ******************************************************************************/
832 float64_t ProjectR(
833  float64_t *x, int32_t n, float64_t lambda, int32_t *a, float64_t b,
834  float64_t *c, float64_t l, float64_t u)
835 {
836  int32_t i;
837  float64_t r = 0.0;
838 
839  for (i = 0; i < n; i++)
840  {
841  x[i] = -c[i] + lambda*(float64_t)a[i];
842  if (x[i] >= u) x[i] = u;
843  else if (x[i] < l) x[i] = l;
844  r += (float64_t)a[i]*x[i];
845  }
846 
847  return (r - b);
848 }
849 
850 /******************************************************************************
851  *** Dai-Fletcher QP projector (Y. Dai and R. Fletcher, "New Algorithms for ***
852  *** Singly Linear Constrained Quadratic Programs Subject to Lower and ***
853  *** Upper Bounds"; Math. Prog. to appear) ***
854  ******************************************************************************/
855 /*** ***
856  *** Solves the problem min x'*x/2 + c'*x ***
857  *** subj to a'*x - b = 0 ***
858  *** l <= x <= u ***
859  ******************************************************************************/
860 int32_t ProjectDai(
861  int32_t n, int32_t *a, float64_t b, float64_t *c, float64_t l, float64_t u,
862  float64_t *x, float64_t &lam_ext)
863 {
864  float64_t lambda, lambdal, lambdau, dlambda, lambda_new, tol_lam;
865  float64_t r, rl, ru, s, tol_r;
866  int32_t iter;
867 
868  tol_lam = 1.0e-11;
869  tol_r = 1.0e-10 * sqrt((u-l)*(float64_t)n);
870  lambda = lam_ext;
871  dlambda = 0.5;
872  iter = 1;
873  b = -b;
874 
875  // Bracketing Phase
876  r = ProjectR(x, n, lambda, a, b, c, l, u);
877  if (fabs(r) < tol_r)
878  return 0;
879 
880  if (r < 0.0)
881  {
882  lambdal = lambda;
883  rl = r;
884  lambda = lambda + dlambda;
885  r = ProjectR(x, n, lambda, a, b, c, l, u);
886  while (r < 0.0)
887  {
888  lambdal = lambda;
889  s = rl/r - 1.0;
890  if (s < 0.1) s = 0.1;
891  dlambda = dlambda + dlambda/s;
892  lambda = lambda + dlambda;
893  rl = r;
894  r = ProjectR(x, n, lambda, a, b, c, l, u);
895  }
896  lambdau = lambda;
897  ru = r;
898  }
899  else
900  {
901  lambdau = lambda;
902  ru = r;
903  lambda = lambda - dlambda;
904  r = ProjectR(x, n, lambda, a, b, c, l, u);
905  while (r > 0.0)
906  {
907  lambdau = lambda;
908  s = ru/r - 1.0;
909  if (s < 0.1) s = 0.1;
910  dlambda = dlambda + dlambda/s;
911  lambda = lambda - dlambda;
912  ru = r;
913  r = ProjectR(x, n, lambda, a, b, c, l, u);
914  }
915  lambdal = lambda;
916  rl = r;
917  }
918 
919 
920  // Secant Phase
921  s = 1.0 - rl/ru;
922  dlambda = dlambda/s;
923  lambda = lambdau - dlambda;
924  r = ProjectR(x, n, lambda, a, b, c, l, u);
925 
926  while ( fabs(r) > tol_r
927  && dlambda > tol_lam * (1.0 + fabs(lambda))
928  && iter < maxprojections )
929  {
930  iter++;
931  if (r > 0.0)
932  {
933  if (s <= 2.0)
934  {
935  lambdau = lambda;
936  ru = r;
937  s = 1.0 - rl/ru;
938  dlambda = (lambdau - lambdal) / s;
939  lambda = lambdau - dlambda;
940  }
941  else
942  {
943  s = ru/r-1.0;
944  if (s < 0.1) s = 0.1;
945  dlambda = (lambdau - lambda) / s;
946  lambda_new = 0.75*lambdal + 0.25*lambda;
947  if (lambda_new < (lambda - dlambda))
948  lambda_new = lambda - dlambda;
949  lambdau = lambda;
950  ru = r;
951  lambda = lambda_new;
952  s = (lambdau - lambdal) / (lambdau - lambda);
953  }
954  }
955  else
956  {
957  if (s >= 2.0)
958  {
959  lambdal = lambda;
960  rl = r;
961  s = 1.0 - rl/ru;
962  dlambda = (lambdau - lambdal) / s;
963  lambda = lambdau - dlambda;
964  }
965  else
966  {
967  s = rl/r - 1.0;
968  if (s < 0.1) s = 0.1;
969  dlambda = (lambda-lambdal) / s;
970  lambda_new = 0.75*lambdau + 0.25*lambda;
971  if (lambda_new > (lambda + dlambda))
972  lambda_new = lambda + dlambda;
973  lambdal = lambda;
974  rl = r;
975  lambda = lambda_new;
976  s = (lambdau - lambdal) / (lambdau-lambda);
977  }
978  }
979  r = ProjectR(x, n, lambda, a, b, c, l, u);
980  }
981 
982  lam_ext = lambda;
983  if (iter >= maxprojections)
984  SG_SERROR("Projector exits after max iterations: %d\n", iter);
985 
986  return (iter);
987 }
988 
989 #define SWAP(a,b) { register float64_t t=(a);(a)=(b);(b)=t; }
990 
991 /*** Median computation using Quick Select ***/
992 float64_t quick_select(float64_t *arr, int32_t n)
993 {
994  int32_t low, high;
995  int32_t median;
996  int32_t middle, l, h;
997 
998  low = 0;
999  high = n-1;
1000  median = (low + high) / 2;
1001 
1002  for (;;)
1003  {
1004  if (high <= low)
1005  return arr[median];
1006 
1007  if (high == low + 1)
1008  {
1009  if (arr[low] > arr[high])
1010  SWAP(arr[low], arr[high]);
1011  return arr[median];
1012  }
1013 
1014  middle = (low + high) / 2;
1015  if (arr[middle] > arr[high]) SWAP(arr[middle], arr[high]);
1016  if (arr[low] > arr[high]) SWAP(arr[low], arr[high]);
1017  if (arr[middle] > arr[low]) SWAP(arr[middle], arr[low]);
1018 
1019  SWAP(arr[middle], arr[low+1]);
1020 
1021  l = low + 1;
1022  h = high;
1023  for (;;)
1024  {
1025  do l++; while (arr[low] > arr[l]);
1026  do h--; while (arr[h] > arr[low]);
1027  if (h < l)
1028  break;
1029  SWAP(arr[l], arr[h]);
1030  }
1031 
1032  SWAP(arr[low], arr[h]);
1033  if (h <= median)
1034  low = l;
1035  if (h >= median)
1036  high = h - 1;
1037  }
1038 }
1039 
1040 /******************************************************************************
1041  *** Pardalos-Kovoor projector (P.M. Pardalos and N. Kovoor, "An Algorithm ***
1042  *** for a Singly Constrained Class of Quadratic Programs Subject to Upper ***
1043  *** and Lower Bounds"; Math. Prog. 46, 1990, 321-328). ***
1044  ******************************************************************************
1045  *** Solves the problem ***
1046  *** min x'*x/2 + qk'*x ***
1047  *** subj to iy'*x + e = 0 ***
1048  *** l <= x <= u ***
1049  *** iy(i) ~= 0 ***
1050  ******************************************************************************/
1051 
1052 int32_t Pardalos(
1053  int32_t n, int32_t *iy, float64_t e, float64_t *qk, float64_t low,
1054  float64_t up, float64_t *x)
1055 {
1056  int32_t i, l, iter; /* conters */
1057  int32_t luv, lxint; /* dimensions */
1058  float64_t d, xmin, xmax, xmold, xmid, xx, ts, sw, s, s1, testsum;
1059 
1060  /*** allocation-dependant settings ***/
1061 #ifdef VARIABLES_ON_STACK
1062  int32_t uv[in], uvt[in];
1063  float64_t xint[2*in+2], xint2[2*in+2], a[in], b[in], at[in], bt[in];
1064  float64_t newdia[in], newdt[in];
1065 #else
1066 
1067  int32_t *uv, *uvt;
1068  float64_t *xint, *xint2, *a, *b, *at, *bt, *newdia, *newdt;
1069 
1070  /*** arrays allocation ***/
1071  uv = SG_MALLOC(int32_t, n);
1072  uvt = SG_MALLOC(int32_t, n);
1073  a = SG_MALLOC(float64_t, n);
1074  b = SG_MALLOC(float64_t, n);
1075  at = SG_MALLOC(float64_t, n);
1076  bt = SG_MALLOC(float64_t, n);
1077  newdia = SG_MALLOC(float64_t, n);
1078  newdt = SG_MALLOC(float64_t, n);
1079  xint = SG_MALLOC(float64_t, (2*n + 2));
1080  xint2 = SG_MALLOC(float64_t, (2*n + 2));
1081 
1082 #endif
1083 
1084  d = 0.0;
1085  for (i = 0; i < n; i++)
1086  d += iy[i] * qk[i];
1087  d = 0.5 * (d-e);
1088 
1089  for (i = 0; i < n; i++)
1090  {
1091  /* The following computations should divide by iy[i] instead *
1092  * of multiply by iy[i], but this is correct for binary classification *
1093  * with labels -1 and 1. */
1094  if (iy[i] > 0)
1095  {
1096  a[i] = ((qk[i] + low) * iy[i]) * 0.5;
1097  b[i] = ((up + qk[i]) * iy[i]) * 0.5;
1098  }
1099  else
1100  {
1101  b[i] = ((qk[i] + low) * iy[i]) * 0.5;
1102  a[i] = ((up + qk[i]) * iy[i]) * 0.5;
1103  }
1104  newdia[i] = (iy[i]*iy[i]);
1105  }
1106 
1107  xmin = -1e33;
1108  xmax = 1e33;
1109 
1110  /* arrays initialization */
1111  for (i = 0; i < n; i++)
1112  {
1113  uv[i] = i; /* contains the unset variables */
1114  xint[i] = a[i];
1115  xint[n+i] = b[i];
1116  }
1117 
1118  xmid = xmin;
1119  xint[2*n] = xmin;
1120  xint[2*n+1] = xmax;
1121  ts = 0.0;
1122  sw = 0.0;
1123  luv = n;
1124  lxint = 2*n+2;
1125 
1126  iter = 0;
1127  do {
1128  for (i = 0; i < luv; i++)
1129  {
1130  at[i] = a[uv[i]];
1131  bt[i] = b[uv[i]];
1132  newdt[i] = newdia[uv[i]];
1133  }
1134 
1135  xmold = xmid;
1136  xmid = quick_select(xint, lxint);
1137  if (xmold == xmid)
1138  xmid = xint[(int32_t)(ThRandPos % lxint)];
1139 
1140  s = ts;
1141  s1 = sw;
1142  for (i = 0; i < luv; i++)
1143  {
1144  if (xmid > bt[i])
1145  s += newdt[i]*bt[i];
1146  else if (xmid < at[i])
1147  s += newdt[i]*at[i];
1148  else
1149  s1 += newdt[i];
1150  }
1151 
1152  testsum = s + s1*xmid;
1153  if (testsum <= (d+(1e-15)))
1154  xmin = xmid;
1155  if (testsum >= (d-(1e-15)))
1156  xmax = xmid;
1157 
1158  l = 0;
1159  for (i = 0; i < lxint; i++)
1160  if((xint[i] >= xmin) && (xint[i] <= xmax))
1161  xint2[l++] = xint[i];
1162  lxint = l;
1163  memcpy(xint, xint2, lxint*sizeof(float64_t));
1164 
1165  l = 0;
1166  for (i = 0; i < luv; i++)
1167  {
1168  if (xmin >= bt[i])
1169  ts += newdt[i]*bt[i];
1170  else if (xmax <= at[i])
1171  ts += newdt[i]*at[i];
1172  else if ((at[i] <= xmin) && (bt[i] >= xmax))
1173  sw += newdt[i];
1174  else
1175  uvt[l++] = uv[i];
1176  }
1177  luv = l;
1178  memcpy(uv, uvt, luv*sizeof(int32_t));
1179  iter++;
1180  } while(luv != 0 && iter < maxprojections);
1181 
1182  if (sw == 0)
1183  xx = xmin;
1184  else
1185  xx = (d-ts) / sw;
1186 
1187  for (i = 0; i < n; i++)
1188  {
1189  if (b[i] <= xmin)
1190  x[i] = b[i];
1191  else if (a[i] >= xmax)
1192  x[i] = a[i];
1193  else if ((a[i]<=xmin) && (xmax<=b[i]))
1194  x[i] = xx;
1195  else
1196  SG_SWARNING("Inner solver troubles...\n");
1197  }
1198 
1199  for (i = 0; i < n; i++)
1200  x[i] = (2.0*x[i]*iy[i]-qk[i]);
1201 
1202 #ifndef VARIABLES_ON_STACK
1203  SG_FREE(newdt);
1204  SG_FREE(newdia);
1205  SG_FREE(a);
1206  SG_FREE(b);
1207  SG_FREE(uvt);
1208  SG_FREE(uv);
1209  SG_FREE(bt);
1210  SG_FREE(at);
1211  SG_FREE(xint2);
1212  SG_FREE(xint);
1213 #endif
1214 
1215  return(iter);
1216 }
1217 
1218 /******************************************************************************/
1219 /*** Wrapper method to call the selected inner projector ***/
1220 /******************************************************************************/
1221 int32_t InnerProjector(
1222  int32_t method, int32_t n, int32_t *iy, float64_t e, float64_t *qk,
1223  float64_t l, float64_t u, float64_t *x, float64_t &lambda)
1224 {
1225  if (method == 0)
1226  return Pardalos(n, iy, e, qk, l, u, x);
1227  else
1228  return ProjectDai(n, iy, e, qk, l, u, x, lambda);
1229 }
1230 }
1231 /******************************************************************************/
1232 /*** End of gpm.cpp file ***/
1233 /******************************************************************************/
SHOGUN Machine Learning Toolbox - Documentation