SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GaussianProcessMachine.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) 2013 Roman Votyakov
8  */
9 
10 #include <shogun/lib/config.h>
11 
12 #ifdef HAVE_EIGEN3
13 
16 #include <shogun/kernel/Kernel.h>
18 
20 
21 using namespace shogun;
22 using namespace Eigen;
23 
25 {
26  init();
27 }
28 
30 {
31  init();
32  set_inference_method(method);
33 }
34 
35 void CGaussianProcessMachine::init()
36 {
37  m_method=NULL;
38 
39  SG_ADD((CSGObject**) &m_method, "inference_method", "Inference method",
40  MS_AVAILABLE);
41 }
42 
44 {
45  SG_UNREF(m_method);
46 }
47 
49 {
50  REQUIRE(m_method, "Inference method should not be NULL\n")
51 
52  CFeatures* feat;
53 
54  // use latent features for FITC inference method
55  if (m_method->get_inference_type()==INF_FITC)
56  {
57  CFITCInferenceMethod* fitc_method=
59  feat=fitc_method->get_latent_features();
60  SG_UNREF(fitc_method);
61  }
62  else
63  feat=m_method->get_features();
64 
65  // get kernel and compute kernel matrix: K(feat, data)*scale^2
66  CKernel* training_kernel=m_method->get_kernel();
67  CKernel* kernel=CKernel::obtain_from_generic(training_kernel->clone());
68  SG_UNREF(training_kernel);
69 
70  kernel->init(feat, data);
71 
72  // get kernel matrix and create eigen representation of it
73  SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix();
74  Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols);
75 
76  // compute Ks=Ks*scale^2
77  eigen_Ks*=CMath::sq(m_method->get_scale());
78 
79  // cleanup
80  SG_UNREF(feat);
81  SG_UNREF(kernel);
82 
83  // get alpha and create eigen representation of it
84  SGVector<float64_t> alpha=m_method->get_alpha();
85  Map<VectorXd> eigen_alpha(alpha.vector, alpha.vlen);
86 
87  // get mean and create eigen representation of it
88  CMeanFunction* mean_function=m_method->get_mean();
89  SGVector<float64_t> m=mean_function->get_mean_vector(data);
90  Map<VectorXd> eigen_m(m.vector, m.vlen);
91  SG_UNREF(mean_function);
92 
93  // compute mean: mu=Ks'*alpha+m
95  Map<VectorXd> eigen_mu(mu.vector, mu.vlen);
96  eigen_mu=eigen_Ks.adjoint()*eigen_alpha+eigen_m;
97 
98  return mu;
99 }
100 
102  CFeatures* data)
103 {
104  REQUIRE(m_method, "Inference method should not be NULL\n")
105 
106  CFeatures* feat;
107 
108  // use latent features for FITC inference method
109  if (m_method->get_inference_type()==INF_FITC)
110  {
111  CFITCInferenceMethod* fitc_method=
113  feat=fitc_method->get_latent_features();
114  SG_UNREF(fitc_method);
115  }
116  else
117  feat=m_method->get_features();
118 
119  SG_REF(data);
120 
121  // get kernel and compute kernel matrix: K(data, data)*scale^2
122  CKernel* training_kernel=m_method->get_kernel();
123  CKernel* kernel=CKernel::obtain_from_generic(training_kernel->clone());
124  SG_UNREF(training_kernel);
125 
126  kernel->init(data, data);
127 
128  // get kernel matrix and create eigen representation of it
129  SGMatrix<float64_t> k_tsts=kernel->get_kernel_matrix();
130  Map<MatrixXd> eigen_Kss(k_tsts.matrix, k_tsts.num_rows, k_tsts.num_cols);
131 
132  // compute Kss=Kss*scale^2
133  eigen_Kss*=CMath::sq(m_method->get_scale());
134 
135  // compute kernel matrix: K(feat, data)*scale^2
136  kernel->init(feat, data);
137 
138  // get kernel matrix and create eigen representation of it
139  SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix();
140  Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols);
141 
142  // compute Ks=Ks*scale^2
143  eigen_Ks*=CMath::sq(m_method->get_scale());
144 
145  // cleanup
146  SG_UNREF(kernel);
147  SG_UNREF(feat);
148  SG_UNREF(data);
149 
150  // get shogun representation of cholesky and create eigen representation
151  SGMatrix<float64_t> L=m_method->get_cholesky();
152  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
153 
154  // result variance vector
155  SGVector<float64_t> s2(k_tsts.num_cols);
156  Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
157 
158  if (eigen_L.isUpperTriangular())
159  {
160  // get shogun of diagonal sigma vector and create eigen representation
161  SGVector<float64_t> sW=m_method->get_diagonal_vector();
162  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
163 
164  // solve L' * V = sW * Ks and compute V.^2
165  MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
166  eigen_sW.asDiagonal()*eigen_Ks);
167  MatrixXd eigen_sV=eigen_V.cwiseProduct(eigen_V);
168 
169  eigen_s2=eigen_Kss.diagonal()-eigen_sV.colwise().sum().adjoint();
170  }
171  else
172  {
173  // M = Ks .* (L * Ks)
174  MatrixXd eigen_M=eigen_Ks.cwiseProduct(eigen_L*eigen_Ks);
175  eigen_s2=eigen_Kss.diagonal()+eigen_M.colwise().sum().adjoint();
176  }
177 
178  return s2;
179 }
180 
181 #endif /* HAVE_EIGEN3 */

SHOGUN Machine Learning Toolbox - Documentation