More Cholesky fixes.

* Cholesky decs are NOT rank revealing so remove all the rank/isPositiveDefinite etc stuff.
* fix bug in LLT: s/return/continue/
* introduce machine_epsilon constants, they are actually needed for Higman's formula determining
  the cutoff in Cholesky. Btw fix the page reference to his book (chat with Keir).
* solve methods always return true, since this isn't a rank revealing dec. Actually... they already did always return true!! Now it's explicit.
* updated dox and unit-test

Eigen/src/Cholesky/LDLT.h

@@ -43,6 +43,9 @@
   * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root
   * on D also stabilizes the computation.
   *
+  * Remember that Cholesky decompositions are not rank-revealing.  Also, do not use a Cholesky decomposition to determine
+  * whether a system of equations has a solution.
+  *
   * \sa MatrixBase::ldlt(), class LLT
   */
  /* THIS PART OF THE DOX IS CURRENTLY DISABLED BECAUSE INACCURATE BECAUSE OF BUG IN THE DECOMPOSITION CODE
@@ -88,25 +91,6 @@ template<typename MatrixType> class LDLT
     /** \returns true if the matrix is negative (semidefinite) */
     inline bool isNegative(void) const { return m_sign == -1; }
 
-    /** \returns true if the matrix is invertible */
-    inline bool isInvertible(void) const { return m_rank == m_matrix.rows(); }
-
-    /** \returns true if the matrix is positive definite */
-    inline bool isPositiveDefinite(void) const { return isPositive() && isInvertible(); }
-
-    /** \returns true if the matrix is negative definite */
-    inline bool isNegativeDefinite(void) const { return isNegative() && isInvertible(); }
-
-    /** \returns the rank of the matrix of which *this is the LDLT decomposition.
-      *
-      * \note This is computed at the time of the construction of the LDLT decomposition. This
-      *       method does not perform any further computation.
-      */
-    inline int rank() const
-    {
-      return m_rank;
-    }
-
     template<typename RhsDerived, typename ResDerived>
     bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
 
@@ -125,7 +109,7 @@ template<typename MatrixType> class LDLT
     MatrixType m_matrix;
     IntColVectorType m_p;
     IntColVectorType m_transpositions;
-    int m_rank, m_sign;
+    int m_sign;
 };
 
 /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
@@ -135,7 +119,6 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
 {
   ei_assert(a.rows()==a.cols());
   const int size = a.rows();
-  m_rank = size;
 
   m_matrix = a;
 
@@ -168,8 +151,8 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
       // to the largest overall, the algorithm bails.  This cutoff is suggested
       // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
       // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
-      // Algorithms" page 208, also by Higham.
+      // Algorithms" page 217, also by Higham.
-      cutoff = ei_abs(precision<RealScalar>() * size * biggest_in_corner);
+      cutoff = ei_abs(machine_epsilon<Scalar>() * size * biggest_in_corner);
 
       m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
     }
@@ -178,7 +161,6 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
     if(biggest_in_corner < cutoff)
     {
       for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i;
-      m_rank = j;
       break;
     }
 
@@ -200,11 +182,9 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
     m_matrix.coeffRef(j,j) = Djj;
 
     // Finish early if the matrix is not full rank.
-    if(ei_abs(Djj) < cutoff) // i made experiments, this is better than isMuchSmallerThan(biggest_in_corner), and of course
+    if(ei_abs(Djj) < cutoff)
-                             // much better than plain sign comparison as used to be done before.
     {
       for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i;
-      m_rank = j;
       break;
     }
 
@@ -230,7 +210,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
 /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
   * The result is stored in \a result
   *
-  * \returns true in case of success, false otherwise.
+  * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
   *
   * In other words, it computes \f$ b = A^{-1} b \f$ with
   * \f$ P^T{L^{*}}^{-1} D^{-1} L^{-1} P b \f$ from right to left.
@@ -252,6 +232,8 @@ bool LDLT<MatrixType>
   *
   * \param bAndX represents both the right-hand side matrix b and result x.
   *
+  * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
+  *
   * This version avoids a copy when the right hand side matrix b is not
   * needed anymore.
   *
@@ -264,8 +246,6 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
   const int size = m_matrix.rows();
   ei_assert(size == bAndX.rows());
 
-  if (m_rank != size) return false;
-
   // z = P b
   for(int i = 0; i < size; ++i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i));
 

Eigen/src/Cholesky/LLT.h

@@ -41,6 +41,10 @@
   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
   * situations like generalised eigen problems with hermitian matrices.
   *
+  * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
+  * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations
+  * has a solution.
+  *
   * \sa MatrixBase::llt(), class LDLT
   */
  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
