MOAB
4.9.3pre
|
Robust Cholesky decomposition of a matrix with pivoting. More...
#include <LDLT.h>
Public Types | |
enum | { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options & ~RowMajorBit, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, UpLo = _UpLo } |
typedef _MatrixType | MatrixType |
typedef MatrixType::Scalar | Scalar |
typedef NumTraits< typename MatrixType::Scalar >::Real | RealScalar |
typedef Eigen::Index | Index |
typedef MatrixType::StorageIndex | StorageIndex |
typedef Matrix< Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1 > | TmpMatrixType |
typedef Transpositions < RowsAtCompileTime, MaxRowsAtCompileTime > | TranspositionType |
typedef PermutationMatrix < RowsAtCompileTime, MaxRowsAtCompileTime > | PermutationType |
typedef internal::LDLT_Traits < MatrixType, UpLo > | Traits |
Public Member Functions | |
LDLT () | |
Default Constructor. | |
LDLT (Index size) | |
Default Constructor with memory preallocation. | |
template<typename InputType > | |
LDLT (const EigenBase< InputType > &matrix) | |
Constructor with decomposition. | |
void | setZero () |
Traits::MatrixU | matrixU () const |
Traits::MatrixL | matrixL () const |
const TranspositionType & | transpositionsP () const |
Diagonal< const MatrixType > | vectorD () const |
bool | isPositive () const |
bool | isNegative (void) const |
template<typename Rhs > | |
const Solve< LDLT, Rhs > | solve (const MatrixBase< Rhs > &b) const |
template<typename Derived > | |
bool | solveInPlace (MatrixBase< Derived > &bAndX) const |
template<typename InputType > | |
LDLT & | compute (const EigenBase< InputType > &matrix) |
template<typename Derived > | |
LDLT & | rankUpdate (const MatrixBase< Derived > &w, const RealScalar &alpha=1) |
const MatrixType & | matrixLDLT () const |
MatrixType | reconstructedMatrix () const |
Index | rows () const |
Index | cols () const |
ComputationInfo | info () const |
Reports whether previous computation was successful. | |
template<typename RhsType , typename DstType > | |
EIGEN_DEVICE_FUNC void | _solve_impl (const RhsType &rhs, DstType &dst) const |
template<typename Derived > | |
LDLT< MatrixType, _UpLo > & | rankUpdate (const MatrixBase< Derived > &w, const typename LDLT< MatrixType, _UpLo >::RealScalar &sigma) |
Static Protected Member Functions | |
static void | check_template_parameters () |
Protected Attributes | |
MatrixType | m_matrix |
TranspositionType | m_transpositions |
TmpMatrixType | m_temporary |
internal::SignMatrix | m_sign |
bool | m_isInitialized |
Robust Cholesky decomposition of a matrix with pivoting.
_MatrixType | the type of the matrix of which to compute the LDL^T Cholesky decomposition |
_UpLo | the triangular part that will be used for the decompositon: Lower (default) or Upper. The other triangular part won't be read. |
Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite matrix such that
, where P is a permutation matrix, L is lower triangular with a unit diagonal and D is a diagonal matrix.
The decomposition uses pivoting to ensure stability, so that L will have zeros in the bottom right rank(A) - n submatrix. Avoiding the square root on D also stabilizes the computation.
Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky decomposition to determine whether a system of equations has a solution.
typedef Eigen::Index Eigen::LDLT< _MatrixType, _UpLo >::Index |
typedef _MatrixType Eigen::LDLT< _MatrixType, _UpLo >::MatrixType |
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> Eigen::LDLT< _MatrixType, _UpLo >::PermutationType |
typedef NumTraits<typename MatrixType::Scalar>::Real Eigen::LDLT< _MatrixType, _UpLo >::RealScalar |
typedef MatrixType::Scalar Eigen::LDLT< _MatrixType, _UpLo >::Scalar |
typedef MatrixType::StorageIndex Eigen::LDLT< _MatrixType, _UpLo >::StorageIndex |
typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> Eigen::LDLT< _MatrixType, _UpLo >::TmpMatrixType |
typedef internal::LDLT_Traits<MatrixType,UpLo> Eigen::LDLT< _MatrixType, _UpLo >::Traits |
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> Eigen::LDLT< _MatrixType, _UpLo >::TranspositionType |
anonymous enum |
RowsAtCompileTime | |
ColsAtCompileTime | |
Options | |
MaxRowsAtCompileTime | |
MaxColsAtCompileTime | |
UpLo |
Definition at line 52 of file LDLT.h.
{ RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here! MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, UpLo = _UpLo };
Eigen::LDLT< _MatrixType, _UpLo >::LDLT | ( | ) | [inline] |
Default Constructor.
The default constructor is useful in cases in which the user intends to perform decompositions via LDLT::compute(const MatrixType&).
Definition at line 76 of file LDLT.h.
: m_matrix(), m_transpositions(), m_sign(internal::ZeroSign), m_isInitialized(false) {}
Eigen::LDLT< _MatrixType, _UpLo >::LDLT | ( | Index | size | ) | [inline, explicit] |
Default Constructor with memory preallocation.
Like the default constructor but with preallocation of the internal data according to the specified problem size.
Definition at line 89 of file LDLT.h.
: m_matrix(size, size), m_transpositions(size), m_temporary(size), m_sign(internal::ZeroSign), m_isInitialized(false) {}
Eigen::LDLT< _MatrixType, _UpLo >::LDLT | ( | const EigenBase< InputType > & | matrix | ) | [inline, explicit] |
Constructor with decomposition.
This calculates the decomposition for the input matrix.
Definition at line 103 of file LDLT.h.
: m_matrix(matrix.rows(), matrix.cols()), m_transpositions(matrix.rows()), m_temporary(matrix.rows()), m_sign(internal::ZeroSign), m_isInitialized(false) { compute(matrix.derived()); }
void Eigen::LDLT< _MatrixType, _UpLo >::_solve_impl | ( | const RhsType & | rhs, |
DstType & | dst | ||
) | const |
Definition at line 488 of file LDLT.h.
{ eigen_assert(rhs.rows() == rows()); // dst = P b dst = m_transpositions * rhs; // dst = L^-1 (P b) matrixL().solveInPlace(dst); // dst = D^-1 (L^-1 P b) // more precisely, use pseudo-inverse of D (see bug 241) using std::abs; const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD()); // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon // as motivated by LAPACK's xGELSS: // RealScalar tolerance = numext::maxi(vecD.array().abs().maxCoeff() * NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest()); // However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest // diagonal element is not well justified and leads to numerical issues in some cases. // Moreover, Lapack's xSYTRS routines use 0 for the tolerance. RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest(); for (Index i = 0; i < vecD.size(); ++i) { if(abs(vecD(i)) > tolerance) dst.row(i) /= vecD(i); else dst.row(i).setZero(); } // dst = L^-T (D^-1 L^-1 P b) matrixU().solveInPlace(dst); // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b dst = m_transpositions.transpose() * dst; }
static void Eigen::LDLT< _MatrixType, _UpLo >::check_template_parameters | ( | ) | [inline, static, protected] |
Index Eigen::LDLT< _MatrixType, _UpLo >::cols | ( | ) | const [inline] |
LDLT< MatrixType, _UpLo > & Eigen::LDLT< MatrixType, _UpLo >::compute | ( | const EigenBase< InputType > & | a | ) |
Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of matrix
Definition at line 433 of file LDLT.h.
{ check_template_parameters(); eigen_assert(a.rows()==a.cols()); const Index size = a.rows(); m_matrix = a.derived(); m_transpositions.resize(size); m_isInitialized = false; m_temporary.resize(size); m_sign = internal::ZeroSign; internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign); m_isInitialized = true; return *this; }
ComputationInfo Eigen::LDLT< _MatrixType, _UpLo >::info | ( | ) | const [inline] |
Reports whether previous computation was successful.
Success
if computation was succesful, NumericalIssue
if the matrix.appears to be negative. Definition at line 218 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); return Success; }
bool Eigen::LDLT< _MatrixType, _UpLo >::isNegative | ( | void | ) | const [inline] |
Definition at line 158 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign; }
bool Eigen::LDLT< _MatrixType, _UpLo >::isPositive | ( | ) | const [inline] |
Definition at line 151 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign; }
Traits::MatrixL Eigen::LDLT< _MatrixType, _UpLo >::matrixL | ( | ) | const [inline] |
Definition at line 129 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); return Traits::getL(m_matrix); }
const MatrixType& Eigen::LDLT< _MatrixType, _UpLo >::matrixLDLT | ( | ) | const [inline] |
TODO: document the storage layout
Definition at line 202 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); return m_matrix; }
Traits::MatrixU Eigen::LDLT< _MatrixType, _UpLo >::matrixU | ( | ) | const [inline] |
Definition at line 122 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); return Traits::getU(m_matrix); }
LDLT& Eigen::LDLT< _MatrixType, _UpLo >::rankUpdate | ( | const MatrixBase< Derived > & | w, |
const RealScalar & | alpha = 1 |
||
) |
LDLT<MatrixType,_UpLo>& Eigen::LDLT< _MatrixType, _UpLo >::rankUpdate | ( | const MatrixBase< Derived > & | w, |
const typename LDLT< MatrixType, _UpLo >::RealScalar & | sigma | ||
) |
Update the LDLT decomposition: given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T.
w | a vector to be incorporated into the decomposition. |
sigma | a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1. |
Definition at line 460 of file LDLT.h.
{ typedef typename TranspositionType::StorageIndex IndexType; const Index size = w.rows(); if (m_isInitialized) { eigen_assert(m_matrix.rows()==size); } else { m_matrix.resize(size,size); m_matrix.setZero(); m_transpositions.resize(size); for (Index i = 0; i < size; i++) m_transpositions.coeffRef(i) = IndexType(i); m_temporary.resize(size); m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef; m_isInitialized = true; } internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma); return *this; }
MatrixType Eigen::LDLT< MatrixType, _UpLo >::reconstructedMatrix | ( | ) | const |
Definition at line 554 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); const Index size = m_matrix.rows(); MatrixType res(size,size); // P res.setIdentity(); res = transpositionsP() * res; // L^* P res = matrixU() * res; // D(L^*P) res = vectorD().real().asDiagonal() * res; // L(DL^*P) res = matrixL() * res; // P^T (LDL^*P) res = transpositionsP().transpose() * res; return res; }
Index Eigen::LDLT< _MatrixType, _UpLo >::rows | ( | ) | const [inline] |
void Eigen::LDLT< _MatrixType, _UpLo >::setZero | ( | ) | [inline] |
Clear any existing decomposition
Definition at line 116 of file LDLT.h.
{ m_isInitialized = false; }
const Solve<LDLT, Rhs> Eigen::LDLT< _MatrixType, _UpLo >::solve | ( | const MatrixBase< Rhs > & | b | ) | const [inline] |
This function also supports in-place solves using the syntax x = decompositionObject.solve(x)
.
More precisely, this method solves using the decomposition
by solving the systems
,
,
,
and
in succession. If the matrix
is singular, then
will also be singular (all the other matrices are invertible). In that case, the least-square solution of
is computed. This does not mean that this function computes the least-square solution of
is
is singular.
Definition at line 181 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); eigen_assert(m_matrix.rows()==b.rows() && "LDLT::solve(): invalid number of rows of the right hand side matrix b"); return Solve<LDLT, Rhs>(*this, b.derived()); }
bool Eigen::LDLT< MatrixType, _UpLo >::solveInPlace | ( | MatrixBase< Derived > & | bAndX | ) | const |
use x = ldlt_object.solve(x);
This is the in-place version of solve().
bAndX | represents both the right-hand side matrix b and result x. |
This version avoids a copy when the right hand side matrix b is not needed anymore.
Definition at line 540 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); eigen_assert(m_matrix.rows() == bAndX.rows()); bAndX = this->solve(bAndX); return true; }
const TranspositionType& Eigen::LDLT< _MatrixType, _UpLo >::transpositionsP | ( | ) | const [inline] |
Definition at line 137 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); return m_transpositions; }
Diagonal<const MatrixType> Eigen::LDLT< _MatrixType, _UpLo >::vectorD | ( | ) | const [inline] |
Definition at line 144 of file LDLT.h.
{ eigen_assert(m_isInitialized && "LDLT is not initialized."); return m_matrix.diagonal(); }
bool Eigen::LDLT< _MatrixType, _UpLo >::m_isInitialized [protected] |
MatrixType Eigen::LDLT< _MatrixType, _UpLo >::m_matrix [protected] |
Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U. The strict upper part is used during the decomposition, the strict lower part correspond to the coefficients of L (its diagonal is equal to 1 and is not stored), and the diagonal entries correspond to D.
internal::SignMatrix Eigen::LDLT< _MatrixType, _UpLo >::m_sign [protected] |
TmpMatrixType Eigen::LDLT< _MatrixType, _UpLo >::m_temporary [protected] |
TranspositionType Eigen::LDLT< _MatrixType, _UpLo >::m_transpositions [protected] |