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

Add a Transpositions class to ease the representation and

Browse Source 
manipulation of permutations as a sequence of transpositions.
Make LDLT use it.
This commit is contained in:
Gael Guennebaud 2010-06-04 23:17:57 +02:00
parent 1ff1bd69ac
commit bfeba41174
5 changed files with 351 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Eigen/Core
View File

@ -260,6 +260,7 @@ using std::size_t;
#include "src/Core/Diagonal.h"
#include "src/Core/DiagonalProduct.h"
#include "src/Core/PermutationMatrix.h"
#include "src/Core/Transpositions.h"
#include "src/Core/Redux.h"
#include "src/Core/Visitor.h"
#include "src/Core/Fuzzy.h"

87
Eigen/src/Cholesky/LDLT.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Keir Mierle <mierle@gmail.com>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
@ -69,9 +69,12 @@ template<typename _MatrixType, int _UpLo> class LDLT
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename MatrixType::Index Index;
typedef typename ei_plain_col_type<MatrixType, Index>::type IntColVectorType;
// typedef typename ei_plain_col_type<MatrixType, Index>::type IntColVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
typedef LDLT_Traits<MatrixType,UpLo> Traits;
/** \brief Default Constructor.
@ -79,7 +82,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
* The default constructor is useful in cases in which the user intends to
* perform decompositions via LDLT::compute(const MatrixType&).
*/
LDLT() : m_matrix(), m_p(), m_transpositions(), m_isInitialized(false) {}
LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {}
/** \brief Default Constructor with memory preallocation
*
@ -89,7 +92,6 @@ template<typename _MatrixType, int _UpLo> class LDLT
*/
LDLT(Index size)
: m_matrix(size, size),
m_p(size),
m_transpositions(size),
m_temporary(size),
m_isInitialized(false)
@ -97,7 +99,6 @@ template<typename _MatrixType, int _UpLo> class LDLT
LDLT(const MatrixType& matrix)
: m_matrix(matrix.rows(), matrix.cols()),
m_p(matrix.rows()),
m_transpositions(matrix.rows()),
m_temporary(matrix.rows()),
m_isInitialized(false)
@ -119,14 +120,12 @@ template<typename _MatrixType, int _UpLo> class LDLT
return Traits::getL(m_matrix);
}
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning,
* see the examples given in the documentation of class FullPivLU.
/** \returns the permutation matrix P as a transposition sequence.
*/
inline const IntColVectorType& permutationP() const
inline const TranspositionType& transpositionsP() const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_p;
return m_transpositions;
}
/** \returns the coefficients of the diagonal matrix D */
@ -195,8 +194,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
* is not stored), and the diagonal entries correspond to D.
*/
MatrixType m_matrix;
IntColVectorType m_p;
IntColVectorType m_transpositions; // FIXME do we really need to store permanently the transpositions?
TranspositionType m_transpositions;
TmpMatrixType m_temporary;
int m_sign;
bool m_isInitialized;
@ -206,18 +204,18 @@ template<int UpLo> struct ei_ldlt_inplace;
template<> struct ei_ldlt_inplace<Lower>
{
template<typename MatrixType, typename Transpositions, typename Workspace>
static bool unblocked(MatrixType& mat, Transpositions& transpositions, Workspace& temp, int* sign=0)
template<typename MatrixType, typename TranspositionType, typename Workspace>
static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
ei_assert(mat.rows()==mat.cols());
const Index size = mat.rows();
if (size <= 1)
{
transpositions.setZero();
transpositions.setIdentity();
if(sign)
*sign = ei_real(mat.coeff(0,0))>0 ? 1:-1;
return true;
@ -272,15 +270,15 @@ template<> struct ei_ldlt_inplace<Lower>
if((rs>0) && (ei_abs(mat.coeffRef(k,k)) > cutoff))
A21 /= mat.coeffRef(k,k);
}
return true;
}
};
template<> struct ei_ldlt_inplace<Upper>
{
template<typename MatrixType, typename Transpositions, typename Workspace>
static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, Transpositions& transpositions, Workspace& temp, int* sign=0)
template<typename MatrixType, typename TranspositionType, typename Workspace>
static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
{
Transpose<MatrixType> matt(mat);
return ei_ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
@ -313,22 +311,12 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
m_matrix = a;
m_p.resize(size);
m_transpositions.resize(size);
m_isInitialized = false;
m_temporary.resize(size);
ei_ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
if(size==0)
m_p.setZero();
// Reverse applied swaps to get P matrix.
for(Index k = 0; k < size; ++k) m_p.coeffRef(k) = k;
for(Index k = size-1; k >= 0; --k) {
std::swap(m_p.coeffRef(k), m_p.coeffRef(m_transpositions.coeff(k)));
}
m_isInitialized = true;
return *this;
}
@ -342,8 +330,21 @@ struct ei_solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dst = rhs();
dec().solveInPlace(dst);
ei_assert(rhs().rows() == dec().matrixLDLT().rows());
// dst = P b
dst = dec().transpositionsP() * rhs();
// dst = L^-1 (P b)
dec().matrixL().solveInPlace(dst);
// dst = D^-1 (L^-1 P b)
dst = dec().vectorD().asDiagonal().inverse() * dst;
// dst = L^-T (D^-1 L^-1 P b)
dec().matrixU().solveInPlace(dst);
// dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b
dst = dec().transpositionsP().transpose() * dst;
}
};
@ -366,21 +367,7 @@ bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
const Index size = m_matrix.rows();
ei_assert(size == bAndX.rows());
// z = P b
for(Index i = 0; i < size; ++i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i));
// y = L^-1 z
//matrixL().solveInPlace(bAndX);
matrixL().solveInPlace(bAndX);
// w = D^-1 y
bAndX = m_matrix.diagonal().asDiagonal().inverse() * bAndX;
// u = L^-T w
matrixU().solveInPlace(bAndX);
// x = P^T u
for (Index i = size-1; i >= 0; --i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i));
bAndX = this->solve(bAndX);
return true;
}
@ -394,10 +381,10 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
ei_assert(m_isInitialized && "LDLT is not initialized.");
const Index size = m_matrix.rows();
MatrixType res(size,size);
res.setIdentity();
// PI
for(Index i = 0; i < size; ++i) res.row(m_transpositions.coeff(i)).swap(res.row(i));
// P
res.setIdentity();
res = transpositionsP() * res;
// L^* P
res = matrixU() * res;
// D(L^*P)
@ -405,7 +392,7 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
// L(DL^*P)
res = matrixL() * res;
// P^T (LDL^*P)
for (Index i = size-1; i >= 0; --i) res.row(m_transpositions.coeff(i)).swap(res.row(i));
res = transpositionsP().transpose() * res;
return res;
}

