![]() |
Eigen
3.3.3
|
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