Skip to content

Commit 84d9bb7

Browse files
Lucas-HaubertLucas-Haubert
committed
completed llt and ldlt
1 parent 1ddfa87 commit 84d9bb7

File tree

9 files changed

+415
-27
lines changed

9 files changed

+415
-27
lines changed

‎include/nanoeigenpy/decompositions.hpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
#pragma once
2+
3+
#include "nanoeigenpy/decompositions/ldlt.hpp"
4+
#include "nanoeigenpy/decompositions/llt.hpp"

‎include/nanoeigenpy/decompositions/ldlt.hpp‎

Lines changed: 104 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -6,31 +6,124 @@
66
namespace nanoeigenpy {
77
namespace nb = nanobind;
88

9+
template <typename MatrixType, typename MatrixOrVector>
10+
MatrixOrVector solve(const Eigen::LDLT<MatrixType> &c, const MatrixOrVector &vec) {
11+
return c.solve(vec);
12+
}
13+
914
template <typename MatrixType>
1015
void exposeLDLTSolver(nb::module_ m, const char *name) {
1116
using Solver = Eigen::LDLT<MatrixType>;
1217
using Scalar = typename MatrixType::Scalar;
1318
using VectorType = Eigen::Matrix<Scalar, -1, 1>;
14-
auto cl = nb::class_<Solver>(m, name)
15-
.def(nb::init<const MatrixType &>(), nb::arg("matrix"))
16-
.def(nb::init<Eigen::DenseIndex>(), nb::arg("size"))
19+
auto cl = nb::class_<Solver>(m, name,
20+
"Robust Cholesky decomposition of a matrix with pivoting.\n\n"
21+
"Perform a robust Cholesky decomposition of a positive semidefinite or "
22+
"negative semidefinite matrix $ A $ such that $ A = P^TLDL^*P $, where "
23+
"P is a permutation matrix, L is lower triangular with a unit diagonal "
24+
"and D is a diagonal matrix.\n\n"
25+
"The decomposition uses pivoting to ensure stability, so that L will "
26+
"have zeros in the bottom right rank(A) - n submatrix. Avoiding the "
27+
"square root on D also stabilizes the computation.")
28+
29+
.def(nb::init<>(),
30+
"Default constructor.")
31+
.def(nb::init<Eigen::DenseIndex>(), nb::arg("size"),
32+
"Default constructor with memory preallocation.")
33+
.def(nb::init<const MatrixType &>(), nb::arg("matrix"),
34+
"Constructs a LLT factorization from a given matrix.")
35+
1736
.def(EigenBaseVisitor())
18-
.def("isNegative", &Solver::isNegative)
19-
.def("isPositive", &Solver::isPositive)
37+
38+
.def("isNegative", &Solver::isNegative,
39+
"Returns true if the matrix is negative (semidefinite).")
40+
.def("isPositive", &Solver::isPositive,
41+
"Returns true if the matrix is positive (semidefinite).")
42+
2043
.def("matrixL",
21-
[](Solver const &c) -> MatrixType { return c.matrixL(); })
44+
[](Solver const &c) -> MatrixType { return c.matrixL(); },
45+
"Returns the lower triangular matrix L.")
2246
.def("matrixU",
23-
[](Solver const &c) -> MatrixType { return c.matrixU(); })
47+
[](Solver const &c) -> MatrixType { return c.matrixU(); },
48+
"Returns the upper triangular matrix U.")
2449
.def("vectorD",
25-
[](Solver const &c) -> VectorType { return c.vectorD(); })
50+
[](Solver const &c) -> VectorType { return c.vectorD(); },
51+
"Returns the coefficients of the diagonal matrix D.")
2652
.def("matrixLDLT", &Solver::matrixLDLT,
53+
"Returns the LDLT decomposition matrix.",
2754
nb::rv_policy::reference_internal)
28-
.def(
29-
"rankUpdate",
55+
56+
.def("transpositionsP",
57+
[](Solver const &c) -> MatrixType { return c.transpositionsP() *
58+
MatrixType::Identity(c.matrixL().rows(), c.matrixL().rows()); },
59+
"Returns the permutation matrix P.")
60+
61+
.def("rankUpdate",
3062
[](Solver &c, VectorType const &w, Scalar sigma) {
3163
return c.rankUpdate(w, sigma);
3264
},
33-
nb::arg("w"), nb::arg("sigma"));
65+
nb::arg("w"), nb::arg("sigma"))
66+
67+
#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
68+
.def("adjoint", &Solver::adjoint,
69+
"Returns the adjoint, that is, a reference to the decomposition "
70+
"itself as if the underlying matrix is self-adjoint.",
71+
nb::rv_policy::reference)
72+
#endif
73+
74+
.def("compute",
75+
[](Solver &c, VectorType const &matrix) {
76+
return c.compute(matrix);
77+
},
78+
nb::arg("matrix"),
79+
"Computes the LDLT of given matrix.",
80+
nb::rv_policy::reference)
81+
.def("info", &Solver::info,
82+
"NumericalIssue if the input contains INF or NaN values or "
83+
"overflow occured. Returns Success otherwise.")
84+
85+
#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
86+
.def("rcond", &Solver::rcond,
87+
"Returns an estimate of the reciprocal condition number of the "
88+
"matrix.")
89+
#endif
90+
91+
.def("reconstructedMatrix", &Solver::reconstructedMatrix,
92+
"Returns the matrix represented by the decomposition, i.e., it "
93+
"returns the product: L L^*. This function is provided for debug "
94+
"purpose.")
95+
96+
.def("solve",
97+
[](Solver const &c, VectorType const &b) -> VectorType { return solve(c, b); },
98+
nb::arg("b"),
99+
"Returns the solution x of A x = b using the current "
100+
"decomposition of A.")
101+
.def("solve",
102+
[](Solver const &c, MatrixType const &B) -> MatrixType { return solve(c, B); },
103+
nb::arg("B"),
104+
"Returns the solution X of A X = B using the current "
105+
"decomposition of A where B is a right hand side matrix.")
106+
107+
.def("setZero", &Solver::setZero,
108+
"Clear any existing decomposition.")
109+
110+
.def("id",
111+
[](Solver const &c) -> int64_t { return reinterpret_cast<int64_t>(&c); },
112+
"Returns the unique identity of an object.\n"
113+
"For objects held in C++, it corresponds to its memory address.");
114+
34115
}
35116

