Merge pull request #2846 from stan-dev/feature/2733-complex-eigendecomposition

Add overloads to eigendecomposition for complex matrices

Author: WardBrian

stan/math/prim/fun/eigenvalues.hpp

@@ -7,14 +7,47 @@
 namespace stan {
 namespace math {
 
-template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr>
-inline auto eigenvalues(const EigMat& m) {
+/**
+ * Return the eigenvalues of a (real-valued) matrix
+ *
+ * @tparam EigMat type of real matrix argument
+ * @param[in] m matrix to find the eigenvectors of. Must be square and have a
+ * non-zero size.
+ * @return a complex-valued column vector with entries the eigenvectors of `m`
+ */
+template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
+          require_not_vt_complex<EigMat>* = nullptr>
+inline Eigen::Matrix<complex_return_t<value_type_t<EigMat>>, 1, -1> eigenvalues(
+    const EigMat& m) {
   using PlainMat = plain_type_t<EigMat>;
   const PlainMat& m_eval = m;
   check_nonzero_size("eigenvalues", "m", m_eval);
   check_square("eigenvalues", "m", m_eval);
 
-  Eigen::EigenSolver<PlainMat> solver(m_eval);
+  Eigen::EigenSolver<PlainMat> solver(m_eval, false);
+  return solver.eigenvalues();
+}
+
+/**
+ * Return the eigenvalues of a (complex-valued) matrix
+ *
+ * @tparam EigCplxMat type of complex matrix argument
+ * @param[in] m matrix to find the eigenvectors of. Must be square and have a
+ * non-zero size.
+ * @return a complex-valued column vector with entries the eigenvectors of `m`
+ */
+template <typename EigCplxMat,
+          require_eigen_matrix_dynamic_t<EigCplxMat>* = nullptr,
+          require_vt_complex<EigCplxMat>* = nullptr>
+inline Eigen::Matrix<complex_return_t<value_type_t<EigCplxMat>>, 1, -1>
+eigenvalues(const EigCplxMat& m) {
+  using PlainMat = Eigen::Matrix<scalar_type_t<EigCplxMat>, -1, -1>;
+  const PlainMat& m_eval = m;
+  check_nonzero_size("eigenvalues", "m", m_eval);
+  check_square("eigenvalues", "m", m_eval);
+
+  Eigen::ComplexEigenSolver<PlainMat> solver(m_eval, false);
+
   return solver.eigenvalues();
 }
 

stan/math/prim/fun/eigenvectors.hpp

@@ -7,8 +7,19 @@
 namespace stan {
 namespace math {
 
-template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr>
-inline auto eigenvectors(const EigMat& m) {
+/**
+ * Return the eigenvectors of a (real-valued) matrix
+ *
+ * @tparam EigMat type of real matrix argument
+ * @param[in] m matrix to find the eigenvectors of. Must be square and have a
+ * non-zero size.
+ * @return a complex-valued matrix with columns representing the eigenvectors of
+ * `m`
+ */
+template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
+          require_not_vt_complex<EigMat>* = nullptr>
+inline Eigen::Matrix<complex_return_t<value_type_t<EigMat>>, -1, -1>
+eigenvectors(const EigMat& m) {
   using PlainMat = plain_type_t<EigMat>;
   const PlainMat& m_eval = m;
   check_nonzero_size("eigenvectors", "m", m_eval);
@@ -18,6 +29,29 @@ inline auto eigenvectors(const EigMat& m) {
   return solver.eigenvectors();
 }
 
+/**
+ * Return the eigenvectors of a (complex-valued) matrix
+ *
+ * @tparam EigCplxMat type of complex matrix argument
+ * @param[in] m matrix to find the eigenvectors of. Must be square and have a
+ * non-zero size.
+ * @return a complex-valued matrix with columns representing the eigenvectors of
+ * `m`
+ */
+template <typename EigCplxMat,
+          require_eigen_matrix_dynamic_t<EigCplxMat>* = nullptr,
+          require_vt_complex<EigCplxMat>* = nullptr>
+inline Eigen::Matrix<complex_return_t<value_type_t<EigCplxMat>>, -1, -1>
+eigenvectors(const EigCplxMat& m) {
+  using PlainMat = Eigen::Matrix<scalar_type_t<EigCplxMat>, -1, -1>;
+  const PlainMat& m_eval = m;
+  check_nonzero_size("eigenvectors", "m", m_eval);
+  check_square("eigenvectors", "m", m_eval);
+
+  Eigen::ComplexEigenSolver<PlainMat> solver(m_eval);
+  return solver.eigenvectors();
+}
+
 }  // namespace math
 }  // namespace stan
 #endif

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

@@ -15,5 +15,19 @@ TEST(mathMixFun, eigenvalues) {
   EXPECT_THROW(f(a32), std::invalid_argument);
 }
 
+TEST(mathMixFun, eigenvaluesComplex) {
+  auto f = [](const auto& x) {
+    using stan::math::eigenvalues;
+    return eigenvalues(stan::math::to_complex(x, 0));
+  };
+  for (const auto& x : stan::test::square_test_matrices(0, 2)) {
+    stan::test::expect_ad(f, x);
+  }
+
+  Eigen::MatrixXd a32(3, 2);
+  a32 << 3, -5, 7, -7.2, 9.1, -6.3;
+  EXPECT_THROW(f(a32), std::invalid_argument);
+}
+
 // see eigenvectors_test.cpp for test of eigenvectors() and eigenvalues()
 // using reconstruction identities

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

@@ -15,19 +15,33 @@ TEST(mathMixFun, eigenvectors) {
   EXPECT_THROW(f(a32), std::invalid_argument);
 }
 
+TEST(mathMixFun, eigenvectorsComplex) {
+  auto f = [](const auto& x) {
+    using stan::math::eigenvectors;
+    return eigenvectors(stan::math::to_complex(x, 0));
+  };
+  for (const auto& x : stan::test::square_test_matrices(0, 2)) {
+    stan::test::expect_ad(f, x);
+  }
+
+  Eigen::MatrixXd a32(3, 2);
+  a32 << 3, -5, 7, -7.2, 9.1, -6.3;
+  EXPECT_THROW(f(a32), std::invalid_argument);
+}
+
 template <typename T>
 void expect_identity_matrix(const T& x) {
-  EXPECT_EQUAL(x.rows(), x.cols());
+  EXPECT_EQ(x.rows(), x.cols());
   for (int j = 0; j < x.cols(); ++j) {
     for (int i = 0; i < x.rows(); ++i) {
-      EXPECT_NEAR(i == j ? 1 : 0, x(i, j), 1e-6);
+      EXPECT_NEAR(i == j ? 1 : 0, stan::math::value_of_rec(x(i, j)), 1e-6);
     }
   }
 }
 
 template <typename T>
 void expectEigenvectorsId() {
-  for (const auto& m_d : stan::test::square_test_matrices(0, 2)) {
+  for (const auto& m_d : stan::test::square_test_matrices(1, 2)) {
     Eigen::Matrix<T, -1, -1> m(m_d);
     auto vecs = eigenvectors(m).eval();
     auto vals = eigenvalues(m).eval();
@@ -36,7 +50,21 @@ void expectEigenvectorsId() {
   }
 }
 
-// THESE TESTS USED TO WORK STANDALONE
+template <typename T>
+void expectComplexEigenvectorsId() {
+  Eigen::Matrix<std::complex<T>, -1, -1> c22(2, 2);
+  c22 << stan::math::to_complex(T(0), T(-1)),
+      stan::math::to_complex(T(0), T(0)), stan::math::to_complex(T(2), T(0)),
+      stan::math::to_complex(T(4), T(0));
+  auto eigenvalues = stan::math::eigenvalues(c22);
+  auto eigenvectors = stan::math::eigenvectors(c22);
+
+  auto I = (eigenvectors.inverse() * c22 * eigenvectors
+            * eigenvalues.asDiagonal().inverse())
+               .real();
+
+  expect_identity_matrix(I);
+}
 
 TEST(mathMixFun, eigenvectorsId) {
   using d_t = double;
@@ -46,9 +74,15 @@ TEST(mathMixFun, eigenvectorsId) {
   using fv_t = stan::math::fvar<stan::math::var>;
   using ffv_t = stan::math::fvar<fv_t>;
 
-  // expectEigenvectorsId<v_t>();
-  // expectEigenvectorsId<fd_t>();
-  // expectEigenvectorsId<ffd_t>();
-  // expectEigenvectorsId<fv_t>();
-  // expectEigenvectorsId<ffv_t>();
+  expectEigenvectorsId<v_t>();
+  expectEigenvectorsId<fd_t>();
+  expectEigenvectorsId<ffd_t>();
+  expectEigenvectorsId<fv_t>();
+  expectEigenvectorsId<ffv_t>();
+
+  expectComplexEigenvectorsId<v_t>();
+  expectComplexEigenvectorsId<fd_t>();
+  expectComplexEigenvectorsId<ffd_t>();
+  expectComplexEigenvectorsId<fv_t>();
+  expectComplexEigenvectorsId<ffv_t>();
 }