SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
pr_loqo.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Purpose: solves quadratic programming problem for pattern recognition
8  * for support vectors
9  *
10  * Written (W) 1997-1998 Alex J. Smola
11  * Written (W) 1999-2009 Soeren Sonnenburg
12  * Written (W) 1999-2008 Gunnar Raetsch
13  * Copyright (C) 1997-2009 Fraunhofer Institute FIRST and Max-Planck-Society
14  */
15 
16 #include <shogun/lib/common.h>
17 #include <shogun/io/SGIO.h>
21 
22 #include <math.h>
23 #include <time.h>
24 #include <stdlib.h>
25 #include <stdio.h>
26 #ifdef HAVE_LAPACK
27 extern "C" {
28 #ifdef HAVE_ATLAS
29 #include <atlas/clapack.h>
30 #endif
31 }
32 #endif
33 
34 namespace shogun
35 {
36 
37 #define PREDICTOR 1
38 #define CORRECTOR 2
39 
40 /*****************************************************************
41  replace this by any other function that will exit gracefully
42  in a larger system
43  ***************************************************************/
44 
45 void nrerror(char error_text[])
46 {
47  SG_SDEBUG("terminating optimizer - %s\n", error_text);
48  // exit(1);
49 }
50 
51 /*****************************************************************
52  taken from numerical recipes and modified to accept pointers
53  moreover numerical recipes code seems to be buggy (at least the
54  ones on the web)
55 
56  cholesky solver and backsubstitution
57  leaves upper right triangle intact (rows first order)
58  ***************************************************************/
59 
60 #ifdef HAVE_LAPACK
61 bool choldc(float64_t* a, int32_t n, float64_t* p)
62 {
63  if (n<=0)
64  return false;
65 
66  float64_t* a2=SG_MALLOC(float64_t, n*n);
67 
68  for (int32_t i=0; i<n; i++)
69  {
70  for (int32_t j=0; j<n; j++)
71  a2[n*i+j]=a[n*i+j];
72  }
73 
74  /* int for calling external lib */
75  int result=clapack_dpotrf(CblasRowMajor, CblasUpper, (int) n, a2, (int) n);
76 
77  for (int32_t i=0; i<n; i++)
78  p[i]=a2[(n+1)*i];
79 
80  for (int32_t i=0; i<n; i++)
81  {
82  for (int32_t j=i+1; j<n; j++)
83  {
84  a[n*j+i]=a2[n*i+j];
85  }
86  }
87 
88  if (result>0)
89  SG_SDEBUG("Choldc failed, matrix not positive definite\n");
90 
91  SG_FREE(a2);
92 
93  return result==0;
94 }
95 #else
96 bool choldc(float64_t a[], int32_t n, float64_t p[])
97 {
98  void nrerror(char error_text[]);
99  int32_t i, j, k;
100  float64_t sum;
101 
102  for (i = 0; i < n; i++)
103  {
104  for (j = i; j < n; j++)
105  {
106  sum=a[n*i + j];
107 
108  for (k=i-1; k>=0; k--)
109  sum -= a[n*i + k]*a[n*j + k];
110 
111  if (i == j)
112  {
113  if (sum <= 0.0)
114  {
115  SG_SDEBUG("Choldc failed, matrix not positive definite");
116  sum = 0.0;
117  return false;
118  }
119 
120  p[i]=sqrt(sum);
121 
122  }
123  else
124  a[n*j + i] = sum/p[i];
125  }
126  }
127 
128  return true;
129 }
130 #endif
131 
132 void cholsb(
133  float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[])
134 {
135  int32_t i, k;
136  float64_t sum;
137 
138  for (i=0; i<n; i++) {
139  sum=b[i];
140  for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
141  x[i]=sum/p[i];
142  }
143 
144  for (i=n-1; i>=0; i--) {
145  sum=x[i];
146  for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
147  x[i]=sum/p[i];
148  }
149 }
150 
151 /*****************************************************************
152  sometimes we only need the forward or backward pass of the
153  backsubstitution, hence we provide these two routines separately
154  ***************************************************************/
155 
157  float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[])
158 {
159  int32_t i, k;
160  float64_t sum;
161 
162  for (i=0; i<n; i++) {
163  sum=b[i];
164  for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
165  x[i]=sum/p[i];
166  }
167 }
168 
170  float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[])
171 {
172  int32_t i, k;
173  float64_t sum;
174 
175  for (i=n-1; i>=0; i--) {
176  sum=b[i];
177  for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
178  x[i]=sum/p[i];
179  }
180 }
181 
182 /*****************************************************************
183  solves the system | -H_x A' | |x_x| = |c_x|
184  | A H_y| |x_y| |c_y|
185 
186  with H_x (and H_y) positive (semidefinite) matrices
187  and n, m the respective sizes of H_x and H_y
188 
189  for variables see pg. 48 of notebook or do the calculations on a
190  sheet of paper again
191 
192  predictor solves the whole thing, corrector assues that H_x didn't
193  change and relies on the results of the predictor. therefore do
194  _not_ modify workspace
195 
196  if you want to speed tune anything in the code here's the right
197  place to do so: about 95% of the time is being spent in
198  here. something like an iterative refinement would be nice,
199  especially when switching from float64_t to single precision. if you
200  have a fast parallel cholesky use it instead of the numrec
201  implementations.
202 
203  side effects: changes H_y (but this is just the unit matrix or zero anyway
204  in our case)
205  ***************************************************************/
206 
208  int32_t n, int32_t m, float64_t h_x[], float64_t h_y[], float64_t a[],
209  float64_t x_x[], float64_t x_y[], float64_t c_x[], float64_t c_y[],
210  float64_t workspace[], int32_t step)
211 {
212  int32_t i,j,k;
213 
214  float64_t *p_x;
215  float64_t *p_y;
216  float64_t *t_a;
217  float64_t *t_c;
218  float64_t *t_y;
219 
220  p_x = workspace; /* together n + m + n*m + n + m = n*(m+2)+2*m */
221  p_y = p_x + n;
222  t_a = p_y + m;
223  t_c = t_a + n*m;
224  t_y = t_c + n;
225 
226  if (step == PREDICTOR) {
227  if (!choldc(h_x, n, p_x)) /* do cholesky decomposition */
228  return false;
229 
230  for (i=0; i<m; i++) /* forward pass for A' */
231  chol_forward(h_x, n, p_x, a+i*n, t_a+i*n);
232 
233  for (i=0; i<m; i++) /* compute (h_y + a h_x^-1A') */
234  for (j=i; j<m; j++)
235  for (k=0; k<n; k++)
236  h_y[m*i + j] += t_a[n*j + k] * t_a[n*i + k];
237 
238  choldc(h_y, m, p_y); /* and cholesky decomposition */
239  }
240 
241  chol_forward(h_x, n, p_x, c_x, t_c);
242  /* forward pass for c */
243 
244  for (i=0; i<m; i++) { /* and solve for x_y */
245  t_y[i] = c_y[i];
246  for (j=0; j<n; j++)
247  t_y[i] += t_a[i*n + j] * t_c[j];
248  }
249 
250  cholsb(h_y, m, p_y, t_y, x_y);
251 
252  for (i=0; i<n; i++) { /* finally solve for x_x */
253  t_c[i] = -t_c[i];
254  for (j=0; j<m; j++)
255  t_c[i] += t_a[j*n + i] * x_y[j];
256  }
257 
258  chol_backward(h_x, n, p_x, t_c, x_x);
259  return true;
260 }
261 
262 /*****************************************************************
263  matrix vector multiplication (symmetric matrix but only one triangle
264  given). computes m*x = y
265  no need to tune it as it's only of O(n^2) but cholesky is of
266  O(n^3). so don't waste your time _here_ although it isn't very
267  elegant.
268  ***************************************************************/
269 
270 void matrix_vector(int32_t n, float64_t m[], float64_t x[], float64_t y[])
271 {
272  int32_t i, j;
273 
274  for (i=0; i<n; i++) {
275  y[i] = m[(n+1) * i] * x[i];
276 
277  for (j=0; j<i; j++)
278  y[i] += m[i + n*j] * x[j];
279 
280  for (j=i+1; j<n; j++)
281  y[i] += m[n*i + j] * x[j];
282  }
283 }
284 
285 /*****************************************************************
286  call only this routine; this is the only one you're interested in
287  for doing quadratical optimization
288 
289  the restart feature exists but it may not be of much use due to the
290  fact that an initial setting, although close but not very close the
291  the actual solution will result in very good starting diagnostics
292  (primal and dual feasibility and small infeasibility gap) but incur
293  later stalling of the optimizer afterwards as we have to enforce
294  positivity of the slacks.
295  ***************************************************************/
296 
297 int32_t pr_loqo(
298  int32_t n, int32_t m, float64_t c[], float64_t h_x[], float64_t a[],
299  float64_t b[], float64_t l[], float64_t u[], float64_t primal[],
300  float64_t dual[], int32_t verb, float64_t sigfig_max, int32_t counter_max,
301  float64_t margin, float64_t bound, int32_t restart)
302 {
303  /* the knobs to be tuned ... */
304  /* float64_t margin = -0.95; we will go up to 95% of the
305  distance between old variables and zero */
306  /* float64_t bound = 10; preset value for the start. small
307  values give good initial
308  feasibility but may result in slow
309  convergence afterwards: we're too
310  close to zero */
311  /* to be allocated */
312  float64_t *workspace;
313  float64_t *diag_h_x;
314  float64_t *h_y;
315  float64_t *c_x;
316  float64_t *c_y;
317  float64_t *h_dot_x;
318  float64_t *rho;
319  float64_t *nu;
320  float64_t *tau;
321  float64_t *sigma;
322  float64_t *gamma_z;
323  float64_t *gamma_s;
324 
325  float64_t *hat_nu;
326  float64_t *hat_tau;
327 
328  float64_t *delta_x;
329  float64_t *delta_y;
330  float64_t *delta_s;
331  float64_t *delta_z;
332  float64_t *delta_g;
333  float64_t *delta_t;
334 
335  float64_t *d;
336 
337  /* from the header - pointers into primal and dual */
338  float64_t *x;
339  float64_t *y;
340  float64_t *g;
341  float64_t *z;
342  float64_t *s;
343  float64_t *t;
344 
345  /* auxiliary variables */
346  float64_t b_plus_1;
347  float64_t c_plus_1;
348 
349  float64_t x_h_x;
350  float64_t primal_inf;
351  float64_t dual_inf;
352 
353  float64_t sigfig;
354  float64_t primal_obj, dual_obj;
355  float64_t mu;
356  float64_t alfa=-1;
357  int32_t counter = 0;
358 
359  int32_t status = STILL_RUNNING;
360 
361  int32_t i,j;
362 
363  /* memory allocation */
364  workspace = SG_MALLOC(float64_t, (n*(m+2)+2*m));
365  diag_h_x = SG_MALLOC(float64_t, n);
366  h_y = SG_MALLOC(float64_t, m*m);
367  c_x = SG_MALLOC(float64_t, n);
368  c_y = SG_MALLOC(float64_t, m);
369  h_dot_x = SG_MALLOC(float64_t, n);
370 
371  rho = SG_MALLOC(float64_t, m);
372  nu = SG_MALLOC(float64_t, n);
373  tau = SG_MALLOC(float64_t, n);
374  sigma = SG_MALLOC(float64_t, n);
375 
376  gamma_z = SG_MALLOC(float64_t, n);
377  gamma_s = SG_MALLOC(float64_t, n);
378 
379  hat_nu = SG_MALLOC(float64_t, n);
380  hat_tau = SG_MALLOC(float64_t, n);
381 
382  delta_x = SG_MALLOC(float64_t, n);
383  delta_y = SG_MALLOC(float64_t, m);
384  delta_s = SG_MALLOC(float64_t, n);
385  delta_z = SG_MALLOC(float64_t, n);
386  delta_g = SG_MALLOC(float64_t, n);
387  delta_t = SG_MALLOC(float64_t, n);
388 
389  d = SG_MALLOC(float64_t, n);
390 
391  /* pointers into the external variables */
392  x = primal; /* n */
393  g = x + n; /* n */
394  t = g + n; /* n */
395 
396  y = dual; /* m */
397  z = y + m; /* n */
398  s = z + n; /* n */
399 
400  /* initial settings */
401  b_plus_1 = 1;
402  c_plus_1 = 0;
403  for (i=0; i<n; i++) c_plus_1 += c[i];
404 
405  /* get diagonal terms */
406  for (i=0; i<n; i++) diag_h_x[i] = h_x[(n+1)*i];
407 
408  /* starting point */
409  if (restart == 1) {
410  /* x, y already preset */
411  for (i=0; i<n; i++) { /* compute g, t for primal feasibility */
412  g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound);
413  t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound);
414  }
415 
416  matrix_vector(n, h_x, x, h_dot_x); /* h_dot_x = h_x * x */
417 
418  for (i=0; i<n; i++) { /* sigma is a dummy variable to calculate z, s */
419  sigma[i] = c[i] + h_dot_x[i];
420  for (j=0; j<m; j++)
421  sigma[i] -= a[n*j + i] * y[j];
422 
423  if (sigma[i] > 0) {
424  s[i] = bound;
425  z[i] = sigma[i] + bound;
426  }
427  else {
428  s[i] = bound - sigma[i];
429  z[i] = bound;
430  }
431  }
432  }
433  else { /* use default start settings */
434  for (i=0; i<m; i++)
435  for (j=i; j<m; j++)
436  h_y[i*m + j] = (i==j) ? 1 : 0;
437 
438  for (i=0; i<n; i++) {
439  c_x[i] = c[i];
440  h_x[(n+1)*i] += 1;
441  }
442 
443  for (i=0; i<m; i++)
444  c_y[i] = b[i];
445 
446  /* and solve the system [-H_x A'; A H_y] [x, y] = [c_x; c_y] */
447  solve_reduced(n, m, h_x, h_y, a, x, y, c_x, c_y, workspace,
448  PREDICTOR);
449 
450  /* initialize the other variables */
451  for (i=0; i<n; i++) {
452  g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound);
453  z[i] = CMath::max(CMath::abs(x[i]), bound);
454  t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound);
455  s[i] = CMath::max(CMath::abs(x[i]), bound);
456  }
457  }
458 
459  for (i=0, mu=0; i<n; i++)
460  mu += z[i] * g[i] + s[i] * t[i];
461  mu = mu / (2*n);
462 
463  /* the main loop */
464  if (verb >= STATUS) {
465  SG_SDEBUG("counter | pri_inf | dual_inf | pri_obj | dual_obj | ");
466  SG_SDEBUG("sigfig | alpha | nu \n");
467  SG_SDEBUG("-------------------------------------------------------");
468  SG_SDEBUG("---------------------------\n");
469  }
470 
471  while (status == STILL_RUNNING) {
472  /* predictor */
473 
474  /* put back original diagonal values */
475  for (i=0; i<n; i++)
476  h_x[(n+1) * i] = diag_h_x[i];
477 
478  matrix_vector(n, h_x, x, h_dot_x); /* compute h_dot_x = h_x * x */
479 
480  for (i=0; i<m; i++) {
481  rho[i] = b[i];
482  for (j=0; j<n; j++)
483  rho[i] -= a[n*i + j] * x[j];
484  }
485 
486  for (i=0; i<n; i++) {
487  nu[i] = l[i] - x[i] + g[i];
488  tau[i] = u[i] - x[i] - t[i];
489 
490  sigma[i] = c[i] - z[i] + s[i] + h_dot_x[i];
491  for (j=0; j<m; j++)
492  sigma[i] -= a[n*j + i] * y[j];
493 
494  gamma_z[i] = - z[i];
495  gamma_s[i] = - s[i];
496  }
497 
498  /* instrumentation */
499  x_h_x = 0;
500  primal_inf = 0;
501  dual_inf = 0;
502 
503  for (i=0; i<n; i++) {
504  x_h_x += h_dot_x[i] * x[i];
505  primal_inf += CMath::sq(tau[i]);
506  primal_inf += CMath::sq(nu[i]);
507  dual_inf += CMath::sq(sigma[i]);
508  }
509  for (i=0; i<m; i++)
510  primal_inf += CMath::sq(rho[i]);
511  primal_inf = sqrt(primal_inf)/b_plus_1;
512  dual_inf = sqrt(dual_inf)/c_plus_1;
513 
514  primal_obj = 0.5 * x_h_x;
515  dual_obj = -0.5 * x_h_x;
516  for (i=0; i<n; i++) {
517  primal_obj += c[i] * x[i];
518  dual_obj += l[i] * z[i] - u[i] * s[i];
519  }
520  for (i=0; i<m; i++)
521  dual_obj += b[i] * y[i];
522 
523  sigfig = log10(CMath::abs(primal_obj) + 1) -
524  log10(CMath::abs(primal_obj - dual_obj));
525  sigfig = CMath::max(sigfig, 0.0);
526 
527  /* the diagnostics - after we computed our results we will
528  analyze them */
529 
530  if (counter > counter_max) status = ITERATION_LIMIT;
531  if (sigfig > sigfig_max) status = OPTIMAL_SOLUTION;
532  if (primal_inf > 10e100) status = PRIMAL_INFEASIBLE;
533  if (dual_inf > 10e100) status = DUAL_INFEASIBLE;
534  if ((primal_inf > 10e100) & (dual_inf > 10e100)) status = PRIMAL_AND_DUAL_INFEASIBLE;
535  if (CMath::abs(primal_obj) > 10e100) status = PRIMAL_UNBOUNDED;
536  if (CMath::abs(dual_obj) > 10e100) status = DUAL_UNBOUNDED;
537 
538  /* write some nice routine to enforce the time limit if you
539  _really_ want, however it's quite useless as you can compute
540  the time from the maximum number of iterations as every
541  iteration costs one cholesky decomposition plus a couple of
542  backsubstitutions */
543 
544  /* generate report */
545  if ((verb >= FLOOD) | ((verb == STATUS) & (status != 0)))
546  SG_SDEBUG("%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n",
547  counter, primal_inf, dual_inf, primal_obj, dual_obj,
548  sigfig, alfa, mu);
549 
550  counter++;
551 
552  if (status == 0) { /* we may keep on going, otherwise
553  it'll cost one loop extra plus a
554  messed up main diagonal of h_x */
555  /* intermediate variables (the ones with hat) */
556  for (i=0; i<n; i++) {
557  hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
558  hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
559  /* diagonal terms */
560  d[i] = z[i] / g[i] + s[i] / t[i];
561  }
562 
563  /* initialization before the cholesky solver */
564  for (i=0; i<n; i++) {
565  h_x[(n+1)*i] = diag_h_x[i] + d[i];
566  c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] -
567  s[i] * hat_tau[i] / t[i];
568  }
569  for (i=0; i<m; i++) {
570  c_y[i] = rho[i];
571  for (j=i; j<m; j++)
572  h_y[m*i + j] = 0;
573  }
574 
575  /* and do it */
576  if (!solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace, PREDICTOR))
577  {
578  status=INCONSISTENT;
579  goto exit_optimizer;
580  }
581 
582  for (i=0; i<n; i++) {
583  /* backsubstitution */
584  delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
585  delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
586 
587  delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
588  delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
589 
590  /* central path (corrector) */
591  gamma_z[i] = mu / g[i] - z[i] - delta_z[i] * delta_g[i] / g[i];
592  gamma_s[i] = mu / t[i] - s[i] - delta_s[i] * delta_t[i] / t[i];
593 
594  /* (some more intermediate variables) the hat variables */
595  hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
596  hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
597 
598  /* initialization before the cholesky */
599  c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - s[i] * hat_tau[i] / t[i];
600  }
601 
602  for (i=0; i<m; i++) { /* comput c_y and rho */
603  c_y[i] = rho[i];
604  for (j=i; j<m; j++)
605  h_y[m*i + j] = 0;
606  }
607 
608  /* and do it */
609  solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,
610  CORRECTOR);
611 
612  for (i=0; i<n; i++) {
613  /* backsubstitution */
614  delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
615  delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
616 
617  delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
618  delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
619  }
620 
621  alfa = -1;
622  for (i=0; i<n; i++) {
623  alfa = CMath::min(alfa, delta_g[i]/g[i]);
624  alfa = CMath::min(alfa, delta_t[i]/t[i]);
625  alfa = CMath::min(alfa, delta_s[i]/s[i]);
626  alfa = CMath::min(alfa, delta_z[i]/z[i]);
627  }
628  alfa = (margin - 1) / alfa;
629 
630  /* compute mu */
631  for (i=0, mu=0; i<n; i++)
632  mu += z[i] * g[i] + s[i] * t[i];
633  mu = mu / (2*n);
634  mu = mu * CMath::sq((alfa - 1) / (alfa + 10));
635 
636  for (i=0; i<n; i++) {
637  x[i] += alfa * delta_x[i];
638  g[i] += alfa * delta_g[i];
639  t[i] += alfa * delta_t[i];
640  z[i] += alfa * delta_z[i];
641  s[i] += alfa * delta_s[i];
642  }
643 
644  for (i=0; i<m; i++)
645  y[i] += alfa * delta_y[i];
646  }
647  }
648 
649 exit_optimizer:
650  if ((status == 1) && (verb >= STATUS)) {
651  SG_SDEBUG("----------------------------------------------------------------------------------\n");
652  SG_SDEBUG("optimization converged\n");
653  }
654 
655  /* free memory */
656  SG_FREE(workspace);
657  SG_FREE(diag_h_x);
658  SG_FREE(h_y);
659  SG_FREE(c_x);
660  SG_FREE(c_y);
661  SG_FREE(h_dot_x);
662 
663  SG_FREE(rho);
664  SG_FREE(nu);
665  SG_FREE(tau);
666  SG_FREE(sigma);
667  SG_FREE(gamma_z);
668  SG_FREE(gamma_s);
669 
670  SG_FREE(hat_nu);
671  SG_FREE(hat_tau);
672 
673  SG_FREE(delta_x);
674  SG_FREE(delta_y);
675  SG_FREE(delta_s);
676  SG_FREE(delta_z);
677  SG_FREE(delta_g);
678  SG_FREE(delta_t);
679 
680  SG_FREE(d);
681 
682  /* and return to sender */
683  return status;
684 }
685 }

SHOGUN Machine Learning Toolbox - Documentation