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