1
Eigen/src/Core/PermutationMatrix.h
View File

@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public

285
Eigen/src/Core/Transpositions.h Normal file
View File

@ -0,0 +1,285 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_TRANSPOSITIONS_H
#define EIGEN_TRANSPOSITIONS_H
/** \class Transpositions
*
* \brief Represents a sequence of transpositions (row/column interchange)
*
* \param SizeAtCompileTime the number of transpositions, or Dynamic
* \param MaxSizeAtCompileTime the maximum number of transpositions, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it.
*
* This class represents a permutation transformation as a sequence of \em n transpositions
* \f$[T_{n-1} \ldots T_{i} \ldots T_{0}]\f$. It is internally stored as a vector of integers \c indices.
* Each transposition \f$ T_{i} \f$ applied on the left of a matrix (\f$ T_{i} M\f$) interchanges
* the rows \c i and \c indices[i] of the matrix \c M.
* A transposition applied on the right (e.g., \f$ M T_{i}\f$) yields a column interchange.
*
* Compared to the class PermutationMatrix, such a sequence of transpositions is what is
* computed during a decomposition with pivoting, and it is faster when applying the permutation in-place.
*
* To apply a sequence of transpositions to a matrix, simply use the operator * as in the following example:
* \code
* Transpositions tr;
* MatrixXf mat;
* mat = tr * mat;
* \endcode
* In this example, we detect that the matrix appears on both side, and so the transpositions
* are applied in-place without any temporary or extra copy.
*
* \sa class PermutationMatrix
*/
template<int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime> class Transpositions;
template<typename TranspositionType, typename MatrixType, int Side, bool Transposed=false> struct ei_transposition_matrix_product_retval;
template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
class Transpositions
{
public:
typedef Matrix<DenseIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
typedef typename IndicesType::Index Index;
inline Transpositions() {}
/** Copy constructor. */
template<int OtherSize, int OtherMaxSize>
inline Transpositions(const Transpositions<OtherSize, OtherMaxSize>& other)
: m_indices(other.indices()) {}
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** Standard copy constructor. Defined only to prevent a default copy constructor
* from hiding the other templated constructor */
inline Transpositions(const Transpositions& other) : m_indices(other.indices()) {}
#endif
/** Generic constructor from expression of the transposition indices. */
template<typename Other>
explicit inline Transpositions(const MatrixBase<Other>& indices) : m_indices(indices)
{}
/** Copies the \a other transpositions into \c *this */
template<int OtherSize, int OtherMaxSize>
Transpositions& operator=(const Transpositions<OtherSize, OtherMaxSize>& other)
{
m_indices = other.indices();
return *this;
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** This is a special case of the templated operator=. Its purpose is to
* prevent a default operator= from hiding the templated operator=.
*/
Transpositions& operator=(const Transpositions& other)
{
m_indices = other.m_indices;
return *this;
}
#endif
/** Constructs an uninitialized permutation matrix of given size.
*/
inline Transpositions(Index size) : m_indices(size)
{}
/** \returns the number of transpositions */
inline Index size() const { return m_indices.size(); }
inline const Index& coeff(Index i) const { return m_indices.coeff(i); }
inline Index& coeffRef(Index i) { return m_indices.coeffRef(i); }
inline const Index& operator()(Index i) const { return m_indices(i); }
inline Index& operator()(Index i) { return m_indices(i); }
/** const version of indices(). */
const IndicesType& indices() const { return m_indices; }
/** \returns a reference to the stored array representing the transpositions. */
IndicesType& indices() { return m_indices; }
/** Resizes to given size. */
inline void resize(int size)
{
m_indices.resize(size);
}
/** Sets \c *this to represents an identity transformation */
void setIdentity()
{
for(int i = 0; i < m_indices.size(); ++i)
m_indices.coeffRef(i) = i;
}
// FIXME: do we want such methods ?
// might be usefull when the target matrix expression is complex, e.g.:
// object.matrix().block(..,..,..,..) = trans * object.matrix().block(..,..,..,..);
/*
template<typename MatrixType>
void applyForwardToRows(MatrixType& mat) const
{
for(Index k=0 ; k<size() ; ++k)
if(m_indices(k)!=k)
mat.row(k).swap(mat.row(m_indices(k)));
}
template<typename MatrixType>
void applyBackwardToRows(MatrixType& mat) const
{
for(Index k=size()-1 ; k>=0 ; --k)
if(m_indices(k)!=k)
mat.row(k).swap(mat.row(m_indices(k)));
}
*/
/** \returns the inverse transformation */
inline Transpose<Transpositions> inverse() const
{ return *this; }
/** \returns the tranpose transformation */
inline Transpose<Transpositions> transpose() const
{ return *this; }
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<int OtherSize, int OtherMaxSize>
Transpositions(const Transpose<Transpositions<OtherSize,OtherMaxSize> >& other)
: m_indices(other.size())
{
Index n = size();
Index j = size-1;
for(Index i=0; i<n;++i,--j)
m_indices.coeffRef(j) = other.nestedTranspositions().indices().coeff(i);
}
#endif
protected:
IndicesType m_indices;
};
/** \returns the \a matrix with the \a transpositions applied to the columns.
*/
template<typename Derived, int SizeAtCompileTime, int MaxSizeAtCompileTime>
inline const ei_transposition_matrix_product_retval<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheRight>
operator*(const MatrixBase<Derived>& matrix,
const Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime> &transpositions)
{
return ei_transposition_matrix_product_retval
<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheRight>
(transpositions, matrix.derived());
}
/** \returns the \a matrix with the \a transpositions applied to the rows.
*/
template<typename Derived, int SizeAtCompileTime, int MaxSizeAtCompileTime>
inline const ei_transposition_matrix_product_retval
<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheLeft>
operator*(const Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime> &transpositions,
const MatrixBase<Derived>& matrix)
{
return ei_transposition_matrix_product_retval
<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheLeft>
(transpositions, matrix.derived());
}
template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
struct ei_traits<ei_transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
{
typedef typename MatrixType::PlainObject ReturnType;
};
template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
struct ei_transposition_matrix_product_retval
: public ReturnByValue<ei_transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
{
typedef typename ei_cleantype<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
typedef typename TranspositionType::Index Index;
ei_transposition_matrix_product_retval(const TranspositionType& tr, const MatrixType& matrix)
: m_transpositions(tr), m_matrix(matrix)
{}
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
template<typename Dest> inline void evalTo(Dest& dst) const
{
const int size = m_transpositions.size();
Index j = 0;
if(!(ei_is_same_type<MatrixTypeNestedCleaned,Dest>::ret && ei_extract_data(dst) == ei_extract_data(m_matrix)))
dst = m_matrix;
for(int k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
if((j=m_transpositions.coeff(k))!=k)
{
if(Side==OnTheLeft)
dst.row(k).swap(dst.row(j));
else if(Side==OnTheRight)
dst.col(k).swap(dst.col(j));
}
}
protected:
const TranspositionType& m_transpositions;
const typename MatrixType::Nested m_matrix;
};
/* Template partial specialization for transposed/inverse transpositions */
template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
class Transpose<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime> >
{
typedef Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime> TranspositionType;
typedef typename TranspositionType::IndicesType IndicesType;
public:
Transpose(const TranspositionType& t) : m_transpositions(t) {}
inline int size() const { return m_transpositions.size(); }
/** \returns the \a matrix with the inverse transpositions applied to the columns.
*/
template<typename Derived> friend
inline const ei_transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>
operator*(const MatrixBase<Derived>& matrix, const Transpose& trt)
{
return ei_transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>(trt.m_transpositions, matrix.derived());
}
/** \returns the \a matrix with the inverse transpositions applied to the rows.
*/
template<typename Derived>
inline const ei_transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>
operator*(const MatrixBase<Derived>& matrix) const
{
return ei_transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>(m_transpositions, matrix.derived());
}
const TranspositionType& nestedTranspositions() const { return m_transpositions; }
protected:
const TranspositionType& m_transpositions;
};
#endif // EIGEN_TRANSPOSITIONS_H

27
test/cholesky.cpp
View File

@ -26,10 +26,21 @@
#define EIGEN_NO_ASSERTION_CHECKING
#endif
static int nb_temporaries;
#define EIGEN_DEBUG_MATRIX_CTOR { if(size!=0) nb_temporaries++; }
#include "main.h"
#include <Eigen/Cholesky>
#include <Eigen/QR>
#define VERIFY_EVALUATION_COUNT(XPR,N) {\
nb_temporaries = 0; \
XPR; \
if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \
VERIFY( (#XPR) && nb_temporaries==N ); \
}
#ifdef HAS_GSL
#include "gsl_helper.h"
#endif
@ -131,6 +142,21 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
VERIFY_IS_APPROX(symm * vecX, vecB);
matX = ldltup.solve(matB);
VERIFY_IS_APPROX(symm * matX, matB);
if(MatrixType::RowsAtCompileTime==Dynamic)
{
// note : each inplace permutation requires a small temporary vector (mask)
// check inplace solve
matX = matB;
VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0);
VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval());
matX = matB;
VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
}
}
}
@ -141,6 +167,7 @@ template<typename MatrixType> void cholesky_verify_assert()
LLT<MatrixType> llt;
VERIFY_RAISES_ASSERT(llt.matrixL())
VERIFY_RAISES_ASSERT(llt.matrixU())
VERIFY_RAISES_ASSERT(llt.solve(tmp))
VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp))