-
Notifications
You must be signed in to change notification settings - Fork 1
/
GPpredmean.cpp
93 lines (77 loc) · 2.98 KB
/
GPpredmean.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
// To run: library(Rcpp); sourceCpp("backsolve_eigen_test.cpp")
// Backsolve triangular lower
// [[Rcpp::depends(RcppEigen)]]
//#include <Rcpp.h>
#include <RcppEigen.h>
using namespace Eigen;
// [[Rcpp::export]]
Rcpp::List GPpredmean_rcpp(const Eigen::Map<Eigen::MatrixXd>& K,
const Eigen::Map<Eigen::MatrixXd>& Kstar,
const Eigen::Map<Eigen::MatrixXd>& Kstarstar,
const Eigen::Map<Eigen::VectorXd>& y,
const Eigen::Map<Eigen::MatrixXd>& H,
const Eigen::Map<Eigen::MatrixXd>& Hstar) {
int m = K.rows(); // m observations
int n = Kstarstar.cols(); // n targets
int k = H.rows(); // k basis function (for mean function)
// int nn = Hstar.cols();
Eigen::MatrixXd L(m, m); // chol(K)
Eigen::MatrixXd d(m, k);
Eigen::MatrixXd d2(k, m);
Eigen::MatrixXd HKHT(k, k); // H K^-1 H^T
Eigen::MatrixXd L2(k, k); // chol(HKHT) = chol( H K^-1 H^T )
Eigen::VectorXd a(m);
Eigen::VectorXd w(k);
Eigen::MatrixXd b(m, n);
Eigen::VectorXd LB(k);
Eigen::VectorXd B(k);
Eigen::MatrixXd R(k, n);
Eigen::MatrixXd e(k, n);
Eigen::MatrixXd btb(n, n);
Eigen::MatrixXd ete(n, n);
Eigen::VectorXd M(n); // Mean
Eigen::MatrixXd C(n, n); // covariance
Eigen::VectorXd logLik1(1);
Eigen::VectorXd logLik2(m);
Eigen::VectorXd logLik3(k);
// cholesky factor L such that K = LL^T
//L = K.adjoint().llt().matrixL();
L = K.llt().matrixL();
// d=L^-1 H^T (forward substitiution)
d = L.triangularView<Lower>().solve(H.adjoint());
//dT = L.triangularView<Upper>().solve(H);
// HKHT = H K^-1 H^T = d^T d
HKHT = MatrixXd(k,k).setZero().selfadjointView<Lower>().rankUpdate(d.adjoint());
// cholesky factor L2 such that HKHT = L2 L2^T
L2 = HKHT.adjoint().llt().matrixL();
//
a = L.triangularView<Lower>().solve(y);
//
b = L.triangularView<Lower>().solve(Kstar);
// Compute B = (H K⁻1 H^T)^-1 (H K^-1 y)
LB = L2.triangularView<Lower>().solve(d.adjoint() * a);
B = L2.transpose().triangularView<Upper>().solve(LB);
// Compute R
R = Hstar - d.adjoint() * b;
//
e = L2.triangularView<Lower>().solve(R);
// mean
M = b.adjoint() * a + R.adjoint() * B;
// Covariance
btb = MatrixXd(n,n).setZero().selfadjointView<Lower>().rankUpdate(b.adjoint());
ete = MatrixXd(n,n).setZero().selfadjointView<Lower>().rankUpdate(e.adjoint());
C = Kstarstar - btb + ete;
// log-likelihood
// second term in the log-likelihood
d2 = L2.triangularView<Lower>().solve(H);
w = d2 * L.adjoint().triangularView<Lower>().solve(a);
logLik1 = 0.5 * a.adjoint()*a;
logLik1 -= 0.5 * w.adjoint()*w;
logLik2 = L.diagonal();
logLik3 = L2.diagonal();
return Rcpp::List::create(Rcpp::Named("mean") = M,
Rcpp::Named("cov") = C,
Rcpp::Named("logLik1") = logLik1,
Rcpp::Named("logLik2") = logLik2,
Rcpp::Named("logLik3") = logLik3);
}