From 9fe6995df16e0c5f231b33051893102950557665 Mon Sep 17 00:00:00 2001 From: nojhan Date: Wed, 9 Nov 2011 15:48:07 +0100 Subject: [PATCH] refactoring of the cholesky decomposition --- edo/src/edoSamplerNormalMulti.h | 160 ++++++++++++++++++-------------- 1 file changed, 90 insertions(+), 70 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 7c08d22b..e5e73a01 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -28,12 +28,19 @@ Authors: #ifndef _edoSamplerNormalMulti_h #define _edoSamplerNormalMulti_h +#include + #include #include #include -//! edoSamplerNormalMulti< EOT > - +/** Sample points in a multi-normal law defined by a mean vector and a covariance matrix. + * + * Given M the mean vector and V the covariance matrix, of order n: + * - draw a vector T in N(0,I) (i.e. each value is drawn in a normal law with mean=0 an stddev=1) + * - compute the Cholesky decomposition L of V (i.e. such as V=LL*) + * - return X = M + LT + */ template< class EOT, typename D = edoNormalMulti< EOT > > class edoSamplerNormalMulti : public edoSampler< D > { @@ -42,135 +49,148 @@ public: edoSamplerNormalMulti( edoRepairer & repairer ) : edoSampler< D >( repairer) {} + /** Cholesky decomposition, given a matrix V, return a matrix L + * such as V = L Lt (Lt being the conjugate transpose of L). + * + * Need a symmetric and positive definite matrix as an input, which + * should be the case of a non-ill-conditionned covariance matrix. + * Thus, expect a (lower) triangular matrix. + */ class Cholesky { - public: - Cholesky( const ublas::symmetric_matrix< AtomType, ublas::lower >& V) - { - unsigned int Vl = V.size1(); + private: + //! The decomposition is a (lower) symetric matrix, just like the covariance matrix + ublas::symmetric_matrix< AtomType, ublas::lower > _L; + public: + //! The decomposition of the covariance matrix + const ublas::symmetric_matrix< AtomType, ublas::lower >& decomposition() const {return _L;} + + /** Computation is made at instanciation and then cached in a member variable, + * use decomposition() to get the result. + */ + Cholesky( const ublas::symmetric_matrix< AtomType, ublas::lower >& V ) + { + factorize( V ); + } + + /** 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 ublas::symmetric_matrix< AtomType, ublas::lower >& V ) + { + unsigned int Vl = V.size1(); // number of lines assert(Vl > 0); - unsigned int Vc = V.size2(); - + unsigned int Vc = V.size2(); // number of columns assert(Vc > 0); - assert( Vl == Vc ); + // FIXME assert definite semi-positive + + // the result goes in _L _L.resize(Vl); - unsigned int i,j,k; + return Vl; + } - // first column - i=0; + /** This standard algorithm makes use of square root and is thus subject + * to round-off errors if the covariance matrix is very ill-conditioned. + */ + void factorize( const ublas::symmetric_matrix< AtomType, ublas::lower >& V) + { + unsigned int N = assert_properties( V ); - // diagonal - j=0; + unsigned int i=0, j=0, k; _L(0, 0) = sqrt( V(0, 0) ); // end of the column - for ( j = 1; j < Vc; ++j ) - { + for ( j = 1; j < N; ++j ) { _L(j, 0) = V(0, j) / _L(0, 0); } // end of the matrix - for ( i = 1; i < Vl; ++i ) // each column - { - + for ( i = 1; i < N; ++i ) { // each column // diagonal double sum = 0.0; - - for ( k = 0; k < i; ++k) - { + for ( k = 0; k < i; ++k) { sum += _L(i, k) * _L(i, k); } - _L(i,i) = sqrt( fabs( V(i,i) - sum) ); + // round-off errors may appear here + assert( V(i,i) - sum >= 0 ); + _L(i,i) = sqrt( V(i,i) - sum ); + //_L(i,i) = sqrt( fabs( V(i,i) - sum) ); - for ( j = i + 1; j < Vl; ++j ) // rows - { + for ( j = i + 1; j < N; ++j ) { // rows // one element sum = 0.0; - - for ( k = 0; k < i; ++k ) - { + for ( k = 0; k < i; ++k ) { sum += _L(j, k) * _L(i, k); } _L(j, i) = (V(j, i) - sum) / _L(i, i); - } - } + + } // for j in ]i,N[ + } // for i in [1,N[ } - const ublas::symmetric_matrix< AtomType, ublas::lower >& get_L() const {return _L;} - private: - ublas::symmetric_matrix< AtomType, ublas::lower > _L; - }; + /** This alternative algorithm does not use square root BUT the covariance + * matrix must be invertible. + * + * Computes L and D such as V = L D Lt + */ + /* + void factorize_robust( const ublas::symmetric_matrix< AtomType, ublas::lower >& V) + { + unsigned int N = assert_properties( V ); + + unsigned int i, j, k; + ublas::symmetric_matrix< AtomType, ublas::lower > D = ublas::zero_matrix(N); + _L(0, 0) = sqrt( V(0, 0) ); + + } + */ + + + }; // class Cholesky + edoSamplerNormalMulti( edoBounder< EOT > & bounder ) : edoSampler< edoNormalMulti< EOT > >( bounder ) {} + EOT sample( edoNormalMulti< EOT >& distrib ) { unsigned int size = distrib.size(); - assert(size > 0); - - //------------------------------------------------------------- // Cholesky factorisation gererating matrix L from covariance // matrix V. - // We must use cholesky.get_L() to get the resulting matrix. + // We must use cholesky.decomposition() to get the resulting matrix. // // L = cholesky decomposition of varcovar - //------------------------------------------------------------- - Cholesky cholesky( distrib.varcovar() ); - ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.get_L(); + ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.decomposition(); - //------------------------------------------------------------- - - - //------------------------------------------------------------- // T = vector of size elements drawn in N(0,1) rng.normal(1.0) - //------------------------------------------------------------- - ublas::vector< AtomType > T( size ); + for ( unsigned int i = 0; i < size; ++i ) { + T( i ) = rng.normal(); + } - for ( unsigned int i = 0; i < size; ++i ) - { - T( i ) = rng.normal( 1.0 ); - } - - //------------------------------------------------------------- - - - //------------------------------------------------------------- - // LT = prod( L, T ) - //------------------------------------------------------------- - + // LT = L * T ublas::vector< AtomType > LT = ublas::prod( L, T ); - //------------------------------------------------------------- - - - //------------------------------------------------------------- // solution = means + LT - //------------------------------------------------------------- - ublas::vector< AtomType > mean = distrib.mean(); - ublas::vector< AtomType > ublas_solution = mean + LT; - EOT solution( size ); - std::copy( ublas_solution.begin(), ublas_solution.end(), solution.begin() ); - //------------------------------------------------------------- - return solution; } };