1 #ifndef STAN_MATH_PRIM_FUN_GP_EXP_QUAD_COV_HPP
2 #define STAN_MATH_PRIM_FUN_GP_EXP_QUAD_COV_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/Eigen.hpp>
7 #include <stan/math/prim/fun/divide_columns.hpp>
8 #include <stan/math/prim/fun/squared_distance.hpp>
9 #include <stan/math/prim/fun/exp.hpp>
10 #include <stan/math/prim/fun/square.hpp>
11 #include <cmath>
12 #include <type_traits>
13 #include <vector>
14
15 namespace stan {
16 namespace math {
17
18 namespace internal {
19
20 /**
21 * Returns a squared exponential kernel.
22 *
23 * @tparam T_x type for each scalar
24 * @tparam T_sigma type of parameter sigma
25 * @tparam T_l type of parameter length scale
26 *
27 * @param x std::vector of scalars that can be used in square distance.
28 * This function assumes each element of x is the same size.
29 * @param sigma_sq square root of the marginal standard deviation or magnitude
30 * @param neg_half_inv_l_sq The half negative inverse of the length scale
31 * @return squared distance
32 */
33 template <typename T_x, typename T_sigma, typename T_l>
34 inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
35 Eigen::Dynamic>
gp_exp_quad_cov(const std::vector<T_x> & x,const T_sigma & sigma_sq,const T_l & neg_half_inv_l_sq)36 gp_exp_quad_cov(const std::vector<T_x> &x, const T_sigma &sigma_sq,
37 const T_l &neg_half_inv_l_sq) {
38 using std::exp;
39 const size_t x_size = x.size();
40 Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
41 Eigen::Dynamic>
42 cov(x_size, x_size);
43 cov.diagonal().array() = sigma_sq;
44 size_t block_size = 10;
45 for (size_t jb = 0; jb < x.size(); jb += block_size) {
46 for (size_t ib = jb; ib < x.size(); ib += block_size) {
47 size_t j_end = std::min(x_size, jb + block_size);
48 for (size_t j = jb; j < j_end; ++j) {
49 size_t i_end = std::min(x_size, ib + block_size);
50 for (size_t i = std::max(ib, j + 1); i < i_end; ++i) {
51 cov.coeffRef(i, j)
52 = sigma_sq
53 * exp(squared_distance(x[i], x[j]) * neg_half_inv_l_sq);
54 }
55 }
56 }
57 }
58 cov.template triangularView<Eigen::Upper>() = cov.transpose();
59 return cov;
60 }
61
62 /**
63 * Returns a squared exponential kernel.
64 *
65 * This function is for the cross covariance matrix
66 * needed to compute posterior predictive density.
67 *
68 * @tparam T_x1 type of first std::vector of scalars
69 * @tparam T_x2 type of second std::vector of scalars
70 * This function assumes each element of x1 and x2 are the same size.
71 * @tparam T_sigma type of sigma
72 * @tparam T_l type of of length scale
73 *
74 * @param x1 std::vector of elements that can be used in square distance
75 * @param x2 std::vector of elements that can be used in square distance
76 * @param sigma_sq square root of the marginal standard deviation or magnitude
77 * @param neg_half_inv_l_sq The half negative inverse of the length scale
78 * @return squared distance
79 */
80 template <typename T_x1, typename T_x2, typename T_sigma, typename T_l>
81 inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l>,
82 Eigen::Dynamic, Eigen::Dynamic>
gp_exp_quad_cov(const std::vector<T_x1> & x1,const std::vector<T_x2> & x2,const T_sigma & sigma_sq,const T_l & neg_half_inv_l_sq)83 gp_exp_quad_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
84 const T_sigma &sigma_sq, const T_l &neg_half_inv_l_sq) {
85 using std::exp;
86 Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l>, Eigen::Dynamic,
87 Eigen::Dynamic>
88 cov(x1.size(), x2.size());
89 size_t block_size = 10;
90
91 for (size_t ib = 0; ib < x1.size(); ib += block_size) {
92 for (size_t jb = 0; jb < x2.size(); jb += block_size) {
93 size_t j_end = std::min(x2.size(), jb + block_size);
94 for (size_t j = jb; j < j_end; ++j) {
95 size_t i_end = std::min(x1.size(), ib + block_size);
96 for (size_t i = ib; i < i_end; ++i) {
97 cov.coeffRef(i, j)
98 = sigma_sq
99 * exp(squared_distance(x1[i], x2[j]) * neg_half_inv_l_sq);
100 }
101 }
102 }
103 }
104 return cov;
105 }
106 } // namespace internal
107
108 /**
109 * Returns a squared exponential kernel.
110 *
111 * @tparam T_x type for each scalar
112 * @tparam T_sigma type of parameter sigma
113 * @tparam T_l type of parameter length scale
114 *
115 * @param x std::vector of scalars that can be used in square distance.
116 * This function assumes each element of x is the same size.
117 * @param sigma marginal standard deviation or magnitude
118 * @param length_scale length scale
119 * @return squared distance
120 * @throw std::domain_error if sigma <= 0, l <= 0, or
121 * x is nan or infinite
122 */
123 template <typename T_x, typename T_sigma, typename T_l>
124 inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
125 Eigen::Dynamic>
gp_exp_quad_cov(const std::vector<T_x> & x,const T_sigma & sigma,const T_l & length_scale)126 gp_exp_quad_cov(const std::vector<T_x> &x, const T_sigma &sigma,
127 const T_l &length_scale) {
128 check_positive("gp_exp_quad_cov", "magnitude", sigma);
129 check_positive("gp_exp_quad_cov", "length scale", length_scale);
130
131 const size_t x_size = x.size();
132 Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
133 Eigen::Dynamic>
134 cov(x_size, x_size);
135
136 if (x_size == 0) {
137 return cov;
138 }
139
140 for (size_t n = 0; n < x_size; ++n) {
141 check_not_nan("gp_exp_quad_cov", "x", x[n]);
142 }
143
144 cov = internal::gp_exp_quad_cov(x, square(sigma),
145 -0.5 / square(length_scale));
146 return cov;
147 }
148
149 /**
150 * Returns a squared exponential kernel.
151 *
152 * @tparam T_x type for each scalar
153 * @tparam T_sigma type of parameter sigma
154 * @tparam T_l type of each length scale parameter
155 *
156 * @param x std::vector of Eigen vectors of scalars.
157 * @param sigma marginal standard deviation or magnitude
158 * @param length_scale std::vector length scale
159 * @return squared distance
160 * @throw std::domain_error if sigma <= 0, l <= 0, or
161 * x is nan or infinite
162 */
163 template <typename T_x, typename T_sigma, typename T_l>
164 inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
165 Eigen::Dynamic>
gp_exp_quad_cov(const std::vector<Eigen::Matrix<T_x,-1,1>> & x,const T_sigma & sigma,const std::vector<T_l> & length_scale)166 gp_exp_quad_cov(const std::vector<Eigen::Matrix<T_x, -1, 1>> &x,
167 const T_sigma &sigma, const std::vector<T_l> &length_scale) {
168 check_positive_finite("gp_exp_quad_cov", "magnitude", sigma);
169 check_positive_finite("gp_exp_quad_cov", "length scale", length_scale);
170
171 size_t x_size = x.size();
172 Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
173 Eigen::Dynamic>
174 cov(x_size, x_size);
175 if (x_size == 0) {
176 return cov;
177 }
178
179 check_size_match("gp_exp_quad_cov", "x dimension", x[0].size(),
180 "number of length scales", length_scale.size());
181 cov = internal::gp_exp_quad_cov(divide_columns(x, length_scale),
182 square(sigma), -0.5);
183 return cov;
184 }
185
186 /**
187 * Returns a squared exponential kernel.
188 *
189 * This function is for the cross covariance matrix
190 * needed to compute posterior predictive density.
191 *
192 * @tparam T_x1 type of first std::vector of scalars
193 * @tparam T_x2 type of second std::vector of scalars
194 * This function assumes each element of x1 and x2 are the same size.
195 * @tparam T_sigma type of sigma
196 * @tparam T_l type of of length scale
197 *
198 * @param x1 std::vector of elements that can be used in square distance
199 * @param x2 std::vector of elements that can be used in square distance
200 * @param sigma standard deviation
201 * @param length_scale length scale
202 * @return squared distance
203 * @throw std::domain_error if sigma <= 0, l <= 0, or
204 * x is nan or infinite
205 */
206 template <typename T_x1, typename T_x2, typename T_sigma, typename T_l>
207 inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l>,
208 Eigen::Dynamic, Eigen::Dynamic>
gp_exp_quad_cov(const std::vector<T_x1> & x1,const std::vector<T_x2> & x2,const T_sigma & sigma,const T_l & length_scale)209 gp_exp_quad_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
210 const T_sigma &sigma, const T_l &length_scale) {
211 const char *function_name = "gp_exp_quad_cov";
212 check_positive(function_name, "magnitude", sigma);
213 check_positive(function_name, "length scale", length_scale);
214
215 const size_t x1_size = x1.size();
216 const size_t x2_size = x2.size();
217 Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l>, Eigen::Dynamic,
218 Eigen::Dynamic>
219 cov(x1_size, x2_size);
220 if (x1_size == 0 || x2_size == 0) {
221 return cov;
222 }
223
224 for (size_t i = 0; i < x1_size; ++i) {
225 check_not_nan(function_name, "x1", x1[i]);
226 }
227 for (size_t i = 0; i < x2_size; ++i) {
228 check_not_nan(function_name, "x2", x2[i]);
229 }
230
231 cov = internal::gp_exp_quad_cov(x1, x2, square(sigma),
232 -0.5 / square(length_scale));
233 return cov;
234 }
235
236 /**
237 * Returns a squared exponential kernel.
238 *
239 * This function is for the cross covariance
240 * matrix needed to compute the posterior predictive density.
241 *
242 * @tparam T_x1 type of first std::vector of elements
243 * @tparam T_x2 type of second std::vector of elements
244 * @tparam T_s type of sigma
245 * @tparam T_l type of length scale
246 *
247 * @param x1 std::vector of Eigen vectors of scalars.
248 * @param x2 std::vector of Eigen vectors of scalars.
249 * @param sigma standard deviation
250 * @param length_scale std::vector of length scale
251 * @return squared distance
252 * @throw std::domain_error if sigma <= 0, l <= 0, or
253 * x is nan or infinite
254 */
255 template <typename T_x1, typename T_x2, typename T_s, typename T_l>
256 inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>,
257 Eigen::Dynamic, Eigen::Dynamic>
gp_exp_quad_cov(const std::vector<Eigen::Matrix<T_x1,-1,1>> & x1,const std::vector<Eigen::Matrix<T_x2,-1,1>> & x2,const T_s & sigma,const std::vector<T_l> & length_scale)258 gp_exp_quad_cov(const std::vector<Eigen::Matrix<T_x1, -1, 1>> &x1,
259 const std::vector<Eigen::Matrix<T_x2, -1, 1>> &x2,
260 const T_s &sigma, const std::vector<T_l> &length_scale) {
261 size_t x1_size = x1.size();
262 size_t x2_size = x2.size();
263 size_t l_size = length_scale.size();
264
265 Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>, Eigen::Dynamic,
266 Eigen::Dynamic>
267 cov(x1_size, x2_size);
268
269 if (x1_size == 0 || x2_size == 0) {
270 return cov;
271 }
272
273 const char *function_name = "gp_exp_quad_cov";
274 for (size_t i = 0; i < x1_size; ++i) {
275 check_not_nan(function_name, "x1", x1[i]);
276 }
277 for (size_t i = 0; i < x2_size; ++i) {
278 check_not_nan(function_name, "x2", x2[i]);
279 }
280 check_positive_finite(function_name, "magnitude", sigma);
281 check_positive_finite(function_name, "length scale", length_scale);
282 check_size_match(function_name, "x dimension", x1[0].size(),
283 "number of length scales", l_size);
284 check_size_match(function_name, "x dimension", x2[0].size(),
285 "number of length scales", l_size);
286 cov = internal::gp_exp_quad_cov(divide_columns(x1, length_scale),
287 divide_columns(x2, length_scale),
288 square(sigma), -0.5);
289 return cov;
290 }
291
292 } // namespace math
293 } // namespace stan
294 #endif
295