Skip to content

Commit c3cfbe6

Browse files
authored
Merge pull request #478 from jcarpent/topic/devel
Expose QR solvers from Eigen
2 parents f65f7ab + 499cd83 commit c3cfbe6

21 files changed

+843
-31
lines changed

CHANGELOG.md

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,11 +6,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
66

77
## [Unreleased]
88

9+
### Added
10+
- Added id() helper to retrieve unique object identifier in Python ([#477](https://github.com/stack-of-tasks/eigenpy/pull/477))
11+
- Expose QR solvers ([#478](https://github.com/stack-of-tasks/eigenpy/pull/478))
12+
913
## [3.6.0] - 2024-06-05
1014

1115
### Added
1216
- Added a deprecation call policy shortcut ([#466](https://github.com/stack-of-tasks/eigenpy/pull/466))
13-
- Added id() helper to retrieve unique object identifier in Python ([#477](https://github.com/stack-of-tasks/eigenpy/pull/477))
1417

1518
### Fixed
1619
- Fix register_symbolic_link_to_registered_type() for multiple successive registrations ([#471](https://github.com/stack-of-tasks/eigenpy/pull/471))

CMakeLists.txt

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -210,6 +210,11 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
210210
include/eigenpy/decompositions/PermutationMatrix.hpp
211211
include/eigenpy/decompositions/LDLT.hpp
212212
include/eigenpy/decompositions/LLT.hpp
213+
include/eigenpy/decompositions/QR.hpp
214+
include/eigenpy/decompositions/HouseholderQR.hpp
215+
include/eigenpy/decompositions/ColPivHouseholderQR.hpp
216+
include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp
217+
include/eigenpy/decompositions/FullPivHouseholderQR.hpp
213218
include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
214219
include/eigenpy/decompositions/minres.hpp)
215220

@@ -282,6 +287,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
282287
src/decompositions/llt-solver.cpp
283288
src/decompositions/ldlt-solver.cpp
284289
src/decompositions/minres-solver.cpp
290+
src/decompositions/qr-solvers.cpp
285291
src/decompositions/eigen-solver.cpp
286292
src/decompositions/self-adjoint-eigen-solver.cpp
287293
src/decompositions/permutation-matrix.cpp
Lines changed: 189 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,189 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decompositions_col_piv_houselholder_qr_hpp__
6+
#define __eigenpy_decompositions_col_piv_houselholder_qr_hpp__
7+
8+
#include "eigenpy/eigenpy.hpp"
9+
#include "eigenpy/utils/scalar-name.hpp"
10+
11+
#include <Eigen/QR>
12+
13+
namespace eigenpy {
14+
15+
template <typename _MatrixType>
16+
struct ColPivHouseholderQRSolverVisitor
17+
: public boost::python::def_visitor<
18+
ColPivHouseholderQRSolverVisitor<_MatrixType> > {
19+
typedef _MatrixType MatrixType;
20+
typedef typename MatrixType::Scalar Scalar;
21+
typedef typename MatrixType::RealScalar RealScalar;
22+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
23+
VectorXs;
24+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
25+
MatrixType::Options>
26+
MatrixXs;
27+
typedef Eigen::ColPivHouseholderQR<MatrixType> Solver;
28+
typedef Solver Self;
29+
30+
template <class PyClass>
31+
void visit(PyClass &cl) const {
32+
cl.def(bp::init<>(bp::arg("self"),
33+
"Default constructor.\n"
34+
"The default constructor is useful in cases in which the "
35+
"user intends to perform decompositions via "
36+
"HouseholderQR.compute(matrix)"))
37+
.def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex>(
38+
bp::args("self", "rows", "cols"),
39+
"Default constructor with memory preallocation.\n"
40+
"Like the default constructor but with preallocation of the "
41+
"internal data according to the specified problem size. "))
42+
.def(bp::init<MatrixType>(
43+
bp::args("self", "matrix"),
44+
"Constructs a QR factorization from a given matrix.\n"
45+
"This constructor computes the QR factorization of the matrix "
46+
"matrix by calling the method compute()."))
47+
48+
.def("absDeterminant", &Self::absDeterminant, bp::arg("self"),
49+
"Returns the absolute value of the determinant of the matrix of "
50+
"which *this is the QR decomposition.\n"
51+
"It has only linear complexity (that is, O(n) where n is the "
52+
"dimension of the square matrix) as the QR decomposition has "
53+
"already been computed.\n"
54+
"Note: This is only for square matrices.")
55+
.def("logAbsDeterminant", &Self::logAbsDeterminant, bp::arg("self"),
56+
"Returns the natural log of the absolute value of the determinant "
57+
"of the matrix of which *this is the QR decomposition.\n"
58+
"It has only linear complexity (that is, O(n) where n is the "
59+
"dimension of the square matrix) as the QR decomposition has "
60+
"already been computed.\n"
61+
"Note: This is only for square matrices. This method is useful to "
62+
"work around the risk of overflow/underflow that's inherent to "
63+
"determinant computation.")
64+
.def("dimensionOfKernel", &Self::dimensionOfKernel, bp::arg("self"),
65+
"Returns the dimension of the kernel of the matrix of which *this "
66+
"is the QR decomposition.")
67+
.def("info", &Self::info, bp::arg("self"),
68+
"Reports whether the QR factorization was successful.\n"
69+
"Note: This function always returns Success. It is provided for "
70+
"compatibility with other factorization routines.")
71+
.def("isInjective", &Self::isInjective, bp::arg("self"),
72+
"Returns true if the matrix associated with this QR decomposition "
73+
"represents an injective linear map, i.e. has trivial kernel; "
74+
"false otherwise.\n"
75+
"\n"
76+
"Note: This method has to determine which pivots should be "
77+
"considered nonzero. For that, it uses the threshold value that "
78+
"you can control by calling setThreshold(threshold).")
79+
.def("isInvertible", &Self::isInvertible, bp::arg("self"),
80+
"Returns true if the matrix associated with the QR decomposition "
81+
"is invertible.\n"
82+
"\n"
83+
"Note: This method has to determine which pivots should be "
84+
"considered nonzero. For that, it uses the threshold value that "
85+
"you can control by calling setThreshold(threshold).")
86+
.def("isSurjective", &Self::isSurjective, bp::arg("self"),
87+
"Returns true if the matrix associated with this QR decomposition "
88+
"represents a surjective linear map; false otherwise.\n"
89+
"\n"
90+
"Note: This method has to determine which pivots should be "
91+
"considered nonzero. For that, it uses the threshold value that "
92+
"you can control by calling setThreshold(threshold).")
93+
.def("maxPivot", &Self::maxPivot, bp::arg("self"),
94+
"Returns the absolute value of the biggest pivot, i.e. the "
95+
"biggest diagonal coefficient of U.")
96+
.def("nonzeroPivots", &Self::nonzeroPivots, bp::arg("self"),
97+
"Returns the number of nonzero pivots in the QR decomposition. "
98+
"Here nonzero is meant in the exact sense, not in a fuzzy sense. "
99+
"So that notion isn't really intrinsically interesting, but it is "
100+
"still useful when implementing algorithms.")
101+
.def("rank", &Self::rank, bp::arg("self"),
102+
"Returns the rank of the matrix associated with the QR "
103+
"decomposition.\n"
104+
"\n"
105+
"Note: This method has to determine which pivots should be "
106+
"considered nonzero. For that, it uses the threshold value that "
107+
"you can control by calling setThreshold(threshold).")
108+
109+
.def("setThreshold",
110+
(Self & (Self::*)(const RealScalar &)) & Self::setThreshold,
111+
bp::args("self", "threshold"),
112+
"Allows to prescribe a threshold to be used by certain methods, "
113+
"such as rank(), who need to determine when pivots are to be "
114+
"considered nonzero. This is not used for the QR decomposition "
115+
"itself.\n"
116+
"\n"
117+
"When it needs to get the threshold value, Eigen calls "
118+
"threshold(). By default, this uses a formula to automatically "
119+
"determine a reasonable threshold. Once you have called the "
120+
"present method setThreshold(const RealScalar&), your value is "
121+
"used instead.\n"
122+
"\n"
123+
"Note: A pivot will be considered nonzero if its absolute value "
124+
"is strictly greater than |pivot| ⩽ threshold×|maxpivot| where "
125+
"maxpivot is the biggest pivot.",
126+
bp::return_self<>())
127+
.def("threshold", &Self::threshold, bp::arg("self"),
128+
"Returns the threshold that will be used by certain methods such "
129+
"as rank().")
130+
131+
.def("matrixQR", &Self::matrixQR, bp::arg("self"),
132+
"Returns the matrix where the Householder QR decomposition is "
133+
"stored in a LAPACK-compatible way.",
134+
bp::return_value_policy<bp::copy_const_reference>())
135+
.def("matrixR", &Self::matrixR, bp::arg("self"),
136+
"Returns the matrix where the result Householder QR is stored.",
137+
bp::return_value_policy<bp::copy_const_reference>())
138+
139+
.def(
140+
"compute",
141+
(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
142+
Solver::compute,
143+
bp::args("self", "matrix"),
144+
"Computes the QR factorization of given matrix.",
145+
bp::return_self<>())
146+
147+
.def("inverse", inverse, bp::arg("self"),
148+
"Returns the inverse of the matrix associated with the QR "
149+
"decomposition..")
150+
151+
.def("solve", &solve<MatrixXs>, bp::args("self", "B"),
152+
"Returns the solution X of A X = B using the current "
153+
"decomposition of A where B is a right hand side matrix.");
154+
}
155+
156+
static void expose() {
157+
static const std::string classname =
158+
"ColPivHouseholderQR" + scalar_name<Scalar>::shortname();
159+
expose(classname);
160+
}
161+
162+
static void expose(const std::string &name) {
163+
bp::class_<Solver>(
164+
name.c_str(),
165+
"This class performs a rank-revealing QR decomposition of a matrix A "
166+
"into matrices P, Q and R such that:\n"
167+
"AP=QR\n"
168+
"by using Householder transformations. Here, P is a permutation "
169+
"matrix, Q a unitary matrix and R an upper triangular matrix.\n"
170+
"\n"
171+
"This decomposition performs column pivoting in order to be "
172+
"rank-revealing and improve numerical stability. It is slower than "
173+
"HouseholderQR, and faster than FullPivHouseholderQR.",
174+
bp::no_init)
175+
.def(ColPivHouseholderQRSolverVisitor())
176+
.def(IdVisitor<Solver>());
177+
}
178+
179+
private:
180+
template <typename MatrixOrVector>
181+
static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
182+
return self.solve(vec);
183+
}
184+
static MatrixXs inverse(const Self &self) { return self.inverse(); }
185+
};
186+
187+
} // namespace eigenpy
188+
189+
#endif // ifndef __eigenpy_decompositions_col_piv_houselholder_qr_hpp__

0 commit comments

Comments
 (0)