SHOGUN
v1.1.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2004-2009 Soeren Sonnenburg, Gunnar Raetsch 00008 * Alexander Zien, Marius Kloft, Chen Guohua 00009 * Copyright (C) 2009 Fraunhofer Institute FIRST and Max-Planck-Society 00010 * Copyright (C) 2010 Ryota Tomioka (University of Tokyo) 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 // -- Delete beta_local for ElasticnetMKL 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; // calling external lib 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 // select dual simplex based optimization 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); /* Problem is minimization */ 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 // -- Initialize subkernel weights for Elasticnet MKL 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); //q-norm = 1 00254 } 00255 00256 kernel->set_subkernel_weights(SGVector<float64_t>(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 //Note that there is currently no safe way to properly resolve this 00297 //situation: the mkl object hands itself as a reference to the svm which 00298 //in turn increases the reference count of mkl and when done decreases 00299 //the count. Thus we have to protect the mkl object from deletion by 00300 //ref()'ing it before we set the callback function and should also 00301 //unref() it afterwards. However, when the reference count was zero 00302 //before this unref() might actually try to destroy this (crash ahead) 00303 //but if we don't actually unref() the object we might leak memory... 00304 //So as a workaround we only unref when the reference count was >1 00305 //before. 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 // -- Direct update of subkernel weights for ElasticnetMKL 00436 // Note that ElasticnetMKL solves a different optimization 00437 // problem from the rest of the solver types 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(SGVector<float64_t>(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; // fraction of root mean squared deviation 00468 float64_t obj; 00469 float64_t Z; 00470 float64_t preR; 00471 int32_t p; 00472 int32_t nofKernelsGood; 00473 00474 // --- optimal beta 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 // --- normalize 00494 CMath::scale_vector(1.0/Z, beta, num_kernels); 00495 00496 // --- regularize & renormalize 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 // --- copy beta into beta_local 00538 for( p=0; p<num_kernels; ++p ) 00539 beta_local[p] = beta[p]; 00540 00541 // --- elastic-net transform 00542 elasticnet_transform(beta, ent_lambda, num_kernels); 00543 00544 // --- objective 00545 obj = -suma; 00546 for (p=0; p<num_kernels; ++p ) 00547 { 00548 //obj += sumw[p] * old_beta[p]*old_beta[p] / beta[p]; 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 // assumes that all constraints are satisfied 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 // Compute Elastic-net dual 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 // SG_PRINT("nm[%d]=%f\n",k,nm[k]); 00612 k++; 00613 00614 00615 SG_UNREF(kn); 00616 kn = ((CCombinedKernel*) kernel)->get_next_kernel(); 00617 } 00618 // initial delta 00619 del = del/CMath::sqrt(2*(1-ent_lambda)); 00620 00621 // Newton's method to optimize delta 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 // SG_PRINT("[%d] fval=%f gg=%f del=%f\n", k, ff, gg, del); 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 // --- optimal beta 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 // --- objective 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; // fraction of root mean squared deviation 00698 float64_t obj; 00699 float64_t Z; 00700 float64_t preR; 00701 int32_t p; 00702 int32_t nofKernelsGood; 00703 00704 // --- optimal beta 00705 nofKernelsGood = num_kernels; 00706 for( p=0; p<num_kernels; ++p ) 00707 { 00708 //SG_PRINT( "MKL-direct: sumw[%3d] = %e ( oldbeta = %e )\n", p, sumw[p], old_beta[p] ); 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 // --- normalize 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 // --- regularize & renormalize 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 // --- objective 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 //float64_t newtStep; 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 // --- check beta 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 // --- compute gamma 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 // should better recompute sumw[] !!! 00845 } 00846 else 00847 { 00848 ASSERT( sumw[p] >= 0 ); 00849 //gamma += CMath::pow( sumw[p] * beta[p]*beta[p], r ); 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 // should better recompute sumw[] !!! 00859 gamma = epsGamma; 00860 } 00861 ASSERT( gamma >= epsGamma ); 00862 //gamma = -gamma; 00863 00864 // --- compute objective 00865 obj = 0.0; 00866 for( p=0; p<num_kernels; ++p ) 00867 { 00868 obj += beta[p] * sumw[p]; 00869 //obj += gamma/mkl_norm * CMath::pow( beta[p], mkl_norm ); 00870 } 00871 if( !( obj >= 0.0 ) ) 00872 SG_WARNING( "negative objective: %e. ", obj ); 00873 //SG_PRINT( "OBJ = %e. \n", obj ); 00874 00875 00876 // === perform Newton steps 00877 for (i = 0; i < nofNewtonSteps; ++i ) 00878 { 00879 // --- compute Newton direction (Hessian is diagonal) 00880 const float64_t gqq1 = mkl_norm * (mkl_norm-1.0) * gamma; 00881 // newtStep = 0.0; 00882 for( p=0; p<num_kernels; ++p ) 00883 { 00884 ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 ); 00885 //const float halfw2p = ( sumw[p] >= 0.0 ) ? sumw[p] : 0.0; 00886 //const float64_t t1 = halfw2p*beta[p] - mkl_norm*gamma*CMath::pow(beta[p],mkl_norm); 00887 //const float64_t t2 = 2.0*halfw2p + gqq1*CMath::pow(beta[p],mkl_norm-1.0); 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 // newtStep += newtDir[p] * grad[p]; 00897 ASSERT( newtDir[p] == newtDir[p] ); 00898 //SG_PRINT( "newtDir[%d] = %6.3f = %e / %e \n", p, newtDir[p], t1, t2 ); 00899 } 00900 //CMath::display_vector( newtDir, num_kernels, "newton direction " ); 00901 //SG_PRINT( "Newton step size = %e\n", Z ); 00902 00903 // --- line search 00904 stepSize = 1.0; 00905 while( stepSize >= epsStep ) 00906 { 00907 // --- perform Newton step 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 // --- normalize new beta (wrt p-norm) 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 //SG_WARNING( "beta[%d] = %e; set to 1. ", p, beta[p] ); 00936 newtBeta[p] = 1.0; 00937 } 00938 ASSERT( 0.0 <= newtBeta[p] && newtBeta[p] <= 1.0 ); 00939 } 00940 00941 // --- objective increased? 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 //SG_PRINT( "step = %.8f: obj = %e -> %e. \n", stepSize, obj, newtObj ); 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 // === return new objective 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]; /* calling external lib */ 00990 double lb[NUMCOLS]; /* calling external lib */ 00991 double ub[NUMCOLS]; /* calling external lib */ 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 // add constraint sum(w)=1; 01015 SG_INFO( "adding the first row\n"); 01016 int initial_rmatbeg[1]; /* calling external lib */ 01017 int initial_rmatind[num_kernels+1]; /* calling external lib */ 01018 double initial_rmatval[num_kernels+1]; /* calling ext lib */ 01019 double initial_rhs[1]; /* calling external lib */ 01020 char initial_sense[1]; 01021 01022 // 1-norm MKL 01023 if (mkl_norm==1) 01024 { 01025 initial_rmatbeg[0] = 0; 01026 initial_rhs[0]=1 ; // rhs=1 ; 01027 initial_sense[0]='E' ; // equality 01028 01029 // sparse matrix 01030 for (int32_t i=0; i<num_kernels; i++) 01031 { 01032 initial_rmatind[i]=i ; //index of non-zero element 01033 initial_rmatval[i]=1 ; //value of ... 01034 } 01035 initial_rmatind[num_kernels]=2*num_kernels ; //number of non-zero elements 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 // 2 and q-norm MKL 01044 { 01045 initial_rmatbeg[0] = 0; 01046 initial_rhs[0]=0 ; // rhs=1 ; 01047 initial_sense[0]='L' ; // <= (inequality) 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 // add constraint w[i]-w[i+1]<s[i]; 01083 // add constraint w[i+1]-w[i]<s[i]; 01084 int rmatbeg[1]; /* calling external lib */ 01085 int rmatind[3]; /* calling external lib */ 01086 double rmatval[3]; /* calling external lib */ 01087 double rhs[1]; /* calling external lib */ 01088 char sense[1]; 01089 01090 rmatbeg[0] = 0; 01091 rhs[0]=0 ; // rhs=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 ; // rhs=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 { // add the new row 01124 //SG_INFO( "add the new row\n") ; 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) // optimize 1 norm MKL 01163 status = CPXlpopt (env, lp_cplex); 01164 else if (mkl_norm==2) // optimize 2-norm MKL 01165 status = CPXbaropt(env, lp_cplex); 01166 else // q-norm MKL 01167 { 01168 float64_t* beta=SG_MALLOC(float64_t, 2*num_kernels+1); 01169 float64_t objval_old=1e-8; //some value to cause the loop to not terminate yet 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 //int rows=CPXgetnumrows(env, lp_cplex); 01178 //int cols=CPXgetnumcols(env, lp_cplex); 01179 //SG_PRINT("rows:%d, cols:%d (kernel:%d)\n", rows, cols, num_kernels); 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 //SG_PRINT("[%d] %f (%f)\n", inner_iters, objval, objval_old); 01202 if ((1-abs(objval/objval_old) < 0.1*mkl_epsilon)) // && (inner_iters>2)) 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 // obtain solution 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 /* calling external lib */ 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 /* 1 norm mkl */ 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++) // skip first 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 // 2-norm or general q-norm 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 // have at most max(100,num_active_rows*2) rows, if not, remove one 01284 if ( (num_rows-start_row>CMath::max(100,2*num_active_rows)) && (max_idx!=-1)) 01285 { 01286 //SG_INFO( "-%i(%i,%i)",max_idx,start_row,num_rows) ; 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 //CMath::display_vector(x, num_kernels, "beta"); 01293 01294 rho = -x[2*num_kernels] ; 01295 SG_FREE(pi); 01296 SG_FREE(slack); 01297 01298 } 01299 else 01300 { 01301 /* then something is wrong and we rather 01302 stop sooner than later */ 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 //set obj function, note glpk indexing is 1-based 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); //unbound 01344 01345 //add first row. sum[w]=1 01346 int row_index = lpx_add_rows(lp_glpk, 1); 01347 int* ind = SG_MALLOC(int, num_kernels+2); 01348 float64_t* val = SG_MALLOC(float64_t, 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 SG_FREE(val); 01359 SG_FREE(ind); 01360 01361 lp_initialized = true; 01362 01363 if (C_mkl!=0.0) 01364 { 01365 for (int32_t q=1; q<num_kernels; q++) 01366 { 01367 int mat_ind[4]; 01368 float64_t mat_val[4]; 01369 int mat_row_index = lpx_add_rows(lp_glpk, 2); 01370 mat_ind[1] = q; 01371 mat_val[1] = 1; 01372 mat_ind[2] = q+1; 01373 mat_val[2] = -1; 01374 mat_ind[3] = num_kernels+q; 01375 mat_val[3] = -1; 01376 lpx_set_mat_row(lp_glpk, mat_row_index, 3, mat_ind, mat_val); 01377 lpx_set_row_bnds(lp_glpk, mat_row_index, LPX_UP, 0, 0); 01378 mat_val[1] = -1; 01379 mat_val[2] = 1; 01380 lpx_set_mat_row(lp_glpk, mat_row_index+1, 3, mat_ind, mat_val); 01381 lpx_set_row_bnds(lp_glpk, mat_row_index+1, LPX_UP, 0, 0); 01382 } 01383 } 01384 } 01385 01386 int* ind=SG_MALLOC(int,num_kernels+2); 01387 float64_t* val=SG_MALLOC(float64_t, num_kernels+2); 01388 int row_index = lpx_add_rows(lp_glpk, 1); 01389 for (int32_t i=1; i<=num_kernels; i++) 01390 { 01391 ind[i] = i; 01392 val[i] = -(sumw[i-1]-suma); 01393 } 01394 ind[num_kernels+1] = 2*num_kernels+1; 01395 val[num_kernels+1] = -1; 01396 lpx_set_mat_row(lp_glpk, row_index, num_kernels+1, ind, val); 01397 lpx_set_row_bnds(lp_glpk, row_index, LPX_UP, 0, 0); 01398 SG_FREE(ind); 01399 SG_FREE(val); 01400 01401 //optimize 01402 lpx_simplex(lp_glpk); 01403 bool res = check_lpx_status(lp_glpk); 01404 if (!res) 01405 SG_ERROR( "Failed to optimize Problem.\n"); 01406 01407 int32_t cur_numrows = lpx_get_num_rows(lp_glpk); 01408 int32_t cur_numcols = lpx_get_num_cols(lp_glpk); 01409 int32_t num_rows=cur_numrows; 01410 ASSERT(cur_numcols<=2*num_kernels+1); 01411 01412 float64_t* col_primal = SG_MALLOC(float64_t, cur_numcols); 01413 float64_t* row_primal = SG_MALLOC(float64_t, cur_numrows); 01414 float64_t* row_dual = SG_MALLOC(float64_t, cur_numrows); 01415 01416 for (int i=0; i<cur_numrows; i++) 01417 { 01418 row_primal[i] = lpx_get_row_prim(lp_glpk, i+1); 01419 row_dual[i] = lpx_get_row_dual(lp_glpk, i+1); 01420 } 01421 for (int i=0; i<cur_numcols; i++) 01422 col_primal[i] = lpx_get_col_prim(lp_glpk, i+1); 01423 01424 obj = -col_primal[2*num_kernels]; 01425 01426 for (int i=0; i<num_kernels; i++) 01427 beta[i] = col_primal[i]; 01428 01429 int32_t num_active_rows=0; 01430 if(res) 01431 { 01432 float64_t max_slack = CMath::INFTY; 01433 int32_t max_idx = -1; 01434 int32_t start_row = 1; 01435 if (C_mkl!=0.0) 01436 start_row += 2*(num_kernels-1); 01437 01438 for (int32_t i= start_row; i<cur_numrows; i++) 01439 { 01440 if (row_dual[i]!=0) 01441 num_active_rows++; 01442 else 01443 { 01444 if (row_primal[i]<max_slack) 01445 { 01446 max_slack = row_primal[i]; 01447 max_idx = i; 01448 } 01449 } 01450 } 01451 01452 if ((num_rows-start_row>CMath::max(100, 2*num_active_rows)) && max_idx!=-1) 01453 { 01454 int del_rows[2]; 01455 del_rows[1] = max_idx+1; 01456 lpx_del_rows(lp_glpk, 1, del_rows); 01457 } 01458 } 01459 01460 SG_FREE(row_dual); 01461 SG_FREE(row_primal); 01462 SG_FREE(col_primal); 01463 #else 01464 SG_ERROR("Glpk not enabled at compile time\n"); 01465 #endif 01466 01467 return obj; 01468 } 01469 01470 void CMKL::compute_sum_beta(float64_t* sumw) 01471 { 01472 ASSERT(sumw); 01473 ASSERT(svm); 01474 01475 int32_t nsv=svm->get_num_support_vectors(); 01476 int32_t num_kernels = kernel->get_num_subkernels(); 01477 float64_t* beta = SG_MALLOC(float64_t, num_kernels); 01478 int32_t nweights=0; 01479 const float64_t* old_beta = kernel->get_subkernel_weights(nweights); 01480 ASSERT(nweights==num_kernels); 01481 ASSERT(old_beta); 01482 01483 for (int32_t i=0; i<num_kernels; i++) 01484 { 01485 beta[i]=0; 01486 sumw[i]=0; 01487 } 01488 01489 for (int32_t n=0; n<num_kernels; n++) 01490 { 01491 beta[n]=1.0; 01492 kernel->set_subkernel_weights(SGVector<float64_t>(beta, num_kernels)); 01493 01494 for (int32_t i=0; i<nsv; i++) 01495 { 01496 int32_t ii=svm->get_support_vector(i); 01497 01498 for (int32_t j=0; j<nsv; j++) 01499 { 01500 int32_t jj=svm->get_support_vector(j); 01501 sumw[n]+=0.5*svm->get_alpha(i)*svm->get_alpha(j)*kernel->kernel(ii,jj); 01502 } 01503 } 01504 beta[n]=0.0; 01505 } 01506 01507 mkl_iterations++; 01508 kernel->set_subkernel_weights(SGVector<float64_t>( (float64_t*) old_beta, num_kernels)); 01509 } 01510 01511 01512 // assumes that all constraints are satisfied 01513 float64_t CMKL::compute_mkl_dual_objective() 01514 { 01515 if (get_solver_type()==ST_ELASTICNET) 01516 { 01517 // -- Compute ElasticnetMKL dual objective 01518 return compute_elasticnet_dual_objective(); 01519 } 01520 01521 int32_t n=get_num_support_vectors(); 01522 float64_t mkl_obj=0; 01523 01524 if (labels && kernel && kernel->get_kernel_type() == K_COMBINED) 01525 { 01526 CKernel* kn = ((CCombinedKernel*)kernel)->get_first_kernel(); 01527 while (kn) 01528 { 01529 float64_t sum=0; 01530 for (int32_t i=0; i<n; i++) 01531 { 01532 int32_t ii=get_support_vector(i); 01533 01534 for (int32_t j=0; j<n; j++) 01535 { 01536 int32_t jj=get_support_vector(j); 01537 sum+=get_alpha(i)*get_alpha(j)*kn->kernel(ii,jj); 01538 } 01539 } 01540 01541 if (mkl_norm==1.0) 01542 mkl_obj = CMath::max(mkl_obj, sum); 01543 else 01544 mkl_obj += CMath::pow(sum, mkl_norm/(mkl_norm-1)); 01545 01546 SG_UNREF(kn); 01547 kn = ((CCombinedKernel*) kernel)->get_next_kernel(); 01548 } 01549 01550 if (mkl_norm==1.0) 01551 mkl_obj=-0.5*mkl_obj; 01552 else 01553 mkl_obj= -0.5*CMath::pow(mkl_obj, (mkl_norm-1)/mkl_norm); 01554 01555 mkl_obj+=compute_sum_alpha(); 01556 } 01557 else 01558 SG_ERROR( "cannot compute objective, labels or kernel not set\n"); 01559 01560 return -mkl_obj; 01561 } 01562 01563 #ifdef USE_CPLEX 01564 void CMKL::set_qnorm_constraints(float64_t* beta, int32_t num_kernels) 01565 { 01566 ASSERT(num_kernels>0); 01567 01568 float64_t* grad_beta=SG_MALLOC(float64_t, num_kernels); 01569 float64_t* hess_beta=SG_MALLOC(float64_t, num_kernels+1); 01570 float64_t* lin_term=SG_MALLOC(float64_t, num_kernels+1); 01571 int* ind=SG_MALLOC(int, num_kernels+1); 01572 01573 //CMath::display_vector(beta, num_kernels, "beta"); 01574 double const_term = 1-CMath::qsq(beta, num_kernels, mkl_norm); 01575 01576 //SG_PRINT("const=%f\n", const_term); 01577 ASSERT(CMath::fequal(const_term, 0.0)); 01578 01579 for (int32_t i=0; i<num_kernels; i++) 01580 { 01581 grad_beta[i]=mkl_norm * pow(beta[i], mkl_norm-1); 01582 hess_beta[i]=0.5*mkl_norm*(mkl_norm-1) * pow(beta[i], mkl_norm-2); 01583 lin_term[i]=grad_beta[i] - 2*beta[i]*hess_beta[i]; 01584 const_term+=grad_beta[i]*beta[i] - CMath::sq(beta[i])*hess_beta[i]; 01585 ind[i]=i; 01586 } 01587 ind[num_kernels]=2*num_kernels; 01588 hess_beta[num_kernels]=0; 01589 lin_term[num_kernels]=0; 01590 01591 int status=0; 01592 int num=CPXgetnumqconstrs (env, lp_cplex); 01593 01594 if (num>0) 01595 { 01596 status = CPXdelqconstrs (env, lp_cplex, 0, 0); 01597 ASSERT(!status); 01598 } 01599 01600 status = CPXaddqconstr (env, lp_cplex, num_kernels+1, num_kernels+1, const_term, 'L', ind, lin_term, 01601 ind, ind, hess_beta, NULL); 01602 ASSERT(!status); 01603 01604 //CPXwriteprob (env, lp_cplex, "prob.lp", NULL); 01605 //CPXqpwrite (env, lp_cplex, "prob.qp"); 01606 01607 SG_FREE(grad_beta); 01608 SG_FREE(hess_beta); 01609 SG_FREE(lin_term); 01610 SG_FREE(ind); 01611 } 01612 #endif // USE_CPLEX