Implement wrapper for diagonally-constrained GMM HMMs.

CMakeLists.txt

@@ -412,8 +412,8 @@ set(MLPACK_LIBRARY_DIRS ${MLPACK_LIBRARY_DIRS} ${Boost_LIBRARY_DIRS})
 add_definitions(-DBOOST_TEST_DYN_LINK)
 
 # Detect OpenMP support in a compiler. If the compiler supports OpenMP, flags
-# to compile with OpenMP are returned and added and the HAS_OPENMP definition is
-# added for compilation.
+# to compile with OpenMP are returned and added and the HAS_OPENMP definition
+# is added for compilation.
 #
 # This way we can skip calls to functions defined in omp.h with code like:
 # #ifdef HAS_OPENMP

HISTORY.md

@@ -1,5 +1,8 @@
 ### mlpack 3.1.0
 ###### ????-??-??
+  * Add DiagonalGaussianDistribution and DiagonalGMM classes to speed up the
+    diagonal covariance computation and deprecate DiagonalConstraint (#1666).
+
   * Add kernel density estimation (KDE) implementation with bindings to other
     languages (#1301).
 

src/mlpack/core.hpp

@@ -281,6 +281,7 @@
 #include <mlpack/core/dists/gaussian_distribution.hpp>
 #include <mlpack/core/dists/laplace_distribution.hpp>
 #include <mlpack/core/dists/gamma_distribution.hpp>
+#include <mlpack/core/dists/diagonal_gaussian_distribution.hpp>
 
 // mlpack::backtrace only for linux
 #ifdef HAS_BFD_DL

src/mlpack/core/dists/CMakeLists.txt

@@ -11,6 +11,8 @@ set(SOURCES
   regression_distribution.cpp
   gamma_distribution.hpp
   gamma_distribution.cpp
+  diagonal_gaussian_distribution.hpp
+  diagonal_gaussian_distribution.cpp
 )
 
 # add directory name to sources

src/mlpack/core/dists/diagonal_gaussian_distribution.cpp

@@ -0,0 +1,148 @@
+/**
+ * @file diagonal_gaussian_distribution.cpp
+ * @author Kim SangYeon
+ *
+ * Implementation of Gaussian distribution class with diagonal covariance.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#include "diagonal_gaussian_distribution.hpp"
+#include <mlpack/methods/gmm/diagonal_constraint.hpp>
+
+using namespace mlpack;
+using namespace mlpack::distribution;
+
+DiagonalGaussianDistribution::DiagonalGaussianDistribution(
+    const arma::vec& mean,
+    const arma::vec& covariance) :
+    mean(mean)
+{
+  Covariance(covariance);
+}
+
+void DiagonalGaussianDistribution::Covariance(const arma::vec& covariance)
+{
+  this->invCov = 1 / covariance;
+  this->logDetCov = arma::accu(log(covariance));
+  this->covariance = covariance;
+}
+
+void DiagonalGaussianDistribution::Covariance(arma::vec&& covariance)
+{
+  this->invCov = 1 / covariance;
+  this->logDetCov = arma::accu(log(covariance));
+  this->covariance = std::move(covariance);
+}
+
+double DiagonalGaussianDistribution::LogProbability(
+    const arma::vec& observation) const
+{
+  const size_t k = observation.n_elem;
+  const arma::vec diff = observation - mean;
+  const arma::vec logExponent = diff.t() * arma::diagmat(invCov) * diff;
+  return -0.5 * k * log2pi - 0.5 * logDetCov - 0.5 * logExponent(0);
+}
+
+void DiagonalGaussianDistribution::LogProbability(
+    const arma::mat& observations,
+    arma::vec& logProbabilities) const
+{
+  const size_t k = observations.n_rows;
+
+  // Column i of 'diffs' is the difference between observations.col(i) and
+  // the mean.
+  arma::mat diffs = observations.each_col() - mean;
+
+  // Calculates log of exponent equation in multivariate gaussian
+  // distribution. We use only diagonal part for faster computation.
+  arma::vec logExponents = -0.5 * arma::trans(diffs % diffs) * invCov;
+
+  logProbabilities = -0.5 * k * log2pi - 0.5 * logDetCov + logExponents;
+}
+
+arma::vec DiagonalGaussianDistribution::Random() const
+{
+  return (arma::sqrt(covariance) % arma::randn<arma::vec>(mean.n_elem)) + mean;
+}
+
+void DiagonalGaussianDistribution::Train(const arma::mat& observations)
+{
+  if (observations.n_cols > 1)
+  {
+    covariance.zeros(observations.n_rows);
+  }
+  else
+  {
+    mean.zeros(0);
+    covariance.zeros(0);
+    return;
+  }
+
+  // Calculate and normalize the mean.
+  mean = arma::sum(observations, 1) / observations.n_cols;
+
+  // Now calculate the covariance.
+  const arma::mat diffs = observations.each_col() - mean;
+  covariance += arma::sum(diffs % diffs, 1);
+
+  // Finish estimating the covariance by normalizing, with the (1 / (n - 1))
+  // to make the estimator unbiased.
+  covariance /= (observations.n_cols - 1);
+  invCov = 1 / covariance;
+  logDetCov = arma::accu(log(covariance));
+}
+
+void DiagonalGaussianDistribution::Train(const arma::mat& observations,
+                                         const arma::vec& probabilities)
+{
+  if (observations.n_cols > 0)
+  {
+    covariance.zeros(observations.n_rows);
+  }
+  else
+  {
+    mean.zeros(0);
+    covariance.zeros(0);
+    return;
+  }
+
+  // We'll normalize the covariance with (v1 - (v2 / v1))
+  // for unbiased estimator in the weighted arithmetic mean. The v1 is the sum
+  // of the weights, and the v2 is the sum of the each weight squared.
+  // If you want to know more detailed description,
+  // please refer to https://en.wikipedia.org/wiki/Weighted_arithmetic_mean.
+  double v1 = arma::accu(probabilities);
+
+  // If their sum is 0, there is nothing in this Gaussian.
+  // At least, set the covariance so that it's invertible.
+  if (v1 == 0)
+  {
+    invCov = 1 / (covariance += 1e-50);
+    logDetCov = arma::accu(log(covariance));
+    return;
+  }
+
+  // Normalize the probabilities.
+  arma::vec normalizedProbs = probabilities / v1;
+
+  // Calculate the mean.
+  mean = observations * normalizedProbs;
+
+  // Now calculate the covariance.
+  const arma::mat diffs = observations.each_col() - mean;
+  covariance += (diffs % diffs) * normalizedProbs;
+
+  // Calculate the sum of each weight squared.
+  const double v2 = arma::accu(normalizedProbs % normalizedProbs);
+
+  // Finish estimating the covariance by normalizing, with
+  // the (1 / (v1 - (v2 / v1))) to make the estimator unbiased.
+  if (v2 != 1)
+    covariance /= (1 - v2);
+
+  invCov = 1 / covariance;
+  logDetCov = arma::accu(log(covariance));
+}

src/mlpack/core/dists/diagonal_gaussian_distribution.hpp

@@ -0,0 +1,156 @@
+/**
+ * @file diagonal_gaussian_distribution.hpp
+ * @author Kim SangYeon
+ *
+ * Implementation of the Gaussian distribution with diagonal covariance.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef MLPACK_CORE_DISTRIBUTIONS_DIAGONAL_GAUSSIAN_DISTRIBUTION_HPP
+#define MLPACK_CORE_DISTRIBUTIONS_DIAGONAL_GAUSSIAN_DISTRIBUTION_HPP
+
+#include <mlpack/prereqs.hpp>
+
+namespace mlpack {
+namespace distribution {
+
+//! A single multivariate Gaussian distribution with diagonal covariance.
+class DiagonalGaussianDistribution
+{
+ private:
+  //! Mean of the distribution.
+  arma::vec mean;
+  //! Diagonal covariance of the distribution.
+  arma::vec covariance;
+  //! Cached inverse of covariance.
+  arma::vec invCov;
+  //! Cached logdet(cov).
+  double logDetCov;
+
+  //! log(2pi)
+  static const constexpr double log2pi = 1.83787706640934533908193770912475883;
+
+ public:
+  //! Default constructor, which creates a Gaussian with zero dimension.
+  DiagonalGaussianDistribution() : logDetCov(0.0) { /* nothing to do. */ }
+
+  /**
+   * Create a Gaussian Distribution with zero mean and diagonal covariance
+   * with the given dimensionality.
+   * 
+   * @param dimension Number of dimensions.
+   */
+  DiagonalGaussianDistribution(const size_t dimension) :
+      mean(arma::zeros<arma::vec>(dimension)),
+      covariance(arma::ones<arma::vec>(dimension)),
+      invCov(arma::ones<arma::vec>(dimension)),
+      logDetCov(0)
+  { /* Nothing to do. */ }
+
+  /**
+   * Create a Gaussian distribution with the given mean and diagonal
+   * covariance.
+   * 
+   * @param mean Mean of distribution.
+   * @param covariance Covariance of distribution.
+   */
+  DiagonalGaussianDistribution(const arma::vec& mean,
+                               const arma::vec& covariance);
+
+  //! Return the dimensionality of this distribution.
+  size_t Dimensionality() const { return mean.n_elem; }
+
+  //! Return the probability of the given observation.
+  double Probability(const arma::vec& observation) const
+  {
+    return exp(LogProbability(observation));
+  }
+
+  //! Return the log probability of the given observation.
+  double LogProbability(const arma::vec& observation) const;
+
+  /**
+   * Calculate the multivariate Gaussian probability density function for each
+   * data point (column) in the given matrix.
+   * 
+   * @param x Matrix of observations.
+   * @param probabilities Output probabilities for each input observation.
+   */
+  void Probability(const arma::mat& x, arma::vec& probabilities) const
+  {
+    arma::vec logProbabilities;
+    LogProbability(x, logProbabilities);
+    probabilities = arma::exp(logProbabilities);
+  }
+
+  /**
+   * Calculate the multivariate Gaussian log probability density function for
+   * each data point (column) in the given matrix.
+   *
+   * @param observations Matrix of observations.
+   * @param probabilities Output log probabilities for each input observation.
+   */
+  void LogProbability(const arma::mat& observations,
+                      arma::vec& logProbabilities) const;
+
+  /**
+   * Return a randomly generated observation according to the probability
+   * distribution defined by this object.
+   * 
+   * @return Random observation from this Diagonal Gaussian distribution.
+   */
+  arma::vec Random() const;
+
+  /**
+   * Estimate the Gaussian distribution directly from the given observations.
+   * 
+   * @param observations Matrix of observations.
+   */
+  void Train(const arma::mat& observations);
+
+  /**
+   * Estimate the Gaussian distribution from the given observations,
+   * taking into account the probability of each observation actually being
+   * from this distribution.
+   * 
+   * @param observations Matrix of observations.
+   * @param probabilities List of probability of the each observation being
+   * from this distribution.
+   */
+  void Train(const arma::mat& observations,
+             const arma::vec& probabilities);
+
+  //! Return the mean.
+  const arma::vec& Mean() const { return mean; }
+
+  //! Return a modifiable copy of the mean.
+  arma::vec& Mean() { return mean; }
+
+  //! Return the covariance matrix.
+  const arma::vec& Covariance() const { return covariance; }
+
+  //! Set the covariance matrix.
+  void Covariance(const arma::vec& covariance);
+
+  //! Set the covariance matrix using move assignment.
+  void Covariance(arma::vec&& covariance);
+
+  //! Serialize the distribution.
+  template<typename Archive>
+  void serialize(Archive& ar, const unsigned int /* version */)
+  {
+    // We just need to serialize each of the members.
+    ar & BOOST_SERIALIZATION_NVP(mean);
+    ar & BOOST_SERIALIZATION_NVP(covariance);
+    ar & BOOST_SERIALIZATION_NVP(invCov);
+    ar & BOOST_SERIALIZATION_NVP(logDetCov);
+  }
+};
+
+} // namespace distribution
+} // namespace mlpack
+
+#endif

