00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "lib/common.h"
00017 #include "lib/io.h"
00018 #include "lib/lapack.h"
00019 #include "lib/Mathematics.h"
00020 #include "classifier/svm/pr_loqo.h"
00021
00022 #include <math.h>
00023 #include <time.h>
00024 #include <stdlib.h>
00025 #include <stdio.h>
00026 #ifdef HAVE_LAPACK
00027 extern "C" {
00028 #ifdef HAVE_ATLAS
00029 #include <clapack.h>
00030 #endif
00031 }
00032 #endif
00033
00034 #define PREDICTOR 1
00035 #define CORRECTOR 2
00036
00037
00038
00039
00040
00041
00042 void nrerror(CHAR error_text[])
00043 {
00044 SG_SDEBUG("terminating optimizer - %s\n", error_text);
00045
00046 }
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 #ifdef HAVE_LAPACK
00058 bool choldc(double* a, int n, double* p)
00059 {
00060 if (n<=0)
00061 return false;
00062
00063 double* a2=new double[n*n];
00064
00065 for (int i=0; i<n; i++)
00066 {
00067 for (int j=0; j<n; j++)
00068 a2[n*i+j]=a[n*i+j];
00069 }
00070
00071 int result=clapack_dpotrf(CblasRowMajor, CblasUpper, n, a2, n);
00072
00073 for (int i=0; i<n; i++)
00074 p[i]=a2[(n+1)*i];
00075
00076 for (int i=0; i<n; i++)
00077 {
00078 for (int j=i+1; j<n; j++)
00079 {
00080 a[n*j+i]=a2[n*i+j];
00081 }
00082 }
00083
00084 if (result>0)
00085 SG_SDEBUG("Choldc failed, matrix not positive definite\n");
00086
00087 delete[] a2;
00088
00089 return result==0;
00090 }
00091 #else
00092 bool choldc(double a[], int n, double p[])
00093 {
00094 void nrerror(CHAR error_text[]);
00095 int i, j, k;
00096 double sum;
00097
00098 for (i = 0; i < n; i++)
00099 {
00100 for (j = i; j < n; j++)
00101 {
00102 sum=a[n*i + j];
00103
00104 for (k=i-1; k>=0; k--)
00105 sum -= a[n*i + k]*a[n*j + k];
00106
00107 if (i == j)
00108 {
00109 if (sum <= 0.0)
00110 {
00111 SG_SERROR("Choldc failed, matrix not positive definite");
00112 sum = 0.0;
00113 return false;
00114 }
00115
00116 p[i]=sqrt(sum);
00117
00118 }
00119 else
00120 a[n*j + i] = sum/p[i];
00121 }
00122 }
00123
00124 return true;
00125 }
00126 #endif
00127
00128 void cholsb(double a[], int n, double p[], double b[], double x[])
00129 {
00130 int i, k;
00131 double sum;
00132
00133 for (i=0; i<n; i++) {
00134 sum=b[i];
00135 for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
00136 x[i]=sum/p[i];
00137 }
00138
00139 for (i=n-1; i>=0; i--) {
00140 sum=x[i];
00141 for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
00142 x[i]=sum/p[i];
00143 }
00144 }
00145
00146
00147
00148
00149
00150
00151 void chol_forward(double a[], int n, double p[], double b[], double x[])
00152 {
00153 int i, k;
00154 double sum;
00155
00156 for (i=0; i<n; i++) {
00157 sum=b[i];
00158 for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
00159 x[i]=sum/p[i];
00160 }
00161 }
00162
00163 void chol_backward(double a[], int n, double p[], double b[], double x[])
00164 {
00165 int i, k;
00166 double sum;
00167
00168 for (i=n-1; i>=0; i--) {
00169 sum=b[i];
00170 for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
00171 x[i]=sum/p[i];
00172 }
00173 }
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200 bool solve_reduced(int n, int m, double h_x[], double h_y[],
00201 double a[], double x_x[], double x_y[],
00202 double c_x[], double c_y[],
00203 double workspace[], int step)
00204 {
00205 int i,j,k;
00206
00207 double *p_x;
00208 double *p_y;
00209 double *t_a;
00210 double *t_c;
00211 double *t_y;
00212
00213 p_x = workspace;
00214 p_y = p_x + n;
00215 t_a = p_y + m;
00216 t_c = t_a + n*m;
00217 t_y = t_c + n;
00218
00219 if (step == PREDICTOR) {
00220 if (!choldc(h_x, n, p_x))
00221 return false;
00222
00223 for (i=0; i<m; i++)
00224 chol_forward(h_x, n, p_x, a+i*n, t_a+i*n);
00225
00226 for (i=0; i<m; i++)
00227 for (j=i; j<m; j++)
00228 for (k=0; k<n; k++)
00229 h_y[m*i + j] += t_a[n*j + k] * t_a[n*i + k];
00230
00231 choldc(h_y, m, p_y);
00232 }
00233
00234 chol_forward(h_x, n, p_x, c_x, t_c);
00235
00236
00237 for (i=0; i<m; i++) {
00238 t_y[i] = c_y[i];
00239 for (j=0; j<n; j++)
00240 t_y[i] += t_a[i*n + j] * t_c[j];
00241 }
00242
00243 cholsb(h_y, m, p_y, t_y, x_y);
00244
00245 for (i=0; i<n; i++) {
00246 t_c[i] = -t_c[i];
00247 for (j=0; j<m; j++)
00248 t_c[i] += t_a[j*n + i] * x_y[j];
00249 }
00250
00251 chol_backward(h_x, n, p_x, t_c, x_x);
00252 return true;
00253 }
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263 void matrix_vector(int n, double m[], double x[], double y[])
00264 {
00265 int i, j;
00266
00267 for (i=0; i<n; i++) {
00268 y[i] = m[(n+1) * i] * x[i];
00269
00270 for (j=0; j<i; j++)
00271 y[i] += m[i + n*j] * x[j];
00272
00273 for (j=i+1; j<n; j++)
00274 y[i] += m[n*i + j] * x[j];
00275 }
00276 }
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290 int pr_loqo(int n, int m, double c[], double h_x[], double a[], double b[],
00291 double l[], double u[], double primal[], double dual[],
00292 int verb, double sigfig_max, int counter_max,
00293 double margin, double bound, int restart)
00294 {
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304 double *workspace;
00305 double *diag_h_x;
00306 double *h_y;
00307 double *c_x;
00308 double *c_y;
00309 double *h_dot_x;
00310 double *rho;
00311 double *nu;
00312 double *tau;
00313 double *sigma;
00314 double *gamma_z;
00315 double *gamma_s;
00316
00317 double *hat_nu;
00318 double *hat_tau;
00319
00320 double *delta_x;
00321 double *delta_y;
00322 double *delta_s;
00323 double *delta_z;
00324 double *delta_g;
00325 double *delta_t;
00326
00327 double *d;
00328
00329
00330 double *x;
00331 double *y;
00332 double *g;
00333 double *z;
00334 double *s;
00335 double *t;
00336
00337
00338 double b_plus_1;
00339 double c_plus_1;
00340
00341 double x_h_x;
00342 double primal_inf;
00343 double dual_inf;
00344
00345 double sigfig;
00346 double primal_obj, dual_obj;
00347 double mu;
00348 double alfa=-1;
00349 int counter = 0;
00350
00351 int status = STILL_RUNNING;
00352
00353 int i,j;
00354
00355
00356 workspace = (double*) malloc((n*(m+2)+2*m)*sizeof(double));
00357 diag_h_x = (double*) malloc(n*sizeof(double));
00358 h_y = (double*) malloc(m*m*sizeof(double));
00359 c_x = (double*) malloc(n*sizeof(double));
00360 c_y = (double*) malloc(m*sizeof(double));
00361 h_dot_x = (double*) malloc(n*sizeof(double));
00362
00363 rho = (double*) malloc(m*sizeof(double));
00364 nu = (double*) malloc(n*sizeof(double));
00365 tau = (double*) malloc(n*sizeof(double));
00366 sigma = (double*) malloc(n*sizeof(double));
00367
00368 gamma_z = (double*) malloc(n*sizeof(double));
00369 gamma_s = (double*) malloc(n*sizeof(double));
00370
00371 hat_nu = (double*) malloc(n*sizeof(double));
00372 hat_tau = (double*) malloc(n*sizeof(double));
00373
00374 delta_x = (double*) malloc(n*sizeof(double));
00375 delta_y = (double*) malloc(m*sizeof(double));
00376 delta_s = (double*) malloc(n*sizeof(double));
00377 delta_z = (double*) malloc(n*sizeof(double));
00378 delta_g = (double*) malloc(n*sizeof(double));
00379 delta_t = (double*) malloc(n*sizeof(double));
00380
00381 d = (double*) malloc(n*sizeof(double));
00382
00383
00384 x = primal;
00385 g = x + n;
00386 t = g + n;
00387
00388 y = dual;
00389 z = y + m;
00390 s = z + n;
00391
00392
00393 b_plus_1 = 1;
00394 c_plus_1 = 0;
00395 for (i=0; i<n; i++) c_plus_1 += c[i];
00396
00397
00398 for (i=0; i<n; i++) diag_h_x[i] = h_x[(n+1)*i];
00399
00400
00401 if (restart == 1) {
00402
00403 for (i=0; i<n; i++) {
00404 g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound);
00405 t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound);
00406 }
00407
00408 matrix_vector(n, h_x, x, h_dot_x);
00409
00410 for (i=0; i<n; i++) {
00411 sigma[i] = c[i] + h_dot_x[i];
00412 for (j=0; j<m; j++)
00413 sigma[i] -= a[n*j + i] * y[j];
00414
00415 if (sigma[i] > 0) {
00416 s[i] = bound;
00417 z[i] = sigma[i] + bound;
00418 }
00419 else {
00420 s[i] = bound - sigma[i];
00421 z[i] = bound;
00422 }
00423 }
00424 }
00425 else {
00426 for (i=0; i<m; i++)
00427 for (j=i; j<m; j++)
00428 h_y[i*m + j] = (i==j) ? 1 : 0;
00429
00430 for (i=0; i<n; i++) {
00431 c_x[i] = c[i];
00432 h_x[(n+1)*i] += 1;
00433 }
00434
00435 for (i=0; i<m; i++)
00436 c_y[i] = b[i];
00437
00438
00439 solve_reduced(n, m, h_x, h_y, a, x, y, c_x, c_y, workspace,
00440 PREDICTOR);
00441
00442
00443 for (i=0; i<n; i++) {
00444 g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound);
00445 z[i] = CMath::max(CMath::abs(x[i]), bound);
00446 t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound);
00447 s[i] = CMath::max(CMath::abs(x[i]), bound);
00448 }
00449 }
00450
00451 for (i=0, mu=0; i<n; i++)
00452 mu += z[i] * g[i] + s[i] * t[i];
00453 mu = mu / (2*n);
00454
00455
00456 if (verb >= STATUS) {
00457 SG_SDEBUG("counter | pri_inf | dual_inf | pri_obj | dual_obj | ");
00458 SG_SDEBUG("sigfig | alpha | nu \n");
00459 SG_SDEBUG("-------------------------------------------------------");
00460 SG_SDEBUG("---------------------------\n");
00461 }
00462
00463 while (status == STILL_RUNNING) {
00464
00465
00466
00467 for (i=0; i<n; i++)
00468 h_x[(n+1) * i] = diag_h_x[i];
00469
00470 matrix_vector(n, h_x, x, h_dot_x);
00471
00472 for (i=0; i<m; i++) {
00473 rho[i] = b[i];
00474 for (j=0; j<n; j++)
00475 rho[i] -= a[n*i + j] * x[j];
00476 }
00477
00478 for (i=0; i<n; i++) {
00479 nu[i] = l[i] - x[i] + g[i];
00480 tau[i] = u[i] - x[i] - t[i];
00481
00482 sigma[i] = c[i] - z[i] + s[i] + h_dot_x[i];
00483 for (j=0; j<m; j++)
00484 sigma[i] -= a[n*j + i] * y[j];
00485
00486 gamma_z[i] = - z[i];
00487 gamma_s[i] = - s[i];
00488 }
00489
00490
00491 x_h_x = 0;
00492 primal_inf = 0;
00493 dual_inf = 0;
00494
00495 for (i=0; i<n; i++) {
00496 x_h_x += h_dot_x[i] * x[i];
00497 primal_inf += CMath::sq(tau[i]);
00498 primal_inf += CMath::sq(nu[i]);
00499 dual_inf += CMath::sq(sigma[i]);
00500 }
00501 for (i=0; i<m; i++)
00502 primal_inf += CMath::sq(rho[i]);
00503 primal_inf = sqrt(primal_inf)/b_plus_1;
00504 dual_inf = sqrt(dual_inf)/c_plus_1;
00505
00506 primal_obj = 0.5 * x_h_x;
00507 dual_obj = -0.5 * x_h_x;
00508 for (i=0; i<n; i++) {
00509 primal_obj += c[i] * x[i];
00510 dual_obj += l[i] * z[i] - u[i] * s[i];
00511 }
00512 for (i=0; i<m; i++)
00513 dual_obj += b[i] * y[i];
00514
00515 sigfig = log10(CMath::abs(primal_obj) + 1) -
00516 log10(CMath::abs(primal_obj - dual_obj));
00517 sigfig = CMath::max(sigfig, 0.0);
00518
00519
00520
00521
00522 if (counter > counter_max) status = ITERATION_LIMIT;
00523 if (sigfig > sigfig_max) status = OPTIMAL_SOLUTION;
00524 if (primal_inf > 10e100) status = PRIMAL_INFEASIBLE;
00525 if (dual_inf > 10e100) status = DUAL_INFEASIBLE;
00526 if ((primal_inf > 10e100) & (dual_inf > 10e100)) status = PRIMAL_AND_DUAL_INFEASIBLE;
00527 if (CMath::abs(primal_obj) > 10e100) status = PRIMAL_UNBOUNDED;
00528 if (CMath::abs(dual_obj) > 10e100) status = DUAL_UNBOUNDED;
00529
00530
00531
00532
00533
00534
00535
00536
00537 if ((verb >= FLOOD) | ((verb == STATUS) & (status != 0)))
00538 SG_SDEBUG("%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n",
00539 counter, primal_inf, dual_inf, primal_obj, dual_obj,
00540 sigfig, alfa, mu);
00541
00542 counter++;
00543
00544 if (status == 0) {
00545
00546
00547
00548 for (i=0; i<n; i++) {
00549 hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
00550 hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
00551
00552 d[i] = z[i] / g[i] + s[i] / t[i];
00553 }
00554
00555
00556 for (i=0; i<n; i++) {
00557 h_x[(n+1)*i] = diag_h_x[i] + d[i];
00558 c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] -
00559 s[i] * hat_tau[i] / t[i];
00560 }
00561 for (i=0; i<m; i++) {
00562 c_y[i] = rho[i];
00563 for (j=i; j<m; j++)
00564 h_y[m*i + j] = 0;
00565 }
00566
00567
00568 if (!solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace, PREDICTOR))
00569 {
00570 status=INCONSISTENT;
00571 goto exit_optimizer;
00572 }
00573
00574 for (i=0; i<n; i++) {
00575
00576 delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
00577 delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
00578
00579 delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
00580 delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
00581
00582
00583 gamma_z[i] = mu / g[i] - z[i] - delta_z[i] * delta_g[i] / g[i];
00584 gamma_s[i] = mu / t[i] - s[i] - delta_s[i] * delta_t[i] / t[i];
00585
00586
00587 hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
00588 hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
00589
00590
00591 c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - s[i] * hat_tau[i] / t[i];
00592 }
00593
00594 for (i=0; i<m; i++) {
00595 c_y[i] = rho[i];
00596 for (j=i; j<m; j++)
00597 h_y[m*i + j] = 0;
00598 }
00599
00600
00601 solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,
00602 CORRECTOR);
00603
00604 for (i=0; i<n; i++) {
00605
00606 delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
00607 delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
00608
00609 delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
00610 delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
00611 }
00612
00613 alfa = -1;
00614 for (i=0; i<n; i++) {
00615 alfa = CMath::min(alfa, delta_g[i]/g[i]);
00616 alfa = CMath::min(alfa, delta_t[i]/t[i]);
00617 alfa = CMath::min(alfa, delta_s[i]/s[i]);
00618 alfa = CMath::min(alfa, delta_z[i]/z[i]);
00619 }
00620 alfa = (margin - 1) / alfa;
00621
00622
00623 for (i=0, mu=0; i<n; i++)
00624 mu += z[i] * g[i] + s[i] * t[i];
00625 mu = mu / (2*n);
00626 mu = mu * CMath::sq((alfa - 1) / (alfa + 10));
00627
00628 for (i=0; i<n; i++) {
00629 x[i] += alfa * delta_x[i];
00630 g[i] += alfa * delta_g[i];
00631 t[i] += alfa * delta_t[i];
00632 z[i] += alfa * delta_z[i];
00633 s[i] += alfa * delta_s[i];
00634 }
00635
00636 for (i=0; i<m; i++)
00637 y[i] += alfa * delta_y[i];
00638 }
00639 }
00640
00641 exit_optimizer:
00642 if ((status == 1) && (verb >= STATUS)) {
00643 SG_SDEBUG("----------------------------------------------------------------------------------\n");
00644 SG_SDEBUG("optimization converged\n");
00645 }
00646
00647
00648 free(workspace);
00649 free(diag_h_x);
00650 free(h_y);
00651 free(c_x);
00652 free(c_y);
00653 free(h_dot_x);
00654
00655 free(rho);
00656 free(nu);
00657 free(tau);
00658 free(sigma);
00659 free(gamma_z);
00660 free(gamma_s);
00661
00662 free(hat_nu);
00663 free(hat_tau);
00664
00665 free(delta_x);
00666 free(delta_y);
00667 free(delta_s);
00668 free(delta_z);
00669 free(delta_g);
00670 free(delta_t);
00671
00672 free(d);
00673
00674
00675 return status;
00676 }