Stan Math Library  2.9.0
reverse mode automatic differentiation
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 
7 #include <stan/math/rev/core.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  (stan::math::ChainableStack::memalloc_
34  .alloc(sizeof(double) * A.rows() * A.cols()))),
35  C_(reinterpret_cast<double*>
36  (stan::math::ChainableStack::memalloc_
37  .alloc(sizeof(double) * B.rows() * B.cols()))),
38  _variRefA(reinterpret_cast<vari**>
39  (stan::math::ChainableStack::memalloc_
40  .alloc(sizeof(vari*) * A.rows() * A.cols()))),
41  _variRefB(reinterpret_cast<vari**>
42  (stan::math::ChainableStack::memalloc_
43  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
44  _variRefC(reinterpret_cast<vari**>
45  (stan::math::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 
95  adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
96  .transpose().colPivHouseholderQr().solve(adjC);
97  adjA.noalias() = -adjB
98  * Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose();
99 
100  pos = 0;
101  for (size_type j = 0; j < adjA.cols(); j++)
102  for (size_type i = 0; i < adjA.rows(); i++)
103  _variRefA[pos++]->adj_ += adjA(i, j);
104 
105  pos = 0;
106  for (size_type j = 0; j < adjB.cols(); j++)
107  for (size_type i = 0; i < adjB.rows(); i++)
108  _variRefB[pos++]->adj_ += adjB(i, j);
109  }
110  };
111 
112  template <int R1, int C1, int R2, int C2>
113  class mdivide_left_dv_vari : public vari {
114  public:
115  int M_; // A.rows() = A.cols() = B.rows()
116  int N_; // B.cols()
117  double* A_;
118  double* C_;
119  vari** _variRefB;
120  vari** _variRefC;
121 
122  mdivide_left_dv_vari(const Eigen::Matrix<double, R1, C1> &A,
123  const Eigen::Matrix<var, R2, C2> &B)
124  : vari(0.0),
125  M_(A.rows()),
126  N_(B.cols()),
127  A_(reinterpret_cast<double*>
128  (stan::math::ChainableStack::memalloc_
129  .alloc(sizeof(double) * A.rows() * A.cols()))),
130  C_(reinterpret_cast<double*>
131  (stan::math::ChainableStack::memalloc_
132  .alloc(sizeof(double) * B.rows() * B.cols()))),
133  _variRefB(reinterpret_cast<vari**>
134  (stan::math::ChainableStack::memalloc_
135  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
136  _variRefC(reinterpret_cast<vari**>
137  (stan::math::ChainableStack::memalloc_
138  .alloc(sizeof(vari*) * B.rows() * B.cols()))) {
139  using Eigen::Matrix;
140  using Eigen::Map;
141 
142  size_t pos = 0;
143  for (size_type j = 0; j < M_; j++) {
144  for (size_type i = 0; i < M_; i++) {
145  A_[pos++] = A(i, j);
146  }
147  }
148 
149  pos = 0;
150  for (size_type j = 0; j < N_; j++) {
151  for (size_type i = 0; i < M_; i++) {
152  _variRefB[pos] = B(i, j).vi_;
153  C_[pos++] = B(i, j).val();
154  }
155  }
156 
157  Matrix<double, R1, C2> C(M_, N_);
158  C = Map<Matrix<double, R1, C2> >(C_, M_, N_);
159 
160  C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
161  .colPivHouseholderQr().solve(C);
162 
163  pos = 0;
164  for (size_type j = 0; j < N_; j++) {
165  for (size_type i = 0; i < M_; i++) {
166  C_[pos] = C(i, j);
167  _variRefC[pos] = new vari(C_[pos], false);
168  pos++;
169  }
170  }
171  }
172 
173  virtual void chain() {
174  using Eigen::Matrix;
175  using Eigen::Map;
176  Eigen::Matrix<double, R2, C2> adjB(M_, N_);
177  Eigen::Matrix<double, R1, C2> adjC(M_, N_);
178 
179  size_t pos = 0;
180  for (size_type j = 0; j < adjC.cols(); j++)
181  for (size_type i = 0; i < adjC.rows(); i++)
182  adjC(i, j) = _variRefC[pos++]->adj_;
183 
184  adjB = Map<Matrix<double, R1, C1> >(A_, M_, M_)
185  .transpose().colPivHouseholderQr().solve(adjC);
186 
187  pos = 0;
188  for (size_type j = 0; j < adjB.cols(); j++)
189  for (size_type i = 0; i < adjB.rows(); i++)
190  _variRefB[pos++]->adj_ += adjB(i, j);
191  }
192  };
193 
194  template <int R1, int C1, int R2, int C2>
195  class mdivide_left_vd_vari : public vari {
196  public:
197  int M_; // A.rows() = A.cols() = B.rows()
198  int N_; // B.cols()
199  double* A_;
200  double* C_;
201  vari** _variRefA;
202  vari** _variRefC;
203 
204  mdivide_left_vd_vari(const Eigen::Matrix<var, R1, C1> &A,
205  const Eigen::Matrix<double, R2, C2> &B)
206  : vari(0.0),
207  M_(A.rows()),
208  N_(B.cols()),
209  A_(reinterpret_cast<double*>
210  (stan::math::ChainableStack::memalloc_
211  .alloc(sizeof(double) * A.rows() * A.cols()))),
212  C_(reinterpret_cast<double*>
213  (stan::math::ChainableStack::memalloc_
214  .alloc(sizeof(double) * B.rows() * B.cols()))),
215  _variRefA(reinterpret_cast<vari**>
216  (stan::math::ChainableStack::memalloc_
217  .alloc(sizeof(vari*) * A.rows() * A.cols()))),
218  _variRefC(reinterpret_cast<vari**>
219  (stan::math::ChainableStack::memalloc_
220  .alloc(sizeof(vari*) * B.rows() * B.cols()))) {
221  using Eigen::Matrix;
222  using Eigen::Map;
223 
224  size_t pos = 0;
225  for (size_type j = 0; j < M_; j++) {
226  for (size_type i = 0; i < M_; i++) {
227  _variRefA[pos] = A(i, j).vi_;
228  A_[pos++] = A(i, j).val();
229  }
230  }
231 
232  Matrix<double, R1, C2> C(M_, N_);
233  C = Map<Matrix<double, R1, C1> >(A_, M_, M_)
234  .colPivHouseholderQr().solve(B);
235 
236  pos = 0;
237  for (size_type j = 0; j < N_; j++) {
238  for (size_type i = 0; i < M_; i++) {
239  C_[pos] = C(i, j);
240  _variRefC[pos] = new vari(C_[pos], false);
241  pos++;
242  }
243  }
244  }
245 
246  virtual void chain() {
247  using Eigen::Matrix;
248  using Eigen::Map;
249  Eigen::Matrix<double, R1, C1> adjA(M_, M_);
250  Eigen::Matrix<double, R1, C2> adjC(M_, N_);
251 
252  size_t pos = 0;
253  for (size_type j = 0; j < adjC.cols(); j++)
254  for (size_type i = 0; i < adjC.rows(); i++)
255  adjC(i, j) = _variRefC[pos++]->adj_;
256 
257  // FIXME: add .noalias() to LHS
258  adjA = -Map<Matrix<double, R1, C1> >(A_, M_, M_)
259  .transpose()
260  .colPivHouseholderQr()
261  .solve(adjC*Map<Matrix<double, R1, C2> >(C_, M_, N_).transpose());
262 
263  pos = 0;
264  for (size_type j = 0; j < adjA.cols(); j++)
265  for (size_type i = 0; i < adjA.rows(); i++)
266  _variRefA[pos++]->adj_ += adjA(i, j);
267  }
268  };
269  }
270 
271  template <int R1, int C1, int R2, int C2>
272  inline
273  Eigen::Matrix<var, R1, C2>
274  mdivide_left(const Eigen::Matrix<var, R1, C1> &A,
275  const Eigen::Matrix<var, R2, C2> &b) {
276  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
277 
278  stan::math::check_square("mdivide_left", "A", A);
279  stan::math::check_multiplicable("mdivide_left",
280  "A", A,
281  "b", b);
282 
283  // NOTE: this is not a memory leak, this vari is used in the
284  // expression graph to evaluate the adjoint, but is not needed
285  // for the returned matrix. Memory will be cleaned up with the
286  // arena allocator.
287  mdivide_left_vv_vari<R1, C1, R2, C2> *baseVari
288  = new mdivide_left_vv_vari<R1, C1, R2, C2>(A, b);
289 
290  size_t pos = 0;
291  for (size_type j = 0; j < res.cols(); j++)
292  for (size_type i = 0; i < res.rows(); i++)
293  res(i, j).vi_ = baseVari->_variRefC[pos++];
294 
295  return res;
296  }
297 
298  template <int R1, int C1, int R2, int C2>
299  inline
300  Eigen::Matrix<var, R1, C2>
301  mdivide_left(const Eigen::Matrix<var, R1, C1> &A,
302  const Eigen::Matrix<double, R2, C2> &b) {
303  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
304 
305  stan::math::check_square("mdivide_left", "A", A);
306  stan::math::check_multiplicable("mdivide_left",
307  "A", A,
308  "b", b);
309 
310  // NOTE: this is not a memory leak, this vari is used in the
311  // expression graph to evaluate the adjoint, but is not needed
312  // for the returned matrix. Memory will be cleaned up with the
313  // arena allocator.
314  mdivide_left_vd_vari<R1, C1, R2, C2> *baseVari
315  = new mdivide_left_vd_vari<R1, C1, R2, C2>(A, b);
316 
317  size_t pos = 0;
318  for (size_type j = 0; j < res.cols(); j++)
319  for (size_type i = 0; i < res.rows(); i++)
320  res(i, j).vi_ = baseVari->_variRefC[pos++];
321 
322  return res;
323  }
324 
325  template <int R1, int C1, int R2, int C2>
326  inline
327  Eigen::Matrix<var, R1, C2>
328  mdivide_left(const Eigen::Matrix<double, R1, C1> &A,
329  const Eigen::Matrix<var, R2, C2> &b) {
330  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
331 
332  stan::math::check_square("mdivide_left", "A", A);
333  stan::math::check_multiplicable("mdivide_left",
334  "A", A,
335  "b", b);
336 
337  // NOTE: this is not a memory leak, this vari is used in the
338  // expression graph to evaluate the adjoint, but is not needed
339  // for the returned matrix. Memory will be cleaned up with the
340  // arena allocator.
341  mdivide_left_dv_vari<R1, C1, R2, C2> *baseVari
342  = new mdivide_left_dv_vari<R1, C1, R2, C2>(A, b);
343 
344  size_t pos = 0;
345  for (size_type j = 0; j < res.cols(); j++)
346  for (size_type i = 0; i < res.rows(); i++)
347  res(i, j).vi_ = baseVari->_variRefC[pos++];
348 
349  return res;
350  }
351 
352  }
353 }
354 #endif
Eigen::Matrix< fvar< T >, R1, C2 > mdivide_left(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
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
double * C_
vari ** _variRefB
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
int M_
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.
vari ** _variRefA
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
double * A_
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:12
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.
AutodiffStackStorage< vari, chainable_alloc > ChainableStack
int N_
vari ** _variRefC

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