SHOGUN  6.1.3
VarianceH1.h
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 - 2017 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 #ifndef VARIANCE_H1__H_
33 #define VARIANCE_H1__H_
34 
35 #include <vector>
36 #include <shogun/lib/common.h>
38 #include <shogun/lib/SGVector.h>
39 #include <shogun/kernel/Kernel.h>
43 
44 using std::vector;
45 
46 namespace shogun
47 {
48 
49 namespace internal
50 {
51 
52 namespace mmd
53 {
54 #ifndef DOXYGEN_SHOULD_SKIP_THIS
55 struct VarianceH1
56 {
57  VarianceH1() : m_lambda(1E-5), m_free_terms(true)
58  {
59  }
60 
61  void init_terms()
62  {
63  m_sum_x=0;
64  m_sum_y=0;
65  m_sum_xy=0;
66  m_sum_sq_x=0;
67  m_sum_sq_y=0;
68  m_sum_sq_xy=0;
69 
70  m_sum_colwise_x.resize(m_n_x);
71  m_sum_colwise_y.resize(m_n_y);
72  m_sum_rowwise_xy.resize(m_n_x);
73  m_sum_colwise_xy.resize(m_n_y);
74  std::fill(m_sum_colwise_x.begin(), m_sum_colwise_x.end(), 0);
75  std::fill(m_sum_colwise_y.begin(), m_sum_colwise_y.end(), 0);
76  std::fill(m_sum_rowwise_xy.begin(), m_sum_rowwise_xy.end(), 0);
77  std::fill(m_sum_colwise_xy.begin(), m_sum_colwise_xy.end(), 0);
78 
79  if (m_second_order_terms.rows()==m_n_x && m_second_order_terms.cols()==m_n_x)
80  m_second_order_terms.setZero();
81  else
82  m_second_order_terms=Eigen::MatrixXd::Zero(m_n_x, m_n_x);
83  }
84 
85  void free_terms()
86  {
87  if (m_free_terms)
88  {
89  m_sum_colwise_x.resize(0);
90  m_sum_colwise_y.resize(0);
91  m_sum_rowwise_xy.resize(0);
92  m_sum_colwise_xy.resize(0);
93  m_second_order_terms=Eigen::MatrixXd::Zero(0, 0);
94  }
95  }
96 
97  template <typename T>
98  void add_terms(T kernel_value, index_t i, index_t j)
99  {
100  if (i<m_n_x && j<m_n_x)
101  {
102  m_sum_x+=2*kernel_value;
103  m_sum_sq_x+=2*kernel_value*kernel_value;
104  m_sum_colwise_x[i]+=kernel_value;
105  m_sum_colwise_x[j]+=kernel_value;
106  m_second_order_terms(i, j)+=kernel_value;
107  m_second_order_terms(j, i)+=kernel_value;
108  }
109  else if (i>=m_n_x && j>=m_n_x)
110  {
111  m_sum_y+=2*kernel_value;
112  m_sum_sq_y+=2*kernel_value*kernel_value;
113  m_sum_colwise_y[i-m_n_x]+=kernel_value;
114  m_sum_colwise_y[j-m_n_x]+=kernel_value;
115  m_second_order_terms(i-m_n_x, j-m_n_x)+=kernel_value;
116  m_second_order_terms(j-m_n_x, i-m_n_x)+=kernel_value;
117  }
118  else if (i<m_n_x && j>=m_n_x)
119  {
120  m_sum_xy+=kernel_value;
121  m_sum_sq_xy+=kernel_value*kernel_value;
122  if (j-i!=m_n_x)
123  {
124  m_second_order_terms(i, j-m_n_x)-=kernel_value;
125  m_second_order_terms(j-m_n_x, i)-=kernel_value;
126  }
127  m_sum_rowwise_xy[i]+=kernel_value;
128  m_sum_colwise_xy[j-m_n_x]+=kernel_value;
129  }
130  }
131 
132  float64_t compute_variance_estimate()
133  {
134  Eigen::Map<Eigen::VectorXd> map_sum_colwise_x(m_sum_colwise_x.data(), m_sum_colwise_x.size());
135  Eigen::Map<Eigen::VectorXd> map_sum_colwise_y(m_sum_colwise_y.data(), m_sum_colwise_y.size());
136  Eigen::Map<Eigen::VectorXd> map_sum_rowwise_xy(m_sum_rowwise_xy.data(), m_sum_rowwise_xy.size());
137  Eigen::Map<Eigen::VectorXd> map_sum_colwise_xy(m_sum_colwise_xy.data(), m_sum_colwise_xy.size());
138 
139  auto t_0=(map_sum_colwise_x.dot(map_sum_colwise_x)-m_sum_sq_x)/m_n_x/(m_n_x-1)/(m_n_x-2);
140  auto t_1=CMath::sq(m_sum_x/m_n_x/(m_n_x-1));
141 
142  auto t_2=map_sum_colwise_x.dot(map_sum_rowwise_xy)*2/m_n_x/(m_n_x-1)/m_n_y;
143  auto t_3=m_sum_x*m_sum_xy*2/m_n_x/m_n_x/(m_n_x-1)/m_n_y;
144 
145  auto t_4=(map_sum_colwise_y.dot(map_sum_colwise_y)-m_sum_sq_y)/m_n_y/(m_n_y-1)/(m_n_y-2);
146  auto t_5=CMath::sq(m_sum_y/m_n_y/(m_n_y-1));
147 
148  auto t_6=map_sum_colwise_y.dot(map_sum_colwise_xy)*2/m_n_y/(m_n_y-1)/m_n_x;
149  auto t_7=m_sum_y*m_sum_xy*2/m_n_y/m_n_y/(m_n_y-1)/m_n_x;
150 
151  auto t_8=(map_sum_rowwise_xy.dot(map_sum_rowwise_xy)-m_sum_sq_xy)/m_n_y/(m_n_y-1)/m_n_x;
152  auto t_9=2*CMath::sq(m_sum_xy/m_n_x/m_n_y);
153  auto t_10=(map_sum_colwise_xy.dot(map_sum_colwise_xy)-m_sum_sq_xy)/m_n_x/(m_n_x-1)/m_n_y;
154 
155  auto var_first=(t_0-t_1)-t_2+t_3+(t_4-t_5)-t_6+t_7+(t_8-t_9+t_10);
156  var_first*=4.0*(m_n_x-2)/m_n_x/(m_n_x-1);
157 
158  auto var_second=2.0/m_n_x/m_n_y/(m_n_x-1)/(m_n_y-1)*m_second_order_terms.array().square().sum();
159 
160  auto variance_estimate=var_first+var_second;
161  if (variance_estimate<0)
162  variance_estimate=var_second;
163 
164  return variance_estimate;
165  }
166 
167  template <class Kernel>
168  float64_t operator()(const Kernel& kernel)
169  {
170  ASSERT(m_n_x>0 && m_n_y>0);
171  ASSERT(m_n_x==m_n_y);
172  const index_t size=m_n_x+m_n_y;
173  init_terms();
174  for (auto j=0; j<size; ++j)
175  {
176  for (auto i=0; i<j; ++i)
177  add_terms(kernel(i, j), i, j);
178  }
179  auto variance_estimate=compute_variance_estimate();
180  free_terms();
181  return variance_estimate;
182  }
183 
184  SGVector<float64_t> operator()(const KernelManager& kernel_mgr)
185  {
186  ASSERT(m_n_x>0 && m_n_y>0);
187  ASSERT(m_n_x==m_n_y);
188  ASSERT(kernel_mgr.num_kernels()>0);
189 
190  const index_t size=m_n_x+m_n_y;
191  SGVector<float64_t> result(kernel_mgr.num_kernels());
192  SelfAdjointPrecomputedKernel kernel_functor(SGVector<float32_t>(size*(size+1)/2));
193  for (auto k=0; k<kernel_mgr.num_kernels(); ++k)
194  {
195  auto kernel=kernel_mgr.kernel_at(k);
196  ASSERT(kernel);
197  kernel_functor.precompute(kernel);
198  init_terms();
199  for (auto i=0; i<size; ++i)
200  {
201  for (auto j=i+1; j<size; ++j)
202  add_terms(kernel_functor(i, j), i, j);
203  }
204  result[k]=compute_variance_estimate();
205  }
206 
207  free_terms();
208  return result;
209  }
210 
211  SGVector<float64_t> test_power(const KernelManager& kernel_mgr)
212  {
213  ASSERT(m_n_x>0 && m_n_y>0);
214  ASSERT(m_n_x==m_n_y);
215  ASSERT(kernel_mgr.num_kernels()>0);
216  ComputeMMD compute_mmd_job;
217  compute_mmd_job.m_n_x=m_n_x;
218  compute_mmd_job.m_n_y=m_n_y;
219  compute_mmd_job.m_stype=ST_UNBIASED_FULL;
220 
221  const index_t size=m_n_x+m_n_y;
222  SGVector<float64_t> result(kernel_mgr.num_kernels());
223  SelfAdjointPrecomputedKernel kernel_functor(SGVector<float32_t>(size*(size+1)/2));
224  for (auto k=0; k<kernel_mgr.num_kernels(); ++k)
225  {
226  auto kernel=kernel_mgr.kernel_at(k);
227  ASSERT(kernel);
228  kernel_functor.precompute(kernel);
229  init_terms();
230  for (auto i=0; i<size; ++i)
231  {
232  for (auto j=i+1; j<size; ++j)
233  add_terms(kernel_functor(i, j), i, j);
234  }
235  auto var_est=compute_variance_estimate();
236  auto mmd_est=compute_mmd_job(kernel_functor);
237  result[k]=mmd_est/CMath::sqrt(var_est+m_lambda);
238  }
239 
240  free_terms();
241  return result;
242  }
243 
244  index_t m_n_x;
245  index_t m_n_y;
246 
247  float64_t m_sum_x;
248  float64_t m_sum_y;
249  float64_t m_sum_xy;
250  float64_t m_sum_sq_x;
251  float64_t m_sum_sq_y;
252  float64_t m_sum_sq_xy;
253  float64_t m_lambda;
254 
255  vector<float64_t> m_sum_colwise_x;
256  vector<float64_t> m_sum_colwise_y;
257  vector<float64_t> m_sum_rowwise_xy;
258  vector<float64_t> m_sum_colwise_xy;
259  Eigen::MatrixXd m_second_order_terms;
260 
261  bool m_free_terms;
262 };
263 #endif // DOXYGEN_SHOULD_SKIP_THIS
264 }
265 
266 }
267 
268 }
269 
270 #endif // VARIANCE_H1__H_
int32_t index_t
Definition: common.h:72
static T sqrt(T x)
Definition: Math.h:428
static T sq(T x)
Definition: Math.h:418
#define ASSERT(x)
Definition: SGIO.h:176
double float64_t
Definition: common.h:60
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18

SHOGUN Machine Learning Toolbox - Documentation