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

SHOGUN Machine Learning Toolbox - Documentation