Automatic Differentiation
 
Loading...
Searching...
No Matches
trace_gen_quad_form.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_REV_FUN_TRACE_GEN_QUAD_FORM_HPP
2#define STAN_MATH_REV_FUN_TRACE_GEN_QUAD_FORM_HPP
3
12#include <type_traits>
13
14namespace stan {
15namespace math {
16namespace internal {
17
18template <typename Td, int Rd, int Cd, typename Ta, int Ra, int Ca, typename Tb,
19 int Rb, int Cb>
21 public:
22 trace_gen_quad_form_vari_alloc(const Eigen::Matrix<Td, Rd, Cd>& D,
23 const Eigen::Matrix<Ta, Ra, Ca>& A,
24 const Eigen::Matrix<Tb, Rb, Cb>& B)
25 : D_(D), A_(A), B_(B) {}
26
27 double compute() {
29 }
30
31 Eigen::Matrix<Td, Rd, Cd> D_;
32 Eigen::Matrix<Ta, Ra, Ca> A_;
33 Eigen::Matrix<Tb, Rb, Cb> B_;
34};
35
36template <typename Td, int Rd, int Cd, typename Ta, int Ra, int Ca, typename Tb,
37 int Rb, int Cb>
39 protected:
40 static inline void computeAdjoints(double adj,
41 const Eigen::Matrix<double, Rd, Cd>& D,
42 const Eigen::Matrix<double, Ra, Ca>& A,
43 const Eigen::Matrix<double, Rb, Cb>& B,
44 Eigen::Matrix<var, Rd, Cd>* varD,
45 Eigen::Matrix<var, Ra, Ca>* varA,
46 Eigen::Matrix<var, Rb, Cb>* varB) {
47 Eigen::Matrix<double, Ca, Cb> AtB;
48 Eigen::Matrix<double, Ra, Cb> BD;
49 if (varB || varA) {
50 BD.noalias() = B * D;
51 }
52 if (varB || varD) {
53 AtB.noalias() = A.transpose() * B;
54 }
55
56 if (varB) {
57 (*varB).adj() += adj * (A * BD + AtB * D.transpose());
58 }
59 if (varA) {
60 (*varA).adj() += adj * (B * BD.transpose());
61 }
62 if (varD) {
63 (*varD).adj() += adj * (B.transpose() * AtB);
64 }
65 }
66
67 public:
70 : vari(impl->compute()), impl_(impl) {}
71
72 virtual void chain() {
74 value_of(impl_->B_),
75 reinterpret_cast<Eigen::Matrix<var, Rd, Cd>*>(
76 std::is_same<Td, var>::value ? (&impl_->D_) : NULL),
77 reinterpret_cast<Eigen::Matrix<var, Ra, Ca>*>(
78 std::is_same<Ta, var>::value ? (&impl_->A_) : NULL),
79 reinterpret_cast<Eigen::Matrix<var, Rb, Cb>*>(
80 std::is_same<Tb, var>::value ? (&impl_->B_) : NULL));
81 }
82
84};
85} // namespace internal
86
87template <typename Td, typename Ta, typename Tb,
91inline var trace_gen_quad_form(const Td& D, const Ta& A, const Tb& B) {
92 using Td_scal = value_type_t<Td>;
93 using Ta_scal = value_type_t<Ta>;
94 using Tb_scal = value_type_t<Tb>;
95 constexpr int Rd = Td::RowsAtCompileTime;
96 constexpr int Cd = Td::ColsAtCompileTime;
97 constexpr int Ra = Ta::RowsAtCompileTime;
98 constexpr int Ca = Ta::ColsAtCompileTime;
99 constexpr int Rb = Tb::RowsAtCompileTime;
100 constexpr int Cb = Tb::ColsAtCompileTime;
101 check_square("trace_gen_quad_form", "A", A);
102 check_square("trace_gen_quad_form", "D", D);
103 check_multiplicable("trace_gen_quad_form", "A", A, "B", B);
104 check_multiplicable("trace_gen_quad_form", "B", B, "D", D);
105
106 auto* baseVari
107 = new internal::trace_gen_quad_form_vari_alloc<Td_scal, Rd, Cd, Ta_scal,
108 Ra, Ca, Tb_scal, Rb, Cb>(
109 D, A, B);
110
111 return var(
112 new internal::trace_gen_quad_form_vari<Td_scal, Rd, Cd, Ta_scal, Ra, Ca,
113 Tb_scal, Rb, Cb>(baseVari));
114}
115
135template <typename Td, typename Ta, typename Tb,
138inline var trace_gen_quad_form(const Td& D, const Ta& A, const Tb& B) {
139 check_square("trace_gen_quad_form", "A", A);
140 check_square("trace_gen_quad_form", "D", D);
141 check_multiplicable("trace_gen_quad_form", "A", A, "B", B);
142 check_multiplicable("trace_gen_quad_form", "B", B, "D", D);
143
144 if constexpr (is_all_autodiff_v<Ta, Tb, Td>) {
148
149 auto arena_BDT = to_arena(arena_B.val_op() * arena_D.val_op().transpose());
150 auto arena_AB = to_arena(arena_A.val_op() * arena_B.val_op());
151
152 var res = (arena_BDT.transpose() * arena_AB).trace();
153
155 [arena_A, arena_B, arena_D, arena_BDT, arena_AB, res]() mutable {
156 double C_adj = res.adj();
157
158 arena_A.adj() += C_adj * arena_BDT * arena_B.val_op().transpose();
159
160 arena_B.adj() += C_adj
161 * (arena_AB * arena_D.val_op()
162 + arena_A.val_op().transpose() * arena_BDT);
163
164 arena_D.adj() += C_adj * (arena_AB.transpose() * arena_B.val_op());
165 });
166
167 return res;
168 } else if constexpr (is_all_autodiff_v<Ta, Tb> && is_constant_v<Td>) {
169 arena_t<promote_scalar_t<double, Td>> arena_D = value_of(D);
170 arena_t<promote_scalar_t<var, Ta>> arena_A = A;
171 arena_t<promote_scalar_t<var, Tb>> arena_B = B;
172
173 auto arena_BDT = to_arena(arena_B.val_op() * arena_D.transpose());
174 auto arena_AB = to_arena(arena_A.val_op() * arena_B.val_op());
175
176 var res = (arena_BDT.transpose() * arena_AB).trace();
177
178 reverse_pass_callback([arena_A, arena_B, arena_D, arena_BDT, arena_AB,
179 res]() mutable {
180 double C_adj = res.adj();
181
182 arena_A.adj() += C_adj * arena_BDT * arena_B.val_op().transpose();
183 arena_B.adj()
184 += C_adj
185 * (arena_AB * arena_D + arena_A.val_op().transpose() * arena_BDT);
186 });
187
188 return res;
189 } else if constexpr (is_all_autodiff_v<Ta, Td> && is_constant_v<Tb>) {
190 arena_t<promote_scalar_t<var, Td>> arena_D = D;
191 arena_t<promote_scalar_t<var, Ta>> arena_A = A;
192 arena_t<promote_scalar_t<double, Tb>> arena_B = value_of(B);
193
194 auto arena_BDT = to_arena(arena_B.val_op() * arena_D.val_op().transpose());
195 auto arena_AB = to_arena(arena_A.val_op() * arena_B.val_op());
196
197 var res = (arena_BDT.transpose() * arena_A.val_op() * arena_B).trace();
198
200 [arena_A, arena_B, arena_D, arena_BDT, arena_AB, res]() mutable {
201 double C_adj = res.adj();
202
203 arena_A.adj() += C_adj * arena_BDT * arena_B.transpose();
204 arena_D.adj() += C_adj * arena_AB.transpose() * arena_B;
205 });
206
207 return res;
208 } else if constexpr (is_autodiff_v<Ta> && is_constant_all_v<Tb, Td>) {
209 arena_t<promote_scalar_t<double, Td>> arena_D = value_of(D);
210 arena_t<promote_scalar_t<var, Ta>> arena_A = A;
211 arena_t<promote_scalar_t<double, Tb>> arena_B = value_of(B);
212
213 auto arena_BDT = to_arena(arena_B * arena_D);
214
215 var res = (arena_BDT.transpose() * arena_A.val_op() * arena_B).trace();
216
217 reverse_pass_callback([arena_A, arena_B, arena_BDT, res]() mutable {
218 arena_A.adj() += res.adj() * arena_BDT * arena_B.val_op().transpose();
219 });
220
221 return res;
222 } else if constexpr (is_constant_v<Ta> && is_all_autodiff_v<Tb, Td>) {
223 arena_t<promote_scalar_t<var, Td>> arena_D = D;
224 arena_t<promote_scalar_t<double, Ta>> arena_A = value_of(A);
225 arena_t<promote_scalar_t<var, Tb>> arena_B = B;
226
227 auto arena_AB = to_arena(arena_A * arena_B.val_op());
228 auto arena_BDT = to_arena(arena_B.val_op() * arena_D.val_op());
229
230 var res = (arena_BDT.transpose() * arena_AB).trace();
231
232 reverse_pass_callback([arena_A, arena_B, arena_D, arena_AB, arena_BDT,
233 res]() mutable {
234 double C_adj = res.adj();
235
236 arena_B.adj()
237 += C_adj
238 * (arena_AB * arena_D.val_op() + arena_A.transpose() * arena_BDT);
239
240 arena_D.adj() += C_adj * (arena_AB.transpose() * arena_B.val_op());
241 });
242
243 return res;
244 } else if constexpr (is_constant_all_v<Ta, Td> && is_autodiff_v<Tb>) {
245 arena_t<promote_scalar_t<double, Td>> arena_D = value_of(D);
246 arena_t<promote_scalar_t<double, Ta>> arena_A = value_of(A);
247 arena_t<promote_scalar_t<var, Tb>> arena_B = B;
248
249 auto arena_AB = to_arena(arena_A * arena_B.val_op());
250 auto arena_BDT = to_arena(arena_B.val_op() * arena_D.val_op());
251
252 var res = (arena_BDT.transpose() * arena_AB).trace();
253
255 [arena_A, arena_B, arena_D, arena_AB, arena_BDT, res]() mutable {
256 arena_B.adj() += res.adj()
257 * (arena_AB * arena_D.val_op()
258 + arena_A.val_op().transpose() * arena_BDT);
259 });
260
261 return res;
262 } else if constexpr (is_constant_all_v<Ta, Tb> && is_autodiff_v<Td>) {
263 arena_t<promote_scalar_t<var, Td>> arena_D = D;
264 arena_t<promote_scalar_t<double, Ta>> arena_A = value_of(A);
265 arena_t<promote_scalar_t<double, Tb>> arena_B = value_of(B);
266
267 auto arena_AB = to_arena(arena_A * arena_B);
268
269 var res = (arena_D.val_op() * arena_B.transpose() * arena_AB).trace();
270
271 reverse_pass_callback([arena_AB, arena_B, arena_D, res]() mutable {
272 arena_D.adj() += res.adj() * (arena_AB.transpose() * arena_B);
273 });
274
275 return res;
276 }
277}
278
279} // namespace math
280} // namespace stan
281#endif
A chainable_alloc is an object which is constructed and destructed normally but the memory lifespan i...
trace_gen_quad_form_vari_alloc(const Eigen::Matrix< Td, Rd, Cd > &D, const Eigen::Matrix< Ta, Ra, Ca > &A, const Eigen::Matrix< Tb, Rb, Cb > &B)
trace_gen_quad_form_vari(trace_gen_quad_form_vari_alloc< Td, Rd, Cd, Ta, Ra, Ca, Tb, Rb, Cb > *impl)
static void computeAdjoints(double adj, const Eigen::Matrix< double, Rd, Cd > &D, const Eigen::Matrix< double, Ra, Ca > &A, const Eigen::Matrix< double, Rb, Cb > &B, Eigen::Matrix< var, Rd, Cd > *varD, Eigen::Matrix< var, Ra, Ca > *varA, Eigen::Matrix< var, Rb, Cb > *varB)
trace_gen_quad_form_vari_alloc< Td, Rd, Cd, Ta, Ra, Ca, Tb, Rb, Cb > * impl_
require_all_t< is_eigen< std::decay_t< Types > >... > require_all_eigen_t
Require all of the types satisfy is_eigen.
Definition is_eigen.hpp:123
require_all_t< is_matrix< std::decay_t< Types > >... > require_all_matrix_t
Require all of the types satisfy is_matrix.
Definition is_matrix.hpp:38
typename value_type< T >::type value_type_t
Helper function for accessing underlying type.
require_any_t< is_var_matrix< std::decay_t< Types > >... > require_any_var_matrix_t
Require any of the types satisfy is_var_matrix.
require_any_t< is_var< std::decay_t< Types > >... > require_any_var_t
Require any of the types satisfy is_var.
Definition is_var.hpp:104
void check_square(const char *function, const char *name, const T_y &y)
Check if the specified matrix is square.
void check_multiplicable(const char *function, const char *name1, const T1 &y1, const char *name2, const T2 &y2)
Check if the matrices can be multiplied.
void reverse_pass_callback(F &&functor)
Puts a callback on the autodiff stack to be called in reverse pass.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
value_type_t< T > trace(const T &m)
Calculates trace (sum of diagonal) of given kernel generator expression.
Definition trace.hpp:22
arena_t< T > to_arena(const T &a)
Converts given argument into a type that either has any dynamic allocation on AD stack or schedules i...
Definition to_arena.hpp:25
auto trace_gen_quad_form(const TD &D, const TA &A, const TB &B)
Return the trace of D times the quadratic form of B and A.
var_value< double > var
Definition var.hpp:1187
typename internal::arena_type_impl< std::decay_t< T > >::type arena_t
Determines a type that can be used in place of T that does any dynamic allocations on the AD stack.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...