Linear algebra: Expose the remaining classes listed in Eigen documentation (decompositions and solvers)

CHANGELOG.md

@@ -8,6 +8,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 ### Added
 
+- Add additional decompositions and solvers from Eigen ([#571](https://github.com/stack-of-tasks/eigenpy/pull/571))
+
+### Added
+
 - Docker images `ghcr.io/stack-of-tasks/eigenpy` ([#575](https://github.com/stack-of-tasks/eigenpy/pull/575))
 
 ### Changed

CMakeLists.txt

@@ -190,10 +190,14 @@ set(${PROJECT_NAME}_SOLVERS_HEADERS
     include/eigenpy/solvers/preconditioners.hpp
     include/eigenpy/solvers/IterativeSolverBase.hpp
     include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp
+    include/eigenpy/solvers/BiCGSTAB.hpp
+    include/eigenpy/solvers/MINRES.hpp
     include/eigenpy/solvers/ConjugateGradient.hpp
     include/eigenpy/solvers/SparseSolverBase.hpp
     include/eigenpy/solvers/BasicPreconditioners.hpp
-    include/eigenpy/solvers/BFGSPreconditioners.hpp)
+    include/eigenpy/solvers/BFGSPreconditioners.hpp
+    include/eigenpy/solvers/IncompleteCholesky.hpp
+    include/eigenpy/solvers/IncompleteLUT.hpp)
 
 set(${PROJECT_NAME}_EIGEN_HEADERS include/eigenpy/eigen/EigenBase.hpp)
 
@@ -205,13 +209,18 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_CHOLMOD_HEADERS
     include/eigenpy/decompositions/sparse/cholmod/CholmodSupernodalLLT.hpp)
 
 set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_ACCELERATE_HEADERS
+    include/eigenpy/decompositions/sparse/accelerate/Accelerate.hpp
     include/eigenpy/decompositions/sparse/accelerate/accelerate.hpp)
 
 set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
-    include/eigenpy/decompositions/sparse/LLT.hpp
-    include/eigenpy/decompositions/sparse/LDLT.hpp
+    include/eigenpy/decompositions/sparse/SimplicialLLT.hpp
+    include/eigenpy/decompositions/sparse/SimplicialLDLT.hpp
+    include/eigenpy/decompositions/sparse/SparseLU.hpp
+    include/eigenpy/decompositions/sparse/SparseQR.hpp
     include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
-    include/eigenpy/decompositions/sparse/SparseSolverBase.hpp)
+    include/eigenpy/decompositions/sparse/SparseSolverBase.hpp
+    include/eigenpy/decompositions/sparse/LDLT.hpp
+    include/eigenpy/decompositions/sparse/LLT.hpp)
 
 if(BUILD_WITH_CHOLMOD_SUPPORT)
   list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
@@ -227,6 +236,16 @@ 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/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
+    include/eigenpy/decompositions/FullPivLU.hpp
+    include/eigenpy/decompositions/PartialPivLU.hpp
     include/eigenpy/decompositions/PermutationMatrix.hpp
     include/eigenpy/decompositions/LDLT.hpp
     include/eigenpy/decompositions/LLT.hpp
@@ -236,6 +255,9 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
     include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp
     include/eigenpy/decompositions/FullPivHouseholderQR.hpp
     include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
