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

SHOGUN Machine Learning Toolbox - Documentation