1 #ifndef STAN_MATH_PRIM_FUN_INVERSE_SPD_HPP
2 #define STAN_MATH_PRIM_FUN_INVERSE_SPD_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/Eigen.hpp>
7
8 namespace stan {
9 namespace math {
10
11 /**
12 * Returns the inverse of the specified symmetric, pos/neg-definite matrix.
13 *
14 * @tparam EigMat type of elements in the matrix
15 * @param m specified matrix
16 * @return Inverse of the matrix (an empty matrix if the specified matrix has
17 * size zero).
18 * @throw std::invalid_argument if the matrix is not symmetric.
19 * @throw std::domain_error if the matrix is not positive definite.
20 */
21 template <typename EigMat>
22 inline Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic>
inverse_spd(const EigMat & m)23 inverse_spd(const EigMat& m) {
24 using Eigen::Dynamic;
25 using Eigen::LDLT;
26 using Eigen::Matrix;
27 using Scalar = value_type_t<EigMat>;
28 if (m.size() == 0) {
29 return {};
30 }
31 const Eigen::Ref<const plain_type_t<EigMat>>& m_ref = m;
32 check_symmetric("inverse_spd", "m", m_ref);
33 plain_type_t<EigMat> mmt = 0.5 * (m_ref + m_ref.transpose());
34 LDLT<plain_type_t<EigMat>> ldlt(mmt);
35 if (ldlt.info() != Eigen::Success) {
36 throw_domain_error("invese_spd", "LDLT factor failed", "", "");
37 }
38 if (!ldlt.isPositive()) {
39 throw_domain_error("invese_spd", "matrix not positive definite", "", "");
40 }
41 Matrix<Scalar, Dynamic, 1> diag_ldlt = ldlt.vectorD();
42 check_positive("inverse_spd", "matrix not positive definite", diag_ldlt);
43
44 return ldlt.solve(
45 Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>::Identity(
46 m.rows(), m.cols()));
47 }
48
49 } // namespace math
50 } // namespace stan
51
52 #endif
53