From 765219aa5180c40be86a380524db691a0fd5861b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 13 Oct 2008 15:53:27 +0000 Subject: [PATCH] Big API change in Cholesky module: * rename Cholesky to LLT * rename CholeskyWithoutSquareRoot to LDLT * rename MatrixBase::cholesky() to llt() * rename MatrixBase::choleskyNoSqrt() to ldlt() * make {LLT,LDLT}::solve() API consistent with other modules Note that we are going to keep a source compatibility untill the next beta release. E.g., the "old" Cholesky* classes, etc are still available for some time. To be clear, Eigen beta2 should be (hopefully) source compatible with beta1, and so beta2 will contain all the deprecated API of beta1. Those features marked as deprecated will be removed in beta3 (or in the final 2.0 if there is no beta 3 !). Also includes various updated in sparse Cholesky. --- Eigen/Cholesky | 6 +- Eigen/src/Cholesky/Cholesky.h | 56 +++-- .../src/Cholesky/CholeskyWithoutSquareRoot.h | 65 +++--- Eigen/src/Cholesky/LDLT.h | 202 ++++++++++++++++++ Eigen/src/Cholesky/LLT.h | 186 ++++++++++++++++ Eigen/src/Core/MatrixBase.h | 3 + Eigen/src/Core/util/ForwardDeclarations.h | 3 + Eigen/src/Core/util/Macros.h | 6 + Eigen/src/LU/LU.h | 3 +- Eigen/src/QR/QrInstantiations.cpp | 2 - Eigen/src/QR/SelfAdjointEigenSolver.h | 2 +- Eigen/src/SVD/SVD.h | 2 +- Eigen/src/Sparse/AmbiVector.h | 12 +- Eigen/src/Sparse/CholmodSupport.h | 22 +- Eigen/src/Sparse/SparseCholesky.h | 72 +++---- Eigen/src/Sparse/TaucsSupport.h | 27 ++- Eigen/src/Sparse/TriangularSolver.h | 2 +- bench/BenchSparseUtil.h | 2 +- bench/benchCholesky.cpp | 30 +-- bench/sparse_product.cpp | 50 +++-- .../{Cholesky_solve.cpp => LLT_solve.cpp} | 6 +- test/cholesky.cpp | 30 +-- test/sparse.cpp | 2 +- 23 files changed, 608 insertions(+), 183 deletions(-) create mode 100644 Eigen/src/Cholesky/LDLT.h create mode 100644 Eigen/src/Cholesky/LLT.h rename doc/snippets/{Cholesky_solve.cpp => LLT_solve.cpp} (57%) diff --git a/Eigen/Cholesky b/Eigen/Cholesky index f5bcec52d..ddde61a31 100644 --- a/Eigen/Cholesky +++ b/Eigen/Cholesky @@ -17,8 +17,8 @@ namespace Eigen { /** \defgroup Cholesky_Module Cholesky module * This module provides two variants of the Cholesky decomposition for selfadjoint (hermitian) matrices. * Those decompositions are accessible via the following MatrixBase methods: - * - MatrixBase::cholesky(), - * - MatrixBase::choleskyNoSqrt() + * - MatrixBase::llt(), + * - MatrixBase::ldlt() * * \code * #include @@ -27,6 +27,8 @@ namespace Eigen { #include "src/Array/CwiseOperators.h" #include "src/Array/Functors.h" +#include "src/Cholesky/LLT.h" +#include "src/Cholesky/LDLT.h" #include "src/Cholesky/Cholesky.h" #include "src/Cholesky/CholeskyWithoutSquareRoot.h" diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h index b94fea8dc..ada413b33 100644 --- a/Eigen/src/Cholesky/Cholesky.h +++ b/Eigen/src/Cholesky/Cholesky.h @@ -29,22 +29,7 @@ * * \class Cholesky * - * \brief Standard Cholesky decomposition of a matrix and associated features - * - * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition - * - * This class performs a standard 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. - * - * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, - * the strict lower part does not have to store correct values. - * - * \sa MatrixBase::cholesky(), class CholeskyWithoutSquareRoot + * \deprecated this class has been renamed LLT */ template class Cholesky { @@ -72,7 +57,10 @@ template class Cholesky inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } template - typename Derived::Eval solve(const MatrixBase &b) const; + typename Derived::Eval solve(const MatrixBase &b) const EIGEN_DEPRECATED; + + template + bool solve(const MatrixBase &b, MatrixBase *result) const; template bool solveInPlace(MatrixBase &bAndX) const; @@ -128,16 +116,7 @@ void Cholesky::compute(const MatrixType& a) } } -/** \returns the solution of \f$ A x = b \f$ using the current decomposition of A. - * In other words, it returns \f$ A^{-1} b \f$ computing - * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left. - * \param b the column vector \f$ b \f$, which can also be a matrix. - * - * Example: \include Cholesky_solve.cpp - * Output: \verbinclude Cholesky_solve.out - * - * \sa MatrixBase::cholesky(), CholeskyWithoutSquareRoot::solve() - */ +/** \deprecated */ template template typename Derived::Eval Cholesky::solve(const MatrixBase &b) const @@ -147,7 +126,6 @@ typename Derived::Eval Cholesky::solve(const MatrixBase &b) typename ei_eval_to_column_major::type x(b); solveInPlace(x); return x; - //return m_matrix.adjoint().template part().solveTriangular(matrixL().solveTriangular(b)); } /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. @@ -162,7 +140,25 @@ typename Derived::Eval Cholesky::solve(const MatrixBase &b) * Example: \include Cholesky_solve.cpp * Output: \verbinclude Cholesky_solve.out * - * \sa MatrixBase::cholesky(), Cholesky::solve() + * \sa MatrixBase::cholesky(), Cholesky::solveInPlace() + */ +template +template +bool Cholesky::solve(const MatrixBase &b, MatrixBase *result) const +{ + const int size = m_matrix.rows(); + ei_assert(size==b.rows() && "Cholesky::solve(): invalid number of rows of the right hand side matrix b"); + return solveInPlace((*result) = b); +} + +/** 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. + * + * \sa Cholesky::solve(), MatrixBase::cholesky() */ template template @@ -178,7 +174,7 @@ bool Cholesky::solveInPlace(MatrixBase &bAndX) const } /** \cholesky_module - * \returns the Cholesky decomposition of \c *this + * \deprecated has been renamed llt() */ template inline const Cholesky::EvalType> diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h index fd111fb1f..af44634a0 100644 --- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h +++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h @@ -25,25 +25,11 @@ #ifndef EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H #define EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H -/** \ingroup Cholesky_Module +/** \deprecated \ingroup Cholesky_Module * * \class CholeskyWithoutSquareRoot * - * \brief Robust Cholesky decomposition of a matrix and associated features - * - * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition - * - * This class performs a Cholesky decomposition without square root of a symmetric, positive definite - * matrix A such that A = L D L^* = U^* D U, where L is lower triangular with a unit diagonal - * and D is a diagonal matrix. - * - * Compared to a standard Cholesky decomposition, avoiding the square roots allows for faster and more - * stable computation. - * - * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, - * the strict lower part does not have to store correct values. - * - * \sa MatrixBase::choleskyNoSqrt(), class Cholesky + * \deprecated this class has been renamed LDLT */ template class CholeskyWithoutSquareRoot { @@ -69,7 +55,10 @@ template class CholeskyWithoutSquareRoot inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } template - typename Derived::Eval solve(const MatrixBase &b) const; + typename Derived::Eval solve(const MatrixBase &b) const EIGEN_DEPRECATED; + + template + bool solve(const MatrixBase &b, MatrixBase *result) const; template bool solveInPlace(MatrixBase &bAndX) const; @@ -141,15 +130,7 @@ void CholeskyWithoutSquareRoot::compute(const MatrixType& a) } } -/** \returns the solution of \f$ A x = b \f$ using the current decomposition of A. - * In other words, it returns \f$ A^{-1} b \f$ computing - * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left. - * \param b the column vector \f$ b \f$, which can also be a matrix. - * - * See Cholesky::solve() for a example. - * - * \sa CholeskyWithoutSquareRoot::solveInPlace(), MatrixBase::choleskyNoSqrt() - */ +/** \deprecated */ template template typename Derived::Eval CholeskyWithoutSquareRoot::solve(const MatrixBase &b) const @@ -173,10 +154,30 @@ typename Derived::Eval CholeskyWithoutSquareRoot::solve(const Matrix * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left. * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$ * - * Example: \include Cholesky_solve.cpp - * Output: \verbinclude Cholesky_solve.out + * Example: \include CholeskyCholeskyWithoutSquareRoot_solve.cpp + * Output: \verbinclude CholeskyCholeskyWithoutSquareRoot_solve.out * - * \sa MatrixBase::cholesky(), CholeskyWithoutSquareRoot::solve() + * \sa CholeskyWithoutSquareRoot::solveInPlace(), MatrixBase::choleskyNoSqrt() + */ +template +template +bool CholeskyWithoutSquareRoot +::solve(const MatrixBase &b, MatrixBase *result) const +{ + const int size = m_matrix.rows(); + ei_assert(size==b.rows() && "Cholesky::solve(): invalid number of rows of the right hand side matrix b"); + *result = b; + return solveInPlace(*result); +} + +/** 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. + * + * \sa CholeskyWithoutSquareRoot::solve(), MatrixBase::choleskyNoSqrt() */ template template @@ -187,13 +188,13 @@ bool CholeskyWithoutSquareRoot::solveInPlace(MatrixBase &bA if (!m_isPositiveDefinite) return false; matrixL().solveTriangularInPlace(bAndX); - bAndX *= m_matrix.cwise().inverse().template part(); + bAndX = (m_matrix.cwise().inverse().template part() * bAndX).lazy(); m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); return true; } -/** \cholesky_module - * \returns the Cholesky decomposition without square root of \c *this +/** \deprecated \cholesky_module + * \deprecated has been renamed ldlt() */ template inline const CholeskyWithoutSquareRoot::EvalType> diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h new file mode 100644 index 000000000..e70a324f6 --- /dev/null +++ b/Eigen/src/Cholesky/LDLT.h @@ -0,0 +1,202 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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 . + +#ifndef EIGEN_LDLT_H +#define EIGEN_LDLT_H + +/** \ingroup cholesky_Module + * + * \class LDLT + * + * \brief Robust Cholesky decomposition of a matrix and associated features + * + * \param MatrixType the type of the matrix of which we are computing the LDL^T Cholesky decomposition + * + * This class performs a Cholesky decomposition without square root of a symmetric, positive definite + * matrix A such that A = L D L^* = U^* D U, where L is lower triangular with a unit diagonal + * and D is a diagonal matrix. + * + * Compared to a standard Cholesky decomposition, avoiding the square roots allows for faster and more + * stable computation. + * + * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, + * the strict lower part does not have to store correct values. + * + * \sa MatrixBase::ldlt(), class LLT + */ +template class LDLT +{ + public: + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix VectorType; + + LDLT(const MatrixType& matrix) + : m_matrix(matrix.rows(), matrix.cols()) + { + compute(matrix); + } + + /** \returns the lower triangular matrix L */ + inline Part matrixL(void) const { return m_matrix; } + + /** \returns the coefficients of the diagonal matrix D */ + inline DiagonalCoeffs vectorD(void) const { return m_matrix.diagonal(); } + + /** \returns true if the matrix is positive definite */ + inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } + + template + bool solve(const MatrixBase &b, MatrixBase *result) const; + + template + bool solveInPlace(MatrixBase &bAndX) const; + + void compute(const MatrixType& matrix); + + protected: + /** \internal + * 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. + */ + MatrixType m_matrix; + + bool m_isPositiveDefinite; +}; + +/** Compute / recompute the LLT decomposition A = L D L^* = U^* D U of \a matrix + */ +template +void LDLT::compute(const MatrixType& a) +{ + assert(a.rows()==a.cols()); + const int size = a.rows(); + m_matrix.resize(size, size); + m_isPositiveDefinite = true; + const RealScalar eps = ei_sqrt(precision()); + + if (size<=1) + { + m_matrix = a; + return; + } + + // Let's preallocate a temporay vector to evaluate the matrix-vector product into it. + // Unlike the standard LLT decomposition, here we cannot evaluate it to the destination + // matrix because it a sub-row which is not compatible suitable for efficient packet evaluation. + // (at least if we assume the matrix is col-major) + Matrix _temporary(size); + + // Note that, in this algorithm the rows of the strict upper part of m_matrix is used to store + // column vector, thus the strange .conjugate() and .transpose()... + + m_matrix.row(0) = a.row(0).conjugate(); + m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0); + for (int j = 1; j < size; ++j) + { + RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0)); + m_matrix.coeffRef(j,j) = tmp; + + if (tmp < eps) + { + m_isPositiveDefinite = false; + return; + } + + int endSize = size-j-1; + if (endSize>0) + { + _temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j) + * m_matrix.col(j).start(j).conjugate() ).lazy(); + + m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate() + - _temporary.end(endSize).transpose(); + + m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp; + } + } +} + +/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. + * The result is stored in \a bAndx + * + * \returns true in case of success, false otherwise. + * + * In other words, it computes \f$ b = A^{-1} b \f$ with + * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left. + * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$ + * + * Example: \include LLTLDLT_solve.cpp + * Output: \verbinclude LLTLDLT_solve.out + * + * \sa LDLT::solveInPlace(), MatrixBase::ldlt() + */ +template +template +bool LDLT +::solve(const MatrixBase &b, MatrixBase *result) const +{ + const int size = m_matrix.rows(); + ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b"); + *result = b; + return solveInPlace(*result); +} + +/** 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. + * + * \sa LDLT::solve(), MatrixBase::ldlt() + */ +template +template +bool LDLT::solveInPlace(MatrixBase &bAndX) const +{ + const int size = m_matrix.rows(); + ei_assert(size==bAndX.rows()); + if (!m_isPositiveDefinite) + return false; + matrixL().solveTriangularInPlace(bAndX); + bAndX = (m_matrix.cwise().inverse().template part() * bAndX).lazy(); + m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); + return true; +} + +/** \cholesky_module + * \returns the Cholesky decomposition without square root of \c *this + */ +template +inline const LDLT::EvalType> +MatrixBase::ldlt() const +{ + return derived(); +} + +#endif // EIGEN_LDLT_H diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h new file mode 100644 index 000000000..8d4a1a752 --- /dev/null +++ b/Eigen/src/Cholesky/LLT.h @@ -0,0 +1,186 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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 . + +#ifndef EIGEN_LLT_H +#define EIGEN_LLT_H + +/** \ingroup cholesky_Module + * + * \class LLT + * + * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features + * + * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition + * + * 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. + * + * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, + * the strict lower part does not have to store correct values. + * + * \sa MatrixBase::llt(), class LDLT + */ +template class LLT +{ + private: + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix VectorType; + + enum { + PacketSize = ei_packet_traits::size, + AlignmentMask = int(PacketSize)-1 + }; + + public: + + LLT(const MatrixType& matrix) + : m_matrix(matrix.rows(), matrix.cols()) + { + compute(matrix); + } + + inline Part matrixL(void) const { return m_matrix; } + + /** \returns true if the matrix is positive definite */ + inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } + + template + bool solve(const MatrixBase &b, MatrixBase *result) const; + + template + bool solveInPlace(MatrixBase &bAndX) const; + + void compute(const MatrixType& matrix); + + protected: + /** \internal + * Used to compute and store L + * The strict upper part is not used and even not initialized. + */ + MatrixType m_matrix; + bool m_isPositiveDefinite; +}; + +/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix + */ +template +void LLT::compute(const MatrixType& a) +{ + assert(a.rows()==a.cols()); + const int size = a.rows(); + m_matrix.resize(size, size); + const RealScalar eps = ei_sqrt(precision()); + + RealScalar x; + x = ei_real(a.coeff(0,0)); + m_isPositiveDefinite = x > eps && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1)); + m_matrix.coeffRef(0,0) = ei_sqrt(x); + m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0)); + for (int j = 1; j < size; ++j) + { + Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).norm2(); + x = ei_real(tmp); + if (x < eps || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1)))) + { + m_isPositiveDefinite = false; + return; + } + m_matrix.coeffRef(j,j) = x = ei_sqrt(x); + + int endSize = size-j-1; + if (endSize>0) { + // Note that when all matrix columns have good alignment, then the following + // product is guaranteed to be optimal with respect to alignment. + m_matrix.col(j).end(endSize) = + (m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy(); + + // FIXME could use a.col instead of a.row + m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint() + - m_matrix.col(j).end(endSize) ) / x; + } + } +} + +/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. + * The result is stored in \a bAndx + * + * \returns true in case of success, false otherwise. + * + * In other words, it computes \f$ b = A^{-1} b \f$ with + * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left. + * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$ + * + * Example: \include LLT_solve.cpp + * Output: \verbinclude LLT_solve.out + * + * \sa LLT::solveInPlace(), MatrixBase::llt() + */ +template +template +bool LLT::solve(const MatrixBase &b, MatrixBase *result) const +{ + const int size = m_matrix.rows(); + ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b"); + return solveInPlace((*result) = b); +} + +/** 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. + * + * \sa LLT::solve(), MatrixBase::llt() + */ +template +template +bool LLT::solveInPlace(MatrixBase &bAndX) const +{ + const int size = m_matrix.rows(); + ei_assert(size==bAndX.rows()); + if (!m_isPositiveDefinite) + return false; + matrixL().solveTriangularInPlace(bAndX); + m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); + return true; +} + +/** \cholesky_module + * \returns the LLT decomposition of \c *this + */ +template +inline const LLT::EvalType> +MatrixBase::llt() const +{ + return LLT::type>(derived()); +} + +#endif // EIGEN_LLT_H diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 944d353d8..8b0129060 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -563,6 +563,9 @@ template class MatrixBase /////////// Cholesky module /////////// + const LLT llt() const; + const LDLT ldlt() const; + // deprecated: const Cholesky cholesky() const; const CholeskyWithoutSquareRoot choleskyNoSqrt() const; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 55016e05d..98d15b415 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -102,6 +102,9 @@ template class PartialRedux; template class LU; template class QR; template class SVD; +template class LLT; +template class LDLT; +// deprecated: template class Cholesky; template class CholeskyWithoutSquareRoot; diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index e3429e146..348f313d2 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -97,6 +97,12 @@ using Eigen::ei_cos; #define EIGEN_DONT_INLINE #endif +#if (defined __GNUC__) +#define EIGEN_DEPRECATED __attribute__((deprecated)) +#else +#define EIGEN_DEPRECATED +#endif + #if (defined __GNUC__) #define EIGEN_ALIGN_128 __attribute__ ((aligned(16))) #else diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index f7ed1b412..554d8bd3c 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -181,7 +181,8 @@ template class LU * \returns true if any solution exists, false if no solution exists. * * \note If there exist more than one solution, this method will arbitrarily choose one. - * If you need a complete analysis of the space of solutions, take the one solution obtained * by this method and add to it elements of the kernel, as determined by kernel(). + * If you need a complete analysis of the space of solutions, take the one solution obtained + * by this method and add to it elements of the kernel, as determined by kernel(). * * Example: \include LU_solve.cpp * Output: \verbinclude LU_solve.out diff --git a/Eigen/src/QR/QrInstantiations.cpp b/Eigen/src/QR/QrInstantiations.cpp index 3b2594e34..dacb05d3d 100644 --- a/Eigen/src/QR/QrInstantiations.cpp +++ b/Eigen/src/QR/QrInstantiations.cpp @@ -26,8 +26,6 @@ #define EIGEN_EXTERN_INSTANTIATIONS #endif #include "../../Core" -// commented because of -pedantic -// #include "../../Cholesky" #undef EIGEN_EXTERN_INSTANTIATIONS #include "../../QR" diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index 9e929620c..fdb2764c4 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -227,7 +227,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); // Compute the cholesky decomposition of matB = L L' - Cholesky cholB(matB); + LLT cholB(matB); // compute C = inv(L) A inv(L') MatrixType matC = matA; diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index c3f3bb235..100ca9310 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -505,7 +505,7 @@ SVD& SVD::sort() /** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A. * The parts of the solution corresponding to zero singular values are ignored. * - * \sa MatrixBase::svd(), LU::solve(), Cholesky::solve() + * \sa MatrixBase::svd(), LU::solve(), LLT::solve() */ template template diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h index 0c4bd1af1..44697bceb 100644 --- a/Eigen/src/Sparse/AmbiVector.h +++ b/Eigen/src/Sparse/AmbiVector.h @@ -28,7 +28,7 @@ /** \internal * Hybrid sparse/dense vector class designed for intensive read-write operations. * - * See BasicSparseCholesky and SparseProduct for usage examples. + * See BasicSparseLLT and SparseProduct for usage examples. */ template class AmbiVector { @@ -269,12 +269,14 @@ class AmbiVector<_Scalar>::Iterator /** Default constructor * \param vec the vector on which we iterate - * \param nonZeroReferenceValue reference value used to prune zero coefficients. - * In practice, the coefficient are compared to \a nonZeroReferenceValue * precision(). + * \param epsilon the minimal value used to prune zero coefficients. + * In practice, all coefficients having a magnitude smaller than \a epsilon + * are skipped. */ - Iterator(const AmbiVector& vec, RealScalar nonZeroReferenceValue = RealScalar(0.1)) : m_vector(vec) + Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*precision()) + : m_vector(vec) { - m_epsilon = nonZeroReferenceValue * precision(); + m_epsilon = epsilon; m_isDense = m_vector.m_mode==IsDense; if (m_isDense) { diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h index ed031b264..ba2161320 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/Eigen/src/Sparse/CholmodSupport.h @@ -99,10 +99,10 @@ SparseMatrix SparseMatrix::Map(cholmod_sparse& cm) } template -class SparseCholesky : public SparseCholesky +class SparseLLT : public SparseLLT { protected: - typedef SparseCholesky Base; + typedef SparseLLT Base; using Base::Scalar; using Base::RealScalar; using Base::MatrixLIsDirty; @@ -113,14 +113,20 @@ class SparseCholesky : public SparseCholesky public: - SparseCholesky(const MatrixType& matrix, int flags = 0) + SparseLLT(int flags = 0) + : Base(flags), m_cholmodFactor(0) + { + cholmod_start(&m_cholmod); + } + + SparseLLT(const MatrixType& matrix, int flags = 0) : Base(matrix, flags), m_cholmodFactor(0) { cholmod_start(&m_cholmod); compute(matrix); } - ~SparseCholesky() + ~SparseLLT() { if (m_cholmodFactor) cholmod_free_factor(&m_cholmodFactor, &m_cholmod); @@ -140,7 +146,7 @@ class SparseCholesky : public SparseCholesky }; template -void SparseCholesky::compute(const MatrixType& a) +void SparseLLT::compute(const MatrixType& a) { if (m_cholmodFactor) { @@ -169,8 +175,8 @@ void SparseCholesky::compute(const MatrixType& a) } template -inline const typename SparseCholesky::CholMatrixType& -SparseCholesky::matrixL() const +inline const typename SparseLLT::CholMatrixType& +SparseLLT::matrixL() const { if (m_status & MatrixLIsDirty) { @@ -187,7 +193,7 @@ SparseCholesky::matrixL() const template template -void SparseCholesky::solveInPlace(MatrixBase &b) const +void SparseLLT::solveInPlace(MatrixBase &b) const { const int size = m_matrix.rows(); ei_assert(size==b.rows()); diff --git a/Eigen/src/Sparse/SparseCholesky.h b/Eigen/src/Sparse/SparseCholesky.h index 839a84bf3..efe005709 100644 --- a/Eigen/src/Sparse/SparseCholesky.h +++ b/Eigen/src/Sparse/SparseCholesky.h @@ -38,25 +38,19 @@ enum { MemoryEfficient = 0x2, SupernodalMultifrontal = 0x4, SupernodalLeftLooking = 0x8 - - -/* - CholUseEigen = 0x0, // Eigen's impl is the default - CholUseTaucs = 0x2, - CholUseCholmod = 0x4*/ }; /** \ingroup Sparse_Module * - * \class SparseCholesky + * \class SparseLLT * - * \brief Standard Cholesky decomposition of a matrix and associated features + * \brief Standard LLT decomposition of a matrix and associated features * - * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition + * \param MatrixType the type of the matrix of which we are computing the LLT decomposition * - * \sa class Cholesky, class CholeskyWithoutSquareRoot + * \sa class LLT, class LDLT */ -template class SparseCholesky +template class SparseLLT { protected: typedef typename MatrixType::Scalar Scalar; @@ -70,66 +64,66 @@ template class SparseCholesky public: - SparseCholesky(const MatrixType& matrix, int flags = 0) + SparseLLT(int flags = 0) + : m_flags(flags), m_status(0) + { + m_precision = RealScalar(0.1) * Eigen::precision(); + } + + SparseLLT(const MatrixType& matrix, int flags = 0) : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0) { + m_precision = RealScalar(0.1) * Eigen::precision(); compute(matrix); } - inline const CholMatrixType& matrixL(void) const { return m_matrix; } + void setPrecision(RealScalar v) { m_precision = v; } + RealScalar precision() const { return m_precision; } - /** \returns true if the matrix is positive definite */ - inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } + void setFlags(int f) { m_flags = f; } + int flags() const { return m_flags; } + + void compute(const MatrixType& matrix); + + inline const CholMatrixType& matrixL(void) const { return m_matrix; } template void solveInPlace(MatrixBase &b) const; - void compute(const MatrixType& matrix); + /** \returns true if the factorization succeeded */ + inline bool succeeded(void) const { return m_succeeded; } protected: - /** \internal - * Used to compute and store L - * The strict upper part is not used and even not initialized. - */ CholMatrixType m_matrix; + RealScalar m_precision; int m_flags; mutable int m_status; - bool m_isPositiveDefinite; + bool m_succeeded; }; -/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix +/** Computes / recomputes the LLT decomposition A = LL^* = U^*U of \a matrix */ template -void SparseCholesky::compute(const MatrixType& a) +void SparseLLT::compute(const MatrixType& a) { assert(a.rows()==a.cols()); const int size = a.rows(); m_matrix.resize(size, size); - const RealScalar eps = ei_sqrt(precision()); +// const RealScalar eps = ei_sqrt(precision()); // allocate a temporary vector for accumulations AmbiVector tempVector(size); + RealScalar density = a.nonZeros()/RealScalar(size*size); // TODO estimate the number of nnz m_matrix.startFill(a.nonZeros()*2); for (int j = 0; j < size; ++j) { -// std::cout << j << "\n"; Scalar x = ei_real(a.coeff(j,j)); int endSize = size-j-1; - // TODO estimate the number of non zero entries -// float ratioLhs = float(lhs.nonZeros())/float(lhs.rows()*lhs.cols()); -// float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols); -// float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f); - - // let's do a more accurate determination of the nnz ratio for the current column j of res - //float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f); - // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector) -// float ratioColRes = ratioRes; -// if (ratioColRes>0.1) -// tempVector.init(IsSparse); - tempVector.init(IsDense); + // TODO better estimate the density ! + tempVector.init(density>0.001? IsDense : IsSparse); tempVector.setBounds(j+1,size); tempVector.setZero(); // init with current matrix a @@ -161,7 +155,7 @@ void SparseCholesky::compute(const MatrixType& a) RealScalar rx = ei_sqrt(ei_real(x)); m_matrix.fill(j,j) = rx; Scalar y = Scalar(1)/rx; - for (typename AmbiVector::Iterator it(tempVector); it; ++it) + for (typename AmbiVector::Iterator it(tempVector, m_precision*rx); it; ++it) { m_matrix.fill(it.index(), j) = it.value() * y; } @@ -171,7 +165,7 @@ void SparseCholesky::compute(const MatrixType& a) template template -void SparseCholesky::solveInPlace(MatrixBase &b) const +void SparseLLT::solveInPlace(MatrixBase &b) const { const int size = m_matrix.rows(); ei_assert(size==b.rows()); diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h index a0a5b4b71..a92a63b57 100644 --- a/Eigen/src/Sparse/TaucsSupport.h +++ b/Eigen/src/Sparse/TaucsSupport.h @@ -79,10 +79,10 @@ SparseMatrix SparseMatrix::Map(taucs_ccs_matrix& tau } template -class SparseCholesky : public SparseCholesky +class SparseLLT : public SparseLLT { protected: - typedef SparseCholesky Base; + typedef SparseLLT Base; using Base::Scalar; using Base::RealScalar; using Base::MatrixLIsDirty; @@ -93,13 +93,18 @@ class SparseCholesky : public SparseCholesky public: - SparseCholesky(const MatrixType& matrix, int flags = 0) + SparseLLT(int flags = 0) + : Base(flags), m_taucsSupernodalFactor(0) + { + } + + SparseLLT(const MatrixType& matrix, int flags = 0) : Base(matrix, flags), m_taucsSupernodalFactor(0) { compute(matrix); } - ~SparseCholesky() + ~SparseLLT() { if (m_taucsSupernodalFactor) taucs_supernodal_factor_free(m_taucsSupernodalFactor); @@ -117,7 +122,7 @@ class SparseCholesky : public SparseCholesky }; template -void SparseCholesky::compute(const MatrixType& a) +void SparseLLT::compute(const MatrixType& a) { if (m_taucsSupernodalFactor) { @@ -128,7 +133,7 @@ void SparseCholesky::compute(const MatrixType& a) if (m_flags & IncompleteFactorization) { taucs_ccs_matrix taucsMatA = const_cast(a).asTaucsMatrix(); - taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0); + taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, Base::m_precision, 0); m_matrix = Base::CholMatrixType::Map(*taucsRes); free(taucsRes); m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty)) @@ -153,8 +158,8 @@ void SparseCholesky::compute(const MatrixType& a) } template -inline const typename SparseCholesky::CholMatrixType& -SparseCholesky::matrixL() const +inline const typename SparseLLT::CholMatrixType& +SparseLLT::matrixL() const { if (m_status & MatrixLIsDirty) { @@ -170,7 +175,7 @@ SparseCholesky::matrixL() const template template -void SparseCholesky::solveInPlace(MatrixBase &b) const +void SparseLLT::solveInPlace(MatrixBase &b) const { const int size = m_matrix.rows(); ei_assert(size==b.rows()); @@ -179,9 +184,9 @@ void SparseCholesky::solveInPlace(MatrixBase &b) cons { // ei_assert(!(m_status & SupernodalFactorIsDirty)); // taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b); - matrixL(); + //matrixL(); } -// else + else { Base::solveInPlace(b); } diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h index 3a1a33b3c..8e71c936d 100644 --- a/Eigen/src/Sparse/TriangularSolver.h +++ b/Eigen/src/Sparse/TriangularSolver.h @@ -102,7 +102,7 @@ struct ei_solve_triangular_selector typename Lhs::InnerIterator it(lhs, i); if(!(Lhs::Flags & UnitDiagBit)) { - std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n"; + // std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n"; ei_assert(it.index()==i); other.coeffRef(i,col) /= it.value(); } diff --git a/bench/BenchSparseUtil.h b/bench/BenchSparseUtil.h index 2c24c29e6..35c9a5263 100644 --- a/bench/BenchSparseUtil.h +++ b/bench/BenchSparseUtil.h @@ -16,7 +16,7 @@ USING_PART_OF_NAMESPACE_EIGEN #endif #ifndef SCALAR -#define SCALAR float +#define SCALAR double #endif typedef SCALAR Scalar; diff --git a/bench/benchCholesky.cpp b/bench/benchCholesky.cpp index f64b61b71..e998d8536 100644 --- a/bench/benchCholesky.cpp +++ b/bench/benchCholesky.cpp @@ -1,5 +1,5 @@ -// g++ -DNDEBUG -O3 -I.. benchCholesky.cpp -o benchCholesky && ./benchCholesky +// g++ -DNDEBUG -O3 -I.. benchLLT.cpp -o benchLLT && ./benchLLT // options: // -DBENCH_GSL -lgsl /usr/lib/libcblas.so.3 // -DEIGEN_DONT_VECTORIZE @@ -9,7 +9,7 @@ // -DSCALAR=double #include -#include +#include #include using namespace Eigen; @@ -24,7 +24,7 @@ using namespace Eigen; typedef float Scalar; template -__attribute__ ((noinline)) void benchCholesky(const MatrixType& m) +__attribute__ ((noinline)) void benchLLT(const MatrixType& m) { int rows = m.rows(); int cols = m.cols(); @@ -54,7 +54,7 @@ __attribute__ ((noinline)) void benchCholesky(const MatrixType& m) timerNoSqrt.start(); for (int k=0; k cholnosqrt(covMat); + LDLT cholnosqrt(covMat); acc += cholnosqrt.matrixL().coeff(r,c); } timerNoSqrt.stop(); @@ -65,7 +65,7 @@ __attribute__ ((noinline)) void benchCholesky(const MatrixType& m) timerSqrt.start(); for (int k=0; k chol(covMat); + LLT chol(covMat); acc += chol.matrixL().coeff(r,c); } timerSqrt.stop(); @@ -124,17 +124,17 @@ int main(int argc, char* argv[]) std::cout << "\n"; for (uint i=0; dynsizes[i]>0; ++i) - benchCholesky(Matrix(dynsizes[i],dynsizes[i])); + benchLLT(Matrix(dynsizes[i],dynsizes[i])); -// benchCholesky(Matrix()); -// benchCholesky(Matrix()); -// benchCholesky(Matrix()); -// benchCholesky(Matrix()); -// benchCholesky(Matrix()); -// benchCholesky(Matrix()); -// benchCholesky(Matrix()); -// benchCholesky(Matrix()); -// benchCholesky(Matrix()); +// benchLLT(Matrix()); +// benchLLT(Matrix()); +// benchLLT(Matrix()); +// benchLLT(Matrix()); +// benchLLT(Matrix()); +// benchLLT(Matrix()); +// benchLLT(Matrix()); +// benchLLT(Matrix()); +// benchLLT(Matrix()); return 0; } diff --git a/bench/sparse_product.cpp b/bench/sparse_product.cpp index 846301fa5..edeb08c5d 100644 --- a/bench/sparse_product.cpp +++ b/bench/sparse_product.cpp @@ -21,6 +21,18 @@ #define MINDENSITY 0.0004 #endif +#ifndef NBTRIES +#define NBTRIES 10 +#endif + +#define BENCH(X) \ + timer.reset(); \ + for (int _j=0; _j DataMatrix; // let's generate some samples on the 3D plane of equation z = 2x+3y (with some noise) DataMatrix samples = DataMatrix::Random(12,2); VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1; -// and let's solve samples * x = elevations in least square sense: -cout << (samples.adjoint() * samples).cholesky().solve((samples.adjoint()*elevations).eval()) << endl; +// and let's solve samples * [x y]^T = elevations in least square sense: +Matrix xy; +(samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations), &xy); +cout << xy << endl; diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 80614346c..f7eb7800e 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -33,7 +33,7 @@ template void cholesky(const MatrixType& m) { /* this test covers the following files: - Cholesky.h CholeskyWithoutSquareRoot.h + LLT.h LDLT.h */ int rows = m.rows(); int cols = m.cols(); @@ -44,8 +44,8 @@ template void cholesky(const MatrixType& m) typedef Matrix VectorType; MatrixType a0 = MatrixType::Random(rows,cols); - VectorType vecB = VectorType::Random(rows); - MatrixType matB = MatrixType::Random(rows,cols); + VectorType vecB = VectorType::Random(rows), vecX(rows); + MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols); SquareMatrixType symm = a0 * a0.adjoint(); // let's make sure the matrix is not singular or near singular MatrixType a1 = MatrixType::Random(rows,cols); @@ -80,28 +80,32 @@ template void cholesky(const MatrixType& m) #endif { - CholeskyWithoutSquareRoot cholnosqrt(symm); - VERIFY(cholnosqrt.isPositiveDefinite()); - VERIFY_IS_APPROX(symm, cholnosqrt.matrixL() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixL().adjoint()); - VERIFY_IS_APPROX(symm * cholnosqrt.solve(vecB), vecB); - VERIFY_IS_APPROX(symm * cholnosqrt.solve(matB), matB); + LDLT ldlt(symm); + VERIFY(ldlt.isPositiveDefinite()); + VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint()); + ldlt.solve(vecB, &vecX); + VERIFY_IS_APPROX(symm * vecX, vecB); + ldlt.solve(matB, &matX); + VERIFY_IS_APPROX(symm * matX, matB); } { - Cholesky chol(symm); + LLT chol(symm); VERIFY(chol.isPositiveDefinite()); VERIFY_IS_APPROX(symm, chol.matrixL() * chol.matrixL().adjoint()); - VERIFY_IS_APPROX(symm * chol.solve(vecB), vecB); - VERIFY_IS_APPROX(symm * chol.solve(matB), matB); + chol.solve(vecB, &vecX); + VERIFY_IS_APPROX(symm * vecX, vecB); + chol.solve(matB, &matX); + VERIFY_IS_APPROX(symm * matX, matB); } // test isPositiveDefinite on non definite matrix if (rows>4) { SquareMatrixType symm = a0.block(0,0,rows,cols-4) * a0.block(0,0,rows,cols-4).adjoint(); - Cholesky chol(symm); + LLT chol(symm); VERIFY(!chol.isPositiveDefinite()); - CholeskyWithoutSquareRoot cholnosqrt(symm); + LDLT cholnosqrt(symm); VERIFY(!cholnosqrt.isPositiveDefinite()); } } diff --git a/test/sparse.cpp b/test/sparse.cpp index 39ea05b8b..f54835972 100644 --- a/test/sparse.cpp +++ b/test/sparse.cpp @@ -217,7 +217,7 @@ template void sparse(int rows, int cols) // TODO test row major } - // test Cholesky + // test LLT { }