Stan Math Library  2.20.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 
4 #include <stan/math/rev/meta.hpp>
7 #include <stan/math/rev/core.hpp>
10 #include <boost/math/tools/promotion.hpp>
11 #include <type_traits>
12 
13 namespace stan {
14 namespace math {
15 
33 template <typename Ta, int Ra, int Ca, typename Tb, int Cb>
34 class multiply_mat_vari : public vari {
35  public:
36  int A_rows_;
37  int A_cols_;
38  int B_cols_;
39  int A_size_;
40  int B_size_;
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()),
67  A_cols_(A.cols()),
68  B_cols_(B.cols()),
69  A_size_(A.size()),
70  B_size_(B.size()),
71  Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(A_size_)),
72  Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(B_size_)),
73  variRefA_(
74  ChainableStack::instance_->memalloc_.alloc_array<vari*>(A_size_)),
75  variRefB_(
76  ChainableStack::instance_->memalloc_.alloc_array<vari*>(B_size_)),
77  variRefAB_(ChainableStack::instance_->memalloc_.alloc_array<vari*>(
78  A_rows_ * B_cols_)) {
79  using Eigen::Map;
80  using Eigen::MatrixXd;
81  for (size_type i = 0; i < A.size(); ++i) {
82  variRefA_[i] = A.coeffRef(i).vi_;
83  Ad_[i] = A.coeffRef(i).val();
84  }
85  for (size_type i = 0; i < B.size(); ++i) {
86  variRefB_[i] = B.coeffRef(i).vi_;
87  Bd_[i] = B.coeffRef(i).val();
88  }
89  MatrixXd AB = Map<MatrixXd>(Ad_, A_rows_, A_cols_)
90  * Map<MatrixXd>(Bd_, A_cols_, B_cols_);
91  for (size_type i = 0; i < AB.size(); ++i)
92  variRefAB_[i] = new vari(AB.coeffRef(i), false);
93  }
94 
95  virtual void chain() {
96  using Eigen::Map;
97  using Eigen::MatrixXd;
98  MatrixXd adjAB(A_rows_, B_cols_);
99  MatrixXd adjA(A_rows_, A_cols_);
100  MatrixXd adjB(A_cols_, B_cols_);
101 
102  for (size_type i = 0; i < adjAB.size(); ++i)
103  adjAB(i) = variRefAB_[i]->adj_;
104  adjA = adjAB * Map<MatrixXd>(Bd_, A_cols_, B_cols_).transpose();
105  adjB = Map<MatrixXd>(Ad_, A_rows_, A_cols_).transpose() * adjAB;
106  for (size_type i = 0; i < A_size_; ++i)
107  variRefA_[i]->adj_ += adjA(i);
108  for (size_type i = 0; i < B_size_; ++i)
109  variRefB_[i]->adj_ += adjB(i);
110  }
111 };
112 
128 template <typename Ta, int Ca, typename Tb>
129 class multiply_mat_vari<Ta, 1, Ca, Tb, 1> : public vari {
130  public:
131  int size_;
132  double* Ad_;
133  double* Bd_;
137 
154  multiply_mat_vari(const Eigen::Matrix<Ta, 1, Ca>& A,
155  const Eigen::Matrix<Tb, Ca, 1>& B)
156  : vari(0.0),
157  size_(A.cols()),
158  Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
159  Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
160  variRefA_(
161  ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)),
162  variRefB_(
163  ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
164  using Eigen::Map;
165  using Eigen::RowVectorXd;
166  using Eigen::VectorXd;
167  for (size_type i = 0; i < size_; ++i) {
168  variRefA_[i] = A.coeffRef(i).vi_;
169  Ad_[i] = A.coeffRef(i).val();
170  }
171  for (size_type i = 0; i < size_; ++i) {
172  variRefB_[i] = B.coeffRef(i).vi_;
173  Bd_[i] = B.coeffRef(i).val();
174  }
175  double AB = Map<RowVectorXd>(Ad_, 1, size_) * Map<VectorXd>(Bd_, size_, 1);
176  variRefAB_ = new vari(AB, false);
177  }
178 
179  virtual void chain() {
180  using Eigen::Map;
181  using Eigen::RowVectorXd;
182  using Eigen::VectorXd;
183  double adjAB;
184 
185  adjAB = variRefAB_->adj_;
186  auto adjA = adjAB * Map<VectorXd>(Bd_, size_, 1).transpose();
187  auto adjB = Map<RowVectorXd>(Ad_, 1, size_).transpose() * adjAB;
188  for (size_type i = 0; i < size_; ++i)
189  variRefA_[i]->adj_ += adjA(i);
190  for (size_type i = 0; i < size_; ++i)
191  variRefB_[i]->adj_ += adjB(i);
192  }
193 };
194 
211 template <int Ra, int Ca, typename Tb, int Cb>
212 class multiply_mat_vari<double, Ra, Ca, Tb, Cb> : public vari {
213  public:
214  int A_rows_;
215  int A_cols_;
216  int B_cols_;
217  int A_size_;
218  int B_size_;
219  double* Ad_;
220  double* Bd_;
223 
240  multiply_mat_vari(const Eigen::Matrix<double, Ra, Ca>& A,
241  const Eigen::Matrix<Tb, Ca, Cb>& B)
242  : vari(0.0),
243  A_rows_(A.rows()),
244  A_cols_(A.cols()),
245  B_cols_(B.cols()),
246  A_size_(A.size()),
247  B_size_(B.size()),
248  Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(A_size_)),
249  Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(B_size_)),
250  variRefB_(
251  ChainableStack::instance_->memalloc_.alloc_array<vari*>(B_size_)),
252  variRefAB_(ChainableStack::instance_->memalloc_.alloc_array<vari*>(
253  A_rows_ * B_cols_)) {
254  using Eigen::Map;
255  using Eigen::MatrixXd;
256  for (size_type i = 0; i < A.size(); ++i)
257  Ad_[i] = A.coeffRef(i);
258  for (size_type i = 0; i < B.size(); ++i) {
259  variRefB_[i] = B.coeffRef(i).vi_;
260  Bd_[i] = B.coeffRef(i).val();
261  }
262  MatrixXd AB = Map<MatrixXd>(Ad_, A_rows_, A_cols_)
263  * Map<MatrixXd>(Bd_, A_cols_, B_cols_);
264  for (size_type i = 0; i < AB.size(); ++i)
265  variRefAB_[i] = new vari(AB.coeffRef(i), false);
266  }
267 
268  virtual void chain() {
269  using Eigen::Map;
270  using Eigen::MatrixXd;
271  MatrixXd adjAB(A_rows_, B_cols_);
272  MatrixXd adjB(A_cols_, B_cols_);
273 
274  for (size_type i = 0; i < adjAB.size(); ++i)
275  adjAB(i) = variRefAB_[i]->adj_;
276  adjB = Map<MatrixXd>(Ad_, A_rows_, A_cols_).transpose() * adjAB;
277  for (size_type i = 0; i < B_size_; ++i)
278  variRefB_[i]->adj_ += adjB(i);
279  }
280 };
281 
297 template <int Ca, typename Tb>
298 class multiply_mat_vari<double, 1, Ca, Tb, 1> : public vari {
299  public:
300  int size_;
301  double* Ad_;
302  double* Bd_;
305 
322  multiply_mat_vari(const Eigen::Matrix<double, 1, Ca>& A,
323  const Eigen::Matrix<Tb, Ca, 1>& B)
324  : vari(0.0),
325  size_(A.cols()),
326  Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
327  Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
328  variRefB_(
329  ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
330  using Eigen::Map;
331  using Eigen::RowVectorXd;
332  using Eigen::VectorXd;
333  for (size_type i = 0; i < size_; ++i)
334  Ad_[i] = A.coeffRef(i);
335  for (size_type i = 0; i < size_; ++i) {
336  variRefB_[i] = B.coeffRef(i).vi_;
337  Bd_[i] = B.coeffRef(i).val();
338  }
339  double AB = Eigen::Map<RowVectorXd>(Ad_, 1, size_)
340  * Eigen::Map<VectorXd>(Bd_, size_, 1);
341  variRefAB_ = new vari(AB, false);
342  }
343 
344  virtual void chain() {
345  using Eigen::Map;
346  using Eigen::RowVectorXd;
347  using Eigen::VectorXd;
348  double adjAB;
349 
350  adjAB = variRefAB_->adj_;
351  auto adjB = Map<RowVectorXd>(Ad_, 1, size_).transpose() * adjAB;
352  for (size_type i = 0; i < size_; ++i)
353  variRefB_[i]->adj_ += adjB(i);
354  }
355 };
356 
373 template <typename Ta, int Ra, int Ca, int Cb>
374 class multiply_mat_vari<Ta, Ra, Ca, double, Cb> : public vari {
375  public:
376  int A_rows_;
377  int A_cols_;
378  int B_cols_;
379  int A_size_;
380  int B_size_;
381  double* Ad_;
382  double* Bd_;
385 
402  multiply_mat_vari(const Eigen::Matrix<Ta, Ra, Ca>& A,
403  const Eigen::Matrix<double, Ca, Cb>& B)
404  : vari(0.0),
405  A_rows_(A.rows()),
406  A_cols_(A.cols()),
407  B_cols_(B.cols()),
408  A_size_(A.size()),
409  B_size_(B.size()),
410  Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(A_size_)),
411  Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(B_size_)),
412  variRefA_(
413  ChainableStack::instance_->memalloc_.alloc_array<vari*>(A_size_)),
414  variRefAB_(ChainableStack::instance_->memalloc_.alloc_array<vari*>(
415  A_rows_ * B_cols_)) {
416  using Eigen::Map;
417  using Eigen::MatrixXd;
418  for (size_type i = 0; i < A_size_; ++i) {
419  variRefA_[i] = A.coeffRef(i).vi_;
420  Ad_[i] = A.coeffRef(i).val();
421  }
422  for (size_type i = 0; i < B_size_; ++i) {
423  Bd_[i] = B.coeffRef(i);
424  }
425  MatrixXd AB = Map<MatrixXd>(Ad_, A_rows_, A_cols_)
426  * Map<MatrixXd>(Bd_, A_cols_, B_cols_);
427  for (size_type i = 0; i < AB.size(); ++i)
428  variRefAB_[i] = new vari(AB.coeffRef(i), false);
429  }
430 
431  virtual void chain() {
432  using Eigen::Map;
433  using Eigen::MatrixXd;
434  MatrixXd adjAB(A_rows_, B_cols_);
435  MatrixXd adjA(A_rows_, A_cols_);
436 
437  for (size_type i = 0; i < adjAB.size(); ++i)
438  adjAB(i) = variRefAB_[i]->adj_;
439  adjA = adjAB * Map<MatrixXd>(Bd_, A_cols_, B_cols_).transpose();
440  for (size_type i = 0; i < A_size_; ++i)
441  variRefA_[i]->adj_ += adjA(i);
442  }
443 };
444 
463 template <typename Ta, int Ca>
464 class multiply_mat_vari<Ta, 1, Ca, double, 1> : public vari {
465  public:
466  int size_;
467  double* Ad_;
468  double* Bd_;
471 
488  multiply_mat_vari(const Eigen::Matrix<Ta, 1, Ca>& A,
489  const Eigen::Matrix<double, Ca, 1>& B)
490  : vari(0.0),
491  size_(A.cols()),
492  Ad_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
493  Bd_(ChainableStack::instance_->memalloc_.alloc_array<double>(size_)),
494  variRefA_(
495  ChainableStack::instance_->memalloc_.alloc_array<vari*>(size_)) {
496  using Eigen::Map;
497  using Eigen::RowVectorXd;
498  using Eigen::VectorXd;
499  for (size_type i = 0; i < size_; ++i) {
500  variRefA_[i] = A.coeffRef(i).vi_;
501  Ad_[i] = A.coeffRef(i).val();
502  }
503  for (size_type i = 0; i < size_; ++i)
504  Bd_[i] = B.coeffRef(i);
505  double AB = Map<RowVectorXd>(Ad_, 1, size_) * Map<VectorXd>(Bd_, size_, 1);
506  variRefAB_ = new vari(AB, false);
507  }
508 
509  virtual void chain() {
510  using Eigen::Map;
511  using Eigen::RowVectorXd;
512  using Eigen::VectorXd;
513  double adjAB;
514 
515  adjAB = variRefAB_->adj_;
516  auto adjA = adjAB * Map<VectorXd>(Bd_, size_, 1).transpose();
517  for (size_type i = 0; i < size_; ++i)
518  variRefA_[i]->adj_ += adjA(i);
519  }
520 };
521 
530 template <typename T1, typename T2>
531 inline typename std::enable_if<
532  (boost::is_scalar<T1>::value || std::is_same<T1, var>::value)
533  && (boost::is_scalar<T2>::value || std::is_same<T2, var>::value),
534  typename boost::math::tools::promote_args<T1, T2>::type>::type
535 multiply(const T1& v, const T2& c) {
536  return v * c;
537 }
538 
549 template <typename T1, typename T2, int R2, int C2>
550 inline Eigen::Matrix<var, R2, C2> multiply(const T1& c,
551  const Eigen::Matrix<T2, R2, C2>& m) {
552  // TODO(trangucci) pull out to eliminate overpromotion of one side
553  // move to matrix.hpp w. promotion?
554  return to_var(m) * to_var(c);
555 }
556 
567 template <typename T1, int R1, int C1, typename T2>
568 inline Eigen::Matrix<var, R1, C1> multiply(const Eigen::Matrix<T1, R1, C1>& m,
569  const T2& c) {
570  // TODO(trangucci) pull out to eliminate overpromotion of one side
571  // move to matrix.hpp w. promotion?
572  return to_var(m) * to_var(c);
573 }
574 
587 template <typename Ta, int Ra, int Ca, typename Tb, int Cb>
588 inline typename std::enable_if<std::is_same<Ta, var>::value
589  || std::is_same<Tb, var>::value,
590  Eigen::Matrix<var, Ra, Cb> >::type
591 multiply(const Eigen::Matrix<Ta, Ra, Ca>& A,
592  const Eigen::Matrix<Tb, Ca, Cb>& B) {
593  check_multiplicable("multiply", "A", A, "B", B);
594  check_not_nan("multiply", "A", A);
595  check_not_nan("multiply", "B", B);
596 
597  // Memory managed with the arena allocator.
600  Eigen::Matrix<var, Ra, Cb> AB_v(A.rows(), B.cols());
601  for (size_type i = 0; i < AB_v.size(); ++i) {
602  AB_v.coeffRef(i).vi_ = baseVari->variRefAB_[i];
603  }
604  return AB_v;
605 }
606 
617 template <typename Ta, int Ca, typename Tb>
618 inline typename std::enable_if<
619  std::is_same<Ta, var>::value || std::is_same<Tb, var>::value, var>::type
620 multiply(const Eigen::Matrix<Ta, 1, Ca>& A, const Eigen::Matrix<Tb, Ca, 1>& B) {
621  check_multiplicable("multiply", "A", A, "B", B);
622  check_not_nan("multiply", "A", A);
623  check_not_nan("multiply", "B", B);
624 
625  // Memory managed with the arena allocator.
628  var AB_v;
629  AB_v.vi_ = baseVari->variRefAB_;
630  return AB_v;
631 }
632 } // namespace math
633 } // namespace stan
634 #endif
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:322
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. ...
Definition: multiply.hpp:509
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:34
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:488
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:14
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:154
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:33
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
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:95
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:344
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
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:240
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:129
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:268
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:431
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
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
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.
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:402
matrix_cl transpose(const matrix_cl &src)
Takes the transpose of the matrix on the OpenCL device.
Definition: transpose.hpp:20
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
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
Definition: multiply.hpp:179
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.
std::vector< var > to_var(const std::vector< double > &v)
Converts argument to an automatic differentiation variable.
Definition: to_var.hpp:19

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