BayesOpt
gaussian_process_normal.cpp
1 /*
2 -------------------------------------------------------------------------
3  This file is part of BayesOpt, an efficient C++ library for
4  Bayesian optimization.
5 
6  Copyright (C) 2011-2015 Ruben Martinez-Cantin <rmcantin@unizar.es>
7 
8  BayesOpt is free software: you can redistribute it and/or modify it
9  under the terms of the GNU Affero General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  BayesOpt is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU Affero General Public License for more details.
17 
18  You should have received a copy of the GNU Affero General Public License
19  along with BayesOpt. If not, see <http://www.gnu.org/licenses/>.
20 ------------------------------------------------------------------------
21 */
22 
23 #include <cstdlib>
24 #include <boost/math/special_functions/fpclassify.hpp>
25 #include <boost/numeric/ublas/banded.hpp>
26 #include "ublas_trace.hpp"
27 #include "ublas_elementwise.hpp"
28 #include "gauss_distribution.hpp"
30 
31 namespace bayesopt
32 {
33 
34  namespace ublas = boost::numeric::ublas;
35 
36  GaussianProcessNormal::GaussianProcessNormal(size_t dim,
37  Parameters params,
38  const Dataset& data,
39  MeanModel& mean,
40  randEngine& eng):
41  HierarchicalGaussianProcess(dim,params,data,mean,eng),
42  mW0(params.mean.coef_mean.size()), mInvVarW(params.mean.coef_mean.size()),
43  mD(params.mean.coef_mean.size(),params.mean.coef_mean.size())
44  {
45  mSigma = params.sigma_s;
46  mW0 = params.mean.coef_mean;
47  for (size_t ii = 0; ii < params.mean.coef_std.size(); ++ii)
48  {
49  double varii = params.mean.coef_std[ii] * params.mean.coef_std[ii];
50  mInvVarW(ii) = 1/varii;
51  }
52  d_ = new GaussianDistribution(eng);
53  } // Constructor
54 
55 
56 
57  GaussianProcessNormal::~GaussianProcessNormal()
58  {
59  delete d_;
60  } // Default destructor
61 
62 
63  ProbabilityDistribution*
64  GaussianProcessNormal::prediction(const vectord &query)
65  {
66  const double kq = computeSelfCorrelation(query);
67  const vectord phi = mMean.getFeatures(query);
68 
69  vectord v = computeCrossCorrelation(query);
70 
71  inplace_solve(mL,v,ublas::lower_tag());
72 
73  vectord rq = phi - prod(v,mKF);
74 
75  vectord rho(rq);
76  inplace_solve(mD,rho,ublas::lower_tag());
77 
78  double yPred = inner_prod(phi,mWMap) + inner_prod(v,mVf);
79  double sPred = sqrt( mSigma * (kq - inner_prod(v,v)
80  + inner_prod(rho,rho)));
81 
82  if ((boost::math::isnan(yPred)) || (boost::math::isnan(sPred)))
83  {
84  throw std::runtime_error("Error in prediction. NaN found.");
85  }
86 
87 
88  d_->setMeanAndStd(yPred,sPred);
89  return d_;
90  }
91 
92 
94  {
95  matrixd KK = computeCorrMatrix();
96  const size_t n = KK.size1();
97  const size_t p = mMean.nFeatures();
98 
99  vectord v0 = mData.mY - prod(trans(mMean.mFeatM),mW0);
100  matrixd WW = zmatrixd(p,p); //TODO: diagonal matrix
101  utils::add_to_diagonal(WW,mInvVarW);
102  matrixd FW = prod(trans(mMean.mFeatM),WW);
103  KK += prod(FW,mMean.mFeatM);
104  matrixd BB(n,n);
106  inplace_solve(BB,v0,ublas::lower_tag());
107  double zz = inner_prod(v0,v0);
108 
109  double lik = 1/(2*mSigma) * zz;
110  lik += utils::log_trace(BB);
111  return lik;
112  }
113 
114 
115 
117  {
118  const size_t p = mMean.nFeatures();
119 
120  mKF = trans(mMean.mFeatM);
121  inplace_solve(mL,mKF,ublas::lower_tag());
122  //TODO: make one line
123  matrixd DD(p,p);
124  DD = prod(trans(mKF),mKF);
125  utils::add_to_diagonal(DD,mInvVarW);
127 
128  vectord vn = mData.mY;
129  inplace_solve(mL,vn,ublas::lower_tag());
130  mWMap = prod(mMean.mFeatM,vn) + utils::ublas_elementwise_prod(mInvVarW,mW0);
131  utils::cholesky_solve(mD,mWMap,ublas::lower());
132 
133  mVf = mData.mY - prod(trans(mMean.mFeatM),mWMap);
134  inplace_solve(mL,mVf,ublas::lower_tag());
135 
136  if (boost::math::isnan(mWMap(0)))
137  {
138  throw std::runtime_error("Error in precomputed prediction. NaN found.");
139  }
140  }
141 
142 } //namespace bayesopt
ProbabilityDistribution * prediction(const vectord &query)
Function that returns the prediction of the GP for a query point in the hypercube [0...
matrixd computeCorrMatrix()
Computes the Correlation (Kernel or Gram) matrix.
void setMeanAndStd(double mean, double std)
Sets the mean and std of the distribution.
void precomputePrediction()
Precompute some values of the prediction that do not depends on the query.
Namespace of the library interface.
Definition: using.dox:1
Elementwise operations for ublas vector/matrix.
matrixd mL
Cholesky decomposition of the Correlation matrix.
GaussianDistribution * d_
Predictive distributions.
Gaussian probability distribution.
v1 ublas_elementwise_prod(const v1 &a, const v2 &b)
Computes the elementwise product of two vectors or matrices.
double negativeLogLikelihood()
Computes the negative log likelihood and its gradient of the data.
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
void cholesky_solve(const TRIA &L, MATRIX &x, ublas::lower)
solve system L L^T x = b inplace
matrixd mFeatM
Value of the mean features at the input points.
vectord mY
Data values.
Definition: dataset.hpp:64
Gaussian process with normal prior on the parameters.