SHOGUN v0.9.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 "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 // -- 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(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 // -- Initialize subkernel weights for Elasticnet MKL 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); //q-norm = 1 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 //Note that there is currently no safe way to properly resolve this 00285 //situation: the mkl object hands itself as a reference to the svm which 00286 //in turn increases the reference count of mkl and when done decreases 00287 //the count. Thus we have to protect the mkl object from deletion by 00288 //ref()'ing it before we set the callback function and should also 00289 //unref() it afterwards. However, when the reference count was zero 00290 //before this unref() might actually try to destroy this (crash ahead) 00291 //but if we don't actually unref() the object we might leak memory... 00292 //So as a workaround we only unref when the reference count was >1 00293 //before. 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 // -- Direct update of subkernel weights for ElasticnetMKL 00395 // Note that ElasticnetMKL solves a different optimization 00396 // problem from the rest of the solver types 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; // fraction of root mean squared deviation 00427 float64_t obj; 00428 float64_t Z; 00429 float64_t preR; 00430 int32_t p; 00431 int32_t nofKernelsGood; 00432 00433 // --- optimal beta 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 // --- normalize 00453 CMath::scale_vector(1.0/Z, beta, num_kernels); 00454 00455 // --- regularize & renormalize 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 // --- copy beta into beta_local 00497 for( p=0; p<num_kernels; ++p ) 00498 beta_local[p] = beta[p]; 00499 00500 // --- elastic-net transform 00501 elasticnet_transform(beta, ent_lambda, num_kernels); 00502 00503 // --- objective 00504 obj = -suma; 00505 for (p=0; p<num_kernels; ++p ) 00506 { 00507 //obj += sumw[p] * old_beta[p]*old_beta[p] / beta[p]; 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 // assumes that all constraints are satisfied 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 // Compute Elastic-net dual 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 // SG_PRINT("nm[%d]=%f\n",k,nm[k]); 00571 k++; 00572 00573 00574 SG_UNREF(kn); 00575 kn = ((CCombinedKernel*) kernel)->get_next_kernel(); 00576 } 00577 // initial delta 00578 del = del/CMath::sqrt(2*(1-ent_lambda)); 00579 00580 // Newton's method to optimize delta 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 // SG_PRINT("[%d] fval=%f gg=%f del=%f\n", k, ff, gg, del); 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; // fraction of root mean squared deviation 00622 float64_t obj; 00623 float64_t Z; 00624 float64_t preR; 00625 int32_t p; 00626 int32_t nofKernelsGood; 00627 00628 // --- optimal beta 00629 nofKernelsGood = num_kernels; 00630 for( p=0; p<num_kernels; ++p ) 00631 { 00632 //SG_PRINT( "MKL-direct: sumw[%3d] = %e ( oldbeta = %e )\n", p, sumw[p], old_beta[p] ); 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 // --- normalize 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 // --- regularize & renormalize 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 // --- objective 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 // --- check beta 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 // --- compute gamma 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 // should better recompute sumw[] !!! 00769 } 00770 else 00771 { 00772 ASSERT( sumw[p] >= 0 ); 00773 //gamma += CMath::pow( sumw[p] * beta[p]*beta[p], r ); 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 // should better recompute sumw[] !!! 00783 gamma = epsGamma; 00784 } 00785 ASSERT( gamma >= epsGamma ); 00786 //gamma = -gamma; 00787 00788 // --- compute objective 00789 obj = 0.0; 00790 for( p=0; p<num_kernels; ++p ) 00791 { 00792 obj += beta[p] * sumw[p]; 00793 //obj += gamma/mkl_norm * CMath::pow( beta[p], mkl_norm ); 00794 } 00795 if( !( obj >= 0.0 ) ) 00796 SG_WARNING( "negative objective: %e. ", obj ); 00797 //SG_PRINT( "OBJ = %e. \n", obj ); 00798 00799 00800 // === perform Newton steps 00801 for (i = 0; i < nofNewtonSteps; ++i ) 00802 { 00803 // --- compute Newton direction (Hessian is diagonal) 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 //const float halfw2p = ( sumw[p] >= 0.0 ) ? sumw[p] : 0.0; 00810 //const float64_t t1 = halfw2p*beta[p] - mkl_norm*gamma*CMath::pow(beta[p],mkl_norm); 00811 //const float64_t t2 = 2.0*halfw2p + gqq1*CMath::pow(beta[p],mkl_norm-1.0); 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 // newtStep += newtDir[p] * grad[p]; 00821 ASSERT( newtDir[p] == newtDir[p] ); 00822 //SG_PRINT( "newtDir[%d] = %6.3f = %e / %e \n", p, newtDir[p], t1, t2 ); 00823 } 00824 //CMath::display_vector( newtDir, num_kernels, "newton direction " ); 00825 //SG_PRINT( "Newton step size = %e\n", Z ); 00826 00827 // --- line search 00828 stepSize = 1.0; 00829 while( stepSize >= epsStep ) 00830 { 00831 // --- perform Newton step 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 // --- normalize new beta (wrt p-norm) 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 //SG_WARNING( "beta[%d] = %e; set to 1. ", p, beta[p] ); 00860 newtBeta[p] = 1.0; 00861 } 00862 ASSERT( 0.0 <= newtBeta[p] && newtBeta[p] <= 1.0 ); 00863 } 00864 00865 // --- objective increased? 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 //SG_PRINT( "step = %.8f: obj = %e -> %e. \n", stepSize, obj, newtObj ); 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 // === return new objective 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]; /* calling external lib */ 00914 double lb[NUMCOLS]; /* calling external lib */ 00915 double ub[NUMCOLS]; /* calling external lib */ 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 // add constraint sum(w)=1; 00939 SG_INFO( "adding the first row\n"); 00940 int initial_rmatbeg[1]; /* calling external lib */ 00941 int initial_rmatind[num_kernels+1]; /* calling external lib */ 00942 double initial_rmatval[num_kernels+1]; /* calling ext lib */ 00943 double initial_rhs[1]; /* calling external lib */ 00944 char initial_sense[1]; 00945 00946 // 1-norm MKL 00947 if (mkl_norm==1) 00948 { 00949 initial_rmatbeg[0] = 0; 00950 initial_rhs[0]=1 ; // rhs=1 ; 00951 initial_sense[0]='E' ; // equality 00952 00953 // sparse matrix 00954 for (int32_t i=0; i<num_kernels; i++) 00955 { 00956 initial_rmatind[i]=i ; //index of non-zero element 00957 initial_rmatval[i]=1 ; //value of ... 00958 } 00959 initial_rmatind[num_kernels]=2*num_kernels ; //number of non-zero elements 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 // 2 and q-norm MKL 00968 { 00969 initial_rmatbeg[0] = 0; 00970 initial_rhs[0]=0 ; // rhs=1 ; 00971 initial_sense[0]='L' ; // <= (inequality) 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 // add constraint w[i]-w[i+1]<s[i]; 01007 // add constraint w[i+1]-w[i]<s[i]; 01008 int rmatbeg[1]; /* calling external lib */ 01009 int rmatind[3]; /* calling external lib */ 01010 double rmatval[3]; /* calling external lib */ 01011 double rhs[1]; /* calling external lib */ 01012 char sense[1]; 01013 01014 rmatbeg[0] = 0; 01015 rhs[0]=0 ; // rhs=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 ; // rhs=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 { // add the new row 01048 //SG_INFO( "add the new row\n") ; 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) // optimize 1 norm MKL 01087 status = CPXlpopt (env, lp_cplex); 01088 else if (mkl_norm==2) // optimize 2-norm MKL 01089 status = CPXbaropt(env, lp_cplex); 01090 else // q-norm MKL 01091 { 01092 float64_t* beta=new float64_t[2*num_kernels+1]; 01093 float64_t objval_old=1e-8; //some value to cause the loop to not terminate yet 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 //int rows=CPXgetnumrows(env, lp_cplex); 01102 //int cols=CPXgetnumcols(env, lp_cplex); 01103 //SG_PRINT("rows:%d, cols:%d (kernel:%d)\n", rows, cols, num_kernels); 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 //SG_PRINT("[%d] %f (%f)\n", inner_iters, objval, objval_old); 01126 if ((1-abs(objval/objval_old) < 0.1*mkl_epsilon)) // && (inner_iters>2)) 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 // obtain solution 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 /* calling external lib */ 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 /* 1 norm mkl */ 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++) // skip first 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 // 2-norm or general q-norm 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 // have at most max(100,num_active_rows*2) rows, if not, remove one 01208 if ( (num_rows-start_row>CMath::max(100,2*num_active_rows)) && (max_idx!=-1)) 01209 { 01210 //SG_INFO( "-%i(%i,%i)",max_idx,start_row,num_rows) ; 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 //CMath::display_vector(x, num_kernels, "beta"); 01217 01218 rho = -x[2*num_kernels] ; 01219 delete[] pi ; 01220 delete[] slack ; 01221 01222 } 01223 else 01224 { 01225 /* then something is wrong and we rather 01226 stop sooner than later */ 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 //set obj function, note glpk indexing is 1-based 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); //unbound 01268 01269 //add first row. sum[w]=1 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 //optimize 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 // assumes that all constraints are satisfied 01433 float64_t CMKL::compute_mkl_dual_objective() 01434 { 01435 if (get_solver_type()==ST_ELASTICNET) 01436 { 01437 // -- Compute ElasticnetMKL dual objective 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 //CMath::display_vector(beta, num_kernels, "beta"); 01494 double const_term = 1-CMath::qsq(beta, num_kernels, mkl_norm); 01495 01496 //SG_PRINT("const=%f\n", const_term); 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 //CPXwriteprob (env, lp_cplex, "prob.lp", NULL); 01525 //CPXqpwrite (env, lp_cplex, "prob.qp"); 01526 01527 delete[] grad_beta; 01528 delete[] hess_beta; 01529 delete[] lin_term; 01530 delete[] ind; 01531 } 01532 #endif // USE_CPLEX