Merge pull request #2577 from spinkney/generalized_inverse_fix

check generalized inverse for full rank symmetric mat

Author: spinkney

stan/math/prim/fun/generalized_inverse.hpp

@@ -31,15 +31,22 @@ template <typename EigMat, require_eigen_t<EigMat>* = nullptr,
 inline Eigen::Matrix<value_type_t<EigMat>, EigMat::ColsAtCompileTime,
                      EigMat::RowsAtCompileTime>
 generalized_inverse(const EigMat& G) {
-  if (G.size() == 0)
+  const auto& G_ref = to_ref(G);
+  if (G_ref.size() == 0)
     return {};
 
-  if (G.rows() == G.cols())
-    return inverse(G);
-
-  const auto& G_ref = to_ref(G);
+  if (G_ref.rows() == G_ref.cols()) {
+    Eigen::CompleteOrthogonalDecomposition<
+        Eigen::Matrix<value_type_t<EigMat>, EigMat::RowsAtCompileTime,
+                      EigMat::ColsAtCompileTime>>
+        complete_ortho_decomp_G = G_ref.completeOrthogonalDecomposition();
+    if (!(complete_ortho_decomp_G.rank() < G_ref.rows()))
+      return inverse(G_ref);
+    else
+      return complete_ortho_decomp_G.pseudoInverse();
+  }
 
-  if (G.rows() < G.cols()) {
+  if (G_ref.rows() < G_ref.cols()) {
     return (G_ref * G_ref.transpose()).ldlt().solve(G_ref).transpose();
   } else {
     return (G_ref.transpose() * G_ref).ldlt().solve(G_ref.transpose());

stan/math/rev/fun/generalized_inverse.hpp

@@ -1,6 +1,7 @@
 #ifndef STAN_MATH_REV_FUN_GENERALIZED_INVERSE_HPP
 #define STAN_MATH_REV_FUN_GENERALIZED_INVERSE_HPP
 
+#include <stan/math/prim/fun/add_diag.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
@@ -66,10 +67,20 @@ inline auto generalized_inverse(const VarMat& G) {
   if (G.size() == 0)
     return ret_type(G);
 
-  if (G.rows() == G.cols())
-    return ret_type(inverse(G));
-
-  if (G.rows() < G.cols()) {
+  if (G.rows() == G.cols()) {
+    arena_t<VarMat> G_arena(G);
+    Eigen::CompleteOrthogonalDecomposition<Eigen::MatrixXd>
+        complete_ortho_decomp_G
+        = G_arena.val().completeOrthogonalDecomposition();
+    if (!(complete_ortho_decomp_G.rank() < G.rows())) {
+      return ret_type(inverse(G_arena));
+    } else {
+      arena_t<ret_type> inv_G(complete_ortho_decomp_G.pseudoInverse());
+      reverse_pass_callback(
+          internal::generalized_inverse_lambda(G_arena, inv_G));
+      return ret_type(inv_G);
+    }
+  } else if (G.rows() < G.cols()) {
     arena_t<VarMat> G_arena(G);
     arena_t<ret_type> inv_G((G_arena.val_op() * G_arena.val_op().transpose())
                                 .ldlt()

test/unit/math/mix/fun/generalized_inverse_test.cpp

@@ -2,7 +2,7 @@
 #include <vector>
 #include <gtest/gtest.h>
 
-TEST(mathMixMatFun, ad_tests) {
+TEST(mathMixMatFun, gen_inv_test_1) {
   using stan::test::expect_ad;
   using stan::test::expect_ad_matvar;
 
@@ -11,25 +11,61 @@ TEST(mathMixMatFun, ad_tests) {
   Eigen::MatrixXd t(0, 0);
   expect_ad(f, t);
   expect_ad_matvar(f, t);
+}
+
+TEST(mathMixMatFun, gen_inv_test_2) {
+  using stan::test::expect_ad;
+  using stan::test::expect_ad_matvar;
+
+  auto f = [](const auto& G) { return stan::math::generalized_inverse(G); };
 
   Eigen::MatrixXd u(1, 1);
   u << 2;
   expect_ad(f, u);
   expect_ad_matvar(f, u);
+}
+
+TEST(mathMixMatFun, gen_inv_test_3) {
+  using stan::test::expect_ad;
+  using stan::test::expect_ad_matvar;
+
+  auto f = [](const auto& G) { return stan::math::generalized_inverse(G); };
 
   Eigen::MatrixXd v(2, 3);
   v << 1, 3, 5, 2, 4, 6;
   expect_ad(f, v);
   expect_ad_matvar(f, v);
+}
+
+TEST(mathMixMatFun, gen_inv_test_4) {
+  using stan::test::expect_ad;
+  using stan::test::expect_ad_matvar;
 
+  auto f = [](const auto& G) { return stan::math::generalized_inverse(G); };
+
+  Eigen::MatrixXd v(2, 3);
   v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1;
   expect_ad(f, v);
   expect_ad_matvar(f, v);
+}
+
+TEST(mathMixMatFun, gen_inv_test_5) {
+  using stan::test::expect_ad;
+  using stan::test::expect_ad_matvar;
+
+  auto f = [](const auto& G) { return stan::math::generalized_inverse(G); };
 
   Eigen::MatrixXd s(2, 4);
   s << 3.4, 2, 5, 1.2, 2, 1, 3.2, 3.1;
   expect_ad(f, s);
   expect_ad_matvar(f, s);
+}
+
+TEST(mathMixMatFun, gen_inv_test_6) {
+  using stan::test::expect_ad;
+  using stan::test::expect_ad_matvar;
+
+  auto f = [](const auto& G) { return stan::math::generalized_inverse(G); };
 
   // issues around zero require looser tolerances for hessians
   stan::test::ad_tolerances tols;
@@ -40,10 +76,24 @@ TEST(mathMixMatFun, ad_tests) {
   w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29;
   expect_ad(tols, f, w);
   expect_ad_matvar(f, w);
+}
+
+TEST(mathMixMatFun, gen_inv_test_8) {
+  using stan::test::expect_ad;
+  using stan::test::expect_ad_matvar;
+
+  auto f = [](const auto& G) { return stan::math::generalized_inverse(G); };
 
   Eigen::MatrixXd z(2, 2);
   z << 1, 2, 5, std::numeric_limits<double>::quiet_NaN();
   EXPECT_NO_THROW(stan::math::generalized_inverse(z));
+}
+
+TEST(mathMixMatFun, gen_inv_test_9) {
+  using stan::test::expect_ad;
+  using stan::test::expect_ad_matvar;
+
+  auto f = [](const auto& G) { return stan::math::generalized_inverse(G); };
 
   // autodiff throws, so following fails (throw behavior must match to pass)
 

test/unit/math/prim/fun/generalized_inverse_test.cpp

@@ -45,3 +45,17 @@ TEST(MathMatrixPrim, Equal2) {
   stan::math::matrix_d m3 = m1 * generalized_inverse(m1);
   EXPECT_MATRIX_NEAR(m2, m3, 1e-9);
 }
+
+TEST(MathMatrixPrim, Equal3) {
+  using stan::math::generalized_inverse;
+
+  stan::math::matrix_d m1(3, 3);
+  m1 << 4, 4, -2, 4, 4, -2, -2, -2, 10;
+
+  stan::math::matrix_d m2(3, 3);
+  m2 << 0.069444444, 0.069444444, 0.0277777777, 0.069444444, 0.069444444,
+      0.0277777777, 0.027777777, 0.027777777, 0.1111111111;
+
+  stan::math::matrix_d m3 = generalized_inverse(m1);
+  EXPECT_MATRIX_NEAR(m2, m3, 1e-9);
+}