Stan Math Library  2.15.0
reverse mode automatic differentiation
cholesky_decompose.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_CHOLESKY_DECOMPOSE_HPP
2 #define STAN_MATH_REV_MAT_FUN_CHOLESKY_DECOMPOSE_HPP
3 
9 #include <stan/math/rev/core.hpp>
14 #include <algorithm>
15 
16 namespace stan {
17  namespace math {
18 
19  class cholesky_block : public vari {
20  public:
21  int M_;
23  typedef Eigen::Block<Eigen::MatrixXd> Block_;
26 
43  cholesky_block(const Eigen::Matrix<var, -1, -1>& A,
44  const Eigen::Matrix<double, -1, -1>& L_A)
45  : vari(0.0),
46  M_(A.rows()),
47  variRefA_(ChainableStack::memalloc_.alloc_array<vari*>
48  (A.rows() * (A.rows() + 1) / 2)),
49  variRefL_(ChainableStack::memalloc_.alloc_array<vari*>
50  (A.rows() * (A.rows() + 1) / 2)) {
51  size_t pos = 0;
52  block_size_ = std::max((M_ / 8 / 16) * 16, 8);
53  block_size_ = std::min(block_size_, 128);
54  for (size_type j = 0; j < M_; ++j) {
55  for (size_type i = j; i < M_; ++i) {
56  variRefA_[pos] = A.coeffRef(i, j).vi_;
57  variRefL_[pos] = new vari(L_A.coeffRef(i, j), false); ++pos;
58  }
59  }
60  }
61 
68  inline void symbolic_rev(Block_& L,
69  Block_& Lbar) {
70  using Eigen::Lower;
71  using Eigen::Upper;
72  using Eigen::StrictlyUpper;
73  L.transposeInPlace();
74  Lbar = (L * Lbar.triangularView<Lower>()).eval();
75  Lbar.triangularView<StrictlyUpper>()
76  = Lbar.adjoint().triangularView<StrictlyUpper>();
77  L.triangularView<Upper>().solveInPlace(Lbar);
78  L.triangularView<Upper>().solveInPlace(Lbar.transpose());
79  }
80 
87  virtual void chain() {
88  using Eigen::MatrixXd;
89  using Eigen::Lower;
90  using Eigen::Block;
91  using Eigen::Upper;
92  using Eigen::StrictlyUpper;
93  MatrixXd Lbar(M_, M_);
94  MatrixXd L(M_, M_);
95 
96  Lbar.setZero();
97  L.setZero();
98  size_t pos = 0;
99  for (size_type j = 0; j < M_; ++j) {
100  for (size_type i = j; i < M_; ++i) {
101  Lbar.coeffRef(i, j) = variRefL_[pos]->adj_;
102  L.coeffRef(i, j) = variRefL_[pos]->val_;
103  ++pos;
104  }
105  }
106 
107  for (int k = M_; k > 0; k -= block_size_) {
108  int j = std::max(0, k - block_size_);
109  Block_ R = L.block(j, 0, k - j, j);
110  Block_ D = L.block(j, j, k - j, k - j);
111  Block_ B = L.block(k, 0, M_ - k, j);
112  Block_ C = L.block(k, j, M_ - k, k - j);
113  Block_ Rbar = Lbar.block(j, 0, k - j, j);
114  Block_ Dbar = Lbar.block(j, j, k - j, k - j);
115  Block_ Bbar = Lbar.block(k, 0, M_ - k, j);
116  Block_ Cbar = Lbar.block(k, j, M_ - k, k - j);
117  if (Cbar.size() > 0) {
118  Cbar
119  = D.transpose().triangularView<Upper>()
120  .solve(Cbar.transpose()).transpose();
121  Bbar.noalias() -= Cbar * R;
122  Dbar.noalias() -= Cbar.transpose() * C;
123  }
124  symbolic_rev(D, Dbar);
125  Rbar.noalias() -= Cbar.transpose() * B;
126  Rbar.noalias() -= Dbar.selfadjointView<Lower>() * R;
127  Dbar.diagonal() *= 0.5;
128  Dbar.triangularView<StrictlyUpper>().setZero();
129  }
130  pos = 0;
131  for (size_type j = 0; j < M_; ++j)
132  for (size_type i = j; i < M_; ++i)
133  variRefA_[pos++]->adj_ += Lbar.coeffRef(i, j);
134  }
135  };
136 
137  class cholesky_scalar : public vari {
138  public:
139  int M_;
142 
157  cholesky_scalar(const Eigen::Matrix<var, -1, -1>& A,
158  const Eigen::Matrix<double, -1, -1>& L_A)
159  : vari(0.0),
160  M_(A.rows()),
161  variRefA_(ChainableStack::memalloc_.alloc_array<vari*>
162  (A.rows() * (A.rows() + 1) / 2)),
163  variRefL_(ChainableStack::memalloc_.alloc_array<vari*>
164  (A.rows() * (A.rows() + 1) / 2)) {
165  size_t accum = 0;
166  size_t accum_i = accum;
167  for (size_type j = 0; j < M_; ++j) {
168  for (size_type i = j; i < M_; ++i) {
169  accum_i += i;
170  size_t pos = j + accum_i;
171  variRefA_[pos] = A.coeffRef(i, j).vi_;
172  variRefL_[pos] = new vari(L_A.coeffRef(i, j), false);
173  }
174  accum += j;
175  accum_i = accum;
176  }
177  }
178 
191  virtual void chain() {
192  using Eigen::Matrix;
193  using Eigen::RowMajor;
194  Matrix<double, -1, -1, RowMajor> adjL(M_, M_);
195  Matrix<double, -1, -1, RowMajor> LA(M_, M_);
196  Matrix<double, -1, -1, RowMajor> adjA(M_, M_);
197  size_t pos = 0;
198  for (size_type i = 0; i < M_; ++i) {
199  for (size_type j = 0; j <= i; ++j) {
200  adjL.coeffRef(i, j) = variRefL_[pos]->adj_;
201  LA.coeffRef(i, j) = variRefL_[pos]->val_;
202  ++pos;
203  }
204  }
205 
206  --pos;
207  for (int i = M_ - 1; i >= 0; --i) {
208  for (int j = i; j >= 0; --j) {
209  if (i == j) {
210  adjA.coeffRef(i, j) = 0.5 * adjL.coeff(i, j)
211  / LA.coeff(i, j);
212  } else {
213  adjA.coeffRef(i, j) = adjL.coeff(i, j)
214  / LA.coeff(j, j);
215  adjL.coeffRef(j, j) -= adjL.coeff(i, j)
216  * LA.coeff(i, j) / LA.coeff(j, j);
217  }
218  for (int k = j - 1; k >=0; --k) {
219  adjL.coeffRef(i, k) -= adjA.coeff(i, j)
220  * LA.coeff(j, k);
221  adjL.coeffRef(j, k) -= adjA.coeff(i, j)
222  * LA.coeff(i, k);
223  }
224  variRefA_[pos--]->adj_ += adjA.coeffRef(i, j);
225  }
226  }
227  }
228  };
229 
243  inline Eigen::Matrix<var, -1, -1>
244  cholesky_decompose(const Eigen::Matrix<var, -1, -1> &A) {
245  check_square("cholesky_decompose", "A", A);
246  check_symmetric("cholesky_decompose", "A", A);
247 
248  Eigen::Matrix<double, -1, -1> L_A(value_of_rec(A));
249  Eigen::LLT<Eigen::MatrixXd> L_factor
250  = L_A.selfadjointView<Eigen::Lower>().llt();
251  check_pos_definite("cholesky_decompose", "m", L_factor);
252  L_A = L_factor.matrixL();
253 
254  // Memory allocated in arena.
255  // cholesky_scalar gradient faster for small matrices compared to
256  // cholesky_block
257  vari* dummy = new vari(0.0, false);
258  Eigen::Matrix<var, -1, -1> L(A.rows(), A.cols());
259  if (L_A.rows() <= 35) {
260  cholesky_scalar *baseVari
261  = new cholesky_scalar(A, L_A);
262  size_t accum = 0;
263  size_t accum_i = accum;
264  for (size_type j = 0; j < L.cols(); ++j) {
265  for (size_type i = j; i < L.cols(); ++i) {
266  accum_i += i;
267  size_t pos = j + accum_i;
268  L.coeffRef(i, j).vi_ = baseVari->variRefL_[pos];
269  }
270  for (size_type k = 0; k < j; ++k)
271  L.coeffRef(k, j).vi_ = dummy;
272  accum += j;
273  accum_i = accum;
274  }
275  } else {
276  cholesky_block *baseVari
277  = new cholesky_block(A, L_A);
278  size_t pos = 0;
279  for (size_type j = 0; j < L.cols(); ++j) {
280  for (size_type i = j; i < L.cols(); ++i) {
281  L.coeffRef(i, j).vi_ = baseVari->variRefL_[pos++];
282  }
283  for (size_type k = 0; k < j; ++k)
284  L.coeffRef(k, j).vi_ = dummy;
285  }
286  }
287  return L;
288  }
289  }
290 }
291 #endif
virtual void chain()
Reverse mode differentiation algorithm refernce:
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
cholesky_block(const Eigen::Matrix< var, -1, -1 > &A, const Eigen::Matrix< double, -1, -1 > &L_A)
Constructor for cholesky function.
int min(const std::vector< int > &x)
Returns the minimum coefficient in the specified column vector.
Definition: min.hpp:20
double value_of_rec(const fvar< T > &v)
Return the value of the specified variable.
cholesky_scalar(const Eigen::Matrix< var, -1, -1 > &A, const Eigen::Matrix< double, -1, -1 > &L_A)
Constructor for cholesky function.
The variable implementation base class.
Definition: vari.hpp:30
friend class var
Definition: vari.hpp:32
const double val_
The value of this variable.
Definition: vari.hpp:38
Empty struct for use in boost::condtional<is_constant_struct<T1>::value, T1, dummy>::type as false co...
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_symmetric(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Check if the specified matrix is symmetric.
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.
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > cholesky_decompose(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Return the lower-triangular Cholesky factor (i.e., matrix square root) of the specified square...
int max(const std::vector< int > &x)
Returns the maximum coefficient in the specified column vector.
Definition: max.hpp:22
void symbolic_rev(Block_ &L, Block_ &Lbar)
Symbolic adjoint calculation for cholesky factor A.
vari(double x)
Construct a variable implementation from a value.
Definition: vari.hpp:58
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
virtual void chain()
Reverse mode differentiation algorithm refernce:
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:12
void check_pos_definite(const char *function, const char *name, const Eigen::Matrix< T_y, -1, -1 > &y)
Check if the specified square, symmetric matrix is positive definite.
Eigen::Block< Eigen::MatrixXd > Block_

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