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
13namespace Eigen {
14
15namespace internal {
16
17template <typename Vector, typename RealVector, bool IsComplex>
18struct 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.
27template <typename Vector>
28struct 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  */
55template <typename Decomposition>
56typename 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  */
157template <typename Decomposition>
158typename Decomposition::RealScalar
159rcond_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