Stan Math Library  2.9.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 
15 namespace stan {
16  namespace math {
17 
19  public:
20  int M_; // A.rows() = A.cols()
23 
24  /* ctor for cholesky function
25  *
26  * Stores varis for A
27  * Instantiates and stores varis for L
28  * Instantiates and stores dummy vari for
29  * upper triangular part of var result returned
30  * in cholesky_decompose function call
31  *
32  * variRefL aren't on the chainable
33  * autodiff stack, only used for storage
34  * and computation. Note that varis for
35  * L are constructed externally in
36  * cholesky_decompose.
37  *
38  * @param matrix A
39  * @param matrix L, cholesky factor of A
40  * */
41  cholesky_decompose_v_vari(const Eigen::Matrix<var, -1, -1>& A,
42  const Eigen::Matrix<double, -1, -1>& L_A)
43  : vari(0.0),
44  M_(A.rows()),
45  variRefA_(ChainableStack::memalloc_.alloc_array<vari*>
46  (A.rows() * (A.rows() + 1) / 2)),
47  variRefL_(ChainableStack::memalloc_.alloc_array<vari*>
48  (A.rows() * (A.rows() + 1) / 2)) {
49  size_t accum = 0;
50  size_t accum_i = accum;
51  for (size_type j = 0; j < M_; ++j) {
52  for (size_type i = j; i < M_; ++i) {
53  accum_i += i;
54  size_t pos = j + accum_i;
55  variRefA_[pos] = A.coeffRef(i, j).vi_;
56  variRefL_[pos] = new vari(L_A.coeffRef(i, j), false);
57  }
58  accum += j;
59  accum_i = accum;
60  }
61  }
62 
63  /* Reverse mode differentiation
64  * algorithm refernce:
65  *
66  * Mike Giles. An extended collection of matrix
67  * derivative results for forward and reverse mode AD.
68  * Jan. 2008.
69  *
70  * Note algorithm as laid out in Giles is
71  * row-major, so Eigen::Matrices are explicitly storage
72  * order RowMajor, whereas Eigen defaults to
73  * ColumnMajor. Also note algorithm
74  * starts by calculating the adjoint for
75  * A(M_ - 1, M_ - 1), hence pos on line 94 is decremented
76  * to start at pos = M_ * (M_ + 1) / 2.
77  * */
78  virtual void chain() {
79  using Eigen::Matrix;
80  using Eigen::RowMajor;
81  Matrix<double, -1, -1, RowMajor> adjL(M_, M_);
82  Matrix<double, -1, -1, RowMajor> LA(M_, M_);
83  Matrix<double, -1, -1, RowMajor> adjA(M_, M_);
84  size_t pos = 0;
85  for (size_type i = 0; i < M_; ++i) {
86  for (size_type j = 0; j <= i; ++j) {
87  adjL.coeffRef(i, j) = variRefL_[pos]->adj_;
88  LA.coeffRef(i, j) = variRefL_[pos]->val_;
89  ++pos;
90  }
91  }
92 
93  --pos;
94  for (int i = M_ - 1; i >= 0; --i) {
95  for (int j = i; j >= 0; --j) {
96  if (i == j) {
97  adjA.coeffRef(i, j) = 0.5 * adjL.coeff(i, j)
98  / LA.coeff(i, j);
99  } else {
100  adjA.coeffRef(i, j) = adjL.coeff(i, j)
101  / LA.coeff(j, j);
102  adjL.coeffRef(j, j) -= adjL.coeff(i, j)
103  * LA.coeff(i, j) / LA.coeff(j, j);
104  }
105  for (int k = j - 1; k >=0; --k) {
106  adjL.coeffRef(i, k) -= adjA.coeff(i, j)
107  * LA.coeff(j, k);
108  adjL.coeffRef(j, k) -= adjA.coeff(i, j)
109  * LA.coeff(i, k);
110  }
111  variRefA_[pos--]->adj_ += adjA.coeffRef(i, j);
112  }
113  }
114  }
115  };
116 
117  /* Reverse mode specialization of
118  * cholesky decomposition
119  *
120  * Internally calls llt rather than using
121  * stan::math::cholesky_decompose in order
122  * to use selfadjointView<Lower> optimization.
123  *
124  * Note chainable stack varis are created
125  * below in Matrix<var, -1, -1>
126  *
127  * @param Matrix A
128  * @return L cholesky factor of A
129  */
130  Eigen::Matrix<var, -1, -1>
131  cholesky_decompose(const Eigen::Matrix<var, -1, -1> &A) {
132  stan::math::check_square("cholesky_decompose", "A", A);
133  stan::math::check_symmetric("cholesky_decompose", "A", A);
134 
135  Eigen::Matrix<double, -1, -1> L_A(value_of_rec(A));
136  Eigen::LLT<Eigen::MatrixXd> L_factor
137  = L_A.selfadjointView<Eigen::Lower>().llt();
138  check_pos_definite("cholesky_decompose", "m", L_factor);
139  L_A = L_factor.matrixL();
140 
141  // NOTE: this is not a memory leak, this vari is used in the
142  // expression graph to evaluate the adjoint, but is not needed
143  // for the returned matrix. Memory will be cleaned up with the
144  // arena allocator.
145  cholesky_decompose_v_vari *baseVari
146  = new cholesky_decompose_v_vari(A, L_A);
147  stan::math::vari dummy(0.0, false);
148  Eigen::Matrix<var, -1, -1> L(A.rows(), A.cols());
149  size_t accum = 0;
150  size_t accum_i = accum;
151  for (size_type j = 0; j < L.cols(); ++j) {
152  for (size_type i = j; i < L.cols(); ++i) {
153  accum_i += i;
154  size_t pos = j + accum_i;
155  L.coeffRef(i, j).vi_ = baseVari->variRefL_[pos];
156  }
157  for (size_type k = 0; k < j; ++k)
158  L.coeffRef(k, j).vi_ = &dummy;
159  accum += j;
160  accum_i = accum;
161  }
162  return L;
163  }
164  }
165 }
166 #endif
vari(const double x)
Construct a variable implementation from a value.
Definition: vari.hpp:58
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. ...
cholesky_decompose_v_vari(const Eigen::Matrix< var,-1,-1 > &A, const Eigen::Matrix< double,-1,-1 > &L_A)
double value_of_rec(const fvar< T > &v)
Return the value of the specified variable.
The variable implementation base class.
Definition: vari.hpp:30
const double val_
The value of this variable.
Definition: vari.hpp:38
Empty struct for use in boost::condtional::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
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...
bool check_pos_definite(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Return true if the specified square, symmetric matrix is positive definite.
bool check_symmetric(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Return true if the specified matrix is symmetric.
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
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.

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