| Stan Math Library 2.12.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

mdivide_left.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_HPP
2 #define STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_HPP
3
4 #include <stan/math/prim/mat/fun/Eigen.hpp>
5 #include <stan/math/prim/mat/fun/typedefs.hpp>
6 #include <stan/math/prim/mat/err/check_multiplicable.hpp>
7 #include <stan/math/rev/core.hpp>
8 #include <stan/math/rev/mat/fun/typedefs.hpp>
9 #include <stan/math/prim/mat/err/check_square.hpp>
10 #include <vector>
11
12 namespace stan {
13 namespace math {
14
15 namespace {
16 template <int R1, int C1, int R2, int C2>
17 class mdivide_left_vv_vari : public vari {
18 public:
19 int M_; // A.rows() = A.cols() = B.rows()
20 int N_; // B.cols()
21 double* A_;
22 double* C_;
23 vari** variRefA_;
24 vari** variRefB_;
25 vari** variRefC_;
26
27 mdivide_left_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::memalloc_
34 .alloc(sizeof(double) * A.rows() * A.cols()))),
35 C_(reinterpret_cast<double*>
36 (ChainableStack::memalloc_
37 .alloc(sizeof(double) * B.rows() * B.cols()))),
38 variRefA_(reinterpret_cast<vari**>
39 (ChainableStack::memalloc_
40 .alloc(sizeof(vari*) * A.rows() * A.cols()))),
41 variRefB_(reinterpret_cast<vari**>
42 (ChainableStack::memalloc_
43 .alloc(sizeof(vari*) * B.rows() * B.cols()))),
44 variRefC_(reinterpret_cast<vari**>
45 (ChainableStack::memalloc_
46 .alloc(sizeof(vari*) * B.rows() * B.cols()))) {
47 using Eigen::Matrix;
48 using Eigen::Map;
49
50 size_t pos = 0;
51 for (size_type j = 0; j < M_; j++) {
52 for (size_type i = 0; i < M_; i++) {
53 variRefA_[pos] = A(i, j).vi_;
54 A_[pos++] = A(i, j).val();
55 }
56 }
57
58 pos = 0;
59 for (size_type j = 0; j < N_; j++) {
60 for (size_type i = 0; i < M_; i++) {
61 variRefB_[pos] = B(i, j).vi_;
62 C_[pos++] = B(i, j).val();
63 }
64 }
65
66 Matrix<double, R1, C2> C(M_, N_);
67 C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
68
69 C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
70 .colPivHouseholderQr().solve(C);
71
72 pos = 0;
73 for (size_type j = 0; j < N_; j++) {
74 for (size_type i = 0; i < M_; i++) {
75 C_[pos] = C(i, j);
76 variRefC_[pos] = new vari(C_[pos], false);
77 pos++;
78 }
79 }
80 }
81
82 virtual void chain() {
83 using Eigen::Matrix;
84 using Eigen::Map;
85 Eigen::Matrix<double, R1, C1> adjA(M_, M_);
86 Eigen::Matrix<double, R2, C2> adjB(M_, N_);
87 Eigen::Matrix<double, R1, C2> adjC(M_, N_);
88
89 size_t pos = 0;
90 for (size_type j = 0; j < adjC.cols(); j++)
91 for (size_type i = 0; i < adjC.rows(); i++)
92 adjC(i, j) = variRefC_[pos++]->adj_;
93
94 adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
95 .transpose().colPivHouseholderQr().solve(adjC);
96 adjA.noalias() = -adjB
97 * Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose();
98
99 pos = 0;
100 for (size_type j = 0; j < adjA.cols(); j++)
101 for (size_type i = 0; i < adjA.rows(); i++)
102 variRefA_[pos++]->adj_ += adjA(i, j);
103
104 pos = 0;
105 for (size_type j = 0; j < adjB.cols(); j++)
106 for (size_type i = 0; i < adjB.rows(); i++)
107 variRefB_[pos++]->adj_ += adjB(i, j);
108 }
109 };
110
111 template <int R1, int C1, int R2, int C2>
112 class mdivide_left_dv_vari : public vari {
113 public:
114 int M_; // A.rows() = A.cols() = B.rows()
115 int N_; // B.cols()
116 double* A_;
117 double* C_;
118 vari** variRefB_;
119 vari** variRefC_;
120
121 mdivide_left_dv_vari(const Eigen::Matrix<double, R1, C1> &A,
122 const Eigen::Matrix<var, R2, C2> &B)
123 : vari(0.0),
124 M_(A.rows()),
125 N_(B.cols()),
126 A_(reinterpret_cast<double*>
127 (ChainableStack::memalloc_
128 .alloc(sizeof(double) * A.rows() * A.cols()))),
129 C_(reinterpret_cast<double*>
130 (ChainableStack::memalloc_
131 .alloc(sizeof(double) * B.rows() * B.cols()))),
132 variRefB_(reinterpret_cast<vari**>
133 (ChainableStack::memalloc_
134 .alloc(sizeof(vari*) * B.rows() * B.cols()))),
135 variRefC_(reinterpret_cast<vari**>
136 (ChainableStack::memalloc_
137 .alloc(sizeof(vari*) * B.rows() * B.cols()))) {
138 using Eigen::Matrix;
139 using Eigen::Map;
140
141 size_t pos = 0;
142 for (size_type j = 0; j < M_; j++) {
143 for (size_type i = 0; i < M_; i++) {
144 A_[pos++] = A(i, j);
145 }
146 }
147
148 pos = 0;
149 for (size_type j = 0; j < N_; j++) {
150 for (size_type i = 0; i < M_; i++) {
151 variRefB_[pos] = B(i, j).vi_;
152 C_[pos++] = B(i, j).val();
153 }
154 }
155
156 Matrix<double, R1, C2> C(M_, N_);
157 C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
158
159 C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
160 .colPivHouseholderQr().solve(C);
161
162 pos = 0;
163 for (size_type j = 0; j < N_; j++) {
164 for (size_type i = 0; i < M_; i++) {
165 C_[pos] = C(i, j);
166 variRefC_[pos] = new vari(C_[pos], false);
167 pos++;
168 }
169 }
170 }
171
172 virtual void chain() {
173 using Eigen::Matrix;
174 using Eigen::Map;
175 Eigen::Matrix<double, R2, C2> adjB(M_, N_);
176 Eigen::Matrix<double, R1, C2> adjC(M_, N_);
177
178 size_t pos = 0;
179 for (size_type j = 0; j < adjC.cols(); j++)
180 for (size_type i = 0; i < adjC.rows(); i++)
181 adjC(i, j) = variRefC_[pos++]->adj_;
182
183 adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
184 .transpose().colPivHouseholderQr().solve(adjC);
185
186 pos = 0;
187 for (size_type j = 0; j < adjB.cols(); j++)
188 for (size_type i = 0; i < adjB.rows(); i++)
189 variRefB_[pos++]->adj_ += adjB(i, j);
190 }
191 };
192
193 template <int R1, int C1, int R2, int C2>
194 class mdivide_left_vd_vari : public vari {
195 public:
196 int M_; // A.rows() = A.cols() = B.rows()
197 int N_; // B.cols()
198 double* A_;
199 double* C_;
200 vari** variRefA_;
201 vari** variRefC_;
202
203 mdivide_left_vd_vari(const Eigen::Matrix<var, R1, C1> &A,
204 const Eigen::Matrix<double, R2, C2> &B)
205 : vari(0.0),
206 M_(A.rows()),
207 N_(B.cols()),
208 A_(reinterpret_cast<double*>
209 (ChainableStack::memalloc_
210 .alloc(sizeof(double) * A.rows() * A.cols()))),
211 C_(reinterpret_cast<double*>
212 (ChainableStack::memalloc_
213 .alloc(sizeof(double) * B.rows() * B.cols()))),
214 variRefA_(reinterpret_cast<vari**>
215 (ChainableStack::memalloc_
216 .alloc(sizeof(vari*) * A.rows() * A.cols()))),
217 variRefC_(reinterpret_cast<vari**>
218 (ChainableStack::memalloc_
219 .alloc(sizeof(vari*) * B.rows() * B.cols()))) {
220 using Eigen::Matrix;
221 using Eigen::Map;
222
223 size_t pos = 0;
224 for (size_type j = 0; j < M_; j++) {
225 for (size_type i = 0; i < M_; i++) {
226 variRefA_[pos] = A(i, j).vi_;
227 A_[pos++] = A(i, j).val();
228 }
229 }
230
231 Matrix<double, R1, C2> C(M_, N_);
232 C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
233 .colPivHouseholderQr().solve(B);
234
235 pos = 0;
236 for (size_type j = 0; j < N_; j++) {
237 for (size_type i = 0; i < M_; i++) {
238 C_[pos] = C(i, j);
239 variRefC_[pos] = new vari(C_[pos], false);
240 pos++;
241 }
242 }
243 }
244
245 virtual void chain() {
246 using Eigen::Matrix;
247 using Eigen::Map;
248 Eigen::Matrix<double, R1, C1> adjA(M_, M_);
249 Eigen::Matrix<double, R1, C2> adjC(M_, N_);
250
251 size_t pos = 0;
252 for (size_type j = 0; j < adjC.cols(); j++)
253 for (size_type i = 0; i < adjC.rows(); i++)
254 adjC(i, j) = variRefC_[pos++]->adj_;
255
256 // FIXME: add .noalias() to LHS
257 adjA = -Map<Matrix<double, R1, C1> >(A_, M_, M_)
258 .transpose()
259 .colPivHouseholderQr()
260 .solve(adjC*Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose());
261
262 pos = 0;
263 for (size_type j = 0; j < adjA.cols(); j++)
264 for (size_type i = 0; i < adjA.rows(); i++)
265 variRefA_[pos++]->adj_ += adjA(i, j);
266 }
267 };
268 }
269
270 template <int R1, int C1, int R2, int C2>
271 inline
272 Eigen::Matrix<var, R1, C2>
273 mdivide_left(const Eigen::Matrix<var, R1, C1> &A,
274 const Eigen::Matrix<var, R2, C2> &b) {
275 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
276
277 check_square("mdivide_left", "A", A);
278 check_multiplicable("mdivide_left",
279 "A", A,
280 "b", b);
281
282 // NOTE: this is not a memory leak, this vari is used in the
283 // expression graph to evaluate the adjoint, but is not needed
284 // for the returned matrix. Memory will be cleaned up with the
285 // arena allocator.
286 mdivide_left_vv_vari<R1, C1, R2, C2> *baseVari
287 = new mdivide_left_vv_vari<R1, C1, R2, C2>(A, b);
288
289 size_t pos = 0;
290 for (size_type j = 0; j < res.cols(); j++)
291 for (size_type i = 0; i < res.rows(); i++)
292 res(i, j).vi_ = baseVari->variRefC_[pos++];
293
294 return res;
295 }
296
297 template <int R1, int C1, int R2, int C2>
298 inline
299 Eigen::Matrix<var, R1, C2>
300 mdivide_left(const Eigen::Matrix<var, R1, C1> &A,
301 const Eigen::Matrix<double, R2, C2> &b) {
302 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
303
304 check_square("mdivide_left", "A", A);
305 check_multiplicable("mdivide_left",
306 "A", A,
307 "b", b);
308
309 // NOTE: this is not a memory leak, this vari is used in the
310 // expression graph to evaluate the adjoint, but is not needed
311 // for the returned matrix. Memory will be cleaned up with the
312 // arena allocator.
313 mdivide_left_vd_vari<R1, C1, R2, C2> *baseVari
314 = new mdivide_left_vd_vari<R1, C1, R2, C2>(A, b);
315
316 size_t pos = 0;
317 for (size_type j = 0; j < res.cols(); j++)
318 for (size_type i = 0; i < res.rows(); i++)
319 res(i, j).vi_ = baseVari->variRefC_[pos++];
320
321 return res;
322 }
323
324 template <int R1, int C1, int R2, int C2>
325 inline
326 Eigen::Matrix<var, R1, C2>
327 mdivide_left(const Eigen::Matrix<double, R1, C1> &A,
328 const Eigen::Matrix<var, R2, C2> &b) {
329 Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
330
331 check_square("mdivide_left", "A", A);
332 check_multiplicable("mdivide_left",
333 "A", A,
334 "b", b);
335
336 // NOTE: this is not a memory leak, this vari is used in the
337 // expression graph to evaluate the adjoint, but is not needed
338 // for the returned matrix. Memory will be cleaned up with the
339 // arena allocator.
340 mdivide_left_dv_vari<R1, C1, R2, C2> *baseVari
341 = new mdivide_left_dv_vari<R1, C1, R2, C2>(A, b);
342
343 size_t pos = 0;
344 for (size_type j = 0; j < res.cols(); j++)
345 for (size_type i = 0; i < res.rows(); i++)
346 res(i, j).vi_ = baseVari->variRefC_[pos++];
347
348 return res;
349 }
350
351 }
352 }
353 #endif
stan::math::mdivide_left
Eigen::Matrix< fvar< T >, R1, C2 > mdivide_left(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
Definition: mdivide_left.hpp:24
stan::math::rows
int rows(const Eigen::Matrix< T, R, C > &m)
Return the number of rows in the specified matrix, vector, or row vector.
Definition: rows.hpp:20
core.hpp
check_multiplicable.hpp
C_
double * C_
Definition: mdivide_left.hpp:22
stan
Definition: log_sum_exp.hpp:8
variRefB_
vari ** variRefB_
Definition: mdivide_left.hpp:24
stan::math::size_type
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Type for sizes and indexes in an Eigen matrix with double e.
Definition: typedefs.hpp:13
M_
int M_
Definition: mdivide_left.hpp:19
typedefs.hpp
stan::math::check_multiplicable
bool check_multiplicable(const char *function, const char *name1, const T1 &y1, const char *name2, const T2 &y2)
Return true if the matrices can be multiplied.
Definition: check_multiplicable.hpp:32
stan::math::cols
int cols(const Eigen::Matrix< T, R, C > &m)
Return the number of columns in the specified matrix, vector, or row vector.
Definition: cols.hpp:20
Eigen.hpp
A_
double * A_
Definition: mdivide_left.hpp:21
variRefA_
vari ** variRefA_
Definition: mdivide_left.hpp:23
variRefC_
vari ** variRefC_
Definition: mdivide_left.hpp:25
stan::math::transpose
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:12
stan::math::check_square
bool check_square(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Return true if the specified matrix is square.
Definition: check_square.hpp:28
stan::math::ChainableStack
AutodiffStackStorage< vari, chainable_alloc > ChainableStack
Definition: chainablestack.hpp:10
N_
int N_
Definition: mdivide_left.hpp:20
typedefs.hpp
check_square.hpp

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