IncompleteCholesky.h

// 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/.
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_INCOMPLETE_CHOlESKY_H
#define EIGEN_INCOMPLETE_CHOlESKY_H
 
#include <vector>
#include <list>
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
template <typename Scalar, int UpLo_ = Lower, typename OrderingType_ = AMDOrdering<int> >
class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar, UpLo_, OrderingType_> > {
 protected:
  using Base = SparseSolverBase<IncompleteCholesky<Scalar, UpLo_, OrderingType_>>;
  using Base::m_isInitialized;
 
 public:
  using RealScalar = typename NumTraits<Scalar>::Real;
  using OrderingType = OrderingType_;
  using PermutationType = typename OrderingType::PermutationType;
  using StorageIndex = typename PermutationType::StorageIndex;
  using FactorType = SparseMatrix<Scalar, ColMajor, StorageIndex>;
  using VectorSx = Matrix<Scalar, Dynamic, 1>;
  using VectorRx = Matrix<RealScalar, Dynamic, 1>;
  using VectorIx = Matrix<StorageIndex, Dynamic, 1>;
  using VectorList = std::vector<std::list<StorageIndex>>;
  enum { UpLo = UpLo_ };
  enum { ColsAtCompileTime = Dynamic, MaxColsAtCompileTime = Dynamic };
 
 public:
  IncompleteCholesky() : m_initialShift(1e-3), m_analysisIsOk(false), m_factorizationIsOk(false) {}

  template <typename MatrixType>
  IncompleteCholesky(const MatrixType& matrix)
      : m_initialShift(1e-3), m_analysisIsOk(false), m_factorizationIsOk(false) {
    compute(matrix);
  }

  constexpr Index rows() const noexcept { return m_L.rows(); }

  constexpr Index cols() const noexcept { return m_L.cols(); }

  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
    return m_info;
  }

  void setInitialShift(RealScalar shift) { m_initialShift = shift; }

  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;
  }

  template <typename MatrixType>
  void factorize(const MatrixType& mat);

  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;
  }

  const FactorType& matrixL() const {
    eigen_assert(m_factorizationIsOk && "factorize() should be called first");
    return m_L;
  }

  const VectorRx& scalingS() const {
    eigen_assert(m_factorizationIsOk && "factorize() should be called first");
    return m_scale;
  }

  const PermutationType& permutationP() const {
    eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
    return m_perm;
  }

  RealScalar shift() const { return m_shift; }
 
 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;
  RealScalar m_shift;  // The final shift parameter.
 
 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_>();
  }
 
  // The algorithm will insert increasingly large shifts on the diagonal until
  // factorization succeeds. Therefore we have to make sure that there is a
  // space in the datastructure to store such values, even if the original
  // matrix has a zero on the diagonal.
  bool modified = false;
  for (Index i = 0; i < mat.cols(); ++i) {
    bool inserted = false;
    m_L.findOrInsertCoeff(i, i, &inserted);
    if (inserted) {
      modified = true;
    }
  }
  if (modified) m_L.makeCompressed();
 
  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 column
  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 the 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 when roughly uniform to speed up solve().
 
  // 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;
 
  m_shift = RealScalar(0);
  if (mindiag <= RealScalar(0.)) m_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]] += m_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
        m_shift = numext::maxi(m_initialShift, RealScalar(2) * m_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