Stan Math Library  2.9.0
reverse mode automatic differentiation
mdivide_left_spd.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_SPD_HPP
2 #define STAN_MATH_REV_MAT_FUN_MDIVIDE_LEFT_SPD_HPP
3 
4 #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_spd_alloc : public chainable_alloc {
18  public:
19  virtual ~mdivide_left_spd_alloc() {}
20 
21  Eigen::LLT< Eigen::Matrix<double, R1, C1> > _llt;
22  Eigen::Matrix<double, R2, C2> C_;
23  };
24 
25  template <int R1, int C1, int R2, int C2>
26  class mdivide_left_spd_vv_vari : public vari {
27  public:
28  int M_; // A.rows() = A.cols() = B.rows()
29  int N_; // B.cols()
30  vari** _variRefA;
31  vari** _variRefB;
32  vari** _variRefC;
33  mdivide_left_spd_alloc<R1, C1, R2, C2> *_alloc;
34 
35  mdivide_left_spd_vv_vari(const Eigen::Matrix<var, R1, C1> &A,
36  const Eigen::Matrix<var, R2, C2> &B)
37  : vari(0.0),
38  M_(A.rows()),
39  N_(B.cols()),
40  _variRefA(reinterpret_cast<vari**>
41  (stan::math::ChainableStack::memalloc_
42  .alloc(sizeof(vari*) * A.rows() * A.cols()))),
43  _variRefB(reinterpret_cast<vari**>
44  (stan::math::ChainableStack::memalloc_
45  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
46  _variRefC(reinterpret_cast<vari**>
47  (stan::math::ChainableStack::memalloc_
48  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
49  _alloc(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
50  using Eigen::Matrix;
51  using Eigen::Map;
52 
53  Matrix<double, R1, C1> Ad(A.rows(), A.cols());
54 
55  size_t pos = 0;
56  for (size_type j = 0; j < M_; j++) {
57  for (size_type i = 0; i < M_; i++) {
58  _variRefA[pos] = A(i, j).vi_;
59  Ad(i, j) = A(i, j).val();
60  pos++;
61  }
62  }
63 
64  pos = 0;
65  _alloc->C_.resize(M_, N_);
66  for (size_type j = 0; j < N_; j++) {
67  for (size_type i = 0; i < M_; i++) {
68  _variRefB[pos] = B(i, j).vi_;
69  _alloc->C_(i, j) = B(i, j).val();
70  pos++;
71  }
72  }
73 
74  _alloc->_llt = Ad.llt();
75  _alloc->_llt.solveInPlace(_alloc->C_);
76 
77  pos = 0;
78  for (size_type j = 0; j < N_; j++) {
79  for (size_type i = 0; i < M_; i++) {
80  _variRefC[pos] = new vari(_alloc->C_(i, j), false);
81  pos++;
82  }
83  }
84  }
85 
86  virtual void chain() {
87  using Eigen::Matrix;
88  using Eigen::Map;
89  Eigen::Matrix<double, R1, C1> adjA(M_, M_);
90  Eigen::Matrix<double, R2, C2> adjB(M_, N_);
91 
92  size_t pos = 0;
93  for (size_type j = 0; j < N_; j++)
94  for (size_type i = 0; i < M_; i++)
95  adjB(i, j) = _variRefC[pos++]->adj_;
96 
97  _alloc->_llt.solveInPlace(adjB);
98  adjA.noalias() = -adjB * _alloc->C_.transpose();
99 
100  pos = 0;
101  for (size_type j = 0; j < M_; j++)
102  for (size_type i = 0; i < M_; i++)
103  _variRefA[pos++]->adj_ += adjA(i, j);
104 
105  pos = 0;
106  for (size_type j = 0; j < N_; j++)
107  for (size_type i = 0; i < M_; 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_spd_dv_vari : public vari {
114  public:
115  int M_; // A.rows() = A.cols() = B.rows()
116  int N_; // B.cols()
117  vari** _variRefB;
118  vari** _variRefC;
119  mdivide_left_spd_alloc<R1, C1, R2, C2> *_alloc;
120 
121  mdivide_left_spd_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  _variRefB(reinterpret_cast<vari**>
127  (stan::math::ChainableStack::memalloc_
128  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
129  _variRefC(reinterpret_cast<vari**>
130  (stan::math::ChainableStack::memalloc_
131  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
132  _alloc(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
133  using Eigen::Matrix;
134  using Eigen::Map;
135 
136  size_t pos = 0;
137  _alloc->C_.resize(M_, N_);
138  for (size_type j = 0; j < N_; j++) {
139  for (size_type i = 0; i < M_; i++) {
140  _variRefB[pos] = B(i, j).vi_;
141  _alloc->C_(i, j) = B(i, j).val();
142  pos++;
143  }
144  }
145 
146  _alloc->_llt = A.llt();
147  _alloc->_llt.solveInPlace(_alloc->C_);
148 
149  pos = 0;
150  for (size_type j = 0; j < N_; j++) {
151  for (size_type i = 0; i < M_; i++) {
152  _variRefC[pos] = new vari(_alloc->C_(i, j), false);
153  pos++;
154  }
155  }
156  }
157 
158  virtual void chain() {
159  using Eigen::Matrix;
160  using Eigen::Map;
161  Eigen::Matrix<double, R2, C2> adjB(M_, N_);
162 
163  size_t pos = 0;
164  for (size_type j = 0; j < adjB.cols(); j++)
165  for (size_type i = 0; i < adjB.rows(); i++)
166  adjB(i, j) = _variRefC[pos++]->adj_;
167 
168  _alloc->_llt.solveInPlace(adjB);
169 
170  pos = 0;
171  for (size_type j = 0; j < adjB.cols(); j++)
172  for (size_type i = 0; i < adjB.rows(); i++)
173  _variRefB[pos++]->adj_ += adjB(i, j);
174  }
175  };
176 
177  template <int R1, int C1, int R2, int C2>
178  class mdivide_left_spd_vd_vari : public vari {
179  public:
180  int M_; // A.rows() = A.cols() = B.rows()
181  int N_; // B.cols()
182  vari** _variRefA;
183  vari** _variRefC;
184  mdivide_left_spd_alloc<R1, C1, R2, C2> *_alloc;
185 
186  mdivide_left_spd_vd_vari(const Eigen::Matrix<var, R1, C1> &A,
187  const Eigen::Matrix<double, R2, C2> &B)
188  : vari(0.0),
189  M_(A.rows()),
190  N_(B.cols()),
191  _variRefA(reinterpret_cast<vari**>
192  (stan::math::ChainableStack::memalloc_
193  .alloc(sizeof(vari*) * A.rows() * A.cols()))),
194  _variRefC(reinterpret_cast<vari**>
195  (stan::math::ChainableStack::memalloc_
196  .alloc(sizeof(vari*) * B.rows() * B.cols()))),
197  _alloc(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
198  using Eigen::Matrix;
199  using Eigen::Map;
200 
201  Matrix<double, R1, C1> Ad(A.rows(), A.cols());
202 
203  size_t pos = 0;
204  for (size_type j = 0; j < M_; j++) {
205  for (size_type i = 0; i < M_; i++) {
206  _variRefA[pos] = A(i, j).vi_;
207  Ad(i, j) = A(i, j).val();
208  pos++;
209  }
210  }
211 
212  _alloc->_llt = Ad.llt();
213  _alloc->C_ = _alloc->_llt.solve(B);
214 
215  pos = 0;
216  for (size_type j = 0; j < N_; j++) {
217  for (size_type i = 0; i < M_; i++) {
218  _variRefC[pos] = new vari(_alloc->C_(i, j), false);
219  pos++;
220  }
221  }
222  }
223 
224  virtual void chain() {
225  using Eigen::Matrix;
226  using Eigen::Map;
227  Eigen::Matrix<double, R1, C1> adjA(M_, M_);
228  Eigen::Matrix<double, R1, C2> adjC(M_, N_);
229 
230  size_t pos = 0;
231  for (size_type j = 0; j < adjC.cols(); j++)
232  for (size_type i = 0; i < adjC.rows(); i++)
233  adjC(i, j) = _variRefC[pos++]->adj_;
234 
235  adjA = -_alloc->_llt.solve(adjC*_alloc->C_.transpose());
236 
237  pos = 0;
238  for (size_type j = 0; j < adjA.cols(); j++)
239  for (size_type i = 0; i < adjA.rows(); i++)
240  _variRefA[pos++]->adj_ += adjA(i, j);
241  }
242  };
243  }
244 
245  template <int R1, int C1, int R2, int C2>
246  inline
247  Eigen::Matrix<var, R1, C2>
248  mdivide_left_spd(const Eigen::Matrix<var, R1, C1> &A,
249  const Eigen::Matrix<var, R2, C2> &b) {
250  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
251 
252  stan::math::check_square("mdivide_left_spd", "A", A);
253  stan::math::check_multiplicable("mdivide_left_spd",
254  "A", A,
255  "b", b);
256 
257  // NOTE: this is not a memory leak, this vari is used in the
258  // expression graph to evaluate the adjoint, but is not needed
259  // for the returned matrix. Memory will be cleaned up with the
260  // arena allocator.
261  mdivide_left_spd_vv_vari<R1, C1, R2, C2> *baseVari
262  = new mdivide_left_spd_vv_vari<R1, C1, R2, C2>(A, b);
263 
264  size_t pos = 0;
265  for (size_type j = 0; j < res.cols(); j++)
266  for (size_type i = 0; i < res.rows(); i++)
267  res(i, j).vi_ = baseVari->_variRefC[pos++];
268 
269  return res;
270  }
271 
272  template <int R1, int C1, int R2, int C2>
273  inline
274  Eigen::Matrix<var, R1, C2>
275  mdivide_left_spd(const Eigen::Matrix<var, R1, C1> &A,
276  const Eigen::Matrix<double, R2, C2> &b) {
277  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
278 
279  stan::math::check_square("mdivide_left_spd", "A", A);
280  stan::math::check_multiplicable("mdivide_left_spd",
281  "A", A,
282  "b", b);
283 
284  // NOTE: this is not a memory leak, this vari is used in the
285  // expression graph to evaluate the adjoint, but is not needed
286  // for the returned matrix. Memory will be cleaned up with the
287  // arena allocator.
288  mdivide_left_spd_vd_vari<R1, C1, R2, C2> *baseVari
289  = new mdivide_left_spd_vd_vari<R1, C1, R2, C2>(A, b);
290 
291  size_t pos = 0;
292  for (size_type j = 0; j < res.cols(); j++)
293  for (size_type i = 0; i < res.rows(); i++)
294  res(i, j).vi_ = baseVari->_variRefC[pos++];
295 
296  return res;
297  }
298 
299  template <int R1, int C1, int R2, int C2>
300  inline
301  Eigen::Matrix<var, R1, C2>
302  mdivide_left_spd(const Eigen::Matrix<double, R1, C1> &A,
303  const Eigen::Matrix<var, R2, C2> &b) {
304  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
305 
306  stan::math::check_square("mdivide_left_spd", "A", A);
307  stan::math::check_multiplicable("mdivide_left_spd",
308  "A", A,
309  "b", b);
310 
311  // NOTE: this is not a memory leak, this vari is used in the
312  // expression graph to evaluate the adjoint, but is not needed
313  // for the returned matrix. Memory will be cleaned up with the
314  // arena allocator.
315  mdivide_left_spd_dv_vari<R1, C1, R2, C2> *baseVari
316  = new mdivide_left_spd_dv_vari<R1, C1, R2, C2>(A, b);
317 
318  size_t pos = 0;
319  for (size_type j = 0; j < res.cols(); j++)
320  for (size_type i = 0; i < res.rows(); i++)
321  res(i, j).vi_ = baseVari->_variRefC[pos++];
322 
323  return res;
324  }
325 
326  }
327 }
328 #endif
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
int N_
vari ** _variRefB
Eigen::LLT< Eigen::Matrix< double, R1, C1 > > _llt
Eigen::Matrix< typename boost::math::tools::promote_args< T1, T2 >::type, R1, C2 > mdivide_left_spd(const Eigen::Matrix< T1, R1, C1 > &A, const Eigen::Matrix< T2, R2, C2 > &b)
Returns the solution of the system Ax=b where A is symmetric positive definite.
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
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.
int M_
mdivide_left_spd_alloc< R1, C1, R2, C2 > * _alloc
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
vari ** _variRefC
vari ** _variRefA
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
Eigen::Matrix< double, R2, C2 > C_

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