From b2b1a964235ea9d15e3b547660ce0a38eb88ddda Mon Sep 17 00:00:00 2001
From: nojhan <nojhan@gmail.com>
Date: Sat, 12 Nov 2011 15:52:18 +0100
Subject: [PATCH] working robust cholesky factorization, with test binary

---
 edo/src/edoSamplerNormalMulti.h | 48 +++++++++-------
 edo/test/CMakeLists.txt         |  1 +
 edo/test/t-cholesky.cpp         | 97 +++++++++++++++++++++++++++++++++
 3 files changed, 126 insertions(+), 20 deletions(-)
 create mode 100644 edo/test/t-cholesky.cpp

diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h
index c16c77ef..997b29c6 100644
--- a/edo/src/edoSamplerNormalMulti.h
+++ b/edo/src/edoSamplerNormalMulti.h
@@ -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
diff --git a/edo/test/CMakeLists.txt b/edo/test/CMakeLists.txt
index 61fdd68f..1ceea4f3 100644
--- a/edo/test/CMakeLists.txt
+++ b/edo/test/CMakeLists.txt
@@ -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
diff --git a/edo/test/t-cholesky.cpp b/edo/test/t-cholesky.cpp
new file mode 100644
index 00000000..11bcaa46
--- /dev/null
+++ b/edo/test/t-cholesky.cpp
@@ -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;
+    
+}