1 #include "MUQ/Approximation/GaussianProcesses/MaternKernel.h"
2
3 #include "MUQ/Modeling/LinearAlgebra/LinearOperator.h"
4 #include "MUQ/Modeling/LinearAlgebra/EigenLinearOperator.h"
5
6 #include "MUQ/Approximation/GaussianProcesses/StateSpaceGP.h"
7
8
9 #include "MUQ/Modeling/LinearAlgebra/CompanionMatrix.h"
10 #include "MUQ/Modeling/LinearAlgebra/LyapunovSolver.h"
11
12 #include "MUQ/Modeling/LinearSDE.h"
13
14
15 #include <unsupported/Eigen/Polynomials>
16 #include <Eigen/SparseCore>
17
18 #include <boost/property_tree/ptree.hpp>
19
20 using namespace muq::Approximation;
21
22
23
MaternKernel(unsigned dimIn,std::vector<unsigned> dimInds,double sigma2In,double lengthIn,double nuIn,Eigen::Vector2d sigmaBounds,Eigen::Vector2d lengthBounds)24 MaternKernel::MaternKernel(unsigned dimIn,
25 std::vector<unsigned> dimInds,
26 double sigma2In,
27 double lengthIn,
28 double nuIn,
29 Eigen::Vector2d sigmaBounds,
30 Eigen::Vector2d lengthBounds) : KernelImpl<MaternKernel>(dimIn, dimInds, 1, 2),
31 nu(nuIn),
32 scale(boost::math::tgamma(nuIn+0.5)/boost::math::tgamma(2.0*nuIn)),
33 weights(BuildWeights(nu-0.5))
34 {
35 CheckNu();
36
37 paramBounds.resize(2,2);
38 paramBounds(0,0) = sigmaBounds(0);
39 paramBounds(1,0) = sigmaBounds(1);
40
41 paramBounds(0,1) = lengthBounds(0);
42 paramBounds(1,1) = lengthBounds(1);
43
44 cachedParams.resize(2);
45 cachedParams(0) = sigma2In;
46 cachedParams(1) = lengthIn;
47 };
48
MaternKernel(unsigned dimIn,double sigma2In,double lengthIn,double nuIn,Eigen::Vector2d sigmaBounds,Eigen::Vector2d lengthBounds)49 MaternKernel::MaternKernel(unsigned dimIn,
50 double sigma2In,
51 double lengthIn,
52 double nuIn,
53 Eigen::Vector2d sigmaBounds,
54 Eigen::Vector2d lengthBounds) : KernelImpl<MaternKernel>(dimIn, 1, 2),
55 nu(nuIn),
56 scale(boost::math::tgamma(nuIn+0.5)/boost::math::tgamma(2.0*nuIn)),
57 weights(BuildWeights(nu-0.5))
58 {
59 CheckNu();
60
61 paramBounds.resize(2,2);
62 paramBounds(0,0) = sigmaBounds(0);
63 paramBounds(1,0) = sigmaBounds(1);
64
65 paramBounds(0,1) = lengthBounds(0);
66 paramBounds(1,1) = lengthBounds(1);
67
68 cachedParams.resize(2);
69 cachedParams(0) = sigma2In;
70 cachedParams(1) = lengthIn;
71 };
72
73
CheckNu() const74 void MaternKernel::CheckNu() const{
75
76 if(nu<=0)
77 throw std::invalid_argument("The value of nu must be greater than 0.");
78
79 if((std::round(nu-0.5)-(nu-0.5)) > 4.0*std::numeric_limits<double>::epsilon())
80 throw std::invalid_argument("The value of nu must take the form nu=i-0.5 for a positive integer i.");
81
82 };
83
84
GetStateSpace(boost::property_tree::ptree sdeOptions) const85 std::tuple<std::shared_ptr<muq::Modeling::LinearSDE>, std::shared_ptr<muq::Modeling::LinearOperator>, Eigen::MatrixXd> MaternKernel::GetStateSpace(boost::property_tree::ptree sdeOptions) const
86 {
87 double sigma2 = cachedParams(0);
88 double length = cachedParams(1);
89 int p = nu-0.5;
90
91 double lambda = sqrt(2.0*nu)/length;
92
93 double q = 2.0*sigma2*boost::math::constants::root_pi<double>() * std::pow(lambda, 2*p+1) * tgamma(p+1) / tgamma(p+0.5);
94
95 Eigen::VectorXd roots = -lambda*Eigen::VectorXd::Ones(p+1);
96 Eigen::VectorXd poly;
97
98 Eigen::roots_to_monicPolynomial( roots, poly );
99
100 poly = poly.head(poly.size()-1).eval();
101
102 auto F = std::make_shared<muq::Modeling::CompanionMatrix>(-1.0*poly);
103
104 std::vector<Eigen::Triplet<double>> Lcoeffs;
105 Lcoeffs.push_back(Eigen::Triplet<double>(poly.size()-1,0,1.0));
106
107 Eigen::SparseMatrix<double> Lmat(poly.size(), 1);
108 Lmat.setFromTriplets(Lcoeffs.begin(), Lcoeffs.end());
109
110 auto L = muq::Modeling::LinearOperator::Create(Lmat);
111
112 Eigen::MatrixXd Q(1,1);
113 Q(0,0) = q;
114
115
116 // Set up the stochastic differential equation
117 auto sde = std::make_shared<muq::Modeling::LinearSDE>(F, L, Q, sdeOptions);
118
119
120 // Define the observation operator, which is just (1,0,...,0) in this case
121 std::vector<Eigen::Triplet<double>> Hcoeffs;
122 Hcoeffs.push_back(Eigen::Triplet<double>(0,0,1.0));
123
124 Eigen::SparseMatrix<double> Hmat(1,poly.size());
125 Hmat.setFromTriplets(Hcoeffs.begin(), Hcoeffs.end());
126
127 auto H = muq::Modeling::LinearOperator::Create(Hmat);
128
129 // Solve the continuous time Lyapunov equation to find the stationary covariance
130 Q = L->Apply(L->Apply(q*Eigen::VectorXd::Ones(1)).transpose());
131
132 Eigen::MatrixXd Pinf = muq::Modeling::LyapunovSolver<double>().compute(F->GetMatrix().transpose(), Q).matrixX().real();
133
134 return std::make_tuple(sde, H, Pinf);
135 }
136
137
BuildWeights(int p)138 Eigen::VectorXd MaternKernel::BuildWeights(int p)
139 {
140 Eigen::VectorXd output(p+1);
141 for(int i=0; i<=p; ++i){
142
143 if(i<p){
144 output(i) = static_cast<double>(Factorial(p+i))/(static_cast<double>(Factorial(i)*Factorial(p-i)));
145 }else{
146 output(i) = static_cast<double>(Factorial(p+i)) / static_cast<double>(Factorial(i));
147 }
148 }
149 return output;
150 }
151