SHOGUN  v3.0.0
 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* kernel=m_method->get_kernel();
67  kernel->init(feat, data);
68 
69  // get kernel matrix and create eigen representation of it
70  SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix();
71  Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols);
72 
73  // compute Ks=Ks*scale^2
74  eigen_Ks*=CMath::sq(m_method->get_scale());
75 
76  // cleanup
77  SG_UNREF(feat);
78  SG_UNREF(kernel);
79 
80  // get alpha and create eigen representation of it
81  SGVector<float64_t> alpha=m_method->get_alpha();
82  Map<VectorXd> eigen_alpha(alpha.vector, alpha.vlen);
83 
84  // get mean and create eigen representation of it
85  CMeanFunction* mean_function=m_method->get_mean();
86  SGVector<float64_t> m=mean_function->get_mean_vector(data);
87  Map<VectorXd> eigen_m(m.vector, m.vlen);
88  SG_UNREF(mean_function);
89 
90  // compute mean: mu=Ks'*alpha+m
92  Map<VectorXd> eigen_mu(mu.vector, mu.vlen);
93  eigen_mu=eigen_Ks.adjoint()*eigen_alpha+eigen_m;
94 
95  return mu;
96 }
97 
99  CFeatures* data)
100 {
101  REQUIRE(m_method, "Inference method should not be NULL\n")
102 
103  CFeatures* feat;
104 
105  // use latent features for FITC inference method
106  if (m_method->get_inference_type()==INF_FITC)
107  {
108  CFITCInferenceMethod* fitc_method=
110  feat=fitc_method->get_latent_features();
111  SG_UNREF(fitc_method);
112  }
113  else
114  feat=m_method->get_features();
115 
116  SG_REF(data);
117 
118  // get kernel and compute kernel matrix: K(data, data)*scale^2
119  CKernel* kernel=m_method->get_kernel();
120  kernel->init(data, data);
121 
122  // get kernel matrix and create eigen representation of it
123  SGMatrix<float64_t> k_tsts=kernel->get_kernel_matrix();
124  Map<MatrixXd> eigen_Kss(k_tsts.matrix, k_tsts.num_rows, k_tsts.num_cols);
125 
126  // compute Kss=Kss*scale^2
127  eigen_Kss*=CMath::sq(m_method->get_scale());
128 
129  kernel->cleanup();
130 
131  // compute kernel matrix: K(feat, data)*scale^2
132  kernel->init(feat, data);
133 
134  // get kernel matrix and create eigen representation of it
135  SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix();
136  Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols);
137 
138  // compute Ks=Ks*scale^2
139  eigen_Ks*=CMath::sq(m_method->get_scale());
140 
141  // cleanup
142  SG_UNREF(kernel);
143  SG_UNREF(feat);
144  SG_UNREF(data);
145 
146  // get shogun representation of cholesky and create eigen representation
147  SGMatrix<float64_t> L=m_method->get_cholesky();
148  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
149 
150  // result variance vector
151  SGVector<float64_t> s2(k_tsts.num_cols);
152  Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
153 
154  if (eigen_L.isUpperTriangular())
155  {
156  // get shogun of diagonal sigma vector and create eigen representation
157  SGVector<float64_t> sW=m_method->get_diagonal_vector();
158  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
159 
160  // solve L' * V = sW * Ks and compute V.^2
161  MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
162  eigen_sW.asDiagonal()*eigen_Ks);
163  MatrixXd eigen_sV=eigen_V.cwiseProduct(eigen_V);
164 
165  eigen_s2=eigen_Kss.diagonal()-eigen_sV.colwise().sum().adjoint();
166  }
167  else
168  {
169  // M = Ks .* (L * Ks)
170  MatrixXd eigen_M=eigen_Ks.cwiseProduct(eigen_L*eigen_Ks);
171  eigen_s2=eigen_Kss.diagonal()+eigen_M.colwise().sum().adjoint();
172  }
173 
174  return s2;
175 }
176 
177 #endif /* HAVE_EIGEN3 */

SHOGUN Machine Learning Toolbox - Documentation