GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Gael Guennebaud 2015-08-04 16:16:02 +02:00
parent 9a4713e505
commit e31fc50280
1 changed files with 103 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

175
unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
View File

@ -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.inverse() * b;
else
x = b;
x = m_scal.asDiagonal() * x;
x = m_L.template triangularView<UnitLower>().solve(x);
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 * 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
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
ScalarType curCol(n); // Store a nonzero values in each column
IndexType irow(n); // Row indices of nonzero elements in each column
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_scale.resize(n);
m_scale.setZero();
for (Index j = 0; j < n; j++)
for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
{
m_scal(j) = m_L.col(j).norm();
m_scal(j) = sqrt(m_scal(j));
}
// 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);
m_scale(j) += numext::abs2(vals(k));
if(rowIdx[k]!=j)
m_scale(rowIdx[k]) += numext::abs2(vals(k));
}
if(mindiag < Scalar(0.)) m_shift = m_shift - mindiag;
// Apply the shift to the diagonal elements of the matrix
for (int j = 0; j < n; j++)
vals[colPtr[j]] += m_shift;
// jki version of the Cholesky factorization
for (int j=0; j < n; ++j)
m_scale = m_scale.cwiseSqrt().cwiseSqrt();
// Scale and compute the shift for the matrix
RealScalar mindiag = NumTraits<RealScalar>::highest();
for (Index j = 0; j < n; j++)
{
//Left-looking factorize the column j
// First, load the jth column into curCol
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);
}
RealScalar shift = 0;
if(mindiag <= RealScalar(0.))
shift = m_initialShift - mindiag;
// 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
for (Index j=0; j < n; ++j)
{
// 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 ;
// 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));
}
}