diff options
author | Rasmus Munk Larsen <rmlarsen@google.com> | 2016-04-01 11:58:17 -0700 |
---|---|---|
committer | Rasmus Munk Larsen <rmlarsen@google.com> | 2016-04-01 11:58:17 -0700 |
commit | 91414e0042779b1b9d312d9255f389e67aa38106 (patch) | |
tree | 7ad8bd061a431e5e66e17ea8d7c87b2872bfe3f6 /Eigen/src/Core/ConditionEstimator.h | |
parent | 1aa89fb85548dc425d54d2cbe7f28915c29db13a (diff) |
Fix comments in ConditionEstimator and minor cleanup.
Diffstat (limited to 'Eigen/src/Core/ConditionEstimator.h')
-rw-r--r-- | Eigen/src/Core/ConditionEstimator.h | 119 |
1 files changed, 63 insertions, 56 deletions
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 <typename Decomposition, bool IsComplex> -struct EstimateInverseL1NormImpl {}; +struct EstimateInverseMatrixL1NormImpl {}; } // namespace internal template <typename Decomposition> @@ -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<Scalar>::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 <typename Decomposition> -struct EstimateInverseL1NormImpl<Decomposition, 0> { +struct EstimateInverseMatrixL1NormImpl<Decomposition, 0> { typedef typename Decomposition::MatrixType MatrixType; typedef typename internal::traits<MatrixType>::Scalar Scalar; typedef typename internal::plain_col_type<MatrixType>::type Vector; @@ -130,8 +137,9 @@ struct EstimateInverseL1NormImpl<Decomposition, 0> { 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<Decomposition, 0> { 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<Decomposition, 0> { // 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<Decomposition, 0> { 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<Scalar>(1) + @@ -194,7 +201,7 @@ struct EstimateInverseL1NormImpl<Decomposition, 0> { // Partial specialization for complex matrices. template <typename Decomposition> -struct EstimateInverseL1NormImpl<Decomposition, 1> { +struct EstimateInverseMatrixL1NormImpl<Decomposition, 1> { typedef typename Decomposition::MatrixType MatrixType; typedef typename internal::traits<MatrixType>::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; @@ -216,8 +223,9 @@ struct EstimateInverseL1NormImpl<Decomposition, 1> { 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<Decomposition, 1> { 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<Decomposition, 1> { // 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<Decomposition, 1> { 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<RealScalar>(1) + |