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