36117
} // namespace nanoeigenpy
118+
119+
120+
// TODO
121+
122+
// Tests that were not done in eigenpy that we could add in nanoeigenpy: (+ those cited in llt.hpp)
123+
// setZero
124+
125+
// Expose supplementary content:
126+
// setZero in LLT decomp too ? (not done in eigenpy)
127+
128+
// Questions about eigenpy itself:
129+
// Relevant to have the reconstructedMatrix method ? (knowing that we are on LDLT, not LLT)

‎include/nanoeigenpy/decompositions/llt.hpp‎

Lines changed: 118 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -6,19 +6,127 @@
66
namespace nanoeigenpy {
77
namespace nb = nanobind;
88

9+
template <typename MatrixType, typename MatrixOrVector>
10+
MatrixOrVector solve(const Eigen::LLT<MatrixType> &c, const MatrixOrVector &vec) {
11+
return c.solve(vec);
12+
}
13+
914
template <typename MatrixType>
1015
void exposeLLTSolver(nb::module_ m, const char *name) {
1116
using Chol = Eigen::LLT<MatrixType>;
12-
auto cl = nb::class_<Chol>(m, name)
13-
.def(nb::init<const MatrixType &>(), nb::arg("matrix"))
14-
.def(nb::init<Eigen::DenseIndex>(), nb::arg("size"))
15-
.def(EigenBaseVisitor())
16-
.def("matrixL",
17-
[](Chol const &c) -> MatrixType { return c.matrixL(); })
18-
.def("matrixU",
19-
[](Chol const &c) -> MatrixType { return c.matrixU(); })
20-
.def("matrixLLT", &Chol::matrixLLT,
21-
nb::rv_policy::reference_internal);
17+
using Scalar = typename MatrixType::Scalar;
18+
using VectorType = Eigen::Matrix<Scalar, -1, 1>;
19+
auto cl = nb::class_<Chol>(m, name,
20+
"Standard Cholesky decomposition (LL^T) of a matrix and associated "
21+
"features.\n\n"
22+
"This class performs a LL^T Cholesky decomposition of a symmetric, "
23+
"positive definite matrix A such that A = LL^* = U^*U, where L is "
24+
"lower triangular.\n\n"
25+
"While the Cholesky decomposition is particularly useful to solve "
26+
"selfadjoint problems like D^*D x = b, for that purpose, we recommend "
27+
"the Cholesky decomposition without square root which is more stable "
28+
"and even faster. Nevertheless, this standard Cholesky decomposition "
29+
"remains useful in many other situations like generalised eigen "
30+
"problems with hermitian matrices.")
31+
32+
.def(nb::init<>(),
33+
"Default constructor.")
34+
.def(nb::init<Eigen::DenseIndex>(), nb::arg("size"),
35+
"Default constructor with memory preallocation.")
36+
.def(nb::init<const MatrixType &>(), nb::arg("matrix"),
37+
"Constructs a LLT factorization from a given matrix.")
38+
39+
.def(EigenBaseVisitor())
40+
41+
.def("matrixL",
42+
[](Chol const &c) -> MatrixType { return c.matrixL(); },
43+
"Returns the lower triangular matrix L.")
44+
.def("matrixU",
45+
[](Chol const &c) -> MatrixType { return c.matrixU(); },
46+
"Returns the upper triangular matrix U.")
47+
.def("matrixLLT", &Chol::matrixLLT,
48+
"Returns the LLT decomposition matrix.",
49+
nb::rv_policy::reference_internal)
50+
51+
#if EIGEN_VERSION_AT_LEAST(3, 3, 90)
52+
.def("rankUpdate",
53+
[](Chol &c, VectorType const &w, Scalar sigma) {
54+
return c.rankUpdate(w, sigma);
55+
},
56+
nb::arg("w"), nb::arg("sigma"),
57+
nb::rv_policy::reference)
58+
#else
59+
.def("rankUpdate",
60+
[](Chol &c, VectorType const &w, Scalar sigma) {
61+
return c.rankUpdate(w, sigma);
62+
},
63+
nb::arg("w"), nb::arg("sigma"))
64+
#endif
65+
66+
#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
67+
.def("adjoint", &Chol::adjoint,
68+
"Returns the adjoint, that is, a reference to the decomposition "
69+
"itself as if the underlying matrix is self-adjoint.",
70+
nb::rv_policy::reference)
71+
#endif
72+
73+
.def("compute",
74+
[](Chol &c, VectorType const &matrix) {
75+
return c.compute(matrix);
76+
},
77+
nb::arg("matrix"),
78+
"Computes the LDLT of given matrix.",
79+
nb::rv_policy::reference)
80+
.def("info", &Chol::info,
81+
"NumericalIssue if the input contains INF or NaN values or "
82+
"overflow occured. Returns Success otherwise.")
83+
84+
#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
85+
.def("rcond", &Chol::rcond,
86+
"Returns an estimate of the reciprocal condition number of the "
87+
"matrix.")
88+
#endif
89+
90+
.def("reconstructedMatrix", &Chol::reconstructedMatrix,
91+
"Returns the matrix represented by the decomposition, i.e., it "
92+
"returns the product: L L^*. This function is provided for debug "
93+
"purpose.")
94+
95+
.def("solve",
96+
[](Chol const &c, VectorType const &b) -> VectorType { return solve(c, b); },
97+
nb::arg("b"),
98+
"Returns the solution x of A x = b using the current "
99+
"decomposition of A.")
100+
.def("solve",
101+
[](Chol const &c, MatrixType const &B) -> MatrixType { return solve(c, B); },
102+
nb::arg("B"),
103+
"Returns the solution X of A X = B using the current "
104+
"decomposition of A where B is a right hand side matrix.")
105+
106+
.def("id",
107+
[](Chol const &c) -> int64_t { return reinterpret_cast<int64_t>(&c); },
108+
"Returns the unique identity of an object.\n"
109+
"For objects held in C++, it corresponds to its memory address.");
110+
22111
}
23112

24113
} // namespace nanoeigenpy
114+
115+
116+
// TODO
117+
118+
// Tests that were not done in eigenpy that we could add in nanoeigenpy:
119+
// Default constructor
120+
// Default constructor with memory preallocation
121+
// matrixLLT
122+
// rankUpdate
123+
// adjoint
124+
// compute
125+
// rcond
126+
// reconstructedMatrix
127+
128+
// Expose supplementary content:
129+
// Expose ComputationInfo to test info
130+
131+
// Assertions for the solve method on vectors x_est and b
132+
// Expose is_approx for vectors too (exposed and tested for ùatrices only in eigenpy)

