12b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// This file is part of Eigen, a lightweight C++ template library 22b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// for linear algebra. 32b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// 42b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2016 Rasmus Munk Larsen (rmlarsen@google.com) 52b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// 62b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// This Source Code Form is subject to the terms of the Mozilla 72b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Public License v. 2.0. If a copy of the MPL was not distributed 82b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 92b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifndef EIGEN_CONDITIONESTIMATOR_H 112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#define EIGEN_CONDITIONESTIMATOR_H 122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace Eigen { 142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <typename Vector, typename RealVector, bool IsComplex> 182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct rcond_compute_sign { 192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static inline Vector run(const Vector& v) { 202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const RealVector v_abs = v.cwiseAbs(); 212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (v_abs.array() == static_cast<typename Vector::RealScalar>(0)) 222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang .select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs)); 232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Partial specialization to avoid elementwise division for real vectors. 272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <typename Vector> 282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct rcond_compute_sign<Vector, Vector, false> { 292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang static inline Vector run(const Vector& v) { 302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (v.array() < static_cast<typename Vector::RealScalar>(0)) 312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang .select(-Vector::Ones(v.size()), Vector::Ones(v.size())); 322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang/** 362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns an estimate of ||inv(matrix)||_1 given a decomposition of 372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \a matrix that implements .solve() and .adjoint().solve() methods. 382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * This function implements Algorithms 4.1 and 5.1 from 402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf 412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * which also forms the basis for the condition number estimators in 422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * LAPACK. Since at most 10 calls to the solve method of dec are 432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * performed, the total cost is O(dims^2), as opposed to O(dims^3) 442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * needed to compute the inverse matrix explicitly. 452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * The most common usage is in estimating the condition number 472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be 482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * computed directly in O(n^2) operations. 492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Supports the following decompositions: FullPivLU, PartialPivLU, LDLT, and 512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * LLT. 522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa FullPivLU, PartialPivLU, LDLT, LLT. 542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <typename Decomposition> 562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtypename Decomposition::RealScalar rcond_invmatrix_L1_norm_estimate(const Decomposition& dec) 572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename Decomposition::MatrixType MatrixType; 592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename Decomposition::Scalar Scalar; 602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename Decomposition::RealScalar RealScalar; 612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename internal::plain_col_type<MatrixType>::type Vector; 622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVector; 632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const bool is_complex = (NumTraits<Scalar>::IsComplex != 0); 642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(dec.rows() == dec.cols()); 662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index n = dec.rows(); 672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (n == 0) 682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return 0; 692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Disable Index to float conversion warning 712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifdef __INTEL_COMPILER 722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang #pragma warning push 732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang #pragma warning ( disable : 2259 ) 742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Vector v = dec.solve(Vector::Ones(n) / Scalar(n)); 762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifdef __INTEL_COMPILER 772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang #pragma warning pop 782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // lower_bound is a lower bound on 812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // ||inv(matrix)||_1 = sup_v ||inv(matrix) v||_1 / ||v||_1 822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // and is the objective maximized by the ("super-") gradient ascent 832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // algorithm below. 842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RealScalar lower_bound = v.template lpNorm<1>(); 852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (n == 1) 862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return lower_bound; 872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Gradient ascent algorithm follows: We know that the optimum is achieved at 892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // one of the simplices v = e_i, so in each iteration we follow a 902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // super-gradient to move towards the optimal one. 912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RealScalar old_lower_bound = lower_bound; 922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Vector sign_vector(n); 932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Vector old_sign_vector; 942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index v_max_abs_index = -1; 952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index old_v_max_abs_index = v_max_abs_index; 962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for (int k = 0; k < 4; ++k) 972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang sign_vector = internal::rcond_compute_sign<Vector, RealVector, is_complex>::run(v); 992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (k > 0 && !is_complex && sign_vector == old_sign_vector) { 1002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Break if the solution stagnated. 1012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang break; 1022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )| 1042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v = dec.adjoint().solve(sign_vector); 1052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v.real().cwiseAbs().maxCoeff(&v_max_abs_index); 1062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (v_max_abs_index == old_v_max_abs_index) { 1072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Break if the solution stagnated. 1082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang break; 1092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Move to the new simplex e_j, where j = v_max_abs_index. 1112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v = dec.solve(Vector::Unit(n, v_max_abs_index)); // v = inv(matrix) * e_j. 1122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang lower_bound = v.template lpNorm<1>(); 1132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (lower_bound <= old_lower_bound) { 1142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Break if the gradient step did not increase the lower_bound. 1152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang break; 1162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (!is_complex) { 1182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang old_sign_vector = sign_vector; 1192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang old_v_max_abs_index = v_max_abs_index; 1212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang old_lower_bound = lower_bound; 1222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // The following calculates an independent estimate of ||matrix||_1 by 1242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // multiplying matrix by a vector with entries of slowly increasing 1252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // magnitude and alternating sign: 1262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1. 1272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // This improvement to Hager's algorithm above is due to Higham. It was 1282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // added to make the algorithm more robust in certain corner cases where 1292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // large elements in the matrix might otherwise escape detection due to 1302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // exact cancellation (especially when op and op_adjoint correspond to a 1312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // sequence of backsubstitutions and permutations), which could cause 1322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Hager's algorithm to vastly underestimate ||matrix||_1. 1332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar alternating_sign(RealScalar(1)); 1342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for (Index i = 0; i < n; ++i) { 1352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // The static_cast is needed when Scalar is a complex and RealScalar implements expression templates 1362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v[i] = alternating_sign * static_cast<RealScalar>(RealScalar(1) + (RealScalar(i) / (RealScalar(n - 1)))); 1372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang alternating_sign = -alternating_sign; 1382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v = dec.solve(v); 1402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const RealScalar alternate_lower_bound = (2 * v.template lpNorm<1>()) / (3 * RealScalar(n)); 1412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return numext::maxi(lower_bound, alternate_lower_bound); 1422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 1432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang/** \brief Reciprocal condition number estimator. 1452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 1462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Computing a decomposition of a dense matrix takes O(n^3) operations, while 1472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * this method estimates the condition number quickly and reliably in O(n^2) 1482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * operations. 1492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 1502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns an estimate of the reciprocal condition number 1512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and 1522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * its decomposition. Supports the following decompositions: FullPivLU, 1532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * PartialPivLU, LDLT, and LLT. 1542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 1552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa FullPivLU, PartialPivLU, LDLT, LLT. 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 1572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <typename Decomposition> 1582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtypename Decomposition::RealScalar 1592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangrcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Decomposition& dec) 1602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 1612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename Decomposition::RealScalar RealScalar; 1622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(dec.rows() == dec.cols()); 1632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (dec.rows() == 0) return RealScalar(1); 1642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (matrix_norm == RealScalar(0)) return RealScalar(0); 1652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (dec.rows() == 1) return RealScalar(1); 1662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec); 1672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0) 1682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : (RealScalar(1) / inverse_matrix_norm) / matrix_norm); 1692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 1702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // namespace internal 1722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // namespace Eigen 1742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 176