From 7c98c04412322e56b3b6f7e235bc7ebb61ab6b43 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 24 Feb 2010 19:16:10 +0100 Subject: [PATCH] add reconstructedMatrix() to LLT, and LUs => they show that some improvements have still to be done for permutations, tr*tr, trapezoidal matrices --- Eigen/src/Cholesky/LDLT.h | 4 ++-- Eigen/src/Cholesky/LLT.h | 12 ++++++++++++ Eigen/src/LU/FullPivLU.h | 29 +++++++++++++++++++++++++++++ Eigen/src/LU/PartialPivLU.h | 20 ++++++++++++++++++++ test/cholesky.cpp | 7 +++---- test/lu.cpp | 21 +++++++++++++++++++++ 6 files changed, 87 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 8cfc256bb..8699fe7e0 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -155,7 +155,7 @@ template class LDLT return m_matrix; } - const MatrixType reconstructedMatrix() const; + MatrixType reconstructedMatrix() const; inline int rows() const { return m_matrix.rows(); } inline int cols() const { return m_matrix.cols(); } @@ -324,7 +324,7 @@ bool LDLT::solveInPlace(MatrixBase &bAndX) const * i.e., it returns the product: P^T L D L^* P. * This function is provided for debug purpose. */ template -const MatrixType LDLT::reconstructedMatrix() const +MatrixType LDLT::reconstructedMatrix() const { ei_assert(m_isInitialized && "LDLT is not initialized."); const int size = m_matrix.rows(); diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 96e1e5f73..2e8df7661 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -133,6 +133,8 @@ template class LLT return m_matrix; } + MatrixType reconstructedMatrix() const; + inline int rows() const { return m_matrix.rows(); } inline int cols() const { return m_matrix.cols(); } @@ -295,6 +297,16 @@ bool LLT::solveInPlace(MatrixBase &bAndX) const return true; } +/** \returns the matrix represented by the decomposition, + * i.e., it returns the product: L L^*. + * This function is provided for debug purpose. */ +template +MatrixType LLT::reconstructedMatrix() const +{ + ei_assert(m_isInitialized && "LLT is not initialized."); + return matrixL() * matrixL().adjoint().toDenseMatrix(); +} + /** \cholesky_module * \returns the LLT decomposition of \c *this */ diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 0a305d52b..cd63b9ec7 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -361,6 +361,8 @@ template class FullPivLU (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols())); } + MatrixType reconstructedMatrix() const; + inline int rows() const { return m_lu.rows(); } inline int cols() const { return m_lu.cols(); } @@ -487,6 +489,33 @@ typename ei_traits::Scalar FullPivLU::determinant() cons return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod()); } +/** \returns the matrix represented by the decomposition, + * i.e., it returns the product: P^{-1} L U Q^{-1}. + * This function is provided for debug purpose. */ +template +MatrixType FullPivLU::reconstructedMatrix() const +{ + ei_assert(m_isInitialized && "LU is not initialized."); + const int smalldim = std::min(m_lu.rows(), m_lu.cols()); + // LU + MatrixType res(m_lu.rows(),m_lu.cols()); + // FIXME the .toDenseMatrix() should not be needed... + res = m_lu.corner(TopLeft,m_lu.rows(),smalldim) + .template triangularView().toDenseMatrix() + * m_lu.corner(TopLeft,smalldim,m_lu.cols()) + .template triangularView().toDenseMatrix(); + + // P^{-1}(LU) + // FIXME implement inplace permutation + res = (m_p.inverse() * res).eval(); + + // (P^{-1}LU)Q^{-1} + // FIXME implement inplace permutation + res = (res * m_q.inverse()).eval(); + + return res; +} + /********* Implementation of kernel() **************************************************/ template diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index ed2354d78..fcffc2458 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -165,6 +165,8 @@ template class PartialPivLU */ typename ei_traits::Scalar determinant() const; + MatrixType reconstructedMatrix() const; + inline int rows() const { return m_lu.rows(); } inline int cols() const { return m_lu.cols(); } @@ -400,6 +402,24 @@ typename ei_traits::Scalar PartialPivLU::determinant() c return Scalar(m_det_p) * m_lu.diagonal().prod(); } +/** \returns the matrix represented by the decomposition, + * i.e., it returns the product: P^{-1} L U. + * This function is provided for debug purpose. */ +template +MatrixType PartialPivLU::reconstructedMatrix() const +{ + ei_assert(m_isInitialized && "LU is not initialized."); + // LU + MatrixType res = m_lu.template triangularView().toDenseMatrix() + * m_lu.template triangularView(); + + // P^{-1}(LU) + // FIXME implement inplace permutation + res = (m_p.inverse() * res).eval(); + + return res; +} + /***** Implementation of solve() *****************************************************/ template diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 1bb808d20..a446f5d73 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -95,7 +95,7 @@ template void cholesky(const MatrixType& m) { LLT chollo(symmLo); - VERIFY_IS_APPROX(symm, chollo.matrixL().toDenseMatrix() * chollo.matrixL().adjoint().toDenseMatrix()); + VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); vecX = chollo.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); matX = chollo.solve(matB); @@ -103,7 +103,7 @@ template void cholesky(const MatrixType& m) // test the upper mode LLT cholup(symmUp); - VERIFY_IS_APPROX(symm, cholup.matrixL().toDenseMatrix() * cholup.matrixL().adjoint().toDenseMatrix()); + VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix()); vecX = cholup.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); matX = cholup.solve(matB); @@ -119,8 +119,7 @@ template void cholesky(const MatrixType& m) { LDLT ldlt(symm); - // TODO(keir): This doesn't make sense now that LDLT pivots. - //VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint()); + VERIFY_IS_APPROX(symm, ldlt.reconstructedMatrix()); vecX = ldlt.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); matX = ldlt.solve(matB); diff --git a/test/lu.cpp b/test/lu.cpp index 442202a33..1ed38cb2b 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -91,6 +91,7 @@ template void lu_non_invertible() KernelMatrixType m1kernel = lu.kernel(); ImageMatrixType m1image = lu.image(m1); + VERIFY_IS_APPROX(m1, lu.reconstructedMatrix()); VERIFY(rank == lu.rank()); VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); VERIFY(!lu.isInjective()); @@ -125,6 +126,7 @@ template void lu_invertible() lu.compute(m1); } while(!lu.isInvertible()); + VERIFY_IS_APPROX(m1, lu.reconstructedMatrix()); VERIFY(0 == lu.dimensionOfKernel()); VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector VERIFY(size == lu.rank()); @@ -138,6 +140,23 @@ template void lu_invertible() VERIFY_IS_APPROX(m2, lu.inverse()*m3); } +template void lu_partial_piv() +{ + /* this test covers the following files: + PartialPivLU.h + */ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + int rows = ei_random(1,4); + int cols = rows; + + MatrixType m1(cols, rows); + m1.setRandom(); + PartialPivLU plu(m1); + + VERIFY_IS_APPROX(m1, plu.reconstructedMatrix()); +} + template void lu_verify_assert() { MatrixType tmp; @@ -180,6 +199,7 @@ void test_lu() CALL_SUBTEST_4( lu_non_invertible() ); CALL_SUBTEST_4( lu_invertible() ); + CALL_SUBTEST_4( lu_partial_piv() ); CALL_SUBTEST_4( lu_verify_assert() ); CALL_SUBTEST_5( lu_non_invertible() ); @@ -188,6 +208,7 @@ void test_lu() CALL_SUBTEST_6( lu_non_invertible() ); CALL_SUBTEST_6( lu_invertible() ); + CALL_SUBTEST_6( lu_partial_piv() ); CALL_SUBTEST_6( lu_verify_assert() ); CALL_SUBTEST_7(( lu_non_invertible >() ));