Skip to content

Commit f105185

Browse files
Lucas-HaubertLucas-Haubert
committed
Add IncompleteCholesky and IncompleteLUT
1 parent e1cb1e9 commit f105185

File tree

8 files changed

+388
-2
lines changed

8 files changed

+388
-2
lines changed

‎CMakeLists.txt‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -213,7 +213,9 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
213213
include/eigenpy/decompositions/sparse/LDLT.hpp
214214
include/eigenpy/decompositions/sparse/LU.hpp
215215
include/eigenpy/decompositions/sparse/SimplicialCholesky.hpp
216-
include/eigenpy/decompositions/sparse/SparseSolverBase.hpp)
216+
include/eigenpy/decompositions/sparse/SparseSolverBase.hpp
217+
include/eigenpy/decompositions/sparse/IncompleteCholesky.hpp
218+
include/eigenpy/decompositions/sparse/IncompleteLUT.hpp)
217219

218220
if(BUILD_WITH_CHOLMOD_SUPPORT)
219221
list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SPARSE_HEADERS
@@ -340,7 +342,9 @@ set(${PROJECT_NAME}_DECOMPOSITIONS_SOURCES
340342
src/decompositions/self-adjoint-eigen-solver.cpp
341343
src/decompositions/permutation-matrix.cpp
342344
src/decompositions/simplicial-llt-solver.cpp
343-
src/decompositions/simplicial-ldlt-solver.cpp)
345+
src/decompositions/simplicial-ldlt-solver.cpp
346+
src/decompositions/incomplete-cholesky.cpp
347+
src/decompositions/incomplete-lut.cpp)
344348

