1 #ifndef STAN_MATH_PRIM_FUN_GENERALIZED_INVERSE_HPP 2 #define STAN_MATH_PRIM_FUN_GENERALIZED_INVERSE_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 #include <stan/math/prim/fun/to_ref.hpp> 8 #include <stan/math/prim/fun/inverse.hpp> 9 10 namespace stan { 11 namespace math { 12 13 /** 14 * Returns the Moore-Penrose generalized inverse of the specified matrix. 15 * 16 * The method is based on the Cholesky computation of the transform as specified 17 * in 18 * 19 * <ul><li> Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse 20 Matrices. 21 * <i>arXiv</i> <b>0804.4809</b> </li></ul> 22 * 23 * @tparam EigMat type of the matrix (must be derived from `Eigen::MatrixBase`) 24 * 25 * @param G specified matrix 26 * @return Generalized inverse of the matrix (an empty matrix if the specified 27 * matrix has size zero). 28 */ 29 template <typename EigMat, require_eigen_t<EigMat>* = nullptr, 30 require_not_vt_var<EigMat>* = nullptr> 31 inline Eigen::Matrix<value_type_t<EigMat>, EigMat::ColsAtCompileTime, 32 EigMat::RowsAtCompileTime> generalized_inverse(const EigMat & G)33generalized_inverse(const EigMat& G) { 34 const auto& G_ref = to_ref(G); 35 if (G_ref.size() == 0) 36 return {}; 37 38 if (G_ref.rows() == G_ref.cols()) { 39 Eigen::CompleteOrthogonalDecomposition< 40 Eigen::Matrix<value_type_t<EigMat>, EigMat::RowsAtCompileTime, 41 EigMat::ColsAtCompileTime>> 42 complete_ortho_decomp_G = G_ref.completeOrthogonalDecomposition(); 43 if (!(complete_ortho_decomp_G.rank() < G_ref.rows())) 44 return inverse(G_ref); 45 else 46 return complete_ortho_decomp_G.pseudoInverse(); 47 } 48 49 if (G_ref.rows() < G_ref.cols()) { 50 return (G_ref * G_ref.transpose()).ldlt().solve(G_ref).transpose(); 51 } else { 52 return (G_ref.transpose() * G_ref).ldlt().solve(G_ref.transpose()); 53 } 54 } 55 56 } // namespace math 57 } // namespace stan 58 59 #endif 60