devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Improve ConditionEstimator docs and tighten test bounds

Browse Source 
libeigen/eigen!2226

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-02-26 22:08:13 -08:00
parent 8525491eb1
commit f64d1e0acc
2 changed files with 34 additions and 43 deletions
Download Patch File Download Diff File Expand all files Collapse all files

61
Eigen/src/Core/ConditionEstimator.h
View File

@@ -40,18 +40,17 @@ struct rcond_compute_sign<Vector, Vector, false> {
* \a matrix that implements .solve() and .adjoint().solve() methods. * \a matrix that implements .solve() and .adjoint().solve() methods.
* *
* This function implements Algorithms 4.1 and 5.1 from * This function implements Algorithms 4.1 and 5.1 from
* http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf * Higham, "Experience with a Matrix Norm Estimator",
* which also forms the basis for the condition number estimators in * SIAM J. Sci. Stat. Comput., 11(4):804-809, 1990.
* LAPACK. Since at most 10 calls to the solve method of dec are * with Higham's alternating-sign safety-net estimate from
* performed, the total cost is O(dims^2), as opposed to O(dims^3) * Higham and Tisseur, "A Block Algorithm for Matrix 1-Norm Estimation,
* needed to compute the inverse matrix explicitly. * 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 * The Hager/Higham gradient ascent uses at most 5 iterations of 2 solves
* ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be * each, giving a total cost of O(n^2).
* computed directly in O(n^2) operations.
* *
* Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, and * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, LLT.
* LLT.
* *
* \sa 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()); eigen_assert(dec.rows() == dec.cols());
const Index n = dec.rows(); const Index n = dec.rows();
if (n == 0) return 0; if (n == 0) return RealScalar(0);
// Disable Index to float conversion warning // Disable Index to float conversion warning
#ifdef __INTEL_COMPILER #ifdef __INTEL_COMPILER
@@ -80,14 +79,12 @@ typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomp
// lower_bound is a lower bound on // lower_bound is a lower bound on
// ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1
// and is the objective maximized by the ("super-") gradient ascent // and is the objective maximized by the supergradient ascent algorithm below.
// algorithm below.
RealScalar lower_bound = v.template lpNorm<1>(); RealScalar lower_bound = v.template lpNorm<1>();
if (n == 1) return lower_bound; if (n == 1) return lower_bound;
// Gradient ascent algorithm follows: We know that the optimum is achieved at // Gradient ascent: the optimum is achieved at a unit vector e_j. Each
// one of the simplices v = e_i, so in each iteration we follow a // iteration follows the supergradient to find which unit vector to probe next.
// super-gradient to move towards the optimal one.
RealScalar old_lower_bound = lower_bound; RealScalar old_lower_bound = lower_bound;
Vector sign_vector(n); Vector sign_vector(n);
Vector old_sign_vector; 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) { for (int k = 0; k < 4; ++k) {
sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v); sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v);
if (k > 0 && !is_complex && sign_vector == old_sign_vector) { if (k > 0 && !is_complex && sign_vector == old_sign_vector) {
// Break if the solution stagnated. // Break if the sign vector stagnated.
break; 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 = dec.adjoint().solve(sign_vector);
v.real().cwiseAbs().maxCoeff(&v_max_abs_index); v.real().cwiseAbs().maxCoeff(&v_max_abs_index);
if (v_max_abs_index == old_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; break;
} }
// Move to the new simplex e_j, where j = v_max_abs_index. // Probe the best unit vector: v = A^{-1} * e_j.
v = dec.solve(Vector::Unit(n, v_max_abs_index)); // v = inv(matrix) * e_j. v = dec.solve(Vector::Unit(n, v_max_abs_index));
lower_bound = v.template lpNorm<1>(); lower_bound = v.template lpNorm<1>();
if (lower_bound <= old_lower_bound) { if (lower_bound <= old_lower_bound) {
// Break if the gradient step did not increase the lower_bound. // No improvement from the gradient step.
break; break;
} }
if (!is_complex) { 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_v_max_abs_index = v_max_abs_index;
old_lower_bound = lower_bound; old_lower_bound = lower_bound;
} }
// The following calculates an independent estimate of ||matrix||_1 by // Higham's alternating-sign estimate: an independent safety-net that catches
// multiplying matrix by a vector with entries of slowly increasing // cases where the gradient ascent converges to a local maximum due to exact
// magnitude and alternating sign: // cancellation patterns (especially with permutations and backsubstitutions).
// v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. // v_i = (-1)^i * (1 + i/(n-1)), then estimate = 2*||A^{-1}*v||_1 / (3*n).
// 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(RealScalar(1)); Scalar alternating_sign(RealScalar(1));
for (Index i = 0; i < n; ++i) { 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>(RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1)))); v[i] = alternating_sign * static_cast<RealScalar>(RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1))));
alternating_sign = -alternating_sign; alternating_sign = -alternating_sign;
} }
v = dec.solve(v); v = dec.solve(v);
const RealScalar alternate_lower_bound = (2 * v.template lpNorm<1>()) / (3 * RealScalar(n)); const RealScalar alt_est = (RealScalar(2) * v.template lpNorm<1>()) / (RealScalar(3) * RealScalar(n));
return numext::maxi(lower_bound, alternate_lower_bound); return numext::maxi(lower_bound, alt_est);
} }
/** \brief Reciprocal condition number estimator. /** \brief Reciprocal condition number estimator.

16
test/condition_estimator.cpp
View File

@@ -29,8 +29,8 @@ void rcond_partial_piv_lu() {
MatrixType m_inverse = lu.inverse(); MatrixType m_inverse = lu.inverse();
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
RealScalar rcond_est = lu.rcond(); RealScalar rcond_est = lu.rcond();
// Verify the estimate is within a factor of 10 of the truth. // Verify the estimate is within a factor of 3 of the truth.
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3);
} }
template <typename MatrixType> template <typename MatrixType>
@@ -47,7 +47,7 @@ void rcond_full_piv_lu() {
MatrixType m_inverse = lu.inverse(); MatrixType m_inverse = lu.inverse();
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
RealScalar rcond_est = lu.rcond(); RealScalar rcond_est = lu.rcond();
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3);
} }
template <typename MatrixType> template <typename MatrixType>
@@ -65,7 +65,7 @@ void rcond_llt() {
MatrixType m_inverse = llt.solve(MatrixType::Identity(size, size)); MatrixType m_inverse = llt.solve(MatrixType::Identity(size, size));
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
RealScalar rcond_est = llt.rcond(); RealScalar rcond_est = llt.rcond();
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3);
} }
template <typename MatrixType> template <typename MatrixType>
@@ -83,7 +83,7 @@ void rcond_ldlt() {
MatrixType m_inverse = ldlt.solve(MatrixType::Identity(size, size)); MatrixType m_inverse = ldlt.solve(MatrixType::Identity(size, size));
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
RealScalar rcond_est = ldlt.rcond(); RealScalar rcond_est = ldlt.rcond();
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3);
} }
template <typename MatrixType> template <typename MatrixType>
@@ -195,21 +195,21 @@ void rcond_2x2() {
Mat2 m_inverse = lu.inverse(); Mat2 m_inverse = lu.inverse();
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
RealScalar rcond_est = lu.rcond(); RealScalar rcond_est = lu.rcond();
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3);
} }
{ {
FullPivLU<Mat2> lu(m); FullPivLU<Mat2> lu(m);
Mat2 m_inverse = lu.inverse(); Mat2 m_inverse = lu.inverse();
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
RealScalar rcond_est = lu.rcond(); RealScalar rcond_est = lu.rcond();
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3);
} }
{ {
LLT<Mat2> llt(m); LLT<Mat2> llt(m);
Mat2 m_inverse = llt.solve(Mat2::Identity()); Mat2 m_inverse = llt.solve(Mat2::Identity());
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse); RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m)) / matrix_l1_norm(m_inverse);
RealScalar rcond_est = llt.rcond(); RealScalar rcond_est = llt.rcond();
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est > rcond / 3 && rcond_est < rcond * 3);
} }
} }