From a9985e3e1ca593a8c2d0baa88dea648d884440de Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Fri, 7 Dec 2018 09:02:25 -0700 Subject: [PATCH] Hoyer sparsity [CI SKIP] --- doc/vector_functionals/vector_functionals.qbk | 55 ++++++++++++++++--- .../boost/math/tools/vector_functionals.hpp | 42 ++++++++++++-- test/vector_functionals_test.cpp | 37 +++++++++++++ 3 files changed, 121 insertions(+), 13 deletions(-) diff --git a/doc/vector_functionals/vector_functionals.qbk b/doc/vector_functionals/vector_functionals.qbk index 046a1bcd3..03a5724e2 100644 --- a/doc/vector_functionals/vector_functionals.qbk +++ b/doc/vector_functionals/vector_functionals.qbk @@ -116,6 +116,8 @@ As a reminder, remember that to actually /get/ vectorization, compile with `-mar We now describe each functional in detail. +[heading Mean] + Compute the mean of a container: std::vector v{1,2,3,4,5}; @@ -125,6 +127,8 @@ The implementation follows [@https://doi.org/10.1137/1.9780898718027 Higham 1.6a The only requirement on the input is that it must be forward iterable, so you can use Eigen vectors, ublas vectors, Armadillo vectors, or a `std::forward_list` to hold your data. +[heading Mean and Sample variance] + Compute the mean and sample variance: std::vector v{1,2,3,4,5}; @@ -132,10 +136,12 @@ Compute the mean and sample variance: The implementation follows [@https://doi.org/10.1137/1.9780898718027 Higham 1.6b]. Note that we do not provide computation of sample variance alone; -we are unaware of any one-pass, numerically stable computation of sample variance which does not require simultaneous computation of the mean. +we are unaware of any one-pass, numerically stable computation of sample variance which does not simultaneously generate the mean. If the mean is not required, simply ignore it. The input datatype must be forward iterable and the range `[first, last)` must contain at least two elements. +[heading Median] + Compute the median of a dataset: std::vector v{1,2,3,4,5}; @@ -145,7 +151,7 @@ Compute the median of a dataset: The calculation of the median is a thin wrapper around the C++11 [@https://en.cppreference.com/w/cpp/algorithm/nth_element nth-element]. Therefore, all requirements of `nth_element` are inherited by the median calculation. - +[heading Sup norm] Compute the sup norm of a dataset: std::vector v{-3, 2, 1}; @@ -158,6 +164,8 @@ Compute the sup norm of a dataset: Note how the calculation of \u2113[super p] norms can be performed in both real and complex arithmetic. +[heading Gini Coefficient] + Compute the Gini coefficient of a dataset: std::vector v{1,0,0,0}; @@ -170,8 +178,13 @@ Compute the Gini coefficient of a dataset: /Nota bene: The input data is altered-in particular, it is sorted./ /Nota bene:/ Different authors use different conventions regarding the overall scale of the Gini coefficient. -We have chosen to follow [@https://arxiv.org/pdf/0811.4706.pdf Hurley and Rickard's definition], which [@https://en.wikipedia.org/wiki/Gini_coefficient Wikipedia] calls a "consistent estimator" of the population Gini coefficient. +We have chosen to follow [@https://arxiv.org/pdf/0811.4706.pdf Hurley and Rickard's definition], which [@https://en.wikipedia.org/wiki/Gini_coefficient Wikipedia] calls a "sample Gini coefficient". Hurley and Rickard's definition places the Gini coefficient in the range [0,1]; Wikipedia's population Gini coefficient is in the range [0, 1 - 1/N]. +If you wish to convert the Boost Gini coefficient to the population Gini coefficient, multiply by (/n/-1)/ /n/. + +/Nota bene:/ There is essentially no reason to pass negative values to the Gini coefficient function. +However, since a single use case (measuring wealth inequality when some people have negative wealth) exists, we do not throw an exception when negative values are encountered. +You should have /very/ good cause to pass negative values to the Gini coefficient calculator. The Gini coefficient, first used to measure wealth inequality, is also one of the best measures of the sparsity of an expansion in a basis. A sparse expansion has most of its norm concentrated in just a few coefficients, making the connection with wealth inequality obvious. @@ -179,18 +192,42 @@ However, for measuring sparsity, the phase of the numbers is irrelevant, so `abs std::vector> v{{0,1}, {0,0}, {0,0}, {0,0}}; double abs_gini = absolute_gini_coefficient(v.begin(), v.end()); - // abs_gini = 1 + // now abs_gini = 1 + std::vector> w{{0,1}, {1,0}, {0,-1}, {-1,0}}; double abs_gini = absolute_gini_coefficient(w.begin(), w.end()); - // abs_gini = 0 + // now abs_gini = 0 + + std::vector u{-1, 1, -1}; + double abs_gini = absolute_gini_coefficient(u.begin(), u.end()); + // now abs_gini = 0 + +Again, Wikipedia denotes our scaling as a "sample Gini coefficient". +We chose this scaling because it always returns unity for a vector which has only one nonzero coefficient. + +If sorting the input data is to much expense for a sparsity measure (is it going to be perfect anyway?), +consider using `hoyer_sparsity`. + +[heading Hoyer Sparsity] + +The Hoyer sparsity measure uses a normalized ratio of the \u2113[super 1] and \u2113[super 2] norms. +As the name suggests, it is used to measure sparsity in an expansion in some basis. + +The Hoyer sparsity computes ([radic]/N/ - \u2113[super 1](v)/\u2113[super 2](v))/([radic]N -1). + +Usage: + + std::vector v{1,0,0}; + Real hs = boost::math::tools::hoyer_sparsity(v.begin(), v.end()); + // hs = 1 + + +[heading \u2113[super /p/] norm] -[heading Examples] +[heading \u2113[super 0] norm] -[heading Performance] - -[heading Caveats] [heading References] diff --git a/include/boost/math/tools/vector_functionals.hpp b/include/boost/math/tools/vector_functionals.hpp index 8e77e97f9..2cc052680 100644 --- a/include/boost/math/tools/vector_functionals.hpp +++ b/include/boost/math/tools/vector_functionals.hpp @@ -211,6 +211,13 @@ auto gini_coefficient(ForwardIterator first, ForwardIterator last) denom += *it; ++i; } + + // If the l1 norm is zero, all elements are zero, so every element is the same. + if (denom == 0) + { + return Real(0); + } + return ((2*num)/denom - i)/(i-2); } @@ -227,19 +234,46 @@ auto absolute_gini_coefficient(ForwardIterator first, ForwardIterator last) decltype(abs(*first)) i = 1; decltype(abs(*first)) num = 0; decltype(abs(*first)) denom = 0; - std::cout << "{"; - for (auto it = first; it != last; ++it) { - std::cout << abs(*it) << ", "; + for (auto it = first; it != last; ++it) + { decltype(abs(*first)) tmp = abs(*it); num += tmp*i; denom += tmp; ++i; } - std::cout << "}\n"; + + // If the l1 norm is zero, all elements are zero, so every element is the same. + if (denom == 0) + { + decltype(abs(*first)) zero = 0; + return zero; + } return ((2*num)/denom - i)/(i-2); } +// The Hoyer sparsity measure is defined in: +// https://arxiv.org/pdf/0811.4706.pdf +template +auto hoyer_sparsity(ForwardIterator first, ForwardIterator last) +{ + using std::abs; + using std::sqrt; + typedef typename std::remove_const())>::type>::type RealOrComplex; + BOOST_ASSERT_MSG(first != last, "Computation of the Hoyer sparsity requires at least one sample."); + decltype(abs(*first)) l1 = 0; + decltype(abs(*first)) l2 = 0; + decltype(abs(*first)) n = 0; + for (auto it = first; it != last; ++it) + { + decltype(abs(*first)) tmp = abs(*it); + l1 += tmp; + l2 += tmp*tmp; + n += 1; + } + decltype(abs(*first)) rootn = sqrt(n); + return (rootn - l1/sqrt(l2) )/ (rootn - 1); +} }}} #endif diff --git a/test/vector_functionals_test.cpp b/test/vector_functionals_test.cpp index c82037932..c535ba702 100644 --- a/test/vector_functionals_test.cpp +++ b/test/vector_functionals_test.cpp @@ -141,6 +141,33 @@ void test_gini_coefficient() v[2] = 1; gini = boost::math::tools::gini_coefficient(v.begin(), v.end()); BOOST_TEST(abs(gini) < tol); + + v[0] = 0; + v[1] = 0; + v[2] = 0; + gini = boost::math::tools::gini_coefficient(v.begin(), v.end()); + BOOST_TEST(abs(gini) < tol); +} + +template +void test_hoyer_sparsity() +{ + using std::sqrt; + Real tol = 5*std::numeric_limits::epsilon(); + std::vector v{1,0,0}; + Real hs = boost::math::tools::hoyer_sparsity(v.begin(), v.end()); + BOOST_TEST(abs(hs - 1) < tol); + + // Does it work with constant iterators? + hs = boost::math::tools::hoyer_sparsity(v.cbegin(), v.cend()); + BOOST_TEST(abs(hs - 1) < tol); + + v[0] = 1; + v[1] = 1; + v[2] = 1; + hs = boost::math::tools::hoyer_sparsity(v.cbegin(), v.cend()); + BOOST_TEST(abs(hs) < tol); + } template @@ -166,6 +193,7 @@ void test_absolute_gini_coefficient() gini = boost::math::tools::absolute_gini_coefficient(w.begin(), w.end()); BOOST_TEST(abs(gini) < tol); + // The Gini index is invariant under "cloning": If w = v \oplus v, then G(w) = G(v). } template @@ -174,6 +202,11 @@ void test_l0_norm() std::vector v{0,0,1}; size_t count = boost::math::tools::l0_norm(v.begin(), v.end()); BOOST_TEST_EQ(count, 1); + + // Compiles with cbegin()/cend()? + count = boost::math::tools::l0_norm(v.cbegin(), v.cend()); + BOOST_TEST_EQ(count, 1); + } int main() @@ -209,6 +242,10 @@ int main() test_absolute_gini_coefficient(); test_absolute_gini_coefficient(); + test_hoyer_sparsity(); + test_hoyer_sparsity(); + test_hoyer_sparsity(); + test_l0_norm(); test_l0_norm(); test_l0_norm();