Automatic Differentiation
 
Loading...
Searching...
No Matches
mdivide_left_tri.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_REV_FUN_MDIVIDE_LEFT_TRI_HPP
2#define STAN_MATH_REV_FUN_MDIVIDE_LEFT_TRI_HPP
3
10
11namespace stan {
12namespace math {
13
14namespace internal {
15
16template <Eigen::UpLoType TriView, int R1, int C1, int R2, int C2>
18 public:
19 int M_; // A.rows() = A.cols() = B.rows()
20 int N_; // B.cols()
21 double *A_;
22 double *C_;
26
27 mdivide_left_tri_vv_vari(const Eigen::Matrix<var, R1, C1> &A,
28 const Eigen::Matrix<var, R2, C2> &B)
29 : vari(0.0),
30 M_(A.rows()),
31 N_(B.cols()),
32 A_(reinterpret_cast<double *>(
33 ChainableStack::instance_->memalloc_.alloc(sizeof(double) * A.rows()
34 * A.cols()))),
35 C_(reinterpret_cast<double *>(
36 ChainableStack::instance_->memalloc_.alloc(sizeof(double) * B.rows()
37 * B.cols()))),
38 variRefA_(reinterpret_cast<vari **>(
39 ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * A.rows()
40 * (A.rows() + 1) / 2))),
41 variRefB_(reinterpret_cast<vari **>(
42 ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
43 * B.cols()))),
44 variRefC_(reinterpret_cast<vari **>(
45 ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
46 * B.cols()))) {
47 using Eigen::Map;
48
49 size_t pos = 0;
50 if (TriView == Eigen::Lower) {
51 for (size_type j = 0; j < M_; j++) {
52 for (size_type i = j; i < M_; i++) {
53 variRefA_[pos++] = A(i, j).vi_;
54 }
55 }
56 } else if (TriView == Eigen::Upper) {
57 for (size_type j = 0; j < M_; j++) {
58 for (size_type i = 0; i < j + 1; i++) {
59 variRefA_[pos++] = A(i, j).vi_;
60 }
61 }
62 }
63
64 Map<matrix_d> c_map(C_, M_, N_);
65 Map<matrix_d> a_map(A_, M_, M_);
66 a_map = A.val();
67 c_map = B.val();
68 Map<matrix_vi>(variRefB_, M_, N_) = B.vi();
69
70 c_map = a_map.template triangularView<TriView>().solve(c_map);
71
72 Map<matrix_vi>(variRefC_, M_, N_)
73 = c_map.unaryExpr([](double x) { return new vari(x, false); });
74 }
75
76 virtual void chain() {
77 using Eigen::Map;
78 matrix_d adjA;
79 matrix_d adjB;
80
81 adjB = Map<matrix_d>(A_, M_, M_)
82 .template triangularView<TriView>()
83 .transpose()
84 .solve(Map<matrix_vi>(variRefC_, M_, N_).adj());
85 adjA = -adjB * Map<matrix_d>(C_, M_, N_).transpose();
86
87 size_t pos = 0;
88 if (TriView == Eigen::Lower) {
89 for (size_type j = 0; j < adjA.cols(); j++) {
90 for (size_type i = j; i < adjA.rows(); i++) {
91 variRefA_[pos++]->adj_ += adjA(i, j);
92 }
93 }
94 } else if (TriView == Eigen::Upper) {
95 for (size_type j = 0; j < adjA.cols(); j++) {
96 for (size_type i = 0; i < j + 1; i++) {
97 variRefA_[pos++]->adj_ += adjA(i, j);
98 }
99 }
100 }
101 Map<matrix_vi>(variRefB_, M_, N_).adj() += adjB;
102 }
103};
104
105template <Eigen::UpLoType TriView, int R1, int C1, int R2, int C2>
107 public:
108 int M_; // A.rows() = A.cols() = B.rows()
109 int N_; // B.cols()
110 double *A_;
111 double *C_;
114
115 mdivide_left_tri_dv_vari(const Eigen::Matrix<double, R1, C1> &A,
116 const Eigen::Matrix<var, R2, C2> &B)
117 : vari(0.0),
118 M_(A.rows()),
119 N_(B.cols()),
120 A_(reinterpret_cast<double *>(
121 ChainableStack::instance_->memalloc_.alloc(sizeof(double) * A.rows()
122 * A.cols()))),
123 C_(reinterpret_cast<double *>(
124 ChainableStack::instance_->memalloc_.alloc(sizeof(double) * B.rows()
125 * B.cols()))),
126 variRefB_(reinterpret_cast<vari **>(
127 ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
128 * B.cols()))),
129 variRefC_(reinterpret_cast<vari **>(
130 ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
131 * B.cols()))) {
132 using Eigen::Map;
133
134 Map<matrix_d>(A_, M_, M_) = A;
135 Map<matrix_vi>(variRefB_, M_, N_) = B.vi();
136 Map<matrix_d> c_map(C_, M_, N_);
137 c_map = B.val();
138
139 c_map = Map<matrix_d>(A_, M_, M_)
140 .template triangularView<TriView>()
141 .solve(c_map);
142
143 Map<matrix_vi>(variRefC_, M_, N_)
144 = c_map.unaryExpr([](double x) { return new vari(x, false); });
145 }
146
147 virtual void chain() {
148 using Eigen::Map;
149
150 Map<matrix_vi>(variRefB_, M_, N_).adj()
151 += Map<matrix_d>(A_, M_, M_)
152 .template triangularView<TriView>()
153 .transpose()
154 .solve(Map<matrix_vi>(variRefC_, M_, N_).adj());
155 }
156};
157
158template <Eigen::UpLoType TriView, int R1, int C1, int R2, int C2>
160 public:
161 int M_; // A.rows() = A.cols() = B.rows()
162 int N_; // B.cols()
163 double *A_;
164 double *C_;
167
168 mdivide_left_tri_vd_vari(const Eigen::Matrix<var, R1, C1> &A,
169 const Eigen::Matrix<double, R2, C2> &B)
170 : vari(0.0),
171 M_(A.rows()),
172 N_(B.cols()),
173 A_(reinterpret_cast<double *>(
174 ChainableStack::instance_->memalloc_.alloc(sizeof(double) * A.rows()
175 * A.cols()))),
176 C_(reinterpret_cast<double *>(
177 ChainableStack::instance_->memalloc_.alloc(sizeof(double) * B.rows()
178 * B.cols()))),
179 variRefA_(reinterpret_cast<vari **>(
180 ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * A.rows()
181 * (A.rows() + 1) / 2))),
182 variRefC_(reinterpret_cast<vari **>(
183 ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
184 * B.cols()))) {
185 using Eigen::Map;
186 using Eigen::Matrix;
187
188 size_t pos = 0;
189 if (TriView == Eigen::Lower) {
190 for (size_type j = 0; j < M_; j++) {
191 for (size_type i = j; i < M_; i++) {
192 variRefA_[pos++] = A(i, j).vi_;
193 }
194 }
195 } else if (TriView == Eigen::Upper) {
196 for (size_type j = 0; j < M_; j++) {
197 for (size_type i = 0; i < j + 1; i++) {
198 variRefA_[pos++] = A(i, j).vi_;
199 }
200 }
201 }
202 Map<matrix_d> Ad(A_, M_, M_);
203 Map<matrix_d> Cd(C_, M_, N_);
204 Ad = A.val();
205
206 Cd = Ad.template triangularView<TriView>().solve(B);
207
208 Map<matrix_vi>(variRefC_, M_, N_)
209 = Cd.unaryExpr([](double x) { return new vari(x, false); });
210 }
211
212 virtual void chain() {
213 using Eigen::Map;
214 using Eigen::Matrix;
215 Matrix<double, R1, C1> adjA(M_, M_);
216 Matrix<double, R1, C2> adjC(M_, N_);
217
218 adjC = Map<matrix_vi>(variRefC_, M_, N_).adj();
219
220 adjA.noalias()
221 = -Map<Matrix<double, R1, C1>>(A_, M_, M_)
222 .template triangularView<TriView>()
223 .transpose()
224 .solve(adjC
225 * Map<Matrix<double, R1, C2>>(C_, M_, N_).transpose());
226
227 size_t pos = 0;
228 if (TriView == Eigen::Lower) {
229 for (size_type j = 0; j < adjA.cols(); j++) {
230 for (size_type i = j; i < adjA.rows(); i++) {
231 variRefA_[pos++]->adj_ += adjA(i, j);
232 }
233 }
234 } else if (TriView == Eigen::Upper) {
235 for (size_type j = 0; j < adjA.cols(); j++) {
236 for (size_type i = 0; i < j + 1; i++) {
237 variRefA_[pos++]->adj_ += adjA(i, j);
238 }
239 }
240 }
241 }
242};
243} // namespace internal
244
245template <Eigen::UpLoType TriView, typename T1, typename T2,
247inline Eigen::Matrix<var, T1::RowsAtCompileTime, T2::ColsAtCompileTime>
248mdivide_left_tri(const T1 &A, const T2 &b) {
249 check_square("mdivide_left_tri", "A", A);
250 check_multiplicable("mdivide_left_tri", "A", A, "b", b);
251 if (A.rows() == 0) {
252 return {0, b.cols()};
253 }
254
255 // NOTE: this is not a memory leak, this vari is used in the
256 // expression graph to evaluate the adjoint, but is not needed
257 // for the returned matrix. Memory will be cleaned up with the
258 // arena allocator.
259 auto *baseVari = new internal::mdivide_left_tri_vv_vari<
260 TriView, T1::RowsAtCompileTime, T1::ColsAtCompileTime,
261 T2::RowsAtCompileTime, T2::ColsAtCompileTime>(A, b);
262
263 Eigen::Matrix<var, T1::RowsAtCompileTime, T2::ColsAtCompileTime> res(
264 b.rows(), b.cols());
265 res.vi()
266 = Eigen::Map<matrix_vi>(&(baseVari->variRefC_[0]), b.rows(), b.cols());
267
268 return res;
269}
270template <Eigen::UpLoType TriView, typename T1, typename T2,
273inline Eigen::Matrix<var, T1::RowsAtCompileTime, T2::ColsAtCompileTime>
274mdivide_left_tri(const T1 &A, const T2 &b) {
275 check_square("mdivide_left_tri", "A", A);
276 check_multiplicable("mdivide_left_tri", "A", A, "b", b);
277 if (A.rows() == 0) {
278 return {0, b.cols()};
279 }
280
281 // NOTE: this is not a memory leak, this vari is used in the
282 // expression graph to evaluate the adjoint, but is not needed
283 // for the returned matrix. Memory will be cleaned up with the
284 // arena allocator.
285 auto *baseVari = new internal::mdivide_left_tri_dv_vari<
286 TriView, T1::RowsAtCompileTime, T1::ColsAtCompileTime,
287 T2::RowsAtCompileTime, T2::ColsAtCompileTime>(A, b);
288
289 Eigen::Matrix<var, T1::RowsAtCompileTime, T2::ColsAtCompileTime> res(
290 b.rows(), b.cols());
291 res.vi()
292 = Eigen::Map<matrix_vi>(&(baseVari->variRefC_[0]), b.rows(), b.cols());
293
294 return res;
295}
296template <Eigen::UpLoType TriView, typename T1, typename T2,
299inline Eigen::Matrix<var, T1::RowsAtCompileTime, T2::ColsAtCompileTime>
300mdivide_left_tri(const T1 &A, const T2 &b) {
301 check_square("mdivide_left_tri", "A", A);
302 check_multiplicable("mdivide_left_tri", "A", A, "b", b);
303 if (A.rows() == 0) {
304 return {0, b.cols()};
305 }
306
307 // NOTE: this is not a memory leak, this vari is used in the
308 // expression graph to evaluate the adjoint, but is not needed
309 // for the returned matrix. Memory will be cleaned up with the
310 // arena allocator.
311 auto *baseVari = new internal::mdivide_left_tri_vd_vari<
312 TriView, T1::RowsAtCompileTime, T1::ColsAtCompileTime,
313 T2::RowsAtCompileTime, T2::ColsAtCompileTime>(A, b);
314
315 Eigen::Matrix<var, T1::RowsAtCompileTime, T2::ColsAtCompileTime> res(
316 b.rows(), b.cols());
317 res.vi()
318 = Eigen::Map<matrix_vi>(&(baseVari->variRefC_[0]), b.rows(), b.cols());
319
320 return res;
321}
322
342template <Eigen::UpLoType TriView, typename T1, typename T2,
343 require_all_matrix_t<T1, T2> * = nullptr,
344 require_any_var_matrix_t<T1, T2> * = nullptr>
345inline auto mdivide_left_tri(const T1 &A, const T2 &B) {
346 using ret_val_type = plain_type_t<decltype(value_of(A) * value_of(B))>;
347 using ret_type = var_value<ret_val_type>;
348
349 if (A.size() == 0) {
350 return ret_type(ret_val_type(0, B.cols()));
351 }
352
353 check_square("mdivide_left_tri", "A", A);
354 check_multiplicable("mdivide_left_tri", "A", A, "B", B);
355
356 if (!is_constant<T1>::value && !is_constant<T2>::value) {
357 arena_t<promote_scalar_t<var, T1>> arena_A = A;
358 arena_t<promote_scalar_t<var, T2>> arena_B = B;
359 auto arena_A_val = to_arena(arena_A.val());
360
361 arena_t<ret_type> res
362 = arena_A_val.template triangularView<TriView>().solve(arena_B.val());
363
364 reverse_pass_callback([arena_A, arena_B, arena_A_val, res]() mutable {
365 promote_scalar_t<double, T2> adjB
366 = arena_A_val.template triangularView<TriView>().transpose().solve(
367 res.adj());
368
369 arena_B.adj() += adjB;
370 arena_A.adj() -= (adjB * res.val().transpose().eval())
371 .template triangularView<TriView>();
372 });
373
374 return ret_type(res);
375 } else if (!is_constant<T1>::value) {
376 arena_t<promote_scalar_t<var, T1>> arena_A = A;
377 auto arena_A_val = to_arena(arena_A.val());
378
379 arena_t<ret_type> res
380 = arena_A_val.template triangularView<TriView>().solve(value_of(B));
381
382 reverse_pass_callback([arena_A, arena_A_val, res]() mutable {
383 promote_scalar_t<double, T2> adjB
384 = arena_A_val.template triangularView<TriView>().transpose().solve(
385 res.adj());
386
387 arena_A.adj() -= (adjB * res.val().transpose().eval())
388 .template triangularView<TriView>();
389 });
390
391 return ret_type(res);
392 } else {
393 arena_t<promote_scalar_t<double, T1>> arena_A = value_of(A);
394 arena_t<promote_scalar_t<var, T2>> arena_B = B;
395
396 arena_t<ret_type> res
397 = arena_A.template triangularView<TriView>().solve(arena_B.val());
398
399 reverse_pass_callback([arena_A, arena_B, res]() mutable {
400 promote_scalar_t<double, T2> adjB
401 = arena_A.template triangularView<TriView>().transpose().solve(
402 res.adj());
403
404 arena_B.adj() += adjB;
405 });
406
407 return ret_type(res);
408 }
409}
410
411} // namespace math
412} // namespace stan
413#endif
mdivide_left_tri_dv_vari(const Eigen::Matrix< double, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)
mdivide_left_tri_vd_vari(const Eigen::Matrix< var, R1, C1 > &A, const Eigen::Matrix< double, R2, C2 > &B)
mdivide_left_tri_vv_vari(const Eigen::Matrix< var, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)
require_all_t< container_type_check_base< is_eigen, value_type_t, TypeCheck, Check >... > require_all_eigen_vt
Require all of the types satisfy is_eigen.
Definition is_eigen.hpp:133
require_t< container_type_check_base< is_eigen, value_type_t, TypeCheck, Check... > > require_eigen_vt
Require type satisfies is_eigen.
Definition is_eigen.hpp:97
int64_t cols(const T_x &x)
Returns the number of columns in the specified kernel generator expression.
Definition cols.hpp:21
int64_t rows(const T_x &x)
Returns the number of rows in the specified kernel generator expression.
Definition rows.hpp:22
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
vari_value< double > vari
Definition vari.hpp:197
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
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Type for sizes and indexes in an Eigen matrix with double elements.
Definition typedefs.hpp:11
auto mdivide_left_tri(const T1 &A, const T2 &b)
Returns the solution of the system Ax=b when A is triangular.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > matrix_d
Type for matrix of double values.
Definition typedefs.hpp:19
typename plain_type< T >::type plain_type_t
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.