// 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/.
/** * \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)
{ typedeftypename Decomposition::MatrixType MatrixType; typedeftypename Decomposition::Scalar Scalar; typedeftypename Decomposition::RealScalar RealScalar; typedeftypename internal::plain_col_type<MatrixType>::type Vector; typedeftypename internal::plain_col_type<MatrixType, RealScalar>::type RealVector; constbool 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)
{ typedeftypename 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);
}
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
