Attachment 341 Details for Bug 608 – A tentative fix

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

[patch] A tentative fix

ldlt_positive_fix.diff (text/plain), 8.12 KB, created by Gael Guennebaud on 2013-06-05 18:10:20 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Gael Guennebaud

Created: 2013-06-05 18:10:20 UTC

Size: 8.12 KB

patch

obsolete

>diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
>--- a/Eigen/src/Cholesky/LDLT.h
>+++ b/Eigen/src/Cholesky/LDLT.h
>@@ -254,17 +254,17 @@ template<> struct ldlt_inplace<Lower>
>     typedef typename MatrixType::Index Index;
>     eigen_assert(mat.rows()==mat.cols());
>     const Index size = mat.rows();
> 
>     if (size <= 1)
>     {
>       transpositions.setIdentity();
>       if(sign)
>-        *sign = real(mat.coeff(0,0))>0 ? 1:-1;
>+        *sign = math::real(mat.coeff(0,0))>0 ? 1:-1;
>       return true;
>     }
> 
>     RealScalar cutoff(0), biggest_in_corner;
> 
>     for (Index k = 0; k < size; ++k)
>     {
>       // Find largest diagonal element
>@@ -273,52 +273,43 @@ template<> struct ldlt_inplace<Lower>
>       index_of_biggest_in_corner += k;
> 
>       if(k == 0)
>       {
>         // The biggest overall is the point of reference to which further diagonals
>         // are compared; if any diagonal is negligible compared
>         // to the largest overall, the algorithm bails.
>         cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
>-
>-        if(sign)
>-          *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
>-      }
>-      else if(sign)
>-      {
>-        // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
>-        int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
>-        if(newSign != *sign)
>-          *sign = 0;
>       }
> 
>       // Finish early if the matrix is not full rank.
>       if(biggest_in_corner < cutoff)
>       {
>         for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
>+        if(sign) *sign = 0;
>         break;
>       }
> 
>       transpositions.coeffRef(k) = index_of_biggest_in_corner;
>       if(k != index_of_biggest_in_corner)
>       {
>         // apply the transposition while taking care to consider only
>         // the lower triangular part
>         Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
>         mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
>         mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
>         std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
>         for(int i=k+1;i<index_of_biggest_in_corner;++i)
>         {
>           Scalar tmp = mat.coeffRef(i,k);
>-          mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i));
>-          mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp);
>+          mat.coeffRef(i,k) = math::conj(mat.coeffRef(index_of_biggest_in_corner,i));
>+          mat.coeffRef(index_of_biggest_in_corner,i) = math::conj(tmp);
>         }
>         if(NumTraits<Scalar>::IsComplex)
>-          mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k));
>+          mat.coeffRef(index_of_biggest_in_corner,k) = math::conj(mat.coeff(index_of_biggest_in_corner,k));
>       }
> 
>       // partition the matrix:
>       //       A00 |  -  |  -
>       // lu  = A10 | A11 |  -
>       //       A20 | A21 | A22
>       Index rs = size - k - 1;
>       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
>@@ -329,16 +320,26 @@ template<> struct ldlt_inplace<Lower>
>       {
>         temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
>         mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
>         if(rs>0)
>           A21.noalias() -= A20 * temp.head(k);
>       }
>       if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
>         A21 /= mat.coeffRef(k,k);
>+      
>+      if(sign)
>+      {
>+        // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
>+        int newSign = math::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
>+        if(k == 0)
>+          *sign = newSign;
>+        else if(*sign != newSign)
>+          *sign = 0;
>+      }
>     }
> 
>     return true;
>   }
> 
>   // Reference for the algorithm: Davis and Hager, "Multiple Rank
>   // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
>   // Trivial rearrangements of their computations (Timothy E. Holy)
>@@ -362,30 +363,30 @@ template<> struct ldlt_inplace<Lower>
>     // Apply the update
>     for (Index j = 0; j < size; j++)
>     {
>       // Check for termination due to an original decomposition of low-rank
>       if (!(isfinite)(alpha))
>         break;
> 
>       // Update the diagonal terms
>-      RealScalar dj = real(mat.coeff(j,j));
>+      RealScalar dj = math::real(mat.coeff(j,j));
>       Scalar wj = w.coeff(j);
>       RealScalar swj2 = sigma*abs2(wj);
>       RealScalar gamma = dj*alpha + swj2;
> 
>       mat.coeffRef(j,j) += swj2/alpha;
>       alpha += swj2/dj;
> 
> 
>       // Update the terms of L
>       Index rs = size-j-1;
>       w.tail(rs) -= wj * mat.col(j).tail(rs);
>       if(gamma != 0)
>-        mat.col(j).tail(rs) += (sigma*conj(wj)/gamma)*w.tail(rs);
>+        mat.col(j).tail(rs) += (sigma*math::conj(wj)/gamma)*w.tail(rs);
>     }
>     return true;
>   }
> 
>   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
>   static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
>   {
>     // Apply the permutation to the input w
>diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
>--- a/Eigen/src/Cholesky/LLT.h
>+++ b/Eigen/src/Cholesky/LLT.h
>@@ -227,17 +227,17 @@ static typename MatrixType::Index llt_ra
>     }
>   }
>   else
>   {
>     temp = vec;
>     RealScalar beta = 1;
>     for(Index j=0; j<n; ++j)
>     {
>-      RealScalar Ljj = real(mat.coeff(j,j));
>+      RealScalar Ljj = math::real(mat.coeff(j,j));
>       RealScalar dj = abs2(Ljj);
>       Scalar wj = temp.coeff(j);
>       RealScalar swj2 = sigma*abs2(wj);
>       RealScalar gamma = dj*beta + swj2;
> 
>       RealScalar x = dj + swj2/beta;
>       if (x<=RealScalar(0))
>         return j;
>@@ -246,17 +246,17 @@ static typename MatrixType::Index llt_ra
>       beta += swj2/dj;
> 
>       // Update the terms of L
>       Index rs = n-j-1;
>       if(rs)
>       {
>         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
>         if(gamma != 0)
>-          mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
>+          mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*math::conj(wj)/gamma)*temp.tail(rs);
>       }
>     }
>   }
>   return -1;
> }
> 
> template<typename Scalar> struct llt_inplace<Scalar, Lower>
> {
>@@ -272,17 +272,17 @@ template<typename Scalar> struct llt_inp
>     for(Index k = 0; k < size; ++k)
>     {
>       Index rs = size-k-1; // remaining size
> 
>       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
>       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
>       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
> 
>-      RealScalar x = real(mat.coeff(k,k));
>+      RealScalar x = math::real(mat.coeff(k,k));
>       if (k>0) x -= A10.squaredNorm();
>       if (x<=RealScalar(0))
>         return k;
>       mat.coeffRef(k,k) = x = sqrt(x);
>       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
>       if (rs>0) A21 *= RealScalar(1)/x;
>     }
>     return -1;
>diff --git a/test/cholesky.cpp b/test/cholesky.cpp
>--- a/test/cholesky.cpp
>+++ b/test/cholesky.cpp
>@@ -263,20 +263,28 @@ template<typename MatrixType> void chole
> 
> // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
> // This test checks that LDLT reports correctly that matrix is indefinite. 
> // See http://forum.kde.org/viewtopic.php?f=74&t=106942
> template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
> {
>   eigen_assert(m.rows() == 2 && m.cols() == 2);
>   MatrixType mat;
>-  mat << 1, 0, 0, -1;
>-  LDLT<MatrixType> ldlt(mat);
>-  VERIFY(!ldlt.isNegative());
>-  VERIFY(!ldlt.isPositive());
>+  {
>+    mat << 1, 0, 0, -1;
>+    LDLT<MatrixType> ldlt(mat);
>+    VERIFY(!ldlt.isNegative());
>+    VERIFY(!ldlt.isPositive());
>+  }
>+  {
>+    mat << 1, 2, 2, 1;
>+    LDLT<MatrixType> ldlt(mat);
>+    VERIFY(!ldlt.isNegative());
>+    VERIFY(!ldlt.isPositive());
>+  }
> }
> 
> template<typename MatrixType> void cholesky_verify_assert()
> {
>   MatrixType tmp;
> 
>   LLT<MatrixType> llt;
>   VERIFY_RAISES_ASSERT(llt.matrixL())

Actions: View | Diff

Attachments on bug 608: 340 | 341