From 6c4ab7f43ed812f71958fc3d17521e03279e9416 Mon Sep 17 00:00:00 2001 From: Stephen Tu Date: Tue, 20 Jan 2015 12:16:40 -0800 Subject: [PATCH 1/8] WIP: gaussian distribution caches covariance factorizations compiles, not tested --- .gitignore | 1 + .../core/dists/gaussian_distribution.cpp | 50 ++++++++++++++++--- .../core/dists/gaussian_distribution.hpp | 36 ++++++++----- .../core/dists/regression_distribution.hpp | 4 +- src/mlpack/methods/gmm/em_fit_impl.hpp | 47 ++++++++++------- src/mlpack/methods/gmm/gmm_convert_main.cpp | 4 +- src/mlpack/methods/hmm/hmm_util_impl.hpp | 8 ++- src/mlpack/tests/distribution_test.cpp | 32 ++++++++---- src/mlpack/tests/gmm_test.cpp | 5 +- src/mlpack/tests/hmm_test.cpp | 10 +++- 10 files changed, 142 insertions(+), 55 deletions(-) diff --git a/.gitignore b/.gitignore index 1b2211df079..726d15396af 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ build* +*.swp diff --git a/src/mlpack/core/dists/gaussian_distribution.cpp b/src/mlpack/core/dists/gaussian_distribution.cpp index cb37e79f1ea..969c33e45aa 100644 --- a/src/mlpack/core/dists/gaussian_distribution.cpp +++ b/src/mlpack/core/dists/gaussian_distribution.cpp @@ -10,21 +10,49 @@ using namespace mlpack; using namespace mlpack::distribution; + +GaussianDistribution::GaussianDistribution(const arma::vec& mean, + const arma::mat& covariance) + : mean(mean) +{ + Covariance(covariance); +} + +void GaussianDistribution::Covariance(const arma::mat& covariance) +{ + this->covariance = covariance; + FactorCovariance(); +} + +void GaussianDistribution::Covariance(arma::mat&& covariance) +{ + this->covariance = std::move(covariance); + FactorCovariance(); +} + +void GaussianDistribution::FactorCovariance() +{ + covLower = arma::chol(covariance, "lower"); + // tell arma that this is lower triangular matrix (for faster inversion) + covLower = arma::trimatl(covLower); + const arma::mat invCovLower = covLower.i(); + invCov = invCovLower.t() * invCovLower; + double sign = 0.; + arma::log_det(logDetCov, sign, covLower); + logDetCov *= 2; +} + double GaussianDistribution::LogProbability(const arma::vec& observation) const { const size_t k = observation.n_elem; - double logdetsigma = 0; - double sign = 0.; - arma::log_det(logdetsigma, sign, covariance); const arma::vec diff = mean - observation; - const arma::vec v = (diff.t() * arma::inv(covariance) * diff); - return -0.5 * k * log2pi - 0.5 * logdetsigma - 0.5 * v(0); + const arma::vec v = (diff.t() * invCov * diff); + return -0.5 * k * log2pi - 0.5 * logDetCov - 0.5 * v(0); } arma::vec GaussianDistribution::Random() const { - // Should we store chol(covariance) for easier calculation later? - return trans(chol(covariance)) * arma::randn(mean.n_elem) + mean; + return covLower * arma::randn(mean.n_elem) + mean; } /** @@ -41,6 +69,7 @@ void GaussianDistribution::Estimate(const arma::mat& observations) } else // This will end up just being empty. { + // TODO(stephentu): why do we allow this case? why not throw an error? mean.zeros(0); covariance.zeros(0); return; @@ -77,6 +106,8 @@ void GaussianDistribution::Estimate(const arma::mat& observations) perturbation *= 10; // Slow, but we don't want to add too much. } } + + FactorCovariance(); } /** @@ -94,6 +125,7 @@ void GaussianDistribution::Estimate(const arma::mat& observations, } else // This will end up just being empty. { + // TODO(stephentu): same as above mean.zeros(0); covariance.zeros(0); return; @@ -114,6 +146,7 @@ void GaussianDistribution::Estimate(const arma::mat& observations, // Nothing in this Gaussian! At least set the covariance so that it's // invertible. covariance.diag() += 1e-50; + // TODO(stephentu): why do we allow this case? return; } @@ -143,6 +176,8 @@ void GaussianDistribution::Estimate(const arma::mat& observations, perturbation *= 10; // Slow, but we don't want to add too much. } } + + FactorCovariance(); } /** @@ -180,4 +215,5 @@ void GaussianDistribution::Load(const util::SaveRestoreUtility& sr) { sr.LoadParameter(mean, "mean"); sr.LoadParameter(covariance, "covariance"); + FactorCovariance(); } diff --git a/src/mlpack/core/dists/gaussian_distribution.hpp b/src/mlpack/core/dists/gaussian_distribution.hpp index 0e7ea2482af..8462aaf34c1 100644 --- a/src/mlpack/core/dists/gaussian_distribution.hpp +++ b/src/mlpack/core/dists/gaussian_distribution.hpp @@ -21,8 +21,14 @@ class GaussianDistribution private: //! Mean of the distribution. arma::vec mean; - //! Covariance of the distribution. + //! Positive definite covariance of the distribution. arma::mat covariance; + //! Lower triangular factor of cov (e.g. cov = LL^T). + arma::mat covLower; + //! Cached inverse of covariance. + arma::mat invCov; + //! Cached logdet(cov). + double logDetCov; //! log(2pi) static const constexpr double log2pi = 1.83787706640934533908193770912475883; @@ -39,14 +45,20 @@ class GaussianDistribution */ GaussianDistribution(const size_t dimension) : mean(arma::zeros(dimension)), - covariance(arma::eye(dimension, dimension)) + covariance(arma::eye(dimension, dimension)), + covLower(arma::eye(dimension, dimension)), + invCov(arma::eye(dimension, dimension)), + logDetCov(0) { /* Nothing to do. */ } /** * Create a Gaussian distribution with the given mean and covariance. + * + * covariance is expected to be positive definite. */ - GaussianDistribution(const arma::vec& mean, const arma::mat& covariance) : - mean(mean), covariance(covariance) { /* Nothing to do. */ } + GaussianDistribution(const arma::vec& mean, const arma::mat& covariance); + + // TODO(stephentu): do we want a (arma::vec&&, arma::mat&&) ctor? //! Return the dimensionality of this distribution. size_t Dimensionality() const { return mean.n_elem; } @@ -119,9 +131,11 @@ class GaussianDistribution const arma::mat& Covariance() const { return covariance; } /** - * Return a modifiable copy of the covariance. + * Set the covariance. */ - arma::mat& Covariance() { return covariance; } + void Covariance(const arma::mat& covariance); + + void Covariance(arma::mat&& covariance); /** * Returns a string representation of this object. @@ -135,7 +149,8 @@ class GaussianDistribution void Load(const util::SaveRestoreUtility& n); static std::string const Type() { return "GaussianDistribution"; } - + private: + void FactorCovariance(); }; @@ -156,17 +171,14 @@ inline void GaussianDistribution::LogProbability(const arma::mat& x, // diffs). We just don't need any of the other elements. We can calculate // the right hand part of the equation (instead of the left side) so that // later we are referencing columns, not rows -- that is faster. - arma::mat rhs = -0.5 * inv(covariance) * diffs; + const arma::mat rhs = -0.5 * invCov * diffs; arma::vec logExponents(diffs.n_cols); // We will now fill this. for (size_t i = 0; i < diffs.n_cols; i++) logExponents(i) = accu(diffs.unsafe_col(i) % rhs.unsafe_col(i)); - double logdetsigma = 0; - double sign = 0.; - arma::log_det(logdetsigma, sign, covariance); const size_t k = x.n_rows; - logProbabilities = -0.5 * k * log2pi - 0.5 * logdetsigma + logExponents; + logProbabilities = -0.5 * k * log2pi - 0.5 * logDetCov + logExponents; } diff --git a/src/mlpack/core/dists/regression_distribution.hpp b/src/mlpack/core/dists/regression_distribution.hpp index 6cda6812ea1..b3186a839d9 100644 --- a/src/mlpack/core/dists/regression_distribution.hpp +++ b/src/mlpack/core/dists/regression_distribution.hpp @@ -48,7 +48,9 @@ class RegressionDistribution rf(regression::LinearRegression(predictors, responses)) { err = GaussianDistribution(1); - err.Covariance() = rf.ComputeError(predictors, responses); + arma::mat cov(1, 1); + cov(0, 0) = rf.ComputeError(predictors, responses); + err.Covariance(std::move(cov)); } /** diff --git a/src/mlpack/methods/gmm/em_fit_impl.hpp b/src/mlpack/methods/gmm/em_fit_impl.hpp index c1457dd4843..dead473e3a8 100644 --- a/src/mlpack/methods/gmm/em_fit_impl.hpp +++ b/src/mlpack/methods/gmm/em_fit_impl.hpp @@ -93,11 +93,12 @@ void EMFit::Estimate( trans(condProb.col(i))); // Don't update if there's no probability of the Gaussian having points. - if (probRowSums[i] != 0.0) - dists[i].Covariance() = (tmp * trans(tmpB)) / probRowSums[i]; - - // Apply covariance constraint. - constraint.ApplyConstraint(dists[i].Covariance()); + if (probRowSums[i] != 0.0) { + arma::mat covariance = (tmp * trans(tmpB)) / probRowSums[i]; + // Apply covariance constraint. + constraint.ApplyConstraint(covariance); + dists[i].Covariance(std::move(covariance)); + } } // Calculate the new values for omega using the updated conditional @@ -180,10 +181,12 @@ void EMFit::Estimate( arma::mat tmpB = tmp % (arma::ones(observations.n_rows) * trans(condProb.col(i) % probabilities)); - dists[i].Covariance() = (tmp * trans(tmpB)) / probRowSums[i]; + arma::mat cov = (tmp * trans(tmpB)) / probRowSums[i]; // Apply covariance constraint. - constraint.ApplyConstraint(dists[i].Covariance()); + constraint.ApplyConstraint(cov); + + dists[i].Covariance(std::move(cov)); } // Calculate the new values for omega using the updated conditional @@ -210,12 +213,16 @@ InitialClustering(const arma::mat& observations, // Run clustering algorithm. clusterer.Cluster(observations, dists.size(), assignments); + std::vector means(dists.size()); + std::vector covs(dists.size()); + // Now calculate the means, covariances, and weights. weights.zeros(); for (size_t i = 0; i < dists.size(); ++i) { - dists[i].Mean().zeros(); - dists[i].Covariance().zeros(); + means[i].zeros(dists[i].Mean().n_elem); + covs[i].zeros(dists[i].Covariance().n_rows, + dists[i].Covariance().n_cols); } // From the assignments, generate our means, covariances, and weights. @@ -224,11 +231,10 @@ InitialClustering(const arma::mat& observations, const size_t cluster = assignments[i]; // Add this to the relevant mean. - dists[cluster].Mean() += observations.col(i); + means[cluster] += observations.col(i); // Add this to the relevant covariance. - dists[cluster].Covariance() += observations.col(i) * - trans(observations.col(i)); + covs[cluster] += observations.col(i) * trans(observations.col(i)); // Now add one to the weights (we will normalize). weights[cluster]++; @@ -237,22 +243,25 @@ InitialClustering(const arma::mat& observations, // Now normalize the mean and covariance. for (size_t i = 0; i < dists.size(); ++i) { - dists[i].Mean() /= (weights[i] > 1) ? weights[i] : 1; + means[i] /= (weights[i] > 1) ? weights[i] : 1; } for (size_t i = 0; i < observations.n_cols; ++i) { const size_t cluster = assignments[i]; - const arma::vec normObs = observations.col(i) - dists[cluster].Mean(); - dists[cluster].Covariance() += normObs * normObs.t(); + const arma::vec normObs = observations.col(i) - means[cluster]; + covs[cluster] += normObs * normObs.t(); } for (size_t i = 0; i < dists.size(); ++i) { - dists[i].Covariance() /= (weights[i] > 1) ? weights[i] : 1; + covs[i] /= (weights[i] > 1) ? weights[i] : 1; // Apply constraints to covariance matrix. - constraint.ApplyConstraint(dists[i].Covariance()); + constraint.ApplyConstraint(covs[i]); + + std::swap(dists[i].Mean(), means[i]); + dists[i].Covariance(std::move(covs[i])); } // Finally, normalize weights. @@ -269,7 +278,7 @@ double EMFit::LogLikelihood( arma::vec phis; arma::mat likelihoods(dists.size(), observations.n_cols); - + for (size_t i = 0; i < dists.size(); ++i) { dists[i].Probability(observations, phis); @@ -283,7 +292,7 @@ double EMFit::LogLikelihood( << "outlier." << std::endl; logLikelihood += log(accu(likelihoods.col(j))); } - + return logLikelihood; } diff --git a/src/mlpack/methods/gmm/gmm_convert_main.cpp b/src/mlpack/methods/gmm/gmm_convert_main.cpp index 7d80ebd0fc0..1f484f90311 100644 --- a/src/mlpack/methods/gmm/gmm_convert_main.cpp +++ b/src/mlpack/methods/gmm/gmm_convert_main.cpp @@ -47,12 +47,14 @@ int main(int argc, char* argv[]) for (size_t i = 0; i < gaussians; ++i) { stringstream o; + arma::mat covariance; o << i; string meanName = "mean" + o.str(); string covName = "covariance" + o.str(); load.LoadParameter(gmm.Component(i).Mean(), meanName); - load.LoadParameter(gmm.Component(i).Covariance(), covName); + load.LoadParameter(covariance, covName); + gmm.Component(i).Covariance(std::move(covariance)); } gmm.Save(CLI::GetParam("output_file")); diff --git a/src/mlpack/methods/hmm/hmm_util_impl.hpp b/src/mlpack/methods/hmm/hmm_util_impl.hpp index d3599abcfac..47fbd9de5a0 100644 --- a/src/mlpack/methods/hmm/hmm_util_impl.hpp +++ b/src/mlpack/methods/hmm/hmm_util_impl.hpp @@ -115,8 +115,10 @@ void ConvertHMM(HMM& hmm, sr.LoadParameter(hmm.Emission()[i].Mean(), s.str()); s.str(""); + arma::mat covariance; s << "hmm_emission_covariance_" << i; - sr.LoadParameter(hmm.Emission()[i].Covariance(), s.str()); + sr.LoadParameter(covariance, s.str()); + hmm.Emission()[i].Covariance(std::move(covariance)); } hmm.Dimensionality() = hmm.Emission()[0].Mean().n_elem; @@ -168,7 +170,9 @@ void ConvertHMM(HMM >& hmm, const util::SaveRestoreUtility& sr) s.str(""); s << "hmm_emission_" << i << "_gaussian_" << g << "_covariance"; - sr.LoadParameter(hmm.Emission()[i].Component(g).Covariance(), s.str()); + arma::mat covariance; + sr.LoadParameter(covariance, s.str()); + hmm.Emission()[i].Component(g).Covariance(std::move(covariance)); } s.str(""); diff --git a/src/mlpack/tests/distribution_test.cpp b/src/mlpack/tests/distribution_test.cpp index 64268fd5064..fc04c1a0403 100644 --- a/src/mlpack/tests/distribution_test.cpp +++ b/src/mlpack/tests/distribution_test.cpp @@ -187,16 +187,23 @@ BOOST_AUTO_TEST_CASE(GaussianUnivariateProbabilityTest) 1e-5); // A few more cases... - g.Covariance().fill(2.0); + arma::mat covariance; + covariance.set_size(g.Covariance().n_rows, g.Covariance().n_cols); + + covariance.fill(2.0); + g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("0.0")), 0.282094791773878, 1e-5); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("1.0")), 0.219695644733861, 1e-5); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("-1.0")), 0.219695644733861, 1e-5); g.Mean().fill(1.0); - g.Covariance().fill(1.0); + covariance.fill(1.0); + g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("1.0")), 0.398942280401433, 1e-5); - g.Covariance().fill(2.0); + + covariance.fill(2.0); + g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("-1.0")), 0.103776874355149, 1e-5); } @@ -215,7 +222,9 @@ BOOST_AUTO_TEST_CASE(GaussianMultivariateProbabilityTest) BOOST_REQUIRE_CLOSE(g.Probability(x), 0.159154943091895, 1e-5); - g.Covariance() = "2 0; 0 2"; + arma::mat covariance; + covariance = "2 0; 0 2"; + g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(x), 0.0795774715459477, 1e-5); @@ -230,7 +239,8 @@ BOOST_AUTO_TEST_CASE(GaussianMultivariateProbabilityTest) BOOST_REQUIRE_CLOSE(g.Probability(-x), 0.0795774715459477, 1e-5); g.Mean() = "1 1"; - g.Covariance() = "2 1.5; 1 4"; + covariance = "2 1.5; 1 4"; + g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(x), 0.0624257046546403, 1e-5); g.Mean() *= -1; @@ -245,11 +255,13 @@ BOOST_AUTO_TEST_CASE(GaussianMultivariateProbabilityTest) // Higher-dimensional case. x = "0 1 2 3 4"; g.Mean() = "5 6 3 3 2"; - g.Covariance() = "6 1 1 0 2;" - "0 7 1 0 1;" - "1 1 4 1 1;" - "1 0 1 7 0;" - "2 0 1 1 6"; + + covariance = "6 1 1 0 2;" + "0 7 1 0 1;" + "1 1 4 1 1;" + "1 0 1 7 0;" + "2 0 1 1 6"; + g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(x), 1.02531207499358e-6, 1e-5); BOOST_REQUIRE_CLOSE(g.Probability(-x), 1.06784794079363e-8, 1e-5); diff --git a/src/mlpack/tests/gmm_test.cpp b/src/mlpack/tests/gmm_test.cpp index e043376be03..d00a50eeec7 100644 --- a/src/mlpack/tests/gmm_test.cpp +++ b/src/mlpack/tests/gmm_test.cpp @@ -490,7 +490,10 @@ BOOST_AUTO_TEST_CASE(GMMLoadSaveTest) for (size_t i = 0; i < gmm.Gaussians(); ++i) { gmm.Component(i).Mean().randu(); - gmm.Component(i).Covariance().randu(); + arma::mat covariance = arma::randu( + gmm.Component(i).Covariance().n_rows, + gmm.Component(i).Covariance().n_cols); + gmm.Component(i).Covariance(std::move(covariance)); } gmm.Save("test-gmm-save.xml"); diff --git a/src/mlpack/tests/hmm_test.cpp b/src/mlpack/tests/hmm_test.cpp index 4368a1fc543..9ae7016a46e 100644 --- a/src/mlpack/tests/hmm_test.cpp +++ b/src/mlpack/tests/hmm_test.cpp @@ -992,7 +992,10 @@ BOOST_AUTO_TEST_CASE(GMMHMMLoadSaveTest) for (size_t i = 0; i < hmm.Emission()[j].Gaussians(); ++i) { hmm.Emission()[j].Component(i).Mean().randu(); - hmm.Emission()[j].Component(i).Covariance().randu(); + arma::mat covariance = arma::randu( + hmm.Emission()[j].Component(i).Covariance().n_rows, + hmm.Emission()[j].Component(i).Covariance().n_cols); + hmm.Emission()[j].Component(i).Covariance(std::move(covariance)); } } @@ -1048,7 +1051,10 @@ BOOST_AUTO_TEST_CASE(GaussianHMMLoadSaveTest) for(size_t j = 0; j < hmm.Emission().size(); ++j) { hmm.Emission()[j].Mean().randu(); - hmm.Emission()[j].Covariance().randu(); + arma::mat covariance = arma::randu( + hmm.Emission()[j].Covariance().n_rows, + hmm.Emission()[j].Covariance().n_cols); + hmm.Emission()[j].Covariance(std::move(covariance)); } util::SaveRestoreUtility sr; From 87d8d991ea0076b0a6f7fb2757aef9ece1f4deb2 Mon Sep 17 00:00:00 2001 From: Stephen Tu Date: Tue, 20 Jan 2015 14:10:48 -0800 Subject: [PATCH 2/8] WIP: fix one test case --- src/mlpack/tests/distribution_test.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mlpack/tests/distribution_test.cpp b/src/mlpack/tests/distribution_test.cpp index fc04c1a0403..2ff1eb0edec 100644 --- a/src/mlpack/tests/distribution_test.cpp +++ b/src/mlpack/tests/distribution_test.cpp @@ -140,6 +140,8 @@ BOOST_AUTO_TEST_CASE(GaussianDistributionDistributionConstructor) mean.randu(); covariance.randu(); + covariance *= covariance.t(); + covariance += arma::eye(3, 3); GaussianDistribution d(mean, covariance); From d30fd427aec35725c134d54381dc8f8bb17efacc Mon Sep 17 00:00:00 2001 From: Stephen Tu Date: Tue, 20 Jan 2015 16:43:27 -0800 Subject: [PATCH 3/8] fix some distribution test cases --- src/mlpack/tests/distribution_test.cpp | 40 +++++++++++++------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/mlpack/tests/distribution_test.cpp b/src/mlpack/tests/distribution_test.cpp index 2ff1eb0edec..9041bb5f8b0 100644 --- a/src/mlpack/tests/distribution_test.cpp +++ b/src/mlpack/tests/distribution_test.cpp @@ -159,20 +159,20 @@ BOOST_AUTO_TEST_CASE(GaussianDistributionDistributionConstructor) BOOST_AUTO_TEST_CASE(GaussianDistributionProbabilityTest) { arma::vec mean("5 6 3 3 2"); - arma::mat cov("6 1 1 0 2;" - "0 7 1 0 1;" + arma::mat cov("6 1 1 1 2;" + "1 7 1 0 0;" "1 1 4 1 1;" "1 0 1 7 0;" - "2 0 1 1 6"); + "2 0 1 0 6"); GaussianDistribution d(mean, cov); - BOOST_REQUIRE_CLOSE(d.Probability("0 1 2 3 4"), 1.02531207499358e-6, 1e-5); - BOOST_REQUIRE_CLOSE(d.Probability("3 2 3 7 8"), 1.82353695848039e-7, 1e-5); - BOOST_REQUIRE_CLOSE(d.Probability("2 2 0 8 1"), 1.29759261892949e-6, 1e-5); - BOOST_REQUIRE_CLOSE(d.Probability("2 1 5 0 1"), 1.33218060268258e-6, 1e-5); - BOOST_REQUIRE_CLOSE(d.Probability("3 0 5 1 0"), 1.12120427975708e-6, 1e-5); - BOOST_REQUIRE_CLOSE(d.Probability("4 0 6 1 0"), 4.57951032485297e-7, 1e-5); + BOOST_REQUIRE_CLOSE(d.LogProbability("0 1 2 3 4"), -13.432076798791542, 1e-5); + BOOST_REQUIRE_CLOSE(d.LogProbability("3 2 3 7 8"), -15.814880322345738, 1e-5); + BOOST_REQUIRE_CLOSE(d.LogProbability("2 2 0 8 1"), -13.754462857772776, 1e-5); + BOOST_REQUIRE_CLOSE(d.LogProbability("2 1 5 0 1"), -13.283283233107898, 1e-5); + BOOST_REQUIRE_CLOSE(d.LogProbability("3 0 5 1 0"), -13.800326511545279, 1e-5); + BOOST_REQUIRE_CLOSE(d.LogProbability("4 0 6 1 0"), -14.900192463287908, 1e-5); } /** @@ -244,33 +244,33 @@ BOOST_AUTO_TEST_CASE(GaussianMultivariateProbabilityTest) covariance = "2 1.5; 1 4"; g.Covariance(std::move(covariance)); - BOOST_REQUIRE_CLOSE(g.Probability(x), 0.0624257046546403, 1e-5); + BOOST_REQUIRE_CLOSE(g.Probability(x), 0.060154914192541771, 1e-5); g.Mean() *= -1; - BOOST_REQUIRE_CLOSE(g.Probability(-x), 0.0624257046546403, 1e-5); + BOOST_REQUIRE_CLOSE(g.Probability(-x), 0.060154914192541771, 1e-5); g.Mean() = "1 1"; x = "-1 4"; - BOOST_REQUIRE_CLOSE(g.Probability(x), 0.00144014867515135, 1e-5); - BOOST_REQUIRE_CLOSE(g.Probability(-x), 0.00133352162064845, 1e-5); + BOOST_REQUIRE_CLOSE(g.Probability(x), 0.0022506270186086271, 1e-5); + BOOST_REQUIRE_CLOSE(g.Probability(-x), 0.0016912950996661141, 1e-5); // Higher-dimensional case. x = "0 1 2 3 4"; g.Mean() = "5 6 3 3 2"; - covariance = "6 1 1 0 2;" - "0 7 1 0 1;" + covariance = "6 1 1 1 2;" + "1 7 1 0 0;" "1 1 4 1 1;" "1 0 1 7 0;" - "2 0 1 1 6"; + "2 0 1 0 6"; g.Covariance(std::move(covariance)); - BOOST_REQUIRE_CLOSE(g.Probability(x), 1.02531207499358e-6, 1e-5); - BOOST_REQUIRE_CLOSE(g.Probability(-x), 1.06784794079363e-8, 1e-5); + BOOST_REQUIRE_CLOSE(g.Probability(x), 1.4673143531128877e-06, 1e-5); + BOOST_REQUIRE_CLOSE(g.Probability(-x), 7.7404143494891786e-09, 1e-8); g.Mean() *= -1; - BOOST_REQUIRE_CLOSE(g.Probability(-x), 1.02531207499358e-6, 1e-5); - BOOST_REQUIRE_CLOSE(g.Probability(x), 1.06784794079363e-8, 1e-5); + BOOST_REQUIRE_CLOSE(g.Probability(-x), 1.4673143531128877e-06, 1e-5); + BOOST_REQUIRE_CLOSE(g.Probability(x), 7.7404143494891786e-09, 1e-8); } From 74f2656684e9e74bcf5463e334e76e04b4b24787 Mon Sep 17 00:00:00 2001 From: Stephen Tu Date: Tue, 20 Jan 2015 17:11:29 -0800 Subject: [PATCH 4/8] some more distribution test fixes --- src/mlpack/tests/distribution_test.cpp | 27 ++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/mlpack/tests/distribution_test.cpp b/src/mlpack/tests/distribution_test.cpp index 9041bb5f8b0..c5d9cedc6a6 100644 --- a/src/mlpack/tests/distribution_test.cpp +++ b/src/mlpack/tests/distribution_test.cpp @@ -190,9 +190,8 @@ BOOST_AUTO_TEST_CASE(GaussianUnivariateProbabilityTest) // A few more cases... arma::mat covariance; - covariance.set_size(g.Covariance().n_rows, g.Covariance().n_cols); - covariance.fill(2.0); + covariance = 2.0; g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("0.0")), 0.282094791773878, 1e-5); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("1.0")), 0.219695644733861, 1e-5); @@ -200,11 +199,11 @@ BOOST_AUTO_TEST_CASE(GaussianUnivariateProbabilityTest) 1e-5); g.Mean().fill(1.0); - covariance.fill(1.0); + covariance = 1.0; g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("1.0")), 0.398942280401433, 1e-5); - covariance.fill(2.0); + covariance = 2.0; g.Covariance(std::move(covariance)); BOOST_REQUIRE_CLOSE(g.Probability(arma::vec("-1.0")), 0.103776874355149, 1e-5); @@ -282,7 +281,11 @@ BOOST_AUTO_TEST_CASE(GaussianMultipointMultivariateProbabilityTest) { // Same case as before. arma::vec mean = "5 6 3 3 2"; - arma::mat cov = "6 1 1 0 2; 0 7 1 0 1; 1 1 4 1 1; 1 0 1 7 0; 2 0 1 1 6"; + arma::mat cov("6 1 1 1 2;" + "1 7 1 0 0;" + "1 1 4 1 1;" + "1 0 1 7 0;" + "2 0 1 0 6"); arma::mat points = "0 3 2 2 3 4;" "1 2 2 1 0 0;" @@ -292,16 +295,16 @@ BOOST_AUTO_TEST_CASE(GaussianMultipointMultivariateProbabilityTest) arma::vec phis; GaussianDistribution g(mean, cov); - g.Probability(points, phis); + g.LogProbability(points, phis); BOOST_REQUIRE_EQUAL(phis.n_elem, 6); - BOOST_REQUIRE_CLOSE(phis(0), 1.02531207499358e-6, 1e-5); - BOOST_REQUIRE_CLOSE(phis(1), 1.82353695848039e-7, 1e-5); - BOOST_REQUIRE_CLOSE(phis(2), 1.29759261892949e-6, 1e-5); - BOOST_REQUIRE_CLOSE(phis(3), 1.33218060268258e-6, 1e-5); - BOOST_REQUIRE_CLOSE(phis(4), 1.12120427975708e-6, 1e-5); - BOOST_REQUIRE_CLOSE(phis(5), 4.57951032485297e-7, 1e-5); + BOOST_REQUIRE_CLOSE(phis(0), -13.432076798791542, 1e-5); + BOOST_REQUIRE_CLOSE(phis(1), -15.814880322345738, 1e-5); + BOOST_REQUIRE_CLOSE(phis(2), -13.754462857772776, 1e-5); + BOOST_REQUIRE_CLOSE(phis(3), -13.283283233107898, 1e-5); + BOOST_REQUIRE_CLOSE(phis(4), -13.800326511545279, 1e-5); + BOOST_REQUIRE_CLOSE(phis(5), -14.900192463287908, 1e-5); } /** From 09e3013749107eca952774beb1881529a1df37ed Mon Sep 17 00:00:00 2001 From: Stephen Tu Date: Tue, 20 Jan 2015 17:16:49 -0800 Subject: [PATCH 5/8] ensure random covariance is positive definite --- src/mlpack/tests/gmm_test.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mlpack/tests/gmm_test.cpp b/src/mlpack/tests/gmm_test.cpp index d00a50eeec7..a83d5e9d512 100644 --- a/src/mlpack/tests/gmm_test.cpp +++ b/src/mlpack/tests/gmm_test.cpp @@ -493,6 +493,8 @@ BOOST_AUTO_TEST_CASE(GMMLoadSaveTest) arma::mat covariance = arma::randu( gmm.Component(i).Covariance().n_rows, gmm.Component(i).Covariance().n_cols); + covariance *= covariance.t(); + covariance += arma::eye(covariance.n_rows, covariance.n_cols); gmm.Component(i).Covariance(std::move(covariance)); } From cbdfab57a0ce0ab1549bbc11b94406f1f1f14b3e Mon Sep 17 00:00:00 2001 From: Stephen Tu Date: Tue, 20 Jan 2015 17:19:17 -0800 Subject: [PATCH 6/8] fix HMM load/save tests by ensuring covariance is positive definite --- src/mlpack/tests/hmm_test.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/mlpack/tests/hmm_test.cpp b/src/mlpack/tests/hmm_test.cpp index 9ae7016a46e..6ec90b295fc 100644 --- a/src/mlpack/tests/hmm_test.cpp +++ b/src/mlpack/tests/hmm_test.cpp @@ -995,6 +995,8 @@ BOOST_AUTO_TEST_CASE(GMMHMMLoadSaveTest) arma::mat covariance = arma::randu( hmm.Emission()[j].Component(i).Covariance().n_rows, hmm.Emission()[j].Component(i).Covariance().n_cols); + covariance *= covariance.t(); + covariance += arma::eye(covariance.n_rows, covariance.n_cols); hmm.Emission()[j].Component(i).Covariance(std::move(covariance)); } } @@ -1054,6 +1056,8 @@ BOOST_AUTO_TEST_CASE(GaussianHMMLoadSaveTest) arma::mat covariance = arma::randu( hmm.Emission()[j].Covariance().n_rows, hmm.Emission()[j].Covariance().n_cols); + covariance *= covariance.t(); + covariance += arma::eye(covariance.n_rows, covariance.n_cols); hmm.Emission()[j].Covariance(std::move(covariance)); } From 10ad21b0a71fbb5b5f1283d4bc8466a585ba31a7 Mon Sep 17 00:00:00 2001 From: Stephen Tu Date: Tue, 20 Jan 2015 17:23:17 -0800 Subject: [PATCH 7/8] remove stray commit --- .gitignore | 1 - 1 file changed, 1 deletion(-) diff --git a/.gitignore b/.gitignore index 726d15396af..1b2211df079 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1 @@ build* -*.swp From 046bdc44ab859d673ca7c013fb487989229f1e34 Mon Sep 17 00:00:00 2001 From: Stephen Tu Date: Thu, 22 Jan 2015 11:59:16 -0800 Subject: [PATCH 8/8] invert the lower triangular matrix the right way thanks to rcurtin for the note --- .../core/dists/gaussian_distribution.cpp | 24 +++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/mlpack/core/dists/gaussian_distribution.cpp b/src/mlpack/core/dists/gaussian_distribution.cpp index 969c33e45aa..cc84b8b86c0 100644 --- a/src/mlpack/core/dists/gaussian_distribution.cpp +++ b/src/mlpack/core/dists/gaussian_distribution.cpp @@ -33,9 +33,25 @@ void GaussianDistribution::Covariance(arma::mat&& covariance) void GaussianDistribution::FactorCovariance() { covLower = arma::chol(covariance, "lower"); - // tell arma that this is lower triangular matrix (for faster inversion) - covLower = arma::trimatl(covLower); - const arma::mat invCovLower = covLower.i(); + + // Comment from rcurtin: + // + // I think the use of the word "interpret" in the Armadillo documentation + // about trimatl and trimatu is somewhat misleading. What the function will + // actually do, when used in that context, is loop over the upper triangular + // part of the matrix and set it all to 0, so this ends up actually just + // burning cycles---also because the operator=() evaluates the expression and + // strips the knowledge that it's a lower triangular matrix. So then the call + // to .i() doesn't actually do anything smarter. + // + // But perusing fn_inv.hpp more closely, there is a specialization that will + // work when called like this: inv(trimatl(covLower)), and will use LAPACK's + // ?trtri functions. However, it will still set the upper triangular part to + // 0 after the method. That last part is unnecessary, but baked into + // Armadillo, so there's not really much that can be done about that without + // discussion with the Armadillo maintainer. + const arma::mat invCovLower = arma::inv(arma::trimatl(covLower)); + invCov = invCovLower.t() * invCovLower; double sign = 0.; arma::log_det(logDetCov, sign, covLower); @@ -146,7 +162,7 @@ void GaussianDistribution::Estimate(const arma::mat& observations, // Nothing in this Gaussian! At least set the covariance so that it's // invertible. covariance.diag() += 1e-50; - // TODO(stephentu): why do we allow this case? + FactorCovariance(); return; }