Automatic Differentiation
 
Loading...
Searching...
No Matches
gp_exp_quad_cov.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_GP_EXP_QUAD_COV_HPP
2#define STAN_MATH_PRIM_FUN_GP_EXP_QUAD_COV_HPP
3
11#include <cmath>
12#include <type_traits>
13#include <vector>
14
15namespace stan {
16namespace math {
17
18namespace internal {
19
33template <typename T_x, typename T_sigma, typename T_l>
34inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
35 Eigen::Dynamic>
36gp_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
80template <typename T_x1, typename T_x2, typename T_sigma, typename T_l>
81inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l>,
82 Eigen::Dynamic, Eigen::Dynamic>
83gp_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
123template <typename T_x, typename T_sigma, typename T_l>
124inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
125 Eigen::Dynamic>
126gp_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
163template <typename T_x, typename T_sigma, typename T_l>
164inline typename Eigen::Matrix<return_type_t<T_x, T_sigma, T_l>, Eigen::Dynamic,
165 Eigen::Dynamic>
166gp_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
206template <typename T_x1, typename T_x2, typename T_sigma, typename T_l>
207inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_sigma, T_l>,
208 Eigen::Dynamic, Eigen::Dynamic>
209gp_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
255template <typename T_x1, typename T_x2, typename T_s, typename T_l>
256inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>,
257 Eigen::Dynamic, Eigen::Dynamic>
258gp_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
void divide_columns(matrix_cl< T1 > &A, const matrix_cl< T2 > &B)
Divides each column of a matrix by a vector.
matrix_cl< return_type_t< T1, T2, T3 > > gp_exp_quad_cov(const matrix_cl< T1 > &x, const T2 sigma, const T3 length_scale)
Squared exponential kernel on the GPU.
int64_t size(const T &m)
Returns the size (number of the elements) of a matrix_cl or var_value<matrix_cl<T>>.
Definition size.hpp:19
Eigen::Matrix< return_type_t< T_x, T_sigma, T_l >, Eigen::Dynamic, 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)
Returns a squared exponential kernel.
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
auto squared_distance(const T_a &a, const T_b &b)
Returns the squared distance.
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
fvar< T > square(const fvar< T > &x)
Definition square.hpp:12
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...