Skip to content

Commit 4553a01

Browse files
jcarpentjcarpent
committed
decompositions: expose CompleteOrthogonalDecomposition
1 parent 3c818e1 commit 4553a01

File tree

4 files changed

+208
-0
lines changed

4 files changed

+208
-0
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -213,6 +213,7 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS
213213
include/eigenpy/decompositions/QR.hpp
214214
include/eigenpy/decompositions/HouseholderQR.hpp
215215
include/eigenpy/decompositions/ColPivHouseholderQR.hpp
216+
include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp
216217
include/eigenpy/decompositions/FullPivHouseholderQR.hpp
217218
include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp
218219
include/eigenpy/decompositions/minres.hpp)

‎include/eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp‎

Lines changed: 204 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,204 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decompositions_complete_orthogonal_decomposition_hpp__
6+
#define __eigenpy_decompositions_complete_orthogonal_decomposition_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 CompleteOrthogonalDecompositionSolverVisitor
17+
: public boost::python::def_visitor<
18+
CompleteOrthogonalDecompositionSolverVisitor<_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::CompleteOrthogonalDecomposition<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+
"CompleteOrthogonalDecomposition.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>(bp::args("self", "matrix"),
43+
"Constructs a complete orthogonal "
44+
"factorization from a given matrix.\n"
45+
"This constructor computes the complete "
46+
"orthogonal factorization of the matrix "
47+
"matrix by calling the method compute()."))
48+
49+
.def("absDeterminant", &Self::absDeterminant, bp::arg("self"),
50+
"Returns the absolute value of the determinant of the matrix "
51+
"associated with the complete orthogonal decomposition.\n"
52+
"It has only linear complexity (that is, O(n) where n is the "
53+
"dimension of the square matrix) as the complete orthogonal "
54+
"decomposition has "
55+
"already been computed.\n"
56+
"Note: This is only for square matrices.")
57+
.def("logAbsDeterminant", &Self::logAbsDeterminant, bp::arg("self"),
58+
"Returns the natural log of the absolute value of the determinant "
59+
"of the matrix of which *this is the complete orthogonal "
60+
"decomposition.\n"
61+
"It has only linear complexity (that is, O(n) where n is the "
62+
"dimension of the square matrix) as the complete orthogonal "
63+
"decomposition has "
64+
"already been computed.\n"
65+
"Note: This is only for square matrices. This method is useful to "
66+
"work around the risk of overflow/underflow that's inherent to "
67+
"determinant computation.")
68+
.def("dimensionOfKernel", &Self::dimensionOfKernel, bp::arg("self"),
69+
"Returns the dimension of the kernel of the matrix of which *this "
70+
"is the complete orthogonal decomposition.")
71+
.def("info", &Self::info, bp::arg("self"),
72+
"Reports whether the complete orthogonal factorization was "
73+
"successful.\n"
74+
"Note: This function always returns Success. It is provided for "
75+
"compatibility with other factorization routines.")
76+
.def("isInjective", &Self::isInjective, bp::arg("self"),
77+
"Returns true if the matrix associated with this complete "
78+
"orthogonal decomposition "
79+
"represents an injective linear map, i.e. has trivial kernel; "
80+
"false otherwise.\n"
81+
"\n"
82+
"Note: This method has to determine which pivots should be "
83+
"considered nonzero. For that, it uses the threshold value that "
84+
"you can control by calling setThreshold(threshold).")
85+
.def("isInvertible", &Self::isInvertible, bp::arg("self"),
86+
"Returns true if the matrix associated with the complete "
87+
"orthogonal decomposition "
88+
"is invertible.\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("isSurjective", &Self::isSurjective, bp::arg("self"),
94+
"Returns true if the matrix associated with this complete "
95+
"orthogonal decomposition "
96+
"represents a surjective linear map; false otherwise.\n"
97+
"\n"
98+
"Note: This method has to determine which pivots should be "
99+
"considered nonzero. For that, it uses the threshold value that "
100+
"you can control by calling setThreshold(threshold).")
101+
.def("maxPivot", &Self::maxPivot, bp::arg("self"),
102+
"Returns the absolute value of the biggest pivot, i.e. the "
103+
"biggest diagonal coefficient of U.")
104+
.def("nonzeroPivots", &Self::nonzeroPivots, bp::arg("self"),
105+
"Returns the number of nonzero pivots in the complete orthogonal "
106+
"decomposition. "
107+
"Here nonzero is meant in the exact sense, not in a fuzzy sense. "
108+
"So that notion isn't really intrinsically interesting, but it is "
109+
"still useful when implementing algorithms.")
110+
.def("rank", &Self::rank, bp::arg("self"),
111+
"Returns the rank of the matrix associated with the complete "
112+
"orthogonal "
113+
"decomposition.\n"
114+
"\n"
115+
"Note: This method has to determine which pivots should be "
116+
"considered nonzero. For that, it uses the threshold value that "
117+
"you can control by calling setThreshold(threshold).")
118+
119+
.def("setThreshold",
120+
(Self & (Self::*)(const RealScalar &)) & Self::setThreshold,
121+
bp::args("self", "threshold"),
122+
"Allows to prescribe a threshold to be used by certain methods, "
123+
"such as rank(), who need to determine when pivots are to be "
124+
"considered nonzero. This is not used for the complete orthogonal "
125+
"decomposition "
126+
"itself.\n"
127+
"\n"
128+
"When it needs to get the threshold value, Eigen calls "
129+
"threshold(). By default, this uses a formula to automatically "
130+
"determine a reasonable threshold. Once you have called the "
131+
"present method setThreshold(const RealScalar&), your value is "
132+
"used instead.\n"
133+
"\n"
134+
"Note: A pivot will be considered nonzero if its absolute value "
135+
"is strictly greater than |pivot| ⩽ threshold×|maxpivot| where "
136+
"maxpivot is the biggest pivot.",
137+
bp::return_self<>())
138+
.def("threshold", &Self::threshold, bp::arg("self"),
139+
"Returns the threshold that will be used by certain methods such "
140+
"as rank().")
141+
142+
.def("matrixQTZ", &Self::matrixQTZ, bp::arg("self"),
143+
"Returns the matrix where the complete orthogonal decomposition "
144+
"is stored.",
145+
bp::return_value_policy<bp::copy_const_reference>())
146+
.def("matrixT", &Self::matrixT, bp::arg("self"),
147+
"Returns the matrix where the complete orthogonal decomposition "
148+
"is stored.",
149+
bp::return_value_policy<bp::copy_const_reference>())
150+
.def("matrixZ", &Self::matrixZ, bp::arg("self"),
151+
"Returns the matrix Z.")
152+
153+
.def(
154+
"compute",
155+
(Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
156+
Solver::compute,
157+
bp::args("self", "matrix"),
158+
"Computes the complete orthogonal factorization of given matrix.",
159+
bp::return_self<>())
160+
161+
.def("pseudoInverse", pseudoInverse, bp::arg("self"),
162+
"Returns the pseudo-inverse of the matrix associated with the "
163+
"complete orthogonal "
164+
"decomposition.")
165+
166+
.def("solve", &solve<MatrixXs>, bp::args("self", "B"),
167+
"Returns the solution X of A X = B using the current "
168+
"decomposition of A where B is a right hand side matrix.");
169+
}
170+
171+
static void expose() {
172+
static const std::string classname =
173+
"CompleteOrthogonalDecomposition" + scalar_name<Scalar>::shortname();
174+
expose(classname);
175+
}
176+
177+
static void expose(const std::string &name) {
178+
bp::class_<Solver>(
179+
name.c_str(),
180+
"This class performs a rank-revealing complete orthogonal "
181+
"decomposition of a matrix A into matrices P, Q, T, and Z such that:\n"
182+
"AP=Q[T000]Z"
183+
"by using Householder transformations. Here, P is a permutation "
184+
"matrix, Q and Z are unitary matrices and T an upper triangular matrix "
185+
"of size rank-by-rank. A may be rank deficient.",
186+
bp::no_init)
187+
.def(CompleteOrthogonalDecompositionSolverVisitor())
188+
.def(IdVisitor<Solver>());
189+
}
190+
191+
private:
192+
template <typename MatrixOrVector>
193+
static MatrixOrVector solve(const Solver &self, const MatrixOrVector &vec) {
194+
return self.solve(vec);
195+
}
196+
static MatrixXs pseudoInverse(const Self &self) {
197+
return self.pseudoInverse();
198+
}
199+
};
200+
201+
} // namespace eigenpy
202+
203+
#endif // ifndef
204+
// __eigenpy_decompositions_complete_orthogonal_decomposition_hpp__

‎include/eigenpy/decompositions/QR.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,5 +8,6 @@
88
#include "eigenpy/decompositions/HouseholderQR.hpp"
99
#include "eigenpy/decompositions/FullPivHouseholderQR.hpp"
1010
#include "eigenpy/decompositions/ColPivHouseholderQR.hpp"
11+
#include "eigenpy/decompositions/CompleteOrthogonalDecomposition.hpp"
1112

1213
#endif // ifndef __eigenpy_decompositions_qr_hpp__

‎src/decompositions/qr-solvers.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,5 +10,7 @@ void exposeQRSolvers() {
1010
HouseholderQRSolverVisitor<MatrixXd>::expose("HouseholderQR");
1111
FullPivHouseholderQRSolverVisitor<MatrixXd>::expose("FullPivHouseholderQR");
1212
ColPivHouseholderQRSolverVisitor<MatrixXd>::expose("ColPivHouseholderQR");
13+
CompleteOrthogonalDecompositionSolverVisitor<MatrixXd>::expose(
14+
"CompleteOrthogonalDecomposition");
1315
}
1416
} // namespace eigenpy

0 commit comments

Comments
 (0)