MKL.cpp

Go to the documentation of this file.
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(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(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[num_kernels+2];
01348         float64_t val[num_kernels+2];
01349         for (int i=1; i<=num_kernels; i++)
01350         {
01351             ind[i] = i;
01352             val[i] = 1;
01353         }
01354         ind[num_kernels+1] = NUMCOLS;
01355         val[num_kernels+1] = 0;
01356         lpx_set_mat_row(lp_glpk, row_index, num_kernels, ind, val);
01357         lpx_set_row_bnds(lp_glpk, row_index, LPX_FX, 1, 1);
01358 
01359         lp_initialized = true;
01360 
01361         if (C_mkl!=0.0)
01362         {
01363             for (int32_t q=1; q<num_kernels; q++)
01364             {
01365                 int mat_ind[4];
01366                 float64_t mat_val[4];
01367                 int mat_row_index = lpx_add_rows(lp_glpk, 2);
01368                 mat_ind[1] = q;
01369                 mat_val[1] = 1;
01370                 mat_ind[2] = q+1;
01371                 mat_val[2] = -1;
01372                 mat_ind[3] = num_kernels+q;
01373                 mat_val[3] = -1;
01374                 lpx_set_mat_row(lp_glpk, mat_row_index, 3, mat_ind, mat_val);
01375                 lpx_set_row_bnds(lp_glpk, mat_row_index, LPX_UP, 0, 0);
01376                 mat_val[1] = -1;
01377                 mat_val[2] = 1;
01378                 lpx_set_mat_row(lp_glpk, mat_row_index+1, 3, mat_ind, mat_val);
01379                 lpx_set_row_bnds(lp_glpk, mat_row_index+1, LPX_UP, 0, 0);
01380             }
01381         }
01382     }
01383 
01384     int ind[num_kernels+2];
01385     float64_t val[num_kernels+2];
01386     int row_index = lpx_add_rows(lp_glpk, 1);
01387     for (int32_t i=1; i<=num_kernels; i++)
01388     {
01389         ind[i] = i;
01390         val[i] = -(sumw[i-1]-suma);
01391     }
01392     ind[num_kernels+1] = 2*num_kernels+1;
01393     val[num_kernels+1] = -1;
01394     lpx_set_mat_row(lp_glpk, row_index, num_kernels+1, ind, val);
01395     lpx_set_row_bnds(lp_glpk, row_index, LPX_UP, 0, 0);
01396 
01397     //optimize
01398     lpx_simplex(lp_glpk);
01399     bool res = check_lpx_status(lp_glpk);
01400     if (!res)
01401         SG_ERROR( "Failed to optimize Problem.\n");
01402 
01403     int32_t cur_numrows = lpx_get_num_rows(lp_glpk);
01404     int32_t cur_numcols = lpx_get_num_cols(lp_glpk);
01405     int32_t num_rows=cur_numrows;
01406     ASSERT(cur_numcols<=2*num_kernels+1);
01407 
01408     float64_t* col_primal = SG_MALLOC(float64_t, cur_numcols);
01409     float64_t* row_primal = SG_MALLOC(float64_t, cur_numrows);
01410     float64_t* row_dual = SG_MALLOC(float64_t, cur_numrows);
01411 
01412     for (int i=0; i<cur_numrows; i++)
01413     {
01414         row_primal[i] = lpx_get_row_prim(lp_glpk, i+1);
01415         row_dual[i] = lpx_get_row_dual(lp_glpk, i+1);
01416     }
01417     for (int i=0; i<cur_numcols; i++)
01418         col_primal[i] = lpx_get_col_prim(lp_glpk, i+1);
01419 
01420     obj = -col_primal[2*num_kernels];
01421 
01422     for (int i=0; i<num_kernels; i++)
01423         beta[i] = col_primal[i];
01424 
01425     int32_t num_active_rows=0;
01426     if(res)
01427     {
01428         float64_t max_slack = CMath::INFTY;
01429         int32_t max_idx = -1;
01430         int32_t start_row = 1;
01431         if (C_mkl!=0.0)
01432             start_row += 2*(num_kernels-1);
01433 
01434         for (int32_t i= start_row; i<cur_numrows; i++)
01435         {
01436             if (row_dual[i]!=0)
01437                 num_active_rows++;
01438             else
01439             {
01440                 if (row_primal[i]<max_slack)
01441                 {
01442                     max_slack = row_primal[i];
01443                     max_idx = i;
01444                 }
01445             }
01446         }
01447 
01448         if ((num_rows-start_row>CMath::max(100, 2*num_active_rows)) && max_idx!=-1)
01449         {
01450             int del_rows[2];
01451             del_rows[1] = max_idx+1;
01452             lpx_del_rows(lp_glpk, 1, del_rows);
01453         }
01454     }
01455 
01456     SG_FREE(row_dual);
01457     SG_FREE(row_primal);
01458     SG_FREE(col_primal);
01459 #else
01460     SG_ERROR("Glpk not enabled at compile time\n");
01461 #endif
01462 
01463     return obj;
01464 }
01465 
01466 void CMKL::compute_sum_beta(float64_t* sumw)
01467 {
01468     ASSERT(sumw);
01469     ASSERT(svm);
01470 
01471     int32_t nsv=svm->get_num_support_vectors();
01472     int32_t num_kernels = kernel->get_num_subkernels();
01473     float64_t* beta = SG_MALLOC(float64_t, num_kernels);
01474     int32_t nweights=0;
01475     const float64_t* old_beta = kernel->get_subkernel_weights(nweights);
01476     ASSERT(nweights==num_kernels);
01477     ASSERT(old_beta);
01478 
01479     for (int32_t i=0; i<num_kernels; i++)
01480     {
01481         beta[i]=0;
01482         sumw[i]=0;
01483     }
01484 
01485     for (int32_t n=0; n<num_kernels; n++)
01486     {
01487         beta[n]=1.0;
01488         kernel->set_subkernel_weights(beta, num_kernels);
01489 
01490         for (int32_t i=0; i<nsv; i++)
01491         {
01492             int32_t ii=svm->get_support_vector(i);
01493 
01494             for (int32_t j=0; j<nsv; j++)
01495             {
01496                 int32_t jj=svm->get_support_vector(j);
01497                 sumw[n]+=0.5*svm->get_alpha(i)*svm->get_alpha(j)*kernel->kernel(ii,jj);
01498             }
01499         }
01500         beta[n]=0.0;
01501     }
01502 
01503     mkl_iterations++;
01504     kernel->set_subkernel_weights( (float64_t*) old_beta, num_kernels);
01505 }
01506 
01507 
01508 // assumes that all constraints are satisfied
01509 float64_t CMKL::compute_mkl_dual_objective()
01510 {
01511     if (get_solver_type()==ST_ELASTICNET)
01512     {
01513         // -- Compute ElasticnetMKL dual objective
01514         return compute_elasticnet_dual_objective();
01515     }
01516 
01517     int32_t n=get_num_support_vectors();
01518     float64_t mkl_obj=0;
01519 
01520     if (labels && kernel && kernel->get_kernel_type() == K_COMBINED)
01521     {
01522         CKernel* kn = ((CCombinedKernel*)kernel)->get_first_kernel();
01523         while (kn)
01524         {
01525             float64_t sum=0;
01526             for (int32_t i=0; i<n; i++)
01527             {
01528                 int32_t ii=get_support_vector(i);
01529 
01530                 for (int32_t j=0; j<n; j++)
01531                 {
01532                     int32_t jj=get_support_vector(j);
01533                     sum+=get_alpha(i)*get_alpha(j)*kn->kernel(ii,jj);
01534                 }
01535             }
01536 
01537             if (mkl_norm==1.0)
01538                 mkl_obj = CMath::max(mkl_obj, sum);
01539             else
01540                 mkl_obj += CMath::pow(sum, mkl_norm/(mkl_norm-1));
01541 
01542             SG_UNREF(kn);
01543             kn = ((CCombinedKernel*) kernel)->get_next_kernel();
01544         }
01545 
01546         if (mkl_norm==1.0)
01547             mkl_obj=-0.5*mkl_obj;
01548         else
01549             mkl_obj= -0.5*CMath::pow(mkl_obj, (mkl_norm-1)/mkl_norm);
01550 
01551         mkl_obj+=compute_sum_alpha();
01552     }
01553     else
01554         SG_ERROR( "cannot compute objective, labels or kernel not set\n");
01555 
01556     return -mkl_obj;
01557 }
01558 
01559 #ifdef USE_CPLEX
01560 void CMKL::set_qnorm_constraints(float64_t* beta, int32_t num_kernels)
01561 {
01562     ASSERT(num_kernels>0);
01563 
01564     float64_t* grad_beta=SG_MALLOC(float64_t, num_kernels);
01565     float64_t* hess_beta=SG_MALLOC(float64_t, num_kernels+1);
01566     float64_t* lin_term=SG_MALLOC(float64_t, num_kernels+1);
01567     int* ind=SG_MALLOC(int, num_kernels+1);
01568 
01569     //CMath::display_vector(beta, num_kernels, "beta");
01570     double const_term = 1-CMath::qsq(beta, num_kernels, mkl_norm);
01571 
01572     //SG_PRINT("const=%f\n", const_term);
01573     ASSERT(CMath::fequal(const_term, 0.0));
01574 
01575     for (int32_t i=0; i<num_kernels; i++)
01576     {
01577         grad_beta[i]=mkl_norm * pow(beta[i], mkl_norm-1);
01578         hess_beta[i]=0.5*mkl_norm*(mkl_norm-1) * pow(beta[i], mkl_norm-2);
01579         lin_term[i]=grad_beta[i] - 2*beta[i]*hess_beta[i];
01580         const_term+=grad_beta[i]*beta[i] - CMath::sq(beta[i])*hess_beta[i];
01581         ind[i]=i;
01582     }
01583     ind[num_kernels]=2*num_kernels;
01584     hess_beta[num_kernels]=0;
01585     lin_term[num_kernels]=0;
01586 
01587     int status=0;
01588     int num=CPXgetnumqconstrs (env, lp_cplex);
01589 
01590     if (num>0)
01591     {
01592         status = CPXdelqconstrs (env, lp_cplex, 0, 0);
01593         ASSERT(!status);
01594     }
01595 
01596     status = CPXaddqconstr (env, lp_cplex, num_kernels+1, num_kernels+1, const_term, 'L', ind, lin_term,
01597             ind, ind, hess_beta, NULL);
01598     ASSERT(!status);
01599 
01600     //CPXwriteprob (env, lp_cplex, "prob.lp", NULL);
01601     //CPXqpwrite (env, lp_cplex, "prob.qp");
01602 
01603     SG_FREE(grad_beta);
01604     SG_FREE(hess_beta);
01605     SG_FREE(lin_term);
01606     SG_FREE(ind);
01607 }
01608 #endif // USE_CPLEX
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation