eigenvalues: HenssenbergDecomposition

Author: Lucas-Haubert

CMakeLists.txt

@@ -238,6 +238,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
     include/eigenpy/decompositions/EigenSolver.hpp
     include/eigenpy/decompositions/GeneralizedEigenSolver.hpp
     include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp
+    include/eigenpy/decompositions/HessenbergDecomposition.hpp
     include/eigenpy/decompositions/ComplexEigenSolver.hpp
     include/eigenpy/decompositions/ComplexSchur.hpp
     include/eigenpy/decompositions/PermutationMatrix.hpp
@@ -331,6 +332,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
     src/decompositions/bdcsvd-solver.cpp
     src/decompositions/jacobisvd-solver.cpp
     src/decompositions/fullpivlu-solver.cpp
+    src/decompositions/hessenberg-decomposition.cpp
     src/decompositions/partialpivlu-solver.cpp
     src/decompositions/minres-solver.cpp
     src/decompositions/sparse-lu-solver.cpp

include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp

@@ -17,30 +17,35 @@ namespace eigenpy {
 
 template <typename _MatrixType>
 struct GeneralizedSelfAdjointEigenSolverVisitor
-    : public boost::python::def_visitor<GeneralizedSelfAdjointEigenSolverVisitor<_MatrixType>> {
+    : 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<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. "))
+            bp::args("matA", "matB", "options"),
+            "Constructor: Computes generalized eigendecomposition of given "
+            "matrix pencil. "))
 
-        .def("compute", &GeneralizedSelfAdjointEigenSolverVisitor::compute_proxy<MatrixType>,
-             bp::args("self", "A", "B"), 
+        .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<>());
+             "Computes generalized eigendecomposition of given matrix pencil.",
+             bp::return_self<>());
   }
 
   static void expose() {
@@ -50,18 +55,21 @@ struct GeneralizedSelfAdjointEigenSolverVisitor
   }
 
   static void expose(const std::string& name) {
-    bp::class_<Solver, bp::bases<Eigen::SelfAdjointEigenSolver<MatrixType>>>(name.c_str(), bp::no_init)
+    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) {
+  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__
+#endif  // ifndef
+        // __eigenpy_decompositions_generalized_self_adjoint_eigen_solver_hpp__

include/eigenpy/decompositions/HessenbergDecomposition.hpp

@@ -0,0 +1,78 @@
+/*
+ * Copyright 2020 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_hessenberg_decomposition_hpp__
+#define __eigenpy_decompositions_hessenberg_decomposition_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 HessenbergDecompositionVisitor
+    : public boost::python::def_visitor<
+          HessenbergDecompositionVisitor<_MatrixType>> {
+  typedef _MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef Eigen::HessenbergDecomposition<MatrixType> Solver;
+
+  template <class PyClass>
+  void visit(PyClass& cl) const {
+    cl.def(bp::init<Eigen::DenseIndex>(
+            bp::arg("size"), "Default constructor; the decomposition will be computed later. "))
+        .def(bp::init<MatrixType>(bp::arg("matrix"),
+            "Constructor; computes Hessenberg decomposition of given matrix. "))
+
+        .def("compute",
+             &HessenbergDecompositionVisitor::compute_proxy<MatrixType>,
+             bp::args("self", "A"),
+             "Computes Hessenberg decomposition of given matrix. ",
+             bp::return_self<>())
+        .def("compute",
+             (Solver &
+              (Solver::*)(const Eigen::EigenBase<MatrixType>& matrix)) &
+                 Solver::compute,
+             bp::args("self", "A"),
+             "Computes Hessenberg decomposition of given matrix. ",
+             bp::return_self<>())
+
+        .def("householderCoefficients", &Solver::householderCoefficients, bp::arg("self"),
+             "Returns the Householder coefficients. ",
+             bp::return_value_policy<bp::copy_const_reference>())
+
+        // matrixH
+        // matrixQ
+
+        .def("packedMatrix", &Solver::packedMatrix, bp::arg("self"),
+             "Returns the internal representation of the decomposition. ",
+             bp::return_value_policy<bp::copy_const_reference>());
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "HessenbergDecomposition" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string& name) {
+    bp::class_<Solver>(name.c_str(), bp::no_init)
+        .def(HessenbergDecompositionVisitor())
+        .def(IdVisitor<Solver>());
+  }
+
+ private:
+  template <typename MatrixType>
+  static Solver& compute_proxy(Solver& self, const Eigen::EigenBase<MatrixType>& matrix) {
+    return self.compute(matrix);
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decompositions_hessenberg_decomposition_hpp__

src/decompositions/decompositions.cpp

@@ -12,6 +12,7 @@ void exposeEigenSolver();
 void exposeGeneralizedEigenSolver();
 void exposeSelfAdjointEigenSolver();
 void exposeGeneralizedSelfAdjointEigenSolver();
+void exposeHessenbergDecomposition();
 void exposeComplexEigenSolver();
 void exposeComplexSchur();
 void exposeLLTSolver();
@@ -35,6 +36,7 @@ void exposeDecompositions() {
   exposeGeneralizedEigenSolver();
   exposeSelfAdjointEigenSolver();
   exposeGeneralizedSelfAdjointEigenSolver();
+  exposeHessenbergDecomposition();
   exposeComplexEigenSolver();
   exposeComplexSchur();
   exposeLLTSolver();

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

@@ -8,6 +8,7 @@
 namespace eigenpy {
 void exposeGeneralizedSelfAdjointEigenSolver() {
   using namespace Eigen;
-  GeneralizedSelfAdjointEigenSolverVisitor<MatrixXd>::expose("GeneralizedSelfAdjointEigenSolver");
+  GeneralizedSelfAdjointEigenSolverVisitor<MatrixXd>::expose(
+      "GeneralizedSelfAdjointEigenSolver");
 }
 }  // namespace eigenpy

src/decompositions/hessenberg-decomposition.cpp

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

unittest/CMakeLists.txt

@@ -155,6 +155,10 @@ add_python_eigenpy_lib_unit_test(
   "py-self-adjoint-eigen-solver"
   "unittest/python/test_self_adjoint_eigen_solver.py")
 
+add_python_eigenpy_lib_unit_test(
+  "py-hessenberg-decomposition"
+  "unittest/python/test_hessenberg_decomposition.py")
+
 add_python_eigenpy_lib_unit_test("py-LLT" "unittest/python/test_LLT.py")
 
 add_python_eigenpy_lib_unit_test("py-LDLT" "unittest/python/test_LDLT.py")

unittest/python/test_generalized_eigen_solver.py

@@ -14,12 +14,13 @@
 V = ges.eigenvectors()
 
 eigenvalues = alphas / betas
-for i in range(dim):
-    v = V[:, i]
-    lam = eigenvalues[i]
+if np.all(np.abs(betas) >= 1e-15):
+    for i in range(dim):
+        v = V[:, i]
+        lam = eigenvalues[i]
 
-    Av = A @ v
-    lam_Bv = lam * (B @ v)
+        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)
+        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

@@ -19,7 +19,7 @@
 for i in range(dim):
     v = V[:, i]
     lam = D[i]
-    
+
     Av = A @ v
     lam_Bv = lam * (B @ v)
 

unittest/python/test_hessenberg_decomposition.py

@@ -0,0 +1,9 @@
+import numpy as np
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+A = rng.random((dim, dim))
+
+ges = eigenpy.HessenbergDecomposition(A)