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

SHOGUN Machine Learning Toolbox - Documentation