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

SHOGUN Machine Learning Toolbox - Documentation