@@ -70,12 +74,6 @@ template<typename MatrixType> class LLT
     /** \returns the lower triangular matrix L */
     inline Part<MatrixType, LowerTriangular> matrixL(void) const { return m_matrix; }
 
-    /** \returns true if the matrix is positive definite */
-    inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
-
-    /** \returns true if the matrix is invertible, which in this context is equivalent to positive definite */
-    inline bool isInvertible(void) const { return m_isPositiveDefinite; }
-
     template<typename RhsDerived, typename ResDerived>
     bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
 
@@ -90,7 +88,6 @@ template<typename MatrixType> class LLT
       * The strict upper part is not used and even not initialized.
       */
     MatrixType m_matrix;
-    bool m_isPositiveDefinite;
 };
 
 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
@@ -101,23 +98,24 @@ void LLT<MatrixType>::compute(const MatrixType& a)
   assert(a.rows()==a.cols());
   const int size = a.rows();
   m_matrix.resize(size, size);
-  const RealScalar reference = size * a.diagonal().cwise().abs().maxCoeff();
+  // The biggest overall is the point of reference to which further diagonals
+  // are compared; if any diagonal is negligible compared
+  // to the largest overall, the algorithm bails.  This cutoff is suggested
+  // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
+  // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
+  // Algorithms" page 217, also by Higham.
+  const RealScalar cutoff = machine_epsilon<Scalar>() * size * a.diagonal().cwise().abs().maxCoeff();
   RealScalar x;
   x = ei_real(a.coeff(0,0));
-  m_isPositiveDefinite = !ei_isMuchSmallerThan(x, reference) && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), reference);
   m_matrix.coeffRef(0,0) = ei_sqrt(x);
   if(size==1)
     return;
   m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
   for (int j = 1; j < size; ++j)
   {
-    Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
+    x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
-    x = ei_real(tmp);
+    if (ei_abs(x) < cutoff) continue;
-    if (ei_isMuchSmallerThan(x, reference) || (!ei_isMuchSmallerThan(ei_imag(tmp), reference)))
+
-    {
-      m_isPositiveDefinite = false;
-      return;
-    }
     m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
 
     int endSize = size-j-1;
@@ -137,7 +135,7 @@ void LLT<MatrixType>::compute(const MatrixType& a)
 /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
   * The result is stored in \a result
   *
-  * \returns true in case of success, false otherwise.
+  * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
   *
   * In other words, it computes \f$ b = A^{-1} b \f$ with
   * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
@@ -160,6 +158,8 @@ bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDeriv
   *
   * \param bAndX represents both the right-hand side matrix b and result x.
   *
+  * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
+  *
   * This version avoids a copy when the right hand side matrix b is not
   * needed anymore.
   *
@@ -171,8 +171,6 @@ bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
 {
   const int size = m_matrix.rows();
   ei_assert(size==bAndX.rows());
-  if (!m_isPositiveDefinite)
-    return false;
   matrixL().solveTriangularInPlace(bAndX);
   m_matrix.adjoint().template part<UpperTriangular>().solveTriangularInPlace(bAndX);
   return true;

Eigen/src/Core/MathFunctions.h

@@ -26,6 +26,7 @@
 #define EIGEN_MATHFUNCTIONS_H
 
 template<typename T> inline typename NumTraits<T>::Real precision();
+template<typename T> inline typename NumTraits<T>::Real machine_epsilon();
 template<typename T> inline T ei_random(T a, T b);
 template<typename T> inline T ei_random();
 template<typename T> inline T ei_random_amplitude()
@@ -50,6 +51,7 @@ template<typename T> inline typename NumTraits<T>::Real ei_hypot(T x, T y)
 **************/
 
 template<> inline int precision<int>() { return 0; }
+template<> inline int machine_epsilon<int>() { return 0; }
 inline int ei_real(int x)  { return x; }
 inline int ei_imag(int)    { return 0; }
 inline int ei_conj(int x)  { return x; }
@@ -102,6 +104,7 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = precision<int>())
 **************/
 
 template<> inline float precision<float>() { return 1e-5f; }
+template<> inline float machine_epsilon<float>() { return 1.192e-07f; }
 inline float ei_real(float x)  { return x; }
 inline float ei_imag(float)    { return 0.f; }
 inline float ei_conj(float x)  { return x; }
@@ -147,6 +150,8 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float
 **************/
 
 template<> inline double precision<double>() { return 1e-11; }
+template<> inline double machine_epsilon<double>() { return 2.220e-16; }
+
 inline double ei_real(double x)  { return x; }
 inline double ei_imag(double)    { return 0.; }
 inline double ei_conj(double x)  { return x; }
@@ -192,6 +197,7 @@ inline bool ei_isApproxOrLessThan(double a, double b, double prec = precision<do
 *********************/
 
 template<> inline float precision<std::complex<float> >() { return precision<float>(); }
+template<> inline float machine_epsilon<std::complex<float> >() { return machine_epsilon<float>(); }
 inline float ei_real(const std::complex<float>& x) { return std::real(x); }
 inline float ei_imag(const std::complex<float>& x) { return std::imag(x); }
 inline std::complex<float> ei_conj(const std::complex<float>& x) { return std::conj(x); }
@@ -225,6 +231,7 @@ inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>&
 **********************/
 
 template<> inline double precision<std::complex<double> >() { return precision<double>(); }
+template<> inline double machine_epsilon<std::complex<double> >() { return machine_epsilon<double>(); }
 inline double ei_real(const std::complex<double>& x) { return std::real(x); }
 inline double ei_imag(const std::complex<double>& x) { return std::imag(x); }
 inline std::complex<double> ei_conj(const std::complex<double>& x) { return std::conj(x); }
@@ -259,6 +266,7 @@ inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double
 ******************/
 
 template<> inline long double precision<long double>() { return precision<double>(); }
+template<> inline long double machine_epsilon<long double>() { return 1.084e-19l; }
 inline long double ei_real(long double x)  { return x; }
 inline long double ei_imag(long double)    { return 0.; }
 inline long double ei_conj(long double x)  { return x; }

test/cholesky.cpp

@@ -86,7 +86,6 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
 
   {
     LLT<SquareMatrixType> chol(symm);
-    VERIFY(chol.isPositiveDefinite());
     VERIFY_IS_APPROX(symm, chol.matrixL() * chol.matrixL().adjoint());
     chol.solve(vecB, &vecX);
     VERIFY_IS_APPROX(symm * vecX, vecB);
@@ -103,18 +102,6 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
 
   {
     LDLT<SquareMatrixType> ldlt(symm);
-    VERIFY(ldlt.isInvertible());
-    if(sign == 1)
-    {
-      VERIFY(ldlt.isPositive());
-      VERIFY(ldlt.isPositiveDefinite());
-    }
-    if(sign == -1)
-    {
-      VERIFY(ldlt.isNegative());
-      VERIFY(ldlt.isNegativeDefinite());
-    }
-
     // TODO(keir): This doesn't make sense now that LDLT pivots.
     //VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint());
     ldlt.solve(vecB, &vecX);
@@ -123,15 +110,6 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     VERIFY_IS_APPROX(symm * matX, matB);
   }
 
-  // test isPositiveDefinite on non definite matrix
-  if (rows>4)
-  {
-    SquareMatrixType symm =  a0.block(0,0,rows,cols-4) * a0.block(0,0,rows,cols-4).adjoint();
-    LLT<SquareMatrixType> chol(symm);
-    VERIFY(!chol.isPositiveDefinite());
-    LDLT<SquareMatrixType> cholnosqrt(symm);
-    VERIFY(!cholnosqrt.isPositiveDefinite());
-  }
 }
 
 template<typename Derived>
@@ -156,29 +134,6 @@ void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m)
   }
 }
 
