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  */
00011 
00012 #include "lib/Signal.h"
00013 #include "classifier/mkl/MKL.h"
00014 #include "classifier/svm/LibSVM.h"
00015 #include "kernel/CombinedKernel.h"
00016 
00017 using namespace shogun;
00018 
00019 CMKL::CMKL(CSVM* s)
00020     : CSVM(), svm(NULL), C_mkl(0), mkl_norm(1),
00021     mkl_iterations(0), mkl_epsilon(1e-5), interleaved_optimization(true),
00022     w_gap(1.0), rho(0)
00023 {
00024     set_constraint_generator(s);
00025 #ifdef USE_CPLEX
00026     lp_cplex = NULL ;
00027     env = NULL ;
00028 #endif
00029 
00030 #ifdef USE_GLPK
00031     lp_glpk = NULL;
00032 #endif
00033 
00034     SG_DEBUG("creating MKL object %p\n", this);
00035     lp_initialized = false ;
00036 }
00037 
00038 CMKL::~CMKL()
00039 {
00040     SG_DEBUG("deleting MKL object %p\n", this);
00041     if (svm)
00042         svm->set_callback_function(NULL, NULL);
00043     SG_UNREF(svm);
00044 }
00045 
00046 void CMKL::init_solver()
00047 {
00048 #ifdef USE_CPLEX
00049     cleanup_cplex();
00050 
00051     if (get_solver_type()==ST_CPLEX)
00052         init_cplex();
00053 #endif
00054 
00055 #ifdef USE_GLPK
00056     cleanup_glpk();
00057 
00058     if (get_solver_type() == ST_GLPK)
00059         init_glpk();
00060 #endif
00061 }
00062 
00063 #ifdef USE_CPLEX
00064 bool CMKL::init_cplex()
00065 {
00066     while (env==NULL)
00067     {
00068         SG_INFO( "trying to initialize CPLEX\n") ;
00069 
00070         int status = 0; // calling external lib
00071         env = CPXopenCPLEX (&status);
00072 
00073         if ( env == NULL )
00074         {
00075             char  errmsg[1024];
00076             SG_WARNING( "Could not open CPLEX environment.\n");
00077             CPXgeterrorstring (env, status, errmsg);
00078             SG_WARNING( "%s", errmsg);
00079             SG_WARNING( "retrying in 60 seconds\n");
00080             sleep(60);
00081         }
00082         else
00083         {
00084             // select dual simplex based optimization
00085             status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, CPX_ALG_DUAL);
00086             if ( status )
00087             {
00088             SG_ERROR( "Failure to select dual lp optimization, error %d.\n", status);
00089             }
00090             else
00091             {
00092                 status = CPXsetintparam (env, CPX_PARAM_DATACHECK, CPX_ON);
00093                 if ( status )
00094                 {
00095                     SG_ERROR( "Failure to turn on data checking, error %d.\n", status);
00096                 }   
00097                 else
00098                 {
00099                     lp_cplex = CPXcreateprob (env, &status, "light");
00100 
00101                     if ( lp_cplex == NULL )
00102                         SG_ERROR( "Failed to create LP.\n");
00103                     else
00104                         CPXchgobjsen (env, lp_cplex, CPX_MIN);  /* Problem is minimization */
00105                 }
00106             }
00107         }
00108     }
00109 
00110     return (lp_cplex != NULL) && (env != NULL);
00111 }
00112 
00113 bool CMKL::cleanup_cplex()
00114 {
00115     bool result=false;
00116 
00117     if (lp_cplex)
00118     {
00119         int32_t status = CPXfreeprob(env, &lp_cplex);
00120         lp_cplex = NULL;
00121         lp_initialized = false;
00122 
00123         if (status)
00124             SG_WARNING( "CPXfreeprob failed, error code %d.\n", status);
00125         else
00126             result = true;
00127     }
00128 
00129     if (env)
00130     {
00131         int32_t status = CPXcloseCPLEX (&env);
00132         env=NULL;
00133         
00134         if (status)
00135         {
00136             char  errmsg[1024];
00137             SG_WARNING( "Could not close CPLEX environment.\n");
00138             CPXgeterrorstring (env, status, errmsg);
00139             SG_WARNING( "%s", errmsg);
00140         }
00141         else
00142             result = true;
00143     }
00144     return result;
00145 }
00146 #endif
00147 
00148 #ifdef USE_GLPK
00149 bool CMKL::init_glpk()
00150 {
00151     lp_glpk = lpx_create_prob();
00152     lpx_set_obj_dir(lp_glpk, LPX_MIN);
00153     lpx_set_int_parm(lp_glpk, LPX_K_DUAL, GLP_ON );
00154     lpx_set_int_parm(lp_glpk, LPX_K_PRESOL, GLP_ON );
00155 
00156     glp_term_out(GLP_OFF);
00157     return (lp_glpk != NULL);
00158 }
00159 
00160 bool CMKL::cleanup_glpk()
00161 {
00162     lp_initialized = false;
00163     if (lp_glpk)
00164         lpx_delete_prob(lp_glpk);
00165     lp_glpk = NULL;
00166     return true;
00167 }
00168 
00169 bool CMKL::check_lpx_status(LPX *lp)
00170 {
00171     int status = lpx_get_status(lp);
00172 
00173     if (status==LPX_INFEAS)
00174     {
00175         SG_PRINT("solution is infeasible!\n");
00176         return false;
00177     }
00178     else if(status==LPX_NOFEAS)
00179     {
00180         SG_PRINT("problem has no feasible solution!\n");
00181         return false;
00182     }
00183     return true;
00184 }
00185 #endif // USE_GLPK
00186 
00187 bool CMKL::train(CFeatures* data)
00188 {
00189     ASSERT(kernel);
00190     ASSERT(labels && labels->get_num_labels());
00191 
00192     if (data)
00193     {
00194         if (labels->get_num_labels() != data->get_num_vectors())
00195             SG_ERROR("Number of training vectors does not match number of labels\n");
00196         kernel->init(data, data);
00197     }
00198 
00199     init_training();
00200     if (!svm)
00201         SG_ERROR("No constraint generator (SVM) set\n");
00202 
00203     int32_t num_label=0;
00204     if (labels)
00205         labels->get_num_labels();
00206 
00207     SG_INFO("%d trainlabels\n", num_label);
00208     if (mkl_epsilon<=0)
00209         mkl_epsilon=1e-2 ;
00210 
00211     SG_INFO("mkl_epsilon = %1.1e\n", mkl_epsilon) ;
00212     SG_INFO("C_mkl = %1.1e\n", C_mkl) ;
00213     SG_INFO("mkl_norm = %1.3e\n", mkl_norm);
00214     SG_INFO("solver = %d\n", get_solver_type());
00215 
00216     int32_t num_weights = -1;
00217     int32_t num_kernels = kernel->get_num_subkernels();
00218     SG_INFO("num_kernels = %d\n", num_kernels);
00219     const float64_t* beta_const   = kernel->get_subkernel_weights(num_weights);
00220     float64_t* beta =  CMath::clone_vector(beta_const, num_weights);
00221     ASSERT(num_weights==num_kernels);
00222     CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm), beta, num_kernels); //q-norm = 1
00223     kernel->set_subkernel_weights(beta, num_kernels);
00224     delete[] beta;
00225 
00226     svm->set_bias_enabled(get_bias_enabled());
00227     svm->set_epsilon(get_epsilon());
00228     svm->set_max_train_time(get_max_train_time());
00229     svm->set_nu(get_nu());
00230     svm->set_C(get_C1(), get_C2());
00231     svm->set_qpsize(get_qpsize());
00232     svm->set_shrinking_enabled(get_shrinking_enabled());
00233     svm->set_linadd_enabled(get_linadd_enabled());
00234     svm->set_batch_computation_enabled(get_batch_computation_enabled());
00235     svm->set_labels(labels);
00236     svm->set_kernel(kernel);
00237 
00238 #ifdef USE_CPLEX
00239     cleanup_cplex();
00240 
00241     if (get_solver_type()==ST_CPLEX)
00242         init_cplex();
00243 #endif
00244 
00245 #ifdef USE_GLPK
00246     if (get_solver_type()==ST_GLPK)
00247         init_glpk();
00248 #endif
00249 
00250     mkl_iterations = 0;
00251     CSignal::clear_cancel();
00252     
00253     if (interleaved_optimization)
00254     {
00255         if (svm->get_classifier_type() != CT_LIGHT && svm->get_classifier_type() != CT_SVRLIGHT)
00256         {
00257             SG_ERROR("Interleaved MKL optimization is currently "
00258                     "only supported with SVMlight\n");
00259         }
00260 
00261         //Note that there is currently no safe way to properly resolve this
00262         //situation: the mkl object hands itself as a reference to the svm which
00263         //in turn increases the reference count of mkl and when done decreases
00264         //the count. Thus we have to protect the mkl object from deletion by
00265         //ref()'ing it before we set the callback function and should also
00266         //unref() it afterwards. However, when the reference count was zero
00267         //before this unref() might actually try to destroy this (crash ahead)
00268         //but if we don't actually unref() the object we might leak memory...
00269         //So as a workaround we only unref when the reference count was >1
00270         //before.
00271         int32_t refs=this->ref();
00272         svm->set_callback_function(this, perform_mkl_step_helper);
00273         svm->train();
00274         SG_DONE();
00275         svm->set_callback_function(NULL, NULL);
00276         if (refs>1)
00277             this->unref();
00278     }
00279     else
00280     {
00281         float64_t* sumw = new float64_t[num_kernels];
00282 
00283         while (true)
00284         {
00285             svm->train();
00286 
00287             float64_t suma=compute_sum_alpha();
00288             compute_sum_beta(sumw);
00289 
00290             mkl_iterations++;
00291             if (perform_mkl_step(sumw, suma) || CSignal::cancel_computations())
00292                 break;
00293         }
00294 
00295         delete[] sumw;
00296     }
00297 #ifdef USE_CPLEX
00298     cleanup_cplex();
00299 #endif
00300 #ifdef USE_GLPK
00301     cleanup_glpk();
00302 #endif
00303 
00304     int32_t nsv=svm->get_num_support_vectors();
00305     create_new_model(nsv);
00306 
00307     set_bias(svm->get_bias());
00308     for (int32_t i=0; i<nsv; i++)
00309     {
00310         set_alpha(i, svm->get_alpha(i));
00311         set_support_vector(i, svm->get_support_vector(i));
00312     }
00313     return true;
00314 }
00315 
00316 bool CMKL::perform_mkl_step(
00317         const float64_t* sumw, float64_t suma)
00318 {
00319     int32_t num_kernels = kernel->get_num_subkernels();
00320     int32_t nweights=0;
00321     const float64_t* old_beta = kernel->get_subkernel_weights(nweights);
00322     ASSERT(nweights==num_kernels);
00323     float64_t* beta = new float64_t[num_kernels];
00324 
00325     int32_t inner_iters=0;
00326     float64_t mkl_objective=0;
00327 
00328     mkl_objective=-suma;
00329     for (int32_t i=0; i<num_kernels; i++)
00330     {
00331         beta[i]=old_beta[i];
00332         mkl_objective+=old_beta[i]*sumw[i];
00333     }
00334 
00335     w_gap = CMath::abs(1-rho/mkl_objective) ;
00336 
00337     if( (w_gap >= mkl_epsilon) ||
00338             (get_solver_type()==ST_AUTO || get_solver_type()==ST_NEWTON || get_solver_type()==ST_DIRECT ) )
00339     {
00340         if (get_solver_type()==ST_AUTO || get_solver_type()==ST_DIRECT)
00341             rho=compute_optimal_betas_directly(beta, old_beta, num_kernels, sumw, suma, mkl_objective);
00342         else if (get_solver_type()==ST_NEWTON)
00343             rho=compute_optimal_betas_newton(beta, old_beta, num_kernels, sumw, suma, mkl_objective);
00344 #ifdef USE_CPLEX
00345         else if (get_solver_type()==ST_CPLEX)
00346             rho=compute_optimal_betas_via_cplex(beta, old_beta, num_kernels, sumw, suma, inner_iters);
00347 #endif
00348 #ifdef USE_GLPK
00349         else if (get_solver_type()==ST_GLPK)
00350             rho=compute_optimal_betas_via_glpk(beta, old_beta, num_kernels, sumw, suma, inner_iters);
00351 #endif
00352         else
00353             SG_ERROR("Solver type not supported (not compiled in?)\n");
00354         w_gap = CMath::abs(1-rho/mkl_objective) ;
00355     }
00356 
00357     kernel->set_subkernel_weights(beta, num_kernels);
00358     delete[] beta;
00359 
00360     return converged();
00361 }
00362 
00363 
00364 float64_t CMKL::compute_optimal_betas_directly(
00365   float64_t* beta, const float64_t* old_beta, const int32_t num_kernels,
00366   const float64_t* sumw, const float64_t suma,
00367   const float64_t mkl_objective )
00368 {
00369   const float64_t epsRegul = 0.01;  // fraction of root mean squared deviation
00370   float64_t obj;
00371   float64_t Z;
00372   float64_t preR;
00373   int32_t p;
00374   int32_t nofKernelsGood;
00375 
00376   // --- optimal beta
00377   nofKernelsGood = num_kernels;
00378   for( p=0; p<num_kernels; ++p ) {
00379     //SG_PRINT( "MKL-direct:  sumw[%3d] = %e  ( oldbeta = %e )\n", p, sumw[p], old_beta[p] );
00380     if( sumw[p] >= 0.0 && old_beta[p] >= 0.0 ) {
00381       beta[p] = sumw[p] * old_beta[p]*old_beta[p] / mkl_norm;
00382       beta[p] = CMath::pow( beta[p], 1.0 / (mkl_norm+1.0) );
00383     } else {
00384       beta[p] = 0.0;
00385       --nofKernelsGood;
00386     }
00387     ASSERT( beta[p] >= 0 );
00388   }
00389 
00390   // --- normalize
00391   Z = 0.0;
00392   for( p=0; p<num_kernels; ++p ) {
00393     Z += CMath::pow( beta[p], mkl_norm );
00394   }
00395   Z = CMath::pow( Z, -1.0/mkl_norm );
00396   ASSERT( Z >= 0 );
00397   for( p=0; p<num_kernels; ++p ) {
00398     beta[p] *= Z;
00399   }
00400 
00401   // --- regularize & renormalize
00402   preR = 0.0;
00403   for( p=0; p<num_kernels; ++p ) {
00404     preR += CMath::pow( old_beta[p] - beta[p], 2.0 );
00405   }
00406   const float64_t R = CMath::sqrt( preR / mkl_norm ) * epsRegul;
00407   if( !( R >= 0 ) ) {
00408     SG_PRINT( "MKL-direct: p = %.3f\n", mkl_norm );
00409     SG_PRINT( "MKL-direct: nofKernelsGood = %d\n", nofKernelsGood );
00410     SG_PRINT( "MKL-direct: Z = %e\n", Z );
00411     SG_PRINT( "MKL-direct: eps = %e\n", epsRegul );
00412     for( p=0; p<num_kernels; ++p ) {
00413       const float64_t t = CMath::pow( old_beta[p] - beta[p], 2.0 );
00414       SG_PRINT( "MKL-direct: t[%3d] = %e  ( diff = %e = %e - %e )\n", p, t, old_beta[p]-beta[p], old_beta[p], beta[p] );
00415     }
00416     SG_PRINT( "MKL-direct: preR = %e\n", preR );
00417     SG_PRINT( "MKL-direct: preR/p = %e\n", preR/mkl_norm );
00418     SG_PRINT( "MKL-direct: sqrt(preR/p) = %e\n", CMath::sqrt(preR/mkl_norm) );
00419     SG_PRINT( "MKL-direct: R = %e\n", R );
00420     SG_ERROR( "Assertion R >= 0 failed!\n" );
00421   }
00422 
00423   Z = 0.0;
00424   for( p=0; p<num_kernels; ++p ) {
00425     beta[p] += R;
00426     Z += CMath::pow( beta[p], mkl_norm );
00427     ASSERT( beta[p] >= 0 );
00428   }
00429   Z = CMath::pow( Z, -1.0/mkl_norm );
00430   ASSERT( Z >= 0 );
00431   for( p=0; p<num_kernels; ++p ) {
00432     beta[p] *= Z;
00433     ASSERT( beta[p] >= 0.0 );
00434     if( beta[p] > 1.0 ) {
00435       beta[p] = 1.0;
00436     }
00437   }
00438 
00439   // --- objective
00440   obj = -suma;
00441   for( p=0; p<num_kernels; ++p ) {
00442     //obj += sumw[p] * old_beta[p]*old_beta[p] / beta[p];
00443     obj += sumw[p] * beta[p];
00444   }
00445   return obj;
00446 }
00447 
00448 float64_t CMKL::compute_optimal_betas_newton(float64_t* beta,
00449         const float64_t* old_beta, int32_t num_kernels,
00450         const float64_t* sumw, float64_t suma,
00451          float64_t mkl_objective)
00452 {
00453   SG_DEBUG("MKL via NEWTON\n");
00454 
00455   if (mkl_norm==1)
00456       SG_ERROR("MKL via NEWTON works only for norms>1\n");
00457 
00458   const double epsBeta = 1e-32;
00459   const double epsGamma = 1e-12;
00460   const double epsWsq = 1e-12;
00461   const double epsNewt = 0.0001;
00462   const double epsStep = 1e-9;
00463   const int nofNewtonSteps = 3;
00464   const double hessRidge = 1e-6;
00465   const int inLogSpace = 0;
00466 
00467   const float64_t r = mkl_norm / ( mkl_norm - 1.0 );
00468   float64_t* newtDir = new float64_t[ num_kernels ];
00469   float64_t* newtBeta = new float64_t[ num_kernels ];
00470   float64_t newtStep;
00471   float64_t stepSize;
00472   float64_t Z;
00473   float64_t obj;
00474   float64_t gamma;
00475   int32_t p;
00476   int i;
00477 
00478 
00479   // === init
00480   for( p=0; p<num_kernels; ++p ) {
00481     //SG_PRINT( "old_beta[%d] = %e;  sumw[%d]=%e.  \n", p, old_beta[p], p, sumw[p] );
00482   }
00483 
00484   // --- check beta
00485   Z = 0.0;
00486   for( p=0; p<num_kernels; ++p ) {
00487     beta[p] = old_beta[p];
00488     if( !( beta[p] >= epsBeta ) ) {
00489       //SG_WARNING( "old_beta[%d] = %e  (sumw[.]=%e);  set to %e.  ", p, beta[p], sumw[p], epsBeta );
00490       beta[p] = epsBeta;
00491     }
00492     ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
00493     Z += CMath::pow( beta[p], mkl_norm );
00494   }
00495   Z = CMath::pow( Z, -1.0/mkl_norm );
00496   if( !( fabs(Z-1.0) <= epsGamma ) ) {
00497     SG_WARNING( "old_beta not normalized (diff=%e);  forcing normalization.  ", Z-1.0 );
00498     for( p=0; p<num_kernels; ++p ) {
00499       beta[p] *= Z;
00500       if( beta[p] > 1.0 ) {
00501     beta[p] = 1.0;
00502       }
00503       ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
00504     }
00505   }
00506 
00507   // --- compute gamma
00508   gamma = 0.0;
00509   for( p=0; p<num_kernels; ++p ) {
00510     if( !( sumw[p] >= 0 ) ) {
00511       if( !( sumw[p] >= -epsWsq ) ) {
00512     SG_WARNING( "sumw[%d] = %e;  treated as 0.  ", p, sumw[p] );
00513       }
00514       // should better recompute sumw[] !!!
00515     } else {
00516       ASSERT( sumw[p] >= 0 );
00517       //gamma += CMath::pow( sumw[p] * beta[p]*beta[p], r );
00518       gamma += CMath::pow( sumw[p] * beta[p]*beta[p] / mkl_norm, r );
00519     }
00520   }
00521   gamma = CMath::pow( gamma, 1.0/r ) / mkl_norm;
00522   ASSERT( gamma > -1e-9 );
00523   if( !( gamma > epsGamma ) ) {
00524     SG_WARNING( "bad gamma: %e;  set to %e.  ", gamma, epsGamma );
00525     // should better recompute sumw[] !!!
00526     gamma = epsGamma;
00527   }
00528   ASSERT( gamma >= epsGamma );
00529   //gamma = -gamma;
00530 
00531   // --- compute objective
00532   obj = 0.0;
00533   for( p=0; p<num_kernels; ++p ) {
00534     obj += beta[p] * sumw[p];
00535     //obj += gamma/mkl_norm * CMath::pow( beta[p], mkl_norm );
00536   }
00537   if( !( obj >= 0.0 ) ) {
00538     SG_WARNING( "negative objective: %e.  ", obj );
00539   }
00540   //SG_PRINT( "OBJ = %e.  \n", obj );
00541 
00542 
00543   // === perform Newton steps
00544   if( nofNewtonSteps > 1 ) {
00545     //SG_DEBUG( "performing %d Newton steps.\n", nofNewtonSteps );
00546   }
00547   for( i = 0; i < nofNewtonSteps; ++i ) {
00548 
00549     // --- compute Newton direction (Hessian is diagonal)
00550     const float64_t gqq1 = mkl_norm * (mkl_norm-1.0) * gamma;
00551     newtStep = 0.0;
00552     for( p=0; p<num_kernels; ++p ) {
00553       ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
00554       //const float halfw2p = ( sumw[p] >= 0.0 ) ? sumw[p] : 0.0;
00555       //const float64_t t1 = halfw2p*beta[p] - mkl_norm*gamma*CMath::pow(beta[p],mkl_norm);
00556       //const float64_t t2 = 2.0*halfw2p + gqq1*CMath::pow(beta[p],mkl_norm-1.0);
00557       const float halfw2p = ( sumw[p] >= 0.0 ) ? (sumw[p]*old_beta[p]*old_beta[p]) : 0.0;
00558       const float64_t t0 = halfw2p*beta[p] - mkl_norm*gamma*CMath::pow(beta[p],mkl_norm+2.0);
00559       const float64_t t1 = ( t0 < 0 ) ? 0.0 : t0;
00560       const float64_t t2 = 2.0*halfw2p + gqq1*CMath::pow(beta[p],mkl_norm+1.0);
00561       if( inLogSpace ) {
00562     newtDir[p] = t1 / ( t1 + t2*beta[p] + hessRidge );
00563       } else {
00564     newtDir[p] = ( t1 == 0.0 ) ? 0.0 : ( t1 / t2 );
00565       }
00566       // newtStep += newtDir[p] * grad[p];
00567       ASSERT( newtDir[p] == newtDir[p] );
00568       //SG_PRINT( "newtDir[%d] = %6.3f = %e / %e \n", p, newtDir[p], t1, t2 );
00569     }
00570     //CMath::display_vector( newtDir, num_kernels, "newton direction  " );
00571     //SG_PRINT( "Newton step size = %e\n", Z );
00572 
00573     // --- line search
00574     stepSize = 1.0;
00575     while( stepSize >= epsStep ) {
00576 
00577       // --- perform Newton step
00578       Z = 0.0;
00579       while( Z == 0.0 ) {
00580     for( p=0; p<num_kernels; ++p ) {
00581       if( inLogSpace ) {
00582         newtBeta[p] = beta[p] * CMath::exp( + stepSize * newtDir[p] );
00583       } else {
00584         newtBeta[p] = beta[p] + stepSize * newtDir[p];
00585       }
00586       if( !( newtBeta[p] >= epsBeta ) ) {
00587         newtBeta[p] = epsBeta;
00588       }
00589       Z += CMath::pow( newtBeta[p], mkl_norm );
00590     }
00591     ASSERT( 0.0 <= Z );
00592     Z = CMath::pow( Z, -1.0/mkl_norm );
00593     if( Z == 0.0 ) {
00594       stepSize /= 2.0;
00595     }
00596       }
00597 
00598       // --- normalize new beta (wrt p-norm)
00599       ASSERT( 0.0 < Z );
00600       ASSERT( Z < CMath::INFTY );
00601       for( p=0; p<num_kernels; ++p ) {
00602     newtBeta[p] *= Z;
00603     if( newtBeta[p] > 1.0 ) {
00604       //SG_WARNING( "beta[%d] = %e;  set to 1.  ", p, beta[p] );
00605       newtBeta[p] = 1.0;
00606     }
00607     ASSERT( 0.0 <= newtBeta[p] && newtBeta[p] <= 1.0 );
00608       }
00609 
00610       // --- objective increased?
00611       float64_t newtObj;
00612       newtObj = 0.0;
00613       for( p=0; p<num_kernels; ++p ) {
00614     newtObj += sumw[p] * old_beta[p]*old_beta[p] / newtBeta[p];
00615       }
00616       //SG_PRINT( "step = %.8f:  obj = %e -> %e.  \n", stepSize, obj, newtObj );
00617       if( newtObj < obj - epsNewt*stepSize*obj ) {
00618     for( p=0; p<num_kernels; ++p ) {
00619       beta[p] = newtBeta[p];
00620     }
00621     obj = newtObj;
00622     break;
00623       }
00624       stepSize /= 2.0;
00625 
00626     }
00627 
00628     if( stepSize < epsStep ) {
00629       break;
00630     }
00631 
00632   }
00633   delete[] newtDir;
00634   delete[] newtBeta;
00635   //CMath::scale_vector( 1.0/CMath::qnorm(beta,num_kernels,mkl_norm), beta, num_kernels );
00636 
00637 
00638   // === return new objective
00639   obj = -suma;
00640   for( p=0; p<num_kernels; ++p ) {
00641     obj += beta[p] * sumw[p];
00642     //obj += sumw[p] * old_beta[p]*old_beta[p] / beta[p];
00643   }
00644   return obj;
00645 }
00646 
00647 
00648 
00649 float64_t CMKL::compute_optimal_betas_via_cplex(float64_t* new_beta, const float64_t* old_beta, int32_t num_kernels,
00650           const float64_t* sumw, float64_t suma, int32_t& inner_iters)
00651 {
00652     SG_DEBUG("MKL via CPLEX\n");
00653 
00654 #ifdef USE_CPLEX
00655     ASSERT(new_beta);
00656     ASSERT(old_beta);
00657 
00658     int32_t NUMCOLS = 2*num_kernels + 1;
00659     double* x=new double[NUMCOLS];
00660 
00661     if (!lp_initialized)
00662     {
00663         SG_INFO( "creating LP\n") ;
00664 
00665         double   obj[NUMCOLS]; /* calling external lib */
00666         double   lb[NUMCOLS]; /* calling external lib */
00667         double   ub[NUMCOLS]; /* calling external lib */
00668 
00669         for (int32_t i=0; i<2*num_kernels; i++)
00670         {
00671             obj[i]=0 ;
00672             lb[i]=0 ;
00673             ub[i]=1 ;
00674         }
00675 
00676         for (int32_t i=num_kernels; i<2*num_kernels; i++)
00677             obj[i]= C_mkl;
00678 
00679         obj[2*num_kernels]=1 ;
00680         lb[2*num_kernels]=-CPX_INFBOUND ;
00681         ub[2*num_kernels]=CPX_INFBOUND ;
00682 
00683         int status = CPXnewcols (env, lp_cplex, NUMCOLS, obj, lb, ub, NULL, NULL);
00684         if ( status ) {
00685             char  errmsg[1024];
00686             CPXgeterrorstring (env, status, errmsg);
00687             SG_ERROR( "%s", errmsg);
00688         }
00689 
00690         // add constraint sum(w)=1;
00691         SG_INFO( "adding the first row\n");
00692         int initial_rmatbeg[1]; /* calling external lib */
00693         int initial_rmatind[num_kernels+1]; /* calling external lib */
00694         double initial_rmatval[num_kernels+1]; /* calling ext lib */
00695         double initial_rhs[1]; /* calling external lib */
00696         char initial_sense[1];
00697 
00698         // 1-norm MKL
00699         if (mkl_norm==1)
00700         {
00701             initial_rmatbeg[0] = 0;
00702             initial_rhs[0]=1 ;     // rhs=1 ;
00703             initial_sense[0]='E' ; // equality
00704 
00705             for (int32_t i=0; i<num_kernels; i++)
00706             {
00707                 initial_rmatind[i]=i ;
00708                 initial_rmatval[i]=1 ;
00709             }
00710             initial_rmatind[num_kernels]=2*num_kernels ;
00711             initial_rmatval[num_kernels]=0 ;
00712 
00713             status = CPXaddrows (env, lp_cplex, 0, 1, num_kernels+1, 
00714                     initial_rhs, initial_sense, initial_rmatbeg,
00715                     initial_rmatind, initial_rmatval, NULL, NULL);
00716 
00717         }
00718         else // 2 and q-norm MKL
00719         {
00720             initial_rmatbeg[0] = 0;
00721             initial_rhs[0]=0 ;     // rhs=1 ;
00722             initial_sense[0]='L' ; // <=  (inequality)
00723 
00724             initial_rmatind[0]=2*num_kernels ;
00725             initial_rmatval[0]=0 ;
00726 
00727             status = CPXaddrows (env, lp_cplex, 0, 1, 1, 
00728                     initial_rhs, initial_sense, initial_rmatbeg,
00729                     initial_rmatind, initial_rmatval, NULL, NULL);
00730 
00731 
00732             if (mkl_norm==2)
00733             {
00734                 for (int32_t i=0; i<num_kernels; i++)
00735                 {
00736                     initial_rmatind[i]=i ;
00737                     initial_rmatval[i]=1 ;
00738                 }
00739                 initial_rmatind[num_kernels]=2*num_kernels ;
00740                 initial_rmatval[num_kernels]=0 ;
00741 
00742                 status = CPXaddqconstr (env, lp_cplex, 0, num_kernels+1, 1.0, 'L', NULL, NULL,
00743                         initial_rmatind, initial_rmatind, initial_rmatval, NULL);
00744             }
00745         }
00746 
00747 
00748         if ( status )
00749             SG_ERROR( "Failed to add the first row.\n");
00750 
00751         lp_initialized = true ;
00752 
00753         if (C_mkl!=0.0)
00754         {
00755             for (int32_t q=0; q<num_kernels-1; q++)
00756             {
00757                 // add constraint w[i]-w[i+1]<s[i];
00758                 // add constraint w[i+1]-w[i]<s[i];
00759                 int rmatbeg[1]; /* calling external lib */
00760                 int rmatind[3]; /* calling external lib */
00761                 double rmatval[3]; /* calling external lib */
00762                 double rhs[1]; /* calling external lib */
00763                 char sense[1];
00764 
00765                 rmatbeg[0] = 0;
00766                 rhs[0]=0 ;     // rhs=0 ;
00767                 sense[0]='L' ; // <=
00768                 rmatind[0]=q ;
00769                 rmatval[0]=1 ;
00770                 rmatind[1]=q+1 ;
00771                 rmatval[1]=-1 ;
00772                 rmatind[2]=num_kernels+q ;
00773                 rmatval[2]=-1 ;
00774                 status = CPXaddrows (env, lp_cplex, 0, 1, 3, 
00775                         rhs, sense, rmatbeg,
00776                         rmatind, rmatval, NULL, NULL);
00777                 if ( status )
00778                     SG_ERROR( "Failed to add a smothness row (1).\n");
00779 
00780                 rmatbeg[0] = 0;
00781                 rhs[0]=0 ;     // rhs=0 ;
00782                 sense[0]='L' ; // <=
00783                 rmatind[0]=q ;
00784                 rmatval[0]=-1 ;
00785                 rmatind[1]=q+1 ;
00786                 rmatval[1]=1 ;
00787                 rmatind[2]=num_kernels+q ;
00788                 rmatval[2]=-1 ;
00789                 status = CPXaddrows (env, lp_cplex, 0, 1, 3, 
00790                         rhs, sense, rmatbeg,
00791                         rmatind, rmatval, NULL, NULL);
00792                 if ( status )
00793                     SG_ERROR( "Failed to add a smothness row (2).\n");
00794             }
00795         }
00796     }
00797 
00798     { // add the new row
00799         //SG_INFO( "add the new row\n") ;
00800 
00801         int rmatbeg[1];
00802         int rmatind[num_kernels+1];
00803         double rmatval[num_kernels+1];
00804         double rhs[1];
00805         char sense[1];
00806 
00807         rmatbeg[0] = 0;
00808         if (mkl_norm==1)
00809             rhs[0]=0 ;
00810         else
00811             rhs[0]=-suma ;
00812 
00813         sense[0]='L' ;
00814 
00815         for (int32_t i=0; i<num_kernels; i++)
00816         {
00817             rmatind[i]=i ;
00818             if (mkl_norm==1)
00819                 rmatval[i]=-(sumw[i]-suma) ;
00820             else
00821                 rmatval[i]=-sumw[i];
00822         }
00823         rmatind[num_kernels]=2*num_kernels ;
00824         rmatval[num_kernels]=-1 ;
00825 
00826         int32_t status = CPXaddrows (env, lp_cplex, 0, 1, num_kernels+1, 
00827                 rhs, sense, rmatbeg,
00828                 rmatind, rmatval, NULL, NULL);
00829         if ( status ) 
00830             SG_ERROR( "Failed to add the new row.\n");
00831     }
00832 
00833     inner_iters=0;
00834     int status;
00835     { 
00836 
00837         if (mkl_norm==1) // optimize 1 norm MKL
00838             status = CPXlpopt (env, lp_cplex);
00839         else if (mkl_norm==2) // optimize 2-norm MKL
00840             status = CPXbaropt(env, lp_cplex);
00841         else // q-norm MKL
00842         {
00843             float64_t* beta=new float64_t[2*num_kernels+1];
00844             float64_t objval_old=1e-8; //some value to cause the loop to not terminate yet
00845             for (int32_t i=0; i<num_kernels; i++)
00846                 beta[i]=old_beta[i];
00847             for (int32_t i=num_kernels; i<2*num_kernels+1; i++)
00848                 beta[i]=0;
00849 
00850             while (true)
00851             {
00852                 //int rows=CPXgetnumrows(env, lp_cplex);
00853                 //int cols=CPXgetnumcols(env, lp_cplex);
00854                 //SG_PRINT("rows:%d, cols:%d (kernel:%d)\n", rows, cols, num_kernels);
00855                 CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm), beta, num_kernels);
00856 
00857                 set_qnorm_constraints(beta, num_kernels);
00858 
00859                 status = CPXbaropt(env, lp_cplex);
00860                 if ( status ) 
00861                     SG_ERROR( "Failed to optimize Problem.\n");
00862 
00863                 int solstat=0;
00864                 double objval=0;
00865                 status=CPXsolution(env, lp_cplex, &solstat, &objval,
00866                         (double*) beta, NULL, NULL, NULL);
00867 
00868                 if ( status )
00869                 {
00870                     CMath::display_vector(beta, num_kernels, "beta");
00871                     SG_ERROR( "Failed to obtain solution.\n");
00872                 }
00873 
00874                 CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm), beta, num_kernels);
00875 
00876                 //SG_PRINT("[%d] %f (%f)\n", inner_iters, objval, objval_old);
00877                 if ((1-abs(objval/objval_old) < 0.1*mkl_epsilon)) // && (inner_iters>2))
00878                     break;
00879 
00880                 objval_old=objval;
00881 
00882                 inner_iters++;
00883             }
00884             delete[] beta;
00885         }
00886 
00887         if ( status ) 
00888             SG_ERROR( "Failed to optimize Problem.\n");
00889 
00890         // obtain solution
00891         int32_t cur_numrows=(int32_t) CPXgetnumrows(env, lp_cplex);
00892         int32_t cur_numcols=(int32_t) CPXgetnumcols(env, lp_cplex);
00893         int32_t num_rows=cur_numrows;
00894         ASSERT(cur_numcols<=2*num_kernels+1);
00895 
00896         float64_t* slack=new float64_t[cur_numrows];
00897         float64_t* pi=new float64_t[cur_numrows];
00898 
00899         /* calling external lib */
00900         int solstat=0;
00901         double objval=0;
00902 
00903         if (mkl_norm==1)
00904         {
00905             status=CPXsolution(env, lp_cplex, &solstat, &objval,
00906                     (double*) x, (double*) pi, (double*) slack, NULL);
00907         }
00908         else
00909         {
00910             status=CPXsolution(env, lp_cplex, &solstat, &objval,
00911                     (double*) x, NULL, (double*) slack, NULL);
00912         }
00913 
00914         int32_t solution_ok = (!status) ;
00915         if ( status )
00916             SG_ERROR( "Failed to obtain solution.\n");
00917 
00918         int32_t num_active_rows=0 ;
00919         if (solution_ok)
00920         {
00921             /* 1 norm mkl */
00922             float64_t max_slack = -CMath::INFTY ;
00923             int32_t max_idx = -1 ;
00924             int32_t start_row = 1 ;
00925             if (C_mkl!=0.0)
00926                 start_row+=2*(num_kernels-1);
00927 
00928             for (int32_t i = start_row; i < cur_numrows; i++)  // skip first
00929             {
00930                 if (mkl_norm==1)
00931                 {
00932                     if ((pi[i]!=0))
00933                         num_active_rows++ ;
00934                     else
00935                     {
00936                         if (slack[i]>max_slack)
00937                         {
00938                             max_slack=slack[i] ;
00939                             max_idx=i ;
00940                         }
00941                     }
00942                 }
00943                 else // 2-norm or general q-norm
00944                 {
00945                     if ((CMath::abs(slack[i])<1e-6))
00946                         num_active_rows++ ;
00947                     else
00948                     {
00949                         if (slack[i]>max_slack)
00950                         {
00951                             max_slack=slack[i] ;
00952                             max_idx=i ;
00953                         }
00954                     }
00955                 }
00956             }
00957 
00958             // have at most max(100,num_active_rows*2) rows, if not, remove one
00959             if ( (num_rows-start_row>CMath::max(100,2*num_active_rows)) && (max_idx!=-1))
00960             {
00961                 //SG_INFO( "-%i(%i,%i)",max_idx,start_row,num_rows) ;
00962                 status = CPXdelrows (env, lp_cplex, max_idx, max_idx) ;
00963                 if ( status ) 
00964                     SG_ERROR( "Failed to remove an old row.\n");
00965             }
00966 
00967             //CMath::display_vector(x, num_kernels, "beta");
00968 
00969             rho = -x[2*num_kernels] ;
00970             delete[] pi ;
00971             delete[] slack ;
00972 
00973         }
00974         else
00975         {
00976             /* then something is wrong and we rather 
00977             stop sooner than later */
00978             rho = 1 ;
00979         }
00980     }
00981     for (int32_t i=0; i<num_kernels; i++)
00982         new_beta[i]=x[i];
00983 
00984     delete[] x;
00985 #else
00986     SG_ERROR("Cplex not enabled at compile time\n");
00987 #endif
00988     return rho;
00989 }
00990 
00991 float64_t CMKL::compute_optimal_betas_via_glpk(float64_t* beta, const float64_t* old_beta,
00992         int num_kernels, const float64_t* sumw, float64_t suma, int32_t& inner_iters)
00993 {
00994     SG_DEBUG("MKL via GLPK\n");
00995 
00996     if (mkl_norm!=1)
00997         SG_ERROR("MKL via GLPK works only for norm=1\n");
00998 
00999     float64_t obj=1.0;
01000 #ifdef USE_GLPK
01001     int32_t NUMCOLS = 2*num_kernels + 1 ;
01002     if (!lp_initialized)
01003     {
01004         SG_INFO( "creating LP\n") ;
01005 
01006         //set obj function, note glpk indexing is 1-based
01007         lpx_add_cols(lp_glpk, NUMCOLS);
01008         for (int i=1; i<=2*num_kernels; i++)
01009         {
01010             lpx_set_obj_coef(lp_glpk, i, 0);
01011             lpx_set_col_bnds(lp_glpk, i, LPX_DB, 0, 1);
01012         }
01013         for (int i=num_kernels+1; i<=2*num_kernels; i++)
01014         {
01015             lpx_set_obj_coef(lp_glpk, i, C_mkl);
01016         }
01017         lpx_set_obj_coef(lp_glpk, NUMCOLS, 1);
01018         lpx_set_col_bnds(lp_glpk, NUMCOLS, LPX_FR, 0,0); //unbound
01019 
01020         //add first row. sum[w]=1
01021         int row_index = lpx_add_rows(lp_glpk, 1);
01022         int ind[num_kernels+2];
01023         float64_t val[num_kernels+2];
01024         for (int i=1; i<=num_kernels; i++)
01025         {
01026             ind[i] = i;
01027             val[i] = 1;
01028         }
01029         ind[num_kernels+1] = NUMCOLS;
01030         val[num_kernels+1] = 0;
01031         lpx_set_mat_row(lp_glpk, row_index, num_kernels, ind, val);
01032         lpx_set_row_bnds(lp_glpk, row_index, LPX_FX, 1, 1);
01033 
01034         lp_initialized = true;
01035 
01036         if (C_mkl!=0.0)
01037         {
01038             for (int32_t q=1; q<num_kernels; q++)
01039             {
01040                 int mat_ind[4];
01041                 float64_t mat_val[4];
01042                 int mat_row_index = lpx_add_rows(lp_glpk, 2);
01043                 mat_ind[1] = q;
01044                 mat_val[1] = 1;
01045                 mat_ind[2] = q+1; 
01046                 mat_val[2] = -1;
01047                 mat_ind[3] = num_kernels+q;
01048                 mat_val[3] = -1;
01049                 lpx_set_mat_row(lp_glpk, mat_row_index, 3, mat_ind, mat_val);
01050                 lpx_set_row_bnds(lp_glpk, mat_row_index, LPX_UP, 0, 0);
01051                 mat_val[1] = -1; 
01052                 mat_val[2] = 1;
01053                 lpx_set_mat_row(lp_glpk, mat_row_index+1, 3, mat_ind, mat_val);
01054                 lpx_set_row_bnds(lp_glpk, mat_row_index+1, LPX_UP, 0, 0);
01055             }
01056         }
01057     }
01058 
01059     int ind[num_kernels+2];
01060     float64_t val[num_kernels+2];
01061     int row_index = lpx_add_rows(lp_glpk, 1);
01062     for (int32_t i=1; i<=num_kernels; i++)
01063     {
01064         ind[i] = i;
01065         val[i] = -(sumw[i-1]-suma);
01066     }
01067     ind[num_kernels+1] = 2*num_kernels+1;
01068     val[num_kernels+1] = -1;
01069     lpx_set_mat_row(lp_glpk, row_index, num_kernels+1, ind, val);
01070     lpx_set_row_bnds(lp_glpk, row_index, LPX_UP, 0, 0);
01071 
01072     //optimize
01073     lpx_simplex(lp_glpk);
01074     bool res = check_lpx_status(lp_glpk);
01075     if (!res)
01076         SG_ERROR( "Failed to optimize Problem.\n");
01077 
01078     int32_t cur_numrows = lpx_get_num_rows(lp_glpk);
01079     int32_t cur_numcols = lpx_get_num_cols(lp_glpk);
01080     int32_t num_rows=cur_numrows;
01081     ASSERT(cur_numcols<=2*num_kernels+1);
01082 
01083     float64_t* col_primal = new float64_t[cur_numcols];
01084     float64_t* row_primal = new float64_t[cur_numrows];
01085     float64_t* row_dual = new float64_t[cur_numrows];
01086 
01087     for (int i=0; i<cur_numrows; i++)
01088     {
01089         row_primal[i] = lpx_get_row_prim(lp_glpk, i+1);
01090         row_dual[i] = lpx_get_row_dual(lp_glpk, i+1);
01091     }
01092     for (int i=0; i<cur_numcols; i++)
01093         col_primal[i] = lpx_get_col_prim(lp_glpk, i+1);
01094 
01095     obj = -col_primal[2*num_kernels];
01096 
01097     for (int i=0; i<num_kernels; i++)
01098         beta[i] = col_primal[i];
01099 
01100     int32_t num_active_rows=0;
01101     if(res)
01102     {
01103         float64_t max_slack = CMath::INFTY;
01104         int32_t max_idx = -1;
01105         int32_t start_row = 1;
01106         if (C_mkl!=0.0)
01107             start_row += 2*(num_kernels-1);
01108 
01109         for (int32_t i= start_row; i<cur_numrows; i++)
01110         {
01111             if (row_dual[i]!=0)
01112                 num_active_rows++;
01113             else
01114             {
01115                 if (row_primal[i]<max_slack)
01116                 {
01117                     max_slack = row_primal[i];
01118                     max_idx = i;
01119                 }
01120             }
01121         }
01122 
01123         if ((num_rows-start_row>CMath::max(100, 2*num_active_rows)) && max_idx!=-1)
01124         {
01125             int del_rows[2];
01126             del_rows[1] = max_idx+1;
01127             lpx_del_rows(lp_glpk, 1, del_rows);
01128         }
01129     }
01130 
01131     delete[] row_dual;
01132     delete[] row_primal;
01133     delete[] col_primal;
01134 #else
01135     SG_ERROR("Glpk not enabled at compile time\n");
01136 #endif
01137 
01138     return obj;
01139 }
01140 
01141 void CMKL::compute_sum_beta(float64_t* sumw)
01142 {
01143     ASSERT(sumw);
01144     ASSERT(svm);
01145 
01146     int32_t nsv=svm->get_num_support_vectors();
01147     int32_t num_kernels = kernel->get_num_subkernels();
01148     float64_t* beta = new float64_t[num_kernels];
01149     int32_t nweights=0;
01150     const float64_t* old_beta = kernel->get_subkernel_weights(nweights);
01151     ASSERT(nweights==num_kernels);
01152     ASSERT(old_beta);
01153 
01154     for (int32_t i=0; i<num_kernels; i++)
01155     {
01156         beta[i]=0;
01157         sumw[i]=0;
01158     }
01159 
01160     for (int32_t n=0; n<num_kernels; n++)
01161     {
01162         beta[n]=1.0;
01163         kernel->set_subkernel_weights(beta, num_kernels);
01164 
01165         for (int32_t i=0; i<nsv; i++)
01166         {   
01167             int32_t ii=svm->get_support_vector(i);
01168 
01169             for (int32_t j=0; j<nsv; j++)
01170             {   
01171                 int32_t jj=svm->get_support_vector(j);
01172                 sumw[n]+=0.5*svm->get_alpha(i)*svm->get_alpha(j)*kernel->kernel(ii,jj);
01173             }
01174         }
01175         beta[n]=0.0;
01176     }
01177 
01178     mkl_iterations++;
01179     kernel->set_subkernel_weights( (float64_t*) old_beta, num_kernels);
01180 }
01181 
01182 
01183 // assumes that all constraints are satisfied
01184 float64_t CMKL::compute_mkl_dual_objective()
01185 {
01186     int32_t n=get_num_support_vectors();
01187     float64_t mkl_obj=0;
01188 
01189     if (labels && kernel && kernel->get_kernel_type() == K_COMBINED)
01190     {
01191         CKernel* kn = ((CCombinedKernel*)kernel)->get_first_kernel();
01192         while (kn)
01193         {
01194             float64_t sum=0;
01195             for (int32_t i=0; i<n; i++)
01196             {
01197                 int32_t ii=get_support_vector(i);
01198 
01199                 for (int32_t j=0; j<n; j++)
01200                 {
01201                     int32_t jj=get_support_vector(j);
01202                     sum+=get_alpha(i)*get_alpha(j)*kn->kernel(ii,jj);
01203                 }
01204             }
01205 
01206             if (mkl_norm==1.0)
01207                 mkl_obj = CMath::max(mkl_obj, sum);
01208             else
01209                 mkl_obj += CMath::pow(sum, mkl_norm/(mkl_norm-1));
01210 
01211             SG_UNREF(kn);
01212             kn = ((CCombinedKernel*) kernel)->get_next_kernel();
01213         }
01214 
01215         if (mkl_norm==1.0)
01216             mkl_obj=-0.5*mkl_obj;
01217         else
01218             mkl_obj= -0.5*CMath::pow(mkl_obj, (mkl_norm-1)/mkl_norm);
01219 
01220         mkl_obj+=compute_sum_alpha();
01221     }
01222     else
01223         SG_ERROR( "cannot compute objective, labels or kernel not set\n");
01224 
01225     return -mkl_obj;
01226 }
01227 
01228 #ifdef USE_CPLEX
01229 void CMKL::set_qnorm_constraints(float64_t* beta, int32_t num_kernels)
01230 {
01231     ASSERT(num_kernels>0);
01232 
01233     float64_t* grad_beta=new float64_t[num_kernels];
01234     float64_t* hess_beta=new float64_t[num_kernels+1];
01235     float64_t* lin_term=new float64_t[num_kernels+1];
01236     int* ind=new int[num_kernels+1];
01237 
01238     //CMath::display_vector(beta, num_kernels, "beta");
01239     double const_term = 1-CMath::qsq(beta, num_kernels, mkl_norm);
01240     //SG_PRINT("const=%f\n", const_term);
01241     ASSERT(CMath::fequal(const_term, 0.0));
01242 
01243     for (int32_t i=0; i<num_kernels; i++)
01244     {
01245         grad_beta[i]=mkl_norm * pow(beta[i], mkl_norm-1);
01246         hess_beta[i]=0.5*mkl_norm*(mkl_norm-1) * pow(beta[i], mkl_norm-2); 
01247         lin_term[i]=grad_beta[i] - 2*beta[i]*hess_beta[i];
01248         const_term+=grad_beta[i]*beta[i] - CMath::sq(beta[i])*hess_beta[i];
01249         ind[i]=i;
01250     }
01251     ind[num_kernels]=2*num_kernels;
01252     hess_beta[num_kernels]=0;
01253     lin_term[num_kernels]=0;
01254 
01255     int status=0;
01256     int num=CPXgetnumqconstrs (env, lp_cplex);
01257 
01258     if (num>0)
01259     {
01260         status = CPXdelqconstrs (env, lp_cplex, 0, 0);
01261         ASSERT(!status);
01262     }
01263 
01264     status = CPXaddqconstr (env, lp_cplex, num_kernels+1, num_kernels+1, const_term, 'L', ind, lin_term,
01265             ind, ind, hess_beta, NULL);
01266     ASSERT(!status);
01267 
01268     //CPXwriteprob (env, lp_cplex, "prob.lp", NULL);
01269     //CPXqpwrite (env, lp_cplex, "prob.qp");
01270 
01271     delete[] grad_beta;
01272     delete[] hess_beta;
01273     delete[] lin_term;
01274     delete[] ind;
01275 }
01276 #endif // USE_CPLEX

SHOGUN Machine Learning Toolbox - Documentation