19 using namespace shogun;
22 mkl_block_norm(1),beta_local(NULL), mkl_iterations(0), mkl_epsilon(1e-5),
23 interleaved_optimization(true), w_gap(1.0), rho(0)
36 SG_DEBUG(
"creating MKL object %p\n",
this)
45 SG_DEBUG(
"deleting MKL object %p\n",
this)
73 SG_INFO(
"trying to initialize CPLEX\n")
76 env = CPXopenCPLEX (&status);
81 SG_WARNING(
"Could not open CPLEX environment.\n")
82 CPXgeterrorstring (
env, status, errmsg);
90 status = CPXsetintparam (
env, CPX_PARAM_LPMETHOD, CPX_ALG_DUAL);
93 SG_ERROR(
"Failure to select dual lp optimization, error %d.\n", status)
97 status = CPXsetintparam (
env, CPX_PARAM_DATACHECK, CPX_ON);
100 SG_ERROR(
"Failure to turn on data checking, error %d.\n", status)
129 SG_WARNING(
"CPXfreeprob failed, error code %d.\n", status)
136 int32_t status = CPXcloseCPLEX (&
env);
142 SG_WARNING(
"Could not close CPLEX environment.\n")
143 CPXgeterrorstring (
env, status, errmsg);
157 glp_set_obj_dir(
lp_glpk, GLP_MIN);
164 glp_term_out(GLP_OFF);
180 int status = glp_get_status(lp);
182 if (status==GLP_INFEAS)
184 SG_PRINT(
"solution is infeasible!\n")
187 else if(status==GLP_NOFEAS)
189 SG_PRINT(
"problem has no feasible solution!\n")
205 SG_ERROR(
"%s::train_machine(): Number of training vectors (%d) does"
206 " not match number of labels (%d)\n",
get_name(),
214 SG_ERROR(
"No constraint generator (SVM) set\n")
231 int32_t num_weights = -1;
233 SG_INFO(
"num_kernels = %d\n", num_kernels)
236 ASSERT(num_weights==num_kernels)
299 SG_ERROR(
"Interleaved MKL optimization is currently "
300 "only supported with SVMlight\n");
313 #ifdef USE_REFERENCE_COUNTING
314 int32_t refs=this->ref();
320 #ifdef USE_REFERENCE_COUNTING
340 SG_SWARNING(
"MKL Algorithm terminates PREMATURELY due to current training time exceeding get_max_train_time ()= %f . It may have not converged yet!\n",
get_max_train_time ())
363 for (int32_t i=0; i<nsv; i++)
376 SG_ERROR(
"Norm must be >= 1, e.g., 1-norm is the standard MKL; norms>1 nonsparse MKL\n")
383 if (lambda>1 || lambda<0)
388 else if (lambda==1.0)
407 SG_SWARNING(
"MKL Algorithm terminates PREMATURELY due to current training time exceeding get_max_train_time ()= %f . It may have not converged yet!\n",
get_max_train_time ())
414 ASSERT(nweights==num_kernels)
417 #if defined(USE_CPLEX) || defined(USE_GLPK)
418 int32_t inner_iters=0;
423 for (int32_t i=0; i<num_kernels; i++)
426 mkl_objective+=old_beta[i]*sumw[i];
460 SG_ERROR(
"Solver type not supported (not compiled in?)\n")
481 int32_t nofKernelsGood;
484 nofKernelsGood = num_kernels;
487 for (p=0; p<num_kernels; ++p )
489 if (sumw[p] >= 0.0 && old_beta[p] >= 0.0 )
491 beta[p] =
CMath::sqrt(sumw[p]*old_beta[p]*old_beta[p]);
508 for( p=0; p<num_kernels; ++p )
513 SG_PRINT(
"MKL-direct: p = %.3f\n", 1.0 )
514 SG_PRINT(
"MKL-direct: nofKernelsGood = %d\n", nofKernelsGood )
515 SG_PRINT(
"MKL-direct: Z = %e\n", Z )
516 SG_PRINT(
"MKL-direct: eps = %e\n", epsRegul )
517 for( p=0; p<num_kernels; ++p )
522 SG_PRINT(
"MKL-direct: preR = %e\n", preR )
523 SG_PRINT(
"MKL-direct: preR/p = %e\n", preR )
525 SG_PRINT(
"MKL-direct: R = %e\n", R )
526 SG_ERROR(
"Assertion R >= 0 failed!\n" )
530 for( p=0; p<num_kernels; ++p )
538 for( p=0; p<num_kernels; ++p )
547 for( p=0; p<num_kernels; ++p )
555 for (p=0; p<num_kernels; ++p )
558 obj += sumw[p] * beta[p];
567 std::list<int32_t> I;
569 for (int32_t i=0; i<len;i++)
579 for (std::list<int32_t>::iterator it=I.begin(); it!=I.end(); it++)
608 for (int32_t i=0; i<n; i++)
612 for (int32_t j=0; j<n; j++)
648 }
while(ff>ff_old+1e-4*gg_old*(del-del_old));
660 SG_ERROR(
"cannot compute objective, labels or kernel not set\n")
675 for( p=0; p<num_kernels; ++p )
689 for( p=0; p<num_kernels; ++p )
694 for( p=0; p<num_kernels; ++p )
695 obj += sumw[p] * beta[p];
711 int32_t nofKernelsGood;
714 nofKernelsGood = num_kernels;
715 for( p=0; p<num_kernels; ++p )
718 if( sumw[p] >= 0.0 && old_beta[p] >= 0.0 )
720 beta[p] = sumw[p] * old_beta[p]*old_beta[p] /
mkl_norm;
733 for( p=0; p<num_kernels; ++p )
738 for( p=0; p<num_kernels; ++p )
743 for( p=0; p<num_kernels; ++p )
744 preR +=
CMath::sq( old_beta[p] - beta[p]);
750 SG_PRINT(
"MKL-direct: nofKernelsGood = %d\n", nofKernelsGood )
751 SG_PRINT(
"MKL-direct: Z = %e\n", Z )
752 SG_PRINT(
"MKL-direct: eps = %e\n", epsRegul )
753 for( p=0; p<num_kernels; ++p )
756 SG_PRINT(
"MKL-direct: t[%3d] = %e ( diff = %e = %e - %e )\n", p, t, old_beta[p]-beta[p], old_beta[p], beta[p] )
758 SG_PRINT(
"MKL-direct: preR = %e\n", preR )
761 SG_PRINT(
"MKL-direct: R = %e\n", R )
762 SG_ERROR(
"Assertion R >= 0 failed!\n" )
766 for( p=0; p<num_kernels; ++p )
774 for( p=0; p<num_kernels; ++p )
784 for( p=0; p<num_kernels; ++p )
785 obj += sumw[p] * beta[p];
791 const float64_t* old_beta, int32_t num_kernels,
798 SG_ERROR(
"MKL via NEWTON works only for norms>1\n")
800 const double epsBeta = 1e-32;
801 const double epsGamma = 1e-12;
802 const double epsWsq = 1e-12;
803 const double epsNewt = 0.0001;
804 const double epsStep = 1e-9;
805 const int nofNewtonSteps = 3;
806 const double hessRidge = 1e-6;
807 const int inLogSpace = 0;
822 for( p=0; p<num_kernels; ++p )
824 beta[p] = old_beta[p];
825 if( !( beta[p] >= epsBeta ) )
828 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 )
833 if( !( fabs(Z-1.0) <= epsGamma ) )
835 SG_WARNING(
"old_beta not normalized (diff=%e); forcing normalization. ", Z-1.0 )
836 for( p=0; p<num_kernels; ++p )
841 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 )
847 for ( p=0; p<num_kernels; ++p )
849 if ( !( sumw[p] >= 0 ) )
851 if( !( sumw[p] >= -epsWsq ) )
852 SG_WARNING(
"sumw[%d] = %e; treated as 0. ", p, sumw[p] )
864 if( !( gamma > epsGamma ) )
866 SG_WARNING(
"bad gamma: %e; set to %e. ", gamma, epsGamma )
870 ASSERT( gamma >= epsGamma )
875 for( p=0; p<num_kernels; ++p )
877 obj += beta[p] * sumw[p];
880 if( !( obj >= 0.0 ) )
886 for (i = 0; i < nofNewtonSteps; ++i )
891 for( p=0; p<num_kernels; ++p )
893 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 )
897 const float halfw2p = ( sumw[p] >= 0.0 ) ? (sumw[p]*old_beta[p]*old_beta[p]) : 0.0;
899 const float64_t t1 = ( t0 < 0 ) ? 0.0 : t0;
902 newtDir[p] = t1 / ( t1 + t2*beta[p] + hessRidge );
904 newtDir[p] = ( t1 == 0.0 ) ? 0.0 : ( t1 / t2 );
906 ASSERT( newtDir[p] == newtDir[p] )
914 while( stepSize >= epsStep )
920 for( p=0; p<num_kernels; ++p )
923 newtBeta[p] = beta[p] *
CMath::exp( + stepSize * newtDir[p] );
925 newtBeta[p] = beta[p] + stepSize * newtDir[p];
926 if( !( newtBeta[p] >= epsBeta ) )
927 newtBeta[p] = epsBeta;
939 for( p=0; p<num_kernels; ++p )
942 if( newtBeta[p] > 1.0 )
947 ASSERT( 0.0 <= newtBeta[p] && newtBeta[p] <= 1.0 )
953 for( p=0; p<num_kernels; ++p )
954 newtObj += sumw[p] * old_beta[p]*old_beta[p] / newtBeta[p];
956 if ( newtObj < obj - epsNewt*stepSize*obj )
958 for( p=0; p<num_kernels; ++p )
959 beta[p] = newtBeta[p];
967 if( stepSize < epsStep )
975 for( p=0; p<num_kernels; ++p )
976 obj += beta[p] * sumw[p];
991 int32_t NUMCOLS = 2*num_kernels + 1;
992 double* x=SG_MALLOC(
double, NUMCOLS);
1002 for (int32_t i=0; i<2*num_kernels; i++)
1009 for (int32_t i=num_kernels; i<2*num_kernels; i++)
1012 obj[2*num_kernels]=1 ;
1013 lb[2*num_kernels]=-CPX_INFBOUND ;
1014 ub[2*num_kernels]=CPX_INFBOUND ;
1016 int status = CPXnewcols (
env,
lp_cplex, NUMCOLS, obj, lb, ub, NULL, NULL);
1019 CPXgeterrorstring (
env, status, errmsg);
1024 SG_INFO(
"adding the first row\n")
1025 int initial_rmatbeg[1];
1026 int initial_rmatind[num_kernels+1];
1027 double initial_rmatval[num_kernels+1];
1028 double initial_rhs[1];
1029 char initial_sense[1];
1034 initial_rmatbeg[0] = 0;
1036 initial_sense[0]=
'E' ;
1039 for (int32_t i=0; i<num_kernels; i++)
1041 initial_rmatind[i]=i ;
1042 initial_rmatval[i]=1 ;
1044 initial_rmatind[num_kernels]=2*num_kernels ;
1045 initial_rmatval[num_kernels]=0 ;
1047 status = CPXaddrows (
env,
lp_cplex, 0, 1, num_kernels+1,
1048 initial_rhs, initial_sense, initial_rmatbeg,
1049 initial_rmatind, initial_rmatval, NULL, NULL);
1054 initial_rmatbeg[0] = 0;
1056 initial_sense[0]=
'L' ;
1058 initial_rmatind[0]=2*num_kernels ;
1059 initial_rmatval[0]=0 ;
1062 initial_rhs, initial_sense, initial_rmatbeg,
1063 initial_rmatind, initial_rmatval, NULL, NULL);
1068 for (int32_t i=0; i<num_kernels; i++)
1070 initial_rmatind[i]=i ;
1071 initial_rmatval[i]=1 ;
1073 initial_rmatind[num_kernels]=2*num_kernels ;
1074 initial_rmatval[num_kernels]=0 ;
1076 status = CPXaddqconstr (
env,
lp_cplex, 0, num_kernels+1, 1.0,
'L', NULL, NULL,
1077 initial_rmatind, initial_rmatind, initial_rmatval, NULL);
1083 SG_ERROR(
"Failed to add the first row.\n")
1089 for (int32_t q=0; q<num_kernels-1; q++)
1106 rmatind[2]=num_kernels+q ;
1109 rhs, sense, rmatbeg,
1110 rmatind, rmatval, NULL, NULL);
1112 SG_ERROR(
"Failed to add a smothness row (1).\n")
1121 rmatind[2]=num_kernels+q ;
1124 rhs, sense, rmatbeg,
1125 rmatind, rmatval, NULL, NULL);
1127 SG_ERROR(
"Failed to add a smothness row (2).\n")
1136 int rmatind[num_kernels+1];
1137 double rmatval[num_kernels+1];
1149 for (int32_t i=0; i<num_kernels; i++)
1153 rmatval[i]=-(sumw[i]-suma) ;
1155 rmatval[i]=-sumw[i];
1157 rmatind[num_kernels]=2*num_kernels ;
1158 rmatval[num_kernels]=-1 ;
1160 int32_t status = CPXaddrows (
env,
lp_cplex, 0, 1, num_kernels+1,
1161 rhs, sense, rmatbeg,
1162 rmatind, rmatval, NULL, NULL);
1164 SG_ERROR(
"Failed to add the new row.\n")
1179 for (int32_t i=0; i<num_kernels; i++)
1180 beta[i]=old_beta[i];
1181 for (int32_t i=num_kernels; i<2*num_kernels+1; i++)
1189 CMath::scale_vector(1/CMath::qnorm(beta, num_kernels,
mkl_norm), beta, num_kernels);
1195 SG_ERROR(
"Failed to optimize Problem.\n")
1199 status=CPXsolution(
env,
lp_cplex, &solstat, &objval,
1200 (
double*) beta, NULL, NULL, NULL);
1204 CMath::display_vector(beta, num_kernels,
"beta");
1205 SG_ERROR(
"Failed to obtain solution.\n")
1208 CMath::scale_vector(1/CMath::qnorm(beta, num_kernels,
mkl_norm), beta, num_kernels);
1222 SG_ERROR(
"Failed to optimize Problem.\n")
1225 int32_t cur_numrows=(int32_t) CPXgetnumrows(
env,
lp_cplex);
1226 int32_t cur_numcols=(int32_t) CPXgetnumcols(
env,
lp_cplex);
1227 int32_t num_rows=cur_numrows;
1228 ASSERT(cur_numcols<=2*num_kernels+1)
1239 status=CPXsolution(
env,
lp_cplex, &solstat, &objval,
1240 (
double*) x, (
double*) pi, (
double*) slack, NULL);
1244 status=CPXsolution(
env,
lp_cplex, &solstat, &objval,
1245 (
double*) x, NULL, (
double*) slack, NULL);
1248 int32_t solution_ok = (!status) ;
1250 SG_ERROR(
"Failed to obtain solution.\n")
1252 int32_t num_active_rows=0 ;
1257 int32_t max_idx = -1 ;
1258 int32_t start_row = 1 ;
1260 start_row+=2*(num_kernels-1);
1262 for (int32_t i = start_row; i < cur_numrows; i++)
1270 if (slack[i]>max_slack)
1272 max_slack=slack[i] ;
1283 if (slack[i]>max_slack)
1285 max_slack=slack[i] ;
1293 if ( (num_rows-start_row>
CMath::max(100,2*num_active_rows)) && (max_idx!=-1))
1296 status = CPXdelrows (
env,
lp_cplex, max_idx, max_idx) ;
1298 SG_ERROR(
"Failed to remove an old row.\n")
1303 rho = -x[2*num_kernels] ;
1315 for (int32_t i=0; i<num_kernels; i++)
1320 SG_ERROR(
"Cplex not enabled at compile time\n")
1331 SG_ERROR(
"MKL via GLPK works only for norm=1\n")
1335 int32_t NUMCOLS = 2*num_kernels + 1 ;
1341 glp_add_cols(
lp_glpk, NUMCOLS);
1342 for (
int i=1; i<=2*num_kernels; i++)
1344 glp_set_obj_coef(
lp_glpk, i, 0);
1345 glp_set_col_bnds(
lp_glpk, i, GLP_DB, 0, 1);
1347 for (
int i=num_kernels+1; i<=2*num_kernels; i++)
1351 glp_set_obj_coef(
lp_glpk, NUMCOLS, 1);
1352 glp_set_col_bnds(
lp_glpk, NUMCOLS, GLP_FR, 0,0);
1355 int row_index = glp_add_rows(
lp_glpk, 1);
1356 int* ind = SG_MALLOC(
int, num_kernels+2);
1358 for (
int i=1; i<=num_kernels; i++)
1363 ind[num_kernels+1] = NUMCOLS;
1364 val[num_kernels+1] = 0;
1365 glp_set_mat_row(
lp_glpk, row_index, num_kernels, ind, val);
1366 glp_set_row_bnds(
lp_glpk, row_index, GLP_FX, 1, 1);
1374 for (int32_t q=1; q<num_kernels; q++)
1378 int mat_row_index = glp_add_rows(
lp_glpk, 2);
1383 mat_ind[3] = num_kernels+q;
1385 glp_set_mat_row(
lp_glpk, mat_row_index, 3, mat_ind, mat_val);
1386 glp_set_row_bnds(
lp_glpk, mat_row_index, GLP_UP, 0, 0);
1389 glp_set_mat_row(
lp_glpk, mat_row_index+1, 3, mat_ind, mat_val);
1390 glp_set_row_bnds(
lp_glpk, mat_row_index+1, GLP_UP, 0, 0);
1395 int* ind=SG_MALLOC(
int,num_kernels+2);
1397 int row_index = glp_add_rows(
lp_glpk, 1);
1398 for (int32_t i=1; i<=num_kernels; i++)
1401 val[i] = -(sumw[i-1]-suma);
1403 ind[num_kernels+1] = 2*num_kernels+1;
1404 val[num_kernels+1] = -1;
1405 glp_set_mat_row(
lp_glpk, row_index, num_kernels+1, ind, val);
1406 glp_set_row_bnds(
lp_glpk, row_index, GLP_UP, 0, 0);
1414 SG_ERROR(
"Failed to optimize Problem.\n")
1416 int32_t cur_numrows = glp_get_num_rows(
lp_glpk);
1417 int32_t cur_numcols = glp_get_num_cols(
lp_glpk);
1418 int32_t num_rows=cur_numrows;
1419 ASSERT(cur_numcols<=2*num_kernels+1)
1425 for (
int i=0; i<cur_numrows; i++)
1427 row_primal[i] = glp_get_row_prim(
lp_glpk, i+1);
1428 row_dual[i] = glp_get_row_dual(
lp_glpk, i+1);
1430 for (
int i=0; i<cur_numcols; i++)
1431 col_primal[i] = glp_get_col_prim(
lp_glpk, i+1);
1433 obj = -col_primal[2*num_kernels];
1435 for (
int i=0; i<num_kernels; i++)
1436 beta[i] = col_primal[i];
1438 int32_t num_active_rows=0;
1442 int32_t max_idx = -1;
1443 int32_t start_row = 1;
1445 start_row += 2*(num_kernels-1);
1447 for (int32_t i= start_row; i<cur_numrows; i++)
1453 if (row_primal[i]<max_slack)
1455 max_slack = row_primal[i];
1461 if ((num_rows-start_row>
CMath::max(100, 2*num_active_rows)) && max_idx!=-1)
1464 del_rows[1] = max_idx+1;
1465 glp_del_rows(
lp_glpk, 1, del_rows);
1470 SG_FREE(row_primal);
1471 SG_FREE(col_primal);
1473 SG_ERROR(
"Glpk not enabled at compile time\n")
1489 ASSERT(nweights==num_kernels)
1492 for (int32_t i=0; i<num_kernels; i++)
1498 for (int32_t n=0; n<num_kernels; n++)
1505 for (int32_t i=0; i<nsv; i++)
1509 for (int32_t j=0; j<nsv; j++)
1541 for (int32_t i=0; i<n; i++)
1545 for (int32_t j=0; j<n; j++)
1561 mkl_obj=-0.5*mkl_obj;
1568 SG_ERROR(
"cannot compute objective, labels or kernel not set\n")
1581 int* ind=SG_MALLOC(
int, num_kernels+1);
1584 double const_term = 1-CMath::qsq(beta, num_kernels,
mkl_norm);
1587 ASSERT(CMath::fequal(const_term, 0.0))
1589 for (int32_t i=0; i<num_kernels; i++)
1593 lin_term[i]=grad_beta[i] - 2*beta[i]*hess_beta[i];
1594 const_term+=grad_beta[i]*beta[i] -
CMath::sq(beta[i])*hess_beta[i];
1597 ind[num_kernels]=2*num_kernels;
1598 hess_beta[num_kernels]=0;
1599 lin_term[num_kernels]=0;
1610 status = CPXaddqconstr (
env,
lp_cplex, num_kernels+1, num_kernels+1, const_term,
'L', ind, lin_term,
1611 ind, ind, hess_beta, NULL);