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

working robust cholesky factorization, with test binary

Browse source
This commit is contained in:
nojhan 2011-11-12 15:52:18 +01:00
commit b2b1a96423
3 changed files with 128 additions and 22 deletions
Download patch file Download diff file Expand all files Collapse all files

48
edo/src/edoSamplerNormalMulti.h
View file

@ -84,7 +84,7 @@ public:
*/
Cholesky(const MatrixType& V, Cholesky::Method use = standard ) : _use(use)
{
factsorize( V );
factorize( V );
}
@ -99,10 +99,19 @@ public:
//! The decomposition of the covariance matrix
const MatrixType & decomposition() const {return _L;}
/** When your using the LDLT robust decomposition (by passing the "robust"
* option to the constructor, @see factorize_LDTL), this is the diagonal
* matrix part.
*/
const MatrixType & diagonal() const {return _D;}
protected:
//! The decomposition is a (lower) symetric matrix, just like the covariance matrix
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.
@ -243,37 +252,36 @@ public:
} // for i in [1,N[
}
/** This alternative algorithm does not use square root BUT the covariance
* matrix must be invertible.
/** This alternative algorithm do not use square root.
*
* Computes L and D such as V = L D Lt
*/
void factorize_LDLT( const MatrixType& V)
{
unsigned int N = assert_properties( 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 );
unsigned int i, j, k;
//MatrixType D = ublas::zero_matrix<AtomType>(N);
std::vector<AtomType> _D(N,0);
_D = ublas::zero_matrix<AtomType>(N,N);
_D(0,0) = V(0,0);
for( int j=0; j<N; ++j ) { // each columns
_L(j, j) = 1;
_D[0] = V(0,0);
_L(0, 0) = 1;
//_L(1,0) = 1/D[0] * V(1,0);
for( j=0; j<N; ++j ) { // each columns
_D[j] = V(j,j);
for( k=0; k<j-1; ++k) { // sum
_D[j] -= _L(j,k) * _L(j,k) * _D[k];
_D(j,j) = V(j,j);
for( int k=0; k<=j-1; ++k) { // sum
_D(j,j) -= _L(j,k) * _L(j,k) * _D(k,k);
}
for( i=j+1; i<N; ++i ) { // remaining rows
for( int i=j+1; i<N; ++i ) { // remaining rows
_L(i,j) = V(i,j);
for( k=0; k<j-1; ++k) { // sum
_L(i,j) -= _L(i,k)*_L(j,k) * _D[k];
for( int k=0; k<=j-1; ++k) { // sum
_L(i,j) -= _L(i,k)*_L(j,k) * _D(k,k);
}
_L(i,j) /= _D[j];
_L(i,j) /= _D(j,j);
} // for i in rows
} // for j in columns

1
edo/test/CMakeLists.txt
View file

@ -33,6 +33,7 @@ LINK_DIRECTORIES(${Boost_LIBRARY_DIRS})
INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/application/common)
SET(SOURCES
t-cholesky
t-edoEstimatorNormalMulti
t-mean-distance
t-bounderno

97
edo/test/t-cholesky.cpp Normal file
View file

@ -0,0 +1,97 @@
/*
The Evolving Distribution Objects framework (EDO) is a template-based,
ANSI-C++ evolutionary computation library which helps you to write your
own estimation of distribution algorithms.
This library 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 2.1 of the License, or (at your option) any later version.
This library 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 for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Copyright (C) 2010 Thales group
*/
/*
Authors:
Johann Dréo <johann.dreo@thalesgroup.com>
*/
#include <vector>
#include <cstdlib>
#include <iostream>
#include <eo>
#include <es.h>
#include <edo>
typedef eoReal< eoMinimizingFitness > EOT;
typedef edoNormalMulti<EOT> EOD;
std::ostream& operator<< (std::ostream& out, const ublas::symmetric_matrix< double, ublas::lower >& mat )
{
for( unsigned int i=0; i<mat.size1(); ++i) {
for( unsigned int j=0; j<=i; ++j) {
out << mat(i,j) << "\t";
} // columns
out << std::endl;
} // rows
return out;
}
int main(int argc, char** argv)
{
unsigned int N = 4;
typedef edoSamplerNormalMulti<EOT,EOD>::Cholesky::MatrixType MatrixType;
// a variance-covariance matrix of size N*N
MatrixType V(N,N);
// random covariance matrix
for( unsigned int i=0; i<N; ++i) {
V(i,i) = 1 + std::pow(rand(),2); // variance should be > 0
for( unsigned int j=i+1; j<N; ++j) {
V(i,j) = rand();
}
}
std::cout << "Covariance matrix" << std::endl << V << std::endl;
std::cout << "-----------------------------------------------------------" << std::endl;
edoSamplerNormalMulti<EOT,EOD>::Cholesky LLT( edoSamplerNormalMulti<EOT,EOD>::Cholesky::standard );
edoSamplerNormalMulti<EOT,EOD>::Cholesky LLTa( edoSamplerNormalMulti<EOT,EOD>::Cholesky::absolute );
edoSamplerNormalMulti<EOT,EOD>::Cholesky LDLT( edoSamplerNormalMulti<EOT,EOD>::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;
MatrixType L2 = LDLT(V);
MatrixType D2 = LDLT.diagonal();
std::cout << "LDLT" << std::endl << L2 << std::endl;
// 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 << "-----------------------------------------------------------" << std::endl;
}