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