Stan Math Library  2.20.0
reverse mode automatic differentiation
gp_exp_quad_cov.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_GP_EXP_QUAD_COV_HPP
2 #define STAN_MATH_REV_MAT_FUN_GP_EXP_QUAD_COV_HPP
3 
4 #include <stan/math/rev/meta.hpp>
10 #include <stan/math/rev/core.hpp>
12 #include <vector>
13 #include <cmath>
14 #include <type_traits>
15 
16 namespace stan {
17 namespace math {
18 
32 template <typename T_x, typename T_sigma, typename T_l>
33 class gp_exp_quad_cov_vari : public vari {
34  public:
35  const size_t size_;
36  const size_t size_ltri_;
37  const double l_d_;
38  const double sigma_d_;
39  const double sigma_sq_d_;
40  double *dist_;
45 
64  gp_exp_quad_cov_vari(const std::vector<T_x> &x, const T_sigma &sigma,
65  const T_l &length_scale)
66  : vari(0.0),
67  size_(x.size()),
68  size_ltri_(size_ * (size_ - 1) / 2),
69  l_d_(value_of(length_scale)),
70  sigma_d_(value_of(sigma)),
71  sigma_sq_d_(sigma_d_ * sigma_d_),
72  dist_(ChainableStack::instance_->memalloc_.alloc_array<double>(
73  size_ltri_)),
74  l_vari_(length_scale.vi_),
75  sigma_vari_(sigma.vi_),
76  cov_lower_(ChainableStack::instance_->memalloc_.alloc_array<vari *>(
77  size_ltri_)),
78  cov_diag_(
79  ChainableStack::instance_->memalloc_.alloc_array<vari *>(size_)) {
80  double inv_half_sq_l_d = 0.5 / (l_d_ * l_d_);
81  size_t pos = 0;
82  for (size_t j = 0; j < size_ - 1; ++j) {
83  for (size_t i = j + 1; i < size_; ++i) {
84  double dist_sq = squared_distance(x[i], x[j]);
85  dist_[pos] = dist_sq;
86  cov_lower_[pos] = new vari(
87  sigma_sq_d_ * std::exp(-dist_sq * inv_half_sq_l_d), false);
88  ++pos;
89  }
90  }
91  for (size_t i = 0; i < size_; ++i)
92  cov_diag_[i] = new vari(sigma_sq_d_, false);
93  }
94 
95  virtual void chain() {
96  double adjl = 0;
97  double adjsigma = 0;
98 
99  for (size_t i = 0; i < size_ltri_; ++i) {
100  vari *el_low = cov_lower_[i];
101  double prod_add = el_low->adj_ * el_low->val_;
102  adjl += prod_add * dist_[i];
103  adjsigma += prod_add;
104  }
105  for (size_t i = 0; i < size_; ++i) {
106  vari *el = cov_diag_[i];
107  adjsigma += el->adj_ * el->val_;
108  }
109  l_vari_->adj_ += adjl / (l_d_ * l_d_ * l_d_);
110  sigma_vari_->adj_ += adjsigma * 2 / sigma_d_;
111  }
112 };
113 
126 template <typename T_x, typename T_l>
127 class gp_exp_quad_cov_vari<T_x, double, T_l> : public vari {
128  public:
129  const size_t size_;
130  const size_t size_ltri_;
131  const double l_d_;
132  const double sigma_d_;
133  const double sigma_sq_d_;
134  double *dist_;
138 
157  gp_exp_quad_cov_vari(const std::vector<T_x> &x, double sigma,
158  const T_l &length_scale)
159  : vari(0.0),
160  size_(x.size()),
161  size_ltri_(size_ * (size_ - 1) / 2),
162  l_d_(value_of(length_scale)),
163  sigma_d_(value_of(sigma)),
164  sigma_sq_d_(sigma_d_ * sigma_d_),
165  dist_(ChainableStack::instance_->memalloc_.alloc_array<double>(
166  size_ltri_)),
167  l_vari_(length_scale.vi_),
168  cov_lower_(ChainableStack::instance_->memalloc_.alloc_array<vari *>(
169  size_ltri_)),
170  cov_diag_(
171  ChainableStack::instance_->memalloc_.alloc_array<vari *>(size_)) {
172  double inv_half_sq_l_d = 0.5 / (l_d_ * l_d_);
173  size_t pos = 0;
174  for (size_t j = 0; j < size_ - 1; ++j) {
175  for (size_t i = j + 1; i < size_; ++i) {
176  double dist_sq = squared_distance(x[i], x[j]);
177  dist_[pos] = dist_sq;
178  cov_lower_[pos] = new vari(
179  sigma_sq_d_ * std::exp(-dist_sq * inv_half_sq_l_d), false);
180  ++pos;
181  }
182  }
183  for (size_t i = 0; i < size_; ++i)
184  cov_diag_[i] = new vari(sigma_sq_d_, false);
185  }
186 
187  virtual void chain() {
188  double adjl = 0;
189 
190  for (size_t i = 0; i < size_ltri_; ++i) {
191  vari *el_low = cov_lower_[i];
192  adjl += el_low->adj_ * el_low->val_ * dist_[i];
193  }
194  l_vari_->adj_ += adjl / (l_d_ * l_d_ * l_d_);
195  }
196 };
197 
209 template <typename T_x>
210 inline typename std::enable_if<
211  std::is_same<typename scalar_type<T_x>::type, double>::value,
212  Eigen::Matrix<var, -1, -1>>::type
213 gp_exp_quad_cov(const std::vector<T_x> &x, const var &sigma,
214  const var &length_scale) {
215  check_positive("gp_exp_quad_cov", "sigma", sigma);
216  check_positive("gp_exp_quad_cov", "length_scale", length_scale);
217  size_t x_size = x.size();
218  for (size_t i = 0; i < x_size; ++i)
219  check_not_nan("gp_exp_quad_cov", "x", x[i]);
220 
221  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
222  if (x_size == 0)
223  return cov;
224 
226  = new gp_exp_quad_cov_vari<T_x, var, var>(x, sigma, length_scale);
227 
228  size_t pos = 0;
229  for (size_t j = 0; j < x_size - 1; ++j) {
230  for (size_t i = (j + 1); i < x_size; ++i) {
231  cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
232  cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
233  ++pos;
234  }
235  cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
236  }
237  cov.coeffRef(x_size - 1, x_size - 1).vi_ = baseVari->cov_diag_[x_size - 1];
238  return cov;
239 }
240 
252 template <typename T_x>
253 inline typename std::enable_if<
254  std::is_same<typename scalar_type<T_x>::type, double>::value,
255  Eigen::Matrix<var, -1, -1>>::type
256 gp_exp_quad_cov(const std::vector<T_x> &x, double sigma,
257  const var &length_scale) {
258  check_positive("gp_exp_quad_cov", "marginal variance", sigma);
259  check_positive("gp_exp_quad_cov", "length-scale", length_scale);
260  size_t x_size = x.size();
261  for (size_t i = 0; i < x_size; ++i)
262  check_not_nan("gp_exp_quad_cov", "x", x[i]);
263 
264  Eigen::Matrix<var, -1, -1> cov(x_size, x_size);
265  if (x_size == 0)
266  return cov;
267 
269  = new gp_exp_quad_cov_vari<T_x, double, var>(x, sigma, length_scale);
270 
271  size_t pos = 0;
272  for (size_t j = 0; j < x_size - 1; ++j) {
273  for (size_t i = (j + 1); i < x_size; ++i) {
274  cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
275  cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
276  ++pos;
277  }
278  cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
279  }
280  cov.coeffRef(x_size - 1, x_size - 1).vi_ = baseVari->cov_diag_[x_size - 1];
281  return cov;
282 }
283 
284 } // namespace math
285 } // namespace stan
286 #endif
gp_exp_quad_cov_vari(const std::vector< T_x > &x, double sigma, const T_l &length_scale)
Constructor for gp_exp_quad_cov.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
Eigen::Matrix< typename stan::return_type< T_x, T_sigma, T_l >::type, Eigen::Dynamic, Eigen::Dynamic > gp_exp_quad_cov(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &length_scale)
Returns a squared exponential kernel.
The variable implementation base class.
Definition: vari.hpp:30
gp_exp_quad_cov_vari(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &length_scale)
Constructor for gp_exp_quad_cov.
friend class var
Definition: vari.hpp:32
const double val_
The value of this variable.
Definition: vari.hpp:38
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
This is a subclass of the vari class for precomputed gradients of gp_exp_quad_cov.
fvar< T > squared_distance(const Eigen::Matrix< fvar< T >, R, C > &v1, const Eigen::Matrix< double, R, C > &v2)
Returns the squared distance between the specified vectors of the same dimensions.
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:11
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
vari(double x)
Construct a variable implementation from a value.
Definition: vari.hpp:58
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
This struct always provides access to the autodiff stack using the singleton pattern.

     [ Stan Home Page ] © 2011–2018, Stan Development Team.