eigenvalues: GeneralizedEigenSolver

Author: Lucas-Haubert

CMakeLists.txt

@@ -236,6 +236,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
     ${${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS}
     include/eigenpy/decompositions/decompositions.hpp
     include/eigenpy/decompositions/EigenSolver.hpp
+    include/eigenpy/decompositions/GeneralizedEigenSolver.hpp
     include/eigenpy/decompositions/ComplexEigenSolver.hpp
     include/eigenpy/decompositions/ComplexSchur.hpp
     include/eigenpy/decompositions/PermutationMatrix.hpp
@@ -320,6 +321,7 @@ set(${PROJECT_NAME}_SOLVERS_SOURCES src/solvers/preconditioners.cpp
 set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
     src/decompositions/decompositions.cpp
     src/decompositions/eigen-solver.cpp
+    src/decompositions/generalized-eigen-solver.cpp
     src/decompositions/complex-eigen-solver.cpp
     src/decompositions/complex-schur.cpp
     src/decompositions/llt-solver.cpp

include/eigenpy/decompositions/GeneralizedEigenSolver.hpp

@@ -0,0 +1,85 @@
+/*
+ * Copyright 2020 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_generalized_eigen_solver_hpp__
+#define __eigenpy_decompositions_generalized_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"
+
+namespace eigenpy {
+
+template <typename _MatrixType>
+struct GeneralizedEigenSolverVisitor
+    : 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<MatrixType, MatrixType, bp::optional<bool>>(
+            bp::args("A", "B", "computeEigenVectors"), 
+            "Computes Schur  of given matrix"))
+
+        .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.",
+             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<>())
+
+     //    .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("info", &Solver::info, bp::arg("self"),
+             "NumericalIssue if the input contains INF or NaN values or "
+             "overflow occured. Returns Success otherwise.")
+
+        .def("setMaxIterations", &Solver::setMaxIterations,
+             bp::args("self", "max_iter"),
+             "Sets the maximum number of iterations allowed.",
+             bp::return_self<>())
+     ;
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "GeneralizedEigenSolver" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string& name) {
+    bp::class_<Solver>(name.c_str(), bp::no_init)
+        .def(GeneralizedEigenSolverVisitor())
+        .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_eigen_solver_hpp__

src/decompositions/decompositions.cpp

@@ -9,6 +9,7 @@
 namespace eigenpy {
 
 void exposeEigenSolver();
+void exposeGeneralizedEigenSolver();
 void exposeComplexEigenSolver();
 void exposeComplexSchur();
 void exposeSelfAdjointEigenSolver();
@@ -30,6 +31,7 @@ void exposeDecompositions() {
   using namespace Eigen;
 
   exposeEigenSolver();
+  exposeGeneralizedEigenSolver();
   exposeComplexEigenSolver();
   exposeComplexSchur();
   exposeSelfAdjointEigenSolver();

src/decompositions/generalized-eigen-solver.cpp

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

unittest/CMakeLists.txt

@@ -137,6 +137,10 @@ add_python_eigenpy_lib_unit_test("py-version" "unittest/python/test_version.py")
 add_python_eigenpy_lib_unit_test("py-eigen-solver"
                                  "unittest/python/test_eigen_solver.py")
 
+add_python_eigenpy_lib_unit_test(
+  "py-generalized-eigen-solver"
+  "unittest/python/test_generalized-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

@@ -0,0 +1,25 @@
+import numpy as np
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+A = rng.random((dim, dim))
+B = rng.random((dim, dim))
+
+ges = eigenpy.GeneralizedEigenSolver(A, B)
+
+alphas = ges.alphas()
+betas = ges.betas()
+V = ges.eigenvectors()
+
+eigenvalues = alphas / betas
+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)