SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
GaussianProcessMachine.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) The Shogun Machine Learning Toolbox
3  * Written (w) 2014 Wu Lin
4  * Written (W) 2013 Roman Votyakov
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  * Code adapted from
32  * Gaussian Process Machine Learning Toolbox
33  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
34  * and
35  * https://gist.github.com/yorkerlin/8a36e8f9b298aa0246a4
36  */
37 
38 #include <shogun/lib/config.h>
39 
40 #ifdef HAVE_EIGEN3
41 
44 #include <shogun/kernel/Kernel.h>
46 
48 
49 using namespace shogun;
50 using namespace Eigen;
51 
53 {
54  init();
55 }
56 
58 {
59  init();
60  set_inference_method(method);
61 }
62 
63 void CGaussianProcessMachine::init()
64 {
65  m_method=NULL;
66 
67  SG_ADD((CSGObject**) &m_method, "inference_method", "Inference method",
68  MS_AVAILABLE);
69 }
70 
72 {
73  SG_UNREF(m_method);
74 }
75 
77 {
78  REQUIRE(m_method, "Inference method should not be NULL\n")
79 
80  CFeatures* feat;
81 
82  // use inducing features for FITC inference method
83  if (m_method->get_inference_type()==INF_FITC_REGRESSION ||
84  m_method->get_inference_type()==INF_FITC_LAPLACIAN_SINGLE)
85  {
86  CSingleFITCLaplacianBase* fitc_method=
87  dynamic_cast<CSingleFITCLaplacianBase *>(m_method);
88  REQUIRE(fitc_method, "Inference method %s does not support FITC inference\n",
89  m_method->get_name());
90  fitc_method->optimize_inducing_features();
91  feat=fitc_method->get_inducing_features();
92  }
93  else
94  feat=m_method->get_features();
95 
96  // get kernel and compute kernel matrix: K(feat, data)*scale^2
97  CKernel* training_kernel=m_method->get_kernel();
98  CKernel* kernel=CKernel::obtain_from_generic(training_kernel->clone());
99  SG_UNREF(training_kernel);
100 
101  kernel->init(feat, data);
102 
103  // get kernel matrix and create eigen representation of it
104  SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix();
105  Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols);
106 
107  // compute Ks=Ks*scale^2
108  eigen_Ks*=CMath::sq(m_method->get_scale());
109 
110  // cleanup
111  SG_UNREF(feat);
112  SG_UNREF(kernel);
113 
114  // get alpha and create eigen representation of it
115  SGVector<float64_t> alpha=m_method->get_alpha();
116  Map<VectorXd> eigen_alpha(alpha.vector, alpha.vlen);
117 
118  // get mean and create eigen representation of it
119  CMeanFunction* mean_function=m_method->get_mean();
120  SGVector<float64_t> mean=mean_function->get_mean_vector(data);
121  Map<VectorXd> eigen_mean(mean.vector, mean.vlen);
122  SG_UNREF(mean_function);
123 
124  const index_t C=alpha.vlen/k_trts.num_rows;
125  const index_t n=k_trts.num_rows;
126  const index_t m=k_trts.num_cols;
127 
128  // compute mean: mu=Ks'*alpha+m
129  SGVector<float64_t> mu(C*m);
130  Map<MatrixXd> eigen_mu_matrix(mu.vector,C,m);
131 
132  for(index_t bl=0; bl<C; bl++)
133  eigen_mu_matrix.block(bl,0,1,m)=(eigen_Ks.adjoint()*eigen_alpha.block(bl*n,0,n,1)+eigen_mean).transpose();
134 
135  return mu;
136 }
137 
139  CFeatures* data)
140 {
141  REQUIRE(m_method, "Inference method should not be NULL\n")
142 
143  CFeatures* feat;
144 
145  bool is_FITC=false;
146  // use inducing features for FITC inference method
147  if (m_method->get_inference_type()==INF_FITC_REGRESSION ||
148  m_method->get_inference_type()==INF_FITC_LAPLACIAN_SINGLE)
149  {
150  CSingleFITCLaplacianBase* fitc_method=
151  dynamic_cast<CSingleFITCLaplacianBase *>(m_method);
152  REQUIRE(fitc_method, "Inference method %s must support FITC inference\n",
153  m_method->get_name());
154  fitc_method->optimize_inducing_features();
155  feat=fitc_method->get_inducing_features();
156  is_FITC=true;
157  }
158  else
159  feat=m_method->get_features();
160 
161  SG_REF(data);
162 
163  // get kernel and compute kernel matrix: K(data, data)*scale^2
164  CKernel* training_kernel=m_method->get_kernel();
165  CKernel* kernel=CKernel::obtain_from_generic(training_kernel->clone());
166  SG_UNREF(training_kernel);
167  kernel->init(data, data);
168 
169  // get kernel matrix and create eigen representation of it
170  SGVector<float64_t> k_tsts=kernel->get_kernel_diagonal();
171  Map<VectorXd> eigen_Kss_diag(k_tsts.vector, k_tsts.vlen);
172 
173  // compute Kss=Kss*scale^2
174  eigen_Kss_diag*=CMath::sq(m_method->get_scale());
175 
176  // compute kernel matrix: K(feat, data)*scale^2
177  kernel->init(feat, data);
178 
179  // get kernel matrix and create eigen representation of it
180  SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix();
181  Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols);
182 
183  // compute Ks=Ks*scale^2
184  eigen_Ks*=CMath::sq(m_method->get_scale());
185 
186  // cleanup
187  SG_UNREF(kernel);
188  SG_UNREF(feat);
189  SG_UNREF(data);
190 
191  // get shogun representation of cholesky and create eigen representation
192  SGMatrix<float64_t> L=m_method->get_cholesky();
193  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
194 
195  SGVector<float64_t> alpha=m_method->get_alpha();
196  const index_t n=k_trts.num_rows;
197  const index_t m=k_tsts.vlen;
198  const index_t C=alpha.vlen/n;
199  // result variance vector
200  SGVector<float64_t> s2(m*C*C);
201  Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
202 
203  if (eigen_L.isUpperTriangular() && !is_FITC)
204  {
205  if (alpha.vlen==L.num_rows)
206  {
207  //binary case
208  // get shogun of diagonal sigma vector and create eigen representation
209  SGVector<float64_t> sW=m_method->get_diagonal_vector();
210  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
211  // solve L' * V = sW * Ks and compute V.^2
212  MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
213  eigen_sW.asDiagonal()*eigen_Ks);
214  MatrixXd eigen_sV=eigen_V.cwiseProduct(eigen_V);
215 
216  eigen_s2=eigen_Kss_diag-eigen_sV.colwise().sum().adjoint();
217  }
218  else
219  {
220  if (m_method->supports_multiclass())
221  {
222  //multiclass case
223  //see the reference code of the gist link, which is based on the algorithm 3.4 of the GPML textbook
224  Map<MatrixXd> &eigen_M=eigen_L;
225  eigen_s2.fill(0);
226 
227  SGMatrix<float64_t> E=m_method->get_multiclass_E();
228  Map<MatrixXd> eigen_E(E.matrix, E.num_rows, E.num_cols);
229  ASSERT(E.num_cols==alpha.vlen);
230  for(index_t bl_i=0; bl_i<C; bl_i++)
231  {
232  //n by m
233  MatrixXd bi=eigen_E.block(0,bl_i*n,n,n)*eigen_Ks;
234  MatrixXd c_cav=eigen_M.triangularView<Upper>().adjoint().solve(bi);
235  c_cav=eigen_M.triangularView<Upper>().solve(c_cav);
236 
237  for(index_t bl_j=0; bl_j<C; bl_j++)
238  {
239  MatrixXd bj=eigen_E.block(0,bl_j*n,n,n)*eigen_Ks;
240  for (index_t idx_m=0; idx_m<m; idx_m++)
241  eigen_s2[bl_j+(bl_i+idx_m*C)*C]=(bj.block(0,idx_m,n,1).array()*c_cav.block(0,idx_m,n,1).array()).sum();
242  }
243  for (index_t idx_m=0; idx_m<m; idx_m++)
244  eigen_s2[bl_i+(bl_i+idx_m*C)*C]+=eigen_Kss_diag(idx_m)-(eigen_Ks.block(0,idx_m,n,1).array()*bi.block(0,idx_m,n,1).array()).sum();
245  }
246  }
247  else
248  {
249  SG_ERROR("Unsupported inference method!\n");
250  return s2;
251  }
252  }
253  }
254  else
255  {
256  // M = Ks .* (L * Ks)
257  MatrixXd eigen_M=eigen_Ks.cwiseProduct(eigen_L*eigen_Ks);
258  eigen_s2=eigen_Kss_diag+eigen_M.colwise().sum().adjoint();
259  }
260 
261  return s2;
262 }
263 
264 #endif /* HAVE_EIGEN3 */
virtual bool init(CFeatures *lhs, CFeatures *rhs)
Definition: Kernel.cpp:98
The Inference Method base class.
int32_t index_t
Definition: common.h:62
virtual CSGObject * clone()
Definition: SGObject.cpp:1361
SGVector< float64_t > get_posterior_variances(CFeatures *data)
static T sq(T x)
Definition: Math.h:450
Definition: SGMatrix.h:20
#define SG_ERROR(...)
Definition: SGIO.h:129
#define REQUIRE(x,...)
Definition: SGIO.h:206
SGVector< float64_t > get_posterior_means(CFeatures *data)
index_t num_cols
Definition: SGMatrix.h:378
virtual SGVector< float64_t > get_mean_vector(const CFeatures *features) const =0
An abstract class of the mean function.
Definition: MeanFunction.h:28
The Fully Independent Conditional Training inference base class for Laplace and regression for 1-D la...
SGMatrix< float64_t > get_kernel_matrix()
Definition: Kernel.h:219
#define SG_REF(x)
Definition: SGObject.h:51
index_t num_rows
Definition: SGMatrix.h:376
index_t vlen
Definition: SGVector.h:494
#define ASSERT(x)
Definition: SGIO.h:201
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:112
virtual const char * get_name() const
static CKernel * obtain_from_generic(CSGObject *kernel)
Definition: Kernel.cpp:885
#define SG_UNREF(x)
Definition: SGObject.h:52
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
The class Features is the base class of all feature objects.
Definition: Features.h:68
SGVector< float64_t > get_kernel_diagonal(SGVector< float64_t > preallocated=SGVector< float64_t >())
Definition: Kernel.h:230
The Kernel base class.
Definition: Kernel.h:158
virtual CFeatures * get_inducing_features()
#define SG_ADD(...)
Definition: SGObject.h:81

SHOGUN Machine Learning Toolbox - Documentation