SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HSIC.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) The Shogun Machine Learning Toolbox
3  * Written (w) 2012-2013 Heiko Strathmann
4  * Written (w) 2014 Soumyajit De
5  * All rights reserved.
6  *
7  * Redistribution and use in source and binary forms, with or without
8  * modification, are permitted provided that the following conditions are met:
9  *
10  * 1. Redistributions of source code must retain the above copyright notice, this
11  * list of conditions and the following disclaimer.
12  * 2. Redistributions in binary form must reproduce the above copyright notice,
13  * this list of conditions and the following disclaimer in the documentation
14  * and/or other materials provided with the distribution.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
20  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  *
27  * The views and conclusions contained in the software and documentation are those
28  * of the authors and should not be interpreted as representing official policies,
29  * either expressed or implied, of the Shogun Development Team.
30  */
31 
32 #include <shogun/statistics/HSIC.h>
35 #include <shogun/kernel/Kernel.h>
37 
38 using namespace shogun;
39 
41 {
42  init();
43 }
44 
45 CHSIC::CHSIC(CKernel* kernel_p, CKernel* kernel_q, CFeatures* p,
46  CFeatures* q) :
47  CKernelIndependenceTest(kernel_p, kernel_q, p, q)
48 {
49  init();
50 
51  if (p && q && p->get_num_vectors()!=q->get_num_vectors())
52  {
53  SG_ERROR("Only features with equal number of vectors are currently "
54  "possible\n");
55  }
56  else
57  m_num_features=p->get_num_vectors();
58 }
59 
61 {
62 }
63 
64 void CHSIC::init()
65 {
66  SG_ADD(&m_num_features, "num_features",
67  "Number of features from each of the distributions",
69 
70  m_num_features=0;
71 }
72 
74 {
75  SG_DEBUG("entering!\n");
76 
77  REQUIRE(m_kernel_p && m_kernel_q, "No or only one kernel specified!\n");
78 
79  REQUIRE(m_p && m_q, "features needed!\n")
80 
81  /* compute kernel matrices */
84 
85  /* center matrices (MATLAB: Kc=H*K*H) */
86  K.center();
87 
88  /* compute MATLAB: sum(sum(Kc' .* (L))), which is biased HSIC */
89  index_t m=m_num_features;
90  SG_DEBUG("Number of samples %d!\n", m);
91 
92  float64_t result=0;
93  for (index_t i=0; i<m; ++i)
94  {
95  for (index_t j=0; j<m; ++j)
96  result+=K(j, i)*L(i, j);
97  }
98 
99  /* return m times statistic */
100  result/=m;
101 
102  SG_DEBUG("leaving!\n");
103 
104  return result;
105 }
106 
108 {
109  float64_t result=0;
111  {
112  case HSIC_GAMMA:
113  {
114  /* fit gamma and return cdf at statistic */
116  result=CStatistics::gamma_cdf(statistic, params[0], params[1]);
117  break;
118  }
119 
120  default:
121  /* sampling null is handled there */
122  result=CIndependenceTest::compute_p_value(statistic);
123  break;
124  }
125 
126  return result;
127 }
128 
130 {
131  float64_t result=0;
133  {
134  case HSIC_GAMMA:
135  {
136  /* fit gamma and return inverse cdf at statistic */
138  result=CStatistics::inverse_gamma_cdf(alpha, params[0], params[1]);
139  break;
140  }
141 
142  default:
143  /* sampling null is handled there */
145  break;
146  }
147 
148  return result;
149 }
150 
152 {
153  REQUIRE(m_kernel_p && m_kernel_q, "No or only one kernel specified!\n");
154 
155  REQUIRE(m_p && m_q, "features needed!\n")
156 
157  index_t m=m_num_features;
158 
159  /* compute kernel matrices */
162 
163  /* compute sum and trace of uncentered kernel matrices, needed later */
164  float64_t trace_K=0;
165  float64_t trace_L=0;
166  float64_t sum_K=0;
167  float64_t sum_L=0;
168  for (index_t i=0; i<m; ++i)
169  {
170  trace_K+=K(i,i);
171  trace_L+=L(i,i);
172  for (index_t j=0; j<m; ++j)
173  {
174  sum_K+=K(i,j);
175  sum_L+=L(i,j);
176  }
177  }
178  SG_DEBUG("sum_K: %f, sum_L: %f, trace_K: %f, trace_L: %f\n", sum_K, sum_L,
179  trace_K, trace_L);
180 
181  /* center both matrices: K=H*K*H, L=H*L*H in MATLAB */
182  K.center();
183  L.center();
184 
185  /* compute the trace of MATLAB: (1/6 * Kc.*Lc).^2 Ü */
186  float64_t trace=0;
187  for (index_t i=0; i<m; ++i)
188  trace+=CMath::pow(K(i,i)*L(i,i), 2);
189 
190  trace/=36.0;
191  SG_DEBUG("trace %f\n", trace)
192 
193  /* compute sum of elements of MATLAB: (1/6 * Kc.*Lc).^2 */
194  float64_t sum=0;
195  for (index_t i=0; i<m; ++i)
196  {
197  for (index_t j=0; j<m; ++j)
198  sum+=CMath::pow(K(i,j)*L(i,j), 2);
199  }
200  sum/=36.0;
201  SG_DEBUG("sum %f\n", sum)
202 
203  /* compute MATLAB: 1/m/(m-1)*(sum(sum(varHSIC)) - sum(diag(varHSIC))),
204  * second term is bias correction */
205  float64_t var_hsic=1.0/m/(m-1)*(sum-trace);
206  SG_DEBUG("1.0/m/(m-1)*(sum-trace): %f\n", var_hsic)
207 
208  /* finally, compute variance of hsic under H0
209  * MATLAB: varHSIC = 72*(m-4)*(m-5)/m/(m-1)/(m-2)/(m-3) * varHSIC */
210  var_hsic=72.0*(m-4)*(m-5)/m/(m-1)/(m-2)/(m-3)*var_hsic;
211  SG_DEBUG("var_hsic: %f\n", var_hsic)
212 
213  /* compute mean of matrices with diagonal elements zero on the base of sums
214  * and trace from K and L which were computed above */
215  float64_t mu_x=1.0/m/(m-1)*(sum_K-trace_K);
216  float64_t mu_y=1.0/m/(m-1)*(sum_L-trace_L);
217  SG_DEBUG("mu_x: %f, mu_y: %f\n", mu_x, mu_y)
218 
219  /* compute mean under H0, MATLAB: 1/m * ( 1 +muX*muY - muX - muY ) */
220  float64_t m_hsic=1.0/m*(1+mu_x*mu_y-mu_x-mu_y);
221  SG_DEBUG("m_hsic: %f\n", m_hsic)
222 
223  /* finally, compute parameters of gamma distirbution */
224  float64_t a=CMath::pow(m_hsic, 2)/var_hsic;
225  float64_t b=var_hsic*m/m_hsic;
226  SG_DEBUG("a: %f, b: %f\n", a, b)
227 
228  SGVector<float64_t> result(2);
229  result[0]=a;
230  result[1]=b;
231 
232  SG_DEBUG("leaving!\n")
233  return result;
234 }
235 
237 {
238  SG_DEBUG("entering!\n");
239 
241 
242  /* distinguish between custom and normal kernels */
244  {
245  /* custom kernels need to to be initialised when a subset is added */
246  CCustomKernel* custom_kernel_p=(CCustomKernel*)m_kernel_p;
247  K=custom_kernel_p->get_kernel_matrix();
248  }
249  else
250  {
251  /* need to init the kernel if kernel is not precomputed - if subsets of
252  * features are in the stack (for permutation), this will handle it */
253  m_kernel_p->init(m_p, m_p);
255  }
256 
257  SG_DEBUG("leaving!\n");
258 
259  return K;
260 }
261 
263 {
264  SG_DEBUG("entering!\n");
265 
267 
268  /* now second half of data for L */
270  {
271  /* custom kernels need to to be initialised - no subsets here */
272  CCustomKernel* custom_kernel_q=(CCustomKernel*)m_kernel_q;
273  L=custom_kernel_q->get_kernel_matrix();
274  }
275  else
276  {
277  /* need to init the kernel if kernel is not precomputed */
278  m_kernel_q->init(m_q, m_q);
280  }
281 
282  SG_DEBUG("leaving!\n");
283 
284  return L;
285 }
286 
288 {
289  SG_DEBUG("entering!\n")
290 
291  /* replace current kernel via precomputed custom kernel and call superclass
292  * method */
293 
294  /* backup references to old kernels */
295  CKernel* kernel_p=m_kernel_p;
296  CKernel* kernel_q=m_kernel_q;
297 
298  /* init kernels before to be sure that everything is fine
299  * kernel function between two samples from different distributions
300  * is never computed - in fact, they may as well have different features */
301  m_kernel_p->init(m_p, m_p);
302  m_kernel_q->init(m_q, m_q);
303 
304  /* precompute kernel matrices */
305  CCustomKernel* precomputed_p=new CCustomKernel(m_kernel_p);
306  CCustomKernel* precomputed_q=new CCustomKernel(m_kernel_q);
307  SG_REF(precomputed_p);
308  SG_REF(precomputed_q);
309 
310  /* temporarily replace own kernels */
311  m_kernel_p=precomputed_p;
312  m_kernel_q=precomputed_q;
313 
314  /* use superclass sample_null which shuffles the entries for one
315  * distribution using index permutation on rows and columns of
316  * kernel matrix from one distribution, while accessing the other
317  * in its original order and then compute statistic */
319 
320  /* restore kernels */
321  m_kernel_p=kernel_p;
322  m_kernel_q=kernel_q;
323 
324  SG_UNREF(precomputed_p);
325  SG_UNREF(precomputed_q);
326 
327  SG_DEBUG("leaving!\n")
328  return null_samples;
329 }
330 
332 {
334  m_num_features=p->get_num_vectors();
335 }
336 
338 {
340  m_num_features=q->get_num_vectors();
341 }
342 

SHOGUN Machine Learning Toolbox - Documentation