nojhan/eodev
Archived
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

the complete robust cholesky factorization -- bugfix in LDLT: return the factorization LD^1/2 instead of L, thus no needs of accessors on D

Browse source
This commit is contained in:
nojhan 2011-11-12 20:13:18 +01:00
commit e0f691c148
2 changed files with 33 additions and 29 deletions
Download patch file Download diff file Expand all files Collapse all files

54
edo/src/edoSamplerNormalMulti.h
View file

@ -41,8 +41,8 @@ Authors:
* - compute the Cholesky decomposition L of V (i.e. such as V=LL*) * - compute the Cholesky decomposition L of V (i.e. such as V=LL*)
* - return X = M + LT * - return X = M + LT
*/ */
template< class EOT, typename D = edoNormalMulti< EOT > > template< class EOT, typename EOD = edoNormalMulti< EOT > >
class edoSamplerNormalMulti : public edoSampler< D > class edoSamplerNormalMulti : public edoSampler< EOD >
{ {
public: public:
typedef typename EOT::AtomType AtomType; typedef typename EOT::AtomType AtomType;
@ -97,23 +97,16 @@ public:
} }
//! The decomposition of the covariance matrix //! The decomposition of the covariance matrix
const MatrixType & decomposition() const {return _L;} const MatrixType & decomposition() const
{
/** When your using the LDLT robust decomposition (by passing the "robust" return _L;
* option to the constructor, @see factorize_LDTL), this is the diagonal }
* matrix part.
*/
const MatrixType & diagonal() const {return _D;}
protected: protected:
//! The decomposition is a (lower) symetric matrix, just like the covariance matrix //! The decomposition is a (lower) symetric matrix, just like the covariance matrix
MatrixType _L; MatrixType _L;
//! The diagonal matrix when using the LDLT factorization
MatrixType _D;
/** Assert that the covariance matrix have the required properties and returns its dimension. /** 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. * Note: if compiled with NDEBUG, will not assert anything and just return the dimension.
@ -222,6 +215,7 @@ public:
unsigned int N = assert_properties( V ); unsigned int N = assert_properties( V );
unsigned int i=0, j=0, k; unsigned int i=0, j=0, k;
_L(0, 0) = sqrt( V(0, 0) ); _L(0, 0) = sqrt( V(0, 0) );
// end of the column // end of the column
@ -264,27 +258,41 @@ public:
// example of an invertible matrix whose decomposition is undefined // example of an invertible matrix whose decomposition is undefined
assert( V(0,0) != 0 ); assert( V(0,0) != 0 );
_D = ublas::zero_matrix<AtomType>(N,N); MatrixType L(N,N);
_D(0,0) = V(0,0); MatrixType D = ublas::zero_matrix<AtomType>(N,N);
D(0,0) = V(0,0);
for( int j=0; j<N; ++j ) { // each columns for( int j=0; j<N; ++j ) { // each columns
_L(j, j) = 1; L(j, j) = 1;
_D(j,j) = V(j,j); D(j,j) = V(j,j);
for( int k=0; k<=j-1; ++k) { // sum for( int k=0; k<=j-1; ++k) { // sum
_D(j,j) -= _L(j,k) * _L(j,k) * _D(k,k); D(j,j) -= L(j,k) * L(j,k) * D(k,k);
} }
for( int i=j+1; i<N; ++i ) { // remaining rows for( int i=j+1; i<N; ++i ) { // remaining rows
_L(i,j) = V(i,j); L(i,j) = V(i,j);
for( int k=0; k<=j-1; ++k) { // sum for( int k=0; k<=j-1; ++k) { // sum
_L(i,j) -= _L(i,k)*_L(j,k) * _D(k,k); L(i,j) -= L(i,k)*L(j,k) * D(k,k);
} }
_L(i,j) /= _D(j,j); L(i,j) /= D(j,j);
} // for i in rows } // for i in rows
} // for j in columns } // for j in columns
// now compute the final symetric matrix: from _LD_LT to _L_LT
// remember that V = (_LD^1/2)(_LD^1/2)^T
// square root of a diagonal matrix is the square root of all its
// scalars
MatrixType D12 = D;
for(int i=0; i<N; ++i) {
D12(i,i) = sqrt(D(i,i));
}
// the factorization is thus _L*D^1/2
_L = ublas::prod( L, D12);
} }
@ -292,11 +300,11 @@ public:
edoSamplerNormalMulti( edoRepairer<EOT> & repairer, typename Cholesky::Method use = Cholesky::absolute ) edoSamplerNormalMulti( edoRepairer<EOT> & repairer, typename Cholesky::Method use = Cholesky::absolute )
: edoSampler< D >( repairer), _cholesky(use) : edoSampler< EOD >( repairer), _cholesky(use)
{} {}
EOT sample( edoNormalMulti< EOT >& distrib ) EOT sample( EOD& distrib )
{ {
unsigned int size = distrib.size(); unsigned int size = distrib.size();
assert(size > 0); assert(size > 0);

8
edo/test/t-cholesky.cpp
View file

@ -85,12 +85,8 @@ int main(int argc, char** argv)
std::cout << "-----------------------------------------------------------" << std::endl; std::cout << "-----------------------------------------------------------" << std::endl;
MatrixType L2 = LDLT(V); MatrixType L2 = LDLT(V);
MatrixType D2 = LDLT.diagonal(); std::cout << "LDLT: L" << std::endl << L2 << std::endl;
std::cout << "LDLT" << std::endl << L2 << std::endl; MatrixType V2 = ublas::prod( L2, ublas::trans(L2) );
// ublas do not allow nested products, we should use a temporary matrix,
// thus the inline instanciation of a MatrixType
// see: http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?Effective_UBLAS
MatrixType V2 = ublas::prod( MatrixType(ublas::prod( L2, D2 )), ublas::trans(L2) );
std::cout << "LDLT covar" << std::endl << V2 << std::endl; std::cout << "LDLT covar" << std::endl << V2 << std::endl;
std::cout << "-----------------------------------------------------------" << std::endl; std::cout << "-----------------------------------------------------------" << std::endl;