1 #ifndef STAN_MATH_REV_FUN_LOG_DETERMINANT_SPD_HPP
2 #define STAN_MATH_REV_FUN_LOG_DETERMINANT_SPD_HPP
3
4 #include <stan/math/rev/meta.hpp>
5 #include <stan/math/rev/core.hpp>
6 #include <stan/math/rev/core/typedefs.hpp>
7 #include <stan/math/prim/err.hpp>
8 #include <stan/math/prim/fun/Eigen.hpp>
9 #include <stan/math/prim/fun/log.hpp>
10 #include <stan/math/prim/fun/sum.hpp>
11 #include <stan/math/prim/fun/typedefs.hpp>
12
13 namespace stan {
14 namespace math {
15
16 /**
17 * Returns the log det of a symmetric, positive-definite matrix
18 *
19 * @tparam EigMat Type of the matrix
20 * @param m a symmetric, positive-definite matrix
21 * @return The log determinant of the specified matrix
22 */
23 template <typename EigMat, require_eigen_vt<is_var, EigMat>* = nullptr>
log_determinant_spd(const EigMat & m)24 inline var log_determinant_spd(const EigMat& m) {
25 if (m.size() == 0) {
26 return 0;
27 }
28
29 matrix_d m_d = m.val();
30 check_symmetric("log_determinant_spd", "m", m_d);
31
32 Eigen::LDLT<matrix_d> ldlt(m_d);
33 if (ldlt.info() != Eigen::Success) {
34 double y = 0;
35 throw_domain_error("log_determinant_spd", "matrix argument", y,
36 "failed LDLT factorization");
37 }
38
39 // compute the inverse of A (needed for the derivative)
40 m_d.setIdentity(m.rows(), m.cols());
41 ldlt.solveInPlace(m_d);
42
43 if (ldlt.isNegative() || (ldlt.vectorD().array() <= 1e-16).any()) {
44 double y = 0;
45 throw_domain_error("log_determinant_spd", "matrix argument", y,
46 "matrix is negative definite");
47 }
48
49 double val = sum(log(ldlt.vectorD()));
50
51 check_finite("log_determinant_spd",
52 "log determininant of the matrix argument", val);
53
54 vari** operands
55 = ChainableStack::instance_->memalloc_.alloc_array<vari*>(m.size());
56 Eigen::Map<matrix_vi>(operands, m.rows(), m.cols()) = m.vi();
57
58 double* gradients
59 = ChainableStack::instance_->memalloc_.alloc_array<double>(m.size());
60 Eigen::Map<matrix_d>(gradients, m.rows(), m.cols()) = m_d;
61
62 return var(
63 new precomputed_gradients_vari(val, m.size(), operands, gradients));
64 }
65
66 /**
67 * Returns the log det of a symmetric, positive-definite matrix
68 *
69 * @tparam EigMat Type of the matrix
70 * @param m a symmetric, positive-definite matrix
71 * @return The log determinant of the specified matrix
72 */
73 template <typename T, require_var_matrix_t<T>* = nullptr>
log_determinant_spd(const T & m)74 inline var log_determinant_spd(const T& m) {
75 check_square("log_determinant_spd", "m", m);
76
77 if (m.size() == 0) {
78 return var(0.0);
79 }
80
81 check_symmetric("log_determinant_spd", "m", m.val());
82
83 auto ldlt = m.val().ldlt();
84 if (ldlt.info() != Eigen::Success) {
85 double y = 0;
86 throw_domain_error("log_determinant_spd", "matrix argument", y,
87 "failed LDLT factorization");
88 }
89
90 if (ldlt.isNegative() || (ldlt.vectorD().array() <= 1e-16).any()) {
91 double y = 0;
92 throw_domain_error("log_determinant_spd", "matrix argument", y,
93 "matrix is negative definite");
94 }
95
96 arena_t<Eigen::MatrixXd> arena_m_inv_transpose
97 = Eigen::MatrixXd::Identity(m.rows(), m.cols());
98 ldlt.solveInPlace(arena_m_inv_transpose);
99
100 return make_callback_var(sum(log(ldlt.vectorD())),
101 [m, arena_m_inv_transpose](const auto& res) mutable {
102 m.adj() += res.adj() * arena_m_inv_transpose;
103 });
104 }
105
106 } // namespace math
107 } // namespace stan
108 #endif
109