123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175 |
- // This file is part of Eigen, a lightweight C++ template library
- // for linear algebra.
- //
- // Copyright (C) 2016 Rasmus Munk Larsen (rmlarsen@google.com)
- //
- // This Source Code Form is subject to the terms of the Mozilla
- // Public License v. 2.0. If a copy of the MPL was not distributed
- // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
- #ifndef EIGEN_CONDITIONESTIMATOR_H
- #define EIGEN_CONDITIONESTIMATOR_H
- namespace Eigen {
- namespace internal {
- template <typename Vector, typename RealVector, bool IsComplex>
- struct rcond_compute_sign {
- static inline Vector run(const Vector& v) {
- const RealVector v_abs = v.cwiseAbs();
- return (v_abs.array() == static_cast<typename Vector::RealScalar>(0))
- .select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs));
- }
- };
- // Partial specialization to avoid elementwise division for real vectors.
- template <typename Vector>
- struct rcond_compute_sign<Vector, Vector, false> {
- static inline Vector run(const Vector& v) {
- return (v.array() < static_cast<typename Vector::RealScalar>(0))
- .select(-Vector::Ones(v.size()), Vector::Ones(v.size()));
- }
- };
- /**
- * \returns an estimate of ||inv(matrix)||_1 given a decomposition of
- * \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.
- *
- * 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.
- *
- * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, and
- * LLT.
- *
- * \sa FullPivLU, PartialPivLU, LDLT, LLT.
- */
- template <typename Decomposition>
- typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomposition& dec)
- {
- typedef typename Decomposition::MatrixType MatrixType;
- typedef typename Decomposition::Scalar Scalar;
- typedef typename Decomposition::RealScalar RealScalar;
- typedef typename internal::plain_col_type<MatrixType>::type Vector;
- typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVector;
- const bool is_complex = (NumTraits<Scalar>::IsComplex != 0);
- eigen_assert(dec.rows() == dec.cols());
- const Index n = dec.rows();
- if (n == 0)
- return 0;
- // Disable Index to float conversion warning
- #ifdef __INTEL_COMPILER
- #pragma warning push
- #pragma warning ( disable : 2259 )
- #endif
- Vector v = dec.solve(Vector::Ones(n) / Scalar(n));
- #ifdef __INTEL_COMPILER
- #pragma warning pop
- #endif
- // 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 = 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.
- RealScalar old_lower_bound = lower_bound;
- Vector sign_vector(n);
- Vector old_sign_vector;
- Index v_max_abs_index = -1;
- Index old_v_max_abs_index = v_max_abs_index;
- for (int k = 0; k < 4; ++k)
- {
- sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v);
- if (k > 0 && !is_complex && sign_vector == old_sign_vector) {
- // Break if the solution stagnated.
- break;
- }
- // 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;
- }
- // 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 = v.template lpNorm<1>();
- if (lower_bound <= old_lower_bound) {
- // Break if the gradient step did not increase the lower_bound.
- break;
- }
- if (!is_complex) {
- 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(RealScalar(1));
- for (Index i = 0; i < n; ++i) {
- // The static_cast is needed when Scalar is a complex and RealScalar implements expression templates
- v[i] = alternating_sign * static_cast<RealScalar>(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);
- }
- /** \brief Reciprocal 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, LDLT, and LLT.
- *
- * \sa FullPivLU, PartialPivLU, LDLT, LLT.
- */
- template <typename Decomposition>
- typename Decomposition::RealScalar
- rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Decomposition& dec)
- {
- typedef typename Decomposition::RealScalar RealScalar;
- eigen_assert(dec.rows() == dec.cols());
- if (dec.rows() == 0) return NumTraits<RealScalar>::infinity();
- if (matrix_norm == RealScalar(0)) return RealScalar(0);
- if (dec.rows() == 1) return RealScalar(1);
- const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec);
- return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0)
- : (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
- }
- } // namespace internal
- } // namespace Eigen
- #endif
|