00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 #include <math.h>
00050 #include <string.h>
00051 #include <limits.h>
00052
00053 #include "lib/config.h"
00054 #include "lib/io.h"
00055 #include "lib/Cplex.h"
00056 #include "lib/Mathematics.h"
00057
00058 #include "classifier/svm/qpbsvmlib.h"
00059 #include "classifier/svm/pr_loqo.h"
00060
00061 using namespace shogun;
00062
00063 #define HISTORY_BUF 1000000
00064
00065 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
00066
00067 CQPBSVMLib::CQPBSVMLib(
00068 float64_t* H, int32_t n, float64_t* f, int32_t m, float64_t UB)
00069 : CSGObject()
00070 {
00071 ASSERT(H && n>0);
00072 m_H=H;
00073 m_dim = n;
00074 m_diag_H=NULL;
00075
00076 m_f=f;
00077 m_UB=UB;
00078 m_tmax = INT_MAX;
00079 m_tolabs = 0;
00080 m_tolrel = 1e-6;
00081 m_tolKKT = 0;
00082 m_solver = QPB_SOLVER_SCA;
00083 }
00084
00085 CQPBSVMLib::~CQPBSVMLib()
00086 {
00087 delete[] m_diag_H;
00088 }
00089
00090 int32_t CQPBSVMLib::solve_qp(float64_t* result, int32_t len)
00091 {
00092 int32_t status = -1;
00093 ASSERT(len==m_dim);
00094 float64_t* Nabla=new float64_t[m_dim];
00095 for (int32_t i=0; i<m_dim; i++)
00096 Nabla[i]=m_f[i];
00097
00098 delete[] m_diag_H;
00099 m_diag_H=new float64_t[m_dim];
00100
00101 for (int32_t i=0; i<m_dim; i++)
00102 m_diag_H[i]=m_H[i*m_dim+i];
00103
00104 float64_t* History=NULL;
00105 int32_t t;
00106 int32_t verb=0;
00107
00108 switch (m_solver)
00109 {
00110 case QPB_SOLVER_GRADDESC:
00111 status = qpbsvm_gradient_descent(result, Nabla, &t, &History, verb );
00112 break;
00113 case QPB_SOLVER_GS:
00114 status = qpbsvm_gauss_seidel(result, Nabla, &t, &History, verb );
00115 break;
00116 case QPB_SOLVER_SCA:
00117 status = qpbsvm_sca(result, Nabla, &t, &History, verb );
00118 break;
00119 case QPB_SOLVER_SCAS:
00120 status = qpbsvm_scas(result, Nabla, &t, &History, verb );
00121 break;
00122 case QPB_SOLVER_SCAMV:
00123 status = qpbsvm_scamv(result, Nabla, &t, &History, verb );
00124 break;
00125 case QPB_SOLVER_PRLOQO:
00126 status = qpbsvm_prloqo(result, Nabla, &t, &History, verb );
00127 break;
00128 #ifdef USE_CPLEX
00129 case QPB_SOLVER_CPLEX:
00130 status = qpbsvm_cplex(result, Nabla, &t, &History, verb );
00131 #else
00132 SG_ERROR("cplex not enabled at compile time - unknow solver\n");
00133 #endif
00134 break;
00135 default:
00136 SG_ERROR("unknown solver\n");
00137 break;
00138 }
00139
00140 delete[] History;
00141 delete[] Nabla;
00142 delete[] m_diag_H;
00143 m_diag_H=NULL;
00144
00145 return status;
00146 }
00147
00148
00149
00150
00151
00152
00153
00154 int32_t CQPBSVMLib::qpbsvm_sca(float64_t *x,
00155 float64_t *Nabla,
00156 int32_t *ptr_t,
00157 float64_t **ptr_History,
00158 int32_t verb)
00159 {
00160 float64_t *History;
00161 float64_t *col_H;
00162 float64_t *tmp_ptr;
00163 float64_t x_old;
00164 float64_t delta_x;
00165 float64_t xHx;
00166 float64_t Q_P;
00167 float64_t Q_D;
00168 float64_t xf;
00169 float64_t xi_sum;
00170 int32_t History_size;
00171 int32_t t;
00172 int32_t i, j;
00173 int32_t exitflag;
00174 int32_t KKTsatisf;
00175
00176
00177
00178
00179
00180 t = 0;
00181
00182 History_size = (m_tmax < HISTORY_BUF ) ? m_tmax+1 : HISTORY_BUF;
00183 History=new float64_t[History_size*2];
00184 memset(History, 0, sizeof(float64_t)*History_size*2);
00185
00186
00187 xHx = 0;
00188 xf = 0;
00189 xi_sum = 0;
00190 for(i = 0; i < m_dim; i++ ) {
00191 xHx += x[i]*(Nabla[i] - m_f[i]);
00192 xf += x[i]*m_f[i];
00193 xi_sum += CMath::max(0.0,-Nabla[i]);
00194 }
00195
00196 Q_P = 0.5*xHx + xf;
00197 Q_D = -0.5*xHx - m_UB*xi_sum;
00198 History[INDEX(0,t,2)] = Q_P;
00199 History[INDEX(1,t,2)] = Q_D;
00200
00201 if( verb > 0 ) {
00202 SG_PRINT("%d: Q_P=%m_f, Q_D=%m_f, Q_P-Q_D=%m_f, (Q_P-Q_D)/|Q_P|=%m_f \n",
00203 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/CMath::abs(Q_P));
00204 }
00205
00206 exitflag = -1;
00207 while( exitflag == -1 )
00208 {
00209 t++;
00210
00211 for(i = 0; i < m_dim; i++ ) {
00212 if( m_diag_H[i] > 0 ) {
00213
00214 x_old = x[i];
00215 x[i] = CMath::min(m_UB,CMath::max(0.0, x[i] - Nabla[i]/m_diag_H[i]));
00216
00217
00218 delta_x = x[i] - x_old;
00219 if( delta_x != 0 ) {
00220 col_H = (float64_t*)get_col(i);
00221 for(j = 0; j < m_dim; j++ ) {
00222 Nabla[j] += col_H[j]*delta_x;
00223 }
00224 }
00225
00226 }
00227 }
00228
00229
00230 xHx = 0;
00231 xf = 0;
00232 xi_sum = 0;
00233 KKTsatisf = 1;
00234 for(i = 0; i < m_dim; i++ ) {
00235 xHx += x[i]*(Nabla[i] - m_f[i]);
00236 xf += x[i]*m_f[i];
00237 xi_sum += CMath::max(0.0,-Nabla[i]);
00238
00239 if((x[i] > 0 && x[i] < m_UB && CMath::abs(Nabla[i]) > m_tolKKT) ||
00240 (x[i] == 0 && Nabla[i] < -m_tolKKT) ||
00241 (x[i] == m_UB && Nabla[i] > m_tolKKT)) KKTsatisf = 0;
00242 }
00243
00244 Q_P = 0.5*xHx + xf;
00245 Q_D = -0.5*xHx - m_UB*xi_sum;
00246
00247
00248 if(t >= m_tmax) exitflag = 0;
00249 else if(Q_P-Q_D <= m_tolabs) exitflag = 1;
00250 else if(Q_P-Q_D <= CMath::abs(Q_P)*m_tolrel) exitflag = 2;
00251 else if(KKTsatisf == 1) exitflag = 3;
00252
00253 if( verb > 0 && (t % verb == 0 || t==1)) {
00254 SG_PRINT("%d: Q_P=%m_f, Q_D=%m_f, Q_P-Q_D=%m_f, (Q_P-Q_D)/|Q_P|=%m_f \n",
00255 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/CMath::abs(Q_P));
00256 }
00257
00258
00259 if( t < History_size ) {
00260 History[INDEX(0,t,2)] = Q_P;
00261 History[INDEX(1,t,2)] = Q_D;
00262 }
00263 else {
00264 tmp_ptr=new float64_t[(History_size+HISTORY_BUF)*2];
00265 memset(tmp_ptr, 0, sizeof(float64_t)*(History_size+HISTORY_BUF)*2);
00266
00267 for( i = 0; i < History_size; i++ ) {
00268 tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
00269 tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
00270 }
00271 tmp_ptr[INDEX(0,t,2)] = Q_P;
00272 tmp_ptr[INDEX(1,t,2)] = Q_D;
00273
00274 History_size += HISTORY_BUF;
00275 delete[] History;
00276 History = tmp_ptr;
00277 }
00278 }
00279
00280 (*ptr_t) = t;
00281 (*ptr_History) = History;
00282
00283 SG_PRINT("QP: %f QD: %f\n", Q_P, Q_D);
00284
00285 return( exitflag );
00286 }
00287
00288
00289
00290
00291
00292
00293
00294
00295 int32_t CQPBSVMLib::qpbsvm_scas(float64_t *x,
00296 float64_t *Nabla,
00297 int32_t *ptr_t,
00298 float64_t **ptr_History,
00299 int32_t verb)
00300 {
00301 float64_t *History;
00302 float64_t *col_H;
00303 float64_t *tmp_ptr;
00304 float64_t x_old;
00305 float64_t x_new;
00306 float64_t delta_x;
00307 float64_t max_x=CMath::INFTY;
00308 float64_t xHx;
00309 float64_t Q_P;
00310 float64_t Q_D;
00311 float64_t xf;
00312 float64_t xi_sum;
00313 float64_t max_update;
00314 float64_t curr_update;
00315 int32_t History_size;
00316 int32_t t;
00317 int32_t i, j;
00318 int32_t max_i=-1;
00319 int32_t exitflag;
00320 int32_t KKTsatisf;
00321
00322
00323
00324
00325
00326 t = 0;
00327
00328 History_size = (m_tmax < HISTORY_BUF ) ? m_tmax+1 : HISTORY_BUF;
00329 History=new float64_t[History_size*2];
00330 memset(History, 0, sizeof(float64_t)*History_size*2);
00331
00332
00333 xHx = 0;
00334 xf = 0;
00335 xi_sum = 0;
00336 for(i = 0; i < m_dim; i++ ) {
00337 xHx += x[i]*(Nabla[i] - m_f[i]);
00338 xf += x[i]*m_f[i];
00339 xi_sum += CMath::max(0.0,-Nabla[i]);
00340 }
00341
00342 Q_P = 0.5*xHx + xf;
00343 Q_D = -0.5*xHx - m_UB*xi_sum;
00344 History[INDEX(0,t,2)] = Q_P;
00345 History[INDEX(1,t,2)] = Q_D;
00346
00347 if( verb > 0 ) {
00348 SG_PRINT("%d: Q_P=%m_f, Q_D=%m_f, Q_P-Q_D=%m_f, (Q_P-Q_D)/|Q_P|=%m_f \n",
00349 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/CMath::abs(Q_P));
00350 }
00351
00352 exitflag = -1;
00353 while( exitflag == -1 )
00354 {
00355 t++;
00356
00357 max_update = -CMath::INFTY;
00358 for(i = 0; i < m_dim; i++ ) {
00359 if( m_diag_H[i] > 0 ) {
00360
00361 x_old = x[i];
00362 x_new = CMath::min(m_UB,CMath::max(0.0, x[i] - Nabla[i]/m_diag_H[i]));
00363
00364 curr_update = -0.5*m_diag_H[i]*(x_new*x_new-x_old*x_old) -
00365 (Nabla[i] - m_diag_H[i]*x_old)*(x_new - x_old);
00366
00367 if( curr_update > max_update ) {
00368 max_i = i;
00369 max_update = curr_update;
00370 max_x = x_new;
00371 }
00372 }
00373 }
00374
00375 x_old = x[max_i];
00376 x[max_i] = max_x;
00377
00378
00379 delta_x = max_x - x_old;
00380 if( delta_x != 0 ) {
00381 col_H = (float64_t*)get_col(max_i);
00382 for(j = 0; j < m_dim; j++ ) {
00383 Nabla[j] += col_H[j]*delta_x;
00384 }
00385 }
00386
00387
00388 xHx = 0;
00389 xf = 0;
00390 xi_sum = 0;
00391 KKTsatisf = 1;
00392 for(i = 0; i < m_dim; i++ ) {
00393 xHx += x[i]*(Nabla[i] - m_f[i]);
00394 xf += x[i]*m_f[i];
00395 xi_sum += CMath::max(0.0,-Nabla[i]);
00396
00397 if((x[i] > 0 && x[i] < m_UB && CMath::abs(Nabla[i]) > m_tolKKT) ||
00398 (x[i] == 0 && Nabla[i] < -m_tolKKT) ||
00399 (x[i] == m_UB && Nabla[i] > m_tolKKT)) KKTsatisf = 0;
00400 }
00401
00402 Q_P = 0.5*xHx + xf;
00403 Q_D = -0.5*xHx - m_UB*xi_sum;
00404
00405
00406 if(t >= m_tmax) exitflag = 0;
00407 else if(Q_P-Q_D <= m_tolabs) exitflag = 1;
00408 else if(Q_P-Q_D <= CMath::abs(Q_P)*m_tolrel) exitflag = 2;
00409 else if(KKTsatisf == 1) exitflag = 3;
00410
00411 if( verb > 0 && (t % verb == 0 || t==1)) {
00412 SG_PRINT("%d: Q_P=%m_f, Q_D=%m_f, Q_P-Q_D=%m_f, (Q_P-Q_D)/|Q_P|=%m_f \n",
00413 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/CMath::abs(Q_P));
00414 }
00415
00416
00417 if( t < History_size ) {
00418 History[INDEX(0,t,2)] = Q_P;
00419 History[INDEX(1,t,2)] = Q_D;
00420 }
00421 else {
00422 tmp_ptr=new float64_t[(History_size+HISTORY_BUF)*2];
00423 memset(tmp_ptr, 0, (History_size+HISTORY_BUF)*2*sizeof(float64_t));
00424 for( i = 0; i < History_size; i++ ) {
00425 tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
00426 tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
00427 }
00428 tmp_ptr[INDEX(0,t,2)] = Q_P;
00429 tmp_ptr[INDEX(1,t,2)] = Q_D;
00430
00431 History_size += HISTORY_BUF;
00432 delete[] History;
00433 History = tmp_ptr;
00434 }
00435 }
00436
00437 (*ptr_t) = t;
00438 (*ptr_History) = History;
00439
00440 return( exitflag );
00441 }
00442
00443
00444
00445
00446
00447
00448
00449 int32_t CQPBSVMLib::qpbsvm_scamv(float64_t *x,
00450 float64_t *Nabla,
00451 int32_t *ptr_t,
00452 float64_t **ptr_History,
00453 int32_t verb)
00454 {
00455 float64_t *History;
00456 float64_t *col_H;
00457 float64_t delta_x;
00458 float64_t x_new;
00459 float64_t max_viol;
00460 float64_t fval;
00461 int32_t t;
00462 int32_t i;
00463 int32_t u=-1;
00464 int32_t exitflag;
00465
00466
00467
00468
00469
00470 t = 0;
00471 exitflag = -1;
00472 while( exitflag == -1 && t <= m_tmax)
00473 {
00474 t++;
00475
00476 max_viol = 0;
00477 for(i = 0; i < m_dim; i++ )
00478 {
00479 if( x[i] == 0 )
00480 {
00481 if( max_viol < -Nabla[i]) { u = i; max_viol = -Nabla[i]; }
00482 }
00483 else if( x[i] > 0 && x[i] < m_UB )
00484 {
00485 if( max_viol < CMath::abs(Nabla[i]) ) { u = i; max_viol = CMath::abs(Nabla[i]); }
00486 }
00487 else if( max_viol < Nabla[i]) { u = i; max_viol = Nabla[i]; }
00488 }
00489
00490
00491
00492 if( max_viol <= m_tolKKT )
00493 {
00494 exitflag = 1;
00495 }
00496 else
00497 {
00498
00499 x_new = CMath::min(m_UB,CMath::max(0.0, x[u] - Nabla[u]/m_diag_H[u]));
00500
00501 delta_x = x_new - x[u];
00502 x[u] = x_new;
00503
00504 col_H = (float64_t*)get_col(u);
00505 for(i = 0; i < m_dim; i++ ) {
00506 Nabla[i] += col_H[i]*delta_x;
00507 }
00508 }
00509 }
00510
00511 History=new float64_t[(t+1)*2];
00512 memset(History, 0, sizeof(float64_t)*(t+1)*2);
00513
00514 fval = 0;
00515 for(fval = 0, i = 0; i < m_dim; i++ ) {
00516 fval += 0.5*x[i]*(Nabla[i]+m_f[i]);
00517 }
00518
00519 History[INDEX(0,t,2)] = fval;
00520 History[INDEX(1,t,2)] = 0;
00521
00522 (*ptr_t) = t;
00523 (*ptr_History) = History;
00524
00525
00526
00527 return( exitflag );
00528 }
00529
00530
00531
00532
00533
00534
00535
00536 int32_t CQPBSVMLib::qpbsvm_prloqo(float64_t *x,
00537 float64_t *Nabla,
00538 int32_t *ptr_t,
00539 float64_t **ptr_History,
00540 int32_t verb)
00541 {
00542 float64_t* lb=new float64_t[m_dim];
00543 float64_t* ub=new float64_t[m_dim];
00544 float64_t* primal=new float64_t[3*m_dim];
00545 float64_t* dual=new float64_t[1+2*m_dim];
00546 float64_t* a=new float64_t[m_dim];
00547
00548 for (int32_t i=0; i<m_dim; i++)
00549 {
00550 a[i]=0.0;
00551 lb[i]=0;
00552 ub[i]=m_UB;
00553 }
00554
00555 float64_t b=0;
00556
00557 CMath::display_vector(m_f, m_dim, "m_f");
00558 int32_t result=pr_loqo(m_dim, 1, m_f, m_H, a, &b, lb, ub, primal, dual,
00559 2, 5, 1, -0.95, 10,0);
00560
00561 delete[] a;
00562 delete[] lb;
00563 delete[] ub;
00564 delete[] primal;
00565 delete[] dual;
00566
00567 *ptr_t=0;
00568 *ptr_History=NULL;
00569 return result;
00570 }
00571
00572 int32_t CQPBSVMLib::qpbsvm_gauss_seidel(float64_t *x,
00573 float64_t *Nabla,
00574 int32_t *ptr_t,
00575 float64_t **ptr_History,
00576 int32_t verb)
00577 {
00578 for (int32_t i=0; i<m_dim; i++)
00579 x[i]=CMath::random(0.0, 1.0);
00580
00581 for (int32_t t=0; t<200; t++)
00582 {
00583 for (int32_t i=0; i<m_dim; i++)
00584 {
00585 x[i]= (-m_f[i]-(CMath::dot(x,&m_H[m_dim*i], m_dim) -
00586 m_H[m_dim*i+i]*x[i]))/m_H[m_dim*i+i];
00587 x[i]=CMath::clamp(x[i], 0.0, 1.0);
00588 }
00589 }
00590
00591 int32_t atbound=0;
00592 for (int32_t i=0; i<m_dim; i++)
00593 {
00594 if (x[i]==0.0 || x[i]==1.0)
00595 atbound++;
00596 }
00597 SG_PRINT("atbound:%d of %d (%2.2f%%)\n", atbound, m_dim, ((float64_t) 100.0*atbound)/m_dim);
00598 *ptr_t=0;
00599 *ptr_History=NULL;
00600 return 0;
00601 }
00602
00603 int32_t CQPBSVMLib::qpbsvm_gradient_descent(float64_t *x,
00604 float64_t *Nabla,
00605 int32_t *ptr_t,
00606 float64_t **ptr_History,
00607 int32_t verb)
00608 {
00609 for (int32_t i=0; i<m_dim; i++)
00610 x[i]=CMath::random(0.0, 1.0);
00611
00612 for (int32_t t=0; t<2000; t++)
00613 {
00614 for (int32_t i=0; i<m_dim; i++)
00615 {
00616 x[i]-=0.001*(CMath::dot(x,&m_H[m_dim*i], m_dim)+m_f[i]);
00617 x[i]=CMath::clamp(x[i], 0.0, 1.0);
00618 }
00619 }
00620
00621 int32_t atbound=0;
00622 for (int32_t i=0; i<m_dim; i++)
00623 {
00624 if (x[i]==0.0 || x[i]==1.0)
00625 atbound++;
00626 }
00627 SG_PRINT("atbound:%d of %d (%2.2f%%)\n", atbound, m_dim, ((float64_t) 100.0*atbound)/m_dim);
00628 *ptr_t=0;
00629 *ptr_History=NULL;
00630 return 0;
00631 }
00632
00633 #ifdef USE_CPLEX
00634
00635
00636
00637
00638
00639
00640 int32_t CQPBSVMLib::qpbsvm_cplex(float64_t *x,
00641 float64_t *Nabla,
00642 int32_t *ptr_t,
00643 float64_t **ptr_History,
00644 int32_t verb)
00645 {
00646 float64_t* lb=new float64_t[m_dim];
00647 float64_t* ub=new float64_t[m_dim];
00648
00649 for (int32_t i=0; i<m_dim; i++)
00650 {
00651 lb[i]=0;
00652 ub[i]=m_UB;
00653 }
00654
00655 CCplex cplex;
00656 cplex.init(E_QP);
00657 cplex.setup_lp(m_f, NULL, 0, m_dim, NULL, lb, ub);
00658 cplex.setup_qp(m_H, m_dim);
00659 cplex.optimize(x);
00660 cplex.cleanup();
00661
00662 delete[] lb;
00663 delete[] ub;
00664
00665 *ptr_t=0;
00666 *ptr_History=NULL;
00667 return 0;
00668 }
00669 #endif