LDLT_factor.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_LDLT_FACTOR_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_LDLT_FACTOR_HPP
3 
7 #include <boost/shared_ptr.hpp>
8 
9 namespace stan {
10 namespace math {
11 
12 /**
13  * LDLT_factor is a thin wrapper on Eigen::LDLT to allow for
14  * reusing factorizations and efficient autodiff of things like
15  * log determinants and solutions to linear systems.
16  *
17  * Memory is allocated in the constructor and stored in a
18  * <code>boost::shared_ptr</code>, which ensures that is freed
19  * when the object is released.
20  *
21  * After the constructor and/or compute() is called users of
22  * LDLT_factor are responsible for calling success() to
23  * check whether the factorization has succeeded. Use of an LDLT_factor
24  * object (e.g., in mdivide_left_ldlt) is undefined if success() is false.
25  *
26  * It's usage pattern is:
27  *
28  * ~~~
29  * Eigen::Matrix<T, R, C> A1, A2;
30  *
31  * LDLT_factor<T, R, C> ldlt_A1(A1);
32  * LDLT_factor<T, R, C> ldlt_A2;
33  * ldlt_A2.compute(A2);
34  * ~~~
35  *
36  * The caller should check that ldlt_A1.success() and ldlt_A2.success()
37  * are true or abort accordingly. Alternatively, call check_ldlt_factor().
38  *
39  * Note that ldlt_A1 and ldlt_A2 are completely equivalent. They simply
40  * demonstrate two different ways to construct the factorization.
41  *
42  * The caller can use the LDLT_factor objects as needed. For
43  * instance
44  *
45  * ~~~
46  * x1 = mdivide_left_ldlt(ldlt_A1, b1);
47  * x2 = mdivide_right_ldlt(b2, ldlt_A2);
48  *
49  * d1 = log_determinant_ldlt(ldlt_A1);
50  * d2 = log_determinant_ldlt(ldlt_A2);
51  * ~~~
52  *
53  * This class is conceptually similar to the corresponding Eigen
54  * class. Any symmetric, positive-definite matrix A can be
55  * decomposed as LDL' where L is unit lower-triangular and D is
56  * diagonal with positive diagonal elements.
57  *
58  * @tparam T scalare type held in the matrix
59  * @tparam R rows (as in Eigen)
60  * @tparam C columns (as in Eigen)
61  */
62 template <typename T, int R, int C>
63 class LDLT_factor {
64  public:
65  typedef Eigen::Matrix<T, Eigen::Dynamic, 1> vector_t;
66  typedef Eigen::Matrix<T, R, C> matrix_t;
67  typedef Eigen::LDLT<matrix_t> ldlt_t;
68  typedef size_t size_type;
69  typedef double value_type;
70 
71  LDLT_factor() : N_(0), ldltP_(new ldlt_t()) {}
72 
73  explicit LDLT_factor(const matrix_t& A) : N_(0), ldltP_(new ldlt_t()) {
74  compute(A);
75  }
76 
77  inline void compute(const matrix_t& A) {
78  check_square("LDLT_factor", "A", A);
79  N_ = A.rows();
80  ldltP_->compute(A);
81  }
82 
83  inline bool success() const {
84  if (ldltP_->info() != Eigen::Success)
85  return false;
86  if (!(ldltP_->isPositive()))
87  return false;
88  vector_t ldltP_diag(ldltP_->vectorD());
89  for (int i = 0; i < ldltP_diag.size(); ++i)
90  if (ldltP_diag(i) <= 0 || is_nan(ldltP_diag(i)))
91  return false;
92  return true;
93  }
94 
95  inline T log_abs_det() const { return ldltP_->vectorD().array().log().sum(); }
96 
97  inline void inverse(matrix_t& invA) const {
98  invA.setIdentity(N_);
99  ldltP_->solveInPlace(invA);
100  }
101 
102 #if EIGEN_VERSION_AT_LEAST(3, 3, 0)
103  template <typename Rhs>
104  inline const Eigen::Solve<ldlt_t, Rhs> solve(
105  const Eigen::MatrixBase<Rhs>& b) const {
106  return ldltP_->solve(b);
107  }
108 #else
109  template <typename Rhs>
110  inline const Eigen::internal::solve_retval<ldlt_t, Rhs> solve(
111  const Eigen::MatrixBase<Rhs>& b) const {
112  return ldltP_->solve(b);
113  }
114 #endif
115 
116  inline matrix_t solveRight(const matrix_t& B) const {
117  return ldltP_->solve(B.transpose()).transpose();
118  }
119 
120  inline vector_t vectorD() const { return ldltP_->vectorD(); }
121 
122  inline ldlt_t matrixLDLT() const { return ldltP_->matrixLDLT(); }
123 
124  inline size_t rows() const { return N_; }
125  inline size_t cols() const { return N_; }
126 
127  size_t N_;
128  boost::shared_ptr<ldlt_t> ldltP_;
129 };
130 
131 } // namespace math
132 } // namespace stan
133 #endif
Eigen::LDLT< matrix_t > ldlt_t
Definition: LDLT_factor.hpp:67
matrix_t solveRight(const matrix_t &B) const
Eigen::Matrix< T, Eigen::Dynamic, 1 > vector_t
Definition: LDLT_factor.hpp:65
ldlt_t matrixLDLT() const
vector_t vectorD() const
void compute(const matrix_t &A)
Definition: LDLT_factor.hpp:77
Eigen::Matrix< T, R, C > matrix_t
Definition: LDLT_factor.hpp:66
boost::shared_ptr< ldlt_t > ldltP_
void check_square(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
static const double A
Definition: Units.h:82
void inverse(matrix_t &invA) const
Definition: LDLT_factor.hpp:97
LDLT_factor(const matrix_t &A)
Definition: LDLT_factor.hpp:73
const hit & b
Definition: hits.cxx:21
double T
Definition: Xdiff_gwt.C:5
int is_nan(const fvar< T > &x)
Definition: is_nan.hpp:19
size_t rows() const
Eigen::Matrix< T, C, R > transpose(const Eigen::Matrix< T, R, C > &m)
Definition: transpose.hpp:10
const Eigen::internal::solve_retval< ldlt_t, Rhs > solve(const Eigen::MatrixBase< Rhs > &b) const
size_t cols() const