Stan Math Library  2.15.0
reverse mode automatic differentiation
cov_exp_quad.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_COV_EXP_QUAD_HPP
2 #define STAN_MATH_REV_MAT_FUN_COV_EXP_QUAD_HPP
3 
4 #include <stan/math/rev/core.hpp>
13 #include <boost/math/tools/promotion.hpp>
14 #include <boost/utility/enable_if.hpp>
15 #include <boost/type_traits.hpp>
16 #include <vector>
17 #include <cmath>
18 
19 namespace stan {
20  namespace math {
21 
35  template <typename T_x, typename T_sigma, typename T_l>
36  class cov_exp_quad_vari : public vari {
37  public:
38  const size_t size_;
39  const size_t size_ltri_;
40  const double l_d_;
41  const double sigma_d_;
42  const double sigma_sq_d_;
43  double* dist_;
48 
67  cov_exp_quad_vari(const std::vector<T_x>& x,
68  const T_sigma& sigma,
69  const T_l& l)
70  : vari(0.0),
71  size_(x.size()),
72  size_ltri_(size_ * (size_ - 1) / 2),
73  l_d_(value_of(l)), sigma_d_(value_of(sigma)),
74  sigma_sq_d_(sigma_d_ * sigma_d_),
75  dist_(ChainableStack::memalloc_.alloc_array<double>(size_ltri_)),
76  l_vari_(l.vi_), sigma_vari_(sigma.vi_),
77  cov_lower_(ChainableStack::memalloc_.alloc_array<vari*>(size_ltri_)),
78  cov_diag_(ChainableStack::memalloc_.alloc_array<vari*>(size_)) {
79  double inv_half_sq_l_d = 0.5 / (l_d_ * l_d_);
80  size_t pos = 0;
81  for (size_t j = 0; j < size_ - 1; ++j) {
82  for (size_t i = j + 1; i < size_; ++i) {
83  double dist_sq = squared_distance(x[i], x[j]);
84  dist_[pos] = dist_sq;
85  cov_lower_[pos] = new vari(sigma_sq_d_
86  * std::exp(-dist_sq
87  * 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 cov_exp_quad_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  cov_exp_quad_vari(const std::vector<T_x>& x,
158  double sigma,
159  const T_l& l)
160  : vari(0.0),
161  size_(x.size()),
162  size_ltri_(size_ * (size_ - 1) / 2),
163  l_d_(value_of(l)), sigma_d_(value_of(sigma)),
164  sigma_sq_d_(sigma_d_ * sigma_d_),
165  dist_(ChainableStack::memalloc_.alloc_array<double>(size_ltri_)),
166  l_vari_(l.vi_),
167  cov_lower_(ChainableStack::memalloc_.alloc_array<vari*>(size_ltri_)),
168  cov_diag_(ChainableStack::memalloc_.alloc_array<vari*>(size_)) {
169  double inv_half_sq_l_d = 0.5 / (l_d_ * l_d_);
170  size_t pos = 0;
171  for (size_t j = 0; j < size_ - 1; ++j) {
172  for (size_t i = j + 1; i < size_; ++i) {
173  double dist_sq = squared_distance(x[i], x[j]);
174  dist_[pos] = dist_sq;
175  cov_lower_[pos] = new vari(sigma_sq_d_
176  * std::exp(-dist_sq
177  * inv_half_sq_l_d), false);
178  ++pos;
179  }
180  }
181  for (size_t i = 0; i < size_; ++i)
182  cov_diag_[i] = new vari(sigma_sq_d_, false);
183  }
184 
185  virtual void chain() {
186  double adjl = 0;
187 
188  for (size_t i = 0; i < size_ltri_; ++i) {
189  vari* el_low = cov_lower_[i];
190  adjl += el_low->adj_ * el_low->val_ * dist_[i];
191  }
192  l_vari_->adj_ += adjl / (l_d_ * l_d_ * l_d_);
193  }
194  };
195 
207  template <typename T_x>
208  inline typename
209  boost::enable_if_c<boost::is_same<typename scalar_type<T_x>::type,
210  double>::value,
211  Eigen::Matrix<var, -1, -1> >::type
212  cov_exp_quad(const std::vector<T_x>& x,
213  const var& sigma,
214  const var& l) {
215  check_positive("cov_exp_quad", "sigma", sigma);
216  check_positive("cov_exp_quad", "l", l);
217  size_t x_size = x.size();
218  for (size_t i = 0; i < x_size; ++i)
219  check_not_nan("cov_exp_quad", "x", x[i]);
220 
221  Eigen::Matrix<var, -1, -1>
222  cov(x_size, x_size);
223  if (x_size == 0)
224  return cov;
225 
227  = new cov_exp_quad_vari<T_x, var, var>(x, sigma, l);
228 
229  size_t pos = 0;
230  for (size_t j = 0; j < x_size - 1; ++j) {
231  for (size_t i = (j + 1); i < x_size; ++i) {
232  cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
233  cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
234  ++pos;
235  }
236  cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
237  }
238  cov.coeffRef(x_size - 1, x_size - 1).vi_
239  = baseVari->cov_diag_[x_size - 1];
240  return cov;
241  }
242 
254  template <typename T_x>
255  inline typename
256  boost::enable_if_c<boost::is_same<typename scalar_type<T_x>::type,
257  double>::value,
258  Eigen::Matrix<var, -1, -1> >::type
259  cov_exp_quad(const std::vector<T_x>& x,
260  double sigma,
261  const var& l) {
262  check_positive("cov_exp_quad", "marginal variance", sigma);
263  check_positive("cov_exp_quad", "length-scale", l);
264  size_t x_size = x.size();
265  for (size_t i = 0; i < x_size; ++i)
266  check_not_nan("cov_exp_quad", "x", x[i]);
267 
268  Eigen::Matrix<var, -1, -1>
269  cov(x_size, x_size);
270  if (x_size == 0)
271  return cov;
272 
274  = new cov_exp_quad_vari<T_x, double, var>(x, sigma, l);
275 
276  size_t pos = 0;
277  for (size_t j = 0; j < x_size - 1; ++j) {
278  for (size_t i = (j + 1); i < x_size; ++i) {
279  cov.coeffRef(i, j).vi_ = baseVari->cov_lower_[pos];
280  cov.coeffRef(j, i).vi_ = cov.coeffRef(i, j).vi_;
281  ++pos;
282  }
283  cov.coeffRef(j, j).vi_ = baseVari->cov_diag_[j];
284  }
285  cov.coeffRef(x_size - 1, x_size - 1).vi_
286  = baseVari->cov_diag_[x_size - 1];
287  return cov;
288  }
289 
290  }
291 }
292 #endif
This is a subclass of the vari class for precomputed gradients of cov_exp_quad.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:16
Eigen::Matrix< typename stan::return_type< T_x, T_sigma, T_l >::type, Eigen::Dynamic, Eigen::Dynamic > cov_exp_quad(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &l)
Returns a squared exponential kernel.
The variable implementation base class.
Definition: vari.hpp:30
friend class var
Definition: vari.hpp:32
const double val_
The value of this variable.
Definition: vari.hpp:38
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:10
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
cov_exp_quad_vari(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &l)
Constructor for cov_exp_quad.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
cov_exp_quad_vari(const std::vector< T_x > &x, double sigma, const T_l &l)
Constructor for cov_exp_quad.
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

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