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