diff --git a/Eigen/src/CMakeLists.txt b/Eigen/src/CMakeLists.txt index ac05e77cd..14a447d86 100644 --- a/Eigen/src/CMakeLists.txt +++ b/Eigen/src/CMakeLists.txt @@ -1,3 +1,4 @@ ADD_SUBDIRECTORY(Core) ADD_SUBDIRECTORY(LU) ADD_SUBDIRECTORY(QR) +ADD_SUBDIRECTORY(Cholesky) diff --git a/Eigen/src/Cholesky/CMakeLists.txt b/Eigen/src/Cholesky/CMakeLists.txt new file mode 100644 index 000000000..fdde84add --- /dev/null +++ b/Eigen/src/Cholesky/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_Cholesky_SRCS "*.h") + +INSTALL(FILES + ${Eigen_Cholesky_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Cholesky + ) diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h new file mode 100644 index 000000000..db7c6e1ed --- /dev/null +++ b/Eigen/src/Cholesky/Cholesky.h @@ -0,0 +1,140 @@ +// 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_CHOLESKY_H +#define EIGEN_CHOLESKY_H + +/** \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 = U'U = LL', where U is upper triangular. + * + * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like A'A 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 + * situation like generalised eigen problem with hermitian matrices. + * + * \sa class CholeskyWithoutSquareRoot + */ +template class Cholesky +{ + public: + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix VectorType; + + Cholesky(const MatrixType& matrix) + : m_matrix(matrix.rows(), matrix.cols()) + { + compute(matrix); + } + + Triangular > > matrixU(void) const + { + return m_matrix.transpose().temporary().upper(); + } + + Triangular matrixL(void) const + { + return m_matrix.lower(); + } + + bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } + + template + typename DerivedVec::Eval solve(MatrixBase &vecB); + + /** Compute / recompute the Cholesky decomposition A = U'U = LL' of \a matrix + */ + void compute(const MatrixType& matrix); + + protected: + /** \internal + * Used to compute and store the cholesky decomposition. + * The strict upper part correspond to the coefficients of the input + * symmetric matrix, while the lower part store U'=L. + */ + MatrixType m_matrix; + bool m_isPositiveDefinite; +}; + +template +void Cholesky::compute(const MatrixType& matrix) +{ + assert(matrix.rows()==matrix.cols()); + const int size = matrix.rows(); + m_matrix = matrix; + + #if 1 + // this version looks faster for large matrices + m_isPositiveDefinite = m_matrix(0,0) > Scalar(0); + m_matrix(0,0) = ei_sqrt(m_matrix(0,0)); + m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0); + for (int j = 1; j < size; ++j) + { + Scalar tmp = m_matrix(j,j) - m_matrix.row(j).start(j).norm2(); + m_isPositiveDefinite = m_isPositiveDefinite && tmp > Scalar(0); + m_matrix(j,j) = ei_sqrt(tmp Scalar(0); + m_matrix(i,i) = ei_sqrt(m_matrix(i,i)); + if (i+1 +template +typename DerivedVec::Eval Cholesky::solve(MatrixBase &vecB) +{ + const int size = m_matrix.rows(); + ei_assert(size==vecB.size()); + + // FIXME .inverseProduct creates a temporary that is not nice since it is called twice + // add a .inverseProductInPlace ?? + return m_matrix.transpose().upper() + .inverseProduct(m_matrix.lower().inverseProduct(vecB)); +} + + +#endif // EIGEN_CHOLESKY_H diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h new file mode 100644 index 000000000..669728fc3 --- /dev/null +++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h @@ -0,0 +1,148 @@ +// 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_CHOLESKY_WITHOUT_SQUARE_ROOT_H +#define EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H + +/** \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 = U' D U = L D L', where U is upper 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. + * + * \todo what about complex matrices ? + * + * \sa class Cholesky + */ +template class CholeskyWithoutSquareRoot +{ + public: + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix VectorType; + + CholeskyWithoutSquareRoot(const MatrixType& matrix) + : m_matrix(matrix.rows(), matrix.cols()) + { + compute(matrix); + } + + Triangular > > matrixU(void) const + { + return m_matrix.transpose().temporary().upperWithUnitDiag(); + } + + Triangular matrixL(void) const + { + return m_matrix.lowerWithUnitDiag(); + } + + DiagonalCoeffs vectorD(void) const + { + return m_matrix.diagonal(); + } + + bool isPositiveDefinite(void) const + { + return m_matrix.diagonal().minCoeff() > Scalar(0); + } + + template + typename DerivedVec::Eval solve(MatrixBase &vecB); + + /** Compute the Cholesky decomposition A = U'DU = LDL' of \a matrix + */ + void compute(const MatrixType& matrix); + + protected: + /** \internal + * Used to compute and store the cholesky decomposition A = U'DU = LDL'. + * The strict upper part is used during the decomposition, the strict lower + * part correspond to the coefficients of U'=L (its diagonal is equal to 1 and + * is not stored), and the diagonal entries correspond to D. + */ + MatrixType m_matrix; +}; + +template +void CholeskyWithoutSquareRoot::compute(const MatrixType& matrix) +{ + assert(matrix.rows()==matrix.cols()); + const int size = matrix.rows(); + m_matrix = matrix; + #if 0 + for (int i = 0; i < size; ++i) + { + Scalar tmp = Scalar(1) / m_matrix(i,i); + for (int j = i+1; j < size; ++j) + { + m_matrix(j,i) = m_matrix(i,j) * tmp; + m_matrix.row(j).end(size-j) -= m_matrix(j,i) * m_matrix.row(i).end(size-j); + } + } + #else + // this version looks faster for large matrices + m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0); + for (int j = 1; j < size; ++j) + { + Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j))(0,0); + m_matrix(j,j) = tmp; + tmp = Scalar(1) / tmp; + for (int i = j+1; i < size; ++i) + { + m_matrix(j,i) = (m_matrix(j,i) - (m_matrix.row(i).start(j) * m_matrix.col(j).start(j))(0,0) ); + m_matrix(i,j) = tmp * m_matrix(j,i); + } + } + #endif +} + +/** Solve A*x = b with A symmeric positive definite using the available Cholesky decomposition. + */ +template +template +typename DerivedVec::Eval CholeskyWithoutSquareRoot::solve(MatrixBase &vecB) +{ + const int size = m_matrix.rows(); + ei_assert(size==vecB.size()); + + // FIXME .inverseProduct creates a temporary that is not nice since it is called twice + // maybe add a .inverseProductInPlace() ?? + return m_matrix.transpose().upperWithUnitDiag() + .inverseProduct( + (m_matrix.lowerWithUnitDiag() + .inverseProduct(vecB)) + .cwiseQuotient(m_matrix.diagonal()) + ); +} + + +#endif // EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 28395e6b2..4ff34f351 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -7,18 +7,19 @@ FIND_PACKAGE(Qt4 REQUIRED) INCLUDE_DIRECTORIES( ${QT_INCLUDE_DIR} ) SET(test_SRCS - triangular.cpp + cholesky.cpp main.cpp - basicstuff.cpp - linearstructure.cpp - product.cpp - adjoint.cpp - submatrices.cpp - miscmatrices.cpp - smallvectors.cpp - map.cpp - cwiseop.cpp - determinant.cpp +# basicstuff.cpp +# linearstructure.cpp +# product.cpp +# adjoint.cpp +# submatrices.cpp +# miscmatrices.cpp +# smallvectors.cpp +# map.cpp +# cwiseop.cpp +# determinant.cpp +# triangular.cpp ) QT4_AUTOMOC(${test_SRCS}) diff --git a/test/cholesky.cpp b/test/cholesky.cpp new file mode 100644 index 000000000..3d8351d8c --- /dev/null +++ b/test/cholesky.cpp @@ -0,0 +1,65 @@ +// 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 . + +#include "main.h" + +#include + +namespace Eigen { + +template void cholesky(const MatrixType& m) +{ + /* this test covers the following files: + Cholesky.h CholeskyWithoutSquareRoot.h + */ + int rows = m.rows(); + int cols = m.cols(); + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix SquareMatrixType; + typedef Matrix VectorType; + + MatrixType a = MatrixType::random(rows,cols).transpose(); + VectorType b = VectorType::random(cols); + SquareMatrixType covMat = a.transpose() * a; + + CholeskyWithoutSquareRoot cholnosqrt(covMat); + VERIFY_IS_APPROX(covMat, cholnosqrt.matrixU().transpose() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixU()); + VERIFY_IS_APPROX(covMat * cholnosqrt.solve(b), b); + + Cholesky chol(covMat); + VERIFY_IS_APPROX(covMat, chol.matrixU().transpose() * chol.matrixU()); + VERIFY_IS_APPROX(covMat * chol.solve(b), b); +} + +void EigenTest::testCholesky() +{ + for(int i = 0; i < 1; i++) { + cholesky(Matrix3f()); + cholesky(Matrix4d()); + cholesky(MatrixXd(17,17)); + } +} + +} // namespace Eigen diff --git a/test/main.h b/test/main.h index cab2cdc01..5d19c3b6b 100644 --- a/test/main.h +++ b/test/main.h @@ -204,17 +204,18 @@ class EigenTest : public QObject EigenTest(int repeat) : m_repeat(repeat) {} private slots: - void testBasicStuff(); - void testLinearStructure(); - void testProduct(); - void testAdjoint(); - void testSubmatrices(); - void testMiscMatrices(); - void testSmallVectors(); - void testMap(); - void testCwiseops(); - void testDeterminant(); - void testTriangular(); +// void testBasicStuff(); +// void testLinearStructure(); +// void testProduct(); +// void testAdjoint(); +// void testSubmatrices(); +// void testMiscMatrices(); +// void testSmallVectors(); +// void testMap(); +// void testCwiseops(); +// void testDeterminant(); +// void testTriangular(); + void testCholesky(); protected: int m_repeat; };