From 1aa89fb85548dc425d54d2cbe7f28915c29db13a Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Fri, 1 Apr 2016 10:27:59 -0700 Subject: Add matrix condition estimator module that implements the Higham/Hager algorithm from http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf used in LPACK. Add rcond() methods to FullPivLU and PartialPivLU. --- test/lu.cpp | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) (limited to 'test') diff --git a/test/lu.cpp b/test/lu.cpp index f14435114..31991520d 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -11,6 +11,11 @@ #include using namespace std; +template +typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) { + return m.cwiseAbs().colwise().sum().maxCoeff(); +} + template void lu_non_invertible() { typedef typename MatrixType::Index Index; @@ -143,7 +148,13 @@ template void lu_invertible() m3 = MatrixType::Random(size,size); m2 = lu.solve(m3); VERIFY_IS_APPROX(m3, m1*m2); - VERIFY_IS_APPROX(m2, lu.inverse()*m3); + MatrixType m1_inverse = lu.inverse(); + VERIFY_IS_APPROX(m2, m1_inverse*m3); + + // Test condition number estimation. + RealScalar rcond = RealScalar(1) / matrix_l1_norm(m1) / matrix_l1_norm(m1_inverse); + // Verify that the estimate is within a factor of 10 of the truth. + VERIFY(lu.rcond() > rcond / 10 && lu.rcond() < rcond * 10); // test solve with transposed lu.template _solve_impl_transposed(m3, m2); @@ -170,6 +181,7 @@ template void lu_partial_piv() PartialPivLU.h */ typedef typename MatrixType::Index Index; + typedef typename NumTraits::Real RealScalar; Index size = internal::random(1,4); MatrixType m1(size, size), m2(size, size), m3(size, size); @@ -181,7 +193,13 @@ template void lu_partial_piv() m3 = MatrixType::Random(size,size); m2 = plu.solve(m3); VERIFY_IS_APPROX(m3, m1*m2); - VERIFY_IS_APPROX(m2, plu.inverse()*m3); + MatrixType m1_inverse = plu.inverse(); + VERIFY_IS_APPROX(m2, m1_inverse*m3); + + // Test condition number estimation. + RealScalar rcond = RealScalar(1) / matrix_l1_norm(m1) / matrix_l1_norm(m1_inverse); + // Verify that the estimate is within a factor of 10 of the truth. + VERIFY(plu.rcond() > rcond / 10 && plu.rcond() < rcond * 10); // test solve with transposed plu.template _solve_impl_transposed(m3, m2); -- cgit v1.2.3 From 91414e0042779b1b9d312d9255f389e67aa38106 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Fri, 1 Apr 2016 11:58:17 -0700 Subject: Fix comments in ConditionEstimator and minor cleanup. --- Eigen/src/Core/ConditionEstimator.h | 119 +++++++++++++++++++----------------- test/lu.cpp | 4 +- 2 files changed, 65 insertions(+), 58 deletions(-) (limited to 'test') diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h index ab6f59319..68e4535aa 100644 --- a/Eigen/src/Core/ConditionEstimator.h +++ b/Eigen/src/Core/ConditionEstimator.h @@ -14,7 +14,7 @@ namespace Eigen { namespace internal { template -struct EstimateInverseL1NormImpl {}; +struct EstimateInverseMatrixL1NormImpl {}; } // namespace internal template @@ -48,7 +48,6 @@ class ConditionEstimator { if (dec.rows() == 0) { return RealScalar(1); } - RealScalar matrix_l1_norm = matrix.cwiseAbs().colwise().sum().maxCoeff(); return rcond(MatrixL1Norm(matrix), dec); } @@ -76,42 +75,50 @@ class ConditionEstimator { if (matrix_norm == 0) { return 0; } - const RealScalar inverse_matrix_norm = EstimateInverseL1Norm(dec); + const RealScalar inverse_matrix_norm = EstimateInverseMatrixL1Norm(dec); return inverse_matrix_norm == 0 ? 0 : (1 / inverse_matrix_norm) / matrix_norm; } - /* - * Fast algorithm for computing a lower bound estimate on the L1 norm of - * the inverse of the matrix using at most 10 calls to the solve method on its - * decomposition. This is an implementation of Algorithm 4.1 in - * http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf - * The most common usage of this algorithm is in estimating the condition - * number ||A||_1 * ||A^{-1}||_1 of a matrix A. While ||A||_1 can be computed - * directly in O(dims^2) operations (see MatrixL1Norm() below), while - * there is no cheap closed-form expression for ||A^{-1}||_1. - * Given a decompostion of A, this algorithm estimates ||A^{-1}|| in O(dims^2) - * operations. This is done by providing operators that use the decomposition - * to solve systems of the form A x = b or A^* z = c by back-substitution, - * each costing O(dims^2) operations. Since at most 10 calls are performed, - * the total cost is O(dims^2), as opposed to O(dims^3) if the inverse matrix - * B^{-1} was formed explicitly. - */ - static RealScalar EstimateInverseL1Norm(const Decomposition& dec) { + /** + * \returns an estimate of ||inv(matrix)||_1 given a decomposition of + * matrix that implements .solve() and .adjoint().solve() methods. + * + * The method implements Algorithms 4.1 and 5.1 from + * http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf + * which also forms the basis for the condition number estimators in + * LAPACK. Since at most 10 calls to the solve method of dec are + * performed, the total cost is O(dims^2), as opposed to O(dims^3) + * needed to compute the inverse matrix explicitly. + * + * The most common usage is in estimating the condition number + * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be + * computed directly in O(n^2) operations. + */ + static RealScalar EstimateInverseMatrixL1Norm(const Decomposition& dec) { eigen_assert(dec.rows() == dec.cols()); - const int n = dec.rows(); - if (n == 0) { + if (dec.rows() == 0) { return 0; } - return internal::EstimateInverseL1NormImpl< + return internal::EstimateInverseMatrixL1NormImpl< Decomposition, NumTraits::IsComplex>::compute(dec); } + + /** + * \returns the induced matrix l1-norm + * ||matrix||_1 = sup ||matrix * v||_1 / ||v||_1, which is equal to + * the greatest absolute column sum. + */ + inline static Scalar MatrixL1Norm(const MatrixType& matrix) { + return matrix.cwiseAbs().colwise().sum().maxCoeff(); + } }; namespace internal { + // Partial specialization for real matrices. template -struct EstimateInverseL1NormImpl { +struct EstimateInverseMatrixL1NormImpl { typedef typename Decomposition::MatrixType MatrixType; typedef typename internal::traits::Scalar Scalar; typedef typename internal::plain_col_type::type Vector; @@ -130,8 +137,9 @@ struct EstimateInverseL1NormImpl { if (n == 1) { return lower_bound; } - // lower_bound is a lower bound on ||inv(A)||_1 = sup_v ||inv(A) v||_1 / - // ||v||_1 and is the objective maximized by the ("super-") gradient ascent + // lower_bound is a lower bound on + // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 + // and is the objective maximized by the ("super-") gradient ascent // algorithm. // Basic idea: We know that the optimum is achieved at one of the simplices // v = e_i, so in each iteration we follow a super-gradient to move towards @@ -143,8 +151,8 @@ struct EstimateInverseL1NormImpl { int v_max_abs_index = -1; int old_v_max_abs_index = v_max_abs_index; for (int k = 0; k < 4; ++k) { - // argmax |inv(A)^T * sign_vector| - v = dec.transpose().solve(sign_vector); + // argmax |inv(matrix)^T * sign_vector| + v = dec.adjoint().solve(sign_vector); v.cwiseAbs().maxCoeff(&v_max_abs_index); if (v_max_abs_index == old_v_max_abs_index) { // Break if the solution stagnated. @@ -153,7 +161,7 @@ struct EstimateInverseL1NormImpl { // Move to the new simplex e_j, where j = v_max_abs_index. v.setZero(); v[v_max_abs_index] = 1; - v = dec.solve(v); // v = inv(A) * e_j. + v = dec.solve(v); // v = inv(matrix) * e_j. lower_bound = VectorL1Norm(v); if (lower_bound <= old_lower_bound) { // Break if the gradient step did not increase the lower_bound. @@ -168,17 +176,16 @@ struct EstimateInverseL1NormImpl { old_v_max_abs_index = v_max_abs_index; old_lower_bound = lower_bound; } - // The following calculates an independent estimate of ||A||_1 by - // multiplying - // A by a vector with entries of slowly increasing magnitude and alternating - // sign: v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. This - // improvement - // to Hager's algorithm above is due to Higham. It was added to make the - // algorithm more robust in certain corner cases where large elements in - // the matrix might otherwise escape detection due to exact cancellation - // (especially when op and op_adjoint correspond to a sequence of - // backsubstitutions and permutations), which could cause Hager's algorithm - // to vastly underestimate ||A||_1. + // The following calculates an independent estimate of ||matrix||_1 by + // multiplying matrix by a vector with entries of slowly increasing + // magnitude and alternating sign: + // v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. + // This improvement to Hager's algorithm above is due to Higham. It was + // added to make the algorithm more robust in certain corner cases where + // large elements in the matrix might otherwise escape detection due to + // exact cancellation (especially when op and op_adjoint correspond to a + // sequence of backsubstitutions and permutations), which could cause + // Hager's algorithm to vastly underestimate ||matrix||_1. Scalar alternating_sign = 1; for (int i = 0; i < n; ++i) { v[i] = alternating_sign * static_cast(1) + @@ -194,7 +201,7 @@ struct EstimateInverseL1NormImpl { // Partial specialization for complex matrices. template -struct EstimateInverseL1NormImpl { +struct EstimateInverseMatrixL1NormImpl { typedef typename Decomposition::MatrixType MatrixType; typedef typename internal::traits::Scalar Scalar; typedef typename NumTraits::Real RealScalar; @@ -216,8 +223,9 @@ struct EstimateInverseL1NormImpl { if (n == 1) { return lower_bound; } - // lower_bound is a lower bound on ||inv(A)||_1 = sup_v ||inv(A) v||_1 / - // ||v||_1 and is the objective maximized by the ("super-") gradient ascent + // lower_bound is a lower bound on + // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 + // and is the objective maximized by the ("super-") gradient ascent // algorithm. // Basic idea: We know that the optimum is achieved at one of the simplices // v = e_i, so in each iteration we follow a super-gradient to move towards @@ -226,7 +234,7 @@ struct EstimateInverseL1NormImpl { int v_max_abs_index = -1; int old_v_max_abs_index = v_max_abs_index; for (int k = 0; k < 4; ++k) { - // argmax |inv(A)^* * sign_vector| + // argmax |inv(matrix)^* * sign_vector| RealVector abs_v = v.cwiseAbs(); const Vector psi = (abs_v.array() == 0).select(v.cwiseQuotient(abs_v), ones); @@ -240,7 +248,7 @@ struct EstimateInverseL1NormImpl { // Move to the new simplex e_j, where j = v_max_abs_index. v.setZero(); v[v_max_abs_index] = 1; - v = dec.solve(v); // v = inv(A) * e_j. + v = dec.solve(v); // v = inv(matrix) * e_j. lower_bound = VectorL1Norm(v); if (lower_bound <= old_lower_bound) { // Break if the gradient step did not increase the lower_bound. @@ -249,17 +257,16 @@ struct EstimateInverseL1NormImpl { old_v_max_abs_index = v_max_abs_index; old_lower_bound = lower_bound; } - // The following calculates an independent estimate of ||A||_1 by - // multiplying - // A by a vector with entries of slowly increasing magnitude and alternating - // sign: v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. This - // improvement - // to Hager's algorithm above is due to Higham. It was added to make the - // algorithm more robust in certain corner cases where large elements in - // the matrix might otherwise escape detection due to exact cancellation - // (especially when op and op_adjoint correspond to a sequence of - // backsubstitutions and permutations), which could cause Hager's algorithm - // to vastly underestimate ||A||_1. + // The following calculates an independent estimate of ||matrix||_1 by + // multiplying matrix by a vector with entries of slowly increasing + // magnitude and alternating sign: + // v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. + // This improvement to Hager's algorithm above is due to Higham. It was + // added to make the algorithm more robust in certain corner cases where + // large elements in the matrix might otherwise escape detection due to + // exact cancellation (especially when op and op_adjoint correspond to a + // sequence of backsubstitutions and permutations), which could cause + // Hager's algorithm to vastly underestimate ||matrix||_1. RealScalar alternating_sign = 1; for (int i = 0; i < n; ++i) { v[i] = alternating_sign * static_cast(1) + diff --git a/test/lu.cpp b/test/lu.cpp index 31991520d..9e8059f58 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -152,7 +152,7 @@ template void lu_invertible() VERIFY_IS_APPROX(m2, m1_inverse*m3); // Test condition number estimation. - RealScalar rcond = RealScalar(1) / matrix_l1_norm(m1) / matrix_l1_norm(m1_inverse); + RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); // Verify that the estimate is within a factor of 10 of the truth. VERIFY(lu.rcond() > rcond / 10 && lu.rcond() < rcond * 10); @@ -197,7 +197,7 @@ template void lu_partial_piv() VERIFY_IS_APPROX(m2, m1_inverse*m3); // Test condition number estimation. - RealScalar rcond = RealScalar(1) / matrix_l1_norm(m1) / matrix_l1_norm(m1_inverse); + RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); // Verify that the estimate is within a factor of 10 of the truth. VERIFY(plu.rcond() > rcond / 10 && plu.rcond() < rcond * 10); -- cgit v1.2.3 From f54137606eb6d68cbafd10d90013e254b26137ed Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Fri, 1 Apr 2016 16:19:45 -0700 Subject: Add condition estimation to Cholesky (LLT) factorization. --- Eigen/src/Cholesky/LLT.h | 53 +++++++++++++++++++++++++++++-------- Eigen/src/Core/ConditionEstimator.h | 37 +++++++++++++++++++------- test/cholesky.cpp | 47 +++++++++++++++++++++++--------- test/lu.cpp | 10 ++++--- 4 files changed, 111 insertions(+), 36 deletions(-) (limited to 'test') diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 74cf5bfe1..b55c5bebf 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -10,7 +10,7 @@ #ifndef EIGEN_LLT_H #define EIGEN_LLT_H -namespace Eigen { +namespace Eigen { namespace internal{ template struct LLT_Traits; @@ -40,7 +40,7 @@ template struct LLT_Traits; * * Example: \include LLT_example.cpp * Output: \verbinclude LLT_example.out - * + * * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT */ /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) @@ -135,6 +135,16 @@ template class LLT template LLT& compute(const EigenBase& matrix); + /** \returns an estimate of the reciprocal condition number of the matrix of + * which *this is the Cholesky decomposition. + */ + RealScalar rcond() const + { + eigen_assert(m_isInitialized && "LLT is not initialized."); + eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative"); + return ConditionEstimator, true >::rcond(m_l1_norm, *this); + } + /** \returns the LLT decomposition matrix * * TODO: document the storage layout @@ -164,7 +174,7 @@ template class LLT template LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); - + #ifndef EIGEN_PARSED_BY_DOXYGEN template EIGEN_DEVICE_FUNC @@ -172,17 +182,18 @@ template class LLT #endif protected: - + static void check_template_parameters() { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); } - + /** \internal * Used to compute and store L * The strict upper part is not used and even not initialized. */ MatrixType m_matrix; + RealScalar m_l1_norm; bool m_isInitialized; ComputationInfo m_info; }; @@ -268,7 +279,7 @@ template struct llt_inplace static Index unblocked(MatrixType& mat) { using std::sqrt; - + eigen_assert(mat.rows()==mat.cols()); const Index size = mat.rows(); for(Index k = 0; k < size; ++k) @@ -328,7 +339,7 @@ template struct llt_inplace return Eigen::internal::llt_rank_update_lower(mat, vec, sigma); } }; - + template struct llt_inplace { typedef typename NumTraits::Real RealScalar; @@ -387,12 +398,32 @@ template LLT& LLT::compute(const EigenBase& a) { check_template_parameters(); - + eigen_assert(a.rows()==a.cols()); const Index size = a.rows(); m_matrix.resize(size, size); m_matrix = a.derived(); + // Compute matrix L1 norm = max abs column sum. + m_l1_norm = RealScalar(0); + if (_UpLo == Lower) { + for (int col = 0; col < size; ++col) { + const RealScalar abs_col_sum = m_matrix.col(col).tail(size - col).cwiseAbs().sum() + + m_matrix.row(col).tail(col).cwiseAbs().sum(); + if (abs_col_sum > m_l1_norm) { + m_l1_norm = abs_col_sum; + } + } + } else { + for (int col = 0; col < a.cols(); ++col) { + const RealScalar abs_col_sum = m_matrix.col(col).tail(col).cwiseAbs().sum() + + m_matrix.row(col).tail(size - col).cwiseAbs().sum(); + if (abs_col_sum > m_l1_norm) { + m_l1_norm = abs_col_sum; + } + } + } + m_isInitialized = true; bool ok = Traits::inplace_decomposition(m_matrix); m_info = ok ? Success : NumericalIssue; @@ -419,7 +450,7 @@ LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, c return *this; } - + #ifndef EIGEN_PARSED_BY_DOXYGEN template template @@ -431,7 +462,7 @@ void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const #endif /** \internal use x = llt_object.solve(x); - * + * * This is the \em in-place version of solve(). * * \param bAndX represents both the right-hand side matrix b and result x. @@ -483,7 +514,7 @@ SelfAdjointView::llt() const return LLT(m_matrix); } #endif // __CUDACC__ - + } // end namespace Eigen #endif // EIGEN_LLT_H diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h index b65306d56..dca3da417 100644 --- a/Eigen/src/Core/ConditionEstimator.h +++ b/Eigen/src/Core/ConditionEstimator.h @@ -13,11 +13,11 @@ namespace Eigen { namespace internal { -template +template struct EstimateInverseMatrixL1NormImpl {}; } // namespace internal -template +template class ConditionEstimator { public: typedef typename Decomposition::MatrixType MatrixType; @@ -101,7 +101,8 @@ class ConditionEstimator { return 0; } return internal::EstimateInverseMatrixL1NormImpl< - Decomposition, NumTraits::IsComplex>::compute(dec); + Decomposition, IsSelfAdjoint, + NumTraits::IsComplex != 0>::compute(dec); } /** @@ -116,9 +117,27 @@ class ConditionEstimator { namespace internal { +template +struct solve_helper { + static inline Vector solve_adjoint(const Decomposition& dec, + const Vector& v) { + return dec.adjoint().solve(v); + } +}; + +// Partial specialization for self_adjoint matrices. +template +struct solve_helper { + static inline Vector solve_adjoint(const Decomposition& dec, + const Vector& v) { + return dec.solve(v); + } +}; + + // Partial specialization for real matrices. -template -struct EstimateInverseMatrixL1NormImpl { +template +struct EstimateInverseMatrixL1NormImpl { typedef typename Decomposition::MatrixType MatrixType; typedef typename internal::traits::Scalar Scalar; typedef typename internal::plain_col_type::type Vector; @@ -152,7 +171,7 @@ struct EstimateInverseMatrixL1NormImpl { int old_v_max_abs_index = v_max_abs_index; for (int k = 0; k < 4; ++k) { // argmax |inv(matrix)^T * sign_vector| - v = dec.adjoint().solve(sign_vector); + v = solve_helper::solve_adjoint(dec, sign_vector); v.cwiseAbs().maxCoeff(&v_max_abs_index); if (v_max_abs_index == old_v_max_abs_index) { // Break if the solution stagnated. @@ -200,8 +219,8 @@ struct EstimateInverseMatrixL1NormImpl { }; // Partial specialization for complex matrices. -template -struct EstimateInverseMatrixL1NormImpl { +template +struct EstimateInverseMatrixL1NormImpl { typedef typename Decomposition::MatrixType MatrixType; typedef typename internal::traits::Scalar Scalar; typedef typename NumTraits::Real RealScalar; @@ -238,7 +257,7 @@ struct EstimateInverseMatrixL1NormImpl { RealVector abs_v = v.cwiseAbs(); const Vector psi = (abs_v.array() == 0).select(v.cwiseQuotient(abs_v), ones); - v = dec.adjoint().solve(psi); + v = solve_helper::solve_adjoint(dec, psi); const RealVector z = v.real(); z.cwiseAbs().maxCoeff(&v_max_abs_index); if (v_max_abs_index == old_v_max_abs_index) { diff --git a/test/cholesky.cpp b/test/cholesky.cpp index d652af5bf..8a21cdbd5 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -17,6 +17,12 @@ #include #include +template +typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) { + MatrixType symm = m.template selfadjointView(); + return symm.cwiseAbs().colwise().sum().maxCoeff(); +} + template class CholType> void test_chol_update(const MatrixType& symm) { typedef typename MatrixType::Scalar Scalar; @@ -77,7 +83,7 @@ template void cholesky(const MatrixType& m) { SquareMatrixType symmUp = symm.template triangularView(); SquareMatrixType symmLo = symm.template triangularView(); - + LLT chollo(symmLo); VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); vecX = chollo.solve(vecB); @@ -85,6 +91,14 @@ template void cholesky(const MatrixType& m) matX = chollo.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); + // Verify that the estimated condition number is within a factor of 10 of the + // truth. + const MatrixType symmLo_inverse = chollo.solve(MatrixType::Identity(rows,cols)); + RealScalar rcond = (RealScalar(1) / matrix_l1_norm(symmLo)) / + matrix_l1_norm(symmLo_inverse); + RealScalar rcond_est = chollo.rcond(); + VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + // test the upper mode LLT cholup(symmUp); VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix()); @@ -93,6 +107,15 @@ template void cholesky(const MatrixType& m) matX = cholup.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); + // Verify that the estimated condition number is within a factor of 10 of the + // truth. + const MatrixType symmUp_inverse = cholup.solve(MatrixType::Identity(rows,cols)); + rcond = (RealScalar(1) / matrix_l1_norm(symmUp)) / + matrix_l1_norm(symmUp_inverse); + rcond_est = cholup.rcond(); + VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + + MatrixType neg = -symmLo; chollo.compute(neg); VERIFY(chollo.info()==NumericalIssue); @@ -101,7 +124,7 @@ template void cholesky(const MatrixType& m) VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL())); VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU())); VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL())); - + // test some special use cases of SelfCwiseBinaryOp: MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); m2 = m1; @@ -167,7 +190,7 @@ template void cholesky(const MatrixType& m) // restore if(sign == -1) symm = -symm; - + // check matrices coming from linear constraints with Lagrange multipliers if(rows>=3) { @@ -183,7 +206,7 @@ template void cholesky(const MatrixType& m) vecX = ldltlo.solve(vecB); VERIFY_IS_APPROX(A * vecX, vecB); } - + // check non-full rank matrices if(rows>=3) { @@ -199,7 +222,7 @@ template void cholesky(const MatrixType& m) vecX = ldltlo.solve(vecB); VERIFY_IS_APPROX(A * vecX, vecB); } - + // check matrices with a wide spectrum if(rows>=3) { @@ -225,7 +248,7 @@ template void cholesky(const MatrixType& m) { RealScalar large_tol = std::sqrt(test_precision()); VERIFY((A * vecX).isApprox(vecB, large_tol)); - + ++g_test_level; VERIFY_IS_APPROX(A * vecX,vecB); --g_test_level; @@ -314,14 +337,14 @@ template void cholesky_bug241(const MatrixType& m) } // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. -// This test checks that LDLT reports correctly that matrix is indefinite. +// This test checks that LDLT reports correctly that matrix is indefinite. // See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736 template void cholesky_definiteness(const MatrixType& m) { eigen_assert(m.rows() == 2 && m.cols() == 2); MatrixType mat; LDLT ldlt(2); - + { mat << 1, 0, 0, -1; ldlt.compute(mat); @@ -384,11 +407,11 @@ void test_cholesky() CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) ); CALL_SUBTEST_4( cholesky(Matrix3f()) ); CALL_SUBTEST_5( cholesky(Matrix4d()) ); - - s = internal::random(1,EIGEN_TEST_MAX_SIZE); + + s = internal::random(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) ); TEST_SET_BUT_UNUSED_VARIABLE(s) - + s = internal::random(1,EIGEN_TEST_MAX_SIZE/2); CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); TEST_SET_BUT_UNUSED_VARIABLE(s) @@ -402,6 +425,6 @@ void test_cholesky() // Test problem size constructors CALL_SUBTEST_9( LLT(10) ); CALL_SUBTEST_9( LDLT(10) ); - + TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) } diff --git a/test/lu.cpp b/test/lu.cpp index 9e8059f58..53b3fcee4 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -151,10 +151,11 @@ template void lu_invertible() MatrixType m1_inverse = lu.inverse(); VERIFY_IS_APPROX(m2, m1_inverse*m3); - // Test condition number estimation. + // Verify that the estimated condition number is within a factor of 10 of the + // truth. RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); - // Verify that the estimate is within a factor of 10 of the truth. - VERIFY(lu.rcond() > rcond / 10 && lu.rcond() < rcond * 10); + const RealScalar rcond_est = lu.rcond(); + VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); // test solve with transposed lu.template _solve_impl_transposed(m3, m2); @@ -199,7 +200,8 @@ template void lu_partial_piv() // Test condition number estimation. RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); // Verify that the estimate is within a factor of 10 of the truth. - VERIFY(plu.rcond() > rcond / 10 && plu.rcond() < rcond * 10); + const RealScalar rcond_est = plu.rcond(); + VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); // test solve with transposed plu.template _solve_impl_transposed(m3, m2); -- cgit v1.2.3 From 9d51f7c457671bfcbab9a1d62d416e1a83e6ad8a Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Fri, 1 Apr 2016 16:48:38 -0700 Subject: Add rcond method to LDLT. --- Eigen/src/Cholesky/LDLT.h | 54 ++++++++++++++++++++++++++++++++++++----------- test/cholesky.cpp | 17 +++++++++++++++ 2 files changed, 59 insertions(+), 12 deletions(-) (limited to 'test') diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index c3cc3746c..9753c84d8 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -13,7 +13,7 @@ #ifndef EIGEN_LDLT_H #define EIGEN_LDLT_H -namespace Eigen { +namespace Eigen { namespace internal { template struct LDLT_Traits; @@ -73,11 +73,11 @@ template class LDLT * The default constructor is useful in cases in which the user intends to * perform decompositions via LDLT::compute(const MatrixType&). */ - LDLT() - : m_matrix(), - m_transpositions(), + LDLT() + : m_matrix(), + m_transpositions(), m_sign(internal::ZeroSign), - m_isInitialized(false) + m_isInitialized(false) {} /** \brief Default Constructor with memory preallocation @@ -168,7 +168,7 @@ template class LDLT * \note_about_checking_solutions * * More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$ - * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$, + * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$, * \f$ L^* y_4 = y_3 \f$ and \f$ P x = y_4 \f$ in succession. If the matrix \f$ A \f$ is singular, then * \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the * least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function @@ -192,6 +192,15 @@ template class LDLT template LDLT& compute(const EigenBase& matrix); + /** \returns an estimate of the reciprocal condition number of the matrix of + * which *this is the LDLT decomposition. + */ + RealScalar rcond() const + { + eigen_assert(m_isInitialized && "LDLT is not initialized."); + return ConditionEstimator, true >::rcond(m_l1_norm, *this); + } + template LDLT& rankUpdate(const MatrixBase& w, const RealScalar& alpha=1); @@ -220,7 +229,7 @@ template class LDLT eigen_assert(m_isInitialized && "LDLT is not initialized."); return Success; } - + #ifndef EIGEN_PARSED_BY_DOXYGEN template EIGEN_DEVICE_FUNC @@ -228,7 +237,7 @@ template class LDLT #endif protected: - + static void check_template_parameters() { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); @@ -241,6 +250,7 @@ template class LDLT * is not stored), and the diagonal entries correspond to D. */ MatrixType m_matrix; + RealScalar m_l1_norm; TranspositionType m_transpositions; TmpMatrixType m_temporary; internal::SignMatrix m_sign; @@ -314,7 +324,7 @@ template<> struct ldlt_inplace if(rs>0) A21.noalias() -= A20 * temp.head(k); } - + // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot // was smaller than the cutoff value. However, since LDLT is not rank-revealing // we should only make sure that we do not introduce INF or NaN values. @@ -433,12 +443,32 @@ template LDLT& LDLT::compute(const EigenBase& a) { check_template_parameters(); - + eigen_assert(a.rows()==a.cols()); const Index size = a.rows(); m_matrix = a.derived(); + // Compute matrix L1 norm = max abs column sum. + m_l1_norm = RealScalar(0); + if (_UpLo == Lower) { + for (int col = 0; col < size; ++col) { + const RealScalar abs_col_sum = m_matrix.col(col).tail(size - col).cwiseAbs().sum() + + m_matrix.row(col).tail(col).cwiseAbs().sum(); + if (abs_col_sum > m_l1_norm) { + m_l1_norm = abs_col_sum; + } + } + } else { + for (int col = 0; col < a.cols(); ++col) { + const RealScalar abs_col_sum = m_matrix.col(col).tail(col).cwiseAbs().sum() + + m_matrix.row(col).tail(size - col).cwiseAbs().sum(); + if (abs_col_sum > m_l1_norm) { + m_l1_norm = abs_col_sum; + } + } + } + m_transpositions.resize(size); m_isInitialized = false; m_temporary.resize(size); @@ -466,7 +496,7 @@ LDLT& LDLT::rankUpdate(const MatrixBase::_solve_impl(const RhsType &rhs, DstType &dst) cons // diagonal element is not well justified and leads to numerical issues in some cases. // Moreover, Lapack's xSYTRS routines use 0 for the tolerance. RealScalar tolerance = RealScalar(1) / NumTraits::highest(); - + for (Index i = 0; i < vecD.size(); ++i) { if(abs(vecD(i)) > tolerance) diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 8a21cdbd5..148a0b388 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -160,6 +160,15 @@ template void cholesky(const MatrixType& m) matX = ldltlo.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); + // Verify that the estimated condition number is within a factor of 10 of the + // truth. + const MatrixType symmLo_inverse = ldltlo.solve(MatrixType::Identity(rows,cols)); + RealScalar rcond = (RealScalar(1) / matrix_l1_norm(symmLo)) / + matrix_l1_norm(symmLo_inverse); + RealScalar rcond_est = ldltlo.rcond(); + VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + + LDLT ldltup(symmUp); VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix()); vecX = ldltup.solve(vecB); @@ -167,6 +176,14 @@ template void cholesky(const MatrixType& m) matX = ldltup.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); + // Verify that the estimated condition number is within a factor of 10 of the + // truth. + const MatrixType symmUp_inverse = ldltup.solve(MatrixType::Identity(rows,cols)); + rcond = (RealScalar(1) / matrix_l1_norm(symmUp)) / + matrix_l1_norm(symmUp_inverse); + rcond_est = ldltup.rcond(); + VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU())); VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL())); VERIFY_IS_APPROX(MatrixType(ldltup.matrixL().transpose().conjugate()), MatrixType(ldltup.matrixU())); -- cgit v1.2.3 From 86e0ed81f8db5a0c9562b62a67a9ba60ec58dec0 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 4 Apr 2016 14:20:01 -0700 Subject: Addresses comments on Eigen pull request PR-174. * Get rid of code-duplication for real vs. complex matrices. * Fix flipped arguments to select. * Make the condition estimation functions free functions. * Use Vector::Unit() to generate canonical unit vectors. * Misc. cleanup. --- Eigen/src/Cholesky/LDLT.h | 12 +- Eigen/src/Cholesky/LLT.h | 12 +- Eigen/src/Core/ConditionEstimator.h | 405 ++++++++++++++---------------------- Eigen/src/LU/FullPivLU.h | 2 +- Eigen/src/LU/PartialPivLU.h | 2 +- test/cholesky.cpp | 8 +- test/lu.cpp | 7 +- 7 files changed, 179 insertions(+), 269 deletions(-) (limited to 'test') diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 9753c84d8..80f18977c 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -198,7 +198,7 @@ template class LDLT RealScalar rcond() const { eigen_assert(m_isInitialized && "LDLT is not initialized."); - return ConditionEstimator, true >::rcond(m_l1_norm, *this); + return ReciprocalConditionNumberEstimate(m_l1_norm, *this); } template @@ -216,6 +216,12 @@ template class LDLT MatrixType reconstructedMatrix() const; + /** \returns the decomposition itself to allow generic code to do + * ldlt.transpose().solve(rhs). + */ + const LDLT& transpose() const { return *this; }; + const LDLT& adjoint() const { return *this; }; + inline Index rows() const { return m_matrix.rows(); } inline Index cols() const { return m_matrix.cols(); } @@ -454,14 +460,14 @@ LDLT& LDLT::compute(const EigenBase m_l1_norm) { m_l1_norm = abs_col_sum; } } } else { for (int col = 0; col < a.cols(); ++col) { - const RealScalar abs_col_sum = m_matrix.col(col).tail(col).cwiseAbs().sum() + + const RealScalar abs_col_sum = m_matrix.col(col).head(col).cwiseAbs().sum() + m_matrix.row(col).tail(size - col).cwiseAbs().sum(); if (abs_col_sum > m_l1_norm) { m_l1_norm = abs_col_sum; diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index b55c5bebf..94da1d52d 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -142,7 +142,7 @@ template class LLT { eigen_assert(m_isInitialized && "LLT is not initialized."); eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative"); - return ConditionEstimator, true >::rcond(m_l1_norm, *this); + return ReciprocalConditionNumberEstimate(m_l1_norm, *this); } /** \returns the LLT decomposition matrix @@ -169,6 +169,12 @@ template class LLT return m_info; } + /** \returns the decomposition itself to allow generic code to do + * llt.transpose().solve(rhs). + */ + const LLT& transpose() const { return *this; }; + const LLT& adjoint() const { return *this; }; + inline Index rows() const { return m_matrix.rows(); } inline Index cols() const { return m_matrix.cols(); } @@ -409,14 +415,14 @@ LLT& LLT::compute(const EigenBase if (_UpLo == Lower) { for (int col = 0; col < size; ++col) { const RealScalar abs_col_sum = m_matrix.col(col).tail(size - col).cwiseAbs().sum() + - m_matrix.row(col).tail(col).cwiseAbs().sum(); + m_matrix.row(col).head(col).cwiseAbs().sum(); if (abs_col_sum > m_l1_norm) { m_l1_norm = abs_col_sum; } } } else { for (int col = 0; col < a.cols(); ++col) { - const RealScalar abs_col_sum = m_matrix.col(col).tail(col).cwiseAbs().sum() + + const RealScalar abs_col_sum = m_matrix.col(col).head(col).cwiseAbs().sum() + m_matrix.row(col).tail(size - col).cwiseAbs().sum(); if (abs_col_sum > m_l1_norm) { m_l1_norm = abs_col_sum; diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h index dca3da417..12b4ae648 100644 --- a/Eigen/src/Core/ConditionEstimator.h +++ b/Eigen/src/Core/ConditionEstimator.h @@ -13,45 +13,35 @@ namespace Eigen { namespace internal { -template -struct EstimateInverseMatrixL1NormImpl {}; -} // namespace internal - -template -class ConditionEstimator { - public: - typedef typename Decomposition::MatrixType MatrixType; - typedef typename internal::traits::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef typename internal::plain_col_type::type Vector; +template +inline typename MatrixType::RealScalar MatrixL1Norm(const MatrixType& matrix) { + return matrix.cwiseAbs().colwise().sum().maxCoeff(); +} + +template +inline typename Vector::RealScalar VectorL1Norm(const Vector& v) { + return v.template lpNorm<1>(); +} + +template +struct SignOrUnity { + static inline Vector run(const Vector& v) { + const RealVector v_abs = v.cwiseAbs(); + return (v_abs.array() == 0).select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs)); + } +}; - /** \class ConditionEstimator - * \ingroup Core_Module - * - * \brief Condition number estimator. - * - * Computing a decomposition of a dense matrix takes O(n^3) operations, while - * this method estimates the condition number quickly and reliably in O(n^2) - * operations. - * - * \returns an estimate of the reciprocal condition number - * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given the matrix and - * its decomposition. Supports the following decompositions: FullPivLU, - * PartialPivLU. - * - * \sa FullPivLU, PartialPivLU. - */ - static RealScalar rcond(const MatrixType& matrix, const Decomposition& dec) { - eigen_assert(matrix.rows() == dec.rows()); - eigen_assert(matrix.cols() == dec.cols()); - eigen_assert(matrix.rows() == matrix.cols()); - if (dec.rows() == 0) { - return RealScalar(1); - } - return rcond(MatrixL1Norm(matrix), dec); +// Partial specialization to avoid elementwise division for real vectors. +template +struct SignOrUnity { + static inline Vector run(const Vector& v) { + return (v.array() < 0).select(-Vector::Ones(v.size()), Vector::Ones(v.size())); } +}; - /** \class ConditionEstimator +} // namespace internal + +/** \class ConditionEstimator * \ingroup Core_Module * * \brief Condition number estimator. @@ -61,245 +51,154 @@ class ConditionEstimator { * operations. * * \returns an estimate of the reciprocal condition number - * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and + * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given the matrix and * its decomposition. Supports the following decompositions: FullPivLU, * PartialPivLU. * * \sa FullPivLU, PartialPivLU. */ - static RealScalar rcond(RealScalar matrix_norm, const Decomposition& dec) { - eigen_assert(dec.rows() == dec.cols()); - if (dec.rows() == 0) { - return 1; - } - if (matrix_norm == 0) { - return 0; - } - const RealScalar inverse_matrix_norm = EstimateInverseMatrixL1Norm(dec); - return inverse_matrix_norm == 0 ? 0 - : (1 / inverse_matrix_norm) / matrix_norm; +template +typename Decomposition::RealScalar ReciprocalConditionNumberEstimate( + const typename Decomposition::MatrixType& matrix, + const Decomposition& dec) { + eigen_assert(matrix.rows() == dec.rows()); + eigen_assert(matrix.cols() == dec.cols()); + eigen_assert(matrix.rows() == matrix.cols()); + if (dec.rows() == 0) { + return Decomposition::RealScalar(1); } - - /** - * \returns an estimate of ||inv(matrix)||_1 given a decomposition of - * matrix that implements .solve() and .adjoint().solve() methods. - * - * The method implements Algorithms 4.1 and 5.1 from - * http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf - * which also forms the basis for the condition number estimators in - * LAPACK. Since at most 10 calls to the solve method of dec are - * performed, the total cost is O(dims^2), as opposed to O(dims^3) - * needed to compute the inverse matrix explicitly. - * - * The most common usage is in estimating the condition number - * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be - * computed directly in O(n^2) operations. - */ - static RealScalar EstimateInverseMatrixL1Norm(const Decomposition& dec) { - eigen_assert(dec.rows() == dec.cols()); - if (dec.rows() == 0) { - return 0; - } - return internal::EstimateInverseMatrixL1NormImpl< - Decomposition, IsSelfAdjoint, - NumTraits::IsComplex != 0>::compute(dec); + return ReciprocalConditionNumberEstimate(MatrixL1Norm(matrix), dec); +} + +/** \class ConditionEstimator + * \ingroup Core_Module + * + * \brief Condition number estimator. + * + * Computing a decomposition of a dense matrix takes O(n^3) operations, while + * this method estimates the condition number quickly and reliably in O(n^2) + * operations. + * + * \returns an estimate of the reciprocal condition number + * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and + * its decomposition. Supports the following decompositions: FullPivLU, + * PartialPivLU. + * + * \sa FullPivLU, PartialPivLU. + */ +template +typename Decomposition::RealScalar ReciprocalConditionNumberEstimate( + typename Decomposition::RealScalar matrix_norm, const Decomposition& dec) { + eigen_assert(dec.rows() == dec.cols()); + if (dec.rows() == 0) { + return 1; } - - /** - * \returns the induced matrix l1-norm - * ||matrix||_1 = sup ||matrix * v||_1 / ||v||_1, which is equal to - * the greatest absolute column sum. - */ - static inline Scalar MatrixL1Norm(const MatrixType& matrix) { - return matrix.cwiseAbs().colwise().sum().maxCoeff(); + if (matrix_norm == 0) { + return 0; } -}; - -namespace internal { - -template -struct solve_helper { - static inline Vector solve_adjoint(const Decomposition& dec, - const Vector& v) { - return dec.adjoint().solve(v); - } -}; - -// Partial specialization for self_adjoint matrices. -template -struct solve_helper { - static inline Vector solve_adjoint(const Decomposition& dec, - const Vector& v) { - return dec.solve(v); - } -}; - - -// Partial specialization for real matrices. -template -struct EstimateInverseMatrixL1NormImpl { + const typename Decomposition::RealScalar inverse_matrix_norm = InverseMatrixL1NormEstimate(dec); + return inverse_matrix_norm == 0 ? 0 : (1 / inverse_matrix_norm) / matrix_norm; +} + +/** + * \returns an estimate of ||inv(matrix)||_1 given a decomposition of + * matrix that implements .solve() and .adjoint().solve() methods. + * + * The method implements Algorithms 4.1 and 5.1 from + * http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf + * which also forms the basis for the condition number estimators in + * LAPACK. Since at most 10 calls to the solve method of dec are + * performed, the total cost is O(dims^2), as opposed to O(dims^3) + * needed to compute the inverse matrix explicitly. + * + * The most common usage is in estimating the condition number + * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be + * computed directly in O(n^2) operations. + */ +template +typename Decomposition::RealScalar InverseMatrixL1NormEstimate( + const Decomposition& dec) { typedef typename Decomposition::MatrixType MatrixType; - typedef typename internal::traits::Scalar Scalar; + typedef typename Decomposition::Scalar Scalar; + typedef typename Decomposition::RealScalar RealScalar; typedef typename internal::plain_col_type::type Vector; + typedef typename internal::plain_col_type::type RealVector; + const bool is_complex = (NumTraits::IsComplex != 0); - // Shorthand for vector L1 norm in Eigen. - static inline Scalar VectorL1Norm(const Vector& v) { - return v.template lpNorm<1>(); + eigen_assert(dec.rows() == dec.cols()); + const int n = dec.rows(); + if (n == 0) { + return 0; } - - static inline Scalar compute(const Decomposition& dec) { - const int n = dec.rows(); - const Vector plus = Vector::Ones(n); - Vector v = plus / n; - v = dec.solve(v); - Scalar lower_bound = VectorL1Norm(v); - if (n == 1) { - return lower_bound; - } - // lower_bound is a lower bound on - // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 - // and is the objective maximized by the ("super-") gradient ascent - // algorithm. - // Basic idea: We know that the optimum is achieved at one of the simplices - // v = e_i, so in each iteration we follow a super-gradient to move towards - // the optimal one. - Scalar old_lower_bound = lower_bound; - const Vector minus = -Vector::Ones(n); - Vector sign_vector = (v.cwiseAbs().array() == 0).select(plus, minus); - Vector old_sign_vector = sign_vector; - int v_max_abs_index = -1; - int old_v_max_abs_index = v_max_abs_index; - for (int k = 0; k < 4; ++k) { - // argmax |inv(matrix)^T * sign_vector| - v = solve_helper::solve_adjoint(dec, sign_vector); - v.cwiseAbs().maxCoeff(&v_max_abs_index); - if (v_max_abs_index == old_v_max_abs_index) { - // Break if the solution stagnated. - break; - } - // Move to the new simplex e_j, where j = v_max_abs_index. - v.setZero(); - v[v_max_abs_index] = 1; - v = dec.solve(v); // v = inv(matrix) * e_j. - lower_bound = VectorL1Norm(v); - if (lower_bound <= old_lower_bound) { - // Break if the gradient step did not increase the lower_bound. - break; - } - sign_vector = (v.array() < 0).select(plus, minus); + Vector v = Vector::Ones(n) / n; + v = dec.solve(v); + + // lower_bound is a lower bound on + // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 + // and is the objective maximized by the ("super-") gradient ascent + // algorithm below. + RealScalar lower_bound = internal::VectorL1Norm(v); + if (n == 1) { + return lower_bound; + } + // Gradient ascent algorithm follows: We know that the optimum is achieved at + // one of the simplices v = e_i, so in each iteration we follow a + // super-gradient to move towards the optimal one. + RealScalar old_lower_bound = lower_bound; + Vector sign_vector(n); + Vector old_sign_vector; + int v_max_abs_index = -1; + int old_v_max_abs_index = v_max_abs_index; + for (int k = 0; k < 4; ++k) { + sign_vector = internal::SignOrUnity::run(v); + if (k > 0 && !is_complex) { if (sign_vector == old_sign_vector) { // Break if the solution stagnated. break; } - old_sign_vector = sign_vector; - old_v_max_abs_index = v_max_abs_index; - old_lower_bound = lower_bound; - } - // The following calculates an independent estimate of ||matrix||_1 by - // multiplying matrix by a vector with entries of slowly increasing - // magnitude and alternating sign: - // v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. - // This improvement to Hager's algorithm above is due to Higham. It was - // added to make the algorithm more robust in certain corner cases where - // large elements in the matrix might otherwise escape detection due to - // exact cancellation (especially when op and op_adjoint correspond to a - // sequence of backsubstitutions and permutations), which could cause - // Hager's algorithm to vastly underestimate ||matrix||_1. - Scalar alternating_sign = 1; - for (int i = 0; i < n; ++i) { - v[i] = alternating_sign * static_cast(1) + - (static_cast(i) / (static_cast(n - 1))); - alternating_sign = -alternating_sign; } - v = dec.solve(v); - const Scalar alternate_lower_bound = - (2 * VectorL1Norm(v)) / (3 * static_cast(n)); - return numext::maxi(lower_bound, alternate_lower_bound); - } -}; - -// Partial specialization for complex matrices. -template -struct EstimateInverseMatrixL1NormImpl { - typedef typename Decomposition::MatrixType MatrixType; - typedef typename internal::traits::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef typename internal::plain_col_type::type Vector; - typedef typename internal::plain_col_type::type - RealVector; - - // Shorthand for vector L1 norm in Eigen. - static inline RealScalar VectorL1Norm(const Vector& v) { - return v.template lpNorm<1>(); - } - - static inline RealScalar compute(const Decomposition& dec) { - const int n = dec.rows(); - const Vector ones = Vector::Ones(n); - Vector v = ones / n; - v = dec.solve(v); - RealScalar lower_bound = VectorL1Norm(v); - if (n == 1) { - return lower_bound; + // v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )| + v = dec.adjoint().solve(sign_vector); + v.real().cwiseAbs().maxCoeff(&v_max_abs_index); + if (v_max_abs_index == old_v_max_abs_index) { + // Break if the solution stagnated. + break; } - // lower_bound is a lower bound on - // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 - // and is the objective maximized by the ("super-") gradient ascent - // algorithm. - // Basic idea: We know that the optimum is achieved at one of the simplices - // v = e_i, so in each iteration we follow a super-gradient to move towards - // the optimal one. - RealScalar old_lower_bound = lower_bound; - int v_max_abs_index = -1; - int old_v_max_abs_index = v_max_abs_index; - for (int k = 0; k < 4; ++k) { - // argmax |inv(matrix)^* * sign_vector| - RealVector abs_v = v.cwiseAbs(); - const Vector psi = - (abs_v.array() == 0).select(v.cwiseQuotient(abs_v), ones); - v = solve_helper::solve_adjoint(dec, psi); - const RealVector z = v.real(); - z.cwiseAbs().maxCoeff(&v_max_abs_index); - if (v_max_abs_index == old_v_max_abs_index) { - // Break if the solution stagnated. - break; - } - // Move to the new simplex e_j, where j = v_max_abs_index. - v.setZero(); - v[v_max_abs_index] = 1; - v = dec.solve(v); // v = inv(matrix) * e_j. - lower_bound = VectorL1Norm(v); - if (lower_bound <= old_lower_bound) { - // Break if the gradient step did not increase the lower_bound. - break; - } - old_v_max_abs_index = v_max_abs_index; - old_lower_bound = lower_bound; + // Move to the new simplex e_j, where j = v_max_abs_index. + v = dec.solve(Vector::Unit(n, v_max_abs_index)); // v = inv(matrix) * e_j. + lower_bound = internal::VectorL1Norm(v); + if (lower_bound <= old_lower_bound) { + // Break if the gradient step did not increase the lower_bound. + break; } - // The following calculates an independent estimate of ||matrix||_1 by - // multiplying matrix by a vector with entries of slowly increasing - // magnitude and alternating sign: - // v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. - // This improvement to Hager's algorithm above is due to Higham. It was - // added to make the algorithm more robust in certain corner cases where - // large elements in the matrix might otherwise escape detection due to - // exact cancellation (especially when op and op_adjoint correspond to a - // sequence of backsubstitutions and permutations), which could cause - // Hager's algorithm to vastly underestimate ||matrix||_1. - RealScalar alternating_sign = 1; - for (int i = 0; i < n; ++i) { - v[i] = alternating_sign * static_cast(1) + - (static_cast(i) / (static_cast(n - 1))); - alternating_sign = -alternating_sign; + if (!is_complex) { + old_sign_vector = sign_vector; } - v = dec.solve(v); - const RealScalar alternate_lower_bound = - (2 * VectorL1Norm(v)) / (3 * static_cast(n)); - return numext::maxi(lower_bound, alternate_lower_bound); + old_v_max_abs_index = v_max_abs_index; + old_lower_bound = lower_bound; } -}; + // The following calculates an independent estimate of ||matrix||_1 by + // multiplying matrix by a vector with entries of slowly increasing + // magnitude and alternating sign: + // v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. + // This improvement to Hager's algorithm above is due to Higham. It was + // added to make the algorithm more robust in certain corner cases where + // large elements in the matrix might otherwise escape detection due to + // exact cancellation (especially when op and op_adjoint correspond to a + // sequence of backsubstitutions and permutations), which could cause + // Hager's algorithm to vastly underestimate ||matrix||_1. + Scalar alternating_sign = 1; + for (int i = 0; i < n; ++i) { + v[i] = alternating_sign * static_cast(1) + + (static_cast(i) / (static_cast(n - 1))); + alternating_sign = -alternating_sign; + } + v = dec.solve(v); + const RealScalar alternate_lower_bound = + (2 * internal::VectorL1Norm(v)) / (3 * static_cast(n)); + return numext::maxi(lower_bound, alternate_lower_bound); +} -} // namespace internal } // namespace Eigen #endif diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index ff0b78c35..978a54eff 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -237,7 +237,7 @@ template class FullPivLU inline RealScalar rcond() const { eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); - return ConditionEstimator >::rcond(m_l1_norm, *this); + return ReciprocalConditionNumberEstimate(m_l1_norm, *this); } /** \returns the determinant of the matrix of which diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 5d71a66d0..b22dd75fe 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -157,7 +157,7 @@ template class PartialPivLU inline RealScalar rcond() const { eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); - return ConditionEstimator >::rcond(m_l1_norm, *this); + return ReciprocalConditionNumberEstimate(m_l1_norm, *this); } /** \returns the inverse of the matrix of which *this is the LU decomposition. diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 148a0b388..b7abc230b 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -91,12 +91,12 @@ template void cholesky(const MatrixType& m) matX = chollo.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); - // Verify that the estimated condition number is within a factor of 10 of the - // truth. const MatrixType symmLo_inverse = chollo.solve(MatrixType::Identity(rows,cols)); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(symmLo)) / matrix_l1_norm(symmLo_inverse); RealScalar rcond_est = chollo.rcond(); + // Verify that the estimated condition number is within a factor of 10 of the + // truth. VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); // test the upper mode @@ -160,12 +160,12 @@ template void cholesky(const MatrixType& m) matX = ldltlo.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); - // Verify that the estimated condition number is within a factor of 10 of the - // truth. const MatrixType symmLo_inverse = ldltlo.solve(MatrixType::Identity(rows,cols)); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(symmLo)) / matrix_l1_norm(symmLo_inverse); RealScalar rcond_est = ldltlo.rcond(); + // Verify that the estimated condition number is within a factor of 10 of the + // truth. VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); diff --git a/test/lu.cpp b/test/lu.cpp index 53b3fcee4..9787f4d86 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -151,10 +151,10 @@ template void lu_invertible() MatrixType m1_inverse = lu.inverse(); VERIFY_IS_APPROX(m2, m1_inverse*m3); - // Verify that the estimated condition number is within a factor of 10 of the - // truth. RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); const RealScalar rcond_est = lu.rcond(); + // Verify that the estimated condition number is within a factor of 10 of the + // truth. VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); // test solve with transposed @@ -197,10 +197,9 @@ template void lu_partial_piv() MatrixType m1_inverse = plu.inverse(); VERIFY_IS_APPROX(m2, m1_inverse*m3); - // Test condition number estimation. RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); - // Verify that the estimate is within a factor of 10 of the truth. const RealScalar rcond_est = plu.rcond(); + // Verify that the estimate is within a factor of 10 of the truth. VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); // test solve with transposed -- cgit v1.2.3