SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
ExactInferenceMethod.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  * Copyright (C) 2012 Jacob Walker
9  * Copyright (C) 2013 Roman Votyakov
10  *
11  * Code adapted from Gaussian Process Machine Learning Toolbox
12  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
13  */
14 
16 
17 #ifdef HAVE_EIGEN3
18 
23 
24 using namespace shogun;
25 using namespace Eigen;
26 
28 {
29 }
30 
32  CMeanFunction* m, CLabels* lab, CLikelihoodModel* mod) :
33  CInferenceMethod(kern, feat, m, lab, mod)
34 {
35 }
36 
38 {
39 }
40 
42 {
44 
45  if (!m_gradient_update)
46  {
47  update_deriv();
48  update_mean();
49  update_cov();
50  m_gradient_update=true;
52  }
53 }
54 
56 {
57  SG_DEBUG("entering\n");
58 
60  update_chol();
61  update_alpha();
62  m_gradient_update=false;
64 
65  SG_DEBUG("leaving\n");
66 }
67 
69 {
71 
73  "Exact inference method can only use Gaussian likelihood function\n")
75  "Labels must be type of CRegressionLabels\n")
76 }
77 
79 {
81  update();
82 
83  // get the sigma variable from the Gaussian likelihood model
85  float64_t sigma=lik->get_sigma();
86  SG_UNREF(lik);
87 
88  // compute diagonal vector: sW=1/sigma
90  result.fill_vector(result.vector, m_features->get_num_vectors(), 1.0/sigma);
91 
92  return result;
93 }
94 
96 {
98  update();
99 
100  // get the sigma variable from the Gaussian likelihood model
102  float64_t sigma=lik->get_sigma();
103  SG_UNREF(lik);
104 
105  // create eigen representation of alpha and L
106  Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
108 
109  // get labels and mean vectors and create eigen representation
110  SGVector<float64_t> y=((CRegressionLabels*) m_labels)->get_labels();
111  Map<VectorXd> eigen_y(y.vector, y.vlen);
113  Map<VectorXd> eigen_m(m.vector, m.vlen);
114 
115  // compute negative log of the marginal likelihood:
116  // nlZ=(y-m)'*alpha/2+sum(log(diag(L)))+n*log(2*pi*sigma^2)/2
117  float64_t result=(eigen_y-eigen_m).dot(eigen_alpha)/2.0+
118  eigen_L.diagonal().array().log().sum()+m_L.num_rows*
119  CMath::log(2*CMath::PI*CMath::sq(sigma))/2.0;
120 
121  return result;
122 }
123 
125 {
127  update();
128 
130 }
131 
133 {
135  update();
136 
137  return SGMatrix<float64_t>(m_L);
138 }
139 
141 {
143 
144  return SGVector<float64_t>(m_mu);
145 }
146 
148 {
150 
151  return SGMatrix<float64_t>(m_Sigma);
152 }
153 
155 {
156  // get the sigma variable from the Gaussian likelihood model
158  float64_t sigma=lik->get_sigma();
159  SG_UNREF(lik);
160 
161  /* check whether to allocate cholesky memory */
164 
165  /* creates views on kernel and cholesky matrix and perform cholesky */
168  LLT<MatrixXd> llt(K*(CMath::sq(m_scale)/CMath::sq(sigma))+
169  MatrixXd::Identity(m_ktrtr.num_rows, m_ktrtr.num_cols));
170  L=llt.matrixU();
171 }
172 
174 {
175  // get the sigma variable from the Gaussian likelihood model
177  float64_t sigma=lik->get_sigma();
178  SG_UNREF(lik);
179 
180  // get labels and mean vector and create eigen representation
181  SGVector<float64_t> y=((CRegressionLabels*) m_labels)->get_labels();
182  Map<VectorXd> eigen_y(y.vector, y.vlen);
184  Map<VectorXd> eigen_m(m.vector, m.vlen);
185 
187 
188  /* creates views on cholesky matrix and alpha and solve system
189  * (L * L^T) * a = y for a */
192 
193  a=L.triangularView<Upper>().adjoint().solve(eigen_y-eigen_m);
194  a=L.triangularView<Upper>().solve(a);
195 
196  a/=CMath::sq(sigma);
197 }
198 
200 {
201  // create eigen representataion of kernel matrix and alpha
203  Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
204 
205  // get mean and create eigen representation of it
207  Map<VectorXd> eigen_m(m.vector, m.vlen);
208 
209  m_mu=SGVector<float64_t>(m.vlen);
210  Map<VectorXd> eigen_mu(m_mu.vector, m_mu.vlen);
211 
212  eigen_mu=eigen_K*CMath::sq(m_scale)*eigen_alpha;
213 }
214 
216 {
217  // create eigen representataion of upper triangular factor L^T and kernel
218  // matrix
221 
223  Map<MatrixXd> eigen_Sigma(m_Sigma.matrix, m_Sigma.num_rows,
224  m_Sigma.num_cols);
225 
226  // compute V = L^(-1) * K, using upper triangular factor L^T
227  MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
228  eigen_K*CMath::sq(m_scale));
229 
231  float64_t sigma=lik->get_sigma();
232  SG_UNREF(lik);
233  eigen_V = eigen_V/sigma;
234 
235  // compute covariance matrix of the posterior: Sigma = K - V^T * V
236  eigen_Sigma=eigen_K*CMath::sq(m_scale)-eigen_V.adjoint()*eigen_V;
237 }
238 
240 {
241  // get the sigma variable from the Gaussian likelihood model
243  float64_t sigma=lik->get_sigma();
244  SG_UNREF(lik);
245 
246  // create eigen representation of derivative matrix and cholesky
248  Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
249 
251  Map<MatrixXd> eigen_Q(m_Q.matrix, m_Q.num_rows, m_Q.num_cols);
252 
253  // solve L * L' * Q = I
254  eigen_Q=eigen_L.triangularView<Upper>().adjoint().solve(
255  MatrixXd::Identity(m_L.num_rows, m_L.num_cols));
256  eigen_Q=eigen_L.triangularView<Upper>().solve(eigen_Q);
257 
258  // divide Q by sigma^2
259  eigen_Q/=CMath::sq(sigma);
260 
261  // create eigen representation of alpha and compute Q=Q-alpha*alpha'
262  eigen_Q-=eigen_alpha*eigen_alpha.transpose();
263 }
264 
266  const TParameter* param)
267 {
268  REQUIRE(!strcmp(param->m_name, "scale"), "Can't compute derivative of "
269  "the nagative log marginal likelihood wrt %s.%s parameter\n",
270  get_name(), param->m_name)
271 
273  Map<MatrixXd> eigen_Q(m_Q.matrix, m_Q.num_rows, m_Q.num_cols);
274 
275  SGVector<float64_t> result(1);
276 
277  // compute derivative wrt kernel scale: dnlZ=sum(Q.*K*scale*2)/2
278  result[0]=(eigen_Q.cwiseProduct(eigen_K)*m_scale*2.0).sum()/2.0;
279 
280  return result;
281 }
282 
284  CInferenceMethod* inference)
285 {
286  if (inference==NULL)
287  return NULL;
288 
289  if (inference->get_inference_type()!=INF_EXACT)
290  SG_SERROR("Provided inference is not of type CExactInferenceMethod!\n")
291 
292  SG_REF(inference);
293  return (CExactInferenceMethod*)inference;
294 }
295 
297  const TParameter* param)
298 {
299  REQUIRE(!strcmp(param->m_name, "sigma"), "Can't compute derivative of "
300  "the nagative log marginal likelihood wrt %s.%s parameter\n",
301  m_model->get_name(), param->m_name)
302 
303  // get the sigma variable from the Gaussian likelihood model
305  float64_t sigma=lik->get_sigma();
306  SG_UNREF(lik);
307 
308  // create eigen representation of the matrix Q
309  Map<MatrixXd> eigen_Q(m_Q.matrix, m_Q.num_rows, m_Q.num_cols);
310 
311  SGVector<float64_t> result(1);
312 
313  // compute derivative wrt likelihood model parameter sigma:
314  // dnlZ=sigma^2*trace(Q)
315  result[0]=CMath::sq(sigma)*eigen_Q.trace();
316 
317  return result;
318 }
319 
321  const TParameter* param)
322 {
323  // create eigen representation of the matrix Q
324  Map<MatrixXd> eigen_Q(m_Q.matrix, m_Q.num_rows, m_Q.num_cols);
325 
326  REQUIRE(param, "Param not set\n");
327  SGVector<float64_t> result;
328  int64_t len=const_cast<TParameter *>(param)->m_datatype.get_num_elements();
329  result=SGVector<float64_t>(len);
330 
331  for (index_t i=0; i<result.vlen; i++)
332  {
334 
335  if (result.vlen==1)
336  dK=m_kernel->get_parameter_gradient(param);
337  else
338  dK=m_kernel->get_parameter_gradient(param, i);
339 
340  Map<MatrixXd> eigen_dK(dK.matrix, dK.num_rows, dK.num_cols);
341 
342  // compute derivative wrt kernel parameter: dnlZ=sum(Q.*dK*scale)/2.0
343  result[i]=(eigen_Q.cwiseProduct(eigen_dK)*CMath::sq(m_scale)).sum()/2.0;
344  }
345 
346  return result;
347 }
348 
350  const TParameter* param)
351 {
352  // create eigen representation of alpha vector
353  Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
354 
355  REQUIRE(param, "Param not set\n");
356  SGVector<float64_t> result;
357  int64_t len=const_cast<TParameter *>(param)->m_datatype.get_num_elements();
358  result=SGVector<float64_t>(len);
359 
360  for (index_t i=0; i<result.vlen; i++)
361  {
363 
364  if (result.vlen==1)
366  else
368 
369  Map<VectorXd> eigen_dmu(dmu.vector, dmu.vlen);
370 
371  // compute derivative wrt mean parameter: dnlZ=-dmu'*alpha
372  result[i]=-eigen_dmu.dot(eigen_alpha);
373  }
374 
375  return result;
376 }
377 
378 #endif /* HAVE_EIGEN3 */
virtual const char * get_name() const =0
static void fill_vector(T *vec, int32_t len, T value)
Definition: SGVector.cpp:223
The Gaussian exact form inference method class.
virtual ELabelType get_label_type() const =0
virtual void update_parameter_hash()
Definition: SGObject.cpp:250
Class that models Gaussian likelihood.
Real Labels are real-valued labels.
SGVector< float64_t > m_alpha
The Inference Method base class.
virtual SGMatrix< float64_t > get_posterior_covariance()
int32_t index_t
Definition: common.h:62
Vector::Scalar dot(Vector a, Vector b)
Definition: Redux.h:56
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
real valued labels (e.g. for regression, classifier outputs)
Definition: LabelTypes.h:22
virtual const char * get_name() const
virtual SGVector< float64_t > get_alpha()
static T sq(T x)
Definition: Math.h:450
Definition: SGMatrix.h:20
virtual ELikelihoodModelType get_model_type() const
parameter struct
virtual int32_t get_num_vectors() const =0
#define REQUIRE(x,...)
Definition: SGIO.h:206
index_t num_cols
Definition: SGMatrix.h:378
virtual SGVector< float64_t > get_mean_vector(const CFeatures *features) const =0
static CExactInferenceMethod * obtain_from_generic(CInferenceMethod *inference)
An abstract class of the mean function.
Definition: MeanFunction.h:28
#define SG_REF(x)
Definition: SGObject.h:51
index_t num_rows
Definition: SGMatrix.h:376
index_t vlen
Definition: SGVector.h:494
SGMatrix< float64_t > m_L
virtual SGMatrix< float64_t > get_cholesky()
virtual SGVector< float64_t > get_derivative_wrt_mean(const TParameter *param)
double float64_t
Definition: common.h:50
virtual SGVector< float64_t > get_derivative_wrt_inference_method(const TParameter *param)
virtual SGVector< float64_t > get_diagonal_vector()
static T sum(T *vec, int32_t len)
Return sum(vec)
Definition: SGVector.h:354
static CGaussianLikelihood * obtain_from_generic(CLikelihoodModel *lik)
virtual SGVector< float64_t > get_parameter_derivative(const CFeatures *features, const TParameter *param, index_t index=-1)
Definition: MeanFunction.h:52
#define SG_UNREF(x)
Definition: SGObject.h:52
#define SG_DEBUG(...)
Definition: SGIO.h:107
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
virtual float64_t get_negative_log_marginal_likelihood()
The class Features is the base class of all feature objects.
Definition: Features.h:68
#define SG_SERROR(...)
Definition: SGIO.h:179
virtual SGVector< float64_t > get_derivative_wrt_kernel(const TParameter *param)
virtual SGMatrix< float64_t > get_parameter_gradient(const TParameter *param, index_t index=-1)
Definition: Kernel.h:850
static float64_t log(float64_t v)
Definition: Math.h:922
virtual void check_members() const
virtual SGVector< float64_t > get_derivative_wrt_likelihood_model(const TParameter *param)
virtual EInferenceType get_inference_type() const
The Kernel base class.
Definition: Kernel.h:158
virtual bool parameter_hash_changed()
Definition: SGObject.cpp:264
The Likelihood model base class.
SGMatrix< float64_t > m_ktrtr
CLikelihoodModel * m_model
virtual SGVector< float64_t > get_posterior_mean()
static const float64_t PI
Definition: Math.h:2055

SHOGUN Machine Learning Toolbox - Documentation