Stan Math Library  2.20.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/meta.hpp>
5 #include <stan/math/rev/core.hpp>
10 #include <vector>
11 
12 namespace stan {
13 namespace math {
14 
15 namespace internal {
16 template <int R1, int C1, int R2, int C2>
18  public:
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>
27  public:
28  int M_; // A.rows() = A.cols() = B.rows()
29  int N_; // B.cols()
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  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * A.rows()
42  * A.cols()))),
43  variRefB_(reinterpret_cast<vari **>(
44  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
45  * B.cols()))),
46  variRefC_(reinterpret_cast<vari **>(
47  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
48  * B.cols()))),
49  alloc_(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
50  using Eigen::Map;
51  using Eigen::Matrix;
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::Map;
88  using Eigen::Matrix;
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>
114  public:
115  int M_; // A.rows() = A.cols() = B.rows()
116  int N_; // B.cols()
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  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  alloc_(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
133  using Eigen::Map;
134  using Eigen::Matrix;
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::Map;
160  using Eigen::Matrix;
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>
179  public:
180  int M_; // A.rows() = A.cols() = B.rows()
181  int N_; // B.cols()
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  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * A.rows()
193  * A.cols()))),
194  variRefC_(reinterpret_cast<vari **>(
195  ChainableStack::instance_->memalloc_.alloc(sizeof(vari *) * B.rows()
196  * B.cols()))),
197  alloc_(new mdivide_left_spd_alloc<R1, C1, R2, C2>()) {
198  using Eigen::Map;
199  using Eigen::Matrix;
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::Map;
226  using Eigen::Matrix;
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 } // namespace internal
244 
245 template <int R1, int C1, int R2, int C2>
246 inline Eigen::Matrix<var, R1, C2> mdivide_left_spd(
247  const Eigen::Matrix<var, R1, C1> &A, const Eigen::Matrix<var, R2, C2> &b) {
248  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
249 
250  check_square("mdivide_left_spd", "A", A);
251  check_multiplicable("mdivide_left_spd", "A", A, "b", b);
252 
253  // NOTE: this is not a memory leak, this vari is used in the
254  // expression graph to evaluate the adjoint, but is not needed
255  // for the returned matrix. Memory will be cleaned up with the
256  // arena allocator.
259 
260  size_t pos = 0;
261  for (size_type j = 0; j < res.cols(); j++)
262  for (size_type i = 0; i < res.rows(); i++)
263  res(i, j).vi_ = baseVari->variRefC_[pos++];
264 
265  return res;
266 }
267 
268 template <int R1, int C1, int R2, int C2>
269 inline Eigen::Matrix<var, R1, C2> mdivide_left_spd(
270  const Eigen::Matrix<var, R1, C1> &A,
271  const Eigen::Matrix<double, R2, C2> &b) {
272  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
273 
274  check_square("mdivide_left_spd", "A", A);
275  check_multiplicable("mdivide_left_spd", "A", A, "b", b);
276 
277  // NOTE: this is not a memory leak, this vari is used in the
278  // expression graph to evaluate the adjoint, but is not needed
279  // for the returned matrix. Memory will be cleaned up with the
280  // arena allocator.
283 
284  size_t pos = 0;
285  for (size_type j = 0; j < res.cols(); j++)
286  for (size_type i = 0; i < res.rows(); i++)
287  res(i, j).vi_ = baseVari->variRefC_[pos++];
288 
289  return res;
290 }
291 
292 template <int R1, int C1, int R2, int C2>
293 inline Eigen::Matrix<var, R1, C2> mdivide_left_spd(
294  const Eigen::Matrix<double, R1, C1> &A,
295  const Eigen::Matrix<var, R2, C2> &b) {
296  Eigen::Matrix<var, R1, C2> res(b.rows(), b.cols());
297 
298  check_square("mdivide_left_spd", "A", A);
299  check_multiplicable("mdivide_left_spd", "A", A, "b", b);
300 
301  // NOTE: this is not a memory leak, this vari is used in the
302  // expression graph to evaluate the adjoint, but is not needed
303  // for the returned matrix. Memory will be cleaned up with the
304  // arena allocator.
307 
308  size_t pos = 0;
309  for (size_type j = 0; j < res.cols(); j++)
310  for (size_type i = 0; i < res.rows(); i++)
311  res(i, j).vi_ = baseVari->variRefC_[pos++];
312 
313  return res;
314 }
315 
316 } // namespace math
317 } // namespace stan
318 #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
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
mdivide_left_spd_alloc< R1, C1, R2, C2 > * alloc_
void check_square(const char *function, const char *name, const matrix_cl &y)
Check if the matrix_cl is square.
The variable implementation base class.
Definition: vari.hpp:30
mdivide_left_spd_vd_vari(const Eigen::Matrix< var, R1, C1 > &A, const Eigen::Matrix< double, R2, C2 > &B)
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:11
mdivide_left_spd_vv_vari(const Eigen::Matrix< var, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Eigen::LLT< Eigen::Matrix< double, R1, C1 > > llt_
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
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.
mdivide_left_spd_alloc< R1, C1, R2, C2 > * alloc_
mdivide_left_spd_alloc< R1, C1, R2, C2 > * alloc_
A chainable_alloc is an object which is constructed and destructed normally but the memory lifespan i...
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
This struct always provides access to the autodiff stack using the singleton pattern.
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
mdivide_left_spd_dv_vari(const Eigen::Matrix< double, R1, C1 > &A, const Eigen::Matrix< var, R2, C2 > &B)

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