SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QuadraticTimeMMD.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  * All rights reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions are met:
8  *
9  * 1. Redistributions of source code must retain the above copyright notice, this
10  * list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright notice,
12  * this list of conditions and the following disclaimer in the documentation
13  * and/or other materials provided with the distribution.
14  *
15  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
19  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25  *
26  * The views and conclusions contained in the software and documentation are those
27  * of the authors and should not be interpreted as representing official policies,
28  * either expressed or implied, of the Shogun Development Team.
29  */
30 
35 #include <shogun/kernel/Kernel.h>
38 
39 using namespace shogun;
40 using namespace Eigen;
41 
43 {
44  init();
45 }
46 
48  index_t m) :
49  CKernelTwoSampleTest(kernel, p_and_q, m)
50 {
51  init();
52 }
53 
55  CFeatures* q) : CKernelTwoSampleTest(kernel, p, q)
56 {
57  init();
58 }
59 
61  CKernelTwoSampleTest(custom_kernel, NULL, m)
62 {
63  init();
64 }
65 
67 {
68 }
69 
70 void CQuadraticTimeMMD::init()
71 {
72  SG_ADD(&m_num_samples_spectrum, "num_samples_spectrum", "Number of samples"
73  " for spectrum method null-distribution approximation",
75  SG_ADD(&m_num_eigenvalues_spectrum, "num_eigenvalues_spectrum", "Number of "
76  " Eigenvalues for spectrum method null-distribution approximation",
78  SG_ADD((machine_int_t*)&m_statistic_type, "statistic_type",
79  "Biased or unbiased MMD", MS_NOT_AVAILABLE);
80 
84 }
85 
87 {
88  /* split computations into three terms from JLMR paper (see documentation )*/
89  index_t m=m_m;
90  index_t n=0;
91 
92  /* check if kernel is precomputed (custom kernel) */
95  else
96  {
97  REQUIRE(m_p_and_q, "The samples are not initialized!\n");
99  }
100 
101  SG_DEBUG("Computing MMD_u with %d samples from p and %d samples from q\n",
102  m, n);
103 
104  /* init kernel with features */
106 
107  /* first term */
108  float64_t first=0;
109  for (index_t i=0; i<m; ++i)
110  {
111  /* ensure i!=j while adding up */
112  for (index_t j=0; j<m; ++j)
113  first+=i==j ? 0 : m_kernel->kernel(i,j);
114  }
115  first/=m*(m-1);
116 
117  /* second term */
118  float64_t second=0;
119  for (index_t i=m; i<m+n; ++i)
120  {
121  /* ensure i!=j while adding up */
122  for (index_t j=m; j<m+n; ++j)
123  second+=i==j ? 0 : m_kernel->kernel(i,j);
124  }
125  second/=n*(n-1);
126 
127  /* third term */
128  float64_t third=0;
129  for (index_t i=0; i<m; ++i)
130  {
131  for (index_t j=m; j<m+n; ++j)
132  third+=m_kernel->kernel(i,j);
133  }
134  third*=2.0/m/n;
135 
136  SG_DEBUG("first=%f, second=%f, third=%f\n", first, second, third);
137 
138  float64_t mmd_u=first+second-third;
139 
140  /* return m*MMD_u when m=n, (m+n)*MMD_u in general case */
141  if (m==n)
142  return m*mmd_u;
143 
144  return (m+n)*mmd_u;
145 }
146 
148 {
149  /* split computations into three terms from JLMR paper (see documentation )*/
150  index_t m=m_m;
151  index_t n=0;
152 
153  /* check if kernel is precomputed (custom kernel) */
155  n=m_kernel->get_num_vec_lhs()-m;
156  else
157  {
158  REQUIRE(m_p_and_q, "The samples are not initialized!\n");
159  n=m_p_and_q->get_num_vectors()-m;
160  }
161 
162  SG_DEBUG("Computing MMD_b with %d samples from p and %d samples from q\n",
163  m, n);
164 
165  /* init kernel with features */
167 
168  /* first term */
169  float64_t first=0;
170  for (index_t i=0; i<m; ++i)
171  {
172  for (index_t j=0; j<m; ++j)
173  first+=m_kernel->kernel(i,j);
174  }
175  first/=m*m;
176 
177  /* second term */
178  float64_t second=0;
179  for (index_t i=m; i<m+n; ++i)
180  {
181  for (index_t j=m; j<m+n; ++j)
182  second+=m_kernel->kernel(i,j);
183  }
184  second/=n*n;
185 
186  /* third term */
187  float64_t third=0;
188  for (index_t i=0; i<m; ++i)
189  {
190  for (index_t j=m; j<m+n; ++j)
191  third+=m_kernel->kernel(i,j);
192  }
193  third*=2.0/m/n;
194 
195  SG_DEBUG("first=%f, second=%f, third=%f\n", first, second, third);
196 
197  float64_t mmd_b=first+second-third;
198 
199  /* return m*MMD_b when m=n, (m+n)*MMD_b in general case */
200  if (m==n)
201  return m*mmd_b;
202 
203  return (m+n)*mmd_b;
204 }
205 
207 {
208  if (!m_kernel)
209  SG_ERROR("No kernel specified!\n")
210 
211  float64_t result=0;
212  switch (m_statistic_type)
213  {
214  case UNBIASED:
216  break;
217  case BIASED:
218  result=compute_biased_statistic();
219  break;
220  default:
221  SG_ERROR("Unknown statistic type!\n");
222  break;
223  }
224 
225  return result;
226 }
227 
229 {
230  float64_t result=0;
231 
233  {
234  case MMD2_SPECTRUM:
235  {
236 #ifdef HAVE_EIGEN3
237  /* get samples from null-distribution and compute p-value of statistic */
240  null_samples.qsort();
241  index_t pos=null_samples.find_position_to_insert(statistic);
242  result=1.0-((float64_t)pos)/null_samples.vlen;
243 #else // HAVE_EIGEN3
244  SG_ERROR("Only possible if shogun is compiled with EIGEN3 enabled\n");
245 #endif // HAVE_EIGEN3
246  break;
247  }
248 
249  case MMD2_GAMMA:
250  {
251  /* fit gamma and return cdf at statistic */
253  result=CStatistics::gamma_cdf(statistic, params[0], params[1]);
254  break;
255  }
256 
257  default:
258  result=CKernelTwoSampleTest::compute_p_value(statistic);
259  break;
260  }
261 
262  return result;
263 }
264 
266  bool multiple_kernels)
267 {
268  SGVector<float64_t> mmds;
269  if (!multiple_kernels)
270  {
271  /* just one mmd result */
272  mmds=SGVector<float64_t>(1);
273  mmds[0]=compute_statistic();
274  }
275  else
276  {
278  "multiple kernels specified, but underlying kernel is not of type "
279  "K_COMBINED\n");
280 
281  /* cast and allocate memory for results */
283  SG_REF(combined);
284  mmds=SGVector<float64_t>(combined->get_num_subkernels());
285 
286  /* iterate through all kernels and compute statistic */
287  /* TODO this might be done in parallel */
288  for (index_t i=0; i<mmds.vlen; ++i)
289  {
290  CKernel* current=combined->get_kernel(i);
291  /* temporarily replace underlying kernel and compute statistic */
292  m_kernel=current;
293  mmds[i]=compute_statistic();
294 
295  SG_UNREF(current);
296  }
297 
298  /* restore combined kernel */
299  m_kernel=combined;
300  SG_UNREF(combined);
301  }
302 
303  return mmds;
304 }
305 
307 {
308  float64_t result=0;
309 
311  {
312  case MMD2_SPECTRUM:
313  {
314 #ifdef HAVE_EIGEN3
315  /* get samples from null-distribution and compute threshold */
318  null_samples.qsort();
319  result=null_samples[index_t(CMath::floor(null_samples.vlen*(1-alpha)))];
320 #else // HAVE_EIGEN3
321  SG_ERROR("Only possible if shogun is compiled with EIGEN3 enabled\n");
322 #endif // HAVE_EIGEN3
323  break;
324  }
325 
326  case MMD2_GAMMA:
327  {
328  /* fit gamma and return inverse cdf at alpha */
330  result=CStatistics::inverse_gamma_cdf(alpha, params[0], params[1]);
331  break;
332  }
333 
334  default:
335  /* sampling null is handled here */
337  break;
338  }
339 
340  return result;
341 }
342 
343 
344 #ifdef HAVE_EIGEN3
346  index_t num_eigenvalues)
347 {
348  REQUIRE(m_kernel, "(%d, %d): No kernel set!\n", num_samples,
349  num_eigenvalues);
351  "(%d, %d): No features set and no custom kernel in use!\n",
352  num_samples, num_eigenvalues);
353 
354  index_t m=m_m;
355  index_t n=0;
356 
357  /* check if kernel is precomputed (custom kernel) */
359  n=m_kernel->get_num_vec_lhs()-m;
360  else
361  {
362  REQUIRE(m_p_and_q, "The samples are not initialized!\n");
363  n=m_p_and_q->get_num_vectors()-m;
364  }
365 
366  if (num_samples<=2)
367  {
368  SG_ERROR("Number of samples has to be at least 2, "
369  "better in the hundreds");
370  }
371 
372  if (num_eigenvalues>m+n-1)
373  SG_ERROR("Number of Eigenvalues too large\n");
374 
375  if (num_eigenvalues<1)
376  SG_ERROR("Number of Eigenvalues too small\n");
377 
378  /* imaginary matrix K=[K KL; KL' L] (MATLAB notation)
379  * K is matrix for XX, L is matrix for YY, KL is XY, LK is YX
380  * works since X and Y are concatenated here */
383 
384  /* center matrix K=H*K*H */
385  K.center();
386 
387  /* compute eigenvalues and select num_eigenvalues largest ones */
388  Map<MatrixXd> c_kernel_matrix(K.matrix, K.num_rows, K.num_cols);
389  SelfAdjointEigenSolver<MatrixXd> eigen_solver(c_kernel_matrix);
390  REQUIRE(eigen_solver.info()==Eigen::Success,
391  "Eigendecomposition failed!\n");
392  index_t max_num_eigenvalues=eigen_solver.eigenvalues().rows();
393 
394  /* precomputing terms with rho_x and rho_y of equation 10 in [1]
395  * (see documentation) */
396  float64_t rho_x=float64_t(m)/(m+n);
397  float64_t rho_y=1-rho_x;
398 
399  /* instead of using two Gaussian rv's ~ N(0,1), we'll use just one rv
400  * ~ N(0, 1/rho_x+1/rho_y) (derived from eq 10 in [1]) */
401  float64_t std_dev=CMath::sqrt(1/rho_x+1/rho_y);
402  float64_t inv_rho_x_y=1/(rho_x*rho_y);
403 
404  SG_DEBUG("Using Gaussian samples ~ N(0,%f)\n", std_dev*std_dev);
405 
406  /* finally, sample from null distribution */
407  SGVector<float64_t> null_samples(num_samples);
408  for (index_t i=0; i<num_samples; ++i)
409  {
410  null_samples[i]=0;
411  for (index_t j=0; j<num_eigenvalues; ++j)
412  {
413  /* compute the right hand multiple of eq. 10 in [1] using one RV.
414  * randn_double() gives a sample from N(0,1), we need samples
415  * from N(0,1/rho_x+1/rho_y) */
416  float64_t z_j=std_dev*CMath::randn_double();
417 
418  SG_DEBUG("z_j=%f\n", z_j);
419 
420  float64_t multiple=CMath::pow(z_j, 2);
421 
422  /* take largest EV, scale by 1/(m+n) on the fly and take abs value*/
423  float64_t eigenvalue_estimate=CMath::abs(1.0/(m+n)
424  *eigen_solver.eigenvalues()[max_num_eigenvalues-1-j]);
425 
427  multiple-=inv_rho_x_y;
428 
429  SG_DEBUG("multiple=%f, eigenvalue=%f\n", multiple,
430  eigenvalue_estimate);
431 
432  null_samples[i]+=eigenvalue_estimate*multiple;
433  }
434  }
435 
436  /* when m=n, return m*MMD^2 instead */
437  if (m==n)
438  null_samples.scale(0.5);
439 
440  return null_samples;
441 }
442 #endif // HAVE_EIGEN3
443 
445 {
446  REQUIRE(m_kernel, "No kernel set!\n");
448  "No features set and no custom kernel in use!\n");
449 
450  index_t num_data;
452  num_data=m_kernel->get_num_vec_rhs();
453  else
454  num_data=m_p_and_q->get_num_vectors();
455 
456  if (m_m!=num_data/2)
457  SG_ERROR("Currently, only equal sample sizes are supported\n");
458 
459  /* evtl. warn user not to use wrong statistic type */
461  {
462  SG_WARNING("Note: provided statistic has "
463  "to be BIASED. Please ensure that! To get rid of warning,"
464  "call %s::set_statistic_type(BIASED)\n", get_name());
465  }
466 
467  /* imaginary matrix K=[K KL; KL' L] (MATLAB notation)
468  * K is matrix for XX, L is matrix for YY, KL is XY, LK is YX
469  * works since X and Y are concatenated here */
471 
472  /* compute mean under H0 of MMD, which is
473  * meanMMD = 2/m * ( 1 - 1/m*sum(diag(KL)) );
474  * in MATLAB.
475  * Remove diagonals on the fly */
476  float64_t mean_mmd=0;
477  for (index_t i=0; i<m_m; ++i)
478  {
479  /* virtual KL matrix is in upper right corner of SHOGUN K matrix
480  * so this sums the diagonal of the matrix between X and Y*/
481  mean_mmd+=m_kernel->kernel(i, m_m+i);
482  }
483  mean_mmd=2.0/m_m*(1.0-1.0/m_m*mean_mmd);
484 
485  /* compute variance under H0 of MMD, which is
486  * varMMD = 2/m/(m-1) * 1/m/(m-1) * sum(sum( (K + L - KL - KL').^2 ));
487  * in MATLAB, so sum up all elements */
488  float64_t var_mmd=0;
489  for (index_t i=0; i<m_m; ++i)
490  {
491  for (index_t j=0; j<m_m; ++j)
492  {
493  /* dont add diagonal of all pairs of imaginary kernel matrices */
494  if (i==j || m_m+i==j || m_m+j==i)
495  continue;
496 
497  float64_t to_add=m_kernel->kernel(i, j);
498  to_add+=m_kernel->kernel(m_m+i, m_m+j);
499  to_add-=m_kernel->kernel(i, m_m+j);
500  to_add-=m_kernel->kernel(m_m+i, j);
501  var_mmd+=CMath::pow(to_add, 2);
502  }
503  }
504  var_mmd*=2.0/m_m/(m_m-1)*1.0/m_m/(m_m-1);
505 
506  /* parameters for gamma distribution */
507  float64_t a=CMath::pow(mean_mmd, 2)/var_mmd;
508  float64_t b=var_mmd*m_m / mean_mmd;
509 
510  SGVector<float64_t> result(2);
511  result[0]=a;
512  result[1]=b;
513 
514  return result;
515 }
516 
518  num_samples_spectrum)
519 {
520  m_num_samples_spectrum=num_samples_spectrum;
521 }
522 
524  index_t num_eigenvalues_spectrum)
525 {
526  m_num_eigenvalues_spectrum=num_eigenvalues_spectrum;
527 }
528 
530  statistic_type)
531 {
532  m_statistic_type=statistic_type;
533 }
534 

SHOGUN Machine Learning Toolbox - Documentation