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