IncompleteCholesky.h 14.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_INCOMPLETE_CHOlESKY_H
#define EIGEN_INCOMPLETE_CHOlESKY_H

#include <vector>
#include <list>

namespace Eigen {  
/** 
  * \brief Modified Incomplete Cholesky with dual threshold
  *
  * References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
  *              Limited memory, SIAM J. Sci. Comput.  21(1), pp. 24-45, 1999
  *
  * \tparam Scalar the scalar type of the input matrices
  * \tparam _UpLo The triangular part that will be used for the computations. It can be Lower
    *               or Upper. Default is Lower.
  * \tparam _OrderingType The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<int>,
  *                       unless EIGEN_MPL2_ONLY is defined, in which case the default is NaturalOrdering<int>.
  *
  * \implsparsesolverconcept
  *
  * It performs the following incomplete factorization: \f$ S P A P' S \approx L L' \f$
  * where L is a lower triangular factor, S is a diagonal scaling matrix, and P is a
  * fill-in reducing permutation as computed by the ordering method.
  *
  * \b Shifting \b strategy: Let \f$ B = S P A P' S \f$  be the scaled matrix on which the factorization is carried out,
  * and \f$ \beta \f$ be the minimum value of the diagonal. If \f$ \beta > 0 \f$ then, the factorization is directly performed
  * on the matrix B. Otherwise, the factorization is performed on the shifted matrix \f$ B + (\sigma+|\beta| I \f$ where
  * \f$ \sigma \f$ is the initial shift value as returned and set by setInitialShift() method. The default value is \f$ \sigma = 10^{-3} \f$.
  * If the factorization fails, then the shift in doubled until it succeed or a maximum of ten attempts. If it still fails, as returned by
  * the info() method, then you can either increase the initial shift, or better use another preconditioning technique.
  *
  */
template <typename Scalar, int _UpLo = Lower, typename _OrderingType =
#ifndef EIGEN_MPL2_ONLY
AMDOrdering<int>
#else
NaturalOrdering<int>
#endif
>
class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
{
  protected:
    typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
    using Base::m_isInitialized;
  public:
    typedef typename NumTraits<Scalar>::Real RealScalar; 
    typedef _OrderingType OrderingType;
    typedef typename OrderingType::PermutationType PermutationType;
    typedef typename PermutationType::StorageIndex StorageIndex; 
    typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
    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 };
    enum {
      ColsAtCompileTime = Dynamic,
      MaxColsAtCompileTime = Dynamic
    };
  public:

    /** Default constructor leaving the object in a partly non-initialized stage.
      *
      * You must call compute() or the pair analyzePattern()/factorize() to make it valid.
      *
      * \sa IncompleteCholesky(const MatrixType&)
      */
    IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {}
    
    /** Constructor computing the incomplete factorization for the given matrix \a matrix.
      */
    template<typename MatrixType>
    IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false)
    {
      compute(matrix);
    }
    
    /** \returns number of rows of the factored matrix */
    Index rows() const { return m_L.rows(); }
    
    /** \returns number of columns of the factored matrix */
    Index cols() const { return m_L.cols(); }
    

    /** \brief Reports whether previous computation was successful.
      *
      * It triggers an assertion if \c *this has not been initialized through the respective constructor,
      * or a call to compute() or analyzePattern().
      *
      * \returns \c Success if computation was successful,
      *          \c NumericalIssue if the matrix appears to be negative.
      */
    ComputationInfo info() const
    {
      eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
      return m_info;
    }
    
    /** \brief Set the initial shift parameter \f$ \sigma \f$.
      */
    void setInitialShift(RealScalar shift) { m_initialShift = shift; }
    
    /** \brief Computes the fill reducing permutation vector using the sparsity pattern of \a mat
      */
    template<typename MatrixType>
    void analyzePattern(const MatrixType& mat)
    {
      OrderingType ord; 
      PermutationType pinv;
      ord(mat.template selfadjointView<UpLo>(), pinv); 
      if(pinv.size()>0) m_perm = pinv.inverse();
      else              m_perm.resize(0);
      m_L.resize(mat.rows(), mat.cols());
      m_analysisIsOk = true;
      m_isInitialized = true;
      m_info = Success;
    }
    
    /** \brief Performs the numerical factorization of the input matrix \a mat
      *
      * The method analyzePattern() or compute() must have been called beforehand
      * with a matrix having the same pattern.
      *
      * \sa compute(), analyzePattern()
      */
    template<typename MatrixType>
    void factorize(const MatrixType& mat);
    
    /** Computes or re-computes the incomplete Cholesky factorization of the input matrix \a mat
      *
      * It is a shortcut for a sequential call to the analyzePattern() and factorize() methods.
      *
      * \sa analyzePattern(), factorize()
      */
    template<typename MatrixType>
    void compute(const MatrixType& mat)
    {
      analyzePattern(mat);
      factorize(mat);
    }
    
    // internal
    template<typename Rhs, typename Dest>
    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() * x;
    }

    /** \returns the sparse lower triangular factor L */
    const FactorType& matrixL() const { eigen_assert("m_factorizationIsOk"); return m_L; }

    /** \returns a vector representing the scaling factor S */
    const VectorRx& scalingS() const { eigen_assert("m_factorizationIsOk"); return m_scale; }

    /** \returns the fill-in reducing permutation P (can be empty for a natural ordering) */
    const PermutationType& permutationP() const { eigen_assert("m_analysisIsOk"); return m_perm; }

  protected:
    FactorType m_L;              // The lower part stored in CSC
    VectorRx m_scale;            // The vector for scaling the matrix 
    RealScalar m_initialShift;   // The initial shift parameter
    bool m_analysisIsOk; 
    bool m_factorizationIsOk; 
    ComputationInfo m_info;
    PermutationType m_perm; 

  private:
    inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol); 
}; 

