Automatic Differentiation
 
Loading...
Searching...
No Matches
cov_exp_quad.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_REV_FUN_COV_EXP_QUAD_HPP
2#define STAN_MATH_REV_FUN_COV_EXP_QUAD_HPP
3
11#include <type_traits>
12#include <vector>
13#include <cmath>
14
15namespace stan {
16namespace math {
17
21template <typename T_x, typename T_sigma, typename T_l>
22class cov_exp_quad_vari : public vari {
23 public:
24 const size_t size_;
25 const size_t size_ltri_;
26 const double l_d_;
27 const double sigma_d_;
28 const double sigma_sq_d_;
29 double* dist_;
34
38 cov_exp_quad_vari(const std::vector<T_x>& x, const T_sigma& sigma,
39 const T_l& l)
40 : vari(0.0),
41 size_(x.size()),
42 size_ltri_(size_ * (size_ - 1) / 2),
43 l_d_(value_of(l)),
44 sigma_d_(value_of(sigma)),
46 dist_(ChainableStack::instance_->memalloc_.alloc_array<double>(
47 size_ltri_)),
48 l_vari_(l.vi_),
49 sigma_vari_(sigma.vi_),
50 cov_lower_(ChainableStack::instance_->memalloc_.alloc_array<vari*>(
51 size_ltri_)),
53 ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
54 double inv_half_sq_l_d = 0.5 / (l_d_ * l_d_);
55 size_t pos = 0;
56 for (size_t j = 0; j < size_ - 1; ++j) {
57 for (size_t i = j + 1; i < size_; ++i) {
58 double dist_sq = squared_distance(x[i], x[j]);
59 dist_[pos] = dist_sq;
60 cov_lower_[pos] = new vari(
61 sigma_sq_d_ * std::exp(-dist_sq * inv_half_sq_l_d), false);
62 ++pos;
63 }
64 }
65 for (size_t i = 0; i < size_; ++i) {
66 cov_diag_[i] = new vari(sigma_sq_d_, false);
67 }
68 }
69
70 virtual void chain() {
71 double adjl = 0;
72 double adjsigma = 0;
73
74 for (size_t i = 0; i < size_ltri_; ++i) {
75 vari* el_low = cov_lower_[i];
76 double prod_add = el_low->adj_ * el_low->val_;
77 adjl += prod_add * dist_[i];
78 adjsigma += prod_add;
79 }
80 for (size_t i = 0; i < size_; ++i) {
81 vari* el = cov_diag_[i];
82 adjsigma += el->adj_ * el->val_;
83 }
84 l_vari_->adj_ += adjl / (l_d_ * l_d_ * l_d_);
85 sigma_vari_->adj_ += adjsigma * 2 / sigma_d_;
86 }
87};
88
92template <typename T_x, typename T_l>
93class cov_exp_quad_vari<T_x, double, T_l> : public vari {
94 public:
95 const size_t size_;
96 const size_t size_ltri_;
97 const double l_d_;
98 const double sigma_d_;
99 const double sigma_sq_d_;
100 double* dist_;
104
108 cov_exp_quad_vari(const std::vector<T_x>& x, double sigma, const T_l& l)
109 : vari(0.0),
110 size_(x.size()),
111 size_ltri_(size_ * (size_ - 1) / 2),
112 l_d_(value_of(l)),
113 sigma_d_(value_of(sigma)),
115 dist_(ChainableStack::instance_->memalloc_.alloc_array<double>(
116 size_ltri_)),
117 l_vari_(l.vi_),
118 cov_lower_(ChainableStack::instance_->memalloc_.alloc_array<vari*>(
119 size_ltri_)),
120 cov_diag_(
121 ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
122 double inv_half_sq_l_d = 0.5 / (l_d_ * l_d_);
123 size_t pos = 0;
124 for (size_t j = 0; j < size_ - 1; ++j) {
125 for (size_t i = j + 1; i < size_; ++i) {
126 double dist_sq = squared_distance(x[i], x[j]);
127 dist_[pos] = dist_sq;
128 cov_lower_[pos] = new vari(
129 sigma_sq_d_ * std::exp(-dist_sq * inv_half_sq_l_d), false);
130 ++pos;
131 }
132 }
133 for (size_t i = 0; i < size_; ++i) {
134 cov_diag_[i] = new vari(sigma_sq_d_, false);
135 }
136 }
137
138 virtual void chain() {
139 double adjl = 0;
140
141 for (size_t i = 0; i < size_ltri_; ++i) {
142 vari* el_low = cov_lower_[i];
143 adjl += el_low->adj_ * el_low->val_ * dist_[i];
144 }
145 l_vari_->adj_ += adjl / (l_d_ * l_d_ * l_d_);
146 }
147};
148
152template <typename T_x,
154inline Eigen::Matrix<var, -1, -1> cov_exp_quad(const std::vector<T_x>& x,
155 const var& sigma, const var& l) {
156 return gp_exp_quad_cov(x, sigma, l);
157}
158
162template <typename T_x,
164inline Eigen::Matrix<var, -1, -1> cov_exp_quad(const std::vector<T_x>& x,
165 double sigma, const var& l) {
166 return gp_exp_quad_cov(x, sigma, l);
167}
168
169} // namespace math
170} // namespace stan
171#endif
cov_exp_quad_vari(const std::vector< T_x > &x, double sigma, const T_l &l)
cov_exp_quad_vari(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &l)
require_t< std::is_arithmetic< std::decay_t< T > > > require_arithmetic_t
Require type satisfies std::is_arithmetic.
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
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
vari_value< double > vari
Definition vari.hpp:197
Eigen::Matrix< return_type_t< T_x, T_sigma, T_l >, Eigen::Dynamic, Eigen::Dynamic > cov_exp_quad(const std::vector< T_x > &x, const T_sigma &sigma, const T_l &length_scale)
var_value< double > var
Definition var.hpp:1187
auto squared_distance(const T_a &a, const T_b &b)
Returns the squared distance.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
This struct always provides access to the autodiff stack using the singleton pattern.