Numerous fixes for IncompleteCholesky. Still have to make it fully exploit the sparse structure of the L factor, and improve robustness to illconditionned problems.

1 files changed, 103 insertions, 74 deletions

diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
index 35cfa315d..a62672201 100644
--- a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
+++ b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
@@ -26,25 +26,29 @@ namespace Eigen {
  * \tparam _OrderingType 
  */
 
-template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> >
+template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> >
 class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
 {
   protected:
     typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
     using Base::m_isInitialized;
   public:
-    typedef SparseMatrix<Scalar,ColMajor> MatrixType;
+    typedef typename NumTraits<Scalar>::Real RealScalar; 
     typedef _OrderingType OrderingType;
-    typedef typename MatrixType::RealScalar RealScalar; 
-    typedef typename MatrixType::Index Index; 
-    typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
-    typedef Matrix<Scalar,Dynamic,1> ScalarType; 
-    typedef Matrix<Index,Dynamic, 1> IndexType;
-    typedef std::vector<std::list<Index> > VectorList; 
+    typedef typename OrderingType::PermutationType PermutationType;
+    typedef typename PermutationType::StorageIndex StorageIndex; 
+    typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
+    typedef FactorType MatrixType;
+    typedef Matrix<Scalar,Dynamic,1> VectorSx;
+    typedef Matrix<RealScalar,Dynamic,1> VectorRx;
+    typedef Matrix<StorageIndex,Dynamic, 1> VectorIx;
+    typedef std::vector<std::list<StorageIndex> > VectorList; 
     enum { UpLo = _UpLo };
   public:
-    IncompleteCholesky() : m_shift(1),m_factorizationIsOk(false) {}
-    IncompleteCholesky(const MatrixType& matrix) : m_shift(1),m_factorizationIsOk(false)
+    IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {}
+    
+    template<typename MatrixType>
+    IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false)
     {
       compute(matrix);
     }
@@ -68,7 +72,7 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up
     /** 
      * \brief Set the initial shift parameter
      */
-    void setShift( Scalar shift) { m_shift = shift; }
+    void setShift( Scalar shift) { m_initialShift = shift; }
     
     /**
     * \brief Computes the fill reducing permutation vector. 
@@ -77,7 +81,10 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up
     void analyzePattern(const MatrixType& mat)
     {
       OrderingType ord; 
-      ord(mat.template selfadjointView<UpLo>(), m_perm); 
+      PermutationType pinv;
+      ord(mat.template selfadjointView<UpLo>(), pinv); 
+      if(pinv.size()>0) m_perm = pinv.inverse();
+      else              m_perm.resize(0);
       m_analysisIsOk = true; 
     }
     
@@ -85,7 +92,7 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up
     void factorize(const MatrixType& amat);
     
     template<typename MatrixType>
-    void compute (const MatrixType& matrix)
+    void compute(const MatrixType& matrix)
     {
       analyzePattern(matrix); 
       factorize(matrix);
@@ -95,30 +102,28 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up
     void _solve_impl(const Rhs& b, Dest& x) const
     {
       eigen_assert(m_factorizationIsOk && "factorize() should be called first");
+      if (m_perm.rows() == b.rows())  x = m_perm * b;
+      else                            x = b;
+      x = m_scale.asDiagonal() * x;
+      x = m_L.template triangularView<Lower>().solve(x);
+      x = m_L.adjoint().template triangularView<Upper>().solve(x);
+      x = m_scale.asDiagonal() * x;
       if (m_perm.rows() == b.rows())
-        x = m_perm.inverse() * b; 
-      else 
-        x = b; 
-      x = m_scal.asDiagonal() * x;
-      x = m_L.template triangularView<UnitLower>().solve(x); 
-      x = m_L.adjoint().template triangularView<Upper>().solve(x); 
-      if (m_perm.rows() == b.rows())
-        x = m_perm * x;
-      x = m_scal.asDiagonal() * x;
+        x = m_perm.inverse() * x;
+      
     }
 
   protected:
-    SparseMatrix<Scalar,ColMajor> m_L;  // The lower part stored in CSC
-    ScalarType m_scal; // The vector for scaling the matrix 
-    Scalar m_shift; //The initial shift parameter
+    FactorType m_L;             // The lower part stored in CSC
+    VectorRx m_scale;            // The vector for scaling the matrix 
+    Scalar m_initialShift;      // The initial shift parameter
     bool m_analysisIsOk; 
     bool m_factorizationIsOk; 
     ComputationInfo m_info;
     PermutationType m_perm; 
     
   private:
-    template <typename IdxType, typename SclType>
-    inline void updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol); 
+    inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol); 
 }; 
 
 template<typename Scalar, int _UpLo, typename OrderingType>
@@ -128,93 +133,118 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
   using std::sqrt;
   eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 
     
-  // Dropping strategies : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
+  // Dropping strategy : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
+  
+  m_L.resize(mat.rows(), mat.cols());
   
   // Apply the fill-reducing permutation computed in analyzePattern()
   if (m_perm.rows() == mat.rows() ) // To detect the null permutation
-    m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
+  {
+    // The temporary is needed to make sure that the diagonal entry is properly sorted
+    FactorType tmp(mat.rows(), mat.cols());
+    tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
+    m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
+  }
   else
+  {
     m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
+  }
   
   Index n = m_L.cols(); 
   Index nnz = m_L.nonZeros();
-  Map<ScalarType> vals(m_L.valuePtr(), nnz); //values
-  Map<IndexType> rowIdx(m_L.innerIndexPtr(), nnz);  //Row indices
-  Map<IndexType> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
-  IndexType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
-  VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
-  ScalarType curCol(n); // Store a  nonzero values in each column
-  IndexType irow(n); // Row indices of nonzero elements in each column
+  Map<VectorSx> vals(m_L.valuePtr(), nnz);         //values
+  Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz);  //Row indices
+  Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
+  VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
+  VectorList listCol(n);  // listCol(j) is a linked list of columns to update column j
+  VectorSx curCol(n);     // Store a  nonzero values in each column
+  VectorIx irow(n);       // Row indices of nonzero elements in each column
   
   
   // Computes the scaling factors 
-  m_scal.resize(n);
-  for (int j = 0; j < n; j++)
-  {
-    m_scal(j) = m_L.col(j).norm();
-    m_scal(j) = sqrt(m_scal(j));
-  }
+  m_scale.resize(n);
+  m_scale.setZero();
+  for (Index j = 0; j < n; j++)
+    for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
+    {
+      m_scale(j) += numext::abs2(vals(k));
+      if(rowIdx[k]!=j)
+        m_scale(rowIdx[k]) += numext::abs2(vals(k));
+    }
+  
+  m_scale = m_scale.cwiseSqrt().cwiseSqrt();
+  
   // Scale and compute the shift for the matrix 
-  Scalar mindiag = vals[0];
-  for (int j = 0; j < n; j++){
-    for (int k = colPtr[j]; k < colPtr[j+1]; k++)
-     vals[k] /= (m_scal(j) * m_scal(rowIdx[k]));
-    mindiag = numext::mini(vals[colPtr[j]], mindiag);
+  RealScalar mindiag = NumTraits<RealScalar>::highest();
+  for (Index j = 0; j < n; j++)
+  {
+    for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
+      vals[k] /= (m_scale(j)*m_scale(rowIdx[k]));
+    eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
+    mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
   }
   
-  if(mindiag < Scalar(0.)) m_shift = m_shift - mindiag;
+  RealScalar shift = 0;
+  if(mindiag <= RealScalar(0.))
+    shift = m_initialShift - mindiag;
+
   // Apply the shift to the diagonal elements of the matrix
-  for (int j = 0; j < n; j++)
-    vals[colPtr[j]] += m_shift;
+  for (Index j = 0; j < n; j++)
+    vals[colPtr[j]] += shift;
+  
   // jki version of the Cholesky factorization 
-  for (int j=0; j < n; ++j)
+  for (Index j=0; j < n; ++j)
   {  
-    //Left-looking factorize the column j 
-    // First, load the jth column into curCol 
+    // Left-looking factorization of the j-th column
+    // First, load the j-th column into curCol 
     Scalar diag = vals[colPtr[j]];  // It is assumed that only the lower part is stored
     curCol.setZero();
-    irow.setLinSpaced(n,0,n-1); 
-    for (int i = colPtr[j] + 1; i < colPtr[j+1]; i++)
+    irow.setLinSpaced(n,0,internal::convert_index<StorageIndex,Index>(n-1)); 
+    for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
     {
       curCol(rowIdx[i]) = vals[i]; 
       irow(rowIdx[i]) = rowIdx[i]; 
     }
-    std::list<int>::iterator k; 
+    typename std::list<StorageIndex>::iterator k; 
     // Browse all previous columns that will update column j
     for(k = listCol[j].begin(); k != listCol[j].end(); k++) 
     {
-      int jk = firstElt(*k); // First element to use in the column 
+      Index jk = firstElt(*k); // First element to use in the column 
+      eigen_internal_assert(rowIdx[jk]==j);
+      Scalar v_j_jk = numext::conj(vals[jk]);
+      
       jk += 1; 
-      for (int i = jk; i < colPtr[*k+1]; i++)
-      {
-        curCol(rowIdx[i]) -= vals[i] * vals[jk] ;
-      }
+      for (Index i = jk; i < colPtr[*k+1]; i++)
+        curCol(rowIdx[i]) -= vals[i] * v_j_jk;
       updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
     }
     
     // Scale the current column
-    if(RealScalar(diag) <= 0) 
+    if(numext::real(diag) <= 0) 
     {
-      std::cerr << "\nNegative diagonal during Incomplete factorization... "<< j << "\n";
+      //std::cerr << "\nNegative diagonal during Incomplete factorization at position " << j << " (value = " << diag << ")\n";
       m_info = NumericalIssue; 
       return; 
     }
-    RealScalar rdiag = sqrt(RealScalar(diag));
+    
+    RealScalar rdiag = sqrt(numext::real(diag));
     vals[colPtr[j]] = rdiag;
-    for (int i = j+1; i < n; i++)
+    // TODO, the following should iterate on the structurally non-zeros only
+    for (Index i = j+1; i < n; i++)
     {
       //Scale 
       curCol(i) /= rdiag;
       //Update the remaining diagonals with curCol
-      vals[colPtr[i]] -= curCol(i) * curCol(i);
+      vals[colPtr[i]] -= numext::abs2(curCol(i));
     }
     // Select the largest p elements
-    //  p is the original number of elements in the column (without the diagonal)
-    int p = colPtr[j+1] - colPtr[j] - 1 ; 
+    // p is the original number of elements in the column (without the diagonal)
+    // TODO, QuickSplit should operate on the structurally non zeros only.
+    Index p = colPtr[j+1] - colPtr[j] - 1 ; 
     internal::QuickSplit(curCol, irow, p); 
     // Insert the largest p elements in the matrix
-    int cpt = 0; 
-    for (int i = colPtr[j]+1; i < colPtr[j+1]; i++)
+    Index cpt = 0; 
+    for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
     {
       vals[i] = curCol(cpt); 
       rowIdx[i] = irow(cpt); 
@@ -230,8 +260,7 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
 }
 
 template<typename Scalar, int _UpLo, typename OrderingType>
-template <typename IdxType, typename SclType>
-inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol)
+inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol)
 {
   if (jk < colPtr(col+1) )
   {
@@ -245,8 +274,8 @@ inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const Idx
       std::swap(rowIdx(jk),rowIdx(minpos));
       std::swap(vals(jk),vals(minpos));
     }
-    firstElt(col) = jk;
-    listCol[rowIdx(jk)].push_back(col);
+    firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
+    listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
   }
 }