345349
if(BUILD_WITH_CHOLMOD_SUPPORT)
346350
list(APPEND ${PROJECT_NAME}_DECOMPOSITIONS_SOURCES

‎include/eigenpy/decompositions/sparse/IncompleteCholesky.hpp‎

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decompositions_sparse_incomplete_cholesky_hpp__
6+
#define __eigenpy_decompositions_sparse_incomplete_cholesky_hpp__
7+
8+
#include "eigenpy/eigenpy.hpp"
9+
#include "eigenpy/utils/scalar-name.hpp"
10+
11+
namespace eigenpy {
12+
13+
template <typename _MatrixType>
14+
struct IncompleteCholeskyVisitor : public boost::python::def_visitor<
15+
IncompleteCholeskyVisitor<_MatrixType>> {
16+
typedef _MatrixType MatrixType;
17+
typedef typename MatrixType::Scalar Scalar;
18+
typedef typename MatrixType::RealScalar RealScalar;
19+
typedef Eigen::IncompleteCholesky<Scalar> Solver;
20+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
21+
DenseVectorXs;
22+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
23+
MatrixType::Options>
24+
DenseMatrixXs;
25+
26+
typedef Eigen::SparseMatrix<Scalar, Eigen::ColMajor> FactorType;
27+
typedef Eigen::Matrix<RealScalar, Eigen::Dynamic, 1> VectorRx;
28+
typedef Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>
29+
PermutationType;
30+
31+
template <class PyClass>
32+
void visit(PyClass& cl) const {
33+
cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
34+
.def(bp::init<MatrixType>(bp::args("self", "matrix"),
35+
"Constructs and performs the LDLT "
36+
"factorization from a given matrix."))
37+
38+
.def("rows", &Solver::rows, bp::arg("self"),
39+
"Returns the number of rows of the matrix.")
40+
.def("cols", &Solver::cols, bp::arg("self"),
41+
"Returns the number of cols of the matrix.")
42+
43+
.def("info", &Solver::info, bp::arg("self"),
44+
"Reports whether previous computation was successful.")
45+
46+
.def("setInitialShift", &Solver::setInitialShift,
47+
bp::args("self", "shift"), "Set the initial shift parameter.")
48+
49+
.def(
50+
"analyzePattern",
51+
+[](Solver& self, const MatrixType& amat) {
52+
self.analyzePattern(amat);
53+
},
54+
bp::arg("matrix"))
55+
.def(
56+
"factorize",
57+
+[](Solver& self, const MatrixType& amat) { self.factorize(amat); },
58+
bp::arg("matrix"))
59+
.def(
60+
"compute",
61+
+[](Solver& self, const MatrixType& amat) { self.compute(amat); },
62+
bp::arg("matrix"))
63+
64+
.def(
65+
"matrixL",
66+
+[](Solver& self) -> const FactorType& { return self.matrixL(); },
67+
bp::return_value_policy<bp::copy_const_reference>())
68+
.def(
69+
"scalingS",
70+
+[](Solver& self) -> const VectorRx& { return self.scalingS(); },
71+
bp::return_value_policy<bp::copy_const_reference>())
72+
.def(
73+
"permutationP",
74+
+[](Solver& self) -> const PermutationType& {
75+
return self.permutationP();
76+
},
77+
bp::return_value_policy<bp::copy_const_reference>())
78+
79+
.def(
80+
"solve",
81+
+[](Solver const& self, const Eigen::Ref<DenseVectorXs const>& b)
82+
-> DenseVectorXs { return self.solve(b); },
83+
bp::arg("b"),
84+
"Returns the solution x of A x = b using the current decomposition "
85+
"of A, where b is a right hand side vector.")
86+
.def(
87+
"solve",
88+
+[](Solver const& self, const Eigen::Ref<DenseMatrixXs const>& B)
89+
-> DenseMatrixXs { return self.solve(B); },
90+
bp::arg("b"),
91+
"Returns the solution X of A X = B using the current decomposition "
92+
"of A where B is a right hand side matrix.")
93+
.def(
94+
"solve",
95+
+[](Solver const& self, const MatrixType& B) -> MatrixType {
96+
DenseMatrixXs B_dense = DenseMatrixXs(B);
97+
DenseMatrixXs X_dense = self.solve(B_dense);
98+
return MatrixType(X_dense.sparseView());
99+
},
100+
bp::arg("b"),
101+
"Returns the solution X of A X = B using the current decomposition "
102+
"of A where B is a right hand side matrix.");
103+
}
104+
105+
static void expose() {
106+
static const std::string classname =
107+
"IncompleteCholesky_" + scalar_name<Scalar>::shortname();
108+
expose(classname);
109+
}
110+
111+
static void expose(const std::string& name) {
112+
bp::class_<Solver, boost::noncopyable>(name.c_str(), "Incomplete Cholesky.",
113+
bp::no_init)
114+
.def(IncompleteCholeskyVisitor())
115+
.def(IdVisitor<Solver>());
116+
}
117+
};
118+
119+
} // namespace eigenpy
120+
121+
#endif // ifndef __eigenpy_decompositions_sparse_incomplete_cholesky_hpp__

‎include/eigenpy/decompositions/sparse/IncompleteLUT.hpp‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#ifndef __eigenpy_decompositions_sparse_incomplete_lut_hpp__
6+
#define __eigenpy_decompositions_sparse_incomplete_lut_hpp__
7+
8+
#include "eigenpy/eigenpy.hpp"
9+
#include "eigenpy/utils/scalar-name.hpp"
10+
11+
namespace eigenpy {
12+
13+
template <typename _MatrixType>
14+
struct IncompleteLUTVisitor
15+
: public boost::python::def_visitor<
16+
IncompleteLUTVisitor<_MatrixType>> {
17+
18+
typedef _MatrixType MatrixType;
19+
typedef typename MatrixType::Scalar Scalar;
20+
typedef typename MatrixType::RealScalar RealScalar;
21+
typedef Eigen::IncompleteLUT<Scalar> Solver;
22+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1, MatrixType::Options>
23+
DenseVectorXs;
24+
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic,
25+
MatrixType::Options>
26+
DenseMatrixXs;
27+
28+
template <class PyClass>
29+
void visit(PyClass &cl) const {
30+
cl.def(bp::init<>(bp::arg("self"), "Default constructor"))
31+
.def(bp::init<MatrixType>(bp::args("self", "matrix"),
32+
"Constructs and performs the LDLT "
33+
"factorization from a given matrix."))
34+
.def(bp::init<const MatrixType&, RealScalar, int>(
35+
(bp::arg("matrix"),
36+
bp::arg("droptol") = Eigen::NumTraits<Scalar>::dummy_precision(),
37+
bp::arg("fillfactor") = 10),
38+
"Constructs an incomplete LU factorization from a given matrix."))
39+
40+
.def("rows", &Solver::rows, bp::arg("self"),
41+
"Returns the number of rows of the matrix.")
42+
.def("cols", &Solver::cols, bp::arg("self"),
43+
"Returns the number of cols of the matrix.")
44+
45+
.def("info", &Solver::info, bp::arg("self"),
46+
"Reports whether previous computation was successful.")
47+
48+
.def(
49+
"analyzePattern",
50+
+[](Solver& self, const MatrixType& amat) {
51+
self.analyzePattern(amat);
52+
},
53+
bp::arg("matrix"))
54+
.def(
55+
"factorize",
56+
+[](Solver& self, const MatrixType& amat) { self.factorize(amat); },
57+
bp::arg("matrix"))
58+
.def(
59+
"compute",
60+
+[](Solver& self, const MatrixType& amat) { self.compute(amat); },
61+
bp::arg("matrix"))
62+
63+
.def("setDroptol", &Solver::setDroptol, bp::arg("self"))
64+
.def("setFillfactor", &Solver::setFillfactor, bp::arg("self"))
65+
66+
.def(
67+
"solve",
68+
+[](Solver const& self, const Eigen::Ref<DenseVectorXs const>& b)
69+
-> DenseVectorXs { return self.solve(b); },
70+
bp::arg("b"),
71+
"Returns the solution x of A x = b using the current decomposition "
72+
"of A, where b is a right hand side vector.")
73+
.def(
74+
"solve",
75+
+[](Solver const& self, const Eigen::Ref<DenseMatrixXs const>& B)
76+
-> DenseMatrixXs { return self.solve(B); },
77+
bp::arg("b"),
78+
"Returns the solution X of A X = B using the current decomposition "
79+
"of A where B is a right hand side matrix.")
80+
.def(
81+
"solve",
82+
+[](Solver const& self, const MatrixType& B) -> MatrixType {
83+
DenseMatrixXs B_dense = DenseMatrixXs(B);
84+
DenseMatrixXs X_dense = self.solve(B_dense);
85+
return MatrixType(X_dense.sparseView());
86+
},
87+
bp::arg("b"),
88+
"Returns the solution X of A X = B using the current decomposition "
89+
"of A where B is a right hand side matrix.")
90+
;
91+
}
92+
93+
static void expose() {
94+
static const std::string classname =
95+
"IncompleteLUT_" + scalar_name<Scalar>::shortname();
96+
expose(classname);
97+
}
98+
99+
static void expose(const std::string &name) {
100+
bp::class_<Solver, boost::noncopyable>(
101+
name.c_str(),
102+
"Incomplete LU factorization with dual-threshold strategy.",
103+
bp::no_init)
104+
.def(IncompleteLUTVisitor())
105+
.def(IdVisitor<Solver>());
106+
}
107+
};
108+
109+
} // namespace eigenpy
110+
111+
#endif // ifndef __eigenpy_decompositions_sparse_incomplete_lut_hpp__

‎src/decompositions/decompositions.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@ void exposeQRSolvers();
2626
void exposeMINRESSolver();
2727
void exposeSimplicialLLTSolver();
2828
void exposeSimplicialLDLTSolver();
29+
void exposeIncompleteCholesky();
30+
void exposeIncompleteLUT();
2931
void exposeSparseLUSolver();
3032
void exposeSparseQRSolver();
3133
void exposePermutationMatrix();
@@ -72,6 +74,8 @@ void exposeDecompositions() {
7274
exposeSimplicialLDLTSolver();
7375
exposeSparseLUSolver();
7476
exposeSparseQRSolver();
77+
exposeIncompleteCholesky();
78+
exposeIncompleteLUT();
7579

7680
exposePermutationMatrix();
7781

‎src/decompositions/incomplete-cholesky.cpp‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#include "eigenpy/decompositions/sparse/IncompleteCholesky.hpp"
6+
7+
namespace eigenpy {
8+
void exposeIncompleteCholesky() {
9+
using namespace Eigen;
10+
typedef SparseMatrix<double, ColMajor> ColMajorSparseMatrix;
11+
IncompleteCholeskyVisitor<ColMajorSparseMatrix>::expose("IncompleteCholesky");
12+
}
13+
} // namespace eigenpy

‎src/decompositions/incomplete-lut.cpp‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
/*
2+
* Copyright 2024 INRIA
3+
*/
4+
5+
#include "eigenpy/decompositions/sparse/IncompleteLUT.hpp"
6+
7+
namespace eigenpy {
8+
void exposeIncompleteLUT() {
9+
using namespace Eigen;
10+
typedef SparseMatrix<double, ColMajor> ColMajorSparseMatrix;
11+
IncompleteLUTVisitor<ColMajorSparseMatrix>::expose("IncompleteLUT");
12+
}
13+
} // namespace eigenpy

‎unittest/python/decompositions/sparse/test_IncompleteCholesky.py‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
import numpy as np
2+
from scipy.sparse import csc_matrix
3+
4+
import eigenpy
5+
6+
dim = 100
7+
rng = np.random.default_rng()
8+
9+
A = rng.random((dim, dim))
10+
A = (A + A.T) * 0.5 + np.diag(5.0 + rng.random(dim))
11+
A = csc_matrix(A)
12+
13+
ichol = eigenpy.IncompleteCholesky(A)
14+
assert ichol.info() == eigenpy.ComputationInfo.Success
15+
assert ichol.rows() == dim
16+
assert ichol.cols() == dim
17+
18+
X = rng.random((dim, 20))
19+
B = A.dot(X)
20+
X_est = ichol.solve(B)
21+
assert isinstance(X_est, np.ndarray)
22+
residual = np.linalg.norm(B - A.dot(X_est)) / np.linalg.norm(B)
23+
assert residual < 0.1
24+
25+
x = rng.random(dim)
26+
b = A.dot(x)
27+
x_est = ichol.solve(b)
28+
assert isinstance(x_est, np.ndarray)
29+
residual = np.linalg.norm(b - A.dot(x_est)) / np.linalg.norm(b)
30+
assert residual < 0.1
31+
32+
X_sparse = csc_matrix(rng.random((dim, 10)))
33+
B_sparse = A.dot(X_sparse).tocsc()
34+
if not B_sparse.has_sorted_indices:
35+
B_sparse.sort_indices()
36+
X_est_sparse = ichol.solve(B_sparse)
37+
assert isinstance(X_est_sparse, csc_matrix)
38+
39+
ichol.analyzePattern(A)
40+
ichol.factorize(A)
41+
ichol.compute(A)
42+
assert ichol.info() == eigenpy.ComputationInfo.Success
43+
44+
L = ichol.matrixL()
45+
S_diag = ichol.scalingS()
46+
perm = ichol.permutationP()
47+
P = perm.toDenseMatrix()
48+
49+
assert isinstance(L, csc_matrix)
50+
assert isinstance(S_diag, np.ndarray)
51+
assert L.shape == (dim, dim)
52+
assert S_diag.shape == (dim,)
53+
54+
L_dense = L.toarray()
55+
upper_part = np.triu(L_dense, k=1)
56+
assert np.allclose(upper_part, 0, atol=1e-12)
57+
58+
assert np.all(S_diag > 0)
59+
60+
S = csc_matrix((S_diag, (range(dim), range(dim))), shape=(dim, dim))
61+
62+
PA = P @ A
63+
PAP = PA @ P.T
64+
SPAP = S @ PAP
65+
SPAPS = SPAP @ S
66+
67+
LLT = L @ L.T
68+
69+
diff = SPAPS - LLT
70+
relative_error = np.linalg.norm(diff.data) / np.linalg.norm(SPAPS.data)
71+
assert relative_error < 0.5

0 commit comments

Comments
 (0)