mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-15 07:10:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix comments in ConditionEstimator and minor cleanup.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2016-04-01 11:58:17 -07:00
parent 1aa89fb855
commit 91414e0042
2 changed files with 65 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -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) +

4
test/lu.cpp
View File

@ -152,7 +152,7 @@ template<typename MatrixType> void lu_invertible()
VERIFY_IS_APPROX(m2, m1_inverse*m3);
// Test condition number estimation.
RealScalar rcond = RealScalar(1) / matrix_l1_norm(m1) / matrix_l1_norm(m1_inverse);
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse);
// Verify that the estimate is within a factor of 10 of the truth.
VERIFY(lu.rcond() > rcond / 10 && lu.rcond() < rcond * 10);
@ -197,7 +197,7 @@ template<typename MatrixType> void lu_partial_piv()
VERIFY_IS_APPROX(m2, m1_inverse*m3);
// Test condition number estimation.
RealScalar rcond = RealScalar(1) / matrix_l1_norm(m1) / matrix_l1_norm(m1_inverse);
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse);
// Verify that the estimate is within a factor of 10 of the truth.
VERIFY(plu.rcond() > rcond / 10 && plu.rcond() < rcond * 10);