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/Core/ConditionEstimator.h | 405 ++++++++++++++---------------------- 1 file changed, 152 insertions(+), 253 deletions(-) (limited to 'Eigen/src/Core/ConditionEstimator.h') 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 -- cgit v1.2.3