diff --git a/HISTORY.md b/HISTORY.md index fcb7568255..ef5bf1d95f 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -33,6 +33,8 @@ ### mlpack 3.3.0 ###### 2020-04-07 + * Added `Normal Distribution` to `ann/dists` (#2382). + * Templated return type of `Forward function` of loss functions (#2339). * Added `R2 Score` regression metric (#2323). diff --git a/src/mlpack/methods/ann/dists/CMakeLists.txt b/src/mlpack/methods/ann/dists/CMakeLists.txt index eea77c413c..fadec3c3c1 100644 --- a/src/mlpack/methods/ann/dists/CMakeLists.txt +++ b/src/mlpack/methods/ann/dists/CMakeLists.txt @@ -3,6 +3,8 @@ set(SOURCES bernoulli_distribution.hpp bernoulli_distribution_impl.hpp + normal_distribution.hpp + normal_distribution_impl.hpp ) # Add directory name to sources. diff --git a/src/mlpack/methods/ann/dists/normal_distribution.hpp b/src/mlpack/methods/ann/dists/normal_distribution.hpp new file mode 100644 index 0000000000..8f9844fd49 --- /dev/null +++ b/src/mlpack/methods/ann/dists/normal_distribution.hpp @@ -0,0 +1,145 @@ +/** + * @file normal_distribution.hpp + * @author xiaohong ji + * @author Nishant Kumar + * + * Definition of the Normal distribution class. + * + * 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_METHODS_ANN_DISTRIBUTIONS_NORMAL_DISTRIBUTION_HPP +#define MLPACK_METHODS_ANN_DISTRIBUTIONS_NORMAL_DISTRIBUTION_HPP + +#include +#include "../activation_functions/logistic_function.hpp" + +namespace mlpack { +namespace ann /** Artificial Neural Network. */ { + +/** + * Implementation of the Normal Distribution function. + * + * Normal distribution is a function which accepts a mean and a standard deviation + * term and creates a probablity distribution out of it. + * + * @tparam DataType Type of the input data. (arma::colvec, arma::mat, + * arma::sp_mat or arma::cube). + */ +template +class NormalDistribution +{ + public: + /** + * Default constructor, which creates a Normal distribution with zero + * dimension. + */ + NormalDistribution(); + + /** + * Create a Normal distribution with the given mean and sigma. + * + * @param mean The mean of the normal distribution. + * @param sigma The standard deviation of the normal distribution. + */ + NormalDistribution(const DataType& mean, const DataType& sigma); + + /** + * Return the probabilities of the given matrix of observations. + * + * @param observation The observation matrix. + */ + DataType Probability(const DataType& observation) const + { + return arma::exp(LogProbability(observation)); + } + + /** + * Return the log probabilities of the given matrix of observations. + * + * @param observation The observation matrix. + */ + DataType LogProbability(const DataType& observation) const; + + /** + * Stores the gradient of the probabilities of the observations + * with respect to mean and standard deviation. + * + * @param observation The observation matrix. + * @param dmu The gradient with respect to mean. + * @param dsigma The gradient with respect to the standard deviation. + */ + void ProbBackward(const DataType& observation, + DataType& dmu, + DataType& dsigma) const; + + /** + * Calculates the normal probability density function for each + * data point (column) in the given matrix. + * + * @param x The observation matrix. + * @param probabilities Output probabilities for each input observation. + */ + void Probability(const DataType& x, DataType& probabilities) const + { + probabilities = Probability(x); + } + + /** + * Calculates the log of normal probability density function for each + * data point (column) in the given matrix. + * + * @param x The observation matrix. + * @param probabilities Output log probabilities for each input observation. + */ + void LogProbability(const DataType& x, DataType& probabilities) const + { + probabilities = LogProbability(x); + } + + /** + * Return a randomly generated observation according to the probability + * distribution defined by this object. + * + * @return Random observation from this Normal distribution. + */ + DataType Sample() const; + + //! Get the mean. + const DataType& Mean() const { return mean; } + + //! Modify the mean. + DataType& Mean() { return mean; } + + //! Get the standard deviation. + const DataType& StandardDeviation() const { return sigma; } + + //! Modify the standard deviation. + DataType& StandardDeviation() { return sigma; } + + //! Return the dimensionality of this distribution. + size_t Dimensionality() const { return mean.n_elem; } + + /** + * Serialize the distribution. + */ + template + void serialize(Archive& ar, const unsigned int /* version */); + + private: + //! Mean of the distribution. + DataType mean; + + //! Standard deviation of the distribution. + DataType sigma; +}; // class NormalDistribution + +} // namespace ann +} // namespace mlpack + +// Include implementation. +#include "normal_distribution_impl.hpp" + +#endif diff --git a/src/mlpack/methods/ann/dists/normal_distribution_impl.hpp b/src/mlpack/methods/ann/dists/normal_distribution_impl.hpp new file mode 100644 index 0000000000..4212f30b70 --- /dev/null +++ b/src/mlpack/methods/ann/dists/normal_distribution_impl.hpp @@ -0,0 +1,78 @@ +/** + * @file normal_distribution_impl.hpp + * @author xiaohong ji + * @author Nishant Kumar + * + * Implementation of the Normal Distribution class. + * + * 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_METHODS_ANN_DISTRIBUTIONS_NORMAL_DISTRIBUTION_IMPL_HPP +#define MLPACK_METHODS_ANN_DISTRIBUTIONS_NORMAL_DISTRIBUTION_IMPL_HPP + +// In case it hasn't yet been included. +#include "normal_distribution.hpp" + +namespace mlpack { +namespace ann /** Artificial Neural Network. */ { + +template +NormalDistribution::NormalDistribution() +{ + // Nothing to do here. +} + +template +NormalDistribution::NormalDistribution( + const DataType& mean, + const DataType& sigma) : + mean(mean), + sigma(sigma) +{ + // Nothing to do here. +} + +template +DataType NormalDistribution::Sample() const +{ + return sigma * arma::randn(mean.n_elem) + mean; +} + +template +DataType NormalDistribution::LogProbability( + const DataType& observation) const +{ + const DataType v1 = arma::log(sigma) + std::log(std::sqrt(2 * M_PI)); + const DataType v2 = arma::square(observation - mean) / + (2 * arma::square(sigma)); + return (-v1 - v2); +} + +template +void NormalDistribution::ProbBackward( + const DataType& observation, + DataType& dmu, + DataType& dsigma) const +{ + dmu = (observation - mean) / (arma::square(sigma)) % Probability(observation); + dsigma = (- 1.0 / sigma + + (arma::square(observation - mean) / arma::pow(sigma, 3))) + % Probability(observation); +} + +template +template +void NormalDistribution::serialize(Archive& ar, + const unsigned int /* version */) +{ + ar & BOOST_SERIALIZATION_NVP(mean); + ar & BOOST_SERIALIZATION_NVP(sigma); +} + +} // namespace ann +} // namespace mlpack + +#endif diff --git a/src/mlpack/tests/ann_dist_test.cpp b/src/mlpack/tests/ann_dist_test.cpp index 6c6dad1371..face8565bd 100644 --- a/src/mlpack/tests/ann_dist_test.cpp +++ b/src/mlpack/tests/ann_dist_test.cpp @@ -1,6 +1,7 @@ /** * @file ann_dist_test.cpp * @author Atharva Khandait + * @author Nishant Kumar * * Tests the ann distributions. * @@ -12,6 +13,7 @@ #include #include +#include #include #include @@ -127,4 +129,174 @@ BOOST_AUTO_TEST_CASE(JacobianBernoulliDistributionLogisticTest) } } +/** + * Normal Distribution module test. + */ +BOOST_AUTO_TEST_CASE(NormalDistributionTest) +{ + arma::vec mu = {1.1, 1.2, 1.5, 1.7}; + arma::vec sigma = {0.1, 0.11, 0.5, 0.23}; + + ann::NormalDistribution<> normalDist(mu, sigma); + + arma::vec x = {1.05, 1.1, 1.7, 2.5}; + + arma::vec prob; + normalDist.LogProbability(x, prob); + + // Testing output of log probability for some random mu, sigma and x. + BOOST_REQUIRE_CLOSE(prob[0], 1.2586464, 1e-3); + BOOST_REQUIRE_CLOSE(prob[1], 0.8751131, 1e-3); + BOOST_REQUIRE_CLOSE(prob[2], -0.30579138, 1e-3); + BOOST_REQUIRE_CLOSE(prob[3], -5.498411, 1e-3); + + arma::vec dmu, dsigma; + normalDist.ProbBackward(x, dmu, dsigma); + + // Testing output of dmu and dsigma for some random mu, sigma and x. + BOOST_REQUIRE_CLOSE(dmu[0], -17.603287, 1e-3); + BOOST_REQUIRE_CLOSE(dsigma[0], -26.40487, 1e-3); + BOOST_REQUIRE_CLOSE(dmu[1], -19.827663, 1e-3); + BOOST_REQUIRE_CLOSE(dsigma[1], -3.7852707, 1e-3); + BOOST_REQUIRE_CLOSE(dmu[2], 0.5892323, 1e-3); + BOOST_REQUIRE_CLOSE(dsigma[2], -1.2373875, 1e-3); + BOOST_REQUIRE_CLOSE(dmu[3], 0.061901994, 1e-3); + BOOST_REQUIRE_CLOSE(dsigma[3], 0.19751444, 1e-3); +} + +/** + * Jacobian Normal Distribution module test for mean. + */ +BOOST_AUTO_TEST_CASE(JacobianNormalDistributionMeanTest) +{ + for (size_t i = 0; i < 5; i++) + { + const size_t targetElements = math::RandInt(2, 1000); + + arma::mat mu; + mu.randn(targetElements, 1); + + arma::mat sigma; + sigma.randu(targetElements, 1); + + arma::mat x; + x.randn(targetElements, 1); + + NormalDistribution<> module(mu, sigma); + + const double perturbation = 1e-6; + arma::mat output, outputA, outputB, jacobianA, jacobianB; + + // Initialize the jacobian matrix. + module.Probability(x, output); + jacobianA = arma::zeros(x.n_elem, output.n_elem); + + for (size_t j = 0; j < x.n_elem; ++j) + { + double original = module.Mean()(j); + module.Mean()(j) = original - perturbation; + module.Probability(x, outputA); + module.Mean()(j) = original + perturbation; + module.Probability(x, outputB); + module.Mean()(j) = original; + + outputB -= outputA; + outputB /= 2 * perturbation; + jacobianA.row(j) = outputB.t(); + } + + // Initialize the derivative parameter. + arma::mat deriv = arma::zeros(output.n_rows, output.n_cols); + + // Share the derivative parameter. + arma::mat derivTemp = arma::mat(deriv.memptr(), deriv.n_rows, deriv.n_cols, + false, false); + + // Initialize the jacobian matrix. + jacobianB = arma::zeros(mu.n_elem, output.n_elem); + + for (size_t k = 0; k < derivTemp.n_elem; ++k) + { + deriv.zeros(); + derivTemp(k) = 1; + + arma::mat deltaMu, deltaSigma; + module.ProbBackward(x, deltaMu, deltaSigma); + + jacobianB.col(k) = deltaMu % deriv; + } + + BOOST_REQUIRE_LE(arma::max(arma::max(arma::abs(jacobianA - jacobianB))), + 5e-3); + } +} + +/** + * Jacobian Normal Distribution module test for standard deviation. + */ +BOOST_AUTO_TEST_CASE(JacobianNormalDistributionStandardDeviationTest) +{ + for (size_t i = 0; i < 5; i++) + { + const size_t targetElements = math::RandInt(2, 1000); + + arma::mat mu; + mu.randn(targetElements, 1); + + arma::mat sigma; + sigma.randu(targetElements, 1); + + arma::mat x; + x.randn(targetElements, 1); + + NormalDistribution<> module(mu, sigma); + + const double perturbation = 1e-6; + arma::mat output, outputA, outputB, jacobianA, jacobianB; + + // Initialize the jacobian matrix. + module.Probability(x, output); + jacobianA = arma::zeros(x.n_elem, output.n_elem); + + for (size_t j = 0; j < x.n_elem; ++j) + { + double original = module.StandardDeviation()(j); + module.StandardDeviation()(j) = original - perturbation; + module.Probability(x, outputA); + module.StandardDeviation()(j) = original + perturbation; + module.Probability(x, outputB); + module.StandardDeviation()(j) = original; + + outputB -= outputA; + outputB /= 2 * perturbation; + jacobianA.row(j) = outputB.t(); + } + + // Initialize the derivative parameter. + arma::mat deriv = arma::zeros(output.n_rows, output.n_cols); + + // Share the derivative parameter. + arma::mat derivTemp = arma::mat(deriv.memptr(), deriv.n_rows, deriv.n_cols, + false, false); + + // Initialize the jacobian matrix. + jacobianB = arma::zeros(sigma.n_elem, output.n_elem); + + for (size_t k = 0; k < derivTemp.n_elem; ++k) + { + deriv.zeros(); + derivTemp(k) = 1; + + arma::mat deltaMu, deltaSigma; + module.ProbBackward(x, deltaMu, deltaSigma); + + jacobianB.col(k) = deltaSigma % deriv; + } + + BOOST_REQUIRE_LE(arma::max(arma::max(arma::abs(jacobianA - jacobianB))), + 5e-3); + } +} + + BOOST_AUTO_TEST_SUITE_END();