decompositions: Completed tests

Author: Lucas-Haubert

include/eigenpy/decompositions/FullPivLU.hpp

@@ -33,8 +33,6 @@ struct FullPivLUSolverVisitor
         .def(bp::init<Eigen::DenseIndex, Eigen::DenseIndex>(
             bp::args("self", "rows", "cols"),
             "Default constructor with memory preallocation"))
-        .def(bp::init<const MatrixType>(bp::args("self", "matrix"),
-                                        "Constructor."))
         .def(bp::init<MatrixType>(
             bp::args("self", "matrix"),
             "Constructs a LU factorization from a given matrix."))
@@ -57,7 +55,7 @@ struct FullPivLUSolverVisitor
              "is the LU decomposition.")
         .def(
             "image",
-            +[](const Solver &self, const MatrixType &mat) -> MatrixType {
+            +[](Solver &self, const MatrixType &mat) -> MatrixType {
               return self.image(mat);
             },
             bp::args("self", "originalMatrix"),
@@ -66,7 +64,7 @@ struct FullPivLUSolverVisitor
             "image (column-space).")
         .def(
             "inverse",
-            +[](const Solver &self) -> MatrixType { return self.inverse(); },
+            +[](Solver &self) -> MatrixType { return self.inverse(); },
             bp::arg("self"),
             "Returns the inverse of the matrix of which *this is the LU "
             "decomposition.")
@@ -76,8 +74,7 @@ struct FullPivLUSolverVisitor
         .def("isSurjective", &Solver::isSurjective, bp::arg("self"))
 
         .def(
-            "kernel",
-            +[](const Solver &self) -> MatrixType { return self.kernel(); },
+            "kernel", +[](Solver &self) -> MatrixType { return self.kernel(); },
             bp::arg("self"),
             "Returns the kernel of the matrix, also called its null-space. "
             "The columns of the returned matrix will form a basis of the "
@@ -90,7 +87,6 @@ struct FullPivLUSolverVisitor
         .def("maxPivot", &Solver::maxPivot, bp::arg("self"))
         .def("nonzeroPivots", &Solver::nonzeroPivots, bp::arg("self"))
 
-        // TODO: Expose so that the return type are convertible to np arrays
         .def("permutationP", &Solver::permutationP, bp::arg("self"),
              "Returns the permutation P.",
              bp::return_value_policy<bp::copy_const_reference>())

include/eigenpy/decompositions/GeneralizedEigenSolver.hpp

@@ -34,9 +34,10 @@ struct GeneralizedEigenSolverVisitor
 
         .def("eigenvectors", &Solver::eigenvectors, bp::arg("self"),
              "Returns an expression of the computed generalized eigenvectors. ")
-        // TODO: Expose so that the return type are convertible to np arrays
-        .def("eigenvalues", &Solver::eigenvalues, bp::arg("self"),
-             "Returns an expression of the computed generalized eigenvalues. ")
+        .def(
+            "eigenvalues",
+            +[](const Solver& c) { return c.eigenvalues().eval(); },
+            "Returns the computed generalized eigenvalues.")
 
         .def("alphas", &Solver::alphas, bp::arg("self"),
              "Returns the vectors containing the alpha values. ")

include/eigenpy/decompositions/HessenbergDecomposition.hpp

@@ -44,9 +44,14 @@ struct HessenbergDecompositionVisitor
              bp::arg("self"), "Returns the Householder coefficients. ",
              bp::return_value_policy<bp::copy_const_reference>())
 
-        // TODO: Expose so that the return type are convertible to np arrays
-        // matrixH
-        // matrixQ
+        .def(
+            "matrixQ",
+            +[](const Solver& c) -> MatrixType { return c.matrixQ(); },
+            "Reconstructs the orthogonal matrix Q in the decomposition.")
+        .def(
+            "matrixH",
+            +[](const Solver& c) -> MatrixType { return c.matrixH(); },
+            "Constructs the Hessenberg matrix H in the decomposition.")
 
         .def("packedMatrix", &Solver::packedMatrix, bp::arg("self"),
              "Returns the internal representation of the decomposition. ",

include/eigenpy/decompositions/JacobiSVD.hpp

@@ -15,11 +15,10 @@
 
 namespace eigenpy {
 
-template <typename _MatrixType>
+template <typename JacobiSVD>
 struct JacobiSVDVisitor
-    : public boost::python::def_visitor<JacobiSVDVisitor<_MatrixType>> {
-  typedef _MatrixType MatrixType;
-  typedef Eigen::JacobiSVD<MatrixType> Solver;
+    : public boost::python::def_visitor<JacobiSVDVisitor<JacobiSVD>> {
+  typedef typename JacobiSVD::MatrixType MatrixType;
   typedef typename MatrixType::Scalar Scalar;
 
   template <class PyClass>
@@ -46,28 +45,29 @@ struct JacobiSVDVisitor
             "and the constructor. If possible, prefer using the Options "
             "template parameter."))
 
-        .def("cols", &Solver::cols, bp::arg("self"),
+        .def("cols", &JacobiSVD::cols, bp::arg("self"),
              "Returns the number of columns. ")
         .def("compute",
-             (Solver & (Solver::*)(const MatrixType &matrix)) & Solver::compute,
+             (JacobiSVD & (JacobiSVD::*)(const MatrixType &matrix)) &
+                 JacobiSVD::compute,
              bp::args("self", "matrix"),
              "Method performing the decomposition of given matrix. Computes "
              "Thin/Full "
              "unitaries U/V if specified using the Options template parameter "
              "or the class constructor. ",
              bp::return_self<>())
         .def("compute",
-             (Solver & (Solver::*)(const MatrixType &matrix,
-                                   unsigned int computationOptions)) &
-                 Solver::compute,
+             (JacobiSVD & (JacobiSVD::*)(const MatrixType &matrix,
+                                         unsigned int computationOptions)) &
+                 JacobiSVD::compute,
              bp::args("self", "matrix", "computationOptions"),
              "Method performing the decomposition of given matrix, as "
              "specified by the computationOptions parameter.  ",
              bp::return_self<>())
-        .def("rows", &Solver::rows, bp::arg("self"),
+        .def("rows", &JacobiSVD::rows, bp::arg("self"),
              "Returns the number of rows. . ")
 
-        .def(SVDBaseVisitor<Solver>());
+        .def(SVDBaseVisitor<JacobiSVD>());
   }
 
   static void expose() {
@@ -77,7 +77,7 @@ struct JacobiSVDVisitor
   }
 
   static void expose(const std::string &name) {
-    bp::class_<Solver, boost::noncopyable>(
+    bp::class_<JacobiSVD, boost::noncopyable>(
         name.c_str(),
         "Two-sided Jacobi SVD decomposition of a rectangular matrix. \n\n"
         "SVD decomposition consists in decomposing any n-by-p matrix A as a "
@@ -109,7 +109,7 @@ struct JacobiSVDVisitor
         "solving.",
         bp::no_init)
         .def(JacobiSVDVisitor())
-        .def(IdVisitor<Solver>());
+        .def(IdVisitor<JacobiSVD>());
   }
 };
 

include/eigenpy/decompositions/PartialPivLU.hpp

@@ -33,8 +33,6 @@ struct PartialPivLUSolverVisitor : public boost::python::def_visitor<
         .def(bp::init<Eigen::DenseIndex>(
             bp::args("self", "size"),
             "Default constructor with memory preallocation"))
-        .def(bp::init<const MatrixType>(bp::args("self", "matrix"),
-                                        "Constructor."))
         .def(bp::init<MatrixType>(
             bp::args("self", "matrix"),
             "Constructs a LU factorization from a given matrix."))
@@ -44,6 +42,13 @@ struct PartialPivLUSolverVisitor : public boost::python::def_visitor<
         .def("determinant", &Solver::determinant, bp::arg("self"),
              "Returns the determinant of the matrix of which *this is the LU "
              "decomposition.")
+        .def(
+            "compute",
+            (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
+                Solver::compute,
+            bp::args("self", "matrix"),
+            "Computes the LU factorization of given matrix.",
+            bp::return_self<>())
         .def(
             "inverse",
             +[](const Solver &self) -> MatrixType { return self.inverse(); },
@@ -54,7 +59,6 @@ struct PartialPivLUSolverVisitor : public boost::python::def_visitor<
              "Returns the LU decomposition matrix.",
              bp::return_internal_reference<>())
 
-        // TODO: Expose so that the return type are convertible to np arrays
         .def("permutationP", &Solver::permutationP, bp::arg("self"),
              "Returns the permutation P.",
              bp::return_value_policy<bp::copy_const_reference>())

include/eigenpy/decompositions/Tridiagonalization.hpp

@@ -20,9 +20,10 @@ struct TridiagonalizationVisitor : public boost::python::def_visitor<
   typedef _MatrixType MatrixType;
   typedef typename MatrixType::Scalar Scalar;
   typedef Eigen::Tridiagonalization<MatrixType> Solver;
+  typedef Eigen::VectorXd VectorType;
 
   template <class PyClass>
-  void visit(PyClass& cl) const {
+  void visit(PyClass &cl) const {
     cl.def(
           bp::init<Eigen::DenseIndex>(bp::arg("size"), "Default constructor. "))
         .def(bp::init<MatrixType>(bp::arg("matrix"),
@@ -31,7 +32,7 @@ struct TridiagonalizationVisitor : public boost::python::def_visitor<
 
         .def(
             "compute",
-            (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType>& matrix)) &
+            (Solver & (Solver::*)(const Eigen::EigenBase<MatrixType> &matrix)) &
                 Solver::compute,
             bp::args("self", "matrix"),
             "Computes tridiagonal decomposition of given matrix. ",
@@ -43,18 +44,28 @@ struct TridiagonalizationVisitor : public boost::python::def_visitor<
              "Returns the internal representation of the decomposition. ",
              bp::return_value_policy<bp::copy_const_reference>())
 
-        // TODO: Expose so that the return type are convertible to np arrays
-        .def("matrixQ", &Solver::matrixQ, bp::arg("self"),
-             "Returns the unitary matrix Q in the decomposition. ")
-        .def("matrixT", &Solver::matrixT, bp::arg("self"),
-             "Returns the unitary matrix T in the decomposition. ")
+        .def(
+            "matrixQ",
+            +[](const Solver &c) -> MatrixType { return c.matrixQ(); },
+            "Returns the unitary matrix Q in the decomposition.")
+        .def(
+            "matrixT", +[](Solver &c) -> MatrixType { return c.matrixT(); },
+            "Returns an expression of the tridiagonal matrix T in the "
+            "decomposition.")
 
-        .def("diagonal", &Solver::diagonal, bp::arg("self"),
-             "Returns the diagonal of the tridiagonal matrix T in the "
-             "decomposition. ")
-        .def("subDiagonal", &Solver::subDiagonal, bp::arg("self"),
-             "Returns the subdiagonal of the tridiagonal matrix T in the "
-             "decomposition.");
+        .def(
+            "diagonal",
+            +[](const Solver &c) -> VectorType { return c.diagonal(); },
+            bp::arg("self"),
+            "Returns the diagonal of the tridiagonal matrix T in the "
+            "decomposition. ")
+
+        .def(
+            "subDiagonal",
+            +[](const Solver &c) -> VectorType { return c.subDiagonal(); },
+            bp::arg("self"),
+            "Returns the subdiagonal of the tridiagonal matrix T in the "
+            "decomposition.");
   }
 
   static void expose() {
@@ -63,7 +74,7 @@ struct TridiagonalizationVisitor : public boost::python::def_visitor<
     expose(classname);
   }
 
-  static void expose(const std::string& name) {
+  static void expose(const std::string &name) {
     bp::class_<Solver>(name.c_str(), bp::no_init)
         .def(TridiagonalizationVisitor())
         .def(IdVisitor<Solver>());

include/eigenpy/decompositions/sparse/LU.hpp

@@ -49,6 +49,16 @@ struct SparseLUVisitor : public boost::python::def_visitor<
              "Returns the natural log of the absolute value of the determinant "
              "of the matrix of which *this is the QR decomposition")
 
+        .def("rows", &Solver::rows, bp::arg("self"),
+             "Returns the number of rows of the matrix.")
+        .def("cols", &Solver::cols, bp::arg("self"),
+             "Returns the number of cols of the matrix.")
+
+        .def("nnzL", &Solver::nnzL, bp::arg("self"),
+             "The number of non zero elements in L")
+        .def("nnzU", &Solver::nnzU, bp::arg("self"),
+             "The number of non zero elements in L")
+
         .def("colsPermutation", &Solver::colsPermutation, bp::arg("self"),
              "Returns a reference to the column matrix permutation PTc such "
              "that Pr A PTc = LU.",
@@ -71,12 +81,6 @@ struct SparseLUVisitor : public boost::python::def_visitor<
              "The given matrix must has the same sparcity than the matrix on "
              "which the symbolic decomposition has been performed.")
 
-        // TODO: Expose so that the return type are convertible to np arrays
-        // transpose()
-        // adjoint()
-        // matrixL()
-        // matrixU()
-
         .def("isSymmetric", &Solver::isSymmetric, bp::args("self", "sym"),
              "Indicate that the pattern of the input matrix is symmetric. ")
 

include/eigenpy/decompositions/sparse/QR.hpp

@@ -57,10 +57,6 @@ struct SparseQRVisitor : public boost::python::def_visitor<
              "The given matrix must has the same sparcity than the matrix on "
              "which the symbolic decomposition has been performed.")
 
-        // TODO: Expose so that the return type are convertible to np arrays
-        // matrixQ
-        // matrixR
-
         .def("colsPermutation", &Solver::colsPermutation, bp::arg("self"),
              "Returns a reference to the column matrix permutation PTc such "
              "that Pr A PTc = LU.",

src/decompositions/jacobisvd-solver.cpp

@@ -7,6 +7,23 @@
 namespace eigenpy {
 void exposeJacobiSVDSolver() {
   using namespace Eigen;
-  JacobiSVDVisitor<MatrixXd>::expose("JacobiSVD");
+  using Eigen::JacobiSVD;
+
+  using Eigen::ColPivHouseholderQRPreconditioner;
+  using Eigen::FullPivHouseholderQRPreconditioner;
+  using Eigen::HouseholderQRPreconditioner;
+  using Eigen::NoQRPreconditioner;
+
+  using ColPivHhJacobiSVD =
+      JacobiSVD<MatrixXd, ColPivHouseholderQRPreconditioner>;
+  using FullPivHhJacobiSVD =
+      JacobiSVD<MatrixXd, FullPivHouseholderQRPreconditioner>;
+  using HhJacobiSVD = JacobiSVD<MatrixXd, HouseholderQRPreconditioner>;
+  using NoPrecondJacobiSVD = JacobiSVD<MatrixXd, NoQRPreconditioner>;
+
+  JacobiSVDVisitor<ColPivHhJacobiSVD>::expose("ColPivHhJacobiSVD");
+  JacobiSVDVisitor<FullPivHhJacobiSVD>::expose("FullPivHhJacobiSVD");
+  JacobiSVDVisitor<HhJacobiSVD>::expose("HhJacobiSVD");
+  JacobiSVDVisitor<NoPrecondJacobiSVD>::expose("NoPrecondJacobiSVD");
 }
 }  // namespace eigenpy

unittest/python/decompositions/sparse/test_SparseLU.py

@@ -1,48 +1,38 @@
 import numpy as np
-import scipy
-from scipy.sparse import csc_matrix
+import scipy.sparse as spa
 
 import eigenpy
 
 dim = 100
 rng = np.random.default_rng()
 
-A = rng.random((dim, dim))
-A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
-A = csc_matrix(A)
+A_fac = spa.random(dim, dim, density=0.25, random_state=rng)
+A = A_fac.T @ A_fac
+A += spa.diags(10.0 * rng.standard_normal(dim) ** 2)
+A = A.tocsc(True)
+A.check_format()
 
-sparselu = eigenpy.SparseLU(A)
+splu = eigenpy.SparseLU(A)
 
-assert sparselu.info() == eigenpy.ComputationInfo.Success
+assert splu.info() == eigenpy.ComputationInfo.Success
 
-# Solve
 X = rng.random((dim, 20))
 B = A.dot(X)
-X_est = sparselu.solve(B)
+X_est = splu.solve(B)
+assert isinstance(X_est, np.ndarray)
 assert eigenpy.is_approx(X, X_est)
 assert eigenpy.is_approx(A.dot(X_est), B)
 
-X_sparse = scipy.sparse.random(dim, 10)
-B_sparse = A.dot(X_sparse)
-B_sparse = B_sparse.tocsc(True)
+splu.analyzePattern(A)
+splu.factorize(A)
 
+X_sparse = spa.random(dim, 10, random_state=rng)
+B_sparse = A.dot(X_sparse)
+B_sparse: spa.csc_matrix = B_sparse.tocsc(True)
 if not B_sparse.has_sorted_indices:
     B_sparse.sort_indices()
 
-X_est = sparselu.solve(B_sparse)
+X_est = splu.solve(B_sparse)
+assert isinstance(X_est, spa.csc_matrix)
 assert eigenpy.is_approx(X_est.toarray(), X_sparse.toarray())
 assert eigenpy.is_approx(A.dot(X_est.toarray()), B_sparse.toarray())
-
-# Others
-det = sparselu.determinant()
-sign_det = sparselu.signDeterminant()
-abs_det = sparselu.absDeterminant()
-log_abs_det = sparselu.logAbsDeterminant()
-
-sparselu.analyzePattern(A)
-sparselu.factorize(A)
-
-cols_permutation = sparselu.colsPermutation()
-rows_permutation = sparselu.rowsPermutation()
-
-sparselu.setPivotThreshold(1e-8)