‎include/nanoeigenpy/eigen-typedef.hpp‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
/// Copyright 2025 INRIA
2+
#pragma once
3+
4+
#include <Eigen/Core>
5+
#include <Eigen/SparseCore>
6+
7+
#define NANOEIGENPY_MAKE_TYPEDEFS(Type, Options, TypeSuffix, Size, SizeSuffix) \
8+
/** \ingroup matrixtypedefs */ \
9+
typedef Eigen::Matrix<Type, Size, Size, Options> \
10+
Matrix##SizeSuffix##TypeSuffix; \
11+
/** \ingroup matrixtypedefs */ \
12+
typedef Eigen::Matrix<Type, Size, 1> Vector##SizeSuffix##TypeSuffix; \
13+
/** \ingroup matrixtypedefs */ \
14+
typedef Eigen::Matrix<Type, 1, Size> RowVector##SizeSuffix##TypeSuffix;

‎include/nanoeigenpy/fwd.hpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
/// Copyright 2025 INRIA
2+
#pragma once
3+
4+
#define NANOEIGENPY_UNUSED_VARIABLE(var) (void)(var)
5+
#define NANOEIGENPY_UNUSED_TYPE(type) NANOEIGENPY_UNUSED_VARIABLE((type *)(NULL))

‎include/nanoeigenpy/utils/is-approx.hpp‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
/// Copyright 2025 INRIA
2+
#pragma once
3+
4+
#include <Eigen/Core>
5+
#include <Eigen/SparseCore>
6+
#include <nanobind/nanobind.h>
7+
8+
#include "nanoeigenpy/eigen-typedef.hpp"
9+
#include "nanoeigenpy/fwd.hpp"
10+
11+
namespace nanoeigenpy {
12+
13+
template <typename MatrixType1, typename MatrixType2>
14+
EIGEN_DONT_INLINE bool is_approx(const Eigen::MatrixBase<MatrixType1>& mat1,
15+
const Eigen::MatrixBase<MatrixType2>& mat2,
16+
const typename MatrixType1::RealScalar& prec) {
17+
return mat1.isApprox(mat2, prec);
18+
}
19+
20+
template <typename MatrixType1, typename MatrixType2>
21+
EIGEN_DONT_INLINE bool is_approx(const Eigen::MatrixBase<MatrixType1>& mat1,
22+
const Eigen::MatrixBase<MatrixType2>& mat2) {
23+
return is_approx(
24+
mat1, mat2,
25+
Eigen::NumTraits<typename MatrixType1::RealScalar>::dummy_precision());
26+
}
27+
28+
template <typename MatrixType1, typename MatrixType2>
29+
EIGEN_DONT_INLINE bool is_approx(
30+
const Eigen::SparseMatrixBase<MatrixType1>& mat1,
31+
const Eigen::SparseMatrixBase<MatrixType2>& mat2,
32+
const typename MatrixType1::RealScalar& prec) {
33+
return mat1.isApprox(mat2, prec);
34+
}
35+
36+
template <typename MatrixType1, typename MatrixType2>
37+
EIGEN_DONT_INLINE bool is_approx(
38+
const Eigen::SparseMatrixBase<MatrixType1>& mat1,
39+
const Eigen::SparseMatrixBase<MatrixType2>& mat2) {
40+
return is_approx(
41+
mat1, mat2,
42+
Eigen::NumTraits<typename MatrixType1::RealScalar>::dummy_precision());
43+
}
44+
45+
namespace nb = nanobind;
46+
47+
template <typename Scalar>
48+
void exposeIsApprox(nb::module_ m) {
49+
enum { Options = 0 };
50+
NANOEIGENPY_MAKE_TYPEDEFS(Scalar, Options, s, Eigen::Dynamic, X);
51+
NANOEIGENPY_UNUSED_TYPE(VectorXs);
52+
NANOEIGENPY_UNUSED_TYPE(RowVectorXs);
53+
typedef typename MatrixXs::RealScalar RealScalar;
54+
55+
using namespace Eigen;
56+
const RealScalar dummy_precision =
57+
Eigen::NumTraits<RealScalar>::dummy_precision();
58+
59+
// Exposition de la fonction is_approx pour matrices denses
60+
m.def("is_approx",
61+
[](const MatrixXs& mat1, const MatrixXs& mat2, RealScalar precision) {
62+
return is_approx(mat1, mat2, precision);
63+
},
64+
nb::arg("mat1"), nb::arg("mat2"),
65+
nb::arg("precision") = dummy_precision,
66+
"Check if two dense matrices are approximately equal.");
67+
}
68+
69+
} // namespace nanoeigenpy
70+

0 commit comments

Comments
 (0)