is_approx: Put is_approx and set the precision

Author: Lucas-Haubert

include/eigenpy/decompositions/SelfAdjointEigenSolver.hpp

@@ -21,6 +21,7 @@ struct SelfAdjointEigenSolverVisitor
   typedef _MatrixType MatrixType;
   typedef typename MatrixType::Scalar Scalar;
   typedef Eigen::SelfAdjointEigenSolver<MatrixType> Solver;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> VectorType;
 
   template <class PyClass>
   void visit(PyClass& cl) const {
@@ -32,12 +33,16 @@ struct SelfAdjointEigenSolverVisitor
             bp::args("self", "matrix", "options"),
             "Computes eigendecomposition of given matrix"))
 
-        .def("eigenvalues", &Solver::eigenvalues, bp::arg("self"),
-             "Returns the eigenvalues of given matrix.",
-             bp::return_internal_reference<>())
-        .def("eigenvectors", &Solver::eigenvectors, bp::arg("self"),
-             "Returns the eigenvectors of given matrix.",
-             bp::return_internal_reference<>())
+        .def(
+            "eigenvalues",
+            +[](Solver& c) -> const VectorType& { return c.eigenvalues(); },
+            "Returns the eigenvalues of given matrix.",
+            bp::return_internal_reference<>())
+        .def(
+            "eigenvectors",
+            +[](Solver& c) -> const MatrixType& { return c.eigenvectors(); },
+            "Returns the eigenvectors of given matrix.",
+            bp::return_internal_reference<>())
 
         .def("compute",
              &SelfAdjointEigenSolverVisitor::compute_proxy<MatrixType>,

unittest/python/test_ComplexSchur.py

@@ -13,8 +13,8 @@
 T = cs.matrixT()
 
 A_complex = A.astype(complex)
-assert eigenpy.is_approx(A_complex, U @ T @ U.conj().T, 1e-10)
-assert eigenpy.is_approx(U @ U.conj().T, np.eye(dim), 1e-10)
+assert eigenpy.is_approx(A_complex, U @ T @ U.conj().T)
+assert eigenpy.is_approx(U @ U.conj().T, np.eye(dim))
 
 for row in range(1, dim):
     for col in range(row):
@@ -33,20 +33,20 @@
 T2 = cs2.matrixT()
 U2 = cs2.matrixU()
 
-assert eigenpy.is_approx(T1, T2, 1e-12)
-assert eigenpy.is_approx(U1, U2, 1e-12)
+assert eigenpy.is_approx(T1, T2)
+assert eigenpy.is_approx(U1, U2)
 
 cs_no_u = eigenpy.ComplexSchur(A, False)
 assert cs_no_u.info() == eigenpy.ComputationInfo.Success
 T_no_u = cs_no_u.matrixT()
 
-assert eigenpy.is_approx(T, T_no_u, 1e-12)
+assert eigenpy.is_approx(T, T_no_u)
 
 cs_compute_no_u = eigenpy.ComplexSchur(dim)
 result_no_u = cs_compute_no_u.compute(A, False)
 assert result_no_u.info() == eigenpy.ComputationInfo.Success
 T_compute_no_u = cs_compute_no_u.matrixT()
-assert eigenpy.is_approx(T, T_compute_no_u, 1e-12)
+assert eigenpy.is_approx(T, T_compute_no_u)
 
 cs_iter = eigenpy.ComplexSchur(dim)
 cs_iter.setMaxIterations(30 * dim)  # m_maxIterationsPerRow * size
@@ -56,8 +56,8 @@
 
 T_iter = cs_iter.matrixT()
 U_iter = cs_iter.matrixU()
-assert eigenpy.is_approx(T, T_iter, 1e-12)
-assert eigenpy.is_approx(U, U_iter, 1e-12)
+assert eigenpy.is_approx(T, T_iter)
+assert eigenpy.is_approx(U, U_iter)
 
 cs_few_iter = eigenpy.ComplexSchur(dim)
 cs_few_iter.setMaxIterations(1)
@@ -74,8 +74,8 @@
 U_triangular = cs_triangular.matrixU()
 
 A_triangular_complex = A_triangular.astype(complex)
-assert eigenpy.is_approx(T_triangular, A_triangular_complex, 1e-10)
-assert eigenpy.is_approx(U_triangular, np.eye(dim, dtype=complex), 1e-10)
+assert eigenpy.is_approx(T_triangular, A_triangular_complex)
+assert eigenpy.is_approx(U_triangular, np.eye(dim, dtype=complex))
 
 hess = eigenpy.HessenbergDecomposition(A)
 H = hess.matrixH()
@@ -89,9 +89,7 @@
 U_from_hess = cs_from_hess.matrixU()
 
 A_complex = A.astype(complex)
-assert eigenpy.is_approx(
-    A_complex, U_from_hess @ T_from_hess @ U_from_hess.conj().T, 1e-10
-)
+assert eigenpy.is_approx(A_complex, U_from_hess @ T_from_hess @ U_from_hess.conj().T)
 
 cs1_id = eigenpy.ComplexSchur(dim)
 cs2_id = eigenpy.ComplexSchur(dim)

unittest/python/test_GeneralizedEigenSolver.py

@@ -35,8 +35,8 @@
 
     Av = A @ v
     lambda_Bv = lambda_k * (B @ v)
-    assert eigenpy.is_approx(Av.real, lambda_Bv.real, 1e-10)
-    assert eigenpy.is_approx(Av.imag, lambda_Bv.imag, 1e-10)
+    assert eigenpy.is_approx(Av.real, lambda_Bv.real, 1e-6)
+    assert eigenpy.is_approx(Av.imag, lambda_Bv.imag, 1e-6)
 
 for k in range(dim):
     v = eigenvectors[:, k]
@@ -45,8 +45,8 @@
 
     alpha_Bv = alpha * (B @ v)
     beta_Av = beta * (A @ v)
-    assert eigenpy.is_approx(alpha_Bv.real, beta_Av.real, 1e-10)
-    assert eigenpy.is_approx(alpha_Bv.imag, beta_Av.imag, 1e-10)
+    assert eigenpy.is_approx(alpha_Bv.real, beta_Av.real, 1e-6)
+    assert eigenpy.is_approx(alpha_Bv.imag, beta_Av.imag, 1e-6)
 
 for k in range(dim):
     if abs(betas[k]) > 1e-12:

unittest/python/test_GeneralizedSelfAdjointEigenSolver.py

@@ -45,10 +45,10 @@ def test_generalized_selfadjoint_eigensolver(options):
                 lam = D[i]
                 Av = A @ v
                 lam_Bv = lam * (B @ v)
-                assert eigenpy.is_approx(Av, lam_Bv, 1e-10)
+                assert eigenpy.is_approx(Av, lam_Bv, 1e-6)
 
             VT_B_V = V.T @ B @ V
-            assert eigenpy.is_approx(VT_B_V, np.eye(dim), 1e-10)
+            assert eigenpy.is_approx(VT_B_V, np.eye(dim), 1e-6)
 
         elif options & eigenpy.DecompositionOptions.ABx_lx:
             AB = A @ B
@@ -57,7 +57,7 @@ def test_generalized_selfadjoint_eigensolver(options):
                 lam = D[i]
                 ABv = AB @ v
                 lam_v = lam * v
-                assert np.allclose(ABv, lam_v)
+                assert eigenpy.is_approx(ABv, lam_v, 1e-6)
 
         elif options & eigenpy.DecompositionOptions.BAx_lx:
             BA = B @ A
@@ -66,7 +66,7 @@ def test_generalized_selfadjoint_eigensolver(options):
                 lam = D[i]
                 BAv = BA @ v
                 lam_v = lam * v
-                assert np.allclose(BAv, lam_v)
+                assert eigenpy.is_approx(BAv, lam_v, 1e-6)
 
     _gsaes_compute = gsaes.compute(A, B)
     _gsaes_compute_options = gsaes.compute(A, B, options)

unittest/python/test_HessenbergDecomposition.py

@@ -15,7 +15,7 @@
     A_reconstructed = Q @ H @ Q.conj().T
 else:
     A_reconstructed = Q @ H @ Q.T
-assert eigenpy.is_approx(A, A_reconstructed, 1e-10)
+assert eigenpy.is_approx(A, A_reconstructed)
 
 for row in range(2, dim):
     for col in range(row - 1):
@@ -25,7 +25,7 @@
     QQ_conj = Q @ Q.conj().T
 else:
     QQ_conj = Q @ Q.T
-assert eigenpy.is_approx(QQ_conj, np.eye(dim), 1e-10)
+assert eigenpy.is_approx(QQ_conj, np.eye(dim))
 
 A_test = rng.random((dim, dim))
 hess1 = eigenpy.HessenbergDecomposition(dim)
@@ -37,8 +37,8 @@
 Q1 = hess1.matrixQ()
 Q2 = hess2.matrixQ()
 
-assert eigenpy.is_approx(H1, H2, 1e-12)
-assert eigenpy.is_approx(Q1, Q2, 1e-12)
+assert eigenpy.is_approx(H1, H2)
+assert eigenpy.is_approx(Q1, Q2)
 
 hCoeffs = hess.householderCoefficients()
 packed = hess.packedMatrix()

unittest/python/test_RealQZ.py

@@ -15,11 +15,11 @@
 Z = realqz.matrixZ()
 T = realqz.matrixT()
 
-assert eigenpy.is_approx(A, Q @ S @ Z, 1e-10)
-assert eigenpy.is_approx(B, Q @ T @ Z, 1e-10)
+assert eigenpy.is_approx(A, Q @ S @ Z)
+assert eigenpy.is_approx(B, Q @ T @ Z)
 
-assert eigenpy.is_approx(Q @ Q.T, np.eye(dim), 1e-10)
-assert eigenpy.is_approx(Z @ Z.T, np.eye(dim), 1e-10)
+assert eigenpy.is_approx(Q @ Q.T, np.eye(dim))
+assert eigenpy.is_approx(Z @ Z.T, np.eye(dim))
 
 for i in range(dim):
     for j in range(i):
@@ -62,8 +62,8 @@
 S_without = realqz_without_qz.matrixS()
 T_without = realqz_without_qz.matrixT()
 
-assert eigenpy.is_approx(S_with, S_without, 1e-12)
-assert eigenpy.is_approx(T_with, T_without, 1e-12)
+assert eigenpy.is_approx(S_with, S_without)
+assert eigenpy.is_approx(T_with, T_without)
 
 iterations = realqz.iterations()
 assert iterations >= 0

unittest/python/test_RealSchur.py

@@ -30,8 +30,8 @@ def verify_is_quasi_triangular(T):
 U = rs.matrixU()
 T = rs.matrixT()
 
-assert eigenpy.is_approx(A, U @ T @ U.T, 1e-10)
-assert eigenpy.is_approx(U @ U.T, np.eye(dim), 1e-10)
+assert eigenpy.is_approx(A, U @ T @ U.T)
+assert eigenpy.is_approx(U @ U.T, np.eye(dim))
 
 verify_is_quasi_triangular(T)
 
@@ -48,20 +48,20 @@ def verify_is_quasi_triangular(T):
 T2 = rs2.matrixT()
 U2 = rs2.matrixU()
 
-assert eigenpy.is_approx(T1, T2, 1e-12)
-assert eigenpy.is_approx(U1, U2, 1e-12)
+assert eigenpy.is_approx(T1, T2)
+assert eigenpy.is_approx(U1, U2)
 
 rs_no_u = eigenpy.RealSchur(A, False)
 assert rs_no_u.info() == eigenpy.ComputationInfo.Success
 T_no_u = rs_no_u.matrixT()
 
-assert eigenpy.is_approx(T, T_no_u, 1e-12)
+assert eigenpy.is_approx(T, T_no_u)
 
 rs_compute_no_u = eigenpy.RealSchur(dim)
 result_no_u = rs_compute_no_u.compute(A, False)
 assert result_no_u.info() == eigenpy.ComputationInfo.Success
 T_compute_no_u = rs_compute_no_u.matrixT()
-assert eigenpy.is_approx(T, T_compute_no_u, 1e-12)
+assert eigenpy.is_approx(T, T_compute_no_u)
 
 rs_iter = eigenpy.RealSchur(dim)
 rs_iter.setMaxIterations(40 * dim)  # m_maxIterationsPerRow * size
@@ -71,8 +71,8 @@ def verify_is_quasi_triangular(T):
 
 T_iter = rs_iter.matrixT()
 U_iter = rs_iter.matrixU()
-assert eigenpy.is_approx(T, T_iter, 1e-12)
-assert eigenpy.is_approx(U, U_iter, 1e-12)
+assert eigenpy.is_approx(T, T_iter)
+assert eigenpy.is_approx(U, U_iter)
 
 if dim > 2:
     rs_few_iter = eigenpy.RealSchur(dim)
@@ -89,8 +89,8 @@ def verify_is_quasi_triangular(T):
 T_triangular = rs_triangular.matrixT()
 U_triangular = rs_triangular.matrixU()
 
-assert eigenpy.is_approx(T_triangular, A_triangular, 1e-10)
-assert eigenpy.is_approx(U_triangular, np.eye(dim), 1e-10)
+assert eigenpy.is_approx(T_triangular, A_triangular)
+assert eigenpy.is_approx(U_triangular, np.eye(dim))
 
 hess = eigenpy.HessenbergDecomposition(A)
 H = hess.matrixH()
@@ -103,7 +103,7 @@ def verify_is_quasi_triangular(T):
 T_from_hess = rs_from_hess.matrixT()
 U_from_hess = rs_from_hess.matrixU()
 
-assert eigenpy.is_approx(A, U_from_hess @ T_from_hess @ U_from_hess.T, 1e-10)
+assert eigenpy.is_approx(A, U_from_hess @ T_from_hess @ U_from_hess.T)
 
 rs1_id = eigenpy.RealSchur(dim)
 rs2_id = eigenpy.RealSchur(dim)

unittest/python/test_Tridiagonalization.py

@@ -12,15 +12,15 @@
 Q = tri.matrixQ()
 T = tri.matrixT()
 
-assert eigenpy.is_approx(A, Q @ T @ Q.T, 1e-10)
-assert eigenpy.is_approx(Q @ Q.T, np.eye(dim), 1e-10)
+assert eigenpy.is_approx(A, Q @ T @ Q.T)
+assert eigenpy.is_approx(Q @ Q.T, np.eye(dim))
 
 for i in range(dim):
     for j in range(dim):
         if abs(i - j) > 1:
             assert abs(T[i, j]) < 1e-12
 
-assert eigenpy.is_approx(T, T.T, 1e-12)
+assert eigenpy.is_approx(T, T.T)
 
 diag = tri.diagonal()
 sub_diag = tri.subDiagonal()
@@ -43,8 +43,8 @@
 Q2 = tri2.matrixQ()
 T2 = tri2.matrixT()
 
-assert eigenpy.is_approx(Q1, Q2, 1e-12)
-assert eigenpy.is_approx(T1, T2, 1e-12)
+assert eigenpy.is_approx(Q1, Q2)
+assert eigenpy.is_approx(T1, T2)
 
 h_coeffs = tri.householderCoefficients()
 packed = tri.packedMatrix()
@@ -66,7 +66,7 @@
 Q_diag = tri_diag.matrixQ()
 T_diag = tri_diag.matrixT()
 
-assert eigenpy.is_approx(A_diag, Q_diag @ T_diag @ Q_diag.T, 1e-10)
+assert eigenpy.is_approx(A_diag, Q_diag @ T_diag @ Q_diag.T)
 for i in range(dim):
     for j in range(dim):
         if i != j:
@@ -84,7 +84,7 @@
 Q_tridiag = tri_tridiag.matrixQ()
 T_tridiag = tri_tridiag.matrixT()
 
-assert eigenpy.is_approx(A_tridiag, Q_tridiag @ T_tridiag @ Q_tridiag.T, 1e-10)
+assert eigenpy.is_approx(A_tridiag, Q_tridiag @ T_tridiag @ Q_tridiag.T)
 
 
 tri1_id = eigenpy.Tridiagonalization(dim)