From fe2cebc0e871dbe35d84e696754195f81ea9384e Mon Sep 17 00:00:00 2001 From: nojhan Date: Sat, 12 Nov 2011 23:44:31 +0100 Subject: [PATCH] BUGFIX: factorized matrix are not symetric, cholesky factorization should process different types for covariance and decomposition + better format output for cholesky test --- edo/src/edoSamplerNormalMulti.h | 58 ++++++++--------- edo/test/t-cholesky.cpp | 106 ++++++++++++++++++++++++-------- 2 files changed, 110 insertions(+), 54 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index e2852447..35200c79 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -49,7 +49,7 @@ public: /** Cholesky decomposition, given a matrix V, return a matrix L - * such as V = L Lt (Lt being the conjugate transpose of L). + * such as V = L L^T (L^T being the transposed of L). * * Need a symmetric and positive definite matrix as an input, which * should be the case of a non-ill-conditionned covariance matrix. @@ -58,7 +58,8 @@ public: class Cholesky { public: - typedef ublas::symmetric_matrix< AtomType, ublas::lower > MatrixType; + typedef ublas::symmetric_matrix< AtomType, ublas::lower > CovarMat; + typedef ublas::matrix< AtomType > FactorMat; enum Method { //! use the standard algorithm, with square root @see factorize_LLT @@ -82,7 +83,8 @@ public: * * Use the standard unstable method by default. */ - Cholesky(const MatrixType& V, Cholesky::Method use = standard ) : _use(use) + Cholesky(const CovarMat& V, Cholesky::Method use = standard ) : + _use(use), _L(ublas::zero_matrix(V.size1(),V.size2())) { factorize( V ); } @@ -90,14 +92,14 @@ public: /** Compute the factorization and return the result */ - const MatrixType& operator()( const MatrixType& V ) + const FactorMat& operator()( const CovarMat& V ) { factorize( V ); return decomposition(); } //! The decomposition of the covariance matrix - const MatrixType & decomposition() const + const FactorMat & decomposition() const { return _L; } @@ -105,18 +107,18 @@ public: protected: //! The decomposition is a (lower) symetric matrix, just like the covariance matrix - MatrixType _L; + FactorMat _L; /** Assert that the covariance matrix have the required properties and returns its dimension. * * Note: if compiled with NDEBUG, will not assert anything and just return the dimension. */ - unsigned assert_properties( const MatrixType& V ) + unsigned assert_properties( const CovarMat& V ) { unsigned int Vl = V.size1(); // number of lines // the result goes in _L - _L.resize(Vl); + _L = ublas::zero_matrix(Vl,Vl); #ifndef NDEBUG assert(Vl > 0); @@ -135,7 +137,6 @@ public: * perform the cholesky factorization * check if all eigenvalues are positives * check if all of the leading principal minors are positive - * */ #endif @@ -146,7 +147,7 @@ public: /** Actually performs the factorization with the method given at * instanciation. Results are cached in _L. */ - void factorize( const MatrixType& V ) + void factorize( const CovarMat& V ) { if( _use == standard ) { factorize_LLT( V ); @@ -161,10 +162,12 @@ public: /** This standard algorithm makes use of square root and is thus subject * to round-off errors if the covariance matrix is very ill-conditioned. * + * Compute L such that V = L L^T + * * When compiled in debug mode and called on ill-conditionned matrix, * will raise an assert before calling the square root on a negative number. */ - void factorize_LLT( const MatrixType& V) + void factorize_LLT( const CovarMat& V) { unsigned int N = assert_properties( V ); @@ -210,7 +213,7 @@ public: * Be aware that this increase round-off errors, this is just a ugly * hack to avoid crash. */ - void factorize_LLT_abs( const MatrixType & V) + void factorize_LLT_abs( const CovarMat & V) { unsigned int N = assert_properties( V ); @@ -247,19 +250,21 @@ public: } - /** This alternative algorithm do not use square root. + /** This alternative algorithm do not use square root in an inner loop, + * but only for some diagonal elements of the matrix D. * - * Computes L and D such as V = L D Lt + * Computes L and D such as V = L D L^T. + * The factorized matrix is (L D^1/2), because V = (L D^1/2) (L D^1/2)^T */ - void factorize_LDLT( const MatrixType& V) + void factorize_LDLT( const CovarMat& V) { // use "int" everywhere, because of the "j-1" operation int N = assert_properties( V ); // example of an invertible matrix whose decomposition is undefined assert( V(0,0) != 0 ); - MatrixType L(N,N); - MatrixType D = ublas::zero_matrix(N,N); + FactorMat L = ublas::zero_matrix(N,N); + FactorMat D = ublas::zero_matrix(N,N); D(0,0) = V(0,0); for( int j=0; j 0); - // Cholesky factorisation gererating matrix L from covariance - // matrix V. - // We must use cholesky.decomposition() to get the resulting matrix. - // // L = cholesky decomposition of varcovar - const typename Cholesky::MatrixType& L = _cholesky( distrib.varcovar() ); + const typename Cholesky::FactorMat& L = _cholesky( distrib.varcovar() ); - // T = vector of size elements drawn in N(0,1) rng.normal(1.0) + // T = vector of size elements drawn in N(0,1) ublas::vector< AtomType > T( size ); for ( unsigned int i = 0; i < size; ++i ) { T( i ) = rng.normal(); diff --git a/edo/test/t-cholesky.cpp b/edo/test/t-cholesky.cpp index c06fd726..5716d848 100644 --- a/edo/test/t-cholesky.cpp +++ b/edo/test/t-cholesky.cpp @@ -28,6 +28,9 @@ Authors: #include #include #include +#include +#include +#include #include #include @@ -36,26 +39,73 @@ Authors: typedef eoReal< eoMinimizingFitness > EOT; typedef edoNormalMulti EOD; -std::ostream& operator<< (std::ostream& out, const ublas::symmetric_matrix< double, ublas::lower >& mat ) + +void setformat( std::ostream& out ) { + out << std::right; + out << std::setfill(' '); + out << std::setw( 5 + std::numeric_limits::digits10); + out << std::setprecision(std::numeric_limits::digits10); + out << std::setiosflags(std::ios_base::showpoint); +} + + +template +std::string format(const MT& mat ) +{ + std::ostringstream out; + setformat(out); + for( unsigned int i=0; i +T round( T val, T prec = 1.0 ) +{ + return (val > 0.0) ? + floor(val * prec + 0.5) / prec : + ceil(val * prec - 0.5) / prec ; +} + + +template +bool equal( const MT& M1, const MT& M2, double prec /* = 1/std::numeric_limits::digits10 ???*/ ) +{ + if( M1.size1() != M2.size1() || M1.size2() != M2.size2() ) { + return false; + } + + for( unsigned int i=0; i::Cholesky LLT( edoSamplerNormalMulti::Cholesky::standard ); + FactorMat L0 = LLT(V); + std::cout << "LLT" << std::endl << format(L0) << std::endl; + CovarMat V0 = ublas::prod( L0, ublas::trans(L0) ); + std::cout << "LLT covar" << std::endl << format(V0) << std::endl; + assert( equal(V0,V,precision) ); + std::cout << linesep << std::endl; + edoSamplerNormalMulti::Cholesky LLTa( edoSamplerNormalMulti::Cholesky::absolute ); - edoSamplerNormalMulti::Cholesky LDLT( edoSamplerNormalMulti::Cholesky::robust ); - - MatrixType L0 = LLT(V); - std::cout << "LLT" << std::endl << L0 << std::endl; - MatrixType V0 = ublas::prod( L0, ublas::trans(L0) ); - std::cout << "LLT covar" << std::endl << V0 << std::endl; - std::cout << "-----------------------------------------------------------" << std::endl; - - MatrixType L1 = LLTa(V); - std::cout << "LLT abs" << std::endl << L1 << std::endl; - MatrixType V1 = ublas::prod( L1, ublas::trans(L1) ); - std::cout << "LLT covar" << std::endl << V1 << std::endl; - std::cout << "-----------------------------------------------------------" << std::endl; + FactorMat L1 = LLTa(V); + std::cout << "LLT abs" << std::endl << format(L1) << std::endl; + CovarMat V1 = ublas::prod( L1, ublas::trans(L1) ); + std::cout << "LLT covar" << std::endl << format(V1) << std::endl; + assert( equal(V1,V,precision) ); + std::cout << linesep << std::endl; - MatrixType L2 = LDLT(V); - std::cout << "LDLT: L" << std::endl << L2 << std::endl; - MatrixType V2 = ublas::prod( L2, ublas::trans(L2) ); - std::cout << "LDLT covar" << std::endl << V2 << std::endl; - std::cout << "-----------------------------------------------------------" << std::endl; + edoSamplerNormalMulti::Cholesky LDLT( edoSamplerNormalMulti::Cholesky::robust ); + FactorMat L2 = LDLT(V); + std::cout << "LDLT" << std::endl << format(L2) << std::endl; + CovarMat V2 = ublas::prod( L2, ublas::trans(L2) ); + std::cout << "LDLT covar" << std::endl << format(V2) << std::endl; + assert( equal(V2,V,precision) ); + std::cout << linesep << std::endl; }