Addresses comments on Eigen pull request PR-174.

* Get rid of code-duplication for real vs. complex matrices.
* Fix flipped arguments to select.
* Make the condition estimation functions free functions.
* Use Vector::Unit() to generate canonical unit vectors.
* Misc. cleanup.

1 files changed, 152 insertions, 253 deletions

diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h
index dca3da417..12b4ae648 100644
--- a/Eigen/src/Core/ConditionEstimator.h
+++ b/Eigen/src/Core/ConditionEstimator.h
@@ -13,45 +13,35 @@
 namespace Eigen {
 
 namespace internal {
-template <typename Decomposition, bool IsSelfAdjoint, bool IsComplex>
-struct EstimateInverseMatrixL1NormImpl {};
-}  // namespace internal
-
-template <typename Decomposition, bool IsSelfAdjoint = false>
-class ConditionEstimator {
- public:
-  typedef typename Decomposition::MatrixType MatrixType;
-  typedef typename internal::traits<MatrixType>::Scalar Scalar;
-  typedef typename NumTraits<Scalar>::Real RealScalar;
-  typedef typename internal::plain_col_type<MatrixType>::type Vector;
+template <typename MatrixType>
+inline typename MatrixType::RealScalar MatrixL1Norm(const MatrixType& matrix) {
+  return matrix.cwiseAbs().colwise().sum().maxCoeff();
+}
+
+template <typename Vector>
+inline typename Vector::RealScalar VectorL1Norm(const Vector& v) {
+  return v.template lpNorm<1>();
+}
+
+template <typename Vector, typename RealVector, bool IsComplex>
+struct SignOrUnity {
+  static inline Vector run(const Vector& v) {
+    const RealVector v_abs = v.cwiseAbs();
+    return (v_abs.array() == 0).select(Vector::Ones(v.size()), v.cwiseQuotient(v_abs));
+  }
+};
 
-  /** \class ConditionEstimator
-    * \ingroup Core_Module
-    *
-    * \brief Condition number estimator.
-    *
-    * Computing a decomposition of a dense matrix takes O(n^3) operations, while
-    * this method estimates the condition number quickly and reliably in O(n^2)
-    * operations.
-    *
-    * \returns an estimate of the reciprocal condition number
-    * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given the matrix and
-    * its decomposition. Supports the following decompositions: FullPivLU,
-    * PartialPivLU.
-    *
-    * \sa FullPivLU, PartialPivLU.
-    */
-  static RealScalar rcond(const MatrixType& matrix, const Decomposition& dec) {
-    eigen_assert(matrix.rows() == dec.rows());
-    eigen_assert(matrix.cols() == dec.cols());
-    eigen_assert(matrix.rows() == matrix.cols());
-    if (dec.rows() == 0) {
-      return RealScalar(1);
-    }
-    return rcond(MatrixL1Norm(matrix), dec);
+// Partial specialization to avoid elementwise division for real vectors.
+template <typename Vector>
+struct SignOrUnity<Vector, Vector, false> {
+  static inline Vector run(const Vector& v) {
+    return (v.array() < 0).select(-Vector::Ones(v.size()), Vector::Ones(v.size()));
   }
+};
 
-  /** \class ConditionEstimator
+}  // namespace internal
+
+/** \class ConditionEstimator
     * \ingroup Core_Module
     *
     * \brief Condition number estimator.
@@ -61,245 +51,154 @@ class ConditionEstimator {
     * operations.
     *
     * \returns an estimate of the reciprocal condition number
-    * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and
+    * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given the matrix and
     * its decomposition. Supports the following decompositions: FullPivLU,
     * PartialPivLU.
     *
     * \sa FullPivLU, PartialPivLU.
     */
-  static RealScalar rcond(RealScalar matrix_norm, const Decomposition& dec) {
-    eigen_assert(dec.rows() == dec.cols());
-    if (dec.rows() == 0) {
-      return 1;
-    }
-    if (matrix_norm == 0) {
-      return 0;
-    }
-    const RealScalar inverse_matrix_norm = EstimateInverseMatrixL1Norm(dec);
-    return inverse_matrix_norm == 0 ? 0
-                                    : (1 / inverse_matrix_norm) / matrix_norm;
+template <typename Decomposition>
+typename Decomposition::RealScalar ReciprocalConditionNumberEstimate(
+    const typename Decomposition::MatrixType& matrix,
+    const Decomposition& dec) {
+  eigen_assert(matrix.rows() == dec.rows());
+  eigen_assert(matrix.cols() == dec.cols());
+  eigen_assert(matrix.rows() == matrix.cols());
+  if (dec.rows() == 0) {
+    return Decomposition::RealScalar(1);
   }
-
-  /**
-   * \returns an estimate of ||inv(matrix)||_1 given a decomposition of
-   * matrix that implements .solve() and .adjoint().solve() methods.
-   *
-   * The method implements Algorithms 4.1 and 5.1 from
-   *   http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf
-   * which also forms the basis for the condition number estimators in
-   * LAPACK. Since at most 10 calls to the solve method of dec are
-   * performed, the total cost is O(dims^2), as opposed to O(dims^3)
-   * needed to compute the inverse matrix explicitly.
-   *
-   * The most common usage is in estimating the condition number
-   * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be
-   * computed directly in O(n^2) operations.
-   */
-  static RealScalar EstimateInverseMatrixL1Norm(const Decomposition& dec) {
-    eigen_assert(dec.rows() == dec.cols());
-    if (dec.rows() == 0) {
-      return 0;
-    }
-    return internal::EstimateInverseMatrixL1NormImpl<
-        Decomposition, IsSelfAdjoint,
-        NumTraits<Scalar>::IsComplex != 0>::compute(dec);
+  return ReciprocalConditionNumberEstimate(MatrixL1Norm(matrix), dec);
+}
+
+/** \class ConditionEstimator
+ * \ingroup Core_Module
+ *
+ * \brief Condition number estimator.
+ *
+ * Computing a decomposition of a dense matrix takes O(n^3) operations, while
+ * this method estimates the condition number quickly and reliably in O(n^2)
+ * operations.
+ *
+ * \returns an estimate of the reciprocal condition number
+ * (1 / (||matrix||_1 * ||inv(matrix)||_1)) of matrix, given ||matrix||_1 and
+ * its decomposition. Supports the following decompositions: FullPivLU,
+ * PartialPivLU.
+ *
+ * \sa FullPivLU, PartialPivLU.
+ */
+template <typename Decomposition>
+typename Decomposition::RealScalar ReciprocalConditionNumberEstimate(
+    typename Decomposition::RealScalar matrix_norm, const Decomposition& dec) {
+  eigen_assert(dec.rows() == dec.cols());
+  if (dec.rows() == 0) {
+    return 1;
   }
-
-  /**
-   * \returns the induced matrix l1-norm
-   * ||matrix||_1 = sup ||matrix * v||_1 / ||v||_1, which is equal to
-   * the greatest absolute column sum.
-   */
-  static inline Scalar MatrixL1Norm(const MatrixType& matrix) {
-    return matrix.cwiseAbs().colwise().sum().maxCoeff();
+  if (matrix_norm == 0) {
+    return 0;
   }
-};
-
-namespace internal {
-
-template <typename Decomposition, typename Vector, bool IsSelfAdjoint = false>
-struct solve_helper {
-  static inline Vector solve_adjoint(const Decomposition& dec,
-                                     const Vector& v) {
-    return dec.adjoint().solve(v);
-  }
-};
-
-// Partial specialization for self_adjoint matrices.
-template <typename Decomposition, typename Vector>
-struct solve_helper<Decomposition, Vector, true> {
-  static inline Vector solve_adjoint(const Decomposition& dec,
-                                     const Vector& v) {
-    return dec.solve(v);
-  }
-};
-
-
-// Partial specialization for real matrices.
-template <typename Decomposition, bool IsSelfAdjoint>
-struct EstimateInverseMatrixL1NormImpl<Decomposition, IsSelfAdjoint, false> {
+  const typename Decomposition::RealScalar inverse_matrix_norm = InverseMatrixL1NormEstimate(dec);
+  return inverse_matrix_norm == 0 ? 0 : (1 / inverse_matrix_norm) / matrix_norm;
+}
+
+/**
+ * \returns an estimate of ||inv(matrix)||_1 given a decomposition of
+ * matrix that implements .solve() and .adjoint().solve() methods.
+ *
+ * The method implements Algorithms 4.1 and 5.1 from
+ *   http://www.maths.manchester.ac.uk/~higham/narep/narep135.pdf
+ * which also forms the basis for the condition number estimators in
+ * LAPACK. Since at most 10 calls to the solve method of dec are
+ * performed, the total cost is O(dims^2), as opposed to O(dims^3)
+ * needed to compute the inverse matrix explicitly.
+ *
+ * The most common usage is in estimating the condition number
+ * ||matrix||_1 * ||inv(matrix)||_1. The first term ||matrix||_1 can be
+ * computed directly in O(n^2) operations.
+ */
+template <typename Decomposition>
+typename Decomposition::RealScalar InverseMatrixL1NormEstimate(
+    const Decomposition& dec) {
   typedef typename Decomposition::MatrixType MatrixType;
-  typedef typename internal::traits<MatrixType>::Scalar Scalar;
+  typedef typename Decomposition::Scalar Scalar;
+  typedef typename Decomposition::RealScalar RealScalar;
   typedef typename internal::plain_col_type<MatrixType>::type Vector;
+  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVector;
+  const bool is_complex = (NumTraits<Scalar>::IsComplex != 0);
 
-  // Shorthand for vector L1 norm in Eigen.
-  static inline Scalar VectorL1Norm(const Vector& v) {
-    return v.template lpNorm<1>();
+  eigen_assert(dec.rows() == dec.cols());
+  const int n = dec.rows();
+  if (n == 0) {
+    return 0;
   }
-
-  static inline Scalar compute(const Decomposition& dec) {
-    const int n = dec.rows();
-    const Vector plus = Vector::Ones(n);
-    Vector v = plus / n;
-    v = dec.solve(v);
-    Scalar lower_bound = VectorL1Norm(v);
-    if (n == 1) {
-      return lower_bound;
-    }
-    // lower_bound is a lower bound on
-    //   ||inv(matrix)||_1  = sup_v ||inv(matrix) v||_1 / ||v||_1
-    // and is the objective maximized by the ("super-") gradient ascent
-    // algorithm.
-    // Basic idea: We know that the optimum is achieved at one of the simplices
-    // v = e_i, so in each iteration we follow a super-gradient to move towards
-    // the optimal one.
-    Scalar old_lower_bound = lower_bound;
-    const Vector minus = -Vector::Ones(n);
-    Vector sign_vector = (v.cwiseAbs().array() == 0).select(plus, minus);
-    Vector old_sign_vector = sign_vector;
-    int v_max_abs_index = -1;
-    int old_v_max_abs_index = v_max_abs_index;
-    for (int k = 0; k < 4; ++k) {
-      // argmax |inv(matrix)^T * sign_vector|
-      v = solve_helper<Decomposition, Vector, IsSelfAdjoint>::solve_adjoint(dec, sign_vector);
-      v.cwiseAbs().maxCoeff(&v_max_abs_index);
-      if (v_max_abs_index == old_v_max_abs_index) {
-        // Break if the solution stagnated.
-        break;
-      }
-      // Move to the new simplex e_j, where j = v_max_abs_index.
-      v.setZero();
-      v[v_max_abs_index] = 1;
-      v = dec.solve(v);  // v = inv(matrix) * e_j.
-      lower_bound = VectorL1Norm(v);
-      if (lower_bound <= old_lower_bound) {
-        // Break if the gradient step did not increase the lower_bound.
-        break;
-      }
-      sign_vector = (v.array() < 0).select(plus, minus);
+  Vector v = Vector::Ones(n) / n;
+  v = dec.solve(v);
+
+  // lower_bound is a lower bound on
+  //   ||inv(matrix)||_1  = sup_v ||inv(matrix) v||_1 / ||v||_1
+  // and is the objective maximized by the ("super-") gradient ascent
+  // algorithm below.
+  RealScalar lower_bound = internal::VectorL1Norm(v);
+  if (n == 1) {
+    return lower_bound;
+  }
+  // Gradient ascent algorithm follows: We know that the optimum is achieved at
+  // one of the simplices v = e_i, so in each iteration we follow a
+  // super-gradient to move towards the optimal one.
+  RealScalar old_lower_bound = lower_bound;
+  Vector sign_vector(n);
+  Vector old_sign_vector;
+  int v_max_abs_index = -1;
+  int old_v_max_abs_index = v_max_abs_index;
+  for (int k = 0; k < 4; ++k) {
+    sign_vector = internal::SignOrUnity<Vector, RealVector, is_complex>::run(v);
+    if (k > 0 && !is_complex) {
       if (sign_vector == old_sign_vector) {
         // Break if the solution stagnated.
         break;
       }
-      old_sign_vector = sign_vector;
-      old_v_max_abs_index = v_max_abs_index;
-      old_lower_bound = lower_bound;
-    }
-    // The following calculates an independent estimate of ||matrix||_1 by
-    // multiplying matrix by a vector with entries of slowly increasing
-    // magnitude and alternating sign:
-    //   v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1.
-    // This improvement to Hager's algorithm above is due to Higham. It was
-    // added to make the algorithm more robust in certain corner cases where
-    // large elements in the matrix might otherwise escape detection due to
-    // exact cancellation (especially when op and op_adjoint correspond to a
-    // sequence of backsubstitutions and permutations), which could cause
-    // Hager's algorithm to vastly underestimate ||matrix||_1.
-    Scalar alternating_sign = 1;
-    for (int i = 0; i < n; ++i) {
-      v[i] = alternating_sign * static_cast<Scalar>(1) +
-             (static_cast<Scalar>(i) / (static_cast<Scalar>(n - 1)));
-      alternating_sign = -alternating_sign;
     }
-    v = dec.solve(v);
-    const Scalar alternate_lower_bound =
-        (2 * VectorL1Norm(v)) / (3 * static_cast<Scalar>(n));
-    return numext::maxi(lower_bound, alternate_lower_bound);
-  }
-};
-
-// Partial specialization for complex matrices.
-template <typename Decomposition, bool IsSelfAdjoint>
-struct EstimateInverseMatrixL1NormImpl<Decomposition, IsSelfAdjoint, true> {
-  typedef typename Decomposition::MatrixType MatrixType;
-  typedef typename internal::traits<MatrixType>::Scalar Scalar;
-  typedef typename NumTraits<Scalar>::Real RealScalar;
-  typedef typename internal::plain_col_type<MatrixType>::type Vector;
-  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type
-      RealVector;
-
-  // Shorthand for vector L1 norm in Eigen.
-  static inline RealScalar VectorL1Norm(const Vector& v) {
-    return v.template lpNorm<1>();
-  }
-
-  static inline RealScalar compute(const Decomposition& dec) {
-    const int n = dec.rows();
-    const Vector ones = Vector::Ones(n);
-    Vector v = ones / n;
-    v = dec.solve(v);
-    RealScalar lower_bound = VectorL1Norm(v);
-    if (n == 1) {
-      return lower_bound;
+    // v_max_abs_index = argmax |real( inv(matrix)^T * sign_vector )|
+    v = dec.adjoint().solve(sign_vector);
+    v.real().cwiseAbs().maxCoeff(&v_max_abs_index);
+    if (v_max_abs_index == old_v_max_abs_index) {
+      // Break if the solution stagnated.
+      break;
     }
-    // lower_bound is a lower bound on
-    //   ||inv(matrix)||_1  = sup_v ||inv(matrix) v||_1 / ||v||_1
-    // and is the objective maximized by the ("super-") gradient ascent
-    // algorithm.
-    // Basic idea: We know that the optimum is achieved at one of the simplices
-    // v = e_i, so in each iteration we follow a super-gradient to move towards
-    // the optimal one.
-    RealScalar old_lower_bound = lower_bound;
-    int v_max_abs_index = -1;
-    int old_v_max_abs_index = v_max_abs_index;
-    for (int k = 0; k < 4; ++k) {
-      // argmax |inv(matrix)^* * sign_vector|
-      RealVector abs_v = v.cwiseAbs();
-      const Vector psi =
-          (abs_v.array() == 0).select(v.cwiseQuotient(abs_v), ones);
-      v = solve_helper<Decomposition, Vector, IsSelfAdjoint>::solve_adjoint(dec, psi);
-      const RealVector z = v.real();
-      z.cwiseAbs().maxCoeff(&v_max_abs_index);
-      if (v_max_abs_index == old_v_max_abs_index) {
-        // Break if the solution stagnated.
-        break;
-      }
-      // Move to the new simplex e_j, where j = v_max_abs_index.
-      v.setZero();
-      v[v_max_abs_index] = 1;
-      v = dec.solve(v);  // v = inv(matrix) * e_j.
-      lower_bound = VectorL1Norm(v);
-      if (lower_bound <= old_lower_bound) {
-        // Break if the gradient step did not increase the lower_bound.
-        break;
-      }
-      old_v_max_abs_index = v_max_abs_index;
-      old_lower_bound = lower_bound;
+    // Move to the new simplex e_j, where j = v_max_abs_index.
+    v = dec.solve(Vector::Unit(n, v_max_abs_index));  // v = inv(matrix) * e_j.
+    lower_bound = internal::VectorL1Norm(v);
+    if (lower_bound <= old_lower_bound) {
+      // Break if the gradient step did not increase the lower_bound.
+      break;
     }
-    // The following calculates an independent estimate of ||matrix||_1 by
-    // multiplying matrix by a vector with entries of slowly increasing
-    // magnitude and alternating sign:
-    //   v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1.
-    // This improvement to Hager's algorithm above is due to Higham. It was
-    // added to make the algorithm more robust in certain corner cases where
-    // large elements in the matrix might otherwise escape detection due to
-    // exact cancellation (especially when op and op_adjoint correspond to a
-    // sequence of backsubstitutions and permutations), which could cause
-    // Hager's algorithm to vastly underestimate ||matrix||_1.
-    RealScalar alternating_sign = 1;
-    for (int i = 0; i < n; ++i) {
-      v[i] = alternating_sign * static_cast<RealScalar>(1) +
-             (static_cast<RealScalar>(i) / (static_cast<RealScalar>(n - 1)));
-      alternating_sign = -alternating_sign;
+    if (!is_complex) {
+      old_sign_vector = sign_vector;
     }
-    v = dec.solve(v);
-    const RealScalar alternate_lower_bound =
-        (2 * VectorL1Norm(v)) / (3 * static_cast<RealScalar>(n));
-    return numext::maxi(lower_bound, alternate_lower_bound);
+    old_v_max_abs_index = v_max_abs_index;
+    old_lower_bound = lower_bound;
   }
-};
+  // The following calculates an independent estimate of ||matrix||_1 by
+  // multiplying matrix by a vector with entries of slowly increasing
+  // magnitude and alternating sign:
+  //   v_i = (-1)^{i} (1 + (i / (dim-1))), i = 0,...,dim-1.
+  // This improvement to Hager's algorithm above is due to Higham. It was
+  // added to make the algorithm more robust in certain corner cases where
+  // large elements in the matrix might otherwise escape detection due to
+  // exact cancellation (especially when op and op_adjoint correspond to a
+  // sequence of backsubstitutions and permutations), which could cause
+  // Hager's algorithm to vastly underestimate ||matrix||_1.
+  Scalar alternating_sign = 1;
+  for (int i = 0; i < n; ++i) {
+    v[i] = alternating_sign * static_cast<RealScalar>(1) +
+           (static_cast<RealScalar>(i) / (static_cast<RealScalar>(n - 1)));
+    alternating_sign = -alternating_sign;
+  }
+  v = dec.solve(v);
+  const RealScalar alternate_lower_bound =
+      (2 * internal::VectorL1Norm(v)) / (3 * static_cast<RealScalar>(n));
+  return numext::maxi(lower_bound, alternate_lower_bound);
+}
 
-}  // namespace internal
 }  // namespace Eigen
 
 #endif