+    include/eigenpy/decompositions/SVDBase.hpp
+    include/eigenpy/decompositions/BDCSVD.hpp
+    include/eigenpy/decompositions/JacobiSVD.hpp
     include/eigenpy/decompositions/minres.hpp)
 
 set(${PROJECT_NAME}_HEADERS
@@ -300,17 +322,36 @@ list(
 # ----------------------------------------------------
 # --- TARGETS ----------------------------------------
 # ----------------------------------------------------
-set(${PROJECT_NAME}_SOLVERS_SOURCES src/solvers/preconditioners.cpp
-                                    src/solvers/solvers.cpp)
+set(${PROJECT_NAME}_SOLVERS_SOURCES
+    src/solvers/preconditioners.cpp
+    src/solvers/solvers.cpp
+    src/solvers/minres.cpp
+    src/solvers/bicgstab.cpp
+    src/solvers/conjugate-gradient.cpp
+    src/solvers/least-squares-conjugate-gradient.cpp
+    src/solvers/incomplete-cholesky.cpp
+    src/solvers/incomplete-lut.cpp)
 
 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
     src/decompositions/ldlt-solver.cpp
-    src/decompositions/minres-solver.cpp
+    src/decompositions/bdcsvd-solver.cpp
+    src/decompositions/jacobisvd-solver.cpp
+    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/sparse-lu-solver.cpp
+    src/decompositions/sparse-qr-solver.cpp
     src/decompositions/qr-solvers.cpp
-    src/decompositions/eigen-solver.cpp
     src/decompositions/self-adjoint-eigen-solver.cpp
     src/decompositions/permutation-matrix.cpp
     src/decompositions/simplicial-llt-solver.cpp

README.md

@@ -19,7 +19,7 @@ EigenPy — Versatile and efficient Python bindings between Numpy and Eigen
 - full support Eigen::Ref avoiding memory allocation
 - full support of the Eigen::Tensor module
 - exposition of the Geometry module of Eigen for easy code prototyping
-- standard matrix decomposion routines of Eigen such as the Cholesky decomposition (SVD and QR decompositions [can be added](#contributing))
+- standard matrix decomposion routines of Eigen such as the Cholesky, SVD and QR decompositions
 - full support of SWIG objects
 - full support of runtime declaration of Numpy scalar types
 - extended API to expose several STL types and some of their Boost equivalents: `optional` types, `std::pair`, maps, variants...
@@ -100,6 +100,7 @@ The following people have been involved in the development of **EigenPy**:
 - [Loïc Estève](https://github.com/lesteve) (Inria): Conda integration
 - [Wilson Jallet](https://manifoldfr.github.io/) (Inria): core developer
 - [Joris Vaillant](https://github.com/jorisv) (Inria): core developer and manager of the project
+- [Lucas Haubert](https://www.linkedin.com/in/lucas-haubert-b668a421a/) (Inria): core developer
 
 If you have taken part in the development of **EigenPy**, feel free to add your name and contribution here.
 

flake.nix

@@ -20,20 +20,29 @@
           devShells.default = pkgs.mkShell { inputsFrom = [ self'.packages.default ]; };
           packages = {
             default = self'.packages.eigenpy;
-            eigenpy = pkgs.python3Packages.eigenpy.overrideAttrs (_: {
-              src = pkgs.lib.fileset.toSource {
-                root = ./.;
-                fileset = pkgs.lib.fileset.unions [
-                  ./CMakeLists.txt
-                  ./doc
-                  ./include
-                  ./package.xml
-                  ./python
-                  ./src
-                  ./unittest
-                ];
-              };
-            });
+            eigen = pkgs.eigen.overrideAttrs {
+              # Apply https://gitlab.com/libeigen/eigen/-/merge_requests/977
+              postPatch = ''
+                substituteInPlace Eigen/src/SVD/BDCSVD.h \
+                  --replace-fail "if (l == 0) {" "if (i >= k && l == 0) {"
+              '';
+            };
+            eigenpy =
+              (pkgs.python3Packages.eigenpy.override { inherit (self'.packages) eigen; }).overrideAttrs
+                (_: {
+                  src = pkgs.lib.fileset.toSource {
+                    root = ./.;
+                    fileset = pkgs.lib.fileset.unions [
+                      ./CMakeLists.txt
+                      ./doc
+                      ./include
+                      ./package.xml
+                      ./python
+                      ./src
+                      ./unittest
+                    ];
+                  };
+                });
           };
         };
     };

include/eigenpy/decompositions/BDCSVD.hpp

@@ -0,0 +1,99 @@
+/*
+ * Copyright 2025 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_bdcsvd_hpp__
+#define __eigenpy_decompositions_bdcsvd_hpp__
+
+#include <Eigen/SVD>
+#include <Eigen/Core>
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
+#include "eigenpy/decompositions/SVDBase.hpp"
+
+namespace eigenpy {
+
+template <typename _MatrixType>
+struct BDCSVDVisitor
+    : public boost::python::def_visitor<BDCSVDVisitor<_MatrixType>> {
+  typedef _MatrixType MatrixType;
+  typedef Eigen::BDCSVD<MatrixType> Solver;
+  typedef typename MatrixType::Scalar Scalar;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
+        .def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex,
+                      bp::optional<unsigned int>>(
+            bp::args("self", "rows", "cols", "computationOptions "),
+            "Default Constructor with memory preallocation. "))
+        .def(bp::init<MatrixType, bp::optional<unsigned int>>(
+            bp::args("self", "matrix", "computationOptions "),
+            "Constructor performing the decomposition of given matrix."))
+
+        .def("cols", &Solver::cols, bp::arg("self"),
+             "Returns the number of columns. ")
+        .def("compute",
+             (Solver & (Solver::*)(const MatrixType &matrix)) & Solver::compute,
+             bp::args("self", "matrix"),
+             "Method performing the decomposition of given matrix. Computes "
+             "Thin/Full "
+             "unitaries U/V if specified using the Options template parameter "
+             "or the class constructor. ",
+             bp::return_self<>())
+        .def("compute",
+             (Solver & (Solver::*)(const MatrixType &matrix,
+                                   unsigned int computationOptions)) &
+                 Solver::compute,
+             bp::args("self", "matrix", "computationOptions"),
+             "Method performing the decomposition of given matrix, as "
+             "specified by the computationOptions parameter.  ",
+             bp::return_self<>())
+        .def("rows", &Solver::rows, bp::arg("self"),
+             "Returns the number of rows. ")
+        .def("setSwitchSize", &Solver::setSwitchSize, bp::args("self", "s"))
+
+        .def(SVDBaseVisitor<Solver>());
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "BDCSVD_" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string &name) {
+    bp::class_<Solver, boost::noncopyable>(
+        name.c_str(),
+        "Class Bidiagonal Divide and Conquer SVD.\n\n"
+        "This class first reduces the input matrix to bi-diagonal form using "
+        "class "
+        "UpperBidiagonalization, and then performs a divide-and-conquer "
+        "diagonalization. "
+        "Small blocks are diagonalized using class JacobiSVD. You can control "
+        "the "
+        "switching size with the setSwitchSize() method, default is 16. For "
+        "small matrice "
+        "(<16), it is thus preferable to directly use JacobiSVD. For larger "
+        "ones, BDCSVD "
+        "is highly recommended and can several order of magnitude faster.\n\n"
+        "Warming: this algorithm is unlikely to provide accurate result when "
+        "compiled with "
+        "unsafe math optimizations. For instance, this concerns Intel's "
+        "compiler (ICC), which "
+        "performs such optimization by default unless you compile with the "
+        "-fp-model precise "
+        "option. Likewise, the -ffast-math option of GCC or clang will "
+        "significantly degrade the "
+        "accuracy.",
+        bp::no_init)
+        .def(BDCSVDVisitor())
+        .def(IdVisitor<Solver>());
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decompositions_bdcsvd_hpp__

include/eigenpy/decompositions/ComplexEigenSolver.hpp

@@ -0,0 +1,86 @@
+/*
+ * Copyright 2025 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_complex_eigen_solver_hpp__
+#define __eigenpy_decompositions_complex_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 ComplexEigenSolverVisitor : public boost::python::def_visitor<
+                                       ComplexEigenSolverVisitor<_MatrixType>> {
+  typedef _MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef Eigen::ComplexEigenSolver<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, bp::optional<bool>>(
+            bp::args("matrix", "compute_eigen_vectors"),
+            "Computes eigendecomposition of given matrix"))
+
+        .def("eigenvalues", &Solver::eigenvalues, bp::arg("self"),
+             "Returns the eigenvalues of given matrix.",
+             bp::return_internal_reference<>())
+        .def("eigenvectors", &Solver::eigenvectors, bp::arg("self"),
+             "Returns the eigenvectors of given matrix.",
+             bp::return_internal_reference<>())
+
+        .def("compute", &ComplexEigenSolverVisitor::compute_proxy<MatrixType>,
+             bp::args("self", "matrix"),
+             "Computes the eigendecomposition of given matrix.",
+             bp::return_self<>())
+        .def("compute",
+             (Solver &
+              (Solver::*)(const Eigen::EigenBase<MatrixType>& matrix, bool)) &
+                 Solver::compute,
+             bp::args("self", "matrix", "compute_eigen_vectors"),
+             "Computes the eigendecomposition of given matrix.",
+             bp::return_self<>())
+
+        .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 =
+        "ComplexEigenSolver" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string& name) {
+    bp::class_<Solver>(name.c_str(), bp::no_init)
+        .def(ComplexEigenSolverVisitor())
+        .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_eigen_solver_hpp__