eigenvalues: Tridiagonalization

Author: Lucas-Haubert

CMakeLists.txt

@@ -240,6 +240,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
     include/eigenpy/decompositions/GeneralizedSelfAdjointEigenSolver.hpp
     include/eigenpy/decompositions/HessenbergDecomposition.hpp
     include/eigenpy/decompositions/RealQZ.hpp
+    include/eigenpy/decompositions/Tridiagonalization.hpp
     include/eigenpy/decompositions/RealSchur.hpp
     include/eigenpy/decompositions/ComplexEigenSolver.hpp
     include/eigenpy/decompositions/ComplexSchur.hpp
@@ -336,6 +337,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
     src/decompositions/fullpivlu-solver.cpp
     src/decompositions/hessenberg-decomposition.cpp
     src/decompositions/real-qz.cpp
+    src/decompositions/tridiagonalization.cpp
     src/decompositions/real-schur.cpp
     src/decompositions/partialpivlu-solver.cpp
     src/decompositions/minres-solver.cpp

include/eigenpy/decompositions/RealSchur.hpp

@@ -16,22 +16,20 @@ namespace eigenpy {
 
 template <typename _MatrixType>
 struct RealSchurVisitor
-    : public boost::python::def_visitor<
-          RealSchurVisitor<_MatrixType>> {
+    : public boost::python::def_visitor<RealSchurVisitor<_MatrixType>> {
   typedef _MatrixType MatrixType;
   typedef typename MatrixType::Scalar Scalar;
   typedef Eigen::RealSchur<MatrixType> Solver;
 
   template <class PyClass>
   void visit(PyClass& cl) const {
-    cl.def(bp::init<Eigen::DenseIndex>(
-            bp::arg("size"), "Default constructor. "))
+    cl.def(
+          bp::init<Eigen::DenseIndex>(bp::arg("size"), "Default constructor. "))
         .def(bp::init<MatrixType, bp::optional<bool>>(
             bp::args("matrix", "computeU"),
             "Constructor; computes real Schur decomposition of given matrix. "))
 
-        .def("compute",
-             &RealSchurVisitor::compute_proxy<MatrixType>,
+        .def("compute", &RealSchurVisitor::compute_proxy<MatrixType>,
              bp::args("self", "matrix"),
              "Computes Schur decomposition of given matrix. ",
              bp::return_self<>())
@@ -45,8 +43,7 @@ struct RealSchurVisitor
 
         .def("computeFromHessenberg",
              (Solver & (Solver::*)(const MatrixType& matrixH,
-                                   const MatrixType& matrixQ,
-                                   bool)) &
+                                   const MatrixType& matrixQ, bool)) &
                  Solver::computeFromHessenberg,
              bp::args("self", "matrixH", "matrixQ", "computeU"),
              "Compute Schur decomposition from a given Hessenberg matrix. ",

include/eigenpy/decompositions/Tridiagonalization.hpp

@@ -0,0 +1,83 @@
+/*
+ * Copyright 2020 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_tridiagonalization_hpp__
+#define __eigenpy_decompositions_tridiagonalization_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 TridiagonalizationVisitor
+    : public boost::python::def_visitor<TridiagonalizationVisitor<_MatrixType>> {
+  typedef _MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef Eigen::Tridiagonalization<MatrixType> Solver;
+
+  template <class PyClass>
+  void visit(PyClass& cl) const {
+    cl.def(
+          bp::init<Eigen::DenseIndex>(bp::arg("size"), "Default constructor. "))
+        .def(bp::init<MatrixType>(
+            bp::arg("matrix"),
+            "Constructor; computes tridiagonal decomposition of given matrix. "))
+
+        .def("compute", &TridiagonalizationVisitor::compute_proxy<MatrixType>,
+             bp::args("self", "matrix"),
+             "Computes tridiagonal decomposition of given matrix. ",
+             bp::return_self<>())
+        .def("compute",
+             (Solver &
+              (Solver::*)(const Eigen::EigenBase<MatrixType>& matrix)) &
+                 Solver::compute,
+             bp::args("self", "matrix"),
+             "Computes tridiagonal decomposition of given matrix. ", bp::return_self<>())
+
+        .def("diagonal", &Solver::diagonal, bp::arg("self"),
+             "Returns the diagonal of the tridiagonal matrix T in the decomposition. ")
+
+        .def("householderCoefficients", &Solver::householderCoefficients,
+             bp::arg("self"), "Returns the Householder coefficients. ")
+
+        .def("matrixQ", &Solver::matrixQ,
+             bp::arg("self"), "Returns the unitary matrix Q in the decomposition. ")
+        .def("matrixT", &Solver::matrixT,
+             bp::arg("self"), "Returns the unitary matrix T in the decomposition. ")
+
+        .def("packedMatrix", &Solver::packedMatrix, bp::arg("self"),
+             "Returns the internal representation of the decomposition. ",
+             bp::return_value_policy<bp::copy_const_reference>())
+
+        .def("subDiagonal", &Solver::subDiagonal, bp::arg("self"),
+             "Returns the subdiagonal of the tridiagonal matrix T in the decomposition.");
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "TridiagonalizationVisitor" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string& name) {
+    bp::class_<Solver>(name.c_str(), bp::no_init)
+        .def(TridiagonalizationVisitor())
+        .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_tridiagonalization_hpp__

src/decompositions/decompositions.cpp

@@ -15,6 +15,7 @@ void exposeGeneralizedSelfAdjointEigenSolver();
 void exposeHessenbergDecomposition();
 void exposeRealQZ();
 void exposeRealSchur();
+void exposeTridiagonalization();
 void exposeComplexEigenSolver();
 void exposeComplexSchur();
 void exposeLLTSolver();
@@ -41,6 +42,7 @@ void exposeDecompositions() {
   exposeHessenbergDecomposition();
   exposeRealQZ();
   exposeRealSchur();
+  exposeTridiagonalization();
   exposeComplexEigenSolver();
   exposeComplexSchur();
   exposeLLTSolver();

src/decompositions/tridiagonalization.cpp

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

unittest/CMakeLists.txt

@@ -164,6 +164,9 @@ add_python_eigenpy_lib_unit_test("py-real-qz" "unittest/python/test_real_qz.py")
 add_python_eigenpy_lib_unit_test("py-real-schur"
                                  "unittest/python/test_real_schur.py")
 
+add_python_eigenpy_lib_unit_test("py-tridiagonalization"
+                                 "unittest/python/test_tridiagonalization.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_tridiagonalization.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+A = rng.random((dim, dim))
+
+tri = eigenpy.Tridiagonalization(A)
+
+# Q = tri.matrixQ()
+# print("Q")
+# print(Q)
+# Q_conj = Q.conj().T
+
+# T = tri.matrixT()
+# print("T")
+# print(T)
+
+# assert eigenpy.is_approx(A, Q @ T @ Q_conj)