cholesky factorization with rounding to zero

edo/src/edoSamplerNormalMulti.h

@@ -70,6 +70,8 @@ public:
             standard,
             //! use the algorithm using absolute value within the square root @see factorize_LLT_abs
             absolute, 
+            //! use the method that set negative square roots to zero @see factorize_LLT_zero
+            zeroing,
             //! use the robust algorithm, without square root @see factorize_LDLT
             robust 
         };
@@ -157,6 +159,8 @@ public:
                 factorize_LLT( V );
             } else if( _use == absolute ) {
                 factorize_LLT_abs( V );
+            } else if( _use == zeroing ) {
+                factorize_LLT_zero( V );
             } else if( _use == robust ) {
                 factorize_LDLT( V );
             }
@@ -254,6 +258,54 @@ public:
         }
 
 
+        /** This standard algorithm makes use of square root but do not fail
+          * if the covariance matrix is very ill-conditioned.
+          * Here, if the diagonal difference ir negative, we set it to zero.
+          *
+          * Be aware that this increase round-off errors, this is just a ugly
+          * hack to avoid crash.
+          */
+        void factorize_LLT_zero( const CovarMat & V)
+        {
+            unsigned int N = assert_properties( V );
+
+            unsigned int i=0, j=0, k;
+
+            _L(0, 0) = sqrt( V(0, 0) );
+
+            // end of the column
+            for ( j = 1; j < N; ++j ) {
+                _L(j, 0) = V(0, j) / _L(0, 0);
+            }
+
+            // end of the matrix
+            for ( i = 1; i < N; ++i ) { // each column
+                // diagonal
+                double sum = 0.0;
+                for ( k = 0; k < i; ++k) {
+                    sum += _L(i, k) * _L(i, k);
+                }
+
+                if(  V(i,i) - sum >= 0 ) {
+                    _L(i,i) = sqrt( V(i,i) - sum);
+                } else {
+                    _L(i,i) = 0;
+                }
+
+                for ( j = i + 1; j < N; ++j ) { // rows
+                    // one element
+                    sum = 0.0;
+                    for ( k = 0; k < i; ++k ) {
+                        sum += _L(j, k) * _L(i, k);
+                    }
+
+                    _L(j, i) = (V(j, i) - sum) / _L(i, i);
+
+                } // for j in ]i,N[
+            } // for i in [1,N[
+        }
+
+
         /** This alternative algorithm do not use square root in an inner loop,
          * but only for some diagonal elements of the matrix D.
          *

edo/test/t-cholesky.cpp

@@ -158,9 +158,10 @@ int main(int argc, char** argv)
 
     edoSamplerNormalMulti<EOT,EOD>::Cholesky LLT( edoSamplerNormalMulti<EOT,EOD>::Cholesky::standard );
     edoSamplerNormalMulti<EOT,EOD>::Cholesky LLTa( edoSamplerNormalMulti<EOT,EOD>::Cholesky::absolute );
+    edoSamplerNormalMulti<EOT,EOD>::Cholesky LLTz( edoSamplerNormalMulti<EOT,EOD>::Cholesky::zeroing );
     edoSamplerNormalMulti<EOT,EOD>::Cholesky LDLT( edoSamplerNormalMulti<EOT,EOD>::Cholesky::robust );
 
-    std::vector<double> s0,s1,s2;
+    std::vector<double> s0,s1,s2,s3;
     for( unsigned int n=0; n<R; ++n ) {
 
         // a variance-covariance matrix of size N*N
@@ -182,15 +183,20 @@ int main(int argc, char** argv)
         CovarMat V1 = ublas::prod( L1, ublas::trans(L1) );
         s1.push_back( trigsum(error(V,V1)) );
 
-        FactorMat L2 = LDLT(V);
+        FactorMat L2 = LLTz(V);
         CovarMat V2 = ublas::prod( L2, ublas::trans(L2) );
         s2.push_back( trigsum(error(V,V2)) );
+
+        FactorMat L3 = LDLT(V);
+        CovarMat V3 = ublas::prod( L3, ublas::trans(L3) );
+        s3.push_back( trigsum(error(V,V3)) );
     }
 
     std::cout << "Average error:" << std::endl;
     std::cout << "\tLLT:  " << sum(s0)/R << std::endl;
     std::cout << "\tLLTa: " << sum(s1)/R << std::endl;
-    std::cout << "\tLDLT: " << sum(s2)/R << std::endl;
+    std::cout << "\tLLTz: " << sum(s2)/R << std::endl;
+    std::cout << "\tLDLT: " << sum(s3)/R << std::endl;
 
 //    double precision = 1e-15;
 //    if( argc >= 4 ) {