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