Stan Math Library  2.15.0
reverse mode automatic differentiation
multiply.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_MULTIPLY_HPP
2 #define STAN_MATH_REV_MAT_FUN_MULTIPLY_HPP
3 
9 #include <stan/math/rev/core.hpp>
13 #include <boost/math/tools/promotion.hpp>
14 #include <boost/utility/enable_if.hpp>
15 #include <boost/type_traits.hpp>
16 
17 namespace stan {
18  namespace math {
19 
37  template <typename TA, int RA, int CA, typename TB, int CB>
38  class multiply_mat_vari : public vari {
39  public:
41  double* Ad_;
42  double* Bd_;
46 
63  multiply_mat_vari(const Eigen::Matrix<TA, RA, CA>& A,
64  const Eigen::Matrix<TB, CA, CB>& B)
65  : vari(0.0),
66  A_rows_(A.rows()), A_cols_(A.cols()),
67  B_cols_(B.cols()), A_size_(A.size()), B_size_(B.size()),
68  Ad_(ChainableStack::memalloc_.alloc_array<double>(A_size_)),
69  Bd_(ChainableStack::memalloc_.alloc_array<double>(B_size_)),
70  variRefA_(ChainableStack::memalloc_.alloc_array<vari*>(A_size_)),
71  variRefB_(ChainableStack::memalloc_.alloc_array<vari*>(B_size_)),
72  variRefAB_(ChainableStack::memalloc_.alloc_array<vari*>(A_rows_
73  * B_cols_)) {
74  using Eigen::Map;
75  using Eigen::MatrixXd;
76  for (size_type i = 0; i < A.size(); ++i) {
77  variRefA_[i] = A.coeffRef(i).vi_;
78  Ad_[i] = A.coeffRef(i).val();
79  }
80  for (size_type i = 0; i < B.size(); ++i) {
81  variRefB_[i] = B.coeffRef(i).vi_;
82  Bd_[i] = B.coeffRef(i).val();
83  }
84  MatrixXd AB
85  = Map<MatrixXd>(Ad_, A_rows_, A_cols_)
86  * Map<MatrixXd>(Bd_, A_cols_, B_cols_);
87  for (size_type i = 0; i < AB.size(); ++i)
88  variRefAB_[i] = new vari(AB.coeffRef(i), false);
89  }
90 
91  virtual void chain() {
92  using Eigen::MatrixXd;
93  using Eigen::Map;
94  MatrixXd adjAB(A_rows_, B_cols_);
95  MatrixXd adjA(A_rows_, A_cols_);
96  MatrixXd adjB(A_cols_, B_cols_);
97 
98  for (size_type i = 0; i < adjAB.size(); ++i)
99  adjAB(i) = variRefAB_[i]->adj_;
100  adjA = adjAB
101  * Map<MatrixXd>(Bd_, A_cols_, B_cols_).transpose();
102  adjB = Map<MatrixXd>(Ad_, A_rows_, A_cols_).transpose()
103  * adjAB;
104  for (size_type i = 0; i < A_size_; ++i)
105  variRefA_[i]->adj_ += adjA(i);
106  for (size_type i = 0; i < B_size_; ++i)
107  variRefB_[i]->adj_ += adjB(i);
108  }
109  };
110 
126  template <typename TA, int CA, typename TB>
127  class multiply_mat_vari<TA, 1, CA, TB, 1> : public vari {
128  public:
129  int size_;
130  double* Ad_;
131  double* Bd_;
135 
152  multiply_mat_vari(const Eigen::Matrix<TA, 1, CA>& A,
153  const Eigen::Matrix<TB, CA, 1>& B)
154  : vari(0.0),
155  size_(A.cols()),
156  Ad_(ChainableStack::memalloc_.alloc_array<double>(size_)),
157  Bd_(ChainableStack::memalloc_.alloc_array<double>(size_)),
158  variRefA_(ChainableStack::memalloc_.alloc_array<vari*>(size_)),
159  variRefB_(ChainableStack::memalloc_.alloc_array<vari*>(size_)) {
160  using Eigen::Map;
161  using Eigen::VectorXd;
162  using Eigen::RowVectorXd;
163  for (size_type i = 0; i < size_; ++i) {
164  variRefA_[i] = A.coeffRef(i).vi_;
165  Ad_[i] = A.coeffRef(i).val();
166  }
167  for (size_type i = 0; i < size_; ++i) {
168  variRefB_[i] = B.coeffRef(i).vi_;
169  Bd_[i] = B.coeffRef(i).val();
170  }
171  double AB = Map<RowVectorXd>(Ad_, 1, size_)
172  * Map<VectorXd>(Bd_, size_, 1);
173  variRefAB_ = new vari(AB, false);
174  }
175 
176  virtual void chain() {
177  using Eigen::VectorXd;
178  using Eigen::RowVectorXd;
179  using Eigen::Map;
180  double adjAB;
181  RowVectorXd adjA(size_);
182  VectorXd adjB(size_);
183 
184  adjAB = variRefAB_->adj_;
185  adjA = adjAB
186  * Map<VectorXd>(Bd_, size_, 1).transpose();
187  adjB = Map<RowVectorXd>(Ad_, 1, size_).transpose()
188  * adjAB;
189  for (size_type i = 0; i < size_; ++i)
190  variRefA_[i]->adj_ += adjA(i);
191  for (size_type i = 0; i < size_; ++i)
192  variRefB_[i]->adj_ += adjB(i);
193  }
194  };
195 
212  template <int RA, int CA, typename TB, int CB>
213  class multiply_mat_vari<double, RA, CA, TB, CB> : public vari {
214  public:
216  double* Ad_;
217  double* Bd_;
220 
237  multiply_mat_vari(const Eigen::Matrix<double, RA, CA>& A,
238  const Eigen::Matrix<TB, CA, CB>& B)
239  : vari(0.0),
240  A_rows_(A.rows()), A_cols_(A.cols()),
241  B_cols_(B.cols()), A_size_(A.size()), B_size_(B.size()),
242  Ad_(ChainableStack::memalloc_.alloc_array<double>(A_size_)),
243  Bd_(ChainableStack::memalloc_.alloc_array<double>(B_size_)),
244  variRefB_(ChainableStack::memalloc_.alloc_array<vari*>(B_size_)),
245  variRefAB_(ChainableStack::memalloc_.alloc_array<vari*>(A_rows_
246  * B_cols_)) {
247  using Eigen::MatrixXd;
248  using Eigen::Map;
249  for (size_type i = 0; i < A.size(); ++i)
250  Ad_[i] = A.coeffRef(i);
251  for (size_type i = 0; i < B.size(); ++i) {
252  variRefB_[i] = B.coeffRef(i).vi_;
253  Bd_[i] = B.coeffRef(i).val();
254  }
255  MatrixXd AB
256  = Map<MatrixXd>(Ad_, A_rows_, A_cols_)
257  * Map<MatrixXd>(Bd_, A_cols_, B_cols_);
258  for (size_type i = 0; i < AB.size(); ++i)
259  variRefAB_[i] = new vari(AB.coeffRef(i), false);
260  }
261 
262  virtual void chain() {
263  using Eigen::MatrixXd;
264  using Eigen::Map;
265  MatrixXd adjAB(A_rows_, B_cols_);
266  MatrixXd adjB(A_cols_, B_cols_);
267 
268  for (size_type i = 0; i < adjAB.size(); ++i)
269  adjAB(i) = variRefAB_[i]->adj_;
270  adjB = Map<MatrixXd>(Ad_, A_rows_, A_cols_).transpose()
271  * adjAB;
272  for (size_type i = 0; i < B_size_; ++i)
273  variRefB_[i]->adj_ += adjB(i);
274  }
275  };
276 
292  template <int CA, typename TB>
293  class multiply_mat_vari<double, 1, CA, TB, 1> : public vari {
294  public:
295  int size_;
296  double* Ad_;
297  double* Bd_;
300 
317  multiply_mat_vari(const Eigen::Matrix<double, 1, CA>& A,
318  const Eigen::Matrix<TB, CA, 1>& B)
319  : vari(0.0),
320  size_(A.cols()),
321  Ad_(ChainableStack::memalloc_.alloc_array<double>(size_)),
322  Bd_(ChainableStack::memalloc_.alloc_array<double>(size_)),
323  variRefB_(ChainableStack::memalloc_.alloc_array<vari*>(size_)) {
324  using Eigen::Map;
325  using Eigen::VectorXd;
326  using Eigen::RowVectorXd;
327  for (size_type i = 0; i < size_; ++i)
328  Ad_[i] = A.coeffRef(i);
329  for (size_type i = 0; i < size_; ++i) {
330  variRefB_[i] = B.coeffRef(i).vi_;
331  Bd_[i] = B.coeffRef(i).val();
332  }
333  double AB
334  = Eigen::Map<RowVectorXd>(Ad_, 1, size_)
335  * Eigen::Map<VectorXd>(Bd_, size_, 1);
336  variRefAB_ = new vari(AB, false);
337  }
338 
339  virtual void chain() {
340  using Eigen::RowVectorXd;
341  using Eigen::VectorXd;
342  using Eigen::Map;
343  double adjAB;
344  VectorXd adjB(size_);
345 
346  adjAB = variRefAB_->adj_;
347  adjB = Map<RowVectorXd>(Ad_, 1, size_).transpose()
348  * adjAB;
349  for (size_type i = 0; i < size_; ++i)
350  variRefB_[i]->adj_ += adjB(i);
351  }
352  };
353 
370  template <typename TA, int RA, int CA, int CB>
371  class multiply_mat_vari<TA, RA, CA, double, CB> : public vari {
372  public:
374  double* Ad_;
375  double* Bd_;
378 
395  multiply_mat_vari(const Eigen::Matrix<TA, RA, CA>& A,
396  const Eigen::Matrix<double, CA, CB>& B)
397  : vari(0.0),
398  A_rows_(A.rows()), A_cols_(A.cols()),
399  B_cols_(B.cols()), A_size_(A.size()), B_size_(B.size()),
400  Ad_(ChainableStack::memalloc_.alloc_array<double>(A_size_)),
401  Bd_(ChainableStack::memalloc_.alloc_array<double>(B_size_)),
402  variRefA_(ChainableStack::memalloc_.alloc_array<vari*>(A_size_)),
403  variRefAB_(ChainableStack::memalloc_.alloc_array<vari*>(A_rows_
404  * B_cols_)) {
405  using Eigen::Map;
406  using Eigen::MatrixXd;
407  for (size_type i = 0; i < A_size_; ++i) {
408  variRefA_[i] = A.coeffRef(i).vi_;
409  Ad_[i] = A.coeffRef(i).val();
410  }
411  for (size_type i = 0; i < B_size_; ++i) {
412  Bd_[i] = B.coeffRef(i);
413  }
414  MatrixXd AB
415  = Map<MatrixXd>(Ad_, A_rows_, A_cols_)
416  * Map<MatrixXd>(Bd_, A_cols_, B_cols_);
417  for (size_type i = 0; i < AB.size(); ++i)
418  variRefAB_[i] = new vari(AB.coeffRef(i), false);
419  }
420 
421  virtual void chain() {
422  using Eigen::MatrixXd;
423  using Eigen::Map;
424  MatrixXd adjAB(A_rows_, B_cols_);
425  MatrixXd adjA(A_rows_, A_cols_);
426 
427  for (size_type i = 0; i < adjAB.size(); ++i)
428  adjAB(i) = variRefAB_[i]->adj_;
429  adjA = adjAB * Map<MatrixXd>(Bd_, A_cols_, B_cols_).transpose();
430  for (size_type i = 0; i < A_size_; ++i)
431  variRefA_[i]->adj_ += adjA(i);
432  }
433  };
434 
453  template <typename TA, int CA>
454  class multiply_mat_vari<TA, 1, CA, double, 1> : public vari {
455  public:
456  int size_;
457  double* Ad_;
458  double* Bd_;
461 
478  multiply_mat_vari(const Eigen::Matrix<TA, 1, CA>& A,
479  const Eigen::Matrix<double, CA, 1>& B)
480  : vari(0.0),
481  size_(A.cols()),
482  Ad_(ChainableStack::memalloc_.alloc_array<double>(size_)),
483  Bd_(ChainableStack::memalloc_.alloc_array<double>(size_)),
484  variRefA_(ChainableStack::memalloc_.alloc_array<vari*>(size_)) {
485  using Eigen::Map;
486  using Eigen::VectorXd;
487  using Eigen::RowVectorXd;
488  for (size_type i = 0; i < size_; ++i) {
489  variRefA_[i] = A.coeffRef(i).vi_;
490  Ad_[i] = A.coeffRef(i).val();
491  }
492  for (size_type i = 0; i < size_; ++i)
493  Bd_[i] = B.coeffRef(i);
494  double AB
495  = Map<RowVectorXd>(Ad_, 1, size_)
496  * Map<VectorXd>(Bd_, size_, 1);
497  variRefAB_ = new vari(AB, false);
498  }
499 
500  virtual void chain() {
501  using Eigen::Map;
502  using Eigen::VectorXd;
503  using Eigen::RowVectorXd;
504  double adjAB;
505  RowVectorXd adjA(size_);
506 
507  adjAB = variRefAB_->adj_;
508  adjA = adjAB
509  * Map<VectorXd>(Bd_, size_, 1).transpose();
510  for (size_type i = 0; i < size_; ++i)
511  variRefA_[i]->adj_ += adjA(i);
512  }
513  };
514 
523  template <typename T1, typename T2>
524  inline typename
525  boost::enable_if_c<
526  (boost::is_scalar<T1>::value || boost::is_same<T1, var>::value)
527  && (boost::is_scalar<T2>::value || boost::is_same<T2, var>::value),
528  typename boost::math::tools::promote_args<T1, T2>::type>::type
529  multiply(const T1& v, const T2& c) {
530  return v * c;
531  }
532 
543  template<typename T1, typename T2, int R2, int C2>
544  inline Eigen::Matrix<var, R2, C2>
545  multiply(const T1& c, const Eigen::Matrix<T2, R2, C2>& m) {
546  // TODO(trangucci) pull out to eliminate overpromotion of one side
547  // move to matrix.hpp w. promotion?
548  return to_var(m) * to_var(c);
549  }
550 
561  template<typename T1, int R1, int C1, typename T2>
562  inline Eigen::Matrix<var, R1, C1>
563  multiply(const Eigen::Matrix<T1, R1, C1>& m, const T2& c) {
564  // TODO(trangucci) pull out to eliminate overpromotion of one side
565  // move to matrix.hpp w. promotion?
566  return to_var(m) * to_var(c);
567  }
568 
581  template <typename TA, int RA, int CA, typename TB, int CB>
582  inline typename
583  boost::enable_if_c<boost::is_same<TA, var>::value
584  || boost::is_same<TB, var>::value,
585  Eigen::Matrix<var, RA, CB> >::type
586  multiply(const Eigen::Matrix<TA, RA, CA> &A,
587  const Eigen::Matrix<TB, CA, CB> &B) {
588  check_multiplicable("multiply", "A", A, "B", B);
589  check_not_nan("multiply", "A", A);
590  check_not_nan("multiply", "B", B);
591 
592  // Memory managed with the arena allocator.
595  Eigen::Matrix<var, RA, CB> AB_v(A.rows(), B.cols());
596  for (size_type i = 0; i < AB_v.size(); ++i) {
597  AB_v.coeffRef(i).vi_ = baseVari->variRefAB_[i];
598  }
599  return AB_v;
600  }
601 
612  template <typename TA, int CA, typename TB>
613  inline typename
614  boost::enable_if_c<boost::is_same<TA, var>::value
615  || boost::is_same<TB, var>::value, var>::type
616  multiply(const Eigen::Matrix<TA, 1, CA> &A,
617  const Eigen::Matrix<TB, CA, 1> &B) {
618  check_multiplicable("multiply", "A", A, "B", B);
619  check_not_nan("multiply", "A", A);
620  check_not_nan("multiply", "B", B);
621 
622  // Memory managed with the arena allocator.
625  var AB_v;
626  AB_v.vi_ = baseVari->variRefAB_;
627  return AB_v;
628  }
629  }
630 }
631 #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
multiply_mat_vari(const Eigen::Matrix< TA, RA, CA > &A, const Eigen::Matrix< double, CA, CB > &B)
Constructor for multiply_mat_vari.
Definition: multiply.hpp:395
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:421
multiply_mat_vari(const Eigen::Matrix< double, 1, CA > &A, const Eigen::Matrix< TB, CA, 1 > &B)
Constructor for multiply_mat_vari.
Definition: multiply.hpp:317
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:91
This is a subclass of the vari class for matrix multiplication A * B where A is N by M and B is M by ...
Definition: multiply.hpp:38
The variable implementation base class.
Definition: vari.hpp:30
Eigen::Matrix< fvar< T >, R1, C1 > multiply(const Eigen::Matrix< fvar< T >, R1, C1 > &m, const fvar< T > &c)
Definition: multiply.hpp:20
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:30
multiply_mat_vari(const Eigen::Matrix< TA, RA, CA > &A, const Eigen::Matrix< TB, CA, CB > &B)
Constructor for multiply_mat_vari.
Definition: multiply.hpp:63
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:339
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:262
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
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
size_t size_
Definition: dot_self.hpp:18
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:176
multiply_mat_vari(const Eigen::Matrix< TA, 1, CA > &A, const Eigen::Matrix< TB, CA, 1 > &B)
Constructor for multiply_mat_vari.
Definition: multiply.hpp:152
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:500
multiply_mat_vari(const Eigen::Matrix< double, RA, CA > &A, const Eigen::Matrix< TB, CA, CB > &B)
Constructor for multiply_mat_vari.
Definition: multiply.hpp:237
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.
vari(double x)
Construct a variable implementation from a value.
Definition: vari.hpp:58
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
multiply_mat_vari(const Eigen::Matrix< TA, 1, CA > &A, const Eigen::Matrix< double, CA, 1 > &B)
Constructor for multiply_mat_vari.
Definition: multiply.hpp:478
This is a subclass of the vari class for matrix multiplication A * B where A is 1 by M and B is M by ...
Definition: multiply.hpp:127
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:12
std::vector< var > to_var(const std::vector< double > &v)
Converts argument to an automatic differentiation variable.
Definition: to_var.hpp:20

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