1 #ifndef STAN_MATH_REV_MAT_FUN_CHOLESKY_DECOMPOSE_HPP 2 #define STAN_MATH_REV_MAT_FUN_CHOLESKY_DECOMPOSE_HPP 23 typedef Eigen::Block<Eigen::MatrixXd>
Block_;
44 const Eigen::Matrix<double, -1, -1>& L_A)
52 block_size_ =
std::max((M_ / 8 / 16) * 16, 8);
53 block_size_ =
std::min(block_size_, 128);
56 variRefA_[pos] = A.coeffRef(i, j).vi_;
57 variRefL_[pos] =
new vari(L_A.coeffRef(i, j),
false); ++pos;
72 using Eigen::StrictlyUpper;
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());
88 using Eigen::MatrixXd;
92 using Eigen::StrictlyUpper;
93 MatrixXd Lbar(M_, M_);
101 Lbar.coeffRef(i, j) = variRefL_[pos]->
adj_;
102 L.coeffRef(i, j) = variRefL_[pos]->
val_;
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) {
119 = D.transpose().triangularView<Upper>()
121 Bbar.noalias() -= Cbar * R;
122 Dbar.noalias() -= Cbar.transpose() * C;
125 Rbar.noalias() -= Cbar.transpose() * B;
126 Rbar.noalias() -= Dbar.selfadjointView<Lower>() * R;
127 Dbar.diagonal() *= 0.5;
128 Dbar.triangularView<StrictlyUpper>().setZero();
133 variRefA_[pos++]->
adj_ += Lbar.coeffRef(i, j);
158 const Eigen::Matrix<double, -1, -1>& L_A)
166 size_t accum_i = accum;
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);
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_);
200 adjL.coeffRef(i, j) = variRefL_[pos]->
adj_;
201 LA.coeffRef(i, j) = variRefL_[pos]->
val_;
207 for (
int i = M_ - 1; i >= 0; --i) {
208 for (
int j = i; j >= 0; --j) {
210 adjA.coeffRef(i, j) = 0.5 * adjL.coeff(i, j)
213 adjA.coeffRef(i, j) = adjL.coeff(i, j)
215 adjL.coeffRef(j, j) -= adjL.coeff(i, j)
216 * LA.coeff(i, j) / LA.coeff(j, j);
218 for (
int k = j - 1; k >=0; --k) {
219 adjL.coeffRef(i, k) -= adjA.coeff(i, j)
221 adjL.coeffRef(j, k) -= adjA.coeff(i, j)
224 variRefA_[pos--]->
adj_ += adjA.coeffRef(i, j);
243 inline Eigen::Matrix<
var, -1, -1>
249 Eigen::LLT<Eigen::MatrixXd> L_factor
250 = L_A.selfadjointView<Eigen::Lower>().llt();
252 L_A = L_factor.matrixL();
258 Eigen::Matrix<
var, -1, -1> L(A.rows(), A.cols());
259 if (L_A.rows() <= 35) {
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) {
267 size_t pos = j + accum_i;
268 L.coeffRef(i, j).vi_ = baseVari->
variRefL_[pos];
271 L.coeffRef(k, j).vi_ = dummy;
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++];
284 L.coeffRef(k, j).vi_ = dummy;
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.
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.
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.
const double val_
The value of this variable.
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.
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.
void symbolic_rev(Block_ &L, Block_ &Lbar)
Symbolic adjoint calculation for cholesky factor A.
vari(double x)
Construct a variable implementation from a value.
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
virtual void chain()
Reverse mode differentiation algorithm refernce:
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
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_