19 using namespace shogun;
22 :
CSVM(), svm(NULL), C_mkl(0), mkl_norm(1), ent_lambda(0), beta_local(NULL),
23 mkl_iterations(0), mkl_epsilon(1e-5), interleaved_optimization(true),
36 SG_DEBUG(
"creating MKL object %p\n",
this);
46 SG_DEBUG(
"deleting MKL object %p\n",
this);
74 SG_INFO(
"trying to initialize CPLEX\n") ;
77 env = CPXopenCPLEX (&status);
82 SG_WARNING(
"Could not open CPLEX environment.\n");
83 CPXgeterrorstring (
env, status, errmsg);
91 status = CPXsetintparam (
env, CPX_PARAM_LPMETHOD, CPX_ALG_DUAL);
94 SG_ERROR(
"Failure to select dual lp optimization, error %d.\n", status);
98 status = CPXsetintparam (
env, CPX_PARAM_DATACHECK, CPX_ON);
101 SG_ERROR(
"Failure to turn on data checking, error %d.\n", status);
108 SG_ERROR(
"Failed to create LP.\n");
130 SG_WARNING(
"CPXfreeprob failed, error code %d.\n", status);
137 int32_t status = CPXcloseCPLEX (&
env);
143 SG_WARNING(
"Could not close CPLEX environment.\n");
144 CPXgeterrorstring (
env, status, errmsg);
158 lpx_set_obj_dir(
lp_glpk, LPX_MIN);
159 lpx_set_int_parm(
lp_glpk, LPX_K_DUAL, GLP_ON );
160 lpx_set_int_parm(
lp_glpk, LPX_K_PRESOL, GLP_ON );
162 glp_term_out(GLP_OFF);
177 int status = lpx_get_status(lp);
179 if (status==LPX_INFEAS)
181 SG_PRINT(
"solution is infeasible!\n");
184 else if(status==LPX_NOFEAS)
186 SG_PRINT(
"problem has no feasible solution!\n");
201 SG_ERROR(
"Number of training vectors does not match number of labels\n");
207 SG_ERROR(
"No constraint generator (SVM) set\n");
224 int32_t num_weights = -1;
226 SG_INFO(
"num_kernels = %d\n", num_kernels);
229 ASSERT(num_weights==num_kernels);
292 SG_ERROR(
"Interleaved MKL optimization is currently "
293 "only supported with SVMlight\n");
306 #ifdef USE_REFERENCE_COUNTING
307 int32_t refs=this->ref();
313 #ifdef USE_REFERENCE_COUNTING
333 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 ());
356 for (int32_t i=0; i<nsv; i++)
369 SG_ERROR(
"Norm must be >= 1, e.g., 1-norm is the standard MKL; norms>1 nonsparse MKL\n");
376 if (lambda>1 || lambda<0)
381 else if (lambda==1.0)
400 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 ());
407 ASSERT(nweights==num_kernels);
410 int32_t inner_iters=0;
414 for (int32_t i=0; i<num_kernels; i++)
417 mkl_objective+=old_beta[i]*sumw[i];
451 SG_ERROR(
"Solver type not supported (not compiled in?)\n");
472 int32_t nofKernelsGood;
475 nofKernelsGood = num_kernels;
478 for (p=0; p<num_kernels; ++p )
480 if (sumw[p] >= 0.0 && old_beta[p] >= 0.0 )
482 beta[p] =
CMath::sqrt(sumw[p]*old_beta[p]*old_beta[p]);
499 for( p=0; p<num_kernels; ++p )
504 SG_PRINT(
"MKL-direct: p = %.3f\n", 1.0 );
505 SG_PRINT(
"MKL-direct: nofKernelsGood = %d\n", nofKernelsGood );
506 SG_PRINT(
"MKL-direct: Z = %e\n", Z );
507 SG_PRINT(
"MKL-direct: eps = %e\n", epsRegul );
508 for( p=0; p<num_kernels; ++p )
513 SG_PRINT(
"MKL-direct: preR = %e\n", preR );
514 SG_PRINT(
"MKL-direct: preR/p = %e\n", preR );
516 SG_PRINT(
"MKL-direct: R = %e\n", R );
517 SG_ERROR(
"Assertion R >= 0 failed!\n" );
521 for( p=0; p<num_kernels; ++p )
529 for( p=0; p<num_kernels; ++p )
538 for( p=0; p<num_kernels; ++p )
546 for (p=0; p<num_kernels; ++p )
549 obj += sumw[p] * beta[p];
558 std::list<int32_t> I;
560 for (int32_t i=0; i<len;i++)
570 for (std::list<int32_t>::iterator it=I.begin(); it!=I.end(); it++)
598 for (int32_t i=0; i<n; i++)
602 for (int32_t j=0; j<n; j++)
639 }
while(ff>ff_old+1e-4*gg_old*(del-del_old));
651 SG_ERROR(
"cannot compute objective, labels or kernel not set\n");
666 for( p=0; p<num_kernels; ++p )
680 for( p=0; p<num_kernels; ++p )
685 for( p=0; p<num_kernels; ++p )
686 obj += sumw[p] * beta[p];
702 int32_t nofKernelsGood;
705 nofKernelsGood = num_kernels;
706 for( p=0; p<num_kernels; ++p )
709 if( sumw[p] >= 0.0 && old_beta[p] >= 0.0 )
711 beta[p] = sumw[p] * old_beta[p]*old_beta[p] /
mkl_norm;
724 for( p=0; p<num_kernels; ++p )
729 for( p=0; p<num_kernels; ++p )
734 for( p=0; p<num_kernels; ++p )
735 preR +=
CMath::sq( old_beta[p] - beta[p]);
741 SG_PRINT(
"MKL-direct: nofKernelsGood = %d\n", nofKernelsGood );
742 SG_PRINT(
"MKL-direct: Z = %e\n", Z );
743 SG_PRINT(
"MKL-direct: eps = %e\n", epsRegul );
744 for( p=0; p<num_kernels; ++p )
747 SG_PRINT(
"MKL-direct: t[%3d] = %e ( diff = %e = %e - %e )\n", p, t, old_beta[p]-beta[p], old_beta[p], beta[p] );
749 SG_PRINT(
"MKL-direct: preR = %e\n", preR );
752 SG_PRINT(
"MKL-direct: R = %e\n", R );
753 SG_ERROR(
"Assertion R >= 0 failed!\n" );
757 for( p=0; p<num_kernels; ++p )
765 for( p=0; p<num_kernels; ++p )
775 for( p=0; p<num_kernels; ++p )
776 obj += sumw[p] * beta[p];
782 const float64_t* old_beta, int32_t num_kernels,
789 SG_ERROR(
"MKL via NEWTON works only for norms>1\n");
791 const double epsBeta = 1e-32;
792 const double epsGamma = 1e-12;
793 const double epsWsq = 1e-12;
794 const double epsNewt = 0.0001;
795 const double epsStep = 1e-9;
796 const int nofNewtonSteps = 3;
797 const double hessRidge = 1e-6;
798 const int inLogSpace = 0;
813 for( p=0; p<num_kernels; ++p )
815 beta[p] = old_beta[p];
816 if( !( beta[p] >= epsBeta ) )
819 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
824 if( !( fabs(Z-1.0) <= epsGamma ) )
826 SG_WARNING(
"old_beta not normalized (diff=%e); forcing normalization. ", Z-1.0 );
827 for( p=0; p<num_kernels; ++p )
832 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
838 for ( p=0; p<num_kernels; ++p )
840 if ( !( sumw[p] >= 0 ) )
842 if( !( sumw[p] >= -epsWsq ) )
843 SG_WARNING(
"sumw[%d] = %e; treated as 0. ", p, sumw[p] );
855 if( !( gamma > epsGamma ) )
857 SG_WARNING(
"bad gamma: %e; set to %e. ", gamma, epsGamma );
861 ASSERT( gamma >= epsGamma );
866 for( p=0; p<num_kernels; ++p )
868 obj += beta[p] * sumw[p];
871 if( !( obj >= 0.0 ) )
872 SG_WARNING(
"negative objective: %e. ", obj );
877 for (i = 0; i < nofNewtonSteps; ++i )
882 for( p=0; p<num_kernels; ++p )
884 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
888 const float halfw2p = ( sumw[p] >= 0.0 ) ? (sumw[p]*old_beta[p]*old_beta[p]) : 0.0;
890 const float64_t t1 = ( t0 < 0 ) ? 0.0 : t0;
893 newtDir[p] = t1 / ( t1 + t2*beta[p] + hessRidge );
895 newtDir[p] = ( t1 == 0.0 ) ? 0.0 : ( t1 / t2 );
897 ASSERT( newtDir[p] == newtDir[p] );
905 while( stepSize >= epsStep )
911 for( p=0; p<num_kernels; ++p )
914 newtBeta[p] = beta[p] *
CMath::exp( + stepSize * newtDir[p] );
916 newtBeta[p] = beta[p] + stepSize * newtDir[p];
917 if( !( newtBeta[p] >= epsBeta ) )
918 newtBeta[p] = epsBeta;
930 for( p=0; p<num_kernels; ++p )
933 if( newtBeta[p] > 1.0 )
938 ASSERT( 0.0 <= newtBeta[p] && newtBeta[p] <= 1.0 );
944 for( p=0; p<num_kernels; ++p )
945 newtObj += sumw[p] * old_beta[p]*old_beta[p] / newtBeta[p];
947 if ( newtObj < obj - epsNewt*stepSize*obj )
949 for( p=0; p<num_kernels; ++p )
950 beta[p] = newtBeta[p];
958 if( stepSize < epsStep )
966 for( p=0; p<num_kernels; ++p )
967 obj += beta[p] * sumw[p];
982 int32_t NUMCOLS = 2*num_kernels + 1;
993 for (int32_t i=0; i<2*num_kernels; i++)
1000 for (int32_t i=num_kernels; i<2*num_kernels; i++)
1003 obj[2*num_kernels]=1 ;
1004 lb[2*num_kernels]=-CPX_INFBOUND ;
1005 ub[2*num_kernels]=CPX_INFBOUND ;
1007 int status = CPXnewcols (
env,
lp_cplex, NUMCOLS, obj, lb, ub, NULL, NULL);
1010 CPXgeterrorstring (
env, status, errmsg);
1015 SG_INFO(
"adding the first row\n");
1016 int initial_rmatbeg[1];
1017 int initial_rmatind[num_kernels+1];
1018 double initial_rmatval[num_kernels+1];
1019 double initial_rhs[1];
1020 char initial_sense[1];
1025 initial_rmatbeg[0] = 0;
1027 initial_sense[0]=
'E' ;
1030 for (int32_t i=0; i<num_kernels; i++)
1032 initial_rmatind[i]=i ;
1033 initial_rmatval[i]=1 ;
1035 initial_rmatind[num_kernels]=2*num_kernels ;
1036 initial_rmatval[num_kernels]=0 ;
1038 status = CPXaddrows (
env,
lp_cplex, 0, 1, num_kernels+1,
1039 initial_rhs, initial_sense, initial_rmatbeg,
1040 initial_rmatind, initial_rmatval, NULL, NULL);
1045 initial_rmatbeg[0] = 0;
1047 initial_sense[0]=
'L' ;
1049 initial_rmatind[0]=2*num_kernels ;
1050 initial_rmatval[0]=0 ;
1053 initial_rhs, initial_sense, initial_rmatbeg,
1054 initial_rmatind, initial_rmatval, NULL, NULL);
1059 for (int32_t i=0; i<num_kernels; i++)
1061 initial_rmatind[i]=i ;
1062 initial_rmatval[i]=1 ;
1064 initial_rmatind[num_kernels]=2*num_kernels ;
1065 initial_rmatval[num_kernels]=0 ;
1067 status = CPXaddqconstr (
env,
lp_cplex, 0, num_kernels+1, 1.0,
'L', NULL, NULL,
1068 initial_rmatind, initial_rmatind, initial_rmatval, NULL);
1074 SG_ERROR(
"Failed to add the first row.\n");
1080 for (int32_t q=0; q<num_kernels-1; q++)
1097 rmatind[2]=num_kernels+q ;
1100 rhs, sense, rmatbeg,
1101 rmatind, rmatval, NULL, NULL);
1103 SG_ERROR(
"Failed to add a smothness row (1).\n");
1112 rmatind[2]=num_kernels+q ;
1115 rhs, sense, rmatbeg,
1116 rmatind, rmatval, NULL, NULL);
1118 SG_ERROR(
"Failed to add a smothness row (2).\n");
1127 int rmatind[num_kernels+1];
1128 double rmatval[num_kernels+1];
1140 for (int32_t i=0; i<num_kernels; i++)
1144 rmatval[i]=-(sumw[i]-suma) ;
1146 rmatval[i]=-sumw[i];
1148 rmatind[num_kernels]=2*num_kernels ;
1149 rmatval[num_kernels]=-1 ;
1151 int32_t status = CPXaddrows (
env,
lp_cplex, 0, 1, num_kernels+1,
1152 rhs, sense, rmatbeg,
1153 rmatind, rmatval, NULL, NULL);
1155 SG_ERROR(
"Failed to add the new row.\n");
1170 for (int32_t i=0; i<num_kernels; i++)
1171 beta[i]=old_beta[i];
1172 for (int32_t i=num_kernels; i<2*num_kernels+1; i++)
1186 SG_ERROR(
"Failed to optimize Problem.\n");
1190 status=CPXsolution(
env,
lp_cplex, &solstat, &objval,
1191 (
double*) beta, NULL, NULL, NULL);
1196 SG_ERROR(
"Failed to obtain solution.\n");
1213 SG_ERROR(
"Failed to optimize Problem.\n");
1216 int32_t cur_numrows=(int32_t) CPXgetnumrows(
env,
lp_cplex);
1217 int32_t cur_numcols=(int32_t) CPXgetnumcols(
env,
lp_cplex);
1218 int32_t num_rows=cur_numrows;
1219 ASSERT(cur_numcols<=2*num_kernels+1);
1230 status=CPXsolution(
env,
lp_cplex, &solstat, &objval,
1231 (
double*) x, (
double*) pi, (
double*) slack, NULL);
1235 status=CPXsolution(
env,
lp_cplex, &solstat, &objval,
1236 (
double*) x, NULL, (
double*) slack, NULL);
1239 int32_t solution_ok = (!status) ;
1241 SG_ERROR(
"Failed to obtain solution.\n");
1243 int32_t num_active_rows=0 ;
1248 int32_t max_idx = -1 ;
1249 int32_t start_row = 1 ;
1251 start_row+=2*(num_kernels-1);
1253 for (int32_t i = start_row; i < cur_numrows; i++)
1261 if (slack[i]>max_slack)
1263 max_slack=slack[i] ;
1274 if (slack[i]>max_slack)
1276 max_slack=slack[i] ;
1284 if ( (num_rows-start_row>
CMath::max(100,2*num_active_rows)) && (max_idx!=-1))
1287 status = CPXdelrows (
env,
lp_cplex, max_idx, max_idx) ;
1289 SG_ERROR(
"Failed to remove an old row.\n");
1294 rho = -x[2*num_kernels] ;
1306 for (int32_t i=0; i<num_kernels; i++)
1311 SG_ERROR(
"Cplex not enabled at compile time\n");
1322 SG_ERROR(
"MKL via GLPK works only for norm=1\n");
1326 int32_t NUMCOLS = 2*num_kernels + 1 ;
1332 lpx_add_cols(
lp_glpk, NUMCOLS);
1333 for (
int i=1; i<=2*num_kernels; i++)
1335 lpx_set_obj_coef(
lp_glpk, i, 0);
1336 lpx_set_col_bnds(
lp_glpk, i, LPX_DB, 0, 1);
1338 for (
int i=num_kernels+1; i<=2*num_kernels; i++)
1342 lpx_set_obj_coef(
lp_glpk, NUMCOLS, 1);
1343 lpx_set_col_bnds(
lp_glpk, NUMCOLS, LPX_FR, 0,0);
1346 int row_index = lpx_add_rows(
lp_glpk, 1);
1347 int* ind =
SG_MALLOC(
int, num_kernels+2);
1349 for (
int i=1; i<=num_kernels; i++)
1354 ind[num_kernels+1] = NUMCOLS;
1355 val[num_kernels+1] = 0;
1356 lpx_set_mat_row(
lp_glpk, row_index, num_kernels, ind, val);
1357 lpx_set_row_bnds(
lp_glpk, row_index, LPX_FX, 1, 1);
1365 for (int32_t q=1; q<num_kernels; q++)
1369 int mat_row_index = lpx_add_rows(
lp_glpk, 2);
1374 mat_ind[3] = num_kernels+q;
1376 lpx_set_mat_row(
lp_glpk, mat_row_index, 3, mat_ind, mat_val);
1377 lpx_set_row_bnds(
lp_glpk, mat_row_index, LPX_UP, 0, 0);
1380 lpx_set_mat_row(
lp_glpk, mat_row_index+1, 3, mat_ind, mat_val);
1381 lpx_set_row_bnds(
lp_glpk, mat_row_index+1, LPX_UP, 0, 0);
1388 int row_index = lpx_add_rows(
lp_glpk, 1);
1389 for (int32_t i=1; i<=num_kernels; i++)
1392 val[i] = -(sumw[i-1]-suma);
1394 ind[num_kernels+1] = 2*num_kernels+1;
1395 val[num_kernels+1] = -1;
1396 lpx_set_mat_row(
lp_glpk, row_index, num_kernels+1, ind, val);
1397 lpx_set_row_bnds(
lp_glpk, row_index, LPX_UP, 0, 0);
1405 SG_ERROR(
"Failed to optimize Problem.\n");
1407 int32_t cur_numrows = lpx_get_num_rows(
lp_glpk);
1408 int32_t cur_numcols = lpx_get_num_cols(
lp_glpk);
1409 int32_t num_rows=cur_numrows;
1410 ASSERT(cur_numcols<=2*num_kernels+1);
1416 for (
int i=0; i<cur_numrows; i++)
1418 row_primal[i] = lpx_get_row_prim(
lp_glpk, i+1);
1419 row_dual[i] = lpx_get_row_dual(
lp_glpk, i+1);
1421 for (
int i=0; i<cur_numcols; i++)
1422 col_primal[i] = lpx_get_col_prim(
lp_glpk, i+1);
1424 obj = -col_primal[2*num_kernels];
1426 for (
int i=0; i<num_kernels; i++)
1427 beta[i] = col_primal[i];
1429 int32_t num_active_rows=0;
1433 int32_t max_idx = -1;
1434 int32_t start_row = 1;
1436 start_row += 2*(num_kernels-1);
1438 for (int32_t i= start_row; i<cur_numrows; i++)
1444 if (row_primal[i]<max_slack)
1446 max_slack = row_primal[i];
1452 if ((num_rows-start_row>
CMath::max(100, 2*num_active_rows)) && max_idx!=-1)
1455 del_rows[1] = max_idx+1;
1456 lpx_del_rows(
lp_glpk, 1, del_rows);
1464 SG_ERROR(
"Glpk not enabled at compile time\n");
1480 ASSERT(nweights==num_kernels);
1483 for (int32_t i=0; i<num_kernels; i++)
1489 for (int32_t n=0; n<num_kernels; n++)
1494 for (int32_t i=0; i<nsv; i++)
1498 for (int32_t j=0; j<nsv; j++)
1530 for (int32_t i=0; i<n; i++)
1534 for (int32_t j=0; j<n; j++)
1551 mkl_obj=-0.5*mkl_obj;
1558 SG_ERROR(
"cannot compute objective, labels or kernel not set\n");
1579 for (int32_t i=0; i<num_kernels; i++)
1583 lin_term[i]=grad_beta[i] - 2*beta[i]*hess_beta[i];
1584 const_term+=grad_beta[i]*beta[i] -
CMath::sq(beta[i])*hess_beta[i];
1587 ind[num_kernels]=2*num_kernels;
1588 hess_beta[num_kernels]=0;
1589 lin_term[num_kernels]=0;
1600 status = CPXaddqconstr (
env,
lp_cplex, num_kernels+1, num_kernels+1, const_term,
'L', ind, lin_term,
1601 ind, ind, hess_beta, NULL);