SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MKL.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2004-2009 Soeren Sonnenburg, Gunnar Raetsch
8  * Alexander Zien, Marius Kloft, Chen Guohua
9  * Copyright (C) 2009 Fraunhofer Institute FIRST and Max-Planck-Society
10  * Copyright (C) 2010 Ryota Tomioka (University of Tokyo)
11  */
12 
13 #include <list>
14 #include <shogun/lib/Signal.h>
18 
19 using namespace shogun;
20 
22  : CSVM(), svm(NULL), C_mkl(0), mkl_norm(1), ent_lambda(0), beta_local(NULL),
23  mkl_iterations(0), mkl_epsilon(1e-5), interleaved_optimization(true),
24  w_gap(1.0), rho(0)
25 {
27 #ifdef USE_CPLEX
28  lp_cplex = NULL ;
29  env = NULL ;
30 #endif
31 
32 #ifdef USE_GLPK
33  lp_glpk = NULL;
34 #endif
35 
36  SG_DEBUG("creating MKL object %p\n", this);
37  lp_initialized = false ;
38 }
39 
41 {
42  // -- Delete beta_local for ElasticnetMKL
43  delete [] beta_local;
44  beta_local = NULL;
45 
46  SG_DEBUG("deleting MKL object %p\n", this);
47  if (svm)
48  svm->set_callback_function(NULL, NULL);
49  SG_UNREF(svm);
50 }
51 
53 {
54 #ifdef USE_CPLEX
55  cleanup_cplex();
56 
58  init_cplex();
59 #endif
60 
61 #ifdef USE_GLPK
62  cleanup_glpk();
63 
64  if (get_solver_type() == ST_GLPK)
65  init_glpk();
66 #endif
67 }
68 
69 #ifdef USE_CPLEX
71 {
72  while (env==NULL)
73  {
74  SG_INFO( "trying to initialize CPLEX\n") ;
75 
76  int status = 0; // calling external lib
77  env = CPXopenCPLEX (&status);
78 
79  if ( env == NULL )
80  {
81  char errmsg[1024];
82  SG_WARNING( "Could not open CPLEX environment.\n");
83  CPXgeterrorstring (env, status, errmsg);
84  SG_WARNING( "%s", errmsg);
85  SG_WARNING( "retrying in 60 seconds\n");
86  sleep(60);
87  }
88  else
89  {
90  // select dual simplex based optimization
91  status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, CPX_ALG_DUAL);
92  if ( status )
93  {
94  SG_ERROR( "Failure to select dual lp optimization, error %d.\n", status);
95  }
96  else
97  {
98  status = CPXsetintparam (env, CPX_PARAM_DATACHECK, CPX_ON);
99  if ( status )
100  {
101  SG_ERROR( "Failure to turn on data checking, error %d.\n", status);
102  }
103  else
104  {
105  lp_cplex = CPXcreateprob (env, &status, "light");
106 
107  if ( lp_cplex == NULL )
108  SG_ERROR( "Failed to create LP.\n");
109  else
110  CPXchgobjsen (env, lp_cplex, CPX_MIN); /* Problem is minimization */
111  }
112  }
113  }
114  }
115 
116  return (lp_cplex != NULL) && (env != NULL);
117 }
118 
120 {
121  bool result=false;
122 
123  if (lp_cplex)
124  {
125  int32_t status = CPXfreeprob(env, &lp_cplex);
126  lp_cplex = NULL;
127  lp_initialized = false;
128 
129  if (status)
130  SG_WARNING( "CPXfreeprob failed, error code %d.\n", status);
131  else
132  result = true;
133  }
134 
135  if (env)
136  {
137  int32_t status = CPXcloseCPLEX (&env);
138  env=NULL;
139 
140  if (status)
141  {
142  char errmsg[1024];
143  SG_WARNING( "Could not close CPLEX environment.\n");
144  CPXgeterrorstring (env, status, errmsg);
145  SG_WARNING( "%s", errmsg);
146  }
147  else
148  result = true;
149  }
150  return result;
151 }
152 #endif
153 
154 #ifdef USE_GLPK
156 {
157  lp_glpk = lpx_create_prob();
158  lpx_set_obj_dir(lp_glpk, LPX_MIN);
159  lpx_set_int_parm(lp_glpk, LPX_K_DUAL, GLP_ON );
160  lpx_set_int_parm(lp_glpk, LPX_K_PRESOL, GLP_ON );
161 
162  glp_term_out(GLP_OFF);
163  return (lp_glpk != NULL);
164 }
165 
167 {
168  lp_initialized = false;
169  if (lp_glpk)
170  lpx_delete_prob(lp_glpk);
171  lp_glpk = NULL;
172  return true;
173 }
174 
176 {
177  int status = lpx_get_status(lp);
178 
179  if (status==LPX_INFEAS)
180  {
181  SG_PRINT("solution is infeasible!\n");
182  return false;
183  }
184  else if(status==LPX_NOFEAS)
185  {
186  SG_PRINT("problem has no feasible solution!\n");
187  return false;
188  }
189  return true;
190 }
191 #endif // USE_GLPK
192 
194 {
195  ASSERT(kernel);
197 
198  if (data)
199  {
200  if (labels->get_num_labels() != data->get_num_vectors())
201  SG_ERROR("Number of training vectors does not match number of labels\n");
202  kernel->init(data, data);
203  }
204 
205  init_training();
206  if (!svm)
207  SG_ERROR("No constraint generator (SVM) set\n");
208 
209  int32_t num_label=0;
210  if (labels)
211  num_label = labels->get_num_labels();
212 
213  SG_INFO("%d trainlabels (%ld)\n", num_label, labels);
214  if (mkl_epsilon<=0)
215  mkl_epsilon=1e-2 ;
216 
217  SG_INFO("mkl_epsilon = %1.1e\n", mkl_epsilon) ;
218  SG_INFO("C_mkl = %1.1e\n", C_mkl) ;
219  SG_INFO("mkl_norm = %1.3e\n", mkl_norm);
220  SG_INFO("solver = %d\n", get_solver_type());
221  SG_INFO("ent_lambda = %f\n", ent_lambda);
222  SG_INFO("mkl_block_norm = %f\n", mkl_block_norm);
223 
224  int32_t num_weights = -1;
225  int32_t num_kernels = kernel->get_num_subkernels();
226  SG_INFO("num_kernels = %d\n", num_kernels);
227  const float64_t* beta_const = kernel->get_subkernel_weights(num_weights);
228  float64_t* beta = CMath::clone_vector(beta_const, num_weights);
229  ASSERT(num_weights==num_kernels);
230 
232  mkl_block_norm>=1 &&
233  mkl_block_norm<=2)
234  {
236  SG_WARNING("Switching to ST_DIRECT method with mkl_norm=%g\n", mkl_norm);
238  }
239 
241  {
242  // -- Initialize subkernel weights for Elasticnet MKL
243  CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, 1.0), beta, num_kernels);
244 
245  if (beta_local) { delete [] beta_local; }
246  beta_local = CMath::clone_vector(beta, num_kernels);
247 
248  elasticnet_transform(beta, ent_lambda, num_kernels);
249  }
250  else
251  {
252  CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm),
253  beta, num_kernels); //q-norm = 1
254  }
255 
256  kernel->set_subkernel_weights(SGVector<float64_t>(beta, num_kernels));
257  SG_FREE(beta);
258 
262  svm->set_nu(get_nu());
263  svm->set_C(get_C1(), get_C2());
270 
271 #ifdef USE_CPLEX
272  cleanup_cplex();
273 
274  if (get_solver_type()==ST_CPLEX)
275  init_cplex();
276 #endif
277 
278 #ifdef USE_GLPK
279  if (get_solver_type()==ST_GLPK)
280  init_glpk();
281 #endif
282 
283  mkl_iterations = 0;
285 
287 
289  {
291  {
292  SG_ERROR("Interleaved MKL optimization is currently "
293  "only supported with SVMlight\n");
294  }
295 
296  //Note that there is currently no safe way to properly resolve this
297  //situation: the mkl object hands itself as a reference to the svm which
298  //in turn increases the reference count of mkl and when done decreases
299  //the count. Thus we have to protect the mkl object from deletion by
300  //ref()'ing it before we set the callback function and should also
301  //unref() it afterwards. However, when the reference count was zero
302  //before this unref() might actually try to destroy this (crash ahead)
303  //but if we don't actually unref() the object we might leak memory...
304  //So as a workaround we only unref when the reference count was >1
305  //before.
306 #ifdef USE_REFERENCE_COUNTING
307  int32_t refs=this->ref();
308 #endif
310  svm->train();
311  SG_DONE();
312  svm->set_callback_function(NULL, NULL);
313 #ifdef USE_REFERENCE_COUNTING
314  if (refs>1)
315  this->unref();
316 #endif
317  }
318  else
319  {
320  float64_t* sumw = SG_MALLOC(float64_t, num_kernels);
321 
322 
323 
324  while (true)
325  {
326  svm->train();
327 
329  compute_sum_beta(sumw);
330 
332  {
333  SG_SWARNING("MKL Algorithm terminates PREMATURELY due to current training time exceeding get_max_train_time ()= %f . It may have not converged yet!\n",get_max_train_time ());
334  break;
335  }
336 
337 
338  mkl_iterations++;
339  if (perform_mkl_step(sumw, suma) || CSignal::cancel_computations())
340  break;
341  }
342 
343  SG_FREE(sumw);
344  }
345 #ifdef USE_CPLEX
346  cleanup_cplex();
347 #endif
348 #ifdef USE_GLPK
349  cleanup_glpk();
350 #endif
351 
352  int32_t nsv=svm->get_num_support_vectors();
353  create_new_model(nsv);
354 
355  set_bias(svm->get_bias());
356  for (int32_t i=0; i<nsv; i++)
357  {
358  set_alpha(i, svm->get_alpha(i));
360  }
361  return true;
362 }
363 
364 
366 {
367 
368  if (norm<1)
369  SG_ERROR("Norm must be >= 1, e.g., 1-norm is the standard MKL; norms>1 nonsparse MKL\n");
370 
371  mkl_norm = norm;
372 }
373 
375 {
376  if (lambda>1 || lambda<0)
377  SG_ERROR("0<=lambda<=1\n");
378 
379  if (lambda==0)
380  lambda = 1e-6;
381  else if (lambda==1.0)
382  lambda = 1.0-1e-6;
383 
384  ent_lambda=lambda;
385 }
386 
388 {
389  if (q<1)
390  SG_ERROR("1<=q<=inf\n");
391 
392  mkl_block_norm=q;
393 }
394 
396  const float64_t* sumw, float64_t suma)
397 {
399  {
400  SG_SWARNING("MKL Algorithm terminates PREMATURELY due to current training time exceeding get_max_train_time ()= %f . It may have not converged yet!\n",get_max_train_time ());
401  return true;
402  }
403 
404  int32_t num_kernels = kernel->get_num_subkernels();
405  int32_t nweights=0;
406  const float64_t* old_beta = kernel->get_subkernel_weights(nweights);
407  ASSERT(nweights==num_kernels);
408  float64_t* beta = SG_MALLOC(float64_t, num_kernels);
409 
410  int32_t inner_iters=0;
411  float64_t mkl_objective=0;
412 
413  mkl_objective=-suma;
414  for (int32_t i=0; i<num_kernels; i++)
415  {
416  beta[i]=old_beta[i];
417  mkl_objective+=old_beta[i]*sumw[i];
418  }
419 
420  w_gap = CMath::abs(1-rho/mkl_objective) ;
421 
422  if( (w_gap >= mkl_epsilon) ||
424  {
426  {
427  rho=compute_optimal_betas_directly(beta, old_beta, num_kernels, sumw, suma, mkl_objective);
428  }
429  else if (get_solver_type()==ST_BLOCK_NORM)
430  {
431  rho=compute_optimal_betas_block_norm(beta, old_beta, num_kernels, sumw, suma, mkl_objective);
432  }
433  else if (get_solver_type()==ST_ELASTICNET)
434  {
435  // -- Direct update of subkernel weights for ElasticnetMKL
436  // Note that ElasticnetMKL solves a different optimization
437  // problem from the rest of the solver types
438  rho=compute_optimal_betas_elasticnet(beta, old_beta, num_kernels, sumw, suma, mkl_objective);
439  }
440  else if (get_solver_type()==ST_NEWTON)
441  rho=compute_optimal_betas_newton(beta, old_beta, num_kernels, sumw, suma, mkl_objective);
442 #ifdef USE_CPLEX
443  else if (get_solver_type()==ST_CPLEX)
444  rho=compute_optimal_betas_via_cplex(beta, old_beta, num_kernels, sumw, suma, inner_iters);
445 #endif
446 #ifdef USE_GLPK
447  else if (get_solver_type()==ST_GLPK)
448  rho=compute_optimal_betas_via_glpk(beta, old_beta, num_kernels, sumw, suma, inner_iters);
449 #endif
450  else
451  SG_ERROR("Solver type not supported (not compiled in?)\n");
452 
453  w_gap = CMath::abs(1-rho/mkl_objective) ;
454  }
455 
456  kernel->set_subkernel_weights(SGVector<float64_t>(beta, num_kernels));
457  SG_FREE(beta);
458 
459  return converged();
460 }
461 
463  float64_t* beta, const float64_t* old_beta, const int32_t num_kernels,
464  const float64_t* sumw, const float64_t suma,
465  const float64_t mkl_objective )
466 {
467  const float64_t epsRegul = 0.01; // fraction of root mean squared deviation
468  float64_t obj;
469  float64_t Z;
470  float64_t preR;
471  int32_t p;
472  int32_t nofKernelsGood;
473 
474  // --- optimal beta
475  nofKernelsGood = num_kernels;
476 
477  Z = 0;
478  for (p=0; p<num_kernels; ++p )
479  {
480  if (sumw[p] >= 0.0 && old_beta[p] >= 0.0 )
481  {
482  beta[p] = CMath::sqrt(sumw[p]*old_beta[p]*old_beta[p]);
483  Z += beta[p];
484  }
485  else
486  {
487  beta[p] = 0.0;
488  --nofKernelsGood;
489  }
490  ASSERT( beta[p] >= 0 );
491  }
492 
493  // --- normalize
494  CMath::scale_vector(1.0/Z, beta, num_kernels);
495 
496  // --- regularize & renormalize
497 
498  preR = 0.0;
499  for( p=0; p<num_kernels; ++p )
500  preR += CMath::pow( beta_local[p] - beta[p], 2.0 );
501  const float64_t R = CMath::sqrt( preR ) * epsRegul;
502  if( !( R >= 0 ) )
503  {
504  SG_PRINT( "MKL-direct: p = %.3f\n", 1.0 );
505  SG_PRINT( "MKL-direct: nofKernelsGood = %d\n", nofKernelsGood );
506  SG_PRINT( "MKL-direct: Z = %e\n", Z );
507  SG_PRINT( "MKL-direct: eps = %e\n", epsRegul );
508  for( p=0; p<num_kernels; ++p )
509  {
510  const float64_t t = CMath::pow( beta_local[p] - beta[p], 2.0 );
511  SG_PRINT( "MKL-direct: t[%3d] = %e ( diff = %e = %e - %e )\n", p, t, beta_local[p]-beta[p], beta_local[p], beta[p] );
512  }
513  SG_PRINT( "MKL-direct: preR = %e\n", preR );
514  SG_PRINT( "MKL-direct: preR/p = %e\n", preR );
515  SG_PRINT( "MKL-direct: sqrt(preR/p) = %e\n", CMath::sqrt(preR) );
516  SG_PRINT( "MKL-direct: R = %e\n", R );
517  SG_ERROR( "Assertion R >= 0 failed!\n" );
518  }
519 
520  Z = 0.0;
521  for( p=0; p<num_kernels; ++p )
522  {
523  beta[p] += R;
524  Z += beta[p];
525  ASSERT( beta[p] >= 0 );
526  }
527  Z = CMath::pow( Z, -1.0 );
528  ASSERT( Z >= 0 );
529  for( p=0; p<num_kernels; ++p )
530  {
531  beta[p] *= Z;
532  ASSERT( beta[p] >= 0.0 );
533  if (beta[p] > 1.0 )
534  beta[p] = 1.0;
535  }
536 
537  // --- copy beta into beta_local
538  for( p=0; p<num_kernels; ++p )
539  beta_local[p] = beta[p];
540 
541  // --- elastic-net transform
542  elasticnet_transform(beta, ent_lambda, num_kernels);
543 
544  // --- objective
545  obj = -suma;
546  for (p=0; p<num_kernels; ++p )
547  {
548  //obj += sumw[p] * old_beta[p]*old_beta[p] / beta[p];
549  obj += sumw[p] * beta[p];
550  }
551  return obj;
552 }
553 
555  const float64_t &del, const float64_t* nm, int32_t len,
556  const float64_t &lambda)
557 {
558  std::list<int32_t> I;
559  float64_t gam = 1.0-lambda;
560  for (int32_t i=0; i<len;i++)
561  {
562  if (nm[i]>=CMath::sqrt(2*del*gam))
563  I.push_back(i);
564  }
565  int32_t M=I.size();
566 
567  *ff=del;
568  *gg=-(M*gam+lambda);
569  *hh=0;
570  for (std::list<int32_t>::iterator it=I.begin(); it!=I.end(); it++)
571  {
572  float64_t nmit = nm[*it];
573 
574  *ff += del*gam*CMath::pow(nmit/CMath::sqrt(2*del*gam)-1,2)/lambda;
575  *gg += CMath::sqrt(gam/(2*del))*nmit;
576  *hh += -0.5*CMath::sqrt(gam/(2*CMath::pow(del,3)))*nmit;
577  }
578 }
579 
580 // assumes that all constraints are satisfied
582 {
583  int32_t n=get_num_support_vectors();
584  int32_t num_kernels = kernel->get_num_subkernels();
585  float64_t mkl_obj=0;
586 
588  {
589  // Compute Elastic-net dual
590  float64_t* nm = SG_MALLOC(float64_t, num_kernels);
591  float64_t del=0;
592  CKernel* kn = ((CCombinedKernel*)kernel)->get_first_kernel();
593 
594  int32_t k=0;
595  while (kn)
596  {
597  float64_t sum=0;
598  for (int32_t i=0; i<n; i++)
599  {
600  int32_t ii=get_support_vector(i);
601 
602  for (int32_t j=0; j<n; j++)
603  {
604  int32_t jj=get_support_vector(j);
605  sum+=get_alpha(i)*get_alpha(j)*kn->kernel(ii,jj);
606  }
607  }
608  nm[k]= CMath::pow(sum, 0.5);
609  del = CMath::max(del, nm[k]);
610 
611  // SG_PRINT("nm[%d]=%f\n",k,nm[k]);
612  k++;
613 
614 
615  SG_UNREF(kn);
616  kn = ((CCombinedKernel*) kernel)->get_next_kernel();
617  }
618  // initial delta
619  del = del/CMath::sqrt(2*(1-ent_lambda));
620 
621  // Newton's method to optimize delta
622  k=0;
623  float64_t ff, gg, hh;
624  elasticnet_dual(&ff, &gg, &hh, del, nm, num_kernels, ent_lambda);
625  while (CMath::abs(gg)>+1e-8 && k<100)
626  {
627  float64_t ff_old = ff;
628  float64_t gg_old = gg;
629  float64_t del_old = del;
630 
631  // SG_PRINT("[%d] fval=%f gg=%f del=%f\n", k, ff, gg, del);
632 
633  float64_t step=1.0;
634  do
635  {
636  del=del_old*CMath::exp(-step*gg/(hh*del+gg));
637  elasticnet_dual(&ff, &gg, &hh, del, nm, num_kernels, ent_lambda);
638  step/=2;
639  } while(ff>ff_old+1e-4*gg_old*(del-del_old));
640 
641  k++;
642  }
643  mkl_obj=-ff;
644 
645  SG_FREE(nm);
646 
647  mkl_obj+=compute_sum_alpha();
648 
649  }
650  else
651  SG_ERROR( "cannot compute objective, labels or kernel not set\n");
652 
653  return -mkl_obj;
654 }
655 
657  float64_t* beta, const float64_t* old_beta, const int32_t num_kernels,
658  const float64_t* sumw, const float64_t suma,
659  const float64_t mkl_objective )
660 {
661  float64_t obj;
662  float64_t Z=0;
663  int32_t p;
664 
665  // --- optimal beta
666  for( p=0; p<num_kernels; ++p )
667  {
668  ASSERT(sumw[p]>=0);
669 
670  beta[p] = CMath::pow( sumw[p], -(2.0-mkl_block_norm)/(2.0-2.0*mkl_block_norm));
671  Z+= CMath::pow( sumw[p], -(mkl_block_norm)/(2.0-2.0*mkl_block_norm));
672 
673  ASSERT( beta[p] >= 0 );
674  }
675 
676  ASSERT(Z>=0);
677 
678  Z=1.0/CMath::pow(Z, (2.0-mkl_block_norm)/mkl_block_norm);
679 
680  for( p=0; p<num_kernels; ++p )
681  beta[p] *= Z;
682 
683  // --- objective
684  obj = -suma;
685  for( p=0; p<num_kernels; ++p )
686  obj += sumw[p] * beta[p];
687 
688  return obj;
689 }
690 
691 
693  float64_t* beta, const float64_t* old_beta, const int32_t num_kernels,
694  const float64_t* sumw, const float64_t suma,
695  const float64_t mkl_objective )
696 {
697  const float64_t epsRegul = 0.01; // fraction of root mean squared deviation
698  float64_t obj;
699  float64_t Z;
700  float64_t preR;
701  int32_t p;
702  int32_t nofKernelsGood;
703 
704  // --- optimal beta
705  nofKernelsGood = num_kernels;
706  for( p=0; p<num_kernels; ++p )
707  {
708  //SG_PRINT( "MKL-direct: sumw[%3d] = %e ( oldbeta = %e )\n", p, sumw[p], old_beta[p] );
709  if( sumw[p] >= 0.0 && old_beta[p] >= 0.0 )
710  {
711  beta[p] = sumw[p] * old_beta[p]*old_beta[p] / mkl_norm;
712  beta[p] = CMath::pow( beta[p], 1.0 / (mkl_norm+1.0) );
713  }
714  else
715  {
716  beta[p] = 0.0;
717  --nofKernelsGood;
718  }
719  ASSERT( beta[p] >= 0 );
720  }
721 
722  // --- normalize
723  Z = 0.0;
724  for( p=0; p<num_kernels; ++p )
725  Z += CMath::pow( beta[p], mkl_norm );
726 
727  Z = CMath::pow( Z, -1.0/mkl_norm );
728  ASSERT( Z >= 0 );
729  for( p=0; p<num_kernels; ++p )
730  beta[p] *= Z;
731 
732  // --- regularize & renormalize
733  preR = 0.0;
734  for( p=0; p<num_kernels; ++p )
735  preR += CMath::sq( old_beta[p] - beta[p]);
736 
737  const float64_t R = CMath::sqrt( preR / mkl_norm ) * epsRegul;
738  if( !( R >= 0 ) )
739  {
740  SG_PRINT( "MKL-direct: p = %.3f\n", mkl_norm );
741  SG_PRINT( "MKL-direct: nofKernelsGood = %d\n", nofKernelsGood );
742  SG_PRINT( "MKL-direct: Z = %e\n", Z );
743  SG_PRINT( "MKL-direct: eps = %e\n", epsRegul );
744  for( p=0; p<num_kernels; ++p )
745  {
746  const float64_t t = CMath::pow( old_beta[p] - beta[p], 2.0 );
747  SG_PRINT( "MKL-direct: t[%3d] = %e ( diff = %e = %e - %e )\n", p, t, old_beta[p]-beta[p], old_beta[p], beta[p] );
748  }
749  SG_PRINT( "MKL-direct: preR = %e\n", preR );
750  SG_PRINT( "MKL-direct: preR/p = %e\n", preR/mkl_norm );
751  SG_PRINT( "MKL-direct: sqrt(preR/p) = %e\n", CMath::sqrt(preR/mkl_norm) );
752  SG_PRINT( "MKL-direct: R = %e\n", R );
753  SG_ERROR( "Assertion R >= 0 failed!\n" );
754  }
755 
756  Z = 0.0;
757  for( p=0; p<num_kernels; ++p )
758  {
759  beta[p] += R;
760  Z += CMath::pow( beta[p], mkl_norm );
761  ASSERT( beta[p] >= 0 );
762  }
763  Z = CMath::pow( Z, -1.0/mkl_norm );
764  ASSERT( Z >= 0 );
765  for( p=0; p<num_kernels; ++p )
766  {
767  beta[p] *= Z;
768  ASSERT( beta[p] >= 0.0 );
769  if( beta[p] > 1.0 )
770  beta[p] = 1.0;
771  }
772 
773  // --- objective
774  obj = -suma;
775  for( p=0; p<num_kernels; ++p )
776  obj += sumw[p] * beta[p];
777 
778  return obj;
779 }
780 
782  const float64_t* old_beta, int32_t num_kernels,
783  const float64_t* sumw, float64_t suma,
784  float64_t mkl_objective)
785 {
786  SG_DEBUG("MKL via NEWTON\n");
787 
788  if (mkl_norm==1)
789  SG_ERROR("MKL via NEWTON works only for norms>1\n");
790 
791  const double epsBeta = 1e-32;
792  const double epsGamma = 1e-12;
793  const double epsWsq = 1e-12;
794  const double epsNewt = 0.0001;
795  const double epsStep = 1e-9;
796  const int nofNewtonSteps = 3;
797  const double hessRidge = 1e-6;
798  const int inLogSpace = 0;
799 
800  const float64_t r = mkl_norm / ( mkl_norm - 1.0 );
801  float64_t* newtDir = SG_MALLOC(float64_t, num_kernels );
802  float64_t* newtBeta = SG_MALLOC(float64_t, num_kernels );
803  //float64_t newtStep;
804  float64_t stepSize;
805  float64_t Z;
806  float64_t obj;
807  float64_t gamma;
808  int32_t p;
809  int i;
810 
811  // --- check beta
812  Z = 0.0;
813  for( p=0; p<num_kernels; ++p )
814  {
815  beta[p] = old_beta[p];
816  if( !( beta[p] >= epsBeta ) )
817  beta[p] = epsBeta;
818 
819  ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
820  Z += CMath::pow( beta[p], mkl_norm );
821  }
822 
823  Z = CMath::pow( Z, -1.0/mkl_norm );
824  if( !( fabs(Z-1.0) <= epsGamma ) )
825  {
826  SG_WARNING( "old_beta not normalized (diff=%e); forcing normalization. ", Z-1.0 );
827  for( p=0; p<num_kernels; ++p )
828  {
829  beta[p] *= Z;
830  if( beta[p] > 1.0 )
831  beta[p] = 1.0;
832  ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
833  }
834  }
835 
836  // --- compute gamma
837  gamma = 0.0;
838  for ( p=0; p<num_kernels; ++p )
839  {
840  if ( !( sumw[p] >= 0 ) )
841  {
842  if( !( sumw[p] >= -epsWsq ) )
843  SG_WARNING( "sumw[%d] = %e; treated as 0. ", p, sumw[p] );
844  // should better recompute sumw[] !!!
845  }
846  else
847  {
848  ASSERT( sumw[p] >= 0 );
849  //gamma += CMath::pow( sumw[p] * beta[p]*beta[p], r );
850  gamma += CMath::pow( sumw[p] * beta[p]*beta[p] / mkl_norm, r );
851  }
852  }
853  gamma = CMath::pow( gamma, 1.0/r ) / mkl_norm;
854  ASSERT( gamma > -1e-9 );
855  if( !( gamma > epsGamma ) )
856  {
857  SG_WARNING( "bad gamma: %e; set to %e. ", gamma, epsGamma );
858  // should better recompute sumw[] !!!
859  gamma = epsGamma;
860  }
861  ASSERT( gamma >= epsGamma );
862  //gamma = -gamma;
863 
864  // --- compute objective
865  obj = 0.0;
866  for( p=0; p<num_kernels; ++p )
867  {
868  obj += beta[p] * sumw[p];
869  //obj += gamma/mkl_norm * CMath::pow( beta[p], mkl_norm );
870  }
871  if( !( obj >= 0.0 ) )
872  SG_WARNING( "negative objective: %e. ", obj );
873  //SG_PRINT( "OBJ = %e. \n", obj );
874 
875 
876  // === perform Newton steps
877  for (i = 0; i < nofNewtonSteps; ++i )
878  {
879  // --- compute Newton direction (Hessian is diagonal)
880  const float64_t gqq1 = mkl_norm * (mkl_norm-1.0) * gamma;
881  // newtStep = 0.0;
882  for( p=0; p<num_kernels; ++p )
883  {
884  ASSERT( 0.0 <= beta[p] && beta[p] <= 1.0 );
885  //const float halfw2p = ( sumw[p] >= 0.0 ) ? sumw[p] : 0.0;
886  //const float64_t t1 = halfw2p*beta[p] - mkl_norm*gamma*CMath::pow(beta[p],mkl_norm);
887  //const float64_t t2 = 2.0*halfw2p + gqq1*CMath::pow(beta[p],mkl_norm-1.0);
888  const float halfw2p = ( sumw[p] >= 0.0 ) ? (sumw[p]*old_beta[p]*old_beta[p]) : 0.0;
889  const float64_t t0 = halfw2p*beta[p] - mkl_norm*gamma*CMath::pow(beta[p],mkl_norm+2.0);
890  const float64_t t1 = ( t0 < 0 ) ? 0.0 : t0;
891  const float64_t t2 = 2.0*halfw2p + gqq1*CMath::pow(beta[p],mkl_norm+1.0);
892  if( inLogSpace )
893  newtDir[p] = t1 / ( t1 + t2*beta[p] + hessRidge );
894  else
895  newtDir[p] = ( t1 == 0.0 ) ? 0.0 : ( t1 / t2 );
896  // newtStep += newtDir[p] * grad[p];
897  ASSERT( newtDir[p] == newtDir[p] );
898  //SG_PRINT( "newtDir[%d] = %6.3f = %e / %e \n", p, newtDir[p], t1, t2 );
899  }
900  //CMath::display_vector( newtDir, num_kernels, "newton direction " );
901  //SG_PRINT( "Newton step size = %e\n", Z );
902 
903  // --- line search
904  stepSize = 1.0;
905  while( stepSize >= epsStep )
906  {
907  // --- perform Newton step
908  Z = 0.0;
909  while( Z == 0.0 )
910  {
911  for( p=0; p<num_kernels; ++p )
912  {
913  if( inLogSpace )
914  newtBeta[p] = beta[p] * CMath::exp( + stepSize * newtDir[p] );
915  else
916  newtBeta[p] = beta[p] + stepSize * newtDir[p];
917  if( !( newtBeta[p] >= epsBeta ) )
918  newtBeta[p] = epsBeta;
919  Z += CMath::pow( newtBeta[p], mkl_norm );
920  }
921  ASSERT( 0.0 <= Z );
922  Z = CMath::pow( Z, -1.0/mkl_norm );
923  if( Z == 0.0 )
924  stepSize /= 2.0;
925  }
926 
927  // --- normalize new beta (wrt p-norm)
928  ASSERT( 0.0 < Z );
929  ASSERT( Z < CMath::INFTY );
930  for( p=0; p<num_kernels; ++p )
931  {
932  newtBeta[p] *= Z;
933  if( newtBeta[p] > 1.0 )
934  {
935  //SG_WARNING( "beta[%d] = %e; set to 1. ", p, beta[p] );
936  newtBeta[p] = 1.0;
937  }
938  ASSERT( 0.0 <= newtBeta[p] && newtBeta[p] <= 1.0 );
939  }
940 
941  // --- objective increased?
942  float64_t newtObj;
943  newtObj = 0.0;
944  for( p=0; p<num_kernels; ++p )
945  newtObj += sumw[p] * old_beta[p]*old_beta[p] / newtBeta[p];
946  //SG_PRINT( "step = %.8f: obj = %e -> %e. \n", stepSize, obj, newtObj );
947  if ( newtObj < obj - epsNewt*stepSize*obj )
948  {
949  for( p=0; p<num_kernels; ++p )
950  beta[p] = newtBeta[p];
951  obj = newtObj;
952  break;
953  }
954  stepSize /= 2.0;
955 
956  }
957 
958  if( stepSize < epsStep )
959  break;
960  }
961  SG_FREE(newtDir);
962  SG_FREE(newtBeta);
963 
964  // === return new objective
965  obj = -suma;
966  for( p=0; p<num_kernels; ++p )
967  obj += beta[p] * sumw[p];
968  return obj;
969 }
970 
971 
972 
973 float64_t CMKL::compute_optimal_betas_via_cplex(float64_t* new_beta, const float64_t* old_beta, int32_t num_kernels,
974  const float64_t* sumw, float64_t suma, int32_t& inner_iters)
975 {
976  SG_DEBUG("MKL via CPLEX\n");
977 
978 #ifdef USE_CPLEX
979  ASSERT(new_beta);
980  ASSERT(old_beta);
981 
982  int32_t NUMCOLS = 2*num_kernels + 1;
983  double* x=SG_MALLOC(double, NUMCOLS);
984 
985  if (!lp_initialized)
986  {
987  SG_INFO( "creating LP\n") ;
988 
989  double obj[NUMCOLS]; /* calling external lib */
990  double lb[NUMCOLS]; /* calling external lib */
991  double ub[NUMCOLS]; /* calling external lib */
992 
993  for (int32_t i=0; i<2*num_kernels; i++)
994  {
995  obj[i]=0 ;
996  lb[i]=0 ;
997  ub[i]=1 ;
998  }
999 
1000  for (int32_t i=num_kernels; i<2*num_kernels; i++)
1001  obj[i]= C_mkl;
1002 
1003  obj[2*num_kernels]=1 ;
1004  lb[2*num_kernels]=-CPX_INFBOUND ;
1005  ub[2*num_kernels]=CPX_INFBOUND ;
1006 
1007  int status = CPXnewcols (env, lp_cplex, NUMCOLS, obj, lb, ub, NULL, NULL);
1008  if ( status ) {
1009  char errmsg[1024];
1010  CPXgeterrorstring (env, status, errmsg);
1011  SG_ERROR( "%s", errmsg);
1012  }
1013 
1014  // add constraint sum(w)=1;
1015  SG_INFO( "adding the first row\n");
1016  int initial_rmatbeg[1]; /* calling external lib */
1017  int initial_rmatind[num_kernels+1]; /* calling external lib */
1018  double initial_rmatval[num_kernels+1]; /* calling ext lib */
1019  double initial_rhs[1]; /* calling external lib */
1020  char initial_sense[1];
1021 
1022  // 1-norm MKL
1023  if (mkl_norm==1)
1024  {
1025  initial_rmatbeg[0] = 0;
1026  initial_rhs[0]=1 ; // rhs=1 ;
1027  initial_sense[0]='E' ; // equality
1028 
1029  // sparse matrix
1030  for (int32_t i=0; i<num_kernels; i++)
1031  {
1032  initial_rmatind[i]=i ; //index of non-zero element
1033  initial_rmatval[i]=1 ; //value of ...
1034  }
1035  initial_rmatind[num_kernels]=2*num_kernels ; //number of non-zero elements
1036  initial_rmatval[num_kernels]=0 ;
1037 
1038  status = CPXaddrows (env, lp_cplex, 0, 1, num_kernels+1,
1039  initial_rhs, initial_sense, initial_rmatbeg,
1040  initial_rmatind, initial_rmatval, NULL, NULL);
1041 
1042  }
1043  else // 2 and q-norm MKL
1044  {
1045  initial_rmatbeg[0] = 0;
1046  initial_rhs[0]=0 ; // rhs=1 ;
1047  initial_sense[0]='L' ; // <= (inequality)
1048 
1049  initial_rmatind[0]=2*num_kernels ;
1050  initial_rmatval[0]=0 ;
1051 
1052  status = CPXaddrows (env, lp_cplex, 0, 1, 1,
1053  initial_rhs, initial_sense, initial_rmatbeg,
1054  initial_rmatind, initial_rmatval, NULL, NULL);
1055 
1056 
1057  if (mkl_norm==2)
1058  {
1059  for (int32_t i=0; i<num_kernels; i++)
1060  {
1061  initial_rmatind[i]=i ;
1062  initial_rmatval[i]=1 ;
1063  }
1064  initial_rmatind[num_kernels]=2*num_kernels ;
1065  initial_rmatval[num_kernels]=0 ;
1066 
1067  status = CPXaddqconstr (env, lp_cplex, 0, num_kernels+1, 1.0, 'L', NULL, NULL,
1068  initial_rmatind, initial_rmatind, initial_rmatval, NULL);
1069  }
1070  }
1071 
1072 
1073  if ( status )
1074  SG_ERROR( "Failed to add the first row.\n");
1075 
1076  lp_initialized = true ;
1077 
1078  if (C_mkl!=0.0)
1079  {
1080  for (int32_t q=0; q<num_kernels-1; q++)
1081  {
1082  // add constraint w[i]-w[i+1]<s[i];
1083  // add constraint w[i+1]-w[i]<s[i];
1084  int rmatbeg[1]; /* calling external lib */
1085  int rmatind[3]; /* calling external lib */
1086  double rmatval[3]; /* calling external lib */
1087  double rhs[1]; /* calling external lib */
1088  char sense[1];
1089 
1090  rmatbeg[0] = 0;
1091  rhs[0]=0 ; // rhs=0 ;
1092  sense[0]='L' ; // <=
1093  rmatind[0]=q ;
1094  rmatval[0]=1 ;
1095  rmatind[1]=q+1 ;
1096  rmatval[1]=-1 ;
1097  rmatind[2]=num_kernels+q ;
1098  rmatval[2]=-1 ;
1099  status = CPXaddrows (env, lp_cplex, 0, 1, 3,
1100  rhs, sense, rmatbeg,
1101  rmatind, rmatval, NULL, NULL);
1102  if ( status )
1103  SG_ERROR( "Failed to add a smothness row (1).\n");
1104 
1105  rmatbeg[0] = 0;
1106  rhs[0]=0 ; // rhs=0 ;
1107  sense[0]='L' ; // <=
1108  rmatind[0]=q ;
1109  rmatval[0]=-1 ;
1110  rmatind[1]=q+1 ;
1111  rmatval[1]=1 ;
1112  rmatind[2]=num_kernels+q ;
1113  rmatval[2]=-1 ;
1114  status = CPXaddrows (env, lp_cplex, 0, 1, 3,
1115  rhs, sense, rmatbeg,
1116  rmatind, rmatval, NULL, NULL);
1117  if ( status )
1118  SG_ERROR( "Failed to add a smothness row (2).\n");
1119  }
1120  }
1121  }
1122 
1123  { // add the new row
1124  //SG_INFO( "add the new row\n") ;
1125 
1126  int rmatbeg[1];
1127  int rmatind[num_kernels+1];
1128  double rmatval[num_kernels+1];
1129  double rhs[1];
1130  char sense[1];
1131 
1132  rmatbeg[0] = 0;
1133  if (mkl_norm==1)
1134  rhs[0]=0 ;
1135  else
1136  rhs[0]=-suma ;
1137 
1138  sense[0]='L' ;
1139 
1140  for (int32_t i=0; i<num_kernels; i++)
1141  {
1142  rmatind[i]=i ;
1143  if (mkl_norm==1)
1144  rmatval[i]=-(sumw[i]-suma) ;
1145  else
1146  rmatval[i]=-sumw[i];
1147  }
1148  rmatind[num_kernels]=2*num_kernels ;
1149  rmatval[num_kernels]=-1 ;
1150 
1151  int32_t status = CPXaddrows (env, lp_cplex, 0, 1, num_kernels+1,
1152  rhs, sense, rmatbeg,
1153  rmatind, rmatval, NULL, NULL);
1154  if ( status )
1155  SG_ERROR( "Failed to add the new row.\n");
1156  }
1157 
1158  inner_iters=0;
1159  int status;
1160  {
1161 
1162  if (mkl_norm==1) // optimize 1 norm MKL
1163  status = CPXlpopt (env, lp_cplex);
1164  else if (mkl_norm==2) // optimize 2-norm MKL
1165  status = CPXbaropt(env, lp_cplex);
1166  else // q-norm MKL
1167  {
1168  float64_t* beta=SG_MALLOC(float64_t, 2*num_kernels+1);
1169  float64_t objval_old=1e-8; //some value to cause the loop to not terminate yet
1170  for (int32_t i=0; i<num_kernels; i++)
1171  beta[i]=old_beta[i];
1172  for (int32_t i=num_kernels; i<2*num_kernels+1; i++)
1173  beta[i]=0;
1174 
1175  while (true)
1176  {
1177  //int rows=CPXgetnumrows(env, lp_cplex);
1178  //int cols=CPXgetnumcols(env, lp_cplex);
1179  //SG_PRINT("rows:%d, cols:%d (kernel:%d)\n", rows, cols, num_kernels);
1180  CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm), beta, num_kernels);
1181 
1182  set_qnorm_constraints(beta, num_kernels);
1183 
1184  status = CPXbaropt(env, lp_cplex);
1185  if ( status )
1186  SG_ERROR( "Failed to optimize Problem.\n");
1187 
1188  int solstat=0;
1189  double objval=0;
1190  status=CPXsolution(env, lp_cplex, &solstat, &objval,
1191  (double*) beta, NULL, NULL, NULL);
1192 
1193  if ( status )
1194  {
1195  CMath::display_vector(beta, num_kernels, "beta");
1196  SG_ERROR( "Failed to obtain solution.\n");
1197  }
1198 
1199  CMath::scale_vector(1/CMath::qnorm(beta, num_kernels, mkl_norm), beta, num_kernels);
1200 
1201  //SG_PRINT("[%d] %f (%f)\n", inner_iters, objval, objval_old);
1202  if ((1-abs(objval/objval_old) < 0.1*mkl_epsilon)) // && (inner_iters>2))
1203  break;
1204 
1205  objval_old=objval;
1206 
1207  inner_iters++;
1208  }
1209  SG_FREE(beta);
1210  }
1211 
1212  if ( status )
1213  SG_ERROR( "Failed to optimize Problem.\n");
1214 
1215  // obtain solution
1216  int32_t cur_numrows=(int32_t) CPXgetnumrows(env, lp_cplex);
1217  int32_t cur_numcols=(int32_t) CPXgetnumcols(env, lp_cplex);
1218  int32_t num_rows=cur_numrows;
1219  ASSERT(cur_numcols<=2*num_kernels+1);
1220 
1221  float64_t* slack=SG_MALLOC(float64_t, cur_numrows);
1222  float64_t* pi=SG_MALLOC(float64_t, cur_numrows);
1223 
1224  /* calling external lib */
1225  int solstat=0;
1226  double objval=0;
1227 
1228  if (mkl_norm==1)
1229  {
1230  status=CPXsolution(env, lp_cplex, &solstat, &objval,
1231  (double*) x, (double*) pi, (double*) slack, NULL);
1232  }
1233  else
1234  {
1235  status=CPXsolution(env, lp_cplex, &solstat, &objval,
1236  (double*) x, NULL, (double*) slack, NULL);
1237  }
1238 
1239  int32_t solution_ok = (!status) ;
1240  if ( status )
1241  SG_ERROR( "Failed to obtain solution.\n");
1242 
1243  int32_t num_active_rows=0 ;
1244  if (solution_ok)
1245  {
1246  /* 1 norm mkl */
1247  float64_t max_slack = -CMath::INFTY ;
1248  int32_t max_idx = -1 ;
1249  int32_t start_row = 1 ;
1250  if (C_mkl!=0.0)
1251  start_row+=2*(num_kernels-1);
1252 
1253  for (int32_t i = start_row; i < cur_numrows; i++) // skip first
1254  {
1255  if (mkl_norm==1)
1256  {
1257  if ((pi[i]!=0))
1258  num_active_rows++ ;
1259  else
1260  {
1261  if (slack[i]>max_slack)
1262  {
1263  max_slack=slack[i] ;
1264  max_idx=i ;
1265  }
1266  }
1267  }
1268  else // 2-norm or general q-norm
1269  {
1270  if ((CMath::abs(slack[i])<1e-6))
1271  num_active_rows++ ;
1272  else
1273  {
1274  if (slack[i]>max_slack)
1275  {
1276  max_slack=slack[i] ;
1277  max_idx=i ;
1278  }
1279  }
1280  }
1281  }
1282 
1283  // have at most max(100,num_active_rows*2) rows, if not, remove one
1284  if ( (num_rows-start_row>CMath::max(100,2*num_active_rows)) && (max_idx!=-1))
1285  {
1286  //SG_INFO( "-%i(%i,%i)",max_idx,start_row,num_rows) ;
1287  status = CPXdelrows (env, lp_cplex, max_idx, max_idx) ;
1288  if ( status )
1289  SG_ERROR( "Failed to remove an old row.\n");
1290  }
1291 
1292  //CMath::display_vector(x, num_kernels, "beta");
1293 
1294  rho = -x[2*num_kernels] ;
1295  SG_FREE(pi);
1296  SG_FREE(slack);
1297 
1298  }
1299  else
1300  {
1301  /* then something is wrong and we rather
1302  stop sooner than later */
1303  rho = 1 ;
1304  }
1305  }
1306  for (int32_t i=0; i<num_kernels; i++)
1307  new_beta[i]=x[i];
1308 
1309  SG_FREE(x);
1310 #else
1311  SG_ERROR("Cplex not enabled at compile time\n");
1312 #endif
1313  return rho;
1314 }
1315 
1317  int num_kernels, const float64_t* sumw, float64_t suma, int32_t& inner_iters)
1318 {
1319  SG_DEBUG("MKL via GLPK\n");
1320 
1321  if (mkl_norm!=1)
1322  SG_ERROR("MKL via GLPK works only for norm=1\n");
1323 
1324  float64_t obj=1.0;
1325 #ifdef USE_GLPK
1326  int32_t NUMCOLS = 2*num_kernels + 1 ;
1327  if (!lp_initialized)
1328  {
1329  SG_INFO( "creating LP\n") ;
1330 
1331  //set obj function, note glpk indexing is 1-based
1332  lpx_add_cols(lp_glpk, NUMCOLS);
1333  for (int i=1; i<=2*num_kernels; i++)
1334  {
1335  lpx_set_obj_coef(lp_glpk, i, 0);
1336  lpx_set_col_bnds(lp_glpk, i, LPX_DB, 0, 1);
1337  }
1338  for (int i=num_kernels+1; i<=2*num_kernels; i++)
1339  {
1340  lpx_set_obj_coef(lp_glpk, i, C_mkl);
1341  }
1342  lpx_set_obj_coef(lp_glpk, NUMCOLS, 1);
1343  lpx_set_col_bnds(lp_glpk, NUMCOLS, LPX_FR, 0,0); //unbound
1344 
1345  //add first row. sum[w]=1
1346  int row_index = lpx_add_rows(lp_glpk, 1);
1347  int* ind = SG_MALLOC(int, num_kernels+2);
1348  float64_t* val = SG_MALLOC(float64_t, num_kernels+2);
1349  for (int i=1; i<=num_kernels; i++)
1350  {
1351  ind[i] = i;
1352  val[i] = 1;
1353  }
1354  ind[num_kernels+1] = NUMCOLS;
1355  val[num_kernels+1] = 0;
1356  lpx_set_mat_row(lp_glpk, row_index, num_kernels, ind, val);
1357  lpx_set_row_bnds(lp_glpk, row_index, LPX_FX, 1, 1);
1358  SG_FREE(val);
1359  SG_FREE(ind);
1360 
1361  lp_initialized = true;
1362 
1363  if (C_mkl!=0.0)
1364  {
1365  for (int32_t q=1; q<num_kernels; q++)
1366  {
1367  int mat_ind[4];
1368  float64_t mat_val[4];
1369  int mat_row_index = lpx_add_rows(lp_glpk, 2);
1370  mat_ind[1] = q;
1371  mat_val[1] = 1;
1372  mat_ind[2] = q+1;
1373  mat_val[2] = -1;
1374  mat_ind[3] = num_kernels+q;
1375  mat_val[3] = -1;
1376  lpx_set_mat_row(lp_glpk, mat_row_index, 3, mat_ind, mat_val);
1377  lpx_set_row_bnds(lp_glpk, mat_row_index, LPX_UP, 0, 0);
1378  mat_val[1] = -1;
1379  mat_val[2] = 1;
1380  lpx_set_mat_row(lp_glpk, mat_row_index+1, 3, mat_ind, mat_val);
1381  lpx_set_row_bnds(lp_glpk, mat_row_index+1, LPX_UP, 0, 0);
1382  }
1383  }
1384  }
1385 
1386  int* ind=SG_MALLOC(int,num_kernels+2);
1387  float64_t* val=SG_MALLOC(float64_t, num_kernels+2);
1388  int row_index = lpx_add_rows(lp_glpk, 1);
1389  for (int32_t i=1; i<=num_kernels; i++)
1390  {
1391  ind[i] = i;
1392  val[i] = -(sumw[i-1]-suma);
1393  }
1394  ind[num_kernels+1] = 2*num_kernels+1;
1395  val[num_kernels+1] = -1;
1396  lpx_set_mat_row(lp_glpk, row_index, num_kernels+1, ind, val);
1397  lpx_set_row_bnds(lp_glpk, row_index, LPX_UP, 0, 0);
1398  SG_FREE(ind);
1399  SG_FREE(val);
1400 
1401  //optimize
1402  lpx_simplex(lp_glpk);
1403  bool res = check_lpx_status(lp_glpk);
1404  if (!res)
1405  SG_ERROR( "Failed to optimize Problem.\n");
1406 
1407  int32_t cur_numrows = lpx_get_num_rows(lp_glpk);
1408  int32_t cur_numcols = lpx_get_num_cols(lp_glpk);
1409  int32_t num_rows=cur_numrows;
1410  ASSERT(cur_numcols<=2*num_kernels+1);
1411 
1412  float64_t* col_primal = SG_MALLOC(float64_t, cur_numcols);
1413  float64_t* row_primal = SG_MALLOC(float64_t, cur_numrows);
1414  float64_t* row_dual = SG_MALLOC(float64_t, cur_numrows);
1415 
1416  for (int i=0; i<cur_numrows; i++)
1417  {
1418  row_primal[i] = lpx_get_row_prim(lp_glpk, i+1);
1419  row_dual[i] = lpx_get_row_dual(lp_glpk, i+1);
1420  }
1421  for (int i=0; i<cur_numcols; i++)
1422  col_primal[i] = lpx_get_col_prim(lp_glpk, i+1);
1423 
1424  obj = -col_primal[2*num_kernels];
1425 
1426  for (int i=0; i<num_kernels; i++)
1427  beta[i] = col_primal[i];
1428 
1429  int32_t num_active_rows=0;
1430  if(res)
1431  {
1432  float64_t max_slack = CMath::INFTY;
1433  int32_t max_idx = -1;
1434  int32_t start_row = 1;
1435  if (C_mkl!=0.0)
1436  start_row += 2*(num_kernels-1);
1437 
1438  for (int32_t i= start_row; i<cur_numrows; i++)
1439  {
1440  if (row_dual[i]!=0)
1441  num_active_rows++;
1442  else
1443  {
1444  if (row_primal[i]<max_slack)
1445  {
1446  max_slack = row_primal[i];
1447  max_idx = i;
1448  }
1449  }
1450  }
1451 
1452  if ((num_rows-start_row>CMath::max(100, 2*num_active_rows)) && max_idx!=-1)
1453  {
1454  int del_rows[2];
1455  del_rows[1] = max_idx+1;
1456  lpx_del_rows(lp_glpk, 1, del_rows);
1457  }
1458  }
1459 
1460  SG_FREE(row_dual);
1461  SG_FREE(row_primal);
1462  SG_FREE(col_primal);
1463 #else
1464  SG_ERROR("Glpk not enabled at compile time\n");
1465 #endif
1466 
1467  return obj;
1468 }
1469 
1471 {
1472  ASSERT(sumw);
1473  ASSERT(svm);
1474 
1475  int32_t nsv=svm->get_num_support_vectors();
1476  int32_t num_kernels = kernel->get_num_subkernels();
1477  float64_t* beta = SG_MALLOC(float64_t, num_kernels);
1478  int32_t nweights=0;
1479  const float64_t* old_beta = kernel->get_subkernel_weights(nweights);
1480  ASSERT(nweights==num_kernels);
1481  ASSERT(old_beta);
1482 
1483  for (int32_t i=0; i<num_kernels; i++)
1484  {
1485  beta[i]=0;
1486  sumw[i]=0;
1487  }
1488 
1489  for (int32_t n=0; n<num_kernels; n++)
1490  {
1491  beta[n]=1.0;
1492  kernel->set_subkernel_weights(SGVector<float64_t>(beta, num_kernels));
1493 
1494  for (int32_t i=0; i<nsv; i++)
1495  {
1496  int32_t ii=svm->get_support_vector(i);
1497 
1498  for (int32_t j=0; j<nsv; j++)
1499  {
1500  int32_t jj=svm->get_support_vector(j);
1501  sumw[n]+=0.5*svm->get_alpha(i)*svm->get_alpha(j)*kernel->kernel(ii,jj);
1502  }
1503  }
1504  beta[n]=0.0;
1505  }
1506 
1507  mkl_iterations++;
1508  kernel->set_subkernel_weights(SGVector<float64_t>( (float64_t*) old_beta, num_kernels));
1509 }
1510 
1511 
1512 // assumes that all constraints are satisfied
1514 {
1516  {
1517  // -- Compute ElasticnetMKL dual objective
1519  }
1520 
1521  int32_t n=get_num_support_vectors();
1522  float64_t mkl_obj=0;
1523 
1524  if (labels && kernel && kernel->get_kernel_type() == K_COMBINED)
1525  {
1526  CKernel* kn = ((CCombinedKernel*)kernel)->get_first_kernel();
1527  while (kn)
1528  {
1529  float64_t sum=0;
1530  for (int32_t i=0; i<n; i++)
1531  {
1532  int32_t ii=get_support_vector(i);
1533 
1534  for (int32_t j=0; j<n; j++)
1535  {
1536  int32_t jj=get_support_vector(j);
1537  sum+=get_alpha(i)*get_alpha(j)*kn->kernel(ii,jj);
1538  }
1539  }
1540 
1541  if (mkl_norm==1.0)
1542  mkl_obj = CMath::max(mkl_obj, sum);
1543  else
1544  mkl_obj += CMath::pow(sum, mkl_norm/(mkl_norm-1));
1545 
1546  SG_UNREF(kn);
1547  kn = ((CCombinedKernel*) kernel)->get_next_kernel();
1548  }
1549 
1550  if (mkl_norm==1.0)
1551  mkl_obj=-0.5*mkl_obj;
1552  else
1553  mkl_obj= -0.5*CMath::pow(mkl_obj, (mkl_norm-1)/mkl_norm);
1554 
1555  mkl_obj+=compute_sum_alpha();
1556  }
1557  else
1558  SG_ERROR( "cannot compute objective, labels or kernel not set\n");
1559 
1560  return -mkl_obj;
1561 }
1562 
1563 #ifdef USE_CPLEX
1564 void CMKL::set_qnorm_constraints(float64_t* beta, int32_t num_kernels)
1565 {
1566  ASSERT(num_kernels>0);
1567 
1568  float64_t* grad_beta=SG_MALLOC(float64_t, num_kernels);
1569  float64_t* hess_beta=SG_MALLOC(float64_t, num_kernels+1);
1570  float64_t* lin_term=SG_MALLOC(float64_t, num_kernels+1);
1571  int* ind=SG_MALLOC(int, num_kernels+1);
1572 
1573  //CMath::display_vector(beta, num_kernels, "beta");
1574  double const_term = 1-CMath::qsq(beta, num_kernels, mkl_norm);
1575 
1576  //SG_PRINT("const=%f\n", const_term);
1577  ASSERT(CMath::fequal(const_term, 0.0));
1578 
1579  for (int32_t i=0; i<num_kernels; i++)
1580  {
1581  grad_beta[i]=mkl_norm * pow(beta[i], mkl_norm-1);
1582  hess_beta[i]=0.5*mkl_norm*(mkl_norm-1) * pow(beta[i], mkl_norm-2);
1583  lin_term[i]=grad_beta[i] - 2*beta[i]*hess_beta[i];
1584  const_term+=grad_beta[i]*beta[i] - CMath::sq(beta[i])*hess_beta[i];
1585  ind[i]=i;
1586  }
1587  ind[num_kernels]=2*num_kernels;
1588  hess_beta[num_kernels]=0;
1589  lin_term[num_kernels]=0;
1590 
1591  int status=0;
1592  int num=CPXgetnumqconstrs (env, lp_cplex);
1593 
1594  if (num>0)
1595  {
1596  status = CPXdelqconstrs (env, lp_cplex, 0, 0);
1597  ASSERT(!status);
1598  }
1599 
1600  status = CPXaddqconstr (env, lp_cplex, num_kernels+1, num_kernels+1, const_term, 'L', ind, lin_term,
1601  ind, ind, hess_beta, NULL);
1602  ASSERT(!status);
1603 
1604  //CPXwriteprob (env, lp_cplex, "prob.lp", NULL);
1605  //CPXqpwrite (env, lp_cplex, "prob.qp");
1606 
1607  SG_FREE(grad_beta);
1608  SG_FREE(hess_beta);
1609  SG_FREE(lin_term);
1610  SG_FREE(ind);
1611 }
1612 #endif // USE_CPLEX

SHOGUN Machine Learning Toolbox - Documentation