diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h index 34dae18bb..df27be3bc 100644 --- a/Eigen/src/Core/ConditionEstimator.h +++ b/Eigen/src/Core/ConditionEstimator.h @@ -40,18 +40,17 @@ struct rcond_compute_sign { * \a matrix that implements .solve() and .adjoint().solve() methods. * * This function 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. + * Higham, "Experience with a Matrix Norm Estimator", + * SIAM J. Sci. Stat. Comput., 11(4):804-809, 1990. + * with Higham's alternating-sign safety-net estimate from + * Higham and Tisseur, "A Block Algorithm for Matrix 1-Norm Estimation, + * with an Application to 1-Norm Pseudospectra", SIAM J. Matrix Anal. Appl., + * 21(4):1185-1201, 2000. * - * 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. + * The Hager/Higham gradient ascent uses at most 5 iterations of 2 solves + * each, giving a total cost of O(n^2). * - * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, and - * LLT. + * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, LLT. * * \sa FullPivLU, PartialPivLU, LDLT, LLT. */ @@ -66,7 +65,7 @@ typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomp eigen_assert(dec.rows() == dec.cols()); const Index n = dec.rows(); - if (n == 0) return 0; + if (n == 0) return RealScalar(0); // Disable Index to float conversion warning #ifdef __INTEL_COMPILER @@ -80,14 +79,12 @@ typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomp // 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. + // and is the objective maximized by the supergradient ascent algorithm below. RealScalar lower_bound = v.template lpNorm<1>(); 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. + // Gradient ascent: the optimum is achieved at a unit vector e_j. Each + // iteration follows the supergradient to find which unit vector to probe next. RealScalar old_lower_bound = lower_bound; Vector sign_vector(n); Vector old_sign_vector; @@ -96,21 +93,21 @@ typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomp for (int k = 0; k < 4; ++k) { sign_vector = internal::rcond_compute_sign::run(v); if (k > 0 && !is_complex && sign_vector == old_sign_vector) { - // Break if the solution stagnated. + // Break if the sign vector stagnated. break; } - // v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )| + // Supergradient: z = A^{-T} * sign(v), pick argmax |z_i|. 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. + // Optimality: supergradient points to the same unit vector. break; } - // 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. + // Probe the best unit vector: v = A^{-1} * e_j. + v = dec.solve(Vector::Unit(n, v_max_abs_index)); lower_bound = v.template lpNorm<1>(); if (lower_bound <= old_lower_bound) { - // Break if the gradient step did not increase the lower_bound. + // No improvement from the gradient step. break; } if (!is_complex) { @@ -119,25 +116,19 @@ typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomp 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. + // Higham's alternating-sign estimate: an independent safety-net that catches + // cases where the gradient ascent converges to a local maximum due to exact + // cancellation patterns (especially with permutations and backsubstitutions). + // v_i = (-1)^i * (1 + i/(n-1)), then estimate = 2*||A^{-1}*v||_1 / (3*n). Scalar alternating_sign(RealScalar(1)); for (Index i = 0; i < n; ++i) { - // The static_cast is needed when Scalar is a complex and RealScalar implements expression templates + // The static_cast is needed when Scalar is complex and RealScalar uses expression templates. v[i] = alternating_sign * static_cast(RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1)))); alternating_sign = -alternating_sign; } v = dec.solve(v); - const RealScalar alternate_lower_bound = (2 * v.template lpNorm<1>()) / (3 * RealScalar(n)); - return numext::maxi(lower_bound, alternate_lower_bound); + const RealScalar alt_est = (RealScalar(2) * v.template lpNorm<1>()) / (RealScalar(3) * RealScalar(n)); + return numext::maxi(lower_bound, alt_est); } /** \brief Reciprocal condition number estimator. diff --git a/test/condition_estimator.cpp b/test/condition_estimator.cpp index 94e9e1917..ce65f6ea6 100644 --- a/test/condition_estimator.cpp +++ b/test/condition_estimator.cpp @@ -29,8 +29,8 @@ void rcond_partial_piv_lu() { MatrixType m_inverse = lu.inverse(); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond_est = lu.rcond(); - // Verify the estimate is within a factor of 10 of the truth. - VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + // Verify the estimate is within a factor of 3 of the truth. + VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3); } template @@ -47,7 +47,7 @@ void rcond_full_piv_lu() { MatrixType m_inverse = lu.inverse(); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond_est = lu.rcond(); - VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3); } template @@ -65,7 +65,7 @@ void rcond_llt() { MatrixType m_inverse = llt.solve(MatrixType::Identity(size, size)); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond_est = llt.rcond(); - VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3); } template @@ -83,7 +83,7 @@ void rcond_ldlt() { MatrixType m_inverse = ldlt.solve(MatrixType::Identity(size, size)); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond_est = ldlt.rcond(); - VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3); } template @@ -195,21 +195,21 @@ void rcond_2x2() { Mat2 m_inverse = lu.inverse(); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond_est = lu.rcond(); - VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3); } { FullPivLU lu(m); Mat2 m_inverse = lu.inverse(); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond_est = lu.rcond(); - VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3); } { LLT llt(m); Mat2 m_inverse = llt.solve(Mat2::Identity()); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond_est = llt.rcond(); - VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3); } }