Updates from review round 2

WardBrian · WardBrian · commit 2bd56dad5ae8 · 2023-08-17T15:53:29.000-04:00
diff --git a/stan/math/prim/fun/complex_schur_decompose.hpp b/stan/math/prim/fun/complex_schur_decompose.hpp
@@ -27,13 +27,13 @@ namespace math {
 template <typename M, require_eigen_dense_dynamic_t<M>* = nullptr>
 inline Eigen::Matrix<complex_return_t<scalar_type_t<M>>, -1, -1>
 complex_schur_decompose_u(const M& m) {
-  if (m.size() == 0)
+  if (unlikely(m.size() == 0)) {
     return m;
+  }
   check_square("complex_schur_decompose_u", "m", m);
   using MatType = Eigen::Matrix<scalar_type_t<M>, -1, -1>;
   // copy because ComplexSchur requires Eigen::Matrix type
-  MatType mv = m;
-  Eigen::ComplexSchur<MatType> cs(mv);
+  Eigen::ComplexSchur<MatType> cs{MatType(m)};
   return cs.matrixU();
 }
 
@@ -51,13 +51,13 @@ complex_schur_decompose_u(const M& m) {
 template <typename M, require_eigen_dense_dynamic_t<M>* = nullptr>
 inline Eigen::Matrix<complex_return_t<scalar_type_t<M>>, -1, -1>
 complex_schur_decompose_t(const M& m) {
-  if (m.size() == 0)
+  if (unlikely(m.size() == 0)) {
     return m;
+  }
   check_square("complex_schur_decompose_t", "m", m);
   using MatType = Eigen::Matrix<scalar_type_t<M>, -1, -1>;
   // copy because ComplexSchur requires Eigen::Matrix type
-  MatType mv = m;
-  Eigen::ComplexSchur<MatType> cs(mv, false);
+  Eigen::ComplexSchur<MatType> cs{MatType(m), false};
   return cs.matrixT();
 }
 
@@ -81,13 +81,13 @@ template <typename M, require_eigen_dense_dynamic_t<M>* = nullptr>
 inline std::tuple<Eigen::Matrix<complex_return_t<scalar_type_t<M>>, -1, -1>,
                   Eigen::Matrix<complex_return_t<scalar_type_t<M>>, -1, -1>>
 complex_schur_decompose(const M& m) {
-  if (m.size() == 0)
+  if (unlikely(m.size() == 0)) {
     return std::make_tuple(m, m);
+  }
   check_square("complex_schur_decompose", "m", m);
   using MatType = Eigen::Matrix<scalar_type_t<M>, -1, -1>;
   // copy because ComplexSchur requires Eigen::Matrix type
-  MatType mv = m;
-  Eigen::ComplexSchur<MatType> cs(mv);
+  Eigen::ComplexSchur<MatType> cs{MatType(m)};
   return std::make_tuple(std::move(cs.matrixU()), std::move(cs.matrixT()));
 }
 
diff --git a/stan/math/prim/fun/eigendecompose_sym.hpp b/stan/math/prim/fun/eigendecompose_sym.hpp
@@ -27,9 +27,9 @@ eigendecompose_sym(const EigMat& m) {
     return std::make_tuple(Eigen::Matrix<value_type_t<EigMat>, -1, -1>(0, 0),
                            Eigen::Matrix<value_type_t<EigMat>, -1, 1>(0, 1));
   }
-  check_symmetric("eigendecompose_sym", "m", m);
   using PlainMat = plain_type_t<EigMat>;
   const PlainMat& m_eval = m;
+  check_symmetric("eigendecompose_sym", "m", m_eval);
 
   Eigen::SelfAdjointEigenSolver<PlainMat> solver(m_eval);
   return std::make_tuple(std::move(solver.eigenvectors()),
diff --git a/stan/math/prim/fun/eigenvalues_sym.hpp b/stan/math/prim/fun/eigenvalues_sym.hpp
@@ -24,9 +24,9 @@ Eigen::Matrix<value_type_t<EigMat>, -1, 1> eigenvalues_sym(const EigMat& m) {
   if (unlikely(m.size() == 0)) {
     return Eigen::Matrix<value_type_t<EigMat>, -1, 1>(0, 1);
   }
-  check_symmetric("eigenvalues_sym", "m", m);
   using PlainMat = plain_type_t<EigMat>;
   const PlainMat& m_eval = m;
+  check_symmetric("eigenvalues_sym", "m", m_eval);
 
   Eigen::SelfAdjointEigenSolver<PlainMat> solver(m_eval,
                                                  Eigen::EigenvaluesOnly);
diff --git a/stan/math/prim/fun/eigenvectors_sym.hpp b/stan/math/prim/fun/eigenvectors_sym.hpp
@@ -15,10 +15,10 @@ eigenvectors_sym(const EigMat& m) {
   if (unlikely(m.size() == 0)) {
     return Eigen::Matrix<value_type_t<EigMat>, -1, -1>(0, 0);
   }
-  check_symmetric("eigenvalues_sym", "m", m);
-
   using PlainMat = plain_type_t<EigMat>;
   const PlainMat& m_eval = m;
+  check_symmetric("eigenvalues_sym", "m", m_eval);
+
 
   Eigen::SelfAdjointEigenSolver<PlainMat> solver(m_eval);
   return solver.eigenvectors();
diff --git a/stan/math/prim/fun/singular_values.hpp b/stan/math/prim/fun/singular_values.hpp
@@ -21,8 +21,9 @@ namespace math {
 template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
           require_not_st_var<EigMat>* = nullptr>
 auto singular_values(const EigMat& m) {
-  check_nonzero_size("singular_values", "m", m);
-
+  if (unlikely(m.size() == 0)) {
+    return Eigen::Matrix<base_type_t<EigMat>, Eigen::Dynamic, 1>(0, 1);
+  }
   return Eigen::JacobiSVD<Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic,
                                         Eigen::Dynamic> >(m)
       .singularValues();
diff --git a/stan/math/prim/fun/svd.hpp b/stan/math/prim/fun/svd.hpp
@@ -22,7 +22,11 @@ std::tuple<Eigen::Matrix<value_type_t<EigMat>, -1, -1>,
            Eigen::Matrix<base_type_t<EigMat>, -1, 1>,
            Eigen::Matrix<value_type_t<EigMat>, -1, -1>>
 svd(const EigMat& m) {
-  check_nonzero_size("svd", "m", m);
+  if (unlikely(m.size() == 0)) {
+    return std::make_tuple(Eigen::Matrix<value_type_t<EigMat>, -1, -1>(0, 0),
+                           Eigen::Matrix<base_type_t<EigMat>, -1, 1>(0, 1),
+                           Eigen::Matrix<value_type_t<EigMat>, -1, -1>(0, 0));
+  }
 
   Eigen::JacobiSVD<Eigen::Matrix<value_type_t<EigMat>, -1, -1>> svd(
       m, Eigen::ComputeThinU | Eigen::ComputeThinV);
diff --git a/stan/math/prim/fun/svd_U.hpp b/stan/math/prim/fun/svd_U.hpp
@@ -18,12 +18,12 @@ template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
           require_not_st_var<EigMat>* = nullptr>
 Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic> svd_U(
     const EigMat& m) {
-  check_nonzero_size("svd_U", "m", m);
-
-  return Eigen::JacobiSVD<Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic,
-                                        Eigen::Dynamic> >(m,
-                                                          Eigen::ComputeThinU)
-      .matrixU();
+  using MatType
+      = Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic>;
+  if (unlikely(m.size() == 0)) {
+    return MatType(0, 0);
+  }
+  return Eigen::JacobiSVD<MatType> (m, Eigen::ComputeThinU).matrixU();
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/svd_V.hpp b/stan/math/prim/fun/svd_V.hpp
@@ -18,12 +18,12 @@ template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
           require_not_st_var<EigMat>* = nullptr>
 Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic> svd_V(
     const EigMat& m) {
-  check_nonzero_size("svd_V", "m", m);
-
-  return Eigen::JacobiSVD<Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic,
-                                        Eigen::Dynamic> >(m,
-                                                          Eigen::ComputeThinV)
-      .matrixV();
+  using MatType
+      = Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic>;
+  if (unlikely(m.size() == 0)) {
+    return MatType(0, 0);
+  }
+  return Eigen::JacobiSVD<MatType>(m, Eigen::ComputeThinV).matrixV();
 }
 
 }  // namespace math
diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp
@@ -22,7 +22,9 @@ namespace math {
 template <typename EigMat, require_rev_matrix_t<EigMat>* = nullptr>
 inline auto singular_values(const EigMat& m) {
   using ret_type = return_var_matrix_t<Eigen::VectorXd, EigMat>;
-  check_nonzero_size("singular_values", "m", m);
+  if (unlikely(m.size() == 0)) {
+    return ret_type(Eigen::VectorXd(0));
+  }
 
   auto arena_m = to_arena(m);
 
diff --git a/stan/math/rev/fun/svd.hpp b/stan/math/rev/fun/svd.hpp
@@ -28,7 +28,11 @@ inline auto svd(const EigMat& m) {
   using mat_ret_type = return_var_matrix_t<Eigen::MatrixXd, EigMat>;
   using vec_ret_type = return_var_matrix_t<Eigen::VectorXd, EigMat>;
 
-  check_nonzero_size("svd", "m", m);
+  if (unlikely(m.size() == 0)) {
+    return std::make_tuple(mat_ret_type(Eigen::MatrixXd(0, 0)),
+                           vec_ret_type(Eigen::VectorXd(0, 1)),
+                           mat_ret_type(Eigen::MatrixXd(0, 0)));
+  }
 
   const int M = std::min(m.rows(), m.cols());
   auto arena_m = to_arena(m);
diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp
@@ -24,7 +24,9 @@ namespace math {
 template <typename EigMat, require_rev_matrix_t<EigMat>* = nullptr>
 inline auto svd_U(const EigMat& m) {
   using ret_type = return_var_matrix_t<Eigen::MatrixXd, EigMat>;
-  check_nonzero_size("svd_U", "m", m);
+  if (unlikely(m.size() == 0)) {
+    return ret_type(Eigen::MatrixXd(0, 0));
+  }
 
   const int M = std::min(m.rows(), m.cols());
   auto arena_m = to_arena(m);
diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp
@@ -24,7 +24,9 @@ namespace math {
 template <typename EigMat, require_rev_matrix_t<EigMat>* = nullptr>
 inline auto svd_V(const EigMat& m) {
   using ret_type = return_var_matrix_t<Eigen::MatrixXd, EigMat>;
-  check_nonzero_size("svd_V", "m", m);
+  if (unlikely(m.size() == 0)) {
+    return ret_type(Eigen::MatrixXd(0, 0));
+  }
 
   const int M = std::min(m.rows(), m.cols());
   auto arena_m = to_arena(m);
diff --git a/test/unit/math/mix/fun/singular_values_test.cpp b/test/unit/math/mix/fun/singular_values_test.cpp
@@ -5,7 +5,7 @@ TEST(MathMixMatFun, singularValues) {
   auto f = [](const auto& x) { return stan::math::singular_values(x); };
 
   Eigen::MatrixXd m00(0, 0);
-  EXPECT_THROW(f(m00), std::invalid_argument);
+  EXPECT_NO_THROW(f(m00));
 
   Eigen::MatrixXd m11(1, 1);
   m11 << 1.1;
diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp
@@ -5,7 +5,7 @@ TEST(MathMixMatFun, svd_U) {
   auto f = [](const auto& x) { return stan::math::svd_U(x); };
 
   Eigen::MatrixXd m00(0, 0);
-  EXPECT_THROW(f(m00), std::invalid_argument);
+  EXPECT_NO_THROW(f(m00));
 
   Eigen::MatrixXd m11(1, 1);
   m11 << 1.1;
diff --git a/test/unit/math/mix/fun/svd_V_test.cpp b/test/unit/math/mix/fun/svd_V_test.cpp
@@ -5,7 +5,7 @@ TEST(MathMixMatFun, svd_V) {
   auto f = [](const auto& x) { return stan::math::svd_V(x); };
 
   Eigen::MatrixXd m00(0, 0);
-  EXPECT_THROW(f(m00), std::invalid_argument);
+  EXPECT_NO_THROW(f(m00));
 
   Eigen::MatrixXd m11(1, 1);
   m11 << 1.1;
diff --git a/test/unit/math/mix/fun/svd_test.cpp b/test/unit/math/mix/fun/svd_test.cpp
@@ -7,9 +7,9 @@ TEST(MathMixMatFun, svd) {
   auto h = [](const auto& x) { return std::get<2>(stan::math::svd(x)); };
 
   Eigen::MatrixXd m00(0, 0);
-  EXPECT_THROW(f(m00), std::invalid_argument);
-  EXPECT_THROW(g(m00), std::invalid_argument);
-  EXPECT_THROW(h(m00), std::invalid_argument);
+  EXPECT_NO_THROW(f(m00));
+  EXPECT_NO_THROW(g(m00));
+  EXPECT_NO_THROW(h(m00));
 
   Eigen::MatrixXd m11(1, 1);
   m11 << 1.1;
diff --git a/test/unit/math/prim/fun/singular_values_test.cpp b/test/unit/math/prim/fun/singular_values_test.cpp
@@ -13,7 +13,7 @@ TEST(MathMatrixPrimMat, singular_values) {
   using matrix_c = Eigen::Matrix<compl_t, Eigen::Dynamic, Eigen::Dynamic>;
 
   matrix_d m0(0, 0);
-  EXPECT_THROW(singular_values(m0), std::invalid_argument);
+  EXPECT_NO_THROW(singular_values(m0));
 
   matrix_d m1(1, 1);
   m1 << 1.0;
diff --git a/test/unit/math/prim/fun/svd_U_test.cpp b/test/unit/math/prim/fun/svd_U_test.cpp
@@ -14,7 +14,7 @@ TEST(MathMatrixPrimMat, svd_U) {
   // Values generated using R base::svd
 
   matrix_d m00(0, 0);
-  EXPECT_THROW(svd_U(m00), std::invalid_argument);
+  EXPECT_NO_THROW(svd_U(m00));
 
   matrix_d m11(1, 1);
   m11 << 5;
diff --git a/test/unit/math/prim/fun/svd_V_test.cpp b/test/unit/math/prim/fun/svd_V_test.cpp
@@ -14,7 +14,7 @@ TEST(MathMatrixPrimMat, svd_V) {
   // Values generated using R base::svd
 
   matrix_d m00(0, 0);
-  EXPECT_THROW(svd_V(m00), std::invalid_argument);
+  EXPECT_NO_THROW(svd_V(m00));
 
   matrix_d m11(1, 1);
   m11 << 5;
diff --git a/test/unit/math/prim/fun/svd_test.cpp b/test/unit/math/prim/fun/svd_test.cpp
@@ -16,7 +16,7 @@ TEST(MathMatrixPrimMat, svd) {
 
   matrix_d m00(0, 0);
 
-  EXPECT_THROW(svd(m00), std::invalid_argument);
+  EXPECT_NO_THROW(svd(m00));
 
   matrix_d m11(1, 1);
   m11 << 5;
