00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "lib/Signal.h"
00013 #include "classifier/mkl/MKL.h"
00014 #include "classifier/svm/LibSVM.h"
00015 #include "kernel/CombinedKernel.h"
00016
00017 using namespace shogun;
00018
00019 CMKL::CMKL(CSVM* s)
00020 : CSVM(), svm(NULL), C_mkl(0), mkl_norm(1),
00021 mkl_iterations(0), mkl_epsilon(1e-5), interleaved_optimization(true),
00022 w_gap(1.0), rho(0)
00023 {
00024 set_constraint_generator(s);
00025 #ifdef USE_CPLEX
00026 lp_cplex = NULL ;
00027 env = NULL ;
00028 #endif
00029
00030 #ifdef USE_GLPK
00031 lp_glpk = NULL;
00032 #endif
00033
00034 SG_DEBUG("creating MKL object %p\n", this);
00035 lp_initialized = false ;
00036 }
00037
00038 CMKL::~CMKL()
00039 {
00040 SG_DEBUG("deleting MKL object %p\n", this);
00041 if (svm)
00042 svm->set_callback_function(NULL, NULL);
00043 SG_UNREF(svm);
00044 }
00045
00046 void CMKL::init_solver()
00047 {
00048 #ifdef USE_CPLEX
00049 cleanup_cplex();
00050
00051 if (get_solver_type()==ST_CPLEX)
00052 init_cplex();
00053 #endif
00054
00055 #ifdef USE_GLPK
00056 cleanup_glpk();
00057
00058 if (get_solver_type() == ST_GLPK)
00059 init_glpk();
00060 #endif
00061 }
00062
00063 #ifdef USE_CPLEX
00064 bool CMKL::init_cplex()
00065 {
00066 while (env==NULL)
00067 {
00068 SG_INFO( "trying to initialize CPLEX\n") ;
00069
00070 int status = 0;
00071 env = CPXopenCPLEX (&status);
00072
00073 if ( env == NULL )
00074 {
00075 char errmsg[1024];
00076 SG_WARNING( "Could not open CPLEX environment.\n");
00077 CPXgeterrorstring (env, status, errmsg);
00078 SG_WARNING( "%s", errmsg);
00079 SG_WARNING( "retrying in 60 seconds\n");
00080 sleep(60);
00081 }
00082 else
00083 {
00084
00085 status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, CPX_ALG_DUAL);
00086 if ( status )
00087 {
00088 SG_ERROR( "Failure to select dual lp optimization, error %d.\n", status);
00089 }
00090 else
00091 {
00092 status = CPXsetintparam (env, CPX_PARAM_DATACHECK, CPX_ON);
00093 if ( status )
00094 {
00095 SG_ERROR( "Failure to turn on data checking, error %d.\n", status);
00096 }
00097 else
00098 {
00099 lp_cplex = CPXcreateprob (env, &status, "light");
00100
00101 if ( lp_cplex == NULL )
00102 SG_ERROR( "Failed to create LP.\n");
00103 else
00104 CPXchgobjsen (env, lp_cplex, CPX_MIN);
00105 }
00106 }
00107 }
00108 }
00109
00110 return (lp_cplex != NULL) && (env != NULL);
00111 }
00112
00113 bool CMKL::cleanup_cplex()
00114 {
00115 bool result=false;
00116
00117 if (lp_cplex)
00118 {
00119 int32_t status = CPXfreeprob(env, &lp_cplex);
00120 lp_cplex = NULL;
00121 lp_initialized = false;
00122
00123 if (status)
00124 SG_WARNING( "CPXfreeprob failed, error code %d.\n", status);
00125 else
00126 result = true;
00127 }
00128
00129 if (env)
00130 {
00131 int32_t status = CPXcloseCPLEX (&env);
00132 env=NULL;
00133
00134 if (status)
00135 {
00136 char errmsg[1024];
00137 SG_WARNING( "Could not close CPLEX environment.\n");
00138 CPXgeterrorstring (env, status, errmsg);
00139 SG_WARNING( "%s", errmsg);
00140 }
00141 else
00142 result = true;
00143 }
00144 return result;
00145 }
00146 #endif
00147
00148 #ifdef USE_GLPK
00149 bool CMKL::init_glpk()
00150 {
00151 lp_glpk = lpx_create_prob();
00152 lpx_set_obj_dir(lp_glpk, LPX_MIN);
00153 lpx_set_int_parm(lp_glpk, LPX_K_DUAL, GLP_ON );
00154 lpx_set_int_parm(lp_glpk, LPX_K_PRESOL, GLP_ON );
00155
00156 glp_term_out(GLP_OFF);
00157 return (lp_glpk != NULL);
00158 }
00159
00160 bool CMKL::cleanup_glpk()
00161 {
00162 lp_initialized = false;
00163 if (lp_glpk)
00164 lpx_delete_prob(lp_glpk);
00165 lp_glpk = NULL;
00166 return true;
00167 }
00168
00169 bool CMKL::check_lpx_status(LPX *lp)
00170 {
00171 int status = lpx_get_status(lp);
00172
00173 if (status==LPX_INFEAS)
00174 {
00175 SG_PRINT("solution is infeasible!\n");
00176 return false;
00177 }
00178 else if(status==LPX_NOFEAS)
00179 {
00180 SG_PRINT("problem has no feasible solution!\n");
00181 return false;
00182 }
00183 return true;
00184 }
00185 #endif // USE_GLPK
00186
00187 bool CMKL::train(CFeatures* data)
00188 {
00189 ASSERT(kernel);
00190 ASSERT(labels && labels->get_num_labels());
00191
00192 if (data)
00193 {
00194 if (labels->get_num_labels() != data->get_num_vectors())
00195 SG_ERROR("Number of training vectors does not match number of labels\n");
00196 kernel->init(data, data);
00197 }
00198
00199 init_training();
00200 if (!svm)
00201 SG_ERROR("No constraint generator (SVM) set\n");
00202
00203 int32_t num_label=0;
00204 if (labels)
00205 labels->get_num_labels();
00206
00207 SG_INFO("%d trainlabels\n", num_label);
00208 if (mkl_epsilon<=0)
00209 mkl_epsilon=1e-2 ;
00210
00211 SG_INFO("mkl_epsilon = %1.1e\n", mkl_epsilon) ;
00212 SG_INFO("C_mkl = %1.1e\n", C_mkl) ;
00213 SG_INFO("mkl_norm = %1.3e\n", mkl_norm);
00214 SG_INFO("solver = %d\n", get_solver_type());
00215
00216 int32_t num_weights = -1;
00217 int32_t num_kernels = kernel->get_num_subkernels();
00218 SG_INFO("num_kernels = %d\n", num_kernels);
00219 const float64_t* beta_const = kernel->get_subkernel_weights(num_weights);
00220 float64_t* beta = CMath::clone_vector(beta_const, num_weights);
00221 ASSERT(num_weights==num_kernels);
00222 CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm), beta, num_kernels);
00223 kernel->set_subkernel_weights(beta, num_kernels);
00224 delete[] beta;
00225
00226 svm->set_bias_enabled(get_bias_enabled());
00227 svm->set_epsilon(get_epsilon());
00228 svm->set_max_train_time(get_max_train_time());
00229 svm->set_nu(get_nu());
00230 svm->set_C(get_C1(), get_C2());
00231 svm->set_qpsize(get_qpsize());
00232 svm->set_shrinking_enabled(get_shrinking_enabled());
00233 svm->set_linadd_enabled(get_linadd_enabled());
00234 svm->set_batch_computation_enabled(get_batch_computation_enabled());
00235 svm->set_labels(labels);
00236 svm->set_kernel(kernel);
00237
00238 #ifdef USE_CPLEX
00239 cleanup_cplex();
00240
00241 if (get_solver_type()==ST_CPLEX)
00242 init_cplex();
00243 #endif
00244
00245 #ifdef USE_GLPK
00246 if (get_solver_type()==ST_GLPK)
00247 init_glpk();
00248 #endif
00249
00250 mkl_iterations = 0;
00251 CSignal::clear_cancel();
00252
00253 if (interleaved_optimization)
00254 {
00255 if (svm->get_classifier_type() != CT_LIGHT && svm->get_classifier_type() != CT_SVRLIGHT)
00256 {
00257 SG_ERROR("Interleaved MKL optimization is currently "
00258 "only supported with SVMlight\n");
00259 }
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271 int32_t refs=this->ref();
00272 svm->set_callback_function(this, perform_mkl_step_helper);
00273 svm->train();
00274 SG_DONE();
00275 svm->set_callback_function(NULL, NULL);
00276 if (refs>1)
00277 this->unref();
00278 }
00279 else
00280 {
00281 float64_t* sumw = new float64_t[num_kernels];
00282
00283 while (true)
00284 {
00285 svm->train();
00286
00287 float64_t suma=compute_sum_alpha();
00288 compute_sum_beta(sumw);
00289
00290 mkl_iterations++;
00291 if (perform_mkl_step(sumw, suma) || CSignal::cancel_computations())
00292 break;
00293 }
00294
00295 delete[] sumw;
00296 }
00297 #ifdef USE_CPLEX
00298 cleanup_cplex();
00299 #endif
00300 #ifdef USE_GLPK
00301 cleanup_glpk();
00302 #endif
00303
00304 int32_t nsv=svm->get_num_support_vectors();
00305 create_new_model(nsv);
00306
00307 set_bias(svm->get_bias());
00308 for (int32_t i=0; i<nsv; i++)
00309 {
00310 set_alpha(i, svm->get_alpha(i));
00311 set_support_vector(i, svm->get_support_vector(i));
00312 }
00313 return true;
00314 }
00315
00316 bool CMKL::perform_mkl_step(
00317 const float64_t* sumw, float64_t suma)
00318 {
00319 int32_t num_kernels = kernel->get_num_subkernels();
00320 int32_t nweights=0;
00321 const float64_t* old_beta = kernel->get_subkernel_weights(nweights);
00322 ASSERT(nweights==num_kernels);
00323 float64_t* beta = new float64_t[num_kernels];
00324
00325 int32_t inner_iters=0;
00326 float64_t mkl_objective=0;
00327
00328 mkl_objective=-suma;
00329 for (int32_t i=0; i<num_kernels; i++)
00330 {
00331 beta[i]=old_beta[i];
00332 mkl_objective+=old_beta[i]*sumw[i];
00333 }
00334
00335 w_gap = CMath::abs(1-rho/mkl_objective) ;
00336
00337 if( (w_gap >= mkl_epsilon) ||
00338 (get_solver_type()==ST_AUTO || get_solver_type()==ST_NEWTON || get_solver_type()==ST_DIRECT ) )
00339 {
00340 if (get_solver_type()==ST_AUTO || get_solver_type()==ST_DIRECT)
00341 rho=compute_optimal_betas_directly(beta, old_beta, num_kernels, sumw, suma, mkl_objective);
00342 else if (get_solver_type()==ST_NEWTON)
00343 rho=compute_optimal_betas_newton(beta, old_beta, num_kernels, sumw, suma, mkl_objective);
00344 #ifdef USE_CPLEX
00345 else if (get_solver_type()==ST_CPLEX)
00346 rho=compute_optimal_betas_via_cplex(beta, old_beta, num_kernels, sumw, suma, inner_iters);
00347 #endif
00348 #ifdef USE_GLPK
00349 else if (get_solver_type()==ST_GLPK)
00350 rho=compute_optimal_betas_via_glpk(beta, old_beta, num_kernels, sumw, suma, inner_iters);
00351 #endif
00352 else
00353 SG_ERROR("Solver type not supported (not compiled in?)\n");
00354 w_gap = CMath::abs(1-rho/mkl_objective) ;
00355 }
00356
00357 kernel->set_subkernel_weights(beta, num_kernels);
00358 delete[] beta;
00359
00360 return converged();
00361 }
00362
00363
00364 float64_t CMKL::compute_optimal_betas_directly(
00365 float64_t* beta, const float64_t* old_beta, const int32_t num_kernels,
00366 const float64_t* sumw, const float64_t suma,
00367 const float64_t mkl_objective )
00368 {
00369 const float64_t epsRegul = 0.01;
00370 float64_t obj;
00371 float64_t Z;
00372 float64_t preR;
00373 int32_t p;
00374 int32_t nofKernelsGood;
00375
00376
00377 nofKernelsGood = num_kernels;
00378 for( p=0; p<num_kernels; ++p ) {
00379
00380 if( sumw[p] >= 0.0 && old_beta[p] >= 0.0 ) {
00381 beta[p] = sumw[p] * old_beta[p]*old_beta[p] / mkl_norm;
00382 beta[p] = CMath::pow( beta[p], 1.0 / (mkl_norm+1.0) );
00383 } else {
00384 beta[p] = 0.0;
00385 --nofKernelsGood;
00386 }
00387 ASSERT( beta[p] >= 0 );
00388 }
00389
00390
00391 Z = 0.0;
00392 for( p=0; p<num_kernels; ++p ) {
00393 Z += CMath::pow( beta[p], mkl_norm );
00394 }
00395 Z = CMath::pow( Z, -1.0/mkl_norm );
00396 ASSERT( Z >= 0 );
00397 for( p=0; p<num_kernels; ++p ) {
00398 beta[p] *= Z;
00399 }
00400
00401
00402 preR = 0.0;
00403 for( p=0; p<num_kernels; ++p ) {
00404 preR += CMath::pow( old_beta[p] - beta[p], 2.0 );
00405 }
00406 const float64_t R = CMath::sqrt( preR / mkl_norm ) * epsRegul;
00407 if( !( R >= 0 ) ) {
00408 SG_PRINT( "MKL-direct: p = %.3f\n", mkl_norm );
00409 SG_PRINT( "MKL-direct: nofKernelsGood = %d\n", nofKernelsGood );
00410 SG_PRINT( "MKL-direct: Z = %e\n", Z );
00411 SG_PRINT( "MKL-direct: eps = %e\n", epsRegul );
00412 for( p=0; p<num_kernels; ++p ) {
00413 const float64_t t = CMath::pow( old_beta[p] - beta[p], 2.0 );
00414 SG_PRINT( "MKL-direct: t[%3d] = %e ( diff = %e = %e - %e )\n", p, t, old_beta[p]-beta[p], old_beta[p], beta[p] );
00415 }
00416 SG_PRINT( "MKL-direct: preR = %e\n", preR );
00417 SG_PRINT( "MKL-direct: preR/p = %e\n", preR/mkl_norm );
00418 SG_PRINT( "MKL-direct: sqrt(preR/p) = %e\n", CMath::sqrt(preR/mkl_norm) );
00419 SG_PRINT( "MKL-direct: R = %e\n", R );
00420 SG_ERROR( "Assertion R >= 0 failed!\n" );
00421 }
00422
00423 Z = 0.0;
00424 for( p=0; p<num_kernels; ++p ) {
00425 beta[p] += R;
00426 Z += CMath::pow( beta[p], mkl_norm );
00427 ASSERT( beta[p] >= 0 );
00428 }
00429 Z = CMath::pow( Z, -1.0/mkl_norm );
00430 ASSERT( Z >= 0 );
00431 for( p=0; p<num_kernels; ++p ) {
00432 beta[p] *= Z;
00433 ASSERT( beta[p] >= 0.0 );
00434 if( beta[p] > 1.0 ) {
00435 beta[p] = 1.0;
00436 }
00437 }
00438
00439
00440 obj = -suma;
00441 for( p=0; p<num_kernels; ++p ) {
00442
00443 obj += sumw[p] * beta[p];
00444 }
00445 return obj;
00446 }
00447
00448 float64_t CMKL::compute_optimal_betas_newton(float64_t* beta,
00449 const float64_t* old_beta, int32_t num_kernels,
00450 const float64_t* sumw, float64_t suma,
00451 float64_t mkl_objective)
00452 {
00453 SG_DEBUG("MKL via NEWTON\n");
00454
00455 if (mkl_norm==1)
00456 SG_ERROR("MKL via NEWTON works only for norms>1\n");
00457
00458 const double epsBeta = 1e-32;
00459 const double epsGamma = 1e-12;
00460 const double epsWsq = 1e-12;
00461 const double epsNewt = 0.0001;
00462 const double epsStep = 1e-9;
00463 const int nofNewtonSteps = 3;
00464 const double hessRidge = 1e-6;
00465 const int inLogSpace = 0;
00466
00467 const float64_t r = mkl_norm / ( mkl_norm - 1.0 );
00468 float64_t* newtDir = new float64_t[ num_kernels ];
00469 float64_t* newtBeta = new float64_t[ num_kernels ];
00470 float64_t newtStep;
00471 float64_t stepSize;
00472 float64_t Z;
00473 float64_t obj;
00474 float64_t gamma;
00475 int32_t p;
00476 int i;
00477
00478
00479
00480 for( p=0; p<num_kernels; ++p ) {
00481
00482 }
00483
00484
00485 Z = 0.0;
00486 for( p=0; p<num_kernels; ++p ) {
00487 beta[p] = old_beta[p];
00488 if( !( beta[p] >= epsBeta ) ) {
00489
00490 beta[p] = epsBeta;
00491 }
00492 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
00493 Z += CMath::pow( beta[p], mkl_norm );
00494 }
00495 Z = CMath::pow( Z, -1.0/mkl_norm );
00496 if( !( fabs(Z-1.0) <= epsGamma ) ) {
00497 SG_WARNING( "old_beta not normalized (diff=%e); forcing normalization. ", Z-1.0 );
00498 for( p=0; p<num_kernels; ++p ) {
00499 beta[p] *= Z;
00500 if( beta[p] > 1.0 ) {
00501 beta[p] = 1.0;
00502 }
00503 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
00504 }
00505 }
00506
00507
00508 gamma = 0.0;
00509 for( p=0; p<num_kernels; ++p ) {
00510 if( !( sumw[p] >= 0 ) ) {
00511 if( !( sumw[p] >= -epsWsq ) ) {
00512 SG_WARNING( "sumw[%d] = %e; treated as 0. ", p, sumw[p] );
00513 }
00514
00515 } else {
00516 ASSERT( sumw[p] >= 0 );
00517
00518 gamma += CMath::pow( sumw[p] * beta[p]*beta[p] / mkl_norm, r );
00519 }
00520 }
00521 gamma = CMath::pow( gamma, 1.0/r ) / mkl_norm;
00522 ASSERT( gamma > -1e-9 );
00523 if( !( gamma > epsGamma ) ) {
00524 SG_WARNING( "bad gamma: %e; set to %e. ", gamma, epsGamma );
00525
00526 gamma = epsGamma;
00527 }
00528 ASSERT( gamma >= epsGamma );
00529
00530
00531
00532 obj = 0.0;
00533 for( p=0; p<num_kernels; ++p ) {
00534 obj += beta[p] * sumw[p];
00535
00536 }
00537 if( !( obj >= 0.0 ) ) {
00538 SG_WARNING( "negative objective: %e. ", obj );
00539 }
00540
00541
00542
00543
00544 if( nofNewtonSteps > 1 ) {
00545
00546 }
00547 for( i = 0; i < nofNewtonSteps; ++i ) {
00548
00549
00550 const float64_t gqq1 = mkl_norm * (mkl_norm-1.0) * gamma;
00551 newtStep = 0.0;
00552 for( p=0; p<num_kernels; ++p ) {
00553 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
00554
00555
00556
00557 const float halfw2p = ( sumw[p] >= 0.0 ) ? (sumw[p]*old_beta[p]*old_beta[p]) : 0.0;
00558 const float64_t t0 = halfw2p*beta[p] - mkl_norm*gamma*CMath::pow(beta[p],mkl_norm+2.0);
00559 const float64_t t1 = ( t0 < 0 ) ? 0.0 : t0;
00560 const float64_t t2 = 2.0*halfw2p + gqq1*CMath::pow(beta[p],mkl_norm+1.0);
00561 if( inLogSpace ) {
00562 newtDir[p] = t1 / ( t1 + t2*beta[p] + hessRidge );
00563 } else {
00564 newtDir[p] = ( t1 == 0.0 ) ? 0.0 : ( t1 / t2 );
00565 }
00566
00567 ASSERT( newtDir[p] == newtDir[p] );
00568
00569 }
00570
00571
00572
00573
00574 stepSize = 1.0;
00575 while( stepSize >= epsStep ) {
00576
00577
00578 Z = 0.0;
00579 while( Z == 0.0 ) {
00580 for( p=0; p<num_kernels; ++p ) {
00581 if( inLogSpace ) {
00582 newtBeta[p] = beta[p] * CMath::exp( + stepSize * newtDir[p] );
00583 } else {
00584 newtBeta[p] = beta[p] + stepSize * newtDir[p];
00585 }
00586 if( !( newtBeta[p] >= epsBeta ) ) {
00587 newtBeta[p] = epsBeta;
00588 }
00589 Z += CMath::pow( newtBeta[p], mkl_norm );
00590 }
00591 ASSERT( 0.0 <= Z );
00592 Z = CMath::pow( Z, -1.0/mkl_norm );
00593 if( Z == 0.0 ) {
00594 stepSize /= 2.0;
00595 }
00596 }
00597
00598
00599 ASSERT( 0.0 < Z );
00600 ASSERT( Z < CMath::INFTY );
00601 for( p=0; p<num_kernels; ++p ) {
00602 newtBeta[p] *= Z;
00603 if( newtBeta[p] > 1.0 ) {
00604
00605 newtBeta[p] = 1.0;
00606 }
00607 ASSERT( 0.0 <= newtBeta[p] && newtBeta[p] <= 1.0 );
00608 }
00609
00610
00611 float64_t newtObj;
00612 newtObj = 0.0;
00613 for( p=0; p<num_kernels; ++p ) {
00614 newtObj += sumw[p] * old_beta[p]*old_beta[p] / newtBeta[p];
00615 }
00616
00617 if( newtObj < obj - epsNewt*stepSize*obj ) {
00618 for( p=0; p<num_kernels; ++p ) {
00619 beta[p] = newtBeta[p];
00620 }
00621 obj = newtObj;
00622 break;
00623 }
00624 stepSize /= 2.0;
00625
00626 }
00627
00628 if( stepSize < epsStep ) {
00629 break;
00630 }
00631
00632 }
00633 delete[] newtDir;
00634 delete[] newtBeta;
00635
00636
00637
00638
00639 obj = -suma;
00640 for( p=0; p<num_kernels; ++p ) {
00641 obj += beta[p] * sumw[p];
00642
00643 }
00644 return obj;
00645 }
00646
00647
00648
00649 float64_t CMKL::compute_optimal_betas_via_cplex(float64_t* new_beta, const float64_t* old_beta, int32_t num_kernels,
00650 const float64_t* sumw, float64_t suma, int32_t& inner_iters)
00651 {
00652 SG_DEBUG("MKL via CPLEX\n");
00653
00654 #ifdef USE_CPLEX
00655 ASSERT(new_beta);
00656 ASSERT(old_beta);
00657
00658 int32_t NUMCOLS = 2*num_kernels + 1;
00659 double* x=new double[NUMCOLS];
00660
00661 if (!lp_initialized)
00662 {
00663 SG_INFO( "creating LP\n") ;
00664
00665 double obj[NUMCOLS];
00666 double lb[NUMCOLS];
00667 double ub[NUMCOLS];
00668
00669 for (int32_t i=0; i<2*num_kernels; i++)
00670 {
00671 obj[i]=0 ;
00672 lb[i]=0 ;
00673 ub[i]=1 ;
00674 }
00675
00676 for (int32_t i=num_kernels; i<2*num_kernels; i++)
00677 obj[i]= C_mkl;
00678
00679 obj[2*num_kernels]=1 ;
00680 lb[2*num_kernels]=-CPX_INFBOUND ;
00681 ub[2*num_kernels]=CPX_INFBOUND ;
00682
00683 int status = CPXnewcols (env, lp_cplex, NUMCOLS, obj, lb, ub, NULL, NULL);
00684 if ( status ) {
00685 char errmsg[1024];
00686 CPXgeterrorstring (env, status, errmsg);
00687 SG_ERROR( "%s", errmsg);
00688 }
00689
00690
00691 SG_INFO( "adding the first row\n");
00692 int initial_rmatbeg[1];
00693 int initial_rmatind[num_kernels+1];
00694 double initial_rmatval[num_kernels+1];
00695 double initial_rhs[1];
00696 char initial_sense[1];
00697
00698
00699 if (mkl_norm==1)
00700 {
00701 initial_rmatbeg[0] = 0;
00702 initial_rhs[0]=1 ;
00703 initial_sense[0]='E' ;
00704
00705 for (int32_t i=0; i<num_kernels; i++)
00706 {
00707 initial_rmatind[i]=i ;
00708 initial_rmatval[i]=1 ;
00709 }
00710 initial_rmatind[num_kernels]=2*num_kernels ;
00711 initial_rmatval[num_kernels]=0 ;
00712
00713 status = CPXaddrows (env, lp_cplex, 0, 1, num_kernels+1,
00714 initial_rhs, initial_sense, initial_rmatbeg,
00715 initial_rmatind, initial_rmatval, NULL, NULL);
00716
00717 }
00718 else
00719 {
00720 initial_rmatbeg[0] = 0;
00721 initial_rhs[0]=0 ;
00722 initial_sense[0]='L' ;
00723
00724 initial_rmatind[0]=2*num_kernels ;
00725 initial_rmatval[0]=0 ;
00726
00727 status = CPXaddrows (env, lp_cplex, 0, 1, 1,
00728 initial_rhs, initial_sense, initial_rmatbeg,
00729 initial_rmatind, initial_rmatval, NULL, NULL);
00730
00731
00732 if (mkl_norm==2)
00733 {
00734 for (int32_t i=0; i<num_kernels; i++)
00735 {
00736 initial_rmatind[i]=i ;
00737 initial_rmatval[i]=1 ;
00738 }
00739 initial_rmatind[num_kernels]=2*num_kernels ;
00740 initial_rmatval[num_kernels]=0 ;
00741
00742 status = CPXaddqconstr (env, lp_cplex, 0, num_kernels+1, 1.0, 'L', NULL, NULL,
00743 initial_rmatind, initial_rmatind, initial_rmatval, NULL);
00744 }
00745 }
00746
00747
00748 if ( status )
00749 SG_ERROR( "Failed to add the first row.\n");
00750
00751 lp_initialized = true ;
00752
00753 if (C_mkl!=0.0)
00754 {
00755 for (int32_t q=0; q<num_kernels-1; q++)
00756 {
00757
00758
00759 int rmatbeg[1];
00760 int rmatind[3];
00761 double rmatval[3];
00762 double rhs[1];
00763 char sense[1];
00764
00765 rmatbeg[0] = 0;
00766 rhs[0]=0 ;
00767 sense[0]='L' ;
00768 rmatind[0]=q ;
00769 rmatval[0]=1 ;
00770 rmatind[1]=q+1 ;
00771 rmatval[1]=-1 ;
00772 rmatind[2]=num_kernels+q ;
00773 rmatval[2]=-1 ;
00774 status = CPXaddrows (env, lp_cplex, 0, 1, 3,
00775 rhs, sense, rmatbeg,
00776 rmatind, rmatval, NULL, NULL);
00777 if ( status )
00778 SG_ERROR( "Failed to add a smothness row (1).\n");
00779
00780 rmatbeg[0] = 0;
00781 rhs[0]=0 ;
00782 sense[0]='L' ;
00783 rmatind[0]=q ;
00784 rmatval[0]=-1 ;
00785 rmatind[1]=q+1 ;
00786 rmatval[1]=1 ;
00787 rmatind[2]=num_kernels+q ;
00788 rmatval[2]=-1 ;
00789 status = CPXaddrows (env, lp_cplex, 0, 1, 3,
00790 rhs, sense, rmatbeg,
00791 rmatind, rmatval, NULL, NULL);
00792 if ( status )
00793 SG_ERROR( "Failed to add a smothness row (2).\n");
00794 }
00795 }
00796 }
00797
00798 {
00799
00800
00801 int rmatbeg[1];
00802 int rmatind[num_kernels+1];
00803 double rmatval[num_kernels+1];
00804 double rhs[1];
00805 char sense[1];
00806
00807 rmatbeg[0] = 0;
00808 if (mkl_norm==1)
00809 rhs[0]=0 ;
00810 else
00811 rhs[0]=-suma ;
00812
00813 sense[0]='L' ;
00814
00815 for (int32_t i=0; i<num_kernels; i++)
00816 {
00817 rmatind[i]=i ;
00818 if (mkl_norm==1)
00819 rmatval[i]=-(sumw[i]-suma) ;
00820 else
00821 rmatval[i]=-sumw[i];
00822 }
00823 rmatind[num_kernels]=2*num_kernels ;
00824 rmatval[num_kernels]=-1 ;
00825
00826 int32_t status = CPXaddrows (env, lp_cplex, 0, 1, num_kernels+1,
00827 rhs, sense, rmatbeg,
00828 rmatind, rmatval, NULL, NULL);
00829 if ( status )
00830 SG_ERROR( "Failed to add the new row.\n");
00831 }
00832
00833 inner_iters=0;
00834 int status;
00835 {
00836
00837 if (mkl_norm==1)
00838 status = CPXlpopt (env, lp_cplex);
00839 else if (mkl_norm==2)
00840 status = CPXbaropt(env, lp_cplex);
00841 else
00842 {
00843 float64_t* beta=new float64_t[2*num_kernels+1];
00844 float64_t objval_old=1e-8;
00845 for (int32_t i=0; i<num_kernels; i++)
00846 beta[i]=old_beta[i];
00847 for (int32_t i=num_kernels; i<2*num_kernels+1; i++)
00848 beta[i]=0;
00849
00850 while (true)
00851 {
00852
00853
00854
00855 CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm), beta, num_kernels);
00856
00857 set_qnorm_constraints(beta, num_kernels);
00858
00859 status = CPXbaropt(env, lp_cplex);
00860 if ( status )
00861 SG_ERROR( "Failed to optimize Problem.\n");
00862
00863 int solstat=0;
00864 double objval=0;
00865 status=CPXsolution(env, lp_cplex, &solstat, &objval,
00866 (double*) beta, NULL, NULL, NULL);
00867
00868 if ( status )
00869 {
00870 CMath::display_vector(beta, num_kernels, "beta");
00871 SG_ERROR( "Failed to obtain solution.\n");
00872 }
00873
00874 CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm), beta, num_kernels);
00875
00876
00877 if ((1-abs(objval/objval_old) < 0.1*mkl_epsilon))
00878 break;
00879
00880 objval_old=objval;
00881
00882 inner_iters++;
00883 }
00884 delete[] beta;
00885 }
00886
00887 if ( status )
00888 SG_ERROR( "Failed to optimize Problem.\n");
00889
00890
00891 int32_t cur_numrows=(int32_t) CPXgetnumrows(env, lp_cplex);
00892 int32_t cur_numcols=(int32_t) CPXgetnumcols(env, lp_cplex);
00893 int32_t num_rows=cur_numrows;
00894 ASSERT(cur_numcols<=2*num_kernels+1);
00895
00896 float64_t* slack=new float64_t[cur_numrows];
00897 float64_t* pi=new float64_t[cur_numrows];
00898
00899
00900 int solstat=0;
00901 double objval=0;
00902
00903 if (mkl_norm==1)
00904 {
00905 status=CPXsolution(env, lp_cplex, &solstat, &objval,
00906 (double*) x, (double*) pi, (double*) slack, NULL);
00907 }
00908 else
00909 {
00910 status=CPXsolution(env, lp_cplex, &solstat, &objval,
00911 (double*) x, NULL, (double*) slack, NULL);
00912 }
00913
00914 int32_t solution_ok = (!status) ;
00915 if ( status )
00916 SG_ERROR( "Failed to obtain solution.\n");
00917
00918 int32_t num_active_rows=0 ;
00919 if (solution_ok)
00920 {
00921
00922 float64_t max_slack = -CMath::INFTY ;
00923 int32_t max_idx = -1 ;
00924 int32_t start_row = 1 ;
00925 if (C_mkl!=0.0)
00926 start_row+=2*(num_kernels-1);
00927
00928 for (int32_t i = start_row; i < cur_numrows; i++)
00929 {
00930 if (mkl_norm==1)
00931 {
00932 if ((pi[i]!=0))
00933 num_active_rows++ ;
00934 else
00935 {
00936 if (slack[i]>max_slack)
00937 {
00938 max_slack=slack[i] ;
00939 max_idx=i ;
00940 }
00941 }
00942 }
00943 else
00944 {
00945 if ((CMath::abs(slack[i])<1e-6))
00946 num_active_rows++ ;
00947 else
00948 {
00949 if (slack[i]>max_slack)
00950 {
00951 max_slack=slack[i] ;
00952 max_idx=i ;
00953 }
00954 }
00955 }
00956 }
00957
00958
00959 if ( (num_rows-start_row>CMath::max(100,2*num_active_rows)) && (max_idx!=-1))
00960 {
00961
00962 status = CPXdelrows (env, lp_cplex, max_idx, max_idx) ;
00963 if ( status )
00964 SG_ERROR( "Failed to remove an old row.\n");
00965 }
00966
00967
00968
00969 rho = -x[2*num_kernels] ;
00970 delete[] pi ;
00971 delete[] slack ;
00972
00973 }
00974 else
00975 {
00976
00977
00978 rho = 1 ;
00979 }
00980 }
00981 for (int32_t i=0; i<num_kernels; i++)
00982 new_beta[i]=x[i];
00983
00984 delete[] x;
00985 #else
00986 SG_ERROR("Cplex not enabled at compile time\n");
00987 #endif
00988 return rho;
00989 }
00990
00991 float64_t CMKL::compute_optimal_betas_via_glpk(float64_t* beta, const float64_t* old_beta,
00992 int num_kernels, const float64_t* sumw, float64_t suma, int32_t& inner_iters)
00993 {
00994 SG_DEBUG("MKL via GLPK\n");
00995
00996 if (mkl_norm!=1)
00997 SG_ERROR("MKL via GLPK works only for norm=1\n");
00998
00999 float64_t obj=1.0;
01000 #ifdef USE_GLPK
01001 int32_t NUMCOLS = 2*num_kernels + 1 ;
01002 if (!lp_initialized)
01003 {
01004 SG_INFO( "creating LP\n") ;
01005
01006
01007 lpx_add_cols(lp_glpk, NUMCOLS);
01008 for (int i=1; i<=2*num_kernels; i++)
01009 {
01010 lpx_set_obj_coef(lp_glpk, i, 0);
01011 lpx_set_col_bnds(lp_glpk, i, LPX_DB, 0, 1);
01012 }
01013 for (int i=num_kernels+1; i<=2*num_kernels; i++)
01014 {
01015 lpx_set_obj_coef(lp_glpk, i, C_mkl);
01016 }
01017 lpx_set_obj_coef(lp_glpk, NUMCOLS, 1);
01018 lpx_set_col_bnds(lp_glpk, NUMCOLS, LPX_FR, 0,0);
01019
01020
01021 int row_index = lpx_add_rows(lp_glpk, 1);
01022 int ind[num_kernels+2];
01023 float64_t val[num_kernels+2];
01024 for (int i=1; i<=num_kernels; i++)
01025 {
01026 ind[i] = i;
01027 val[i] = 1;
01028 }
01029 ind[num_kernels+1] = NUMCOLS;
01030 val[num_kernels+1] = 0;
01031 lpx_set_mat_row(lp_glpk, row_index, num_kernels, ind, val);
01032 lpx_set_row_bnds(lp_glpk, row_index, LPX_FX, 1, 1);
01033
01034 lp_initialized = true;
01035
01036 if (C_mkl!=0.0)
01037 {
01038 for (int32_t q=1; q<num_kernels; q++)
01039 {
01040 int mat_ind[4];
01041 float64_t mat_val[4];
01042 int mat_row_index = lpx_add_rows(lp_glpk, 2);
01043 mat_ind[1] = q;
01044 mat_val[1] = 1;
01045 mat_ind[2] = q+1;
01046 mat_val[2] = -1;
01047 mat_ind[3] = num_kernels+q;
01048 mat_val[3] = -1;
01049 lpx_set_mat_row(lp_glpk, mat_row_index, 3, mat_ind, mat_val);
01050 lpx_set_row_bnds(lp_glpk, mat_row_index, LPX_UP, 0, 0);
01051 mat_val[1] = -1;
01052 mat_val[2] = 1;
01053 lpx_set_mat_row(lp_glpk, mat_row_index+1, 3, mat_ind, mat_val);
01054 lpx_set_row_bnds(lp_glpk, mat_row_index+1, LPX_UP, 0, 0);
01055 }
01056 }
01057 }
01058
01059 int ind[num_kernels+2];
01060 float64_t val[num_kernels+2];
01061 int row_index = lpx_add_rows(lp_glpk, 1);
01062 for (int32_t i=1; i<=num_kernels; i++)
01063 {
01064 ind[i] = i;
01065 val[i] = -(sumw[i-1]-suma);
01066 }
01067 ind[num_kernels+1] = 2*num_kernels+1;
01068 val[num_kernels+1] = -1;
01069 lpx_set_mat_row(lp_glpk, row_index, num_kernels+1, ind, val);
01070 lpx_set_row_bnds(lp_glpk, row_index, LPX_UP, 0, 0);
01071
01072
01073 lpx_simplex(lp_glpk);
01074 bool res = check_lpx_status(lp_glpk);
01075 if (!res)
01076 SG_ERROR( "Failed to optimize Problem.\n");
01077
01078 int32_t cur_numrows = lpx_get_num_rows(lp_glpk);
01079 int32_t cur_numcols = lpx_get_num_cols(lp_glpk);
01080 int32_t num_rows=cur_numrows;
01081 ASSERT(cur_numcols<=2*num_kernels+1);
01082
01083 float64_t* col_primal = new float64_t[cur_numcols];
01084 float64_t* row_primal = new float64_t[cur_numrows];
01085 float64_t* row_dual = new float64_t[cur_numrows];
01086
01087 for (int i=0; i<cur_numrows; i++)
01088 {
01089 row_primal[i] = lpx_get_row_prim(lp_glpk, i+1);
01090 row_dual[i] = lpx_get_row_dual(lp_glpk, i+1);
01091 }
01092 for (int i=0; i<cur_numcols; i++)
01093 col_primal[i] = lpx_get_col_prim(lp_glpk, i+1);
01094
01095 obj = -col_primal[2*num_kernels];
01096
01097 for (int i=0; i<num_kernels; i++)
01098 beta[i] = col_primal[i];
01099
01100 int32_t num_active_rows=0;
01101 if(res)
01102 {
01103 float64_t max_slack = CMath::INFTY;
01104 int32_t max_idx = -1;
01105 int32_t start_row = 1;
01106 if (C_mkl!=0.0)
01107 start_row += 2*(num_kernels-1);
01108
01109 for (int32_t i= start_row; i<cur_numrows; i++)
01110 {
01111 if (row_dual[i]!=0)
01112 num_active_rows++;
01113 else
01114 {
01115 if (row_primal[i]<max_slack)
01116 {
01117 max_slack = row_primal[i];
01118 max_idx = i;
01119 }
01120 }
01121 }
01122
01123 if ((num_rows-start_row>CMath::max(100, 2*num_active_rows)) && max_idx!=-1)
01124 {
01125 int del_rows[2];
01126 del_rows[1] = max_idx+1;
01127 lpx_del_rows(lp_glpk, 1, del_rows);
01128 }
01129 }
01130
01131 delete[] row_dual;
01132 delete[] row_primal;
01133 delete[] col_primal;
01134 #else
01135 SG_ERROR("Glpk not enabled at compile time\n");
01136 #endif
01137
01138 return obj;
01139 }
01140
01141 void CMKL::compute_sum_beta(float64_t* sumw)
01142 {
01143 ASSERT(sumw);
01144 ASSERT(svm);
01145
01146 int32_t nsv=svm->get_num_support_vectors();
01147 int32_t num_kernels = kernel->get_num_subkernels();
01148 float64_t* beta = new float64_t[num_kernels];
01149 int32_t nweights=0;
01150 const float64_t* old_beta = kernel->get_subkernel_weights(nweights);
01151 ASSERT(nweights==num_kernels);
01152 ASSERT(old_beta);
01153
01154 for (int32_t i=0; i<num_kernels; i++)
01155 {
01156 beta[i]=0;
01157 sumw[i]=0;
01158 }
01159
01160 for (int32_t n=0; n<num_kernels; n++)
01161 {
01162 beta[n]=1.0;
01163 kernel->set_subkernel_weights(beta, num_kernels);
01164
01165 for (int32_t i=0; i<nsv; i++)
01166 {
01167 int32_t ii=svm->get_support_vector(i);
01168
01169 for (int32_t j=0; j<nsv; j++)
01170 {
01171 int32_t jj=svm->get_support_vector(j);
01172 sumw[n]+=0.5*svm->get_alpha(i)*svm->get_alpha(j)*kernel->kernel(ii,jj);
01173 }
01174 }
01175 beta[n]=0.0;
01176 }
01177
01178 mkl_iterations++;
01179 kernel->set_subkernel_weights( (float64_t*) old_beta, num_kernels);
01180 }
01181
01182
01183
01184 float64_t CMKL::compute_mkl_dual_objective()
01185 {
01186 int32_t n=get_num_support_vectors();
01187 float64_t mkl_obj=0;
01188
01189 if (labels && kernel && kernel->get_kernel_type() == K_COMBINED)
01190 {
01191 CKernel* kn = ((CCombinedKernel*)kernel)->get_first_kernel();
01192 while (kn)
01193 {
01194 float64_t sum=0;
01195 for (int32_t i=0; i<n; i++)
01196 {
01197 int32_t ii=get_support_vector(i);
01198
01199 for (int32_t j=0; j<n; j++)
01200 {
01201 int32_t jj=get_support_vector(j);
01202 sum+=get_alpha(i)*get_alpha(j)*kn->kernel(ii,jj);
01203 }
01204 }
01205
01206 if (mkl_norm==1.0)
01207 mkl_obj = CMath::max(mkl_obj, sum);
01208 else
01209 mkl_obj += CMath::pow(sum, mkl_norm/(mkl_norm-1));
01210
01211 SG_UNREF(kn);
01212 kn = ((CCombinedKernel*) kernel)->get_next_kernel();
01213 }
01214
01215 if (mkl_norm==1.0)
01216 mkl_obj=-0.5*mkl_obj;
01217 else
01218 mkl_obj= -0.5*CMath::pow(mkl_obj, (mkl_norm-1)/mkl_norm);
01219
01220 mkl_obj+=compute_sum_alpha();
01221 }
01222 else
01223 SG_ERROR( "cannot compute objective, labels or kernel not set\n");
01224
01225 return -mkl_obj;
01226 }
01227
01228 #ifdef USE_CPLEX
01229 void CMKL::set_qnorm_constraints(float64_t* beta, int32_t num_kernels)
01230 {
01231 ASSERT(num_kernels>0);
01232
01233 float64_t* grad_beta=new float64_t[num_kernels];
01234 float64_t* hess_beta=new float64_t[num_kernels+1];
01235 float64_t* lin_term=new float64_t[num_kernels+1];
01236 int* ind=new int[num_kernels+1];
01237
01238
01239 double const_term = 1-CMath::qsq(beta, num_kernels, mkl_norm);
01240
01241 ASSERT(CMath::fequal(const_term, 0.0));
01242
01243 for (int32_t i=0; i<num_kernels; i++)
01244 {
01245 grad_beta[i]=mkl_norm * pow(beta[i], mkl_norm-1);
01246 hess_beta[i]=0.5*mkl_norm*(mkl_norm-1) * pow(beta[i], mkl_norm-2);
01247 lin_term[i]=grad_beta[i] - 2*beta[i]*hess_beta[i];
01248 const_term+=grad_beta[i]*beta[i] - CMath::sq(beta[i])*hess_beta[i];
01249 ind[i]=i;
01250 }
01251 ind[num_kernels]=2*num_kernels;
01252 hess_beta[num_kernels]=0;
01253 lin_term[num_kernels]=0;
01254
01255 int status=0;
01256 int num=CPXgetnumqconstrs (env, lp_cplex);
01257
01258 if (num>0)
01259 {
01260 status = CPXdelqconstrs (env, lp_cplex, 0, 0);
01261 ASSERT(!status);
01262 }
01263
01264 status = CPXaddqconstr (env, lp_cplex, num_kernels+1, num_kernels+1, const_term, 'L', ind, lin_term,
01265 ind, ind, hess_beta, NULL);
01266 ASSERT(!status);
01267
01268
01269
01270
01271 delete[] grad_beta;
01272 delete[] hess_beta;
01273 delete[] lin_term;
01274 delete[] ind;
01275 }
01276 #endif // USE_CPLEX