Expose QR solvers from Eigen

CHANGELOG.md

@@ -11,6 +11,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 ### Added
 - Added a deprecation call policy shortcut ([#466](https://github.com/stack-of-tasks/eigenpy/pull/466))
 - Added id() helper to retrieve unique object identifier in Python ([#477](https://github.com/stack-of-tasks/eigenpy/pull/477))
+- Expose QR solvers ([#478](https://github.com/stack-of-tasks/eigenpy/pull/478))
 
 ### Fixed
 - Fix register_symbolic_link_to_registered_type() for multiple successive registrations ([#471](https://github.com/stack-of-tasks/eigenpy/pull/471))

CMakeLists.txt

@@ -210,6 +210,8 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
     include/eigenpy/decompositions/PermutationMatrix.hpp
     include/eigenpy/decompositions/LDLT.hpp
     include/eigenpy/decompositions/LLT.hpp
+    include/eigenpy/decompositions/QR.hpp
+    include/eigenpy/decompositions/HouseholderQR.hpp
     include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
     include/eigenpy/decompositions/minres.hpp)
 
@@ -282,6 +284,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
     src/decompositions/llt-solver.cpp
     src/decompositions/ldlt-solver.cpp
     src/decompositions/minres-solver.cpp
+    src/decompositions/qr-solvers.cpp
     src/decompositions/eigen-solver.cpp
     src/decompositions/self-adjoint-eigen-solver.cpp
     src/decompositions/permutation-matrix.cpp

include/eigenpy/decompositions/EigenSolver.hpp

@@ -2,8 +2,8 @@
  * Copyright 2020 INRIA
  */
 
-#ifndef __eigenpy_decomposition_eigen_solver_hpp__
-#define __eigenpy_decomposition_eigen_solver_hpp__
+#ifndef __eigenpy_decompositions_eigen_solver_hpp__
+#define __eigenpy_decompositions_eigen_solver_hpp__
 
 #include <Eigen/Core>
 #include <Eigen/Eigenvalues>
@@ -90,4 +90,4 @@ struct EigenSolverVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_eigen_solver_hpp__
+#endif  // ifndef __eigenpy_decompositions_eigen_solver_hpp__

include/eigenpy/decompositions/HouseholderQR.hpp

@@ -0,0 +1,119 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_houselholder_qr_hpp__
+#define __eigenpy_decompositions_houselholder_qr_hpp__
+
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/utils/scalar-name.hpp"
+#include "eigenpy/eigen/EigenBase.hpp"
+
+#include <Eigen/QR>
+
+namespace eigenpy {
+
+template <typename _MatrixType>
+struct HouseholderQRSolverVisitor
+    : public boost::python::def_visitor<
+          HouseholderQRSolverVisitor<_MatrixType> > {
+  typedef _MatrixType MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
+      VectorXs;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
+                        MatrixType::Options>
+      MatrixXs;
+  typedef Eigen::HouseholderQR<MatrixType> Solver;
+  typedef Solver Self;
+
+  template <class PyClass>
+  void visit(PyClass &cl) const {
+    cl.def(bp::init<>(bp::arg("self"),
+                      "Default constructor.\n"
+                      "The default constructor is useful in cases in which the "
+                      "user intends to perform decompositions via "
+                      "HouseholderQR.compute(matrix)"))
+        .def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex>(
+            bp::args("self", "rows", "cols"),
+            "Default constructor with memory preallocation.\n"
+            "Like the default constructor but with preallocation of the "
+            "internal data according to the specified problem size. "))
+        .def(bp::init<MatrixType>(
+            bp::args("self", "matrix"),
+            "Constructs a QR factorization from a given matrix.\n"
+            "This constructor computes the QR factorization of the matrix "
+            "matrix by calling the method compute()."))
+
+        .def("absDeterminant", &Self::absDeterminant, bp::arg("self"),
+             "Returns the absolute value of the determinant of the matrix of "
+             "which *this is the QR decomposition.\n"
+             "It has only linear complexity (that is, O(n) where n is the "
+             "dimension of the square matrix) as the QR decomposition has "
+             "already been computed.\n"
+             "Note: This is only for square matrices.")
+        .def("logAbsDeterminant", &Self::logAbsDeterminant, bp::arg("self"),
+             "Returns the natural log of the absolute value of the determinant "
+             "of the matrix of which *this is the QR decomposition.\n"
+             "It has only linear complexity (that is, O(n) where n is the "
+             "dimension of the square matrix) as the QR decomposition has "
+             "already been computed.\n"
+             "Note: This is only for square matrices. This method is useful to "
+             "work around the risk of overflow/underflow that's inherent to "
+             "determinant computation.")
+
+        .def("matrixQR", &Self::matrixQR, bp::arg("self"),
+             "Returns the matrix where the Householder QR decomposition is "
+             "stored in a LAPACK-compatible way.",
+             bp::return_value_policy<bp::copy_const_reference>())
+
+        .def(
+            "compute",
+            (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
+                Solver::compute,
+            bp::args("self", "matrix"),
+            "Computes the QR factorization of given matrix.",
+            bp::return_self<>())
+
+        .def("solve", &solve<MatrixXs>, bp::args("self", "B"),
+             "Returns the solution X of A X = B using the current "
+             "decomposition of A where B is a right hand side matrix.");
+  }
+
+  static void expose() {
+    static const std::string classname =
+        "LDLT" + scalar_name<Scalar>::shortname();
+    expose(classname);
+  }
+
+  static void expose(const std::string &name) {
+    bp::class_<Solver>(
+        name.c_str(),
+        "This class performs a QR decomposition of a matrix A into matrices Q "
+        "and R such that A=QR by using Householder transformations.\n"
+        "Here, Q a unitary matrix and R an upper triangular matrix. The result "
+        "is stored in a compact way compatible with LAPACK.\n"
+        "\n"
+        "Note that no pivoting is performed. This is not a rank-revealing "
+        "decomposition. If you want that feature, use FullPivHouseholderQR or "
+        "ColPivHouseholderQR instead.\n"
+        "\n"
+        "This Householder QR decomposition is faster, but less numerically "
+        "stable and less feature-full than FullPivHouseholderQR or "
+        "ColPivHouseholderQR.",
+        bp::no_init)
+        .def(HouseholderQRSolverVisitor())
+        .def(IdVisitor<Solver>());
+  }
+
+ private:
+  template <typename MatrixOrVector>
+  static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
+    return self.solve(vec);
+  }
+};
+
+}  // namespace eigenpy
+
+#endif  // ifndef __eigenpy_decompositions_houselholder_qr_hpp__

include/eigenpy/decompositions/LDLT.hpp

@@ -2,8 +2,8 @@
  * Copyright 2020-2024 INRIA
  */
 
-#ifndef __eigenpy_decomposition_ldlt_hpp__
-#define __eigenpy_decomposition_ldlt_hpp__
+#ifndef __eigenpy_decompositions_ldlt_hpp__
+#define __eigenpy_decompositions_ldlt_hpp__
 
 #include <Eigen/Cholesky>
 #include <Eigen/Core>
@@ -141,4 +141,4 @@ struct LDLTSolverVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_ldlt_hpp__
+#endif  // ifndef __eigenpy_decompositions_ldlt_hpp__

include/eigenpy/decompositions/LLT.hpp

@@ -2,8 +2,8 @@
  * Copyright 2020-2024 INRIA
  */
 
-#ifndef __eigenpy_decomposition_llt_hpp__
-#define __eigenpy_decomposition_llt_hpp__
+#ifndef __eigenpy_decompositions_llt_hpp__
+#define __eigenpy_decompositions_llt_hpp__
 
 #include <Eigen/Cholesky>
 #include <Eigen/Core>
@@ -131,4 +131,4 @@ struct LLTSolverVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_llt_hpp__
+#endif  // ifndef __eigenpy_decompositions_llt_hpp__

include/eigenpy/decompositions/PermutationMatrix.hpp

@@ -2,8 +2,8 @@
  * Copyright 2024 INRIA
  */
 
-#ifndef __eigenpy_decomposition_permutation_matrix_hpp__
-#define __eigenpy_decomposition_permutation_matrix_hpp__
+#ifndef __eigenpy_decompositions_permutation_matrix_hpp__
+#define __eigenpy_decompositions_permutation_matrix_hpp__
 
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/eigen/EigenBase.hpp"
@@ -104,4 +104,4 @@ struct PermutationMatrixVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_permutation_matrix_hpp__
+#endif  // ifndef __eigenpy_decompositions_permutation_matrix_hpp__

include/eigenpy/decompositions/QR.hpp

@@ -0,0 +1,10 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#ifndef __eigenpy_decompositions_qr_hpp__
+#define __eigenpy_decompositions_qr_hpp__
+
+#include "eigenpy/decompositions/HouseholderQR.hpp"
+
+#endif  // ifndef __eigenpy_decompositions_qr_hpp__

include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp

@@ -2,8 +2,8 @@
  * Copyright 2020-2024 INRIA
  */
 
-#ifndef __eigenpy_decomposition_self_adjoint_eigen_solver_hpp__
-#define __eigenpy_decomposition_self_adjoint_eigen_solver_hpp__
+#ifndef __eigenpy_decompositions_self_adjoint_eigen_solver_hpp__
+#define __eigenpy_decompositions_self_adjoint_eigen_solver_hpp__
 
 #include <Eigen/Core>
 #include <Eigen/Eigenvalues>
@@ -102,4 +102,4 @@ struct SelfAdjointEigenSolverVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_self_adjoint_eigen_solver_hpp__
+#endif  // ifndef __eigenpy_decompositions_self_adjoint_eigen_solver_hpp__

include/eigenpy/decompositions/minres.hpp

@@ -2,8 +2,8 @@
  * Copyright 2021 INRIA
  */
 
-#ifndef __eigenpy_decomposition_minres_hpp__
-#define __eigenpy_decomposition_minres_hpp__
+#ifndef __eigenpy_decompositions_minres_hpp__
+#define __eigenpy_decompositions_minres_hpp__
 
 #include <Eigen/Core>
 #include <iostream>
@@ -178,4 +178,4 @@ struct MINRESSolverVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_minres_hpp__
+#endif  // ifndef __eigenpy_decompositions_minres_hpp__

include/eigenpy/decompositions/sparse/LDLT.hpp

@@ -2,8 +2,8 @@
  * Copyright 2024 INRIA
  */
 
-#ifndef __eigenpy_decomposition_sparse_ldlt_hpp__
-#define __eigenpy_decomposition_sparse_ldlt_hpp__
+#ifndef __eigenpy_decompositions_sparse_ldlt_hpp__
+#define __eigenpy_decompositions_sparse_ldlt_hpp__
 
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp"
@@ -69,4 +69,4 @@ struct SimplicialLDLTVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_sparse_ldlt_hpp__
+#endif  // ifndef __eigenpy_decompositions_sparse_ldlt_hpp__

include/eigenpy/decompositions/sparse/LLT.hpp

@@ -2,8 +2,8 @@
  * Copyright 2024 INRIA
  */
 
-#ifndef __eigenpy_decomposition_sparse_llt_hpp__
-#define __eigenpy_decomposition_sparse_llt_hpp__
+#ifndef __eigenpy_decompositions_sparse_llt_hpp__
+#define __eigenpy_decompositions_sparse_llt_hpp__
 
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/decompositions/sparse/SimplicialCholesky.hpp"
@@ -64,4 +64,4 @@ struct SimplicialLLTVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_sparse_llt_hpp__
+#endif  // ifndef __eigenpy_decompositions_sparse_llt_hpp__

include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp

@@ -2,8 +2,8 @@
  * Copyright 2024 INRIA
  */
 
-#ifndef __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
-#define __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
+#ifndef __eigenpy_decompositions_sparse_simplicial_cholesky_hpp__
+#define __eigenpy_decompositions_sparse_simplicial_cholesky_hpp__
 
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/eigen/EigenBase.hpp"
@@ -91,4 +91,4 @@ struct SimplicialCholeskyVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_sparse_simplicial_cholesky_hpp__
+#endif  // ifndef __eigenpy_decompositions_sparse_simplicial_cholesky_hpp__

include/eigenpy/decompositions/sparse/SparseSolverBase.hpp

@@ -2,8 +2,8 @@
  * Copyright 2024 INRIA
  */
 
-#ifndef __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
-#define __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
+#ifndef __eigenpy_decompositions_sparse_sparse_solver_base_hpp__
+#define __eigenpy_decompositions_sparse_sparse_solver_base_hpp__
 
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/eigen/EigenBase.hpp"
@@ -51,4 +51,4 @@ struct SparseSolverBaseVisitor
 
 }  // namespace eigenpy
 
-#endif  // ifndef __eigenpy_decomposition_sparse_sparse_solver_base_hpp__
+#endif  // ifndef __eigenpy_decompositions_sparse_sparse_solver_base_hpp__

src/decompositions/decompositions.cpp

@@ -12,6 +12,7 @@ void exposeEigenSolver();
 void exposeSelfAdjointEigenSolver();
 void exposeLLTSolver();
 void exposeLDLTSolver();
+void exposeQRSolvers();
 void exposeMINRESSolver();
 void exposeSimplicialLLTSolver();
 void exposeSimplicialLDLTSolver();
@@ -24,6 +25,7 @@ void exposeDecompositions() {
   exposeSelfAdjointEigenSolver();
   exposeLLTSolver();
   exposeLDLTSolver();
+  exposeQRSolvers();
   exposeMINRESSolver();
 
   {

src/decompositions/qr-solvers.cpp

@@ -0,0 +1,12 @@
+/*
+ * Copyright 2024 INRIA
+ */
+
+#include "eigenpy/decompositions/QR.hpp"
+
+namespace eigenpy {
+void exposeQRSolvers() {
+  using namespace Eigen;
+  HouseholderQRSolverVisitor<MatrixXd>::expose("HouseholderQR");
+}
+}  // namespace eigenpy

unittest/CMakeLists.txt

@@ -136,6 +136,8 @@ add_python_eigenpy_lib_unit_test("py-LDLT" "unittest/python/test_LDLT.py")
 
 add_python_eigenpy_lib_unit_test("py-id" "unittest/python/test_id.py")
 
+add_python_eigenpy_lib_unit_test("py-QR" "unittest/python/test_QR.py")
+
 if(NOT WIN32)
   add_python_eigenpy_lib_unit_test("py-MINRES" "unittest/python/test_MINRES.py")
 endif(NOT WIN32)

unittest/python/test_QR.py

@@ -0,0 +1,22 @@
+import numpy as np
+
+import eigenpy
+
+rows = 20
+cols = 100
+rng = np.random.default_rng()
+
+A = rng.random((rows, cols))
+
+# Test HouseholderQR decomposition
+householder_qr = eigenpy.HouseholderQR()
+householder_qr = eigenpy.HouseholderQR(rows, cols)
+householder_qr = eigenpy.HouseholderQR(A)
+
+householder_qr_eye = eigenpy.HouseholderQR(np.eye(rows, rows))
+X = rng.random((rows, 20))
+assert householder_qr_eye.absDeterminant() == 1.0
+assert householder_qr_eye.logAbsDeterminant() == 0.0
+
+Y = householder_qr_eye.solve(X)
+assert (X == Y).all()