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