59 #include <shogun/lib/external/pr_loqo.h>
63 #define HISTORY_BUF 1000000
65 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
67 CQPBSVMLib::CQPBSVMLib()
81 m_solver = QPB_SOLVER_SCA;
84 CQPBSVMLib::CQPBSVMLib(
99 m_solver = QPB_SOLVER_SCA;
102 CQPBSVMLib::~CQPBSVMLib()
107 int32_t CQPBSVMLib::solve_qp(
float64_t* result, int32_t len)
112 for (int32_t i=0; i<m_dim; i++)
116 m_diag_H=SG_MALLOC(float64_t, m_dim);
118 for (int32_t i=0; i<m_dim; i++)
119 m_diag_H[i]=m_H[i*m_dim+i];
121 float64_t* History=NULL;
127 case QPB_SOLVER_GRADDESC:
128 status = qpbsvm_gradient_descent(result, Nabla, &t, &History, verb );
131 status = qpbsvm_gauss_seidel(result, Nabla, &t, &History, verb );
134 status = qpbsvm_sca(result, Nabla, &t, &History, verb );
136 case QPB_SOLVER_SCAS:
137 status = qpbsvm_scas(result, Nabla, &t, &History, verb );
139 case QPB_SOLVER_SCAMV:
140 status = qpbsvm_scamv(result, Nabla, &t, &History, verb );
142 case QPB_SOLVER_PRLOQO:
143 status = qpbsvm_prloqo(result, Nabla, &t, &History, verb );
146 case QPB_SOLVER_CPLEX:
147 status = qpbsvm_cplex(result, Nabla, &t, &History, verb );
149 SG_ERROR(
"cplex not enabled at compile time - unknow solver\n")
171 int32_t CQPBSVMLib::qpbsvm_sca(float64_t *x,
174 float64_t **ptr_History,
187 int32_t History_size;
200 History=SG_MALLOC(float64_t, History_size*2);
201 memset(History, 0,
sizeof(float64_t)*History_size*2);
207 for(i = 0; i < m_dim; i++ ) {
208 xHx += x[i]*(Nabla[i] - m_f[i]);
214 Q_D = -0.5*xHx - m_UB*xi_sum;
215 History[
INDEX(0,t,2)] = Q_P;
216 History[
INDEX(1,t,2)] = Q_D;
219 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",
220 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/
CMath::abs(Q_P));
224 while( exitflag == -1 )
228 for(i = 0; i < m_dim; i++ ) {
229 if( m_diag_H[i] > 0 ) {
235 delta_x = x[i] - x_old;
237 col_H = (float64_t*)get_col(i);
238 for(j = 0; j < m_dim; j++ ) {
239 Nabla[j] += col_H[j]*delta_x;
251 for(i = 0; i < m_dim; i++ ) {
252 xHx += x[i]*(Nabla[i] - m_f[i]);
256 if((x[i] > 0 && x[i] < m_UB &&
CMath::abs(Nabla[i]) > m_tolKKT) ||
257 (x[i] == 0 && Nabla[i] < -m_tolKKT) ||
258 (x[i] == m_UB && Nabla[i] > m_tolKKT)) KKTsatisf = 0;
262 Q_D = -0.5*xHx - m_UB*xi_sum;
265 if(t >= m_tmax) exitflag = 0;
266 else if(Q_P-Q_D <= m_tolabs) exitflag = 1;
267 else if(Q_P-Q_D <=
CMath::abs(Q_P)*m_tolrel) exitflag = 2;
268 else if(KKTsatisf == 1) exitflag = 3;
270 if( verb > 0 && (t % verb == 0 || t==1)) {
271 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",
272 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/
CMath::abs(Q_P));
276 if( t < History_size ) {
277 History[
INDEX(0,t,2)] = Q_P;
278 History[
INDEX(1,t,2)] = Q_D;
281 tmp_ptr=SG_MALLOC(float64_t, (History_size+
HISTORY_BUF)*2);
282 memset(tmp_ptr, 0,
sizeof(float64_t)*(History_size+
HISTORY_BUF)*2);
284 for( i = 0; i < History_size; i++ ) {
285 tmp_ptr[
INDEX(0,i,2)] = History[
INDEX(0,i,2)];
286 tmp_ptr[
INDEX(1,i,2)] = History[
INDEX(1,i,2)];
288 tmp_ptr[
INDEX(0,t,2)] = Q_P;
289 tmp_ptr[
INDEX(1,t,2)] = Q_D;
298 (*ptr_History) = History;
300 SG_PRINT(
"QP: %f QD: %f\n", Q_P, Q_D)
312 int32_t CQPBSVMLib::qpbsvm_scas(float64_t *x,
315 float64_t **ptr_History,
330 float64_t max_update;
331 float64_t curr_update;
332 int32_t History_size;
346 History=SG_MALLOC(float64_t, History_size*2);
347 memset(History, 0,
sizeof(float64_t)*History_size*2);
353 for(i = 0; i < m_dim; i++ ) {
354 xHx += x[i]*(Nabla[i] - m_f[i]);
360 Q_D = -0.5*xHx - m_UB*xi_sum;
361 History[
INDEX(0,t,2)] = Q_P;
362 History[
INDEX(1,t,2)] = Q_D;
365 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",
366 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/
CMath::abs(Q_P));
370 while( exitflag == -1 )
375 for(i = 0; i < m_dim; i++ ) {
376 if( m_diag_H[i] > 0 ) {
381 curr_update = -0.5*m_diag_H[i]*(x_new*x_new-x_old*x_old) -
382 (Nabla[i] - m_diag_H[i]*x_old)*(x_new - x_old);
384 if( curr_update > max_update ) {
386 max_update = curr_update;
396 delta_x = max_x - x_old;
398 col_H = (float64_t*)get_col(max_i);
399 for(j = 0; j < m_dim; j++ ) {
400 Nabla[j] += col_H[j]*delta_x;
409 for(i = 0; i < m_dim; i++ ) {
410 xHx += x[i]*(Nabla[i] - m_f[i]);
414 if((x[i] > 0 && x[i] < m_UB &&
CMath::abs(Nabla[i]) > m_tolKKT) ||
415 (x[i] == 0 && Nabla[i] < -m_tolKKT) ||
416 (x[i] == m_UB && Nabla[i] > m_tolKKT)) KKTsatisf = 0;
420 Q_D = -0.5*xHx - m_UB*xi_sum;
423 if(t >= m_tmax) exitflag = 0;
424 else if(Q_P-Q_D <= m_tolabs) exitflag = 1;
425 else if(Q_P-Q_D <=
CMath::abs(Q_P)*m_tolrel) exitflag = 2;
426 else if(KKTsatisf == 1) exitflag = 3;
428 if( verb > 0 && (t % verb == 0 || t==1)) {
429 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",
430 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/
CMath::abs(Q_P));
434 if( t < History_size ) {
435 History[
INDEX(0,t,2)] = Q_P;
436 History[
INDEX(1,t,2)] = Q_D;
439 tmp_ptr=SG_MALLOC(float64_t, (History_size+
HISTORY_BUF)*2);
440 memset(tmp_ptr, 0, (History_size+
HISTORY_BUF)*2*
sizeof(float64_t));
441 for( i = 0; i < History_size; i++ ) {
442 tmp_ptr[
INDEX(0,i,2)] = History[
INDEX(0,i,2)];
443 tmp_ptr[
INDEX(1,i,2)] = History[
INDEX(1,i,2)];
445 tmp_ptr[
INDEX(0,t,2)] = Q_P;
446 tmp_ptr[
INDEX(1,t,2)] = Q_D;
455 (*ptr_History) = History;
466 int32_t CQPBSVMLib::qpbsvm_scamv(float64_t *x,
469 float64_t **ptr_History,
489 while( exitflag == -1 && t <= m_tmax)
494 for(i = 0; i < m_dim; i++ )
498 if( max_viol < -Nabla[i]) { u = i; max_viol = -Nabla[i]; }
500 else if( x[i] > 0 && x[i] < m_UB )
504 else if( max_viol < Nabla[i]) { u = i; max_viol = Nabla[i]; }
509 if( max_viol <= m_tolKKT )
518 delta_x = x_new - x[u];
521 col_H = (float64_t*)get_col(u);
522 for(i = 0; i < m_dim; i++ ) {
523 Nabla[i] += col_H[i]*delta_x;
528 History=SG_MALLOC(float64_t, (t+1)*2);
529 memset(History, 0,
sizeof(float64_t)*(t+1)*2);
532 for(fval = 0, i = 0; i < m_dim; i++ ) {
533 fval += 0.5*x[i]*(Nabla[i]+m_f[i]);
536 History[
INDEX(0,t,2)] = fval;
537 History[
INDEX(1,t,2)] = 0;
540 (*ptr_History) = History;
553 int32_t CQPBSVMLib::qpbsvm_prloqo(float64_t *x,
556 float64_t **ptr_History,
559 float64_t* lb=SG_MALLOC(float64_t, m_dim);
560 float64_t* ub=SG_MALLOC(float64_t, m_dim);
561 float64_t* primal=SG_MALLOC(float64_t, 3*m_dim);
562 float64_t* dual=SG_MALLOC(float64_t, 1+2*m_dim);
563 float64_t* a=SG_MALLOC(float64_t, m_dim);
565 for (int32_t i=0; i<m_dim; i++)
575 int32_t result=pr_loqo(m_dim, 1, m_f, m_H, a, &b, lb, ub, primal, dual,
576 2, 5, 1, -0.95, 10,0);
589 int32_t CQPBSVMLib::qpbsvm_gauss_seidel(float64_t *x,
592 float64_t **ptr_History,
595 for (int32_t i=0; i<m_dim; i++)
598 for (int32_t t=0; t<200; t++)
600 for (int32_t i=0; i<m_dim; i++)
602 x[i]= (-m_f[i]-(
CMath::dot(x,&m_H[m_dim*i], m_dim) -
603 m_H[m_dim*i+i]*x[i]))/m_H[m_dim*i+i];
609 for (int32_t i=0; i<m_dim; i++)
611 if (x[i]==0.0 || x[i]==1.0)
614 SG_PRINT(
"atbound:%d of %d (%2.2f%%)\n", atbound, m_dim, ((float64_t) 100.0*atbound)/m_dim)
620 int32_t CQPBSVMLib::qpbsvm_gradient_descent(float64_t *x,
623 float64_t **ptr_History,
626 for (int32_t i=0; i<m_dim; i++)
629 for (int32_t t=0; t<2000; t++)
631 for (int32_t i=0; i<m_dim; i++)
633 x[i]-=0.001*(
CMath::dot(x,&m_H[m_dim*i], m_dim)+m_f[i]);
639 for (int32_t i=0; i<m_dim; i++)
641 if (x[i]==0.0 || x[i]==1.0)
644 SG_PRINT(
"atbound:%d of %d (%2.2f%%)\n", atbound, m_dim, ((float64_t) 100.0*atbound)/m_dim)
657 int32_t CQPBSVMLib::qpbsvm_cplex(float64_t *x,
660 float64_t **ptr_History,
663 float64_t* lb=SG_MALLOC(float64_t, m_dim);
664 float64_t* ub=SG_MALLOC(float64_t, m_dim);
666 for (int32_t i=0; i<m_dim; i++)
674 cplex.
setup_lp(m_f, NULL, 0, m_dim, NULL, lb, ub);
687 #endif //USE_GPL_SHOGUN
Class CCplex to encapsulate access to the commercial cplex general purpose optimizer.
bool setup_qp(float64_t *H, int32_t dim)
bool init(E_PROB_TYPE t, int32_t timeout=60)
init cplex with problem type t and retry timeout 60 seconds
static const float64_t INFTY
infinity
#define INDEX(ROW, COL, DIM)
bool optimize(float64_t *sol, float64_t *lambda=NULL)
void display_vector(const char *name="vector", const char *prefix="") const
Class SGObject is the base class of all shogun objects.
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
all of classes and functions are contained in the shogun namespace
bool setup_lp(float64_t *objective, float64_t *constraints_mat, int32_t rows, int32_t cols, float64_t *rhs, float64_t *lb, float64_t *ub)
#define SG_UNSTABLE(func,...)
static T clamp(T value, T lb, T ub)