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  float64_t result=0;
91  for (index_t i=0; i<m; ++i)
92  {
93  for (index_t j=0; j<m; ++j)
94  result+=K(j, i)*L(i, j);
95  }
96 
97  /* return m times statistic */
98  result/=m;
99 
100  SG_DEBUG("leaving!\n");
101 
102  return result;
103 }
104 
106 {
107  float64_t result=0;
109  {
110  case HSIC_GAMMA:
111  {
112  /* fit gamma and return cdf at statistic */
114  result=CStatistics::gamma_cdf(statistic, params[0], params[1]);
115  break;
116  }
117 
118  default:
119  /* sampling null is handled there */
120  result=CIndependenceTest::compute_p_value(statistic);
121  break;
122  }
123 
124  return result;
125 }
126 
128 {
129  float64_t result=0;
131  {
132  case HSIC_GAMMA:
133  {
134  /* fit gamma and return inverse cdf at statistic */
136  result=CStatistics::inverse_gamma_cdf(alpha, params[0], params[1]);
137  break;
138  }
139 
140  default:
141  /* sampling null is handled there */
143  break;
144  }
145 
146  return result;
147 }
148 
150 {
151  REQUIRE(m_kernel_p && m_kernel_q, "No or only one kernel specified!\n");
152 
153  REQUIRE(m_p && m_q, "features needed!\n")
154 
155  index_t m=m_num_features;
156 
157  /* compute kernel matrices */
160 
161  /* compute sum and trace of uncentered kernel matrices, needed later */
162  float64_t trace_K=0;
163  float64_t trace_L=0;
164  float64_t sum_K=0;
165  float64_t sum_L=0;
166  for (index_t i=0; i<m; ++i)
167  {
168  trace_K+=K(i,i);
169  trace_L+=L(i,i);
170  for (index_t j=0; j<m; ++j)
171  {
172  sum_K+=K(i,j);
173  sum_L+=L(i,j);
174  }
175  }
176  SG_DEBUG("sum_K: %f, sum_L: %f, trace_K: %f, trace_L: %f\n", sum_K, sum_L,
177  trace_K, trace_L);
178 
179  /* center both matrices: K=H*K*H, L=H*L*H in MATLAB */
180  K.center();
181  L.center();
182 
183  /* compute the trace of MATLAB: (1/6 * Kc.*Lc).^2 Ü */
184  float64_t trace=0;
185  for (index_t i=0; i<m; ++i)
186  trace+=CMath::pow(K(i,i)*L(i,i), 2);
187 
188  trace/=36.0;
189  SG_DEBUG("trace %f\n", trace)
190 
191  /* compute sum of elements of MATLAB: (1/6 * Kc.*Lc).^2 */
192  float64_t sum=0;
193  for (index_t i=0; i<m; ++i)
194  {
195  for (index_t j=0; j<m; ++j)
196  sum+=CMath::pow(K(i,j)*L(i,j), 2);
197  }
198  sum/=36.0;
199  SG_DEBUG("sum %f\n", sum)
200 
201  /* compute MATLAB: 1/m/(m-1)*(sum(sum(varHSIC)) - sum(diag(varHSIC))),
202  * second term is bias correction */
203  float64_t var_hsic=1.0/m/(m-1)*(sum-trace);
204  SG_DEBUG("1.0/m/(m-1)*(sum-trace): %f\n", var_hsic)
205 
206  /* finally, compute variance of hsic under H0
207  * MATLAB: varHSIC = 72*(m-4)*(m-5)/m/(m-1)/(m-2)/(m-3) * varHSIC */
208  var_hsic=72.0*(m-4)*(m-5)/m/(m-1)/(m-2)/(m-3)*var_hsic;
209  SG_DEBUG("var_hsic: %f\n", var_hsic)
210 
211  /* compute mean of matrices with diagonal elements zero on the base of sums
212  * and trace from K and L which were computed above */
213  float64_t mu_x=1.0/m/(m-1)*(sum_K-trace_K);
214  float64_t mu_y=1.0/m/(m-1)*(sum_L-trace_L);
215  SG_DEBUG("mu_x: %f, mu_y: %f\n", mu_x, mu_y)
216 
217  /* compute mean under H0, MATLAB: 1/m * ( 1 +muX*muY - muX - muY ) */
218  float64_t m_hsic=1.0/m*(1+mu_x*mu_y-mu_x-mu_y);
219  SG_DEBUG("m_hsic: %f\n", m_hsic)
220 
221  /* finally, compute parameters of gamma distirbution */
222  float64_t a=CMath::pow(m_hsic, 2)/var_hsic;
223  float64_t b=var_hsic*m/m_hsic;
224  SG_DEBUG("a: %f, b: %f\n", a, b)
225 
226  SGVector<float64_t> result(2);
227  result[0]=a;
228  result[1]=b;
229 
230  SG_DEBUG("leaving!\n")
231  return result;
232 }
233 
235 {
236  SG_DEBUG("entering!\n");
237 
239 
240  /* distinguish between custom and normal kernels */
242  {
243  /* custom kernels need to to be initialised when a subset is added */
244  CCustomKernel* custom_kernel_p=(CCustomKernel*)m_kernel_p;
245  K=custom_kernel_p->get_kernel_matrix();
246  }
247  else
248  {
249  /* need to init the kernel if kernel is not precomputed - if subsets of
250  * features are in the stack (for permutation), this will handle it */
251  m_kernel_p->init(m_p, m_p);
253  }
254 
255  SG_DEBUG("leaving!\n");
256 
257  return K;
258 }
259 
261 {
262  SG_DEBUG("entering!\n");
263 
265 
266  /* now second half of data for L */
268  {
269  /* custom kernels need to to be initialised - no subsets here */
270  CCustomKernel* custom_kernel_q=(CCustomKernel*)m_kernel_q;
271  L=custom_kernel_q->get_kernel_matrix();
272  }
273  else
274  {
275  /* need to init the kernel if kernel is not precomputed */
276  m_kernel_q->init(m_q, m_q);
278  }
279 
280  SG_DEBUG("leaving!\n");
281 
282  return L;
283 }
284 
286 {
287  SG_DEBUG("entering!\n")
288 
289  /* replace current kernel via precomputed custom kernel and call superclass
290  * method */
291 
292  /* backup references to old kernels */
293  CKernel* kernel_p=m_kernel_p;
294  CKernel* kernel_q=m_kernel_q;
295 
296  /* init kernels before to be sure that everything is fine
297  * kernel function between two samples from different distributions
298  * is never computed - in fact, they may as well have different features */
299  m_kernel_p->init(m_p, m_p);
300  m_kernel_q->init(m_q, m_q);
301 
302  /* precompute kernel matrices */
303  CCustomKernel* precomputed_p=new CCustomKernel(m_kernel_p);
304  CCustomKernel* precomputed_q=new CCustomKernel(m_kernel_q);
305  SG_REF(precomputed_p);
306  SG_REF(precomputed_q);
307 
308  /* temporarily replace own kernels */
309  m_kernel_p=precomputed_p;
310  m_kernel_q=precomputed_q;
311 
312  /* use superclass sample_null which shuffles the entries for one
313  * distribution using index permutation on rows and columns of
314  * kernel matrix from one distribution, while accessing the other
315  * in its original order and then compute statistic */
317 
318  /* restore kernels */
319  m_kernel_p=kernel_p;
320  m_kernel_q=kernel_q;
321 
322  SG_UNREF(precomputed_p);
323  SG_UNREF(precomputed_q);
324 
325  SG_DEBUG("leaving!\n")
326  return null_samples;
327 }

SHOGUN Machine Learning Toolbox - Documentation