backporting LLT accuracy fixes

Eigen/src/Cholesky/LLT.h

@@ -1,5 +1,5 @@
 // This file is part of Eigen, a lightweight C++ template library
-// for linear algebra. Eigen itself is part of the KDE project.
+// for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
@@ -41,11 +41,16 @@
   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
   * situations like generalised eigen problems with hermitian matrices.
   *
-  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+  * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
-  * the strict lower part does not have to store correct values.
+  * 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)
+  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+  * the strict lower part does not have to store correct values.
+  */
 template<typename MatrixType> class LLT
 {
   private:
@@ -60,17 +65,30 @@ template<typename MatrixType> class LLT
 
   public:
 
+    /** 
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via LLT::compute(const MatrixType&).
+    */
+    LLT() : m_matrix(), m_isInitialized(false) {}
+
     LLT(const MatrixType& matrix)
-      : m_matrix(matrix.rows(), matrix.cols())
+      : m_matrix(matrix.rows(), matrix.cols()),
+        m_isInitialized(false)
     {
       compute(matrix);
     }
 
     /** \returns the lower triangular matrix L */
-    inline Part<MatrixType, LowerTriangular> matrixL(void) const { return m_matrix; }
+    inline Part<MatrixType, LowerTriangular> matrixL(void) const 
-
+    { 
-    /** \returns true if the matrix is positive definite */
+      ei_assert(m_isInitialized && "LLT is not initialized.");
-    inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
+      return m_matrix; 
+    }
+    
+    /** \deprecated */
+    inline bool isPositiveDefinite(void) const { return m_isInitialized && m_isPositiveDefinite; }
 
     template<typename RhsDerived, typename ResDerived>
     bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
@@ -86,6 +104,7 @@ template<typename MatrixType> class LLT
       * The strict upper part is not used and even not initialized.
       */
     MatrixType m_matrix;
+    bool m_isInitialized;
     bool m_isPositiveDefinite;
 };
 
@@ -95,24 +114,34 @@ template<typename MatrixType>
 void LLT<MatrixType>::compute(const MatrixType& a)
 {
   assert(a.rows()==a.cols());
+  m_isPositiveDefinite = true;
   const int size = a.rows();
   m_matrix.resize(size, size);
-  const RealScalar eps = ei_sqrt(precision<Scalar>());
+  // 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 = x > eps && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1));
   m_matrix.coeffRef(0,0) = ei_sqrt(x);
+  if(size==1)
+  {
+    m_isInitialized = true;
+    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 (x < cutoff)
-    if (x < eps || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1))))
     {
       m_isPositiveDefinite = false;
-      return;
+      continue;
     }
+
     m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
 
     int endSize = size-j-1;
@@ -127,12 +156,14 @@ void LLT<MatrixType>::compute(const MatrixType& a)
         - m_matrix.col(j).end(endSize) ) / x;
     }
   }
+
+  m_isInitialized = true;
 }
 
 /** 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.
@@ -146,6 +177,7 @@ template<typename MatrixType>
 template<typename RhsDerived, typename ResDerived>
 bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
 {
+  ei_assert(m_isInitialized && "LLT is not initialized.");
   const int size = m_matrix.rows();
   ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
   return solveInPlace((*result) = b);
@@ -155,6 +187,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.
   *
@@ -164,10 +198,9 @@ template<typename MatrixType>
 template<typename Derived>
 bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
 {
+  ei_assert(m_isInitialized && "LLT is not initialized.");
   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()
@@ -49,6 +50,7 @@ template<typename T> inline T 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; }
@@ -93,6 +95,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; }
@@ -138,6 +141,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; }
@@ -183,6 +188,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); }
@@ -216,6 +222,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); }
@@ -250,6 +257,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; }