// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008 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/.
template<typename MatrixType, int UpLo> struct LLT_Traits;
}
/** \ingroup Cholesky_Module * * \class LLT * * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features * * \tparam _MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition * \tparam _UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper. * The other triangular part won't be read. * * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite * matrix A such that A = LL^* = U^*U, where L is lower triangular. * * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like D^*D x = b, * for that purpose, we recommend the Cholesky decomposition without square root which is more stable * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other * situations like generalised eigen problems with hermitian matrices. * * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices, * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations * has a solution. * * Example: \include LLT_example.cpp * Output: \verbinclude LLT_example.out * * \b Performance: for best performance, it is recommended to use a column-major storage format * with the Lower triangular part (the default), or, equivalently, a row-major storage format * with the Upper triangular part. Otherwise, you might get a 20% slowdown for the full factorization * step, and rank-updates can be up to 3 times slower. * * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. * * Note that during the decomposition, only the lower (or upper, as defined by _UpLo) triangular part of A is considered. * Therefore, the strict lower part does not have to store correct values. * * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
*/ template<typename _MatrixType, int _UpLo> class LLT
: public SolverBase<LLT<_MatrixType, _UpLo> >
{ public: typedef _MatrixType MatrixType; typedef SolverBase<LLT> Base; friendclass SolverBase<LLT>;
/** * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to * perform decompositions via LLT::compute(const MatrixType&).
*/
LLT() : m_matrix(), m_isInitialized(false) {}
/** \brief Default Constructor with memory preallocation * * Like the default constructor but with preallocation of the internal data * according to the specified problem \a size. * \sa LLT()
*/ explicit LLT(Index size) : m_matrix(size, size),
m_isInitialized(false) {}
/** \brief Constructs a LLT factorization from a given matrix * * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when * \c MatrixType is a Eigen::Ref. * * \sa LLT(const EigenBase&)
*/ template<typename InputType> explicit LLT(EigenBase<InputType>& matrix)
: m_matrix(matrix.derived()),
m_isInitialized(false)
{
compute(matrix.derived());
}
/** \returns a view of the upper triangular matrix U */ inlinetypename Traits::MatrixU matrixU() const
{
eigen_assert(m_isInitialized && "LLT is not initialized."); return Traits::getU(m_matrix);
}
/** \returns a view of the lower triangular matrix L */ inlinetypename Traits::MatrixL matrixL() const
{
eigen_assert(m_isInitialized && "LLT is not initialized."); return Traits::getL(m_matrix);
}
#ifdef EIGEN_PARSED_BY_DOXYGEN /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * Since this LLT class assumes anyway that the matrix A is invertible, the solution * theoretically exists and is unique regardless of b. * * Example: \include LLT_solve.cpp * Output: \verbinclude LLT_solve.out * * \sa solveInPlace(), MatrixBase::llt(), SelfAdjointView::llt()
*/ template<typename Rhs> inlineconst Solve<LLT, Rhs>
solve(const MatrixBase<Rhs>& b) const; #endif
/** \returns an estimate of the reciprocal condition number of the matrix of * which \c *this is the Cholesky decomposition.
*/
RealScalar rcond() const
{
eigen_assert(m_isInitialized && "LLT is not initialized.");
eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative"); return internal::rcond_estimate_helper(m_l1_norm, *this);
}
/** \returns the LLT decomposition matrix * * TODO: document the storage layout
*/ inlineconst MatrixType& matrixLLT() const
{
eigen_assert(m_isInitialized && "LLT is not initialized."); return m_matrix;
}
MatrixType reconstructedMatrix() const;
/** \brief Reports whether previous computation was successful. * * \returns \c Success if computation was successful, * \c NumericalIssue if the matrix.appears not to be positive definite.
*/
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "LLT is not initialized."); return m_info;
}
/** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix is self-adjoint. * * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as: * \code x = decomposition.adjoint().solve(b) \endcode
*/ const LLT& adjoint() const EIGEN_NOEXCEPT { return *this; };
inline EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); } inline EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); }
/** \internal * Used to compute and store L * The strict upper part is not used and even not initialized.
*/
MatrixType m_matrix;
RealScalar m_l1_norm; bool m_isInitialized;
ComputationInfo m_info;
};
namespace internal {
template<typename Scalar, int UpLo> struct llt_inplace;
Index n = mat.cols();
eigen_assert(mat.rows()==n && vec.size()==n);
TempVectorType temp;
if(sigma>0)
{ // This version is based on Givens rotations. // It is faster than the other one below, but only works for updates, // i.e., for sigma > 0
temp = sqrt(sigma) * vec;
RealScalar x = numext::real(mat.coeff(k,k)); if (k>0) x -= A10.squaredNorm(); if (x<=RealScalar(0)) return k;
mat.coeffRef(k,k) = x = sqrt(x); if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); if (rs>0) A21 /= x;
} return -1;
}
template<typename MatrixType> static Index blocked(MatrixType& m)
{
eigen_assert(m.rows()==m.cols());
Index size = m.rows(); if(size<32) return unblocked(m);
/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix * * \returns a reference to *this * * Example: \include TutorialLinAlgComputeTwice.cpp * Output: \verbinclude TutorialLinAlgComputeTwice.out
*/ template<typename MatrixType, int _UpLo> template<typename InputType>
LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& ;a)
{
check_template_parameters();
eigen_assert(a.rows()==a.cols()); const Index size = a.rows();
m_matrix.resize(size, size); if (!internal::is_same_dense(m_matrix, a.derived()))
m_matrix = a.derived();
// Compute matrix L1 norm = max abs column sum.
m_l1_norm = RealScalar(0); // TODO move this code to SelfAdjointView for (Index col = 0; col < size; ++col) {
RealScalar abs_col_sum; if (_UpLo == Lower)
abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>(); else
abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>(); if (abs_col_sum > m_l1_norm)
m_l1_norm = abs_col_sum;
}
m_isInitialized = true; bool ok = Traits::inplace_decomposition(m_matrix);
m_info = ok ? Success : NumericalIssue;
return *this;
}
/** Performs a rank one update (or dowdate) of the current decomposition. * If A = LL^* before the rank one update, * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector * of same dimension.
*/ template<typename _MatrixType, int _UpLo> template<typename VectorType>
LLT<_MatrixType,_UpLo> & LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
eigen_assert(v.size()==m_matrix.cols());
eigen_assert(m_isInitialized); if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
m_info = NumericalIssue; else
m_info = Success;
/** \internal use x = llt_object.solve(x); * * This is the \em in-place version of solve(). * * \param 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. * * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here. * This function will const_cast it, so constness isn't honored here. * * \sa LLT::solve(), MatrixBase::llt()
*/ template<typename MatrixType, int _UpLo> template<typename Derived> void LLT<MatrixType,_UpLo>::solveInPlace(const MatrixBase<Derived> &bAndX) const
{
eigen_assert(m_isInitialized && "LLT is not initialized.");
eigen_assert(m_matrix.rows()==bAndX.rows());
matrixL().solveInPlace(bAndX);
matrixU().solveInPlace(bAndX);
}
/** \returns the matrix represented by the decomposition, * i.e., it returns the product: L L^*.
* This function is provided for debug purpose. */ template<typename MatrixType, int _UpLo>
MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
{
eigen_assert(m_isInitialized && "LLT is not initialized."); return matrixL() * matrixL().adjoint().toDenseMatrix();
}
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