// Based on the following paper:
//   C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
//   Limited memory, SIAM J. Sci. Comput.  21(1), pp. 24-45, 1999
//   http://ftp.mcs.anl.gov/pub/tech_reports/reports/P682.pdf
template<typename Scalar, int _UpLo, typename OrderingType>
template<typename _MatrixType>
void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
{
  using std::sqrt;
  eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 
    
  // 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
  
  // Apply the fill-reducing permutation computed in analyzePattern()
  if (m_perm.rows() == mat.rows() ) // To detect the null permutation
  {
    // 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<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 col_vals(n);   // Store a  nonzero values in each column
  VectorIx col_irow(n);   // Row indices of nonzero elements in each column
  VectorIx col_pattern(n);
  col_pattern.fill(-1);
  StorageIndex col_nnz;
  
  
  // Computes the scaling factors 
  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();

  for (Index j = 0; j < n; ++j)
    if(m_scale(j)>(std::numeric_limits<RealScalar>::min)())
      m_scale(j) = RealScalar(1)/m_scale(j);
    else
      m_scale(j) = 1;

  // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
  
  // Scale and compute the shift for the matrix 
  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);
  }

  FactorType L_save = m_L;
  
  RealScalar shift = 0;
  if(mindiag <= RealScalar(0.))
    shift = m_initialShift - mindiag;

  m_info = NumericalIssue;

  // Try to perform the incomplete factorization using the current shift
  int iter = 0;
  do
  {
    // Apply the shift to the diagonal elements of the matrix
    for (Index j = 0; j < n; j++)
      vals[colPtr[j]] += shift;

    // jki version of the Cholesky factorization
    Index j=0;
    for (; j < n; ++j)
    {
      // Left-looking factorization of the j-th column
      // First, load the j-th column into col_vals
      Scalar diag = vals[colPtr[j]];  // It is assumed that only the lower part is stored
      col_nnz = 0;
      for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
      {
        StorageIndex l = rowIdx[i];
        col_vals(col_nnz) = vals[i];
        col_irow(col_nnz) = l;
        col_pattern(l) = col_nnz;
        col_nnz++;
      }
      {
        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++)
        {
          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 (Index i = jk; i < colPtr[*k+1]; i++)
          {
            StorageIndex l = rowIdx[i];
            if(col_pattern[l]<0)
            {
              col_vals(col_nnz) = vals[i] * v_j_jk;
              col_irow[col_nnz] = l;
              col_pattern(l) = col_nnz;
              col_nnz++;
            }
            else
              col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
          }
          updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
        }
      }

      // Scale the current column
      if(numext::real(diag) <= 0)
      {
        if(++iter>=10)
          return;

        // increase shift
        shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
        // restore m_L, col_pattern, and listCol
        vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
        rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
        colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
        col_pattern.fill(-1);
        for(Index i=0; i<n; ++i)
          listCol[i].clear();

        break;
      }

      RealScalar rdiag = sqrt(numext::real(diag));
      vals[colPtr[j]] = rdiag;
      for (Index k = 0; k<col_nnz; ++k)
      {
        Index i = col_irow[k];
        //Scale
        col_vals(k) /= rdiag;
        //Update the remaining diagonals with col_vals
        vals[colPtr[i]] -= numext::abs2(col_vals(k));
      }
      // Select the largest p elements
      // p is the original number of elements in the column (without the diagonal)
      Index p = colPtr[j+1] - colPtr[j] - 1 ;
      Ref<VectorSx> cvals = col_vals.head(col_nnz);
      Ref<VectorIx> cirow = col_irow.head(col_nnz);
      internal::QuickSplit(cvals,cirow, p);
      // Insert the largest p elements in the matrix
      Index cpt = 0;
      for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
      {
        vals[i] = col_vals(cpt);
        rowIdx[i] = col_irow(cpt);
        // restore col_pattern:
        col_pattern(col_irow(cpt)) = -1;
        cpt++;
      }
      // Get the first smallest row index and put it after the diagonal element
      Index jk = colPtr(j)+1;
      updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
    }

    if(j==n)
    {
      m_factorizationIsOk = true;
      m_info = Success;
    }
  } while(m_info!=Success);
}

template<typename Scalar, int _UpLo, typename OrderingType>
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) )
  {
    Index p = colPtr(col+1) - jk;
    Index minpos; 
    rowIdx.segment(jk,p).minCoeff(&minpos);
    minpos += jk;
    if (rowIdx(minpos) != rowIdx(jk))
    {
      //Swap
      std::swap(rowIdx(jk),rowIdx(minpos));
      std::swap(vals(jk),vals(minpos));
    }
    firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
    listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
  }
}

} // end namespace Eigen 

#endif