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

SHOGUN Machine Learning Toolbox - Documentation