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