src/Core/ConditionEstimator.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2016 Rasmus Munk Larsen (rmlarsen@google.com)
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_CONDITIONESTIMATOR_H
0011 #define EIGEN_CONDITIONESTIMATOR_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 template <typename Vector, typename RealVector, bool IsComplex>
0018 struct rcond_compute_sign {
0019   static inline Vector run(const Vector& v) {
0020     const RealVector v_abs = v.cwiseAbs();
0021     return (v_abs.array() == static_cast<typename Vector::RealScalar>(0))
0022             .select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs));
0023   }
0024 };
0025
0026 // Partial specialization to avoid elementwise division for real vectors.
0027 template <typename Vector>
0028 struct rcond_compute_sign<Vector, Vector, false> {
0029   static inline Vector run(const Vector& v) {
0030     return (v.array() < static_cast<typename Vector::RealScalar>(0))
0031            .select(-Vector::Ones(v.size()), Vector::Ones(v.size()));
0032   }
0033 };
0034
0035 /**
0036   * \returns an estimate of ||inv(matrix)||_1 given a decomposition of
0037   * \a matrix that implements .solve() and .adjoint().solve() methods.
0038   *
0039   * This function implements Algorithms 4.1 and 5.1 from
0040   *   http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf
0041   * which also forms the basis for the condition number estimators in
0042   * LAPACK. Since at most 10 calls to the solve method of dec are
0043   * performed, the total cost is O(dims^2), as opposed to O(dims^3)
0044   * needed to compute the inverse matrix explicitly.
0045   *
0046   * The most common usage is in estimating the condition number
0047   * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be
0048   * computed directly in O(n^2) operations.
0049   *
0050   * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, and
0051   * LLT.
0052   *
0053   * \sa FullPivLU, PartialPivLU, LDLT, LLT.
0054   */
0055 template <typename Decomposition>
0056 typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomposition& dec)
0057 {
0058   typedef typename Decomposition::MatrixType MatrixType;
0059   typedef typename Decomposition::Scalar Scalar;
0060   typedef typename Decomposition::RealScalar RealScalar;
0061   typedef typename internal::plain_col_type<MatrixType>::type Vector;
0062   typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVector;
0063   const bool is_complex = (NumTraits<Scalar>::IsComplex != 0);
0064
0065   eigen_assert(dec.rows() == dec.cols());
0066   const Index n = dec.rows();
0067   if (n == 0)
0068     return 0;
0069
0070   // Disable Index to float conversion warning
0071 #ifdef __INTEL_COMPILER
0072   #pragma warning push
0073   #pragma warning ( disable : 2259 )
0074 #endif
0075   Vector v = dec.solve(Vector::Ones(n) / Scalar(n));
0076 #ifdef __INTEL_COMPILER
0077   #pragma warning pop
0078 #endif
0079
0080   // lower_bound is a lower bound on
0081   //   ||inv(matrix)||_1  = sup_v ||inv(matrix) v||_1 / ||v||_1
0082   // and is the objective maximized by the ("super-") gradient ascent
0083   // algorithm below.
0084   RealScalar lower_bound = v.template lpNorm<1>();
0085   if (n == 1)
0086     return lower_bound;
0087
0088   // Gradient ascent algorithm follows: We know that the optimum is achieved at
0089   // one of the simplices v = e_i, so in each iteration we follow a
0090   // super-gradient to move towards the optimal one.
0091   RealScalar old_lower_bound = lower_bound;
0092   Vector sign_vector(n);
0093   Vector old_sign_vector;
0094   Index v_max_abs_index = -1;
0095   Index old_v_max_abs_index = v_max_abs_index;
0096   for (int k = 0; k < 4; ++k)
0097   {
0098     sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v);
0099     if (k > 0 && !is_complex && sign_vector == old_sign_vector) {
0100       // Break if the solution stagnated.
0101       break;
0102     }
0103     // v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )|
0104     v = dec.adjoint().solve(sign_vector);
0105     v.real().cwiseAbs().maxCoeff(&v_max_abs_index);
0106     if (v_max_abs_index == old_v_max_abs_index) {
0107       // Break if the solution stagnated.
0108       break;
0109     }
0110     // Move to the new simplex e_j, where j = v_max_abs_index.
0111     v = dec.solve(Vector::Unit(n, v_max_abs_index));  // v = inv(matrix) * e_j.
0112     lower_bound = v.template lpNorm<1>();
0113     if (lower_bound <= old_lower_bound) {
0114       // Break if the gradient step did not increase the lower_bound.
0115       break;
0116     }
0117     if (!is_complex) {
0118       old_sign_vector = sign_vector;
0119     }
0120     old_v_max_abs_index = v_max_abs_index;
0121     old_lower_bound = lower_bound;
0122   }
0123   // The following calculates an independent estimate of ||matrix||_1 by
0124   // multiplying matrix by a vector with entries of slowly increasing
0125   // magnitude and alternating sign:
0126   //   v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1.
0127   // This improvement to Hager's algorithm above is due to Higham. It was
0128   // added to make the algorithm more robust in certain corner cases where
0129   // large elements in the matrix might otherwise escape detection due to
0130   // exact cancellation (especially when op and op_adjoint correspond to a
0131   // sequence of backsubstitutions and permutations), which could cause
0132   // Hager's algorithm to vastly underestimate ||matrix||_1.
0133   Scalar alternating_sign(RealScalar(1));
0134   for (Index i = 0; i < n; ++i) {
0135     // The static_cast is needed when Scalar is a complex and RealScalar implements expression templates
0136     v[i] = alternating_sign * static_cast<RealScalar>(RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1))));
0137     alternating_sign = -alternating_sign;
0138   }
0139   v = dec.solve(v);
0140   const RealScalar alternate_lower_bound = (2 * v.template lpNorm<1>()) / (3 * RealScalar(n));
0141   return numext::maxi(lower_bound, alternate_lower_bound);
0142 }
0143
0144 /** \brief Reciprocal condition number estimator.
0145   *
0146   * Computing a decomposition of a dense matrix takes O(n^3) operations, while
0147   * this method estimates the condition number quickly and reliably in O(n^2)
0148   * operations.
0149   *
0150   * \returns an estimate of the reciprocal condition number
0151   * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and
0152   * its decomposition. Supports the following decompositions: FullPivLU,
0153   * PartialPivLU, LDLT, and LLT.
0154   *
0155   * \sa FullPivLU, PartialPivLU, LDLT, LLT.
0156   */
0157 template <typename Decomposition>
0158 typename Decomposition::RealScalar
0159 rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Decomposition& dec)
0160 {
0161   typedef typename Decomposition::RealScalar RealScalar;
0162   eigen_assert(dec.rows() == dec.cols());
0163   if (dec.rows() == 0)              return NumTraits<RealScalar>::infinity();
0164   if (matrix_norm == RealScalar(0)) return RealScalar(0);
0165   if (dec.rows() == 1)              return RealScalar(1);
0166   const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec);
0167   return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0)
0168                                                : (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
0169 }
0170
0171 }  // namespace internal
0172
0173 }  // namespace Eigen
0174
0175 #endif