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@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 {
runrcond_compute_sign19   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