[mlpack] 24/40: A better positive definite constraint strategy.

[Barak A. Pearlmutter]

This is an automated email from the git hooks/post-receive script.

bap pushed a commit to branch master
in repository mlpack.

commit 9935250ba2d5e9fde591717e763bb35a4a1cbbc2
Author: Ryan Curtin <ryan at ratml.org>
Date:   Thu Jan 28 10:23:03 2016 -0500

    A better positive definite constraint strategy.
    
    Instead of just checking if the Cholesky factorization succeeds (which it might
    and it may still be inaccurate), use the condition number of the matrix to
    prevent numerical instability.  I suspect this may be somewhat slower than
    before, though, but this is the price of stability...
---
 .../methods/gmm/positive_definite_constraint.hpp   | 68 ++++++++--------------
 1 file changed, 24 insertions(+), 44 deletions(-)

diff --git a/src/mlpack/methods/gmm/positive_definite_constraint.hpp b/src/mlpack/methods/gmm/positive_definite_constraint.hpp
index d13566f..762feb2 100644
--- a/src/mlpack/methods/gmm/positive_definite_constraint.hpp
+++ b/src/mlpack/methods/gmm/positive_definite_constraint.hpp
@@ -22,7 +22,9 @@ namespace gmm {
 /**
  * Given a covariance matrix, force the matrix to be positive definite.  Also
  * force a minimum value on the diagonal, so that even if the matrix is
- * invertible, it doesn't 
+ * invertible, it doesn't cause problems with Cholesky decompositions.  The
+ * forcing here is also done in order to bring the condition number of the
+ * matrix under 1e5 (10k), which should help with numerical stability.
  */
 class PositiveDefiniteConstraint
 {
@@ -35,53 +37,31 @@ class PositiveDefiniteConstraint
    */
   static void ApplyConstraint(arma::mat& covariance)
   {
-    // Make sure each diagonal element is at least 1e-50.
-    for (size_t i = 0; i < covariance.n_cols; ++i)
-      if (std::abs(covariance(i, i)) < 1e-50)
-        covariance(i, i) = 1e-50;
+    // What we want to do is make sure that the matrix is positive definite and
+    // that the condition number isn't too large.  We also need to ensure that
+    // the covariance matrix is not too close to zero (hence, we ensure that all
+    // eigenvalues are at least 1e-50).
+    arma::vec eigval;
+    arma::mat eigvec;
+    arma::eig_sym(eigval, eigvec, covariance);
 
-    // Realistically, all we care about is that we can perform a Cholesky
-    // decomposition of the matrix, so that FactorCovariance() doesn't fail
-    // later.  Therefore, that's what we'll do to check for positive
-    // definiteness...
-    //
-    // Note that other techniques like checking the determinant *could* work,
-    // but floating-point errors mean that various decompositions may start to
-    // fail when the matrix gets close to being indefinite.  This is why we test
-    // with chol() and not something else, since that's what will be used later.
-    //
-    // We also need to make sure that the errors go to nowhere, so we have to
-    // call set_stream_err2()...
-    std::ostringstream oss;
-    std::ostream& originalStream = arma::get_stream_err2();
-    arma::set_stream_err2(oss); // Thus, errors won't be displayed.
-
-    arma::mat covLower;
-    #if (ARMA_VERSION_MAJOR < 4) || \
-        ((ARMA_VERSION_MAJOR == 4) && (ARMA_VERSION_MINOR < 500))
-    if (!arma::chol(covLower, covariance))
-    #else
-    if (!arma::chol(covLower, covariance, "lower"))
-    #endif
+    // If the matrix is not positive definite or if the condition number is
+    // large, we must project it back onto the cone of positive definite
+    // matrices with reasonable condition number (I'm picking 1e5 here, not for
+    // any particular reason).
+    if ((eigval[0] < 0.0) || ((eigval[eigval.n_elem - 1] / eigval[0]) > 1e5) ||
+        (eigval[eigval.n_elem - 1] < 1e-50))
     {
-      Log::Debug << "Covariance matrix is not positive definite.  Adding "
-          << "perturbation." << std::endl;
+      // Project any negative eigenvalues back to non-negative, and project any
+      // too-small eigenvalues to a large enough value.  Make them as small as
+      // possible to satisfy our constraint on the condition number.
+      const double minEigval = std::max(eigval[eigval.n_elem - 1] / 1e5, 1e-50);
+      for (size_t i = 0; i < eigval.n_elem; ++i)
+        eigval[i] = std::max(eigval[i], minEigval);
 
-      double perturbation = 1e-15;
-      #if (ARMA_VERSION_MAJOR < 4) || \
-          ((ARMA_VERSION_MAJOR == 4) && (ARMA_VERSION_MAJOR < 500))
-      while (!arma::chol(covLower, covariance))
-      #else
-      while (!arma::chol(covLower, covariance, "lower"))
-      #endif
-      {
-        covariance.diag() += perturbation;
-        perturbation *= 10;
-      }
+      // Now reassemble the covariance matrix.
+      covariance = eigvec * arma::diagmat(eigval) * eigvec.t();
     }
-
-    // Restore the original stream state.
-    arma::set_stream_err2(originalStream);
   }
 
   //! Serialize the constraint (which stores nothing, so, nothing to do).

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/mlpack.git