positive_definite_constraint.hpp
Go to the documentation of this file.
1 
12 #ifndef MLPACK_METHODS_GMM_POSITIVE_DEFINITE_CONSTRAINT_HPP
13 #define MLPACK_METHODS_GMM_POSITIVE_DEFINITE_CONSTRAINT_HPP
14 
15 #include <mlpack/prereqs.hpp>
16 
17 namespace mlpack {
18 namespace gmm {
19 
28 {
29  public:
36  static void ApplyConstraint(arma::mat& covariance)
37  {
38  // What we want to do is make sure that the matrix is positive definite and
39  // that the condition number isn't too large. We also need to ensure that
40  // the covariance matrix is not too close to zero (hence, we ensure that all
41  // eigenvalues are at least 1e-50).
42  arma::vec eigval;
43  arma::mat eigvec;
44  covariance = arma::symmatu(covariance);
45  if (!arma::eig_sym(eigval, eigvec, covariance))
46  {
47  Log::Fatal << "applying to constraint could not be accomplished."
48  << std::endl;
49  }
50 
51  // If the matrix is not positive definite or if the condition number is
52  // large, we must project it back onto the cone of positive definite
53  // matrices with reasonable condition number (I'm picking 1e5 here, not for
54  // any particular reason).
55  if ((eigval[0] < 0.0) || ((eigval[eigval.n_elem - 1] / eigval[0]) > 1e5) ||
56  (eigval[eigval.n_elem - 1] < 1e-50))
57  {
58  // Project any negative eigenvalues back to non-negative, and project any
59  // too-small eigenvalues to a large enough value. Make them as small as
60  // possible to satisfy our constraint on the condition number.
61  const double minEigval = std::max(eigval[eigval.n_elem - 1] / 1e5, 1e-50);
62  for (size_t i = 0; i < eigval.n_elem; ++i)
63  eigval[i] = std::max(eigval[i], minEigval);
64 
65  // Now reassemble the covariance matrix.
66  covariance = eigvec * arma::diagmat(eigval) * eigvec.t();
67  }
68  }
69 
75  static void ApplyConstraint(arma::vec& diagCovariance)
76  {
77  // If the matrix is not positive definite or if the condition number is
78  // large, we must project it back onto the cone of positive definite
79  // matrices with reasonable condition number (I'm picking 1e5 here, not for
80  // any particular reason).
81  double maxEigval = -DBL_MAX;
82  for (size_t i = 0; i < diagCovariance.n_elem; ++i)
83  {
84  if (diagCovariance[i] > maxEigval)
85  maxEigval = diagCovariance[i];
86  }
87 
88  for (size_t i = 0; i < diagCovariance.n_elem; ++i)
89  {
90  if ((diagCovariance[i] < 0.0) || ((maxEigval / diagCovariance[i]) > 1e5)
91  || (maxEigval < 1e-50))
92  {
93  diagCovariance[i] = std::max(maxEigval / 1e5, 1e-50);
94  }
95  }
96  }
97 
99  template<typename Archive>
100  static void serialize(Archive& /* ar */, const uint32_t /* version */) { }
101 };
102 
103 } // namespace gmm
104 } // namespace mlpack
105 
106 #endif
107 
Linear algebra utility functions, generally performed on matrices or vectors.
static void ApplyConstraint(arma::mat &covariance)
Apply the positive definiteness constraint to the given covariance matrix, and ensure each value on t...
static void ApplyConstraint(arma::vec &diagCovariance)
Apply the positive definiteness constraint to the given diagonal covariance matrix (which is represen...
The core includes that mlpack expects; standard C++ includes, Armadillo, cereal, and a few basic mlpa...
constexpr T const & max(T const &lhs, T const &rhs)
Definition: algorithm.hpp:79
static void serialize(Archive &, const uint32_t)
Serialize the constraint (which stores nothing, so, nothing to do).
static util::PrefixedOutStream Fatal
Definition: log.hpp:106
Given a covariance matrix, force the matrix to be positive definite.