diff options
| -rw-r--r-- | Eigen/Core | 9 | ||||
| -rw-r--r-- | Eigen/src/Core/ConditionEstimator.h | 279 | ||||
| -rw-r--r-- | Eigen/src/LU/FullPivLU.h | 13 | ||||
| -rw-r--r-- | Eigen/src/LU/PartialPivLU.h | 19 | ||||
| -rw-r--r-- | test/lu.cpp | 22 |
5 files changed, 332 insertions, 10 deletions
diff --git a/Eigen/Core b/Eigen/Core index 8428c51e4..0a196c814 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -33,13 +33,13 @@ #ifdef EIGEN_EXCEPTIONS #undef EIGEN_EXCEPTIONS #endif - + // All functions callable from CUDA code must be qualified with __device__ #define EIGEN_DEVICE_FUNC __host__ __device__ - + #else #define EIGEN_DEVICE_FUNC - + #endif #if defined(__CUDA_ARCH__) @@ -282,7 +282,7 @@ inline static const char *SimdInstructionSetsInUse(void) { // we use size_t frequently and we'll never remember to prepend it with std:: everytime just to // ensure QNX/QCC support using std::size_t; -// gcc 4.6.0 wants std:: for ptrdiff_t +// gcc 4.6.0 wants std:: for ptrdiff_t using std::ptrdiff_t; /** \defgroup Core_Module Core module @@ -422,6 +422,7 @@ using std::ptrdiff_t; #include "src/Core/products/TriangularSolverVector.h" #include "src/Core/BandMatrix.h" #include "src/Core/CoreIterators.h" +#include "src/Core/ConditionEstimator.h" #include "src/Core/BooleanRedux.h" #include "src/Core/Select.h" diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h new file mode 100644 index 000000000..ab6f59319 --- /dev/null +++ b/Eigen/src/Core/ConditionEstimator.h @@ -0,0 +1,279 @@ +// 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 Decomposition, bool IsComplex> +struct EstimateInverseL1NormImpl {}; +} // namespace internal + +template <typename Decomposition> +class ConditionEstimator { + public: + typedef typename Decomposition::MatrixType MatrixType; + typedef typename internal::traits<MatrixType>::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef typename internal::plain_col_type<MatrixType>::type Vector; + + /** \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); + } + RealScalar matrix_l1_norm = matrix.cwiseAbs().colwise().sum().maxCoeff(); + return rcond(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. + */ + 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 = EstimateInverseL1Norm(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) { + eigen_assert(dec.rows() == dec.cols()); + const int n = dec.rows(); + if (n == 0) { + return 0; + } + return internal::EstimateInverseL1NormImpl< + Decomposition, NumTraits<Scalar>::IsComplex>::compute(dec); + } +}; + +namespace internal { +// Partial specialization for real matrices. +template <typename Decomposition> +struct EstimateInverseL1NormImpl<Decomposition, 0> { + typedef typename Decomposition::MatrixType MatrixType; + typedef typename internal::traits<MatrixType>::Scalar Scalar; + typedef typename internal::plain_col_type<MatrixType>::type Vector; + + // Shorthand for vector L1 norm in Eigen. + inline static Scalar VectorL1Norm(const Vector& v) { + return v.template lpNorm<1>(); + } + + 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(A)||_1 = sup_v ||inv(A) 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(A)^T * sign_vector| + v = dec.transpose().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. + 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(A) * 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); + 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 ||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. + Scalar alternating_sign = 1; + for (int i = 0; i < n; ++i) { + v[i] = alternating_sign * static_cast<Scalar>(1) + + (static_cast<Scalar>(i) / (static_cast<Scalar>(n - 1))); + alternating_sign = -alternating_sign; + } + v = dec.solve(v); + const Scalar alternate_lower_bound = + (2 * VectorL1Norm(v)) / (3 * static_cast<Scalar>(n)); + return numext::maxi(lower_bound, alternate_lower_bound); + } +}; + +// Partial specialization for complex matrices. +template <typename Decomposition> +struct EstimateInverseL1NormImpl<Decomposition, 1> { + typedef typename Decomposition::MatrixType MatrixType; + typedef typename internal::traits<MatrixType>::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef typename internal::plain_col_type<MatrixType>::type Vector; + typedef typename internal::plain_col_type<MatrixType, RealScalar>::type + RealVector; + + // Shorthand for vector L1 norm in Eigen. + inline static 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; + } + // 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 + // 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(A)^* * sign_vector| + RealVector abs_v = v.cwiseAbs(); + const Vector psi = + (abs_v.array() == 0).select(v.cwiseQuotient(abs_v), ones); + v = dec.adjoint().solve(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(A) * 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; + } + // 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. + RealScalar alternating_sign = 1; + for (int i = 0; i < n; ++i) { + v[i] = alternating_sign * static_cast<RealScalar>(1) + + (static_cast<RealScalar>(i) / (static_cast<RealScalar>(n - 1))); + alternating_sign = -alternating_sign; + } + v = dec.solve(v); + const RealScalar alternate_lower_bound = + (2 * VectorL1Norm(v)) / (3 * static_cast<RealScalar>(n)); + return numext::maxi(lower_bound, alternate_lower_bound); + } +}; + +} // namespace internal +} // namespace Eigen + +#endif diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 1721213d6..ff0b78c35 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -231,6 +231,15 @@ template<typename _MatrixType> class FullPivLU return Solve<FullPivLU, Rhs>(*this, b.derived()); } + /** \returns an estimate of the reciprocal condition number of the matrix of which *this is + the LU decomposition. + */ + inline RealScalar rcond() const + { + eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); + return ConditionEstimator<FullPivLU<_MatrixType> >::rcond(m_l1_norm, *this); + } + /** \returns the determinant of the matrix of which * *this is the LU decomposition. It has only linear complexity * (that is, O(n) where n is the dimension of the square matrix) @@ -410,6 +419,7 @@ template<typename _MatrixType> class FullPivLU IntColVectorType m_rowsTranspositions; IntRowVectorType m_colsTranspositions; Index m_det_pq, m_nonzero_pivots; + RealScalar m_l1_norm; RealScalar m_maxpivot, m_prescribedThreshold; bool m_isInitialized, m_usePrescribedThreshold; }; @@ -455,11 +465,12 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const EigenBase<InputType> // the permutations are stored as int indices, so just to be sure: eigen_assert(matrix.rows()<=NumTraits<int>::highest() && matrix.cols()<=NumTraits<int>::highest()); - m_isInitialized = true; m_lu = matrix.derived(); + m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff(); computeInPlace(); + m_isInitialized = true; return *this; } diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index ab7797d2a..5d71a66d0 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -76,7 +76,6 @@ template<typename _MatrixType> class PartialPivLU typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; typedef typename MatrixType::PlainObject PlainObject; - /** * \brief Default Constructor. * @@ -152,6 +151,15 @@ template<typename _MatrixType> class PartialPivLU return Solve<PartialPivLU, Rhs>(*this, b.derived()); } + /** \returns an estimate of the reciprocal condition number of the matrix of which *this is + the LU decomposition. + */ + inline RealScalar rcond() const + { + eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); + return ConditionEstimator<PartialPivLU<_MatrixType> >::rcond(m_l1_norm, *this); + } + /** \returns the inverse of the matrix of which *this is the LU decomposition. * * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for @@ -178,7 +186,7 @@ template<typename _MatrixType> class PartialPivLU * * \sa MatrixBase::determinant() */ - typename internal::traits<MatrixType>::Scalar determinant() const; + Scalar determinant() const; MatrixType reconstructedMatrix() const; @@ -247,6 +255,7 @@ template<typename _MatrixType> class PartialPivLU PermutationType m_p; TranspositionType m_rowsTranspositions; Index m_det_p; + RealScalar m_l1_norm; bool m_isInitialized; }; @@ -256,6 +265,7 @@ PartialPivLU<MatrixType>::PartialPivLU() m_p(), m_rowsTranspositions(), m_det_p(0), + m_l1_norm(0), m_isInitialized(false) { } @@ -266,6 +276,7 @@ PartialPivLU<MatrixType>::PartialPivLU(Index size) m_p(size), m_rowsTranspositions(size), m_det_p(0), + m_l1_norm(0), m_isInitialized(false) { } @@ -277,6 +288,7 @@ PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix) m_p(matrix.rows()), m_rowsTranspositions(matrix.rows()), m_det_p(0), + m_l1_norm(0), m_isInitialized(false) { compute(matrix.derived()); @@ -467,6 +479,7 @@ PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const EigenBase<Inpu eigen_assert(matrix.rows()<NumTraits<int>::highest()); m_lu = matrix.derived(); + m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff(); eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); const Index size = matrix.rows(); @@ -484,7 +497,7 @@ PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const EigenBase<Inpu } template<typename MatrixType> -typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const +typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const { eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); return Scalar(m_det_p) * m_lu.diagonal().prod(); diff --git a/test/lu.cpp b/test/lu.cpp index f14435114..31991520d 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -11,6 +11,11 @@ #include <Eigen/LU> using namespace std; +template<typename MatrixType> +typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) { + return m.cwiseAbs().colwise().sum().maxCoeff(); +} + template<typename MatrixType> void lu_non_invertible() { typedef typename MatrixType::Index Index; @@ -143,7 +148,13 @@ template<typename MatrixType> void lu_invertible() m3 = MatrixType::Random(size,size); m2 = lu.solve(m3); VERIFY_IS_APPROX(m3, m1*m2); - VERIFY_IS_APPROX(m2, lu.inverse()*m3); + MatrixType m1_inverse = lu.inverse(); + VERIFY_IS_APPROX(m2, m1_inverse*m3); + + // Test condition number estimation. + 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); // test solve with transposed lu.template _solve_impl_transposed<false>(m3, m2); @@ -170,6 +181,7 @@ template<typename MatrixType> void lu_partial_piv() PartialPivLU.h */ typedef typename MatrixType::Index Index; + typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; Index size = internal::random<Index>(1,4); MatrixType m1(size, size), m2(size, size), m3(size, size); @@ -181,7 +193,13 @@ template<typename MatrixType> void lu_partial_piv() m3 = MatrixType::Random(size,size); m2 = plu.solve(m3); VERIFY_IS_APPROX(m3, m1*m2); - VERIFY_IS_APPROX(m2, plu.inverse()*m3); + MatrixType m1_inverse = plu.inverse(); + VERIFY_IS_APPROX(m2, m1_inverse*m3); + + // Test condition number estimation. + 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); // test solve with transposed plu.template _solve_impl_transposed<false>(m3, m2); |