Eigen  3.3.3
ConditionEstimator.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2016 Rasmus Munk Larsen (rmlarsen@google.com)
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_CONDITIONESTIMATOR_H
00011 #define EIGEN_CONDITIONESTIMATOR_H
00012 
00013 namespace Eigen {
00014 
00015 namespace internal {
00016 
00017 template <typename Vector, typename RealVector, bool IsComplex>
00018 struct rcond_compute_sign {
00019   static inline Vector run(const Vector& v) {
00020     const RealVector v_abs = v.cwiseAbs();
00021     return (v_abs.array() == static_cast<typename Vector::RealScalar>(0))
00022             .select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs));
00023   }
00024 };
00025 
00026 // Partial specialization to avoid elementwise division for real vectors.
00027 template <typename Vector>
00028 struct rcond_compute_sign<Vector, Vector, false> {
00029   static inline Vector run(const Vector& v) {
00030     return (v.array() < static_cast<typename Vector::RealScalar>(0))
00031            .select(-Vector::Ones(v.size()), Vector::Ones(v.size()));
00032   }
00033 };
00034 
00055 template <typename Decomposition>
00056 typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomposition& dec)
00057 {
00058   typedef typename Decomposition::MatrixType MatrixType;
00059   typedef typename Decomposition::Scalar Scalar;
00060   typedef typename Decomposition::RealScalar RealScalar;
00061   typedef typename internal::plain_col_type<MatrixType>::type Vector;
00062   typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVector;
00063   const bool is_complex = (NumTraits<Scalar>::IsComplex != 0);
00064 
00065   eigen_assert(dec.rows() == dec.cols());
00066   const Index n = dec.rows();
00067   if (n == 0)
00068     return 0;
00069 
00070   // Disable Index to float conversion warning
00071 #ifdef __INTEL_COMPILER
00072   #pragma warning push
00073   #pragma warning ( disable : 2259 )
00074 #endif
00075   Vector v = dec.solve(Vector::Ones(n) / Scalar(n));
00076 #ifdef __INTEL_COMPILER
00077   #pragma warning pop
00078 #endif
00079 
00080   // lower_bound is a lower bound on
00081   //   ||inv(matrix)||_1  = sup_v ||inv(matrix) v||_1 / ||v||_1
00082   // and is the objective maximized by the ("super-") gradient ascent
00083   // algorithm below.
00084   RealScalar lower_bound = v.template lpNorm<1>();
00085   if (n == 1)
00086     return lower_bound;
00087 
00088   // Gradient ascent algorithm follows: We know that the optimum is achieved at
00089   // one of the simplices v = e_i, so in each iteration we follow a
00090   // super-gradient to move towards the optimal one.
00091   RealScalar old_lower_bound = lower_bound;
00092   Vector sign_vector(n);
00093   Vector old_sign_vector;
00094   Index v_max_abs_index = -1;
00095   Index old_v_max_abs_index = v_max_abs_index;
00096   for (int k = 0; k < 4; ++k)
00097   {
00098     sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v);
00099     if (k > 0 && !is_complex && sign_vector == old_sign_vector) {
00100       // Break if the solution stagnated.
00101       break;
00102     }
00103     // v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )|
00104     v = dec.adjoint().solve(sign_vector);
00105     v.real().cwiseAbs().maxCoeff(&v_max_abs_index);
00106     if (v_max_abs_index == old_v_max_abs_index) {
00107       // Break if the solution stagnated.
00108       break;
00109     }
00110     // Move to the new simplex e_j, where j = v_max_abs_index.
00111     v = dec.solve(Vector::Unit(n, v_max_abs_index));  // v = inv(matrix) * e_j.
00112     lower_bound = v.template lpNorm<1>();
00113     if (lower_bound <= old_lower_bound) {
00114       // Break if the gradient step did not increase the lower_bound.
00115       break;
00116     }
00117     if (!is_complex) {
00118       old_sign_vector = sign_vector;
00119     }
00120     old_v_max_abs_index = v_max_abs_index;
00121     old_lower_bound = lower_bound;
00122   }
00123   // The following calculates an independent estimate of ||matrix||_1 by
00124   // multiplying matrix by a vector with entries of slowly increasing
00125   // magnitude and alternating sign:
00126   //   v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1.
00127   // This improvement to Hager's algorithm above is due to Higham. It was
00128   // added to make the algorithm more robust in certain corner cases where
00129   // large elements in the matrix might otherwise escape detection due to
00130   // exact cancellation (especially when op and op_adjoint correspond to a
00131   // sequence of backsubstitutions and permutations), which could cause
00132   // Hager's algorithm to vastly underestimate ||matrix||_1.
00133   Scalar alternating_sign(RealScalar(1));
00134   for (Index i = 0; i < n; ++i) {
00135     // The static_cast is needed when Scalar is a complex and RealScalar implements expression templates
00136     v[i] = alternating_sign * static_cast<RealScalar>(RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1))));
00137     alternating_sign = -alternating_sign;
00138   }
00139   v = dec.solve(v);
00140   const RealScalar alternate_lower_bound = (2 * v.template lpNorm<1>()) / (3 * RealScalar(n));
00141   return numext::maxi(lower_bound, alternate_lower_bound);
00142 }
00143 
00157 template <typename Decomposition>
00158 typename Decomposition::RealScalar
00159 rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Decomposition& dec)
00160 {
00161   typedef typename Decomposition::RealScalar RealScalar;
00162   eigen_assert(dec.rows() == dec.cols());
00163   if (dec.rows() == 0)              return RealScalar(1);
00164   if (matrix_norm == RealScalar(0)) return RealScalar(0);
00165   if (dec.rows() == 1)              return RealScalar(1);
00166   const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec);
00167   return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0)
00168                                                : (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
00169 }
00170 
00171 }  // namespace internal
00172 
00173 }  // namespace Eigen
00174 
00175 #endif
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends