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

SHOGUN Machine Learning Toolbox - Documentation