58 #include <shogun/lib/external/pr_loqo.h>
62 #define HISTORY_BUF 1000000
64 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
111 for (int32_t i=0; i<
m_dim; i++)
117 for (int32_t i=0; i<
m_dim; i++)
126 case QPB_SOLVER_GRADDESC:
133 status =
qpbsvm_sca(result, Nabla, &t, &History, verb );
135 case QPB_SOLVER_SCAS:
136 status =
qpbsvm_scas(result, Nabla, &t, &History, verb );
138 case QPB_SOLVER_SCAMV:
139 status =
qpbsvm_scamv(result, Nabla, &t, &History, verb );
141 case QPB_SOLVER_PRLOQO:
145 case QPB_SOLVER_CPLEX:
146 status =
qpbsvm_cplex(result, Nabla, &t, &History, verb );
148 SG_ERROR(
"cplex not enabled at compile time - unknow solver\n")
186 int32_t History_size;
199 History=SG_MALLOC(
float64_t, History_size*2);
200 memset(History, 0,
sizeof(
float64_t)*History_size*2);
206 for(i = 0; i <
m_dim; i++ ) {
207 xHx += x[i]*(Nabla[i] -
m_f[i]);
213 Q_D = -0.5*xHx -
m_UB*xi_sum;
214 History[
INDEX(0,t,2)] = Q_P;
215 History[
INDEX(1,t,2)] = Q_D;
218 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",
219 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/
CMath::abs(Q_P));
223 while( exitflag == -1 )
227 for(i = 0; i <
m_dim; i++ ) {
234 delta_x = x[i] - x_old;
237 for(j = 0; j <
m_dim; j++ ) {
238 Nabla[j] += col_H[j]*delta_x;
250 for(i = 0; i <
m_dim; i++ ) {
251 xHx += x[i]*(Nabla[i] -
m_f[i]);
256 (x[i] == 0 && Nabla[i] < -
m_tolKKT) ||
261 Q_D = -0.5*xHx -
m_UB*xi_sum;
264 if(t >=
m_tmax) exitflag = 0;
265 else if(Q_P-Q_D <=
m_tolabs) exitflag = 1;
267 else if(KKTsatisf == 1) exitflag = 3;
269 if( verb > 0 && (t % verb == 0 || t==1)) {
270 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",
271 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/
CMath::abs(Q_P));
275 if( t < History_size ) {
276 History[
INDEX(0,t,2)] = Q_P;
277 History[
INDEX(1,t,2)] = Q_D;
283 for( i = 0; i < History_size; i++ ) {
284 tmp_ptr[
INDEX(0,i,2)] = History[
INDEX(0,i,2)];
285 tmp_ptr[
INDEX(1,i,2)] = History[
INDEX(1,i,2)];
287 tmp_ptr[
INDEX(0,t,2)] = Q_P;
288 tmp_ptr[
INDEX(1,t,2)] = Q_D;
297 (*ptr_History) = History;
299 SG_PRINT(
"QP: %f QD: %f\n", Q_P, Q_D)
331 int32_t History_size;
345 History=SG_MALLOC(
float64_t, History_size*2);
346 memset(History, 0,
sizeof(
float64_t)*History_size*2);
352 for(i = 0; i <
m_dim; i++ ) {
353 xHx += x[i]*(Nabla[i] -
m_f[i]);
359 Q_D = -0.5*xHx -
m_UB*xi_sum;
360 History[
INDEX(0,t,2)] = Q_P;
361 History[
INDEX(1,t,2)] = Q_D;
364 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",
365 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/
CMath::abs(Q_P));
369 while( exitflag == -1 )
374 for(i = 0; i <
m_dim; i++ ) {
380 curr_update = -0.5*
m_diag_H[i]*(x_new*x_new-x_old*x_old) -
381 (Nabla[i] -
m_diag_H[i]*x_old)*(x_new - x_old);
383 if( curr_update > max_update ) {
385 max_update = curr_update;
395 delta_x = max_x - x_old;
398 for(j = 0; j <
m_dim; j++ ) {
399 Nabla[j] += col_H[j]*delta_x;
408 for(i = 0; i <
m_dim; i++ ) {
409 xHx += x[i]*(Nabla[i] -
m_f[i]);
414 (x[i] == 0 && Nabla[i] < -
m_tolKKT) ||
419 Q_D = -0.5*xHx -
m_UB*xi_sum;
422 if(t >=
m_tmax) exitflag = 0;
423 else if(Q_P-Q_D <=
m_tolabs) exitflag = 1;
425 else if(KKTsatisf == 1) exitflag = 3;
427 if( verb > 0 && (t % verb == 0 || t==1)) {
428 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",
429 t, Q_P, Q_D, Q_P-Q_D,(Q_P-Q_D)/
CMath::abs(Q_P));
433 if( t < History_size ) {
434 History[
INDEX(0,t,2)] = Q_P;
435 History[
INDEX(1,t,2)] = Q_D;
440 for( i = 0; i < History_size; i++ ) {
441 tmp_ptr[
INDEX(0,i,2)] = History[
INDEX(0,i,2)];
442 tmp_ptr[
INDEX(1,i,2)] = History[
INDEX(1,i,2)];
444 tmp_ptr[
INDEX(0,t,2)] = Q_P;
445 tmp_ptr[
INDEX(1,t,2)] = Q_D;
454 (*ptr_History) = History;
488 while( exitflag == -1 && t <=
m_tmax)
493 for(i = 0; i <
m_dim; i++ )
497 if( max_viol < -Nabla[i]) { u = i; max_viol = -Nabla[i]; }
499 else if( x[i] > 0 && x[i] <
m_UB )
503 else if( max_viol < Nabla[i]) { u = i; max_viol = Nabla[i]; }
517 delta_x = x_new - x[u];
521 for(i = 0; i <
m_dim; i++ ) {
522 Nabla[i] += col_H[i]*delta_x;
528 memset(History, 0,
sizeof(
float64_t)*(t+1)*2);
531 for(fval = 0, i = 0; i <
m_dim; i++ ) {
532 fval += 0.5*x[i]*(Nabla[i]+
m_f[i]);
535 History[
INDEX(0,t,2)] = fval;
536 History[
INDEX(1,t,2)] = 0;
539 (*ptr_History) = History;
564 for (int32_t i=0; i<
m_dim; i++)
574 int32_t result=pr_loqo(m_dim, 1,
m_f,
m_H, a, &b, lb, ub, primal, dual,
575 2, 5, 1, -0.95, 10,0);
594 for (int32_t i=0; i<
m_dim; i++)
597 for (int32_t t=0; t<200; t++)
599 for (int32_t i=0; i<
m_dim; i++)
602 m_H[m_dim*i+i]*x[i]))/
m_H[m_dim*i+i];
608 for (int32_t i=0; i<
m_dim; i++)
610 if (x[i]==0.0 || x[i]==1.0)
613 SG_PRINT(
"atbound:%d of %d (%2.2f%%)\n", atbound, m_dim, ((
float64_t) 100.0*atbound)/m_dim)
625 for (int32_t i=0; i<
m_dim; i++)
628 for (int32_t t=0; t<2000; t++)
630 for (int32_t i=0; i<
m_dim; i++)
638 for (int32_t i=0; i<
m_dim; i++)
640 if (x[i]==0.0 || x[i]==1.0)
643 SG_PRINT(
"atbound:%d of %d (%2.2f%%)\n", atbound, m_dim, ((
float64_t) 100.0*atbound)/m_dim)
665 for (int32_t i=0; i<
m_dim; i++)
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
int32_t qpbsvm_scamv(float64_t *x, float64_t *Nabla, int32_t *ptr_t, float64_t **ptr_History, int32_t verb)
float64_t * get_col(int32_t col)
static const float64_t INFTY
infinity
bool optimize(float64_t *sol, float64_t *lambda=NULL)
int32_t qpbsvm_sca(float64_t *x, float64_t *Nabla, int32_t *ptr_t, float64_t **ptr_History, int32_t verb)
void display_vector(const char *name="vector", const char *prefix="") const
#define INDEX(ROW, COL, DIM)
int32_t qpbsvm_gauss_seidel(float64_t *x, float64_t *Nabla, int32_t *ptr_t, float64_t **ptr_History, int32_t verb)
Class SGObject is the base class of all shogun objects.
int32_t qpbsvm_prloqo(float64_t *x, float64_t *Nabla, int32_t *ptr_t, float64_t **ptr_History, int32_t verb)
int32_t qpbsvm_gradient_descent(float64_t *x, float64_t *Nabla, int32_t *ptr_t, float64_t **ptr_History, int32_t verb)
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
Compute dot product between v1 and v2 (blas optimized)
int32_t qpbsvm_scas(float64_t *x, float64_t *Nabla, int32_t *ptr_t, float64_t **ptr_History, int32_t verb)
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)
int32_t qpbsvm_cplex(float64_t *x, float64_t *Nabla, int32_t *ptr_t, float64_t **ptr_History, int32_t verb)
int32_t solve_qp(float64_t *result, int32_t len)
result has to be allocated & zeroed
#define SG_UNSTABLE(func,...)
static T clamp(T value, T lb, T ub)