eigenvalues: ComplexSchur

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/ComplexEigenSolver.hpp
+    include/eigenpy/decompositions/ComplexSchur.hpp
     include/eigenpy/decompositions/PermutationMatrix.hpp
     include/eigenpy/decompositions/LDLT.hpp
     include/eigenpy/decompositions/LLT.hpp
@@ -320,6 +321,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
     src/decompositions/decompositions.cpp
     src/decompositions/eigen-solver.cpp
     src/decompositions/complex-eigen-solver.cpp
+    src/decompositions/complex-schur.cpp
     src/decompositions/llt-solver.cpp
     src/decompositions/ldlt-solver.cpp
     src/decompositions/bdcsvd-solver.cpp

include/eigenpy/decompositions/ComplexSchur.hpp

@@ -0,0 +1,90 @@
+/*
+ * Copyright 2020 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_complex_schur_hpp__
+#define __eigenpy_decompositions_complex_schur_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 ComplexSchurVisitor
+    : public boost::python::def_visitor<ComplexSchurVisitor<_MatrixType>> {
+  typedef _MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef Eigen::ComplexSchur<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::optional<bool>>(
+            bp::args("matrix", "computeU"), "Computes Schur  of given matrix"))
+
+        .def("compute", &ComplexSchurVisitor::compute_proxy<MatrixType>,
+             bp::args("self", "matrix"), "Computes the Schur of given matrix.",
+             bp::return_self<>())
+        .def("compute",
+             (Solver &
+              (Solver::*)(const Eigen::EigenBase<MatrixType>& matrix, bool)) &
+                 Solver::compute,
+             bp::args("self", "matrix", "computeU"),
+             "Computes the Schur of given matrix.", bp::return_self<>())
+
+        .def("computeFromHessenberg",
+             (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType>& matrixH,
+                                   const Eigen::EigenBase<MatrixType>& matrixQ,
+                                   bool)) &
+                 Solver::computeFromHessenberg,
+             bp::args("self", "matrix", "computeU"),
+             "Compute Schur decomposition from a given Hessenberg matrix. ",
+             bp::return_self<>())
+
+        .def("matrixT", &Solver::matrixT, bp::arg("self"),
+             "Returns the triangular matrix in the Schur decomposition. ",
+             bp::return_value_policy<bp::copy_const_reference>())
+        .def("matrixU", &Solver::matrixU, bp::arg("self"),
+             "Returns the unitary matrix in the Schur decomposition. ",
+             bp::return_value_policy<bp::copy_const_reference>())
+
+        .def("info", &Solver::info, bp::arg("self"),
+             "NumericalIssue if the input contains INF or NaN values or "
+             "overflow occured. Returns Success otherwise.")
+
+        .def("getMaxIterations", &Solver::getMaxIterations, bp::arg("self"),
+             "Returns the maximum number of iterations.")
+        .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 =
+        "ComplexSchur" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string& name) {
+    bp::class_<Solver>(name.c_str(), bp::no_init)
+        .def(ComplexSchurVisitor())
+        .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_complex_schur_hpp__

src/decompositions/complex-schur.cpp

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

src/decompositions/decompositions.cpp

@@ -10,6 +10,7 @@ namespace eigenpy {
 
 void exposeEigenSolver();
 void exposeComplexEigenSolver();
+void exposeComplexSchur();
 void exposeSelfAdjointEigenSolver();
 void exposeLLTSolver();
 void exposeLDLTSolver();
@@ -30,6 +31,7 @@ void exposeDecompositions() {
 
   exposeEigenSolver();
   exposeComplexEigenSolver();
+  exposeComplexSchur();
   exposeSelfAdjointEigenSolver();
   exposeLLTSolver();
   exposeLDLTSolver();

unittest/CMakeLists.txt

@@ -140,6 +140,9 @@ add_python_eigenpy_lib_unit_test("py-eigen-solver"
 add_python_eigenpy_lib_unit_test("py-complex-eigen-solver"
                                  "unittest/python/test_complex_eigen_solver.py")
 
+add_python_eigenpy_lib_unit_test("py-complex-schur"
+                                 "unittest/python/test_complex_schur.py")
+
 add_python_eigenpy_lib_unit_test(
   "py-self-adjoint-eigen-solver"
   "unittest/python/test_self_adjoint_eigen_solver.py")

unittest/python/test_complex_eigen_solver.py

@@ -6,10 +6,10 @@
 rng = np.random.default_rng()
 A = rng.random((dim, dim))
 
-es = eigenpy.ComplexEigenSolver(A)
+ces = eigenpy.ComplexEigenSolver(A)
 
-V = es.eigenvectors()
-D = es.eigenvalues()
+V = ces.eigenvectors()
+D = ces.eigenvalues()
 
 assert eigenpy.is_approx(A.dot(V).real, V.dot(np.diag(D)).real)
 assert eigenpy.is_approx(A.dot(V).imag, V.dot(np.diag(D)).imag)

unittest/python/test_complex_schur.py

@@ -0,0 +1,16 @@
+import numpy as np
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+A = rng.random((dim, dim))
+
+cs = eigenpy.ComplexSchur(A)
+
+U = cs.matrixU()
+T = cs.matrixT()
+U_star = U.conj().T
+
+assert eigenpy.is_approx(A.real, (U @ T @ U_star).real)
+assert np.allclose(A.imag, (U @ T @ U_star).imag, atol=1e-10)