SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
EPInferenceMethod.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  * Based on ideas from GAUSSIAN PROCESS REGRESSION AND CLASSIFICATION Toolbox
10  * Copyright (C) 2005-2013 by Carl Edward Rasmussen & Hannes Nickisch under the
11  * FreeBSD License
12  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
13  */
14 
16 
17 #ifdef HAVE_EIGEN3
18 
24 
26 
27 using namespace shogun;
28 using namespace Eigen;
29 
30 // try to use previously allocated memory for SGVector
31 #define CREATE_SGVECTOR(vec, len, sg_type) \
32  { \
33  if (!vec.vector || vec.vlen!=len) \
34  vec=SGVector<sg_type>(len); \
35  }
36 
37 // try to use previously allocated memory for SGMatrix
38 #define CREATE_SGMATRIX(mat, rows, cols, sg_type) \
39  { \
40  if (!mat.matrix || mat.num_rows!=rows || mat.num_cols!=cols) \
41  mat=SGMatrix<sg_type>(rows, cols); \
42  }
43 
45 {
46  init();
47 }
48 
50  CMeanFunction* mean, CLabels* labels, CLikelihoodModel* model)
51  : CInferenceMethod(kernel, features, mean, labels, model)
52 {
53  init();
54 }
55 
57 {
58 }
59 
60 void CEPInferenceMethod::init()
61 {
62  m_max_sweep=15;
63  m_min_sweep=2;
64  m_tol=1e-4;
65 }
66 
68 {
70  update();
71 
72  return m_nlZ;
73 }
74 
76 {
78  update();
79 
81 }
82 
84 {
86  update();
87 
88  return SGMatrix<float64_t>(m_L);
89 }
90 
92 {
94  update();
95 
96  return SGVector<float64_t>(m_sttau);
97 }
98 
100 {
101  if (update_parameter_hash())
102  update();
103 
104  return SGVector<float64_t>(m_mu);
105 }
106 
108 {
109  if (update_parameter_hash())
110  update();
111 
112  return SGMatrix<float64_t>(m_Sigma);
113 }
114 
116 {
117  // update kernel and feature matrix
119 
120  // get number of labels (trainig examples)
122 
123  // try to use tilde values from previous call
124  if (m_ttau.vlen==n)
125  {
126  update_chol();
130  }
131 
132  // get mean vector
134 
135  // get and scale diagonal of the kernel matrix
137  ktrtr_diag.scale(CMath::sq(m_scale));
138 
139  // marginal likelihood for ttau = tnu = 0
141  mean, ktrtr_diag, m_labels));
142 
143  // use zero values if we have no better guess or it's better
144  if (m_ttau.vlen!=n || m_nlZ>nlZ0)
145  {
146  CREATE_SGVECTOR(m_ttau, n, float64_t);
147  m_ttau.zero();
148 
149  CREATE_SGVECTOR(m_sttau, n, float64_t);
150  m_sttau.zero();
151 
152  CREATE_SGVECTOR(m_tnu, n, float64_t);
153  m_tnu.zero();
154 
156 
157  // copy data manually, since we don't have appropriate method
158  for (index_t i=0; i<m_ktrtr.num_rows; i++)
159  for (index_t j=0; j<m_ktrtr.num_cols; j++)
160  m_Sigma(i,j)=m_ktrtr(i,j)*CMath::sq(m_scale);
161 
162  CREATE_SGVECTOR(m_mu, n, float64_t);
163  m_mu.zero();
164 
165  // set marginal likelihood
166  m_nlZ=nlZ0;
167  }
168 
169  // create vector of the random permutation
171 
172  // cavity tau and nu vectors
173  SGVector<float64_t> tau_n(n);
174  SGVector<float64_t> nu_n(n);
175 
176  // cavity mu and s2 vectors
177  SGVector<float64_t> mu_n(n);
178  SGVector<float64_t> s2_n(n);
179 
180  float64_t nlZ_old=CMath::INFTY;
181  uint32_t sweep=0;
182 
183  while ((CMath::abs(m_nlZ-nlZ_old)>m_tol && sweep<m_max_sweep) ||
184  sweep<m_min_sweep)
185  {
186  nlZ_old=m_nlZ;
187  sweep++;
188 
189  // shuffle random permutation
190  randperm.permute();
191 
192  for (index_t j=0; j<n; j++)
193  {
194  index_t i=randperm[j];
195 
196  // find cavity paramters
197  tau_n[i]=1.0/m_Sigma(i,i)-m_ttau[i];
198  nu_n[i]=m_mu[i]/m_Sigma(i,i)+mean[i]*tau_n[i]-m_tnu[i];
199 
200  // compute cavity mean and variance
201  mu_n[i]=nu_n[i]/tau_n[i];
202  s2_n[i]=1.0/tau_n[i];
203 
204  // get moments
205  float64_t mu=m_model->get_first_moment(mu_n, s2_n, m_labels, i);
206  float64_t s2=m_model->get_second_moment(mu_n, s2_n, m_labels, i);
207 
208  // save old value of ttau
209  float64_t ttau_old=m_ttau[i];
210 
211  // compute ttau and sqrt(ttau)
212  m_ttau[i]=CMath::max(1.0/s2-tau_n[i], 0.0);
213  m_sttau[i]=CMath::sqrt(m_ttau[i]);
214 
215  // compute tnu
216  m_tnu[i]=mu/s2-nu_n[i];
217 
218  // compute difference ds2=ttau_new-ttau_old
219  float64_t ds2=m_ttau[i]-ttau_old;
220 
221  // create eigen representation of Sigma, tnu and mu
222  Map<MatrixXd> eigen_Sigma(m_Sigma.matrix, m_Sigma.num_rows,
223  m_Sigma.num_cols);
224  Map<VectorXd> eigen_tnu(m_tnu.vector, m_tnu.vlen);
225  Map<VectorXd> eigen_mu(m_mu.vector, m_mu.vlen);
226 
227  VectorXd eigen_si=eigen_Sigma.col(i);
228 
229  // rank-1 update Sigma
230  eigen_Sigma=eigen_Sigma-ds2/(1.0+ds2*eigen_si(i))*eigen_si*
231  eigen_si.adjoint();
232 
233  // update mu
234  eigen_mu=eigen_Sigma*eigen_tnu;
235  }
236 
237  // update upper triangular factor (L^T) of Cholesky decomposition of
238  // matrix B, approximate posterior covariance and mean, negative
239  // marginal likelihood
240  update_chol();
244  }
245 
246  if (sweep==m_max_sweep && CMath::abs(m_nlZ-nlZ_old)>m_tol)
247  {
248  SG_ERROR("Maximum number (%d) of sweeps reached, but tolerance (%f) was "
249  "not yet reached. You can manually set maximum number of sweeps "
250  "or tolerance to fix this problem.\n", m_max_sweep, m_tol);
251  }
252 
253  // update vector alpha
254  update_alpha();
255 
256  // update matrices to compute derivatives
257  update_deriv();
258 }
259 
261 {
262  // create eigen representations kernel matrix, L^T, sqrt(ttau) and tnu
263  Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
264  Map<VectorXd> eigen_tnu(m_tnu.vector, m_tnu.vlen);
265  Map<VectorXd> eigen_sttau(m_sttau.vector, m_sttau.vlen);
266  Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
267 
268  // create shogun and eigen representation of the alpha vector
270  Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
271 
272  // solve LL^T * v = tS^(1/2) * K * tnu
273  VectorXd eigen_v=eigen_L.triangularView<Upper>().adjoint().solve(
274  eigen_sttau.cwiseProduct(eigen_K*CMath::sq(m_scale)*eigen_tnu));
275  eigen_v=eigen_L.triangularView<Upper>().solve(eigen_v);
276 
277  // compute alpha = (I - tS^(1/2) * B^(-1) * tS(1/2) * K) * tnu =
278  // tnu - tS(1/2) * (L^T)^(-1) * L^(-1) * tS^(1/2) * K * tnu =
279  // tnu - tS(1/2) * v
280  eigen_alpha=eigen_tnu-eigen_sttau.cwiseProduct(eigen_v);
281 }
282 
284 {
285  // create eigen representations of kernel matrix and sqrt(ttau)
286  Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
287  Map<VectorXd> eigen_sttau(m_sttau.vector, m_sttau.vlen);
288 
289  // create shogun and eigen representation of the upper triangular factor
290  // (L^T) of the Cholesky decomposition of the matrix B
292  Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
293 
294  // compute upper triangular factor L^T of the Cholesky decomposion of the
295  // matrix: B = tS^(1/2) * K * tS^(1/2) + I
296  LLT<MatrixXd> eigen_chol((eigen_sttau*eigen_sttau.adjoint()).cwiseProduct(
297  eigen_K*CMath::sq(m_scale))+
298  MatrixXd::Identity(m_L.num_rows, m_L.num_cols));
299 
300  eigen_L=eigen_chol.matrixU();
301 }
302 
304 {
305  // create eigen representations of kernel matrix, L^T matrix and sqrt(ttau)
306  Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
307  Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
308  Map<VectorXd> eigen_sttau(m_sttau.vector, m_sttau.vlen);
309 
310  // create shogun and eigen representation of the approximate covariance
311  // matrix
313  Map<MatrixXd> eigen_Sigma(m_Sigma.matrix, m_Sigma.num_rows, m_Sigma.num_cols);
314 
315  // compute V = L^(-1) * tS^(1/2) * K, using upper triangular factor L^T
316  MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
317  eigen_sttau.asDiagonal()*eigen_K*CMath::sq(m_scale));
318 
319  // compute covariance matrix of the posterior:
320  // Sigma = K - K * tS^(1/2) * (L * L^T)^(-1) * tS^(1/2) * K =
321  // K - (K * tS^(1/2)) * (L^T)^(-1) * L^(-1) * tS^(1/2) * K =
322  // K - (tS^(1/2) * K)^T * (L^(-1))^T * L^(-1) * tS^(1/2) * K = K - V^T * V
323  eigen_Sigma=eigen_K*CMath::sq(m_scale)-eigen_V.adjoint()*eigen_V;
324 }
325 
327 {
328  // create eigen representation of posterior covariance matrix and tnu
329  Map<MatrixXd> eigen_Sigma(m_Sigma.matrix, m_Sigma.num_rows, m_Sigma.num_cols);
330  Map<VectorXd> eigen_tnu(m_tnu.vector, m_tnu.vlen);
331 
332  // create shogun and eigen representation of the approximate mean vector
333  CREATE_SGVECTOR(m_mu, m_tnu.vlen, float64_t);
334  Map<VectorXd> eigen_mu(m_mu.vector, m_mu.vlen);
335 
336  // compute mean vector of the approximate posterior: mu = Sigma * tnu
337  eigen_mu=eigen_Sigma*eigen_tnu;
338 }
339 
341 {
342  // create eigen representation of Sigma, L, mu, tnu, ttau
343  Map<MatrixXd> eigen_Sigma(m_Sigma.matrix, m_Sigma.num_rows, m_Sigma.num_cols);
344  Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
345  Map<VectorXd> eigen_mu(m_mu.vector, m_mu.vlen);
346  Map<VectorXd> eigen_tnu(m_tnu.vector, m_tnu.vlen);
347  Map<VectorXd> eigen_ttau(m_ttau.vector, m_ttau.vlen);
348 
349  // get and create eigen representation of the mean vector
351  Map<VectorXd> eigen_m(m.vector, m.vlen);
352 
353  // compute vector of cavity parameter tau_n
354  VectorXd eigen_tau_n=(VectorXd::Ones(m_ttau.vlen)).cwiseQuotient(
355  eigen_Sigma.diagonal())-eigen_ttau;
356 
357  // compute vector of cavity parameter nu_n
358  VectorXd eigen_nu_n=eigen_mu.cwiseQuotient(eigen_Sigma.diagonal())-
359  eigen_tnu+eigen_m.cwiseProduct(eigen_tau_n);
360 
361  // compute cavity mean: mu_n=nu_n/tau_n
362  SGVector<float64_t> mu_n(m_ttau.vlen);
363  Map<VectorXd> eigen_mu_n(mu_n.vector, mu_n.vlen);
364 
365  eigen_mu_n=eigen_nu_n.cwiseQuotient(eigen_tau_n);
366 
367  // compute cavity variance: s2_n=1.0/tau_n
368  SGVector<float64_t> s2_n(m_ttau.vlen);
369  Map<VectorXd> eigen_s2_n(s2_n.vector, s2_n.vlen);
370 
371  eigen_s2_n=(VectorXd::Ones(m_ttau.vlen)).cwiseQuotient(eigen_tau_n);
372 
374  m_model->get_log_zeroth_moments(mu_n, s2_n, m_labels));
375 
376  // compute nlZ_part1=sum(log(diag(L)))-sum(lZ)-tnu'*Sigma*tnu/2
377  float64_t nlZ_part1=eigen_L.diagonal().array().log().sum()-lZ-
378  (eigen_tnu.adjoint()*eigen_Sigma).dot(eigen_tnu)/2.0;
379 
380  // compute nlZ_part2=sum(tnu.^2./(tau_n+ttau))/2-sum(log(1+ttau./tau_n))/2
381  float64_t nlZ_part2=(eigen_tnu.array().square()/
382  (eigen_tau_n+eigen_ttau).array()).sum()/2.0-(1.0+eigen_ttau.array()/
383  eigen_tau_n.array()).log().sum()/2.0;
384 
385  // compute nlZ_part3=-(nu_n-m.*tau_n)'*((ttau./tau_n.*(nu_n-m.*tau_n)-2*tnu)
386  // ./(ttau+tau_n))/2
387  float64_t nlZ_part3=-(eigen_nu_n-eigen_m.cwiseProduct(eigen_tau_n)).dot(
388  ((eigen_ttau.array()/eigen_tau_n.array()*(eigen_nu_n.array()-
389  eigen_m.array()*eigen_tau_n.array())-2*eigen_tnu.array())/
390  (eigen_ttau.array()+eigen_tau_n.array())).matrix())/2.0;
391 
392  // compute nlZ=nlZ_part1+nlZ_part2+nlZ_part3
393  m_nlZ=nlZ_part1+nlZ_part2+nlZ_part3;
394 }
395 
397 {
398  // create eigen representation of L, sstau, alpha
399  Map<MatrixXd> eigen_L(m_L.matrix, m_L.num_rows, m_L.num_cols);
400  Map<VectorXd> eigen_sttau(m_sttau.vector, m_sttau.vlen);
401  Map<VectorXd> eigen_alpha(m_alpha.vector, m_alpha.vlen);
402 
403  // create shogun and eigen representation of F
405  Map<MatrixXd> eigen_F(m_F.matrix, m_F.num_rows, m_F.num_cols);
406 
407  // solve L*L^T * V = diag(sqrt(ttau))
408  MatrixXd V=eigen_L.triangularView<Upper>().adjoint().solve(
409  MatrixXd(eigen_sttau.asDiagonal()));
410  V=eigen_L.triangularView<Upper>().solve(V);
411 
412  // compute F=alpha*alpha'-repmat(sW,1,n).*solve_chol(L,diag(sW))
413  eigen_F=eigen_alpha*eigen_alpha.adjoint()-eigen_sttau.asDiagonal()*V;
414 }
415 
417  const TParameter* param)
418 {
419  REQUIRE(!strcmp(param->m_name, "scale"), "Can't compute derivative of "
420  "the nagative log marginal likelihood wrt %s.%s parameter\n",
421  get_name(), param->m_name)
422 
423  Map<MatrixXd> eigen_K(m_ktrtr.matrix, m_ktrtr.num_rows, m_ktrtr.num_cols);
424  Map<MatrixXd> eigen_F(m_F.matrix, m_F.num_rows, m_F.num_cols);
425 
426  SGVector<float64_t> result(1);
427 
428  // compute derivative wrt kernel scale: dnlZ=-sum(F.*K*scale*2)/2
429  result[0]=-(eigen_F.cwiseProduct(eigen_K)*m_scale*2.0).sum()/2.0;
430 
431  return result;
432 }
433 
435  const TParameter* param)
436 {
438  return SGVector<float64_t>();
439 }
440 
442  const TParameter* param)
443 {
444  // create eigen representation of the matrix Q
445  Map<MatrixXd> eigen_F(m_F.matrix, m_F.num_rows, m_F.num_cols);
446 
447  SGVector<float64_t> result;
448 
449  if (param->m_datatype.m_ctype==CT_VECTOR ||
450  param->m_datatype.m_ctype==CT_SGVECTOR)
451  {
453  "Length of the parameter %s should not be NULL\n", param->m_name)
454  result=SGVector<float64_t>(*(param->m_datatype.m_length_y));
455  }
456  else
457  {
458  result=SGVector<float64_t>(1);
459  }
460 
461  for (index_t i=0; i<result.vlen; i++)
462  {
464 
465  if (result.vlen==1)
466  dK=m_kernel->get_parameter_gradient(param);
467  else
468  dK=m_kernel->get_parameter_gradient(param, i);
469 
470  Map<MatrixXd> eigen_dK(dK.matrix, dK.num_rows, dK.num_cols);
471 
472  // compute derivative wrt kernel parameter: dnlZ=-sum(F.*dK*scale^2)/2.0
473  result[i]=-(eigen_F.cwiseProduct(eigen_dK)*CMath::sq(m_scale)).sum()/2.0;
474  }
475 
476  return result;
477 }
478 
480  const TParameter* param)
481 {
483  return SGVector<float64_t>();
484 }
485 
486 #endif /* HAVE_EIGEN3 */
virtual SGVector< float64_t > get_diagonal_vector()
virtual SGVector< float64_t > get_alpha()
SGVector< float64_t > m_alpha
The Inference Method base class.
int32_t index_t
Definition: common.h:60
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:35
static const float64_t INFTY
infinity
Definition: Math.h:1330
virtual SGMatrix< float64_t > get_posterior_covariance()
virtual int32_t get_num_labels() const =0
static void permute(T *vec, int32_t n)
Definition: SGVector.cpp:710
static T sq(T x)
x^2
Definition: Math.h:239
parameter struct
Definition: Parameter.h:26
#define SG_ERROR(...)
Definition: SGIO.h:131
virtual SGVector< float64_t > get_log_zeroth_moments(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab) const =0
#define REQUIRE(x,...)
Definition: SGIO.h:208
#define SG_NOTIMPLEMENTED
Definition: SGIO.h:141
index_t num_cols
Definition: SGMatrix.h:303
virtual SGVector< float64_t > get_mean_vector(const CFeatures *features) const =0
#define CREATE_SGMATRIX(mat, rows, cols, sg_type)
virtual float64_t get_negative_log_marginal_likelihood()
virtual float64_t get_second_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const =0
An abstract class of the mean function.
Definition: MeanFunction.h:26
void scale(T alpha)
scale vector inplace
Definition: SGVector.cpp:1200
virtual SGVector< float64_t > get_derivative_wrt_kernel(const TParameter *param)
index_t num_rows
Definition: SGMatrix.h:301
virtual bool update_parameter_hash()
Definition: SGObject.cpp:187
virtual SGVector< float64_t > get_posterior_mean()
SGVector< T > get_diagonal_vector() const
Definition: SGMatrix.cpp:978
TSGDataType m_datatype
Definition: Parameter.h:156
index_t vlen
Definition: SGVector.h:706
SGMatrix< float64_t > m_L
virtual SGVector< float64_t > get_derivative_wrt_mean(const TParameter *param)
static SGVector< T > randperm_vec(int32_t n)
Definition: SGVector.cpp:664
double float64_t
Definition: common.h:48
#define CREATE_SGVECTOR(vec, len, sg_type)
static T sum(T *vec, int32_t len)
return sum(vec)
Definition: SGVector.h:506
static T max(T a, T b)
return the maximum of two integers
Definition: Math.h:160
virtual SGVector< float64_t > get_derivative_wrt_likelihood_model(const TParameter *param)
index_t * m_length_y
Definition: DataType.h:77
EContainerType m_ctype
Definition: DataType.h:70
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:16
virtual SGVector< float64_t > get_derivative_wrt_inference_method(const TParameter *param)
The class Features is the base class of all feature objects.
Definition: Features.h:62
virtual SGMatrix< float64_t > get_parameter_gradient(const TParameter *param, index_t index=-1)
Definition: Kernel.h:574
The Kernel base class.
Definition: Kernel.h:150
virtual SGMatrix< float64_t > get_cholesky()
static float32_t sqrt(float32_t x)
x^0.5
Definition: Math.h:245
virtual float64_t get_first_moment(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab, index_t i) const =0
virtual const char * get_name() const
The Likelihood model base class.
SGMatrix< float64_t > m_ktrtr
CLikelihoodModel * m_model
static T abs(T a)
return the absolute value of a number
Definition: Math.h:179

SHOGUN Machine Learning Toolbox - Documentation