Stan Math Library  2.14.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  (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*>
128  .alloc(sizeof(double) * A.rows() * A.cols()))),
129  C_(reinterpret_cast<double*>
131  .alloc(sizeof(double) * B.rows() * B.cols()))),
132  variRefB_(reinterpret_cast<vari**>
134  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
135  variRefC_(reinterpret_cast<vari**>
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*>
210  .alloc(sizeof(double) * A.rows() * A.cols()))),
211  C_(reinterpret_cast<double*>
213  .alloc(sizeof(double) * B.rows() * B.cols()))),
214  variRefA_(reinterpret_cast<vari**>
216  .alloc(sizeof(vari*) * A.rows() * A.cols()))),
217  variRefC_(reinterpret_cast<vari**>
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
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_
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_
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 check_square(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Check if the specified matrix is square.
vari ** variRefA_
vari ** variRefC_
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:12
AutodiffStackStorage< vari, chainable_alloc > ChainableStack
int N_

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