diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h index ab6f593190..68e4535aad 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 31991520df..9e8059f583 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);