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

SHOGUN Machine Learning Toolbox - Documentation