-template<typename MatrixType> void ldlt_rank()
-{
-  // NOTE there seems to be a problem with too small sizes -- could easily lie in the doSomeRankPreservingOperations function
-  int rows = ei_random<int>(50,200);
-  int rank = ei_random<int>(40, rows-1);
-
-
-  // generate a random positive matrix a of given rank
-  MatrixType m = MatrixType::Random(rows,rows);
-  QR<MatrixType> qr(m);
-  typedef typename MatrixType::Scalar Scalar;
-  typedef typename NumTraits<Scalar>::Real RealScalar;
-  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> DiagVectorType;
-  DiagVectorType d(rows);
-  d.setZero();
-  for(int i = 0; i < rank; i++) d(i)=RealScalar(1);
-  MatrixType a = qr.matrixQ() * d.asDiagonal() * qr.matrixQ().adjoint();
-
-  LDLT<MatrixType> ldlt(a);
-
-  VERIFY( ei_abs(ldlt.rank() - rank) <= rank / 20 ); // yes, LDLT::rank is a bit inaccurate...
-}
-
 
 void test_cholesky()
 {
@@ -191,9 +146,4 @@ void test_cholesky()
     CALL_SUBTEST( cholesky(MatrixXd(17,17)) );
     CALL_SUBTEST( cholesky(MatrixXf(200,200)) );
   }
-  for(int i = 0; i < g_repeat/3; i++) {
-    CALL_SUBTEST( ldlt_rank<MatrixXd>() );
-    CALL_SUBTEST( ldlt_rank<MatrixXf>() );
-    CALL_SUBTEST( ldlt_rank<MatrixXcd>() );
-  }
 }