SHOGUN  4.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
MMDKernelSelectionComb.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) 2012-2013 Heiko Strathmann
8  */
9 
10 #include <shogun/lib/config.h>
11 #ifdef USE_GPL_SHOGUN
12 
16 
17 using namespace shogun;
18 
19 CMMDKernelSelectionComb::CMMDKernelSelectionComb() :
21 {
22  init();
23 }
24 
25 CMMDKernelSelectionComb::CMMDKernelSelectionComb(
27 {
28  init();
29 }
30 
31 CMMDKernelSelectionComb::~CMMDKernelSelectionComb()
32 {
33 }
34 
35 void CMMDKernelSelectionComb::init()
36 {
37  SG_ADD(&m_opt_max_iterations, "opt_max_iterations", "Maximum number of "
38  "iterations for qp solver", MS_NOT_AVAILABLE);
39  SG_ADD(&m_opt_epsilon, "opt_epsilon", "Stopping criterion for qp solver",
41  SG_ADD(&m_opt_low_cut, "opt_low_cut", "Low cut value for optimization "
42  "kernel weights", MS_NOT_AVAILABLE);
43 
44  /* sensible values for optimization */
45  m_opt_max_iterations=10000;
46  m_opt_epsilon=10E-15;
47  m_opt_low_cut=10E-7;
48 }
49 
50 CKernel* CMMDKernelSelectionComb::select_kernel()
51 {
52  /* cast is safe due to assertion in constructor */
53  CCombinedKernel* combined=(CCombinedKernel*)m_estimator->get_kernel();
54 
55  /* optimise for kernel weights and set them */
56  SGVector<float64_t> weights=compute_measures();
57  combined->set_subkernel_weights(weights);
58 
59  /* note that kernel is SG_REF'ed from getter above */
60  return combined;
61 }
62 
63 /* no reference counting, use the static context constructor of SGMatrix */
64 SGMatrix<float64_t> CMMDKernelSelectionComb::m_Q=SGMatrix<float64_t>(false);
65 
66 const float64_t* CMMDKernelSelectionComb::get_Q_col(uint32_t i)
67 {
68  return &m_Q[m_Q.num_rows*i];
69 }
70 
72 void CMMDKernelSelectionComb::print_state(libqp_state_T state)
73 {
74  SG_SDEBUG("libqp state: primal=%f\n", state.QP);
75 }
76 
77 SGVector<float64_t> CMMDKernelSelectionComb::solve_optimization(
79 {
80  /* readability */
81  index_t num_kernels=mmds.vlen;
82 
83  /* compute sum of mmds to generate feasible point for convex program */
84  float64_t sum_mmds=0;
85  for (index_t i=0; i<mmds.vlen; ++i)
86  sum_mmds+=mmds[i];
87 
88  /* QP: 0.5*x'*Q*x + f'*x
89  * subject to
90  * mmds'*x = b
91  * LB[i] <= x[i] <= UB[i] for all i=1..n */
92  SGVector<float64_t> Q_diag(num_kernels);
93  SGVector<float64_t> f(num_kernels);
94  SGVector<float64_t> lb(num_kernels);
95  SGVector<float64_t> ub(num_kernels);
96  SGVector<float64_t> weights(num_kernels);
97 
98  /* init everything, there are two cases possible: i) at least one mmd is
99  * is positive, ii) all mmds are negative */
100  bool one_pos=false;
101  for (index_t i=0; i<mmds.vlen; ++i)
102  {
103  if (mmds[i]>0)
104  {
105  SG_DEBUG("found at least one positive MMD\n")
106  one_pos=true;
107  break;
108  }
109  }
110 
111  if (!one_pos)
112  {
113  SG_WARNING("All mmd estimates are negative. This is techically possible,"
114  "although extremely rare. Consider using different kernels. "
115  "This combination will lead to a bad two-sample test. Since any"
116  "combination is bad, will now just return equally distributed "
117  "kernel weights\n");
118 
119  /* if no element is positive, we can choose arbritary weights since
120  * the results will be bad anyway */
121  weights.set_const(1.0/num_kernels);
122  }
123  else
124  {
125  SG_DEBUG("one MMD entry is positive, performing optimisation\n")
126  /* do optimisation, init vectors */
127  for (index_t i=0; i<num_kernels; ++i)
128  {
129  Q_diag[i]=m_Q(i,i);
130  f[i]=0;
131  lb[i]=0;
132  ub[i]=CMath::INFTY;
133 
134  /* initial point has to be feasible, i.e. mmds'*x = b */
135  weights[i]=1.0/sum_mmds;
136  }
137 
138  /* start libqp solver with desired parameters */
139  SG_DEBUG("starting libqp optimization\n")
140  libqp_state_T qp_exitflag=libqp_gsmo_solver(&get_Q_col, Q_diag.vector,
141  f.vector, mmds.vector,
142  one_pos ? 1 : -1,
143  lb.vector, ub.vector,
144  weights.vector, num_kernels, m_opt_max_iterations,
145  m_opt_epsilon, &(CMMDKernelSelectionComb::print_state));
146 
147  SG_DEBUG("libqp returns: nIts=%d, exit_flag: %d\n", qp_exitflag.nIter,
148  qp_exitflag.exitflag);
149 
150  /* set really small entries to zero and sum up for normalization */
151  float64_t sum_weights=0;
152  for (index_t i=0; i<weights.vlen; ++i)
153  {
154  if (weights[i]<m_opt_low_cut)
155  {
156  SG_DEBUG("lowcut: weight[%i]=%f<%f setting to zero\n", i, weights[i],
157  m_opt_low_cut);
158  weights[i]=0;
159  }
160 
161  sum_weights+=weights[i];
162  }
163 
164  /* normalize (allowed since problem is scale invariant) */
165  for (index_t i=0; i<weights.vlen; ++i)
166  weights[i]/=sum_weights;
167  }
168 
169  return weights;
170 }
171 #endif //USE_GPL_SHOGUN
virtual void set_subkernel_weights(SGVector< float64_t > weights)
int32_t index_t
Definition: common.h:62
static const float64_t INFTY
infinity
Definition: Math.h:2048
Base class for kernel selection for MMD-based two-sample test statistic implementations. Provides abstract methods for selecting kernels and computing criteria or kernel weights for the implemented method. In order to implement new methods for kernel selection, simply write a new implementation of this class.
Kernel two sample test base class. Provides an interface for performing a two-sample test using a ker...
index_t vlen
Definition: SGVector.h:492
CKernel * get_kernel(int32_t idx)
double float64_t
Definition: common.h:50
The Combined kernel is used to combine a number of kernels into a single CombinedKernel object by lin...
#define SG_DEBUG(...)
Definition: SGIO.h:107
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
#define SG_SDEBUG(...)
Definition: SGIO.h:168
The Kernel base class.
Definition: Kernel.h:159
#define SG_WARNING(...)
Definition: SGIO.h:128
#define SG_ADD(...)
Definition: SGObject.h:81
void set_const(T const_elem)
Definition: SGVector.cpp:150

SHOGUN Machine Learning Toolbox - Documentation