From 2733b4cb1bf7faefb7d6f7379416132764e3a084 Mon Sep 17 00:00:00 2001 From: Justin Carpentier Date: Thu, 5 Aug 2021 14:28:56 +0200 Subject: [PATCH 1/5] solvers: add MINRES --- CMakeLists.txt | 1 + include/eigenpy/decompositions/minres.hpp | 157 ++++++++++++++++++++++ src/decompositions/decompositions.cpp | 5 +- 3 files changed, 162 insertions(+), 1 deletion(-) create mode 100644 include/eigenpy/decompositions/minres.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 689dc3dba..97b46d60e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -97,6 +97,7 @@ SET(${PROJECT_NAME}_DECOMPOSITIONS_HEADERS include/eigenpy/decompositions/LDLT.hpp include/eigenpy/decompositions/LLT.hpp include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp + include/eigenpy/decompositions/minres.hpp ) SET(${PROJECT_NAME}_HEADERS diff --git a/include/eigenpy/decompositions/minres.hpp b/include/eigenpy/decompositions/minres.hpp new file mode 100644 index 000000000..d2d09a65f --- /dev/null +++ b/include/eigenpy/decompositions/minres.hpp @@ -0,0 +1,157 @@ +/* + * Copyright 2021 INRIA + */ + +#ifndef __eigenpy_decomposition_minres_hpp__ +#define __eigenpy_decomposition_minres_hpp__ + +#include "eigenpy/eigenpy.hpp" + +#include +#include + +#include "eigenpy/utils/scalar-name.hpp" + +namespace eigenpy +{ + template + struct IterativeSolverBaseVisitor + : public boost::python::def_visitor< IterativeSolverBaseVisitor<_Solver> > + { + typedef _Solver Solver; + typedef typename Solver::MatrixType MatrixType; + typedef typename Solver::Preconditioner Preconditioner; + typedef typename Solver::Scalar Scalar; + typedef typename Solver::RealScalar RealScalar; + + typedef Eigen::Matrix MatrixXs; + + template + void visit(PyClass& cl) const + { + namespace bp = boost::python; + cl + .def("analyzePattern", + (Solver & (Solver::*)(const Eigen::EigenBase & matrix))&Solver::analyzePattern, + bp::args("self","A"), + "Initializes the iterative solver for the sparsity pattern of the matrix A for further solving Ax=b problems.", + bp::return_self<>()) + .def("factorize", + (Solver & (Solver::*)(const Eigen::EigenBase & matrix))&Solver::factorize, + bp::args("self","A"), + "Initializes the iterative solver with the numerical values of the matrix A for further solving Ax=b problems.", + bp::return_self<>()) + .def("compute", + (Solver & (Solver::*)(const Eigen::EigenBase & matrix))&Solver::compute, + bp::args("self","A"), + "Initializes the iterative solver with the matrix A for further solving Ax=b problems.", + bp::return_self<>()) + + .def("rows",&Solver::rows,bp::arg("self"),"Returns the number of rows.") + .def("cols",&Solver::cols,bp::arg("self"),"Returns the number of columns.") + .def("tolerance",&Solver::tolerance,bp::arg("self"), + "Returns the tolerance threshold used by the stopping criteria.") + .def("setTolerance",&Solver::setTolerance,bp::args("self","tolerance"), + "Sets the tolerance threshold used by the stopping criteria.\n" + "This value is used as an upper bound to the relative residual error: |Ax-b|/|b|.\n" + "The default value is the machine precision given by NumTraits::epsilon().", + bp::return_self<>()) + .def("preconditioner",(Preconditioner & (Solver::*)())&Solver::preconditioner,bp::arg("self"), + "Returns a read-write reference to the preconditioner for custom configuration.", + bp::return_internal_reference<>()) + + .def("maxIterations",&Solver::maxIterations,bp::arg("self"), + "Returns the max number of iterations.\n" + "It is either the value setted by setMaxIterations or, by default, twice the number of columns of the matrix.") + .def("setMaxIterations",&Solver::setMaxIterations,bp::args("self","max_iterations"), + "Sets the max number of iterations.\n" + "Default is twice the number of columns of the matrix.", + bp::return_self<>()) + + .def("iterations",&Solver::iterations,bp::arg("self"), + "Returns the number of iterations performed during the last solve.") + .def("error",&Solver::error,bp::arg("self"), + "Returns the tolerance error reached during the last solve.\n" + "It is a close approximation of the true relative residual error |Ax-b|/|b|.") + .def("info",&Solver::error,bp::arg("info"), + "Returns Success if the iterations converged, and NoConvergence otherwise.") + + .def("solveWithGuess",&solveWithGuess,bp::args("self","b","x0"), + "Returns the solution x of A x = b using the current decomposition of A and x0 as an initial solution.") + + .def("solve",&solve,bp::args("self","b"), + "Returns the solution x of A x = b using the current decomposition of A where b is a right hand side matrix or vector.") + ; + } + + private: + + template + static MatrixOrVector1 solveWithGuess(const Solver & self, const MatrixOrVector1 & b, const MatrixOrVector2 & guess) + { + return self.solveWithGuess(b,guess); + } + + template + static MatrixOrVector solve(const Solver & self, const MatrixOrVector & mat_or_vec) + { + MatrixOrVector res = self.solve(mat_or_vec); + return res; + } + }; + + template + struct MINRESSolverVisitor + : public boost::python::def_visitor< MINRESSolverVisitor<_MatrixType> > + { + + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Eigen::Matrix VectorXs; + typedef Eigen::Matrix MatrixXs; + typedef Eigen::MINRES Solver; + + template + void visit(PyClass& cl) const + { + namespace bp = boost::python; + cl + .def(bp::init<>("Default constructor")) + .def(bp::init(bp::arg("matrix"), + "Initialize the solver with matrix A for further Ax=b solving.\n" + "This constructor is a shortcut for the default constructor followed by a call to compute().")) + + .def(IterativeSolverBaseVisitor()) + ; + } + + static void expose() + { + static const std::string classname = "MINRES" + scalar_name::shortname(); + expose(classname); + } + + static void expose(const std::string & name) + { + namespace bp = boost::python; + bp::class_(name.c_str(), + "A minimal residual solver for sparse symmetric problems.\n" + "This class allows to solve for A.x = b sparse linear problems using the MINRES algorithm of Paige and Saunders (1975). The sparse matrix A must be symmetric (possibly indefinite). The vectors x and b can be either dense or sparse.\n" + "The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations and NumTraits::epsilon() for the tolerance.\n", + bp::no_init) + .def(MINRESSolverVisitor()); + } + + private: + + template + static MatrixOrVector solve(const Solver & self, const MatrixOrVector & vec) + { + return self.solve(vec); + } + }; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_decomposition_minres_hpp__ diff --git a/src/decompositions/decompositions.cpp b/src/decompositions/decompositions.cpp index a3bbd36e8..93ed609d6 100644 --- a/src/decompositions/decompositions.cpp +++ b/src/decompositions/decompositions.cpp @@ -1,5 +1,5 @@ /* - * Copyright 2020 INRIA + * Copyright 2020-2021 INRIA */ #include "eigenpy/fwd.hpp" @@ -10,6 +10,7 @@ #include "eigenpy/decompositions/SelfAdjointEigenSolver.hpp" #include "eigenpy/decompositions/LLT.hpp" #include "eigenpy/decompositions/LDLT.hpp" +#include "eigenpy/decompositions/minres.hpp" namespace eigenpy { @@ -22,6 +23,8 @@ namespace eigenpy SelfAdjointEigenSolverVisitor::expose("SelfAdjointEigenSolver"); LLTSolverVisitor::expose("LLT"); LDLTSolverVisitor::expose("LDLT"); + + MINRESSolverVisitor::expose("MINRES"); { bp::enum_("DecompositionOptions") From a82e0fe0a09e9582436ce5cf3beedea60aa332bf Mon Sep 17 00:00:00 2001 From: Justin Carpentier Date: Thu, 5 Aug 2021 14:29:07 +0200 Subject: [PATCH 2/5] test: add MINRES test --- unittest/CMakeLists.txt | 5 ++++- unittest/python/test_MINRES.py | 17 +++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) create mode 100644 unittest/python/test_MINRES.py diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt index ba18c9c60..d8730c43c 100644 --- a/unittest/CMakeLists.txt +++ b/unittest/CMakeLists.txt @@ -1,6 +1,6 @@ # # Copyright (c) 2014-2019 CNRS -# Copyright (c) 2018-2020 INRIA +# Copyright (c) 2018-2021 INRIA # MACRO(ADD_LIB_UNIT_TEST test) @@ -62,3 +62,6 @@ SET_TESTS_PROPERTIES("py-LLT" PROPERTIES DEPENDS ${PYWRAP}) ADD_PYTHON_UNIT_TEST("py-LDLT" "unittest/python/test_LDLT.py" "python/eigenpy;unittest") SET_TESTS_PROPERTIES("py-LDLT" PROPERTIES DEPENDS ${PYWRAP}) + +ADD_PYTHON_UNIT_TEST("py-MINRES" "unittest/python/test_MINRES.py" "python/eigenpy;unittest") +SET_TESTS_PROPERTIES("py-MINRES" PROPERTIES DEPENDS ${PYWRAP}) diff --git a/unittest/python/test_MINRES.py b/unittest/python/test_MINRES.py new file mode 100644 index 000000000..3b1daad79 --- /dev/null +++ b/unittest/python/test_MINRES.py @@ -0,0 +1,17 @@ +import eigenpy + +import numpy as np +import numpy.linalg as la + +dim = 100 +A = np.random.rand(dim,dim) + +A = (A + A.T)*0.5 + np.diag(10. + np.random.rand(dim)) + +minres = eigenpy.MINRES(A) + +X = np.random.rand(dim,20) +B = A.dot(X) +X_est = minres.solve(B) +assert eigenpy.is_approx(X,X_est) +assert eigenpy.is_approx(A.dot(X_est),B) From 7b489df1f11a868e9893363ba6714df7a870f373 Mon Sep 17 00:00:00 2001 From: Justin Carpentier Date: Thu, 5 Aug 2021 14:53:55 +0200 Subject: [PATCH 3/5] core: add missing include Related to Eigen --- include/eigenpy/decompositions/minres.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/eigenpy/decompositions/minres.hpp b/include/eigenpy/decompositions/minres.hpp index d2d09a65f..45123ce2e 100644 --- a/include/eigenpy/decompositions/minres.hpp +++ b/include/eigenpy/decompositions/minres.hpp @@ -7,6 +7,7 @@ #include "eigenpy/eigenpy.hpp" +#include #include #include From 7de35b3d7afa241b53a192d8e101df2b1b112d39 Mon Sep 17 00:00:00 2001 From: Justin Carpentier Date: Thu, 5 Aug 2021 16:14:15 +0200 Subject: [PATCH 4/5] test: relax --- unittest/python/test_MINRES.py | 1 - 1 file changed, 1 deletion(-) diff --git a/unittest/python/test_MINRES.py b/unittest/python/test_MINRES.py index 3b1daad79..80f0d0a0e 100644 --- a/unittest/python/test_MINRES.py +++ b/unittest/python/test_MINRES.py @@ -13,5 +13,4 @@ X = np.random.rand(dim,20) B = A.dot(X) X_est = minres.solve(B) -assert eigenpy.is_approx(X,X_est) assert eigenpy.is_approx(A.dot(X_est),B) From aad265bb197c25f5f4e1e6de226bf74af4908f4a Mon Sep 17 00:00:00 2001 From: Justin Carpentier Date: Thu, 5 Aug 2021 17:27:42 +0200 Subject: [PATCH 5/5] test: change A matrix --- unittest/python/test_MINRES.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/unittest/python/test_MINRES.py b/unittest/python/test_MINRES.py index 80f0d0a0e..67cbfac96 100644 --- a/unittest/python/test_MINRES.py +++ b/unittest/python/test_MINRES.py @@ -4,13 +4,11 @@ import numpy.linalg as la dim = 100 -A = np.random.rand(dim,dim) - -A = (A + A.T)*0.5 + np.diag(10. + np.random.rand(dim)) +A = np.eye(dim) minres = eigenpy.MINRES(A) X = np.random.rand(dim,20) B = A.dot(X) X_est = minres.solve(B) -assert eigenpy.is_approx(A.dot(X_est),B) +assert eigenpy.is_approx(A.dot(X_est),B,1e-6)