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

massive cleaning

Browse source
This commit is contained in:
nojhan 2011-12-15 21:38:53 +01:00
commit ff97af692d
2 changed files with 130 additions and 55 deletions
Download patch file Download diff file Expand all files Collapse all files

72
cholesky.h
View file

@ -21,8 +21,20 @@ Authors:
Caner Candan <caner.candan@thalesgroup.com> Caner Candan <caner.candan@thalesgroup.com>
*/ */
#include <stdexcept>
#include <boost/numeric/ublas/symmetric.hpp>
using namespace boost::numeric;
namespace cholesky { namespace cholesky {
class NotDefinitePositive : public std::runtime_error
{
public:
NotDefinitePositive( const std::string & what ) : std::runtime_error( what ) {}
};
/** Cholesky decomposition, given a matrix V, return a matrix L /** Cholesky decomposition, given a matrix V, return a matrix L
* such as V = L L^T (L^T being the transposed of L). * such as V = L L^T (L^T being the transposed of L).
* *
@ -31,7 +43,7 @@ namespace cholesky {
* Thus, expect a (lower) triangular matrix. * Thus, expect a (lower) triangular matrix.
*/ */
template< typename T > template< typename T >
class CholeskyBase class Cholesky
{ {
public: public:
//! The covariance-matrix is symetric //! The covariance-matrix is symetric
@ -58,12 +70,12 @@ public:
} }
/** Compute the factorization and cache the result */ /** Compute the factorization and cache the result */
virtual void operator()( const CovarMat& V ) = 0; virtual void factorize( const CovarMat& V ) = 0;
/** Compute the factorization and return the result */ /** Compute the factorization and return the result */
virtual const FactorMat& operator()( const CovarMat& V ) const FactorMat& operator()( const CovarMat& V )
{ {
(*this)( V ); this->factorize(V);
return decomposition(); return decomposition();
} }
@ -123,19 +135,19 @@ protected:
* will raise an assert before calling the square root on a negative number. * will raise an assert before calling the square root on a negative number.
*/ */
template< typename T > template< typename T >
class CholeskyLLT : public CholeskyBase<T> class LLT : public Cholesky<T>
{ {
public: public:
virtual void operator()( const CovarMat& V ) virtual void factorize( const typename Cholesky<T>::CovarMat& V )
{ {
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) ); this->_L(0, 0) = sqrt( V(0, 0) );
// end of the column // end of the column
for ( j = 1; j < N; ++j ) { for ( j = 1; j < N; ++j ) {
_L(j, 0) = V(0, j) / _L(0, 0); this->_L(j, 0) = V(0, j) / this->_L(0, 0);
} }
// end of the matrix // end of the matrix
@ -143,32 +155,36 @@ public:
// diagonal // diagonal
double sum = 0.0; double sum = 0.0;
for ( k = 0; k < i; ++k) { for ( k = 0; k < i; ++k) {
sum += _L(i, k) * _L(i, k); sum += this->_L(i, k) * this->_L(i, k);
} }
_L(i,i) = L_i_i( V, i, sum ); this->_L(i,i) = L_i_i( V, i, sum );
for ( j = i + 1; j < N; ++j ) { // rows for ( j = i + 1; j < N; ++j ) { // rows
// one element // one element
sum = 0.0; sum = 0.0;
for ( k = 0; k < i; ++k ) { for ( k = 0; k < i; ++k ) {
sum += _L(j, k) * _L(i, k); sum += this->_L(j, k) * this->_L(i, k);
} }
_L(j, i) = (V(j, i) - sum) / _L(i, i); this->_L(j, i) = (V(j, i) - sum) / this->_L(i, i);
} // for j in ]i,N[ } // for j in ]i,N[
} // for i in [1,N[ } // for i in [1,N[
} }
/** The step of the standard LLT algorithm where round off errors may appear */ /** The step of the standard LLT algorithm where round off errors may appear */
inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const inline virtual T L_i_i( const typename Cholesky<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
{ {
// round-off errors may appear here // round-off errors may appear here
assert( V(i,i) - sum >= 0 ); if( V(i,i) - sum < 0 ) {
std::ostringstream oss;
oss << "V(" << i << "/" << V.size1() << ")=" << V(i,i) << " - sum=" << sum << "\t== " << V(i,i)-sum << " < 0 ";
throw NotDefinitePositive(oss.str());
}
return sqrt( V(i,i) - sum ); return sqrt( V(i,i) - sum );
} }
}; };
@ -181,10 +197,10 @@ public:
* hack to avoid crash. * hack to avoid crash.
*/ */
template< typename T > template< typename T >
class CholeskyLLTabs : public CholeskyLLT<T> class LLTabs : public LLT<T>
{ {
public: public:
inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const inline virtual T L_i_i( const typename Cholesky<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
{ {
/***** ugly hack *****/ /***** ugly hack *****/
return sqrt( fabs( V(i,i) - sum) ); return sqrt( fabs( V(i,i) - sum) );
@ -200,10 +216,10 @@ public:
* hack to avoid crash. * hack to avoid crash.
*/ */
template< typename T > template< typename T >
class CholeskyLLTzero : public CholeskyLLT<T> class LLTzero : public LLT<T>
{ {
public: public:
inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const inline virtual T L_i_i( const typename Cholesky<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
{ {
T Lii; T Lii;
if( V(i,i) - sum >= 0 ) { if( V(i,i) - sum >= 0 ) {
@ -224,18 +240,18 @@ public:
* The factorized matrix is (L D^1/2), because V = (L D^1/2) (L D^1/2)^T * The factorized matrix is (L D^1/2), because V = (L D^1/2) (L D^1/2)^T
*/ */
template< typename T > template< typename T >
class CholeskyLDLT : public CholeskyBase<T> class LDLT : public Cholesky<T>
{ {
public: public:
virtual void operator()( const CovarMat& V ) virtual void factorize( const typename Cholesky<T>::CovarMat& V )
{ {
// use "int" everywhere, because of the "j-1" operation // use "int" everywhere, because of the "j-1" operation
int N = assert_properties( V ); int N = assert_properties( V );
// 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 );
FactorMat L = ublas::zero_matrix<T>(N,N); typename Cholesky<T>::FactorMat L = ublas::zero_matrix<T>(N,N);
FactorMat D = ublas::zero_matrix<T>(N,N); typename Cholesky<T>::FactorMat D = ublas::zero_matrix<T>(N,N);
D(0,0) = V(0,0); D(0,0) = V(0,0);
for( int j=0; j<N; ++j ) { // each columns for( int j=0; j<N; ++j ) { // each columns
@ -257,19 +273,19 @@ public:
} // for i in rows } // for i in rows
} // for j in columns } // for j in columns
_L = root( L, D ); this->_L = root( L, D );
} }
inline FactorMat root( FactorMat& L, FactorMat& D ) inline typename Cholesky<T>::FactorMat root( typename Cholesky<T>::FactorMat& L, typename Cholesky<T>::FactorMat& D )
{ {
// now compute the final symetric matrix: _L = L D^1/2 // now compute the final symetric matrix: this->_L = L D^1/2
// remember that V = ( L D^1/2) ( L D^1/2)^T // remember that V = ( L D^1/2) ( L D^1/2)^T
// fortunately, the square root of a diagonal matrix is the square // fortunately, the square root of a diagonal matrix is the square
// root of all its elements // root of all its elements
FactorMat sqrt_D = D; typename Cholesky<T>::FactorMat sqrt_D = D;
for(int i=0; i<N; ++i) { for(unsigned int i=0; i<D.size1(); ++i) {
sqrt_D(i,i) = sqrt(D(i,i)); sqrt_D(i,i) = sqrt(D(i,i));
} }

115
test.cpp
View file

@ -1,3 +1,29 @@
/*
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 <iostream>
#include <limits>
#include <iomanip>
#include <numeric>
#include "cholesky.h" #include "cholesky.h"
@ -79,7 +105,7 @@ MT error( const MT& M1, const MT& M2 )
template<typename MT > template<typename MT >
double trigsum( const MT& M ) double trigsum( const MT& M )
{ {
double sum; double sum = 0;
for( unsigned int i=0; i<M.size1(); ++i ) { for( unsigned int i=0; i<M.size1(); ++i ) {
for( unsigned int j=i; j<M.size2(); ++j ) { // triangular browsing for( unsigned int j=i; j<M.size2(); ++j ) { // triangular browsing
sum += fabs( M(i,j) ); // absolute deviation sum += fabs( M(i,j) ); // absolute deviation
@ -96,52 +122,77 @@ double sum( const T& c )
} }
typedef double real;
int main(int argc, char** argv) int main(int argc, char** argv)
{ {
srand(time(0)); srand(time(0));
unsigned int N = 4; // size of matrix unsigned int M = 10; // sample size
unsigned int R = 1000; // nb of repetitions unsigned int N = 12; // variables number
unsigned int F = 10; // range factor
unsigned int R = 1; // nb of repetitions
if( argc >= 2 ) { if( argc >= 2 ) {
N = std::atoi(argv[1]); M = std::atoi(argv[1]);
} }
if( argc >= 3 ) { if( argc >= 3 ) {
R = std::atoi(argv[2]); N = std::atoi(argv[2]);
}
if( argc >= 4 ) {
F = std::atoi(argv[3]);
}
if( argc >= 5 ) {
R = std::atoi(argv[4]);
} }
std::cout << "Usage: t-cholesky [matrix size] [repetitions]" << std::endl; std::clog << "Usage: test [sample size] [variables number] [random range] [repetitions]" << std::endl;
std::cout << "matrix size = " << N << std::endl; std::clog << "\tsample size = " << M << std::endl;
std::cout << "repetitions = " << R << std::endl; std::clog << "\tmatrix size = " << N << std::endl;
std::clog << "\trandom range = " << F << std::endl;
std::clog << "\trepetitions = " << R << std::endl;
typedef cholesky::Cholesky::CovarMat CovarMat; typedef double real;
typedef cholesky::Cholesky::FactorMat FactorMat; typedef cholesky::Cholesky<real>::CovarMat CovarMat;
typedef cholesky::Cholesky<real>::FactorMat FactorMat;
cholesky::LLT llt(); cholesky::LLT<real> llt;
cholesky::LLTabs llta(); cholesky::LLTabs<real> llta;
cholesky::LLTzero lltz(); cholesky::LLTzero<real> lltz;
cholesky::LDLT ldlt(); cholesky::LDLT<real> ldlt;
unsigned int llt_fail = 0;
std::vector<double> s0,s1,s2,s3; std::vector<double> s0,s1,s2,s3;
for( unsigned int n=0; n<R; ++n ) { for( unsigned int n=0; n<R; ++n ) {
// a variance-covariance matrix of size N*N // a random sample matrix
CovarMat V(N,N); ublas::matrix<real> S(M,N);
for( unsigned int i=0; i<M; ++i) {
// random covariance matrix for( unsigned int j=0; j<N; ++j) {
for( unsigned int i=0; i<N; ++i) { S(i,j) = F * static_cast<real>(rand())/RAND_MAX;
V(i,i) = std::pow(rand(),2); // variance should be >= 0
for( unsigned int j=i+1; j<N; ++j) {
V(i,j) = rand();
} }
} }
FactorMat L0 = llt(V); // a variance-covariance matrix of size N*N
CovarMat V0 = ublas::prod( L0, ublas::trans(L0) ); CovarMat V = ublas::prod( ublas::trans(S), S );
s0.push_back( trigsum(error(V,V0)) ); assert( V.size1() == N && V.size2() == N );
if( R == 1 ) {
std::cout << std::endl << "Covariance matrix:" << std::endl;
std::cout << format(V) << std::endl;
}
FactorMat L0;
try {
L0 = llt(V);
CovarMat V0 = ublas::prod( L0, ublas::trans(L0) );
s0.push_back( trigsum(error(V,V0)) );
} catch( cholesky::NotDefinitePositive & error ) {
llt_fail++;
#ifndef NDEBUG
std::cout << "LLT FAILED:\t" << error.what() << std::endl;
#endif
}
FactorMat L1 = llta(V); FactorMat L1 = llta(V);
CovarMat V1 = ublas::prod( L1, ublas::trans(L1) ); CovarMat V1 = ublas::prod( L1, ublas::trans(L1) );
@ -157,8 +208,16 @@ int main(int argc, char** argv)
} }
std::cout << "Average error:" << std::endl; std::cout << "Average error:" << std::endl;
std::cout << "\tLLT: " << sum(s0)/R << std::endl;
std::cout << "\tLLT: ";
if( s0.size() == 0 ) {
std::cout << "NAN";
} else {
std::cout << sum(s0)/R;
}
std::cout << "\t" << llt_fail << "/" << R << " failures" << std::endl;
std::cout << "\tLLTa: " << sum(s1)/R << std::endl; std::cout << "\tLLTa: " << sum(s1)/R << std::endl;
std::cout << "\tLLTz: " << sum(s2)/R << std::endl; std::cout << "\tLLTz: " << sum(s2)/R << std::endl;
std::cout << "\tLDLT: " << sum(s3)/R << std::endl; std::cout << "\tLDLT: " << sum(s3)/R << std::endl;
}