1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2016 Rasmus Munk Larsen (rmlarsen (at) google.com) 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_CONDITIONESTIMATOR_H 11 #define EIGEN_CONDITIONESTIMATOR_H 12 13 namespace Eigen { 14 15 namespace internal { 16 17 template <typename Vector, typename RealVector, bool IsComplex> 18 struct rcond_compute_sign { 19 static inline Vector run(const Vector& v) { 20 const RealVector v_abs = v.cwiseAbs(); 21 return (v_abs.array() == static_cast<typename Vector::RealScalar>(0)) 22 .select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs)); 23 } 24 }; 25 26 // Partial specialization to avoid elementwise division for real vectors. 27 template <typename Vector> 28 struct rcond_compute_sign<Vector, Vector, false> { 29 static inline Vector run(const Vector& v) { 30 return (v.array() < static_cast<typename Vector::RealScalar>(0)) 31 .select(-Vector::Ones(v.size()), Vector::Ones(v.size())); 32 } 33 }; 34 35 /** 36 * \returns an estimate of ||inv(matrix)||_1 given a decomposition of 37 * \a matrix that implements .solve() and .adjoint().solve() methods. 38 * 39 * This function implements Algorithms 4.1 and 5.1 from 40 * http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf 41 * which also forms the basis for the condition number estimators in 42 * LAPACK. Since at most 10 calls to the solve method of dec are 43 * performed, the total cost is O(dims^2), as opposed to O(dims^3) 44 * needed to compute the inverse matrix explicitly. 45 * 46 * The most common usage is in estimating the condition number 47 * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be 48 * computed directly in O(n^2) operations. 49 * 50 * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, and 51 * LLT. 52 * 53 * \sa FullPivLU, PartialPivLU, LDLT, LLT. 54 */ 55 template <typename Decomposition> 56 typename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomposition& dec) 57 { 58 typedef typename Decomposition::MatrixType MatrixType; 59 typedef typename Decomposition::Scalar Scalar; 60 typedef typename Decomposition::RealScalar RealScalar; 61 typedef typename internal::plain_col_type<MatrixType>::type Vector; 62 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVector; 63 const bool is_complex = (NumTraits<Scalar>::IsComplex != 0); 64 65 eigen_assert(dec.rows() == dec.cols()); 66 const Index n = dec.rows(); 67 if (n == 0) 68 return 0; 69 70 // Disable Index to float conversion warning 71 #ifdef __INTEL_COMPILER 72 #pragma warning push 73 #pragma warning ( disable : 2259 ) 74 #endif 75 Vector v = dec.solve(Vector::Ones(n) / Scalar(n)); 76 #ifdef __INTEL_COMPILER 77 #pragma warning pop 78 #endif 79 80 // lower_bound is a lower bound on 81 // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 82 // and is the objective maximized by the ("super-") gradient ascent 83 // algorithm below. 84 RealScalar lower_bound = v.template lpNorm<1>(); 85 if (n == 1) 86 return lower_bound; 87 88 // Gradient ascent algorithm follows: We know that the optimum is achieved at 89 // one of the simplices v = e_i, so in each iteration we follow a 90 // super-gradient to move towards the optimal one. 91 RealScalar old_lower_bound = lower_bound; 92 Vector sign_vector(n); 93 Vector old_sign_vector; 94 Index v_max_abs_index = -1; 95 Index old_v_max_abs_index = v_max_abs_index; 96 for (int k = 0; k < 4; ++k) 97 { 98 sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v); 99 if (k > 0 && !is_complex && sign_vector == old_sign_vector) { 100 // Break if the solution stagnated. 101 break; 102 } 103 // v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )| 104 v = dec.adjoint().solve(sign_vector); 105 v.real().cwiseAbs().maxCoeff(&v_max_abs_index); 106 if (v_max_abs_index == old_v_max_abs_index) { 107 // Break if the solution stagnated. 108 break; 109 } 110 // Move to the new simplex e_j, where j = v_max_abs_index. 111 v = dec.solve(Vector::Unit(n, v_max_abs_index)); // v = inv(matrix) * e_j. 112 lower_bound = v.template lpNorm<1>(); 113 if (lower_bound <= old_lower_bound) { 114 // Break if the gradient step did not increase the lower_bound. 115 break; 116 } 117 if (!is_complex) { 118 old_sign_vector = sign_vector; 119 } 120 old_v_max_abs_index = v_max_abs_index; 121 old_lower_bound = lower_bound; 122 } 123 // The following calculates an independent estimate of ||matrix||_1 by 124 // multiplying matrix by a vector with entries of slowly increasing 125 // magnitude and alternating sign: 126 // v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. 127 // This improvement to Hager's algorithm above is due to Higham. It was 128 // added to make the algorithm more robust in certain corner cases where 129 // large elements in the matrix might otherwise escape detection due to 130 // exact cancellation (especially when op and op_adjoint correspond to a 131 // sequence of backsubstitutions and permutations), which could cause 132 // Hager's algorithm to vastly underestimate ||matrix||_1. 133 Scalar alternating_sign(RealScalar(1)); 134 for (Index i = 0; i < n; ++i) { 135 // The static_cast is needed when Scalar is a complex and RealScalar implements expression templates 136 v[i] = alternating_sign * static_cast<RealScalar>(RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1)))); 137 alternating_sign = -alternating_sign; 138 } 139 v = dec.solve(v); 140 const RealScalar alternate_lower_bound = (2 * v.template lpNorm<1>()) / (3 * RealScalar(n)); 141 return numext::maxi(lower_bound, alternate_lower_bound); 142 } 143 144 /** \brief Reciprocal condition number estimator. 145 * 146 * Computing a decomposition of a dense matrix takes O(n^3) operations, while 147 * this method estimates the condition number quickly and reliably in O(n^2) 148 * operations. 149 * 150 * \returns an estimate of the reciprocal condition number 151 * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and 152 * its decomposition. Supports the following decompositions: FullPivLU, 153 * PartialPivLU, LDLT, and LLT. 154 * 155 * \sa FullPivLU, PartialPivLU, LDLT, LLT. 156 */ 157 template <typename Decomposition> 158 typename Decomposition::RealScalar 159 rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Decomposition& dec) 160 { 161 typedef typename Decomposition::RealScalar RealScalar; 162 eigen_assert(dec.rows() == dec.cols()); 163 if (dec.rows() == 0) return RealScalar(1); 164 if (matrix_norm == RealScalar(0)) return RealScalar(0); 165 if (dec.rows() == 1) return RealScalar(1); 166 const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec); 167 return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0) 168 : (RealScalar(1) / inverse_matrix_norm) / matrix_norm); 169 } 170 171 } // namespace internal 172 173 } // namespace Eigen 174 175 #endif 176
