1 #ifndef MATERNKERNEL_H
2 #define MATERNKERNEL_H
3 
4 #include "MUQ/Approximation/GaussianProcesses/KernelImpl.h"
5 
6 #include <cmath>
7 #include <stdexcept>
8 
9 #include <boost/property_tree/ptree_fwd.hpp>
10 
11 #include <boost/math/constants/constants.hpp>
12 
13 
14 
15 namespace muq
16 {
17 namespace Approximation
18 {
19 
20 
21 /**
22 
23 @class MaternKernel
24 @ingroup CovarianceKernels
25 This class implements a kernel of the form
26 \f[
27 k(x,y) = \sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}}{l}\tau\right)^\nu K_\nu\left(\frac{\sqrt{2\nu}}{l}\tau\right),
28 \f]
29 where \f$\Gamma(\cdot)\f$ is the gamma function, \f$K_\nu(\cdot)\f$ is a modified Bessel function of the second kind, \f$\nu\f$ is a smoothness parameter, \f$l\f$ is a lengthscale, and \f$\sigma^2\f$ is the variance.  Note that we only allow values of \f$\nu\f$ that take the form \f$i-0.5\f$ for some positive integer \f$i\f$.  Typical choices are \f$\nu=1/2\f$, \f$\nu=3/2\f$, and \f$\nu=5/2\f$.   In the limit, \f$\nu\rightarrow\infty\f$, this Matern kernel converges to the squared exponential kernel (muq::Approximation::SquaredExpKernel) and for \f$\nu=1/2\f$, the Matern kernel is equivalent to the exponential kernel associated with the Ornstein-Uhlenbeck process.
30 
31 Note: the smoothness parameter \f$\nu\f$ is not optimized as a hyperparameter.
32  */
33 class MaternKernel : public KernelImpl<MaternKernel>
34 {
35 
36 public:
37 
MaternKernel(unsigned dimIn,std::vector<unsigned> dimInds,double sigma2In,double lengthIn,double nuIn)38   MaternKernel(unsigned              dimIn,
39                std::vector<unsigned> dimInds,
40                double                sigma2In,
41                double                lengthIn,
42                double                nuIn) : MaternKernel(dimIn,
43                                                           dimInds,
44                                                           sigma2In,
45                                                           lengthIn,
46                                                           nuIn,
47                                                           {0.0, std::numeric_limits<double>::infinity()},
48                                                           {1e-10, std::numeric_limits<double>::infinity()})
49   {};
50 
51     MaternKernel(unsigned              dimIn,
52                  std::vector<unsigned> dimInds,
53                  double                sigma2In,
54                  double                lengthIn,
55                  double                nuIn,
56                  Eigen::Vector2d       sigmaBounds,
57                  Eigen::Vector2d       lengthBounds);
58 
59 
MaternKernel(unsigned dimIn,double sigma2In,double lengthIn,double nuIn)60     MaternKernel(unsigned        dimIn,
61                  double          sigma2In,
62                  double          lengthIn,
63                  double          nuIn) : MaternKernel(dimIn,
64                                                       sigma2In,
65                                                       lengthIn,
66                                                       nuIn,
67                                                       {0.0, std::numeric_limits<double>::infinity()},
68                                                       {1e-10, std::numeric_limits<double>::infinity()})
69     {};
70 
71     MaternKernel(unsigned        dimIn,
72                  double          sigma2In,
73                  double          lengthIn,
74                  double          nuIn,
75                  Eigen::Vector2d sigmaBounds,
76                  Eigen::Vector2d lengthBounds);
77 
78 
~MaternKernel()79     virtual ~MaternKernel(){};
80 
81 
82     template<typename ScalarType1, typename ScalarType2, typename ScalarType3>
FillBlockImpl(Eigen::Ref<const Eigen::Matrix<ScalarType1,Eigen::Dynamic,1>> const & x1,Eigen::Ref<const Eigen::Matrix<ScalarType1,Eigen::Dynamic,1>> const & x2,Eigen::Ref<const Eigen::Matrix<ScalarType2,Eigen::Dynamic,1>> const & params,Eigen::Ref<Eigen::Matrix<ScalarType3,Eigen::Dynamic,Eigen::Dynamic>> block)83     void FillBlockImpl(Eigen::Ref<const Eigen::Matrix<ScalarType1, Eigen::Dynamic, 1>> const& x1,
84                        Eigen::Ref<const Eigen::Matrix<ScalarType1, Eigen::Dynamic, 1>> const& x2,
85                        Eigen::Ref<const Eigen::Matrix<ScalarType2, Eigen::Dynamic, 1>> const& params,
86                        Eigen::Ref<Eigen::Matrix<ScalarType3,Eigen::Dynamic, Eigen::Dynamic>>  block) const
87     {
88       int p = round(nu-0.5);
89 
90       ScalarType1 dist = (x1-x2).norm();
91 
92       block(0,0) = 0.0;
93       for(int i=0; i<=p; ++i)
94         block(0,0) +=  weights(i) * pow(sqrt(8.0*nu)*dist/params(1), p-i);
95 
96       block(0,0) *= params(0)*exp(-sqrt(2.0*nu)*dist / params(1)) * scale;
97     }
98 
99     virtual std::tuple<std::shared_ptr<muq::Modeling::LinearSDE>, std::shared_ptr<muq::Modeling::LinearOperator>, Eigen::MatrixXd> GetStateSpace(boost::property_tree::ptree sdeOptions = boost::property_tree::ptree()) const override;
100 
101 private:
102 
103     const double nu;
104     const double scale; // std::pow(2.0, 1.0-nu)/boost::math::tgamma(nu)
105 
106     const Eigen::VectorXd weights;
107 
108     void CheckNu() const;
109 
110     static Eigen::VectorXd BuildWeights(int p);
111 
Factorial(int n)112     static inline int Factorial(int n){
113       return (n == 1 || n == 0) ? 1 : Factorial(n - 1) * n;
114     }
115 };
116 
117 }
118 }
119 
120 
121 #endif
122