src/mlpack/methods/gmm/CMakeLists.txt

@@ -4,6 +4,9 @@ set(SOURCES
   gmm.hpp
   gmm.cpp
   gmm_impl.hpp
+  diagonal_gmm.hpp
+  diagonal_gmm.cpp
+  diagonal_gmm_impl.hpp
   em_fit.hpp
   em_fit_impl.hpp
   no_constraint.hpp

src/mlpack/methods/gmm/diagonal_constraint.hpp

@@ -30,6 +30,18 @@ class DiagonalConstraint
     covariance = arma::diagmat(arma::clamp(covariance.diag(), 1e-10, DBL_MAX));
   }
 
+  /**
+   * Apply the diagonal constraint to the given diagonal covariance matrix
+   * (which is represented as a vector), and ensure each value on the diagonal
+   * is at least 1e-10.
+   */
+  static void ApplyConstraint(arma::vec& diagCovariance)
+  {
+    // Although the covariance is already diagonal, clamp it to ensure each
+    // value is at least 1e-10.
+    diagCovariance = arma::clamp(diagCovariance, 1e-10, DBL_MAX);
+  }
+
   //! Serialize the constraint (which holds nothing, so, nothing to do).
   template<typename Archive>
   static void serialize(Archive& /* ar */, const unsigned int /* version */) { }