Kstars

gaussian_process.cpp
Go to the documentation of this file.
1/*
2 SPDX-FileCopyrightText: 2014-2017 Max Planck Society.
3 All rights reserved.
4
5 SPDX-License-Identifier: BSD-3-Clause
6*/
7
8/**
9 * @file
10 * @date 2014-2017
11 * @copyright Max Planck Society
12 *
13 * @author Edgar D. Klenske <edgar.klenske@tuebingen.mpg.de>
14 * @author Stephan Wenninger <stephan.wenninger@tuebingen.mpg.de>
15 * @author Raffi Enficiaud <raffi.enficiaud@tuebingen.mpg.de>
16 *
17 * @brief The GP class implements the Gaussian Process functionality.
18 */
19
20#include <cstdint>
21
22#include "gaussian_process.h"
23#include "math_tools.h"
24#include "covariance_functions.h"
25
26// A functor for special orderings
27struct covariance_ordering
28{
29 covariance_ordering(Eigen::VectorXd const &cov) : covariance_(cov) {}
30 bool operator()(int a, int b) const
31 {
32 return (covariance_[a] > covariance_[b]);
33 }
34 Eigen::VectorXd const &covariance_;
35};
36
37GP::GP() : covFunc_(nullptr), // initialize pointer to null
38 covFuncProj_(nullptr), // initialize pointer to null
39 data_loc_(Eigen::VectorXd()),
40 data_out_(Eigen::VectorXd()),
41 data_var_(Eigen::VectorXd()),
42 gram_matrix_(Eigen::MatrixXd()),
43 alpha_(Eigen::VectorXd()),
44 chol_gram_matrix_(Eigen::LDLT<Eigen::MatrixXd>()),
45 log_noise_sd_(-1E20),
46 use_explicit_trend_(false),
47 feature_vectors_(Eigen::MatrixXd()),
48 feature_matrix_(Eigen::MatrixXd()),
49 chol_feature_matrix_(Eigen::LDLT<Eigen::MatrixXd>()),
50 beta_(Eigen::VectorXd())
51{ }
52
53GP::GP(const covariance_functions::CovFunc &covFunc) :
54 covFunc_(covFunc.clone()),
55 covFuncProj_(nullptr),
56 data_loc_(Eigen::VectorXd()),
57 data_out_(Eigen::VectorXd()),
58 data_var_(Eigen::VectorXd()),
59 gram_matrix_(Eigen::MatrixXd()),
60 alpha_(Eigen::VectorXd()),
61 chol_gram_matrix_(Eigen::LDLT<Eigen::MatrixXd>()),
62 log_noise_sd_(-1E20),
63 use_explicit_trend_(false),
64 feature_vectors_(Eigen::MatrixXd()),
65 feature_matrix_(Eigen::MatrixXd()),
66 chol_feature_matrix_(Eigen::LDLT<Eigen::MatrixXd>()),
67 beta_(Eigen::VectorXd())
68{ }
69
70GP::GP(const double noise_variance,
71 const covariance_functions::CovFunc &covFunc) :
72 covFunc_(covFunc.clone()),
73 covFuncProj_(nullptr),
74 data_loc_(Eigen::VectorXd()),
75 data_out_(Eigen::VectorXd()),
76 data_var_(Eigen::VectorXd()),
77 gram_matrix_(Eigen::MatrixXd()),
78 alpha_(Eigen::VectorXd()),
79 chol_gram_matrix_(Eigen::LDLT<Eigen::MatrixXd>()),
80 log_noise_sd_(std::log(noise_variance)),
81 use_explicit_trend_(false),
82 feature_vectors_(Eigen::MatrixXd()),
83 feature_matrix_(Eigen::MatrixXd()),
84 chol_feature_matrix_(Eigen::LDLT<Eigen::MatrixXd>()),
85 beta_(Eigen::VectorXd())
86{ }
87
88GP::~GP()
89{
90 delete this->covFunc_; // tidy up since we are responsible for the covFunc.
91 delete this->covFuncProj_; // tidy up since we are responsible for the covFuncProj.
92}
93
94GP::GP(const GP &that) :
95 covFunc_(nullptr), // initialize to nullptr, clone later
96 covFuncProj_(nullptr), // initialize to nullptr, clone later
97 data_loc_(that.data_loc_),
98 data_out_(that.data_out_),
99 data_var_(that.data_var_),
100 gram_matrix_(that.gram_matrix_),
101 alpha_(that.alpha_),
102 chol_gram_matrix_(that.chol_gram_matrix_),
103 log_noise_sd_(that.log_noise_sd_),
104 use_explicit_trend_(that.use_explicit_trend_),
105 feature_vectors_(that.feature_vectors_),
106 feature_matrix_(that.feature_matrix_),
107 chol_feature_matrix_(that.chol_feature_matrix_),
108 beta_(that.beta_)
109{
110 covFunc_ = that.covFunc_->clone();
111 covFuncProj_ = that.covFuncProj_->clone();
112}
113
114bool GP::setCovarianceFunction(const covariance_functions::CovFunc &covFunc)
115{
116 // can only set the covariance function if training dataset is empty
117 if (data_loc_.size() != 0 || data_out_.size() != 0)
118 return false;
119 delete covFunc_; // initialized to zero, so delete is safe
120 covFunc_ = covFunc.clone();
121
122 return true;
123}
124
125void GP::enableOutputProjection(const covariance_functions::CovFunc &covFuncProj)
126{
127 delete covFuncProj_;// initialized to zero, so delete is safe
128 covFuncProj_ = covFuncProj.clone();
129}
130
131void GP::disableOutputProjection()
132{
133 delete covFuncProj_; // initialized to zero, so delete is safe
134 covFuncProj_ = nullptr;
135}
136
137GP &GP::operator=(const GP &that)
138{
139 if (this != &that)
140 {
141 covariance_functions::CovFunc* temp = covFunc_; // store old pointer...
142 covFunc_ = that.covFunc_->clone(); // ... first clone ...
143 delete temp; // ... and then delete.
144
145 // copy the rest
146 data_loc_ = that.data_loc_;
147 data_out_ = that.data_out_;
148 data_var_ = that.data_var_;
149 gram_matrix_ = that.gram_matrix_;
150 alpha_ = that.alpha_;
151 chol_gram_matrix_ = that.chol_gram_matrix_;
152 log_noise_sd_ = that.log_noise_sd_;
153 }
154 return *this;
155}
156
157Eigen::VectorXd GP::drawSample(const Eigen::VectorXd &locations) const
158{
159 return drawSample(locations,
160 math_tools::generate_normal_random_matrix(locations.rows(), locations.cols()));
161}
162
163Eigen::VectorXd GP::drawSample(const Eigen::VectorXd &locations,
164 const Eigen::VectorXd &random_vector) const
165{
166 Eigen::MatrixXd prior_covariance;
167 Eigen::MatrixXd kernel_matrix;
168 Eigen::LLT<Eigen::MatrixXd> chol_kernel_matrix;
169 Eigen::MatrixXd samples;
170
171 // we need the prior covariance for both, prior and posterior samples.
172 prior_covariance = covFunc_->evaluate(locations, locations);
173 kernel_matrix = prior_covariance;
174
175 if (gram_matrix_.cols() == 0) // no data, i.e. only a prior
176 {
177 kernel_matrix = prior_covariance + JITTER * Eigen::MatrixXd::Identity(
178 prior_covariance.rows(), prior_covariance.cols());
179 }
180 else // we have some data
181 {
182 Eigen::MatrixXd mixed_covariance;
183 mixed_covariance = covFunc_->evaluate(locations, data_loc_);
184 Eigen::MatrixXd posterior_covariance;
185 posterior_covariance = prior_covariance - mixed_covariance *
186 (chol_gram_matrix_.solve(mixed_covariance.transpose()));
187 kernel_matrix = posterior_covariance + JITTER * Eigen::MatrixXd::Identity(
188 posterior_covariance.rows(), posterior_covariance.cols());
189 }
190 chol_kernel_matrix = kernel_matrix.llt();
191
192 // Draw sample: s = chol(K)*x, where x is a random vector
193 samples = chol_kernel_matrix.matrixL() * random_vector;
194
195 // add the measurement noise on return
196 return samples + std::exp(log_noise_sd_) *
197 math_tools::generate_normal_random_matrix(samples.rows(), samples.cols());
198}
199
200void GP::infer()
201{
202 assert(data_loc_.rows() > 0 && "Error: the GP is not yet initialized!");
203
204 // The data covariance matrix
205 Eigen::MatrixXd data_cov = covFunc_->evaluate(data_loc_, data_loc_);
206
207 // swapping in the Gram matrix is faster than directly assigning it
208 gram_matrix_.swap(data_cov); // store the new data_cov as gram matrix
209 if (data_var_.rows() == 0) // homoscedastic
210 {
211 gram_matrix_ += (std::exp(2 * log_noise_sd_) + JITTER) *
212 Eigen::MatrixXd::Identity(gram_matrix_.rows(), gram_matrix_.cols());
213 }
214 else // heteroscedastic
215 {
216 gram_matrix_ += data_var_.asDiagonal();
217 }
218
219 // compute the Cholesky decomposition of the Gram matrix
220 chol_gram_matrix_ = gram_matrix_.ldlt();
221
222 // pre-compute the alpha, which is the solution of the chol to the data
223 alpha_ = chol_gram_matrix_.solve(data_out_);
224
225 if (use_explicit_trend_)
226 {
227 feature_vectors_ = Eigen::MatrixXd(2, data_loc_.rows());
228 // precompute necessary matrices for the explicit trend function
229 feature_vectors_.row(0) = Eigen::MatrixXd::Ones(1, data_loc_.rows()); // instead of pow(0)
230 feature_vectors_.row(1) = data_loc_.array(); // instead of pow(1)
231
232 feature_matrix_ = feature_vectors_ * chol_gram_matrix_.solve(feature_vectors_.transpose());
233 chol_feature_matrix_ = feature_matrix_.ldlt();
234
235 beta_ = chol_feature_matrix_.solve(feature_vectors_) * alpha_;
236 }
237}
238
239void GP::infer(const Eigen::VectorXd &data_loc,
240 const Eigen::VectorXd &data_out,
241 const Eigen::VectorXd &data_var /* = EigenVectorXd() */)
242{
243 data_loc_ = data_loc;
244 data_out_ = data_out;
245 if (data_var.rows() > 0)
246 {
247 data_var_ = data_var;
248 }
249 infer(); // updates the Gram matrix and its Cholesky decomposition
250}
251
252void GP::inferSD(const Eigen::VectorXd &data_loc,
253 const Eigen::VectorXd &data_out,
254 const int n, const Eigen::VectorXd &data_var /* = EigenVectorXd() */,
255 const double prediction_point /*= std::numeric_limits<double>::quiet_NaN()*/)
256{
257 Eigen::VectorXd covariance;
258
259 Eigen::VectorXd prediction_loc(1);
260 if ( math_tools::isNaN(prediction_point) )
261 {
262 // if none given, use the last datapoint as prediction reference
263 prediction_loc = data_loc.tail(1);
264 }
265 else
266 {
267 prediction_loc << prediction_point;
268 }
269
270 // calculate covariance between data and prediction point for point selection
271 covariance = covFunc_->evaluate(data_loc, prediction_loc);
272
273 // generate index vector
274 std::vector<int> index(covariance.size(), 0);
275 for (size_t i = 0 ; i != index.size() ; i++)
276 {
277 index[i] = i;
278 }
279
280 // sort indices with respect to covariance value
281 std::sort(index.begin(), index.end(),
282 covariance_ordering(covariance)
283 );
284
285 bool use_var = data_var.rows() > 0; // true means heteroscedastic noise
286
287 if (n < data_loc.rows())
288 {
289 std::vector<double> loc_arr(n);
290 std::vector<double> out_arr(n);
291 std::vector<double> var_arr(n);
292
293 for (int i = 0; i < n; ++i)
294 {
295 loc_arr[i] = data_loc[index[i]];
296 out_arr[i] = data_out[index[i]];
297 if (use_var)
298 {
299 var_arr[i] = data_var[index[i]];
300 }
301 }
302
303 data_loc_ = Eigen::Map<Eigen::VectorXd>(loc_arr.data(), n, 1);
304 data_out_ = Eigen::Map<Eigen::VectorXd>(out_arr.data(), n, 1);
305 if (use_var)
306 {
307 data_var_ = Eigen::Map<Eigen::VectorXd>(var_arr.data(), n, 1);
308 }
309 }
310 else // we can use all points and don't neet to select
311 {
312 data_loc_ = data_loc;
313 data_out_ = data_out;
314 if (use_var)
315 {
316 data_var_ = data_var;
317 }
318 }
319 infer();
320}
321
322void GP::clearData()
323{
324 gram_matrix_ = Eigen::MatrixXd();
325 chol_gram_matrix_ = Eigen::LDLT<Eigen::MatrixXd>();
326 data_loc_ = Eigen::VectorXd();
327 data_out_ = Eigen::VectorXd();
328}
329
330Eigen::VectorXd GP::predict(const Eigen::VectorXd &locations, Eigen::VectorXd* variances /*=nullptr*/) const
331{
332 // The prior covariance matrix (evaluated on test points)
333 Eigen::MatrixXd prior_cov = covFunc_->evaluate(locations, locations);
334
335 if (data_loc_.rows() == 0) // check if the data is empty
336 {
337 if (variances != nullptr)
338 {
339 (*variances) = prior_cov.diagonal();
340 }
341 return Eigen::VectorXd::Zero(locations.size());
342 }
343 else
344 {
345 // Calculate mixed covariance matrix (test and data points)
346 Eigen::MatrixXd mixed_cov = covFunc_->evaluate(locations, data_loc_);
347
348 // Calculate feature matrix for linear feature
349 Eigen::MatrixXd phi(2, locations.rows());
350 if (use_explicit_trend_)
351 {
352 phi.row(0) = Eigen::MatrixXd::Ones(1, locations.rows()); // instead of pow(0)
353 phi.row(1) = locations.array(); // instead of pow(1)
354
355 return predict(prior_cov, mixed_cov, phi, variances);
356 }
357 return predict(prior_cov, mixed_cov, Eigen::MatrixXd(), variances);
358 }
359}
360
361Eigen::VectorXd GP::predictProjected(const Eigen::VectorXd &locations, Eigen::VectorXd* variances /*=nullptr*/) const
362{
363 // use the suitable covariance function, depending on whether an
364 // output projection is used or not.
365 covariance_functions::CovFunc* covFunc = nullptr;
366 if (covFuncProj_ == nullptr)
367 {
368 covFunc = covFunc_;
369 }
370 else
371 {
372 covFunc = covFuncProj_;
373 }
374
375 assert(covFunc != nullptr);
376
377 // The prior covariance matrix (evaluated on test points)
378 Eigen::MatrixXd prior_cov = covFunc->evaluate(locations, locations);
379
380 if (data_loc_.rows() == 0) // check if the data is empty
381 {
382 if (variances != nullptr)
383 {
384 (*variances) = prior_cov.diagonal();
385 }
386 return Eigen::VectorXd::Zero(locations.size());
387 }
388 else
389 {
390 // The mixed covariance matrix (test and data points)
391 Eigen::MatrixXd mixed_cov = covFunc->evaluate(locations, data_loc_);
392
393 Eigen::MatrixXd phi(2, locations.rows());
394 if (use_explicit_trend_)
395 {
396 // calculate the feature vectors for linear regression
397 phi.row(0) = Eigen::MatrixXd::Ones(1, locations.rows()); // instead of pow(0)
398 phi.row(1) = locations.array(); // instead of pow(1)
399
400 return predict(prior_cov, mixed_cov, phi, variances);
401 }
402 return predict(prior_cov, mixed_cov, Eigen::MatrixXd(), variances);
403 }
404}
405
406Eigen::VectorXd GP::predict(const Eigen::MatrixXd &prior_cov, const Eigen::MatrixXd &mixed_cov,
407 const Eigen::MatrixXd &phi /*=Eigen::MatrixXd()*/, Eigen::VectorXd* variances /*=nullptr*/)
408const
409{
410
411 // calculate GP mean from precomputed alpha vector
412 Eigen::VectorXd m = mixed_cov * alpha_;
413
414 // precompute K^{-1} * mixed_cov
415 Eigen::MatrixXd gamma = chol_gram_matrix_.solve(mixed_cov.transpose());
416
417 Eigen::MatrixXd R;
418
419 // include fixed-features in the calculations
420 if (use_explicit_trend_)
421 {
422 R = phi - feature_vectors_ * gamma;
423
424 m += R.transpose() * beta_;
425 }
426
427 if (variances != nullptr)
428 {
429 // calculate GP variance
430 Eigen::MatrixXd v = prior_cov - mixed_cov * gamma;
431
432 // include fixed-features in the calculations
433 if (use_explicit_trend_)
434 {
435 assert(R.size() > 0);
436 v += R.transpose() * chol_feature_matrix_.solve(R);
437 }
438
439 (*variances) = v.diagonal();
440 }
441 return m;
442}
443
444void GP::setHyperParameters(const Eigen::VectorXd &hyperParameters)
445{
446 assert(hyperParameters.rows() == covFunc_->getParameterCount() + covFunc_->getExtraParameterCount() + 1 &&
447 "Wrong number of hyperparameters supplied to setHyperParameters()!");
448 log_noise_sd_ = hyperParameters[0];
449 covFunc_->setParameters(hyperParameters.segment(1, covFunc_->getParameterCount()));
450 covFunc_->setExtraParameters(hyperParameters.tail(covFunc_->getExtraParameterCount()));
451 if (data_loc_.rows() > 0)
452 {
453 infer();
454 }
455
456 // if the projection kernel is set, set the parameters there as well.
457 if (covFuncProj_ != nullptr)
458 {
459 covFuncProj_->setParameters(hyperParameters.segment(1, covFunc_->getParameterCount()));
460 covFuncProj_->setExtraParameters(hyperParameters.tail(covFunc_->getExtraParameterCount()));
461 }
462}
463
464Eigen::VectorXd GP::getHyperParameters() const
465{
466 Eigen::VectorXd hyperParameters(covFunc_->getParameterCount() + covFunc_->getExtraParameterCount() + 1);
467 hyperParameters << log_noise_sd_, covFunc_->getParameters(), covFunc_->getExtraParameters();
468 return hyperParameters;
469}
470
471void GP::enableExplicitTrend()
472{
473 use_explicit_trend_ = true;
474}
475
476void GP::disableExplicitTrend()
477{
478 use_explicit_trend_ = false;
479}
covariance_functions::CovFunc
Base class definition for covariance functions.
Definition covariance_functions.h:36
covariance_functions::CovFunc::evaluate
virtual Eigen::MatrixXd evaluate(const Eigen::VectorXd &x1, const Eigen::VectorXd &x2)=0
covariance_functions::CovFunc::setParameters
virtual void setParameters(const Eigen::VectorXd &params)=0
Method to set the hyper-parameters.
covariance_functions::CovFunc::getParameterCount
virtual int getParameterCount() const =0
Returns the number of hyper-parameters.
covariance_functions::CovFunc::getParameters
virtual const Eigen::VectorXd & getParameters() const =0
Returns the hyper-parameters.
covariance_functions::CovFunc::clone
virtual CovFunc * clone() const =0
Produces a clone to be able to copy the object.
covariance_functions.h
The file holds the covariance functions that can be used with the GP class.
gaussian_process.h
The GP class implements the Gaussian Process functionality.
math_tools::generate_normal_random_matrix
Eigen::MatrixXd generate_normal_random_matrix(const size_t n, const size_t m)
Definition src/math_tools.cpp:138
math_tools.h
Provides mathematical tools needed for the Gaussian process toolbox.
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri Jul 26 2024 11:59:51 by doxygen 1.11.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.