SHOGUN  4.1.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  * 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  feat=fitc_method->get_inducing_features();
91  }
92  else
93  feat=m_method->get_features();
94 
95  // get kernel and compute kernel matrix: K(feat, data)*scale^2
96  CKernel* training_kernel=m_method->get_kernel();
97  CKernel* kernel=CKernel::obtain_from_generic(training_kernel->clone());
98  SG_UNREF(training_kernel);
99 
100  kernel->init(feat, data);
101 
102  // get kernel matrix and create eigen representation of it
103  SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix();
104  Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols);
105 
106  // compute Ks=Ks*scale^2
107  eigen_Ks*=CMath::sq(m_method->get_scale());
108 
109  // cleanup
110  SG_UNREF(feat);
111  SG_UNREF(kernel);
112 
113  // get alpha and create eigen representation of it
114  SGVector<float64_t> alpha=m_method->get_alpha();
115  Map<VectorXd> eigen_alpha(alpha.vector, alpha.vlen);
116 
117  // get mean and create eigen representation of it
118  CMeanFunction* mean_function=m_method->get_mean();
119  SGVector<float64_t> mean=mean_function->get_mean_vector(data);
120  Map<VectorXd> eigen_mean(mean.vector, mean.vlen);
121  SG_UNREF(mean_function);
122 
123  const index_t C=alpha.vlen/k_trts.num_rows;
124  const index_t n=k_trts.num_rows;
125  const index_t m=k_trts.num_cols;
126 
127  // compute mean: mu=Ks'*alpha+m
128  SGVector<float64_t> mu(C*m);
129  Map<MatrixXd> eigen_mu_matrix(mu.vector,C,m);
130 
131  for(index_t bl=0; bl<C; bl++)
132  eigen_mu_matrix.block(bl,0,1,m)=(eigen_Ks.adjoint()*eigen_alpha.block(bl*n,0,n,1)+eigen_mean).transpose();
133 
134  return mu;
135 }
136 
138  CFeatures* data)
139 {
140  REQUIRE(m_method, "Inference method should not be NULL\n")
141 
142  CFeatures* feat;
143 
144  // use inducing features for FITC inference method
145  if (m_method->get_inference_type()==INF_FITC_REGRESSION ||
146  m_method->get_inference_type()==INF_FITC_LAPLACIAN_SINGLE)
147  {
148  CSingleFITCLaplacianBase* fitc_method=
149  dynamic_cast<CSingleFITCLaplacianBase *>(m_method);
150  REQUIRE(fitc_method, "Inference method %s must support FITC inference\n",
151  m_method->get_name());
152  feat=fitc_method->get_inducing_features();
153  }
154  else
155  feat=m_method->get_features();
156 
157  SG_REF(data);
158 
159  // get kernel and compute kernel matrix: K(data, data)*scale^2
160  CKernel* training_kernel=m_method->get_kernel();
161  CKernel* kernel=CKernel::obtain_from_generic(training_kernel->clone());
162  SG_UNREF(training_kernel);
163 
164  kernel->init(data, data);
165 
166  // get kernel matrix and create eigen representation of it
167  SGVector<float64_t> k_tsts=kernel->get_kernel_diagonal();
168  Map<VectorXd> eigen_Kss_diag(k_tsts.vector, k_tsts.vlen);
169 
170  // compute Kss=Kss*scale^2
171  eigen_Kss_diag*=CMath::sq(m_method->get_scale());
172 
173  // compute kernel matrix: K(feat, data)*scale^2
174  kernel->init(feat, data);
175 
176  // get kernel matrix and create eigen representation of it
177  SGMatrix<float64_t> k_trts=kernel->get_kernel_matrix();
178  Map<MatrixXd> eigen_Ks(k_trts.matrix, k_trts.num_rows, k_trts.num_cols);
179 
180  // compute Ks=Ks*scale^2
181  eigen_Ks*=CMath::sq(m_method->get_scale());
182 
183  // cleanup
184  SG_UNREF(kernel);
185  SG_UNREF(feat);
186  SG_UNREF(data);
187 
188  // get shogun representation of cholesky and create eigen representation
189  SGMatrix<float64_t> L=m_method->get_cholesky();
190  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
191 
192  SGVector<float64_t> alpha=m_method->get_alpha();
193  const index_t n=k_trts.num_rows;
194  const index_t m=k_tsts.vlen;
195  const index_t C=alpha.vlen/n;
196  // result variance vector
197  SGVector<float64_t> s2(m*C*C);
198  Map<VectorXd> eigen_s2(s2.vector, s2.vlen);
199 
200  if (eigen_L.isUpperTriangular())
201  {
202  if (alpha.vlen==L.num_rows)
203  {
204  //binary case
205  // get shogun of diagonal sigma vector and create eigen representation
206  SGVector<float64_t> sW=m_method->get_diagonal_vector();
207  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
208 
209  // solve L' * V = sW * Ks and compute V.^2
210  MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
211  eigen_sW.asDiagonal()*eigen_Ks);
212  MatrixXd eigen_sV=eigen_V.cwiseProduct(eigen_V);
213 
214  eigen_s2=eigen_Kss_diag-eigen_sV.colwise().sum().adjoint();
215  }
216  else
217  {
218  if (m_method->supports_multiclass())
219  {
220  //multiclass case
221  //see the reference code of the gist link, which is based on the algorithm 3.4 of the GPML textbook
222  Map<MatrixXd> &eigen_M=eigen_L;
223  eigen_s2.fill(0);
224 
225  SGMatrix<float64_t> E=m_method->get_multiclass_E();
226  Map<MatrixXd> eigen_E(E.matrix, E.num_rows, E.num_cols);
227  ASSERT(E.num_cols==alpha.vlen);
228  for(index_t bl_i=0; bl_i<C; bl_i++)
229  {
230  //n by m
231  MatrixXd bi=eigen_E.block(0,bl_i*n,n,n)*eigen_Ks;
232  MatrixXd c_cav=eigen_M.triangularView<Upper>().adjoint().solve(bi);
233  c_cav=eigen_M.triangularView<Upper>().solve(c_cav);
234 
235  for(index_t bl_j=0; bl_j<C; bl_j++)
236  {
237  MatrixXd bj=eigen_E.block(0,bl_j*n,n,n)*eigen_Ks;
238  for (index_t idx_m=0; idx_m<m; idx_m++)
239  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();
240  }
241  for (index_t idx_m=0; idx_m<m; idx_m++)
242  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();
243  }
244  }
245  else
246  {
247  SG_ERROR("Unsupported inference method!\n");
248  return s2;
249  }
250  }
251  }
252  else
253  {
254  // M = Ks .* (L * Ks)
255  MatrixXd eigen_M=eigen_Ks.cwiseProduct(eigen_L*eigen_Ks);
256  eigen_s2=eigen_Kss_diag+eigen_M.colwise().sum().adjoint();
257  }
258 
259  return s2;
260 }
261 
262 #endif /* HAVE_EIGEN3 */

SHOGUN Machine Learning Toolbox - Documentation