diff --git a/cholesky.h b/cholesky.h index 3573cdb..6765062 100644 --- a/cholesky.h +++ b/cholesky.h @@ -21,7 +21,9 @@ Authors: Caner Candan */ +#include #include +#include #include using namespace boost::numeric; @@ -277,11 +279,12 @@ public: } + /** Compute the final symetric matrix: _L = L D^1/2 + * remember that V = ( L D^1/2) ( L D^1/2)^T + * the factorization is thus L*D^1/2 + */ inline typename Cholesky::FactorMat root( typename Cholesky::FactorMat& L, typename Cholesky::FactorMat& D ) { - // now compute the final symetric matrix: this->_L = L D^1/2 - // remember that V = ( L D^1/2) ( L D^1/2)^T - // fortunately, the square root of a diagonal matrix is the square // root of all its elements typename Cholesky::FactorMat sqrt_D = D; @@ -289,7 +292,6 @@ public: sqrt_D(i,i) = sqrt(D(i,i)); } - // the factorization is thus L*D^1/2 return ublas::prod( L, sqrt_D ); } }; diff --git a/test.cpp b/test.cpp index b903adb..a8c3684 100644 --- a/test.cpp +++ b/test.cpp @@ -161,6 +161,8 @@ void test( unsigned int M, unsigned int N, unsigned int F, unsigned int R, unsig } // a variance-covariance matrix of size N*N + // Note: a covariance matrix is necessarily semi-positive definite + // thus, any failure in the Cholesky factorization is due to round-off errors CovarMat V = ublas::prod( ublas::trans(S), S ); assert( V.size1() == N && V.size2() == N ); #ifndef NDEBUG