Skip to content

Commit 6c85efc

Browse files
jcarpentjcarpent
committed
decompositions: expose FullPivHouseholderQR
1 parent d6fb4e9 commit 6c85efc

File tree

4 files changed

+186
-0
lines changed

4 files changed

+186
-0
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -212,6 +212,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
212212
include/eigenpy/decompositions/LLT.hpp
213213
include/eigenpy/decompositions/QR.hpp
214214
include/eigenpy/decompositions/HouseholderQR.hpp
215+
include/eigenpy/decompositions/FullPivHouseholderQR.hpp
215216
include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
216217
include/eigenpy/decompositions/minres.hpp)
217218

‎include/eigenpy/decompositions/FullPivHouseholderQR.hpp‎

Lines changed: 183 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,183 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decompositions_full_piv_houselholder_qr_hpp__
6+
#define __eigenpy_decompositions_full_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 FullPivHouseholderQRSolverVisitor
17+
: public boost::python::def_visitor<
18+
FullPivHouseholderQRSolverVisitor<_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::FullPivHouseholderQR<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("isInjective", &Self::isInjective, bp::arg("self"),
68+
"Returns true if the matrix associated with this QR decomposition "
69+
"represents an injective linear map, i.e. has trivial kernel; "
70+
"false otherwise.\n"
71+
"\n"
72+
"Note: This method has to determine which pivots should be "
73+
"considered nonzero. For that, it uses the threshold value that "
74+
"you can control by calling setThreshold(threshold).")
75+
.def("isInvertible", &Self::isInvertible, bp::arg("self"),
76+
"Returns true if the matrix associated with the QR decomposition "
77+
"is invertible.\n"
78+
"\n"
79+
"Note: This method has to determine which pivots should be "
80+
"considered nonzero. For that, it uses the threshold value that "
81+
"you can control by calling setThreshold(threshold).")
82+
.def("isSurjective", &Self::isSurjective, bp::arg("self"),
83+
"Returns true if the matrix associated with this QR decomposition "
84+
"represents a surjective linear map; false otherwise.\n"
85+
"\n"
86+
"Note: This method has to determine which pivots should be "
87+
"considered nonzero. For that, it uses the threshold value that "
88+
"you can control by calling setThreshold(threshold).")
89+
.def("maxPivot", &Self::maxPivot, bp::arg("self"),
90+
"Returns the absolute value of the biggest pivot, i.e. the "
91+
"biggest diagonal coefficient of U.")
92+
.def("nonzeroPivots", &Self::nonzeroPivots, bp::arg("self"),
93+
"Returns the number of nonzero pivots in the QR decomposition. "
94+
"Here nonzero is meant in the exact sense, not in a fuzzy sense. "
95+
"So that notion isn't really intrinsically interesting, but it is "
96+
"still useful when implementing algorithms.")
97+
.def("rank", &Self::rank, bp::arg("self"),
98+
"Returns the rank of the matrix associated with the QR "
99+
"decomposition.\n"
100+
"\n"
101+
"Note: This method has to determine which pivots should be "
102+
"considered nonzero. For that, it uses the threshold value that "
103+
"you can control by calling setThreshold(threshold).")
104+
105+
.def("setThreshold",
106+
(Self & (Self::*)(const RealScalar &)) & Self::setThreshold,
107+
bp::args("self", "threshold"),
108+
"Allows to prescribe a threshold to be used by certain methods, "
109+
"such as rank(), who need to determine when pivots are to be "
110+
"considered nonzero. This is not used for the QR decomposition "
111+
"itself.\n"
112+
"\n"
113+
"When it needs to get the threshold value, Eigen calls "
114+
"threshold(). By default, this uses a formula to automatically "
115+
"determine a reasonable threshold. Once you have called the "
116+
"present method setThreshold(const RealScalar&), your value is "
117+
"used instead.\n"
118+
"\n"
119+
"Note: A pivot will be considered nonzero if its absolute value "
120+
"is strictly greater than |pivot| ⩽ threshold×|maxpivot| where "
121+
"maxpivot is the biggest pivot.",
122+
bp::return_self<>())
123+
.def("threshold", &Self::threshold, bp::arg("self"),
124+
"Returns the threshold that will be used by certain methods such "
125+
"as rank().")
126+
127+
.def("matrixQR", &Self::matrixQR, bp::arg("self"),
128+
"Returns the matrix where the Householder QR decomposition is "
129+
"stored in a LAPACK-compatible way.",
130+
bp::return_value_policy<bp::copy_const_reference>())
131+
132+
.def(
133+
"compute",
134+
(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
135+
Solver::compute,
136+
bp::args("self", "matrix"),
137+
"Computes the QR factorization of given matrix.",
138+
bp::return_self<>())
139+
140+
.def("inverse", inverse, bp::arg("self"),
141+
"Returns the inverse of the matrix associated with the QR "
142+
"decomposition..")
143+
144+
.def("solve", &solve<MatrixXs>, bp::args("self", "B"),
145+
"Returns the solution X of A X = B using the current "
146+
"decomposition of A where B is a right hand side matrix.");
147+
}
148+
149+
static void expose() {
150+
static const std::string classname =
151+
"FullPivHouseholderQR" + scalar_name<Scalar>::shortname();
152+
expose(classname);
153+
}
154+
155+
static void expose(const std::string &name) {
156+
bp::class_<Solver>(
157+
name.c_str(),
158+
"This class performs a rank-revealing QR decomposition of a matrix A "
159+
"into matrices P, P', Q and R such that:\n"
160+
"PAP′=QR\n"
161+
"by using Householder transformations. Here, P and P' are permutation "
162+
"matrices, Q a unitary matrix and R an upper triangular matrix.\n"
163+
"\n"
164+
"This decomposition performs a very prudent full pivoting in order to "
165+
"be rank-revealing and achieve optimal numerical stability. The "
166+
"trade-off is that it is slower than HouseholderQR and "
167+
"ColPivHouseholderQR.",
168+
bp::no_init)
169+
.def(FullPivHouseholderQRSolverVisitor())
170+
.def(IdVisitor<Solver>());
171+
}
172+
173+
private:
174+
template <typename MatrixOrVector>
175+
static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
176+
return self.solve(vec);
177+
}
178+
static MatrixXs inverse(const Self &self) { return self.inverse(); }
179+
};
180+
181+
} // namespace eigenpy
182+
183+
#endif // ifndef __eigenpy_decompositions_full_piv_houselholder_qr_hpp__

‎include/eigenpy/decompositions/QR.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,5 +6,6 @@
66
#define __eigenpy_decompositions_qr_hpp__
77

88
#include "eigenpy/decompositions/HouseholderQR.hpp"
9+
#include "eigenpy/decompositions/FullPivHouseholderQR.hpp"
910

1011
#endif // ifndef __eigenpy_decompositions_qr_hpp__

‎src/decompositions/qr-solvers.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,5 +8,6 @@ namespace eigenpy {
88
void exposeQRSolvers() {
99
using namespace Eigen;
1010
HouseholderQRSolverVisitor<MatrixXd>::expose("HouseholderQR");
11+
FullPivHouseholderQRSolverVisitor<MatrixXd>::expose("FullPivHouseholderQR");
1112
}
1213
} // namespace eigenpy

0 commit comments

Comments
 (0)