eigenvalues: GeneralizedSelfAdjointSolver

Author: Lucas-Haubert

CMakeLists.txt

@@ -237,6 +237,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
     include/eigenpy/decompositions/decompositions.hpp
     include/eigenpy/decompositions/EigenSolver.hpp
     include/eigenpy/decompositions/GeneralizedEigenSolver.hpp
+    include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp
     include/eigenpy/decompositions/ComplexEigenSolver.hpp
     include/eigenpy/decompositions/ComplexSchur.hpp
     include/eigenpy/decompositions/PermutationMatrix.hpp
@@ -322,6 +323,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
     src/decompositions/decompositions.cpp
     src/decompositions/eigen-solver.cpp
     src/decompositions/generalized-eigen-solver.cpp
+    src/decompositions/generalized-self-adjoint-eigen-solver.cpp
     src/decompositions/complex-eigen-solver.cpp
     src/decompositions/complex-schur.cpp
     src/decompositions/llt-solver.cpp

include/eigenpy/decompositions/GeneralizedEigenSolver.hpp

@@ -16,39 +16,42 @@ namespace eigenpy {
 
 template <typename _MatrixType>
 struct GeneralizedEigenSolverVisitor
-    : public boost::python::def_visitor<GeneralizedEigenSolverVisitor<_MatrixType>> {
+    : public boost::python::def_visitor<
+          GeneralizedEigenSolverVisitor<_MatrixType>> {
   typedef _MatrixType MatrixType;
   typedef typename MatrixType::Scalar Scalar;
   typedef Eigen::GeneralizedEigenSolver<MatrixType> Solver;
 
   template <class PyClass>
   void visit(PyClass& cl) const {
     cl.def(bp::init<>("Default constructor"))
-        .def(bp::init<Eigen::DenseIndex>(bp::arg("size"), 
-        "Default constructor with memory preallocation. "))
+        .def(bp::init<Eigen::DenseIndex>(
+            bp::arg("size"), "Default constructor with memory preallocation. "))
         .def(bp::init<MatrixType, MatrixType, bp::optional<bool>>(
-            bp::args("A", "B", "computeEigenVectors"), 
-            "Computes Schur  of given matrix"))
+            bp::args("A", "B", "computeEigenVectors"),
+            "Computes the generalized eigendecomposition of given matrix "
+            "pair. "))
 
         .def("alphas", &Solver::alphas, bp::arg("self"),
              "Returns the vectors containing the alpha values ")
         .def("betas", &Solver::betas, bp::arg("self"),
              "Returns the vectors containing the beta values ")
 
-        .def("compute", &GeneralizedEigenSolverVisitor::compute_proxy<MatrixType>,
-             bp::args("self", "A", "B"), "Computes the Schur of given matrix.",
+        .def("compute",
+             &GeneralizedEigenSolverVisitor::compute_proxy<MatrixType>,
+             bp::args("self", "A", "B"),
+             "Computes generalized eigendecomposition of given matrix. ",
              bp::return_self<>())
         .def("compute",
              (Solver &
               (Solver::*)(const MatrixType& A, const MatrixType& B, bool)) &
                  Solver::compute,
              bp::args("self", "A", "B", "computeEigenvectors"),
-             "Computes generalized eigendecomposition of given matrix. .", bp::return_self<>())
+             "Computes generalized eigendecomposition of given matrix. .",
+             bp::return_self<>())
 
-     //    .def("eigenvalues", &Solver::eigenvalues, bp::arg("self"),
-     //         "Returns an expression of the computed generalized eigenvalues. ")
-     //    .def("eigenvectors", &Solver::eigenvectors, bp::arg("self"),
-     //         "Returns an expression of the computed generalized eigenvectors. ")
+        .def("eigenvectors", &Solver::eigenvectors, bp::arg("self"),
+             "Returns an expression of the computed generalized eigenvectors. ")
 
         .def("info", &Solver::info, bp::arg("self"),
              "NumericalIssue if the input contains INF or NaN values or "
@@ -57,8 +60,7 @@ struct GeneralizedEigenSolverVisitor
         .def("setMaxIterations", &Solver::setMaxIterations,
              bp::args("self", "max_iter"),
              "Sets the maximum number of iterations allowed.",
-             bp::return_self<>())
-     ;
+             bp::return_self<>());
   }
 
   static void expose() {
@@ -75,7 +77,8 @@ struct GeneralizedEigenSolverVisitor
 
  private:
   template <typename MatrixType>
-  static Solver& compute_proxy(Solver& self, const MatrixType& A, const MatrixType& B) {
+  static Solver& compute_proxy(Solver& self, const MatrixType& A,
+                               const MatrixType& B) {
     return self.compute(A, B);
   }
 };

include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp

@@ -0,0 +1,67 @@
+/*
+ * Copyright 2020 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_generalized_self_adjoint_eigen_solver_hpp__
+#define __eigenpy_decompositions_generalized_self_adjoint_eigen_solver_hpp__
+
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
+
+#include "eigenpy/eigen-to-python.hpp"
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+#include "eigenpy/decompositions/SelfAdjointEigenSolver.hpp"
+
+namespace eigenpy {
+
+template <typename _MatrixType>
+struct GeneralizedSelfAdjointEigenSolverVisitor
+    : public boost::python::def_visitor<GeneralizedSelfAdjointEigenSolverVisitor<_MatrixType>> {
+  typedef _MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef Eigen::GeneralizedSelfAdjointEigenSolver<MatrixType> Solver;
+
+  template <class PyClass>
+  void visit(PyClass& cl) const {
+    cl.def(bp::init<>("Default constructor"))
+        .def(bp::init<Eigen::DenseIndex>(bp::arg("size"), 
+        "Default constructor with memory preallocation. "))
+        .def(bp::init<MatrixType, MatrixType, bp::optional<int>>(
+            bp::args("matA", "matB", "options"), 
+            "Constructor: Computes generalized eigendecomposition of given matrix pencil. "))
+
+        .def("compute", &GeneralizedSelfAdjointEigenSolverVisitor::compute_proxy<MatrixType>,
+             bp::args("self", "A", "B"), 
+             "Computes generalized eigendecomposition of given matrix pencil. ",
+             bp::return_self<>())
+        .def("compute",
+             (Solver &
+              (Solver::*)(const MatrixType& A, const MatrixType& B, int)) &
+                 Solver::compute,
+             bp::args("self", "A", "B", "options"),
+             "Computes generalized eigendecomposition of given matrix pencil.", bp::return_self<>());
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "GeneralizedSelfAdjointEigenSolver" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string& name) {
+    bp::class_<Solver, bp::bases<Eigen::SelfAdjointEigenSolver<MatrixType>>>(name.c_str(), bp::no_init)
+        .def(GeneralizedSelfAdjointEigenSolverVisitor())
+        .def(IdVisitor<Solver>());
+  }
+
+ private:
+  template <typename MatrixType>
+  static Solver& compute_proxy(Solver& self, const MatrixType& A, const MatrixType& B) {
+    return self.compute(A, B);
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decompositions_generalized_self_adjoint_eigen_solver_hpp__

src/decompositions/decompositions.cpp

@@ -10,9 +10,10 @@ namespace eigenpy {
 
 void exposeEigenSolver();
 void exposeGeneralizedEigenSolver();
+void exposeSelfAdjointEigenSolver();
+void exposeGeneralizedSelfAdjointEigenSolver();
 void exposeComplexEigenSolver();
 void exposeComplexSchur();
-void exposeSelfAdjointEigenSolver();
 void exposeLLTSolver();
 void exposeLDLTSolver();
 void exposeFullPivLUSolver();
@@ -32,9 +33,10 @@ void exposeDecompositions() {
 
   exposeEigenSolver();
   exposeGeneralizedEigenSolver();
+  exposeSelfAdjointEigenSolver();
+  exposeGeneralizedSelfAdjointEigenSolver();
   exposeComplexEigenSolver();
   exposeComplexSchur();
-  exposeSelfAdjointEigenSolver();
   exposeLLTSolver();
   exposeLDLTSolver();
   exposeFullPivLUSolver();

src/decompositions/generalized-self-adjoint-eigen-solver.cpp

@@ -0,0 +1,13 @@
+
+/*
+ * Copyright 2024 INRIA
+ */
+
+#include "eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp"
+
+namespace eigenpy {
+void exposeGeneralizedSelfAdjointEigenSolver() {
+  using namespace Eigen;
+  GeneralizedSelfAdjointEigenSolverVisitor<MatrixXd>::expose("GeneralizedSelfAdjointEigenSolver");
+}
+}  // namespace eigenpy

unittest/CMakeLists.txt

@@ -139,7 +139,11 @@ add_python_eigenpy_lib_unit_test("py-eigen-solver"
 
 add_python_eigenpy_lib_unit_test(
   "py-generalized-eigen-solver"
-  "unittest/python/test_generalized-eigen_solver.py")
+  "unittest/python/test_generalized_eigen_solver.py")
+
+add_python_eigenpy_lib_unit_test(
+  "py-generalized-self-adjoint-eigen-solver"
+  "unittest/python/test_generalized_self_adjoint_eigen_solver.py")
 
 add_python_eigenpy_lib_unit_test("py-complex-eigen-solver"
                                  "unittest/python/test_complex_eigen_solver.py")

unittest/python/test_generalized_eigen_solver.py

@@ -17,9 +17,9 @@
 for i in range(dim):
     v = V[:, i]
     lam = eigenvalues[i]
-    
+
     Av = A @ v
     lam_Bv = lam * (B @ v)
-    
+
     assert eigenpy.is_approx(Av.real, lam_Bv.real)
     assert eigenpy.is_approx(Av.imag, lam_Bv.imag)

unittest/python/test_generalized_self_adjoint_eigen_solver.py

@@ -0,0 +1,26 @@
+import numpy as np
+
+import eigenpy
+
+dim = 5
+rng = np.random.default_rng()
+
+A = rng.random((dim, dim))
+A = (A + A.T) * 0.5
+
+B = rng.random((dim, dim))
+B = B @ B.T + 0.1 * np.eye(dim)
+
+gsaes = eigenpy.GeneralizedSelfAdjointEigenSolver(A, B)
+
+V = gsaes.eigenvectors()
+D = gsaes.eigenvalues()
+
+for i in range(dim):
+    v = V[:, i]
+    lam = D[i]
+    
+    Av = A @ v
+    lam_Bv = lam * (B @ v)
+
+    assert np.allclose(Av, lam_Bv, atol=1e-6)