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

SHOGUN Machine Learning Toolbox - Documentation