# A collection of vector functionals

I would like some feedback for a collection of what I call "vector functionals", by which I mean maps $$\\ell\colon \mathbb{K}^{n} \to \mathbb{K}\$$, where the field $$\\mathbb{K} = \mathbb{R}\$$ or $$\\mathbb{C}\$$.

Many of these maps have trivial implementations, but still they are well-defined and must be put into a function somewhere, and they still take time and effort and hence it makes sense to place them into some sort of canonical point of authority, which I am trying to do. In addition, a few of the maps (like the Hoyer sparsity and the Gini coefficient) could use popularizing.

I would like feedback on the following:

• Is the design clean? Can it be improved? I have a particular antipathy towards the typedef typename std::remove_const<... everywhere, and I'd like some way to make the iterators for the non-modifying functions "const" in some sense.
• What other "vector functionals" do you think should be added? Or is the notion of "vector functionals" simply too broad, and this should be split into those functions traditionally found in statistics, split into numerical analysis, split into signal processing?
• Have I missed any edge cases? For instance, one common edge case is overflow in $$\\ell^{p}\$$ norms. Another is not using Welford's algorithm for computing variance. In particular, I could not find any literature on rapid and numerically stable evaluation of the Gini coefficient.
• I checked for generation of avx instructions on godbolt for a few of the functions, but is there anything I'm doing which will prevent generation of the AVX instructions in general?

Without further ado, the PR is here, the documentation here, and the complete code below:

#ifndef BOOST_MATH_TOOLS_VECTOR_FUNCTIONALS_HPP
#define BOOST_MATH_TOOLS_VECTOR_FUNCTIONALS_HPP

#include <algorithm>
#include <iterator>
#include <boost/type_traits.hpp>
#include <boost/assert.hpp>
#include <boost/multiprecision/detail/number_base.hpp>

/*
* A set of tools for computing scalar quantities associated with lists of numbers.
*/

namespace boost{ namespace math{ namespace tools {

template<class ForwardIterator>
auto
mean(ForwardIterator first, ForwardIterator last)
{
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type Real;
BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");
Real mu = 0;
Real i = 1;
for(auto it = first; it != last; ++it) {
mu = mu + (*it - mu)/i;
i += 1;
}
return mu;
}

template<class ForwardIterator>
auto
mean_and_population_variance(ForwardIterator first, ForwardIterator last)
{
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type Real;
BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute mean and variance.");
// Higham, Accuracy and Stability, equation 1.6a and 1.6b:
Real M = *first;
Real Q = 0;
Real k = 2;
for (auto it = first + 1; it != last; ++it)
{
Real tmp = *it - M;
Q = Q + ((k-1)*tmp*tmp)/k;
M = M + tmp/k;
k += 1;
}

return std::make_pair(M, Q/(k-1));
}

template<class RandomAccessIterator>
auto median(RandomAccessIterator first, RandomAccessIterator last)
{
size_t num_elems = std::distance(first, last);
BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined.");
if (num_elems & 1)
{
auto middle = first + (num_elems - 1)/2;
nth_element(first, middle, last);
return *middle;
}
else
{
auto middle = first + num_elems/2 - 1;
nth_element(first, middle, last);
nth_element(middle, middle+1, last);
return (*middle + *(middle+1))/2;
}
}

template<class RandomAccessIterator>
auto absolute_median(RandomAccessIterator first, RandomAccessIterator last)
{
using std::abs;
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<RandomAccessIterator>())>::type>::type RealOrComplex;
size_t num_elems = std::distance(first, last);
BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined.");
auto comparator = [](RealOrComplex a, RealOrComplex b) { return abs(a) < abs(b);};
if (num_elems & 1)
{
auto middle = first + (num_elems - 1)/2;
nth_element(first, middle, last, comparator);
return abs(*middle);
}
else
{
auto middle = first + num_elems/2 - 1;
nth_element(first, middle, last, comparator);
nth_element(middle, middle+1, last, comparator);
return (abs(*middle) + abs(*(middle+1)))/abs(static_cast<RealOrComplex>(2));
}
}

// Mallat, "A Wavelet Tour of Signal Processing", equation 2.60:
template<class ForwardIterator>
auto total_variation(ForwardIterator first, ForwardIterator last)
{
using std::abs;
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type Real;
BOOST_ASSERT_MSG(first != last && std::next(first) != last, "At least two samples are required to compute the total variation.");
Real tv = 0;
auto it = first;
Real tmp = *it;
while (++it != last)
{
tv += abs(*it - tmp);
tmp = *it;
}
return tv;
}

// Mallat, equation 10.4 uses the base-2 logarithm.
template<class ForwardIterator>
auto shannon_entropy(ForwardIterator first, ForwardIterator last)
{
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type Real;
using std::log;
Real entropy = 0;
for (auto it = first; it != last; ++it)
{
if (*it != 0)
{
entropy += (*it)*log(*it);
}
}
return -entropy;
}

template<class ForwardIterator>
auto sup_norm(ForwardIterator first, ForwardIterator last)
{
BOOST_ASSERT_MSG(first != last, "At least one value is required to compute the sup norm.");
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type RealOrComplex;
using std::abs;
if constexpr (boost::is_complex<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_complex)
{
auto it = std::max_element(first, last, [](RealOrComplex a, RealOrComplex b) { return abs(b) > abs(a); });
return abs(*it);
}
else
{
auto pair = std::minmax_element(first, last);
if (abs(*pair.first) > abs(*pair.second))
{
return abs(*pair.first);
}
else
{
return abs(*pair.second);
}
}
}

template<class ForwardIterator>
auto l1_norm(ForwardIterator first, ForwardIterator last)
{
using std::abs;
decltype(abs(*first)) l1 = 0;
for (auto it = first; it != last; ++it)
{
l1 += abs(*first);
}
return l1;
}

template<class ForwardIterator>
auto l2_norm(ForwardIterator first, ForwardIterator last)
{
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type RealOrComplex;
using std::abs;
using std::norm;
using std::sqrt;
using std::is_floating_point;
if constexpr (boost::is_complex<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_complex)
{
typedef typename RealOrComplex::value_type Real;
Real l2 = 0;
for (auto it = first; it != last; ++it)
{
l2 += norm(*it);
}
Real result = sqrt(l2);
if (!isfinite(result))
{
Real a = sup_norm(first, last);
l2 = 0;
for (auto it = first; it != last; ++it)
{
l2 += norm(*it/a);
}
return a*sqrt(l2);
}
return result;
}
else if constexpr (is_floating_point<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_floating_point)
{
RealOrComplex l2 = 0;
for (auto it = first; it != last; ++it)
{
l2 += (*it)*(*it);
}
RealOrComplex result = sqrt(l2);
if (!isfinite(result))
{
RealOrComplex a = sup_norm(first, last);
l2 = 0;
for (auto it = first; it != last; ++it)
{
RealOrComplex tmp = *it/a;
l2 += tmp*tmp;
}
return a*sqrt(l2);
}
return result;
}
}

template<class ForwardIterator>
size_t l0_pseudo_norm(ForwardIterator first, ForwardIterator last)
{
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type RealOrComplex;
size_t count = 0;
for (auto it = first; it != last; ++it)
{
if (*it != RealOrComplex(0))
{
++count;
}
}
return count;
}

template<class ForwardIterator>
auto lp_norm(ForwardIterator first, ForwardIterator last, typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type p)
{
using std::pow;
using std::is_floating_point;
using std::isfinite;
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type RealOrComplex;
if constexpr (boost::is_complex<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_complex)
{
BOOST_ASSERT_MSG(p.real() >= 0, "For p < 0, the lp norm is not a norm.");
BOOST_ASSERT_MSG(p.imag() == 0, "For imaginary p, the lp norm is not a norm.");
using std::norm;
decltype(p.real()) lp = 0;
for (auto it = first; it != last; ++it)
{
lp += pow(norm(*it), p.real()/2);
}

auto result = pow(lp, 1/p.real());
if (!isfinite(result))
{
auto a = boost::math::tools::sup_norm(first, last);
decltype(p.real()) lp = 0;
for (auto it = first; it != last; ++it)
{
lp += pow(abs(*it)/a, p.real());
}
result = a*pow(lp, 1/p.real());
}
return result;
}
else if constexpr (is_floating_point<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_floating_point)
{
BOOST_ASSERT_MSG(p >= 0, "For p < 0, the lp norm is not a norm");
RealOrComplex lp = 0;

for (auto it = first; it != last; ++it)
{
lp += pow(abs(*it), p);
}

RealOrComplex result = pow(lp, 1/p);
if (!isfinite(result))
{
RealOrComplex a = boost::math::tools::sup_norm(first, last);
lp = 0;
for (auto it = first; it != last; ++it)
{
lp += pow(abs(*it)/a, p);
}
result = a*pow(lp, 1/p);
}
return result;
}
else
{
BOOST_ASSERT_MSG(false, "Unable to determine if the input type is real or complex.");
}
}

template<class ForwardIterator>
auto gini_coefficient(ForwardIterator first, ForwardIterator last)
{
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type Real;
BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");

std::sort(first, last);

Real i = 1;
Real num = 0;
Real denom = 0;
for (auto it = first; it != last; ++it) {
num += *it*i;
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);
}

template<class ForwardIterator>
auto absolute_gini_coefficient(ForwardIterator first, ForwardIterator last)
{
typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type RealOrComplex;
BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");

std::sort(first, last,  [](RealOrComplex a, RealOrComplex b) { return abs(b) > abs(a); });

decltype(abs(*first)) i = 1;
decltype(abs(*first)) num = 0;
decltype(abs(*first)) denom = 0;
for (auto it = first; it != last; ++it)
{
decltype(abs(*first)) tmp = abs(*it);
num += tmp*i;
denom += tmp;
++i;
}

// 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<class ForwardIterator>
auto hoyer_sparsity(const ForwardIterator first, const ForwardIterator last)
{
using std::abs;
using std::sqrt;
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


What is all this for?

typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type Real;


That's a really complex way of doing something.

Let's break this down.

using T1 = std::declval<ForwardIterator>  // Returns the reference version of ForwardIterator

using T2 = decltype(*T1())                // Creates a temporary object.
// De-reference the temporary
// Gets the type of the de-referenced value.2

using T3 = std::remove_reference<T2>::type; // removes references.
using T4 = std::remove_const<T3>::type;     // removes an const ness.


Simpler to write:

typedef typename std::iterator_traits<ForwardIterator>::value_type Real;


Or if you want to use the modern syntax.

using Real = typename std::iterator_traits<ForwardIterator>::value_type;


Or if you are like me and hate using the typename everywhere you can create your own template using to make this simpler.

template<typename I>
using IterValue = typename std::iterator_traits<I>::value_type;


Now in your code you can use:

template<class ForwardIterator, ReturnType = IterValue<ForwardIterator>>
ReturnType mean(ForwardIterator first, ForwardIterator last)
{
BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");
int i = 1;            // The i value is always an int no matter what the data type is

ReturnType mu = 0;
for(auto it = first; it != last; ++it) {
mu = mu + (*it - mu)/i;             // Sure!
i += 1;
}
return mu;
}

• The using Real = typename ForwardIterator::value_type; is such a clean syntax, but I can't get it to compile with a std::array input type on g++-8 and clang: fatal error: type 'float *' cannot be used prior to '::' because it has no members – user14717 Dec 8 '18 at 20:47
• I kept i as a Real because I was worried it would generate an int -> Real cast on every iteration before the division. However, adding your change didn't generate any meaningful performance differences either way, so your simpler option is probably the best. – user14717 Dec 8 '18 at 20:59
• @user14717: Fixed the issue with iterators that are pointers. – Martin York Dec 8 '18 at 21:24
• You could use auto for most type deductions (instead of manually deducing Real). Only caveat: The algorithms might not support iterators returning proxy objects (looking at you, std::vector<bool>::iterator). Technically, that's still standard library compatible, though, since that kind of iterator can at most fulfill the InputIterator requirements, and the algorithms seem to require ForwardIterators, for which reference must be equal to value_type& or const value_type&. – hoffmale Dec 9 '18 at 0:15
• @hoffmale Yes. Auto would be even better. Tried updating the answer. But the problem here is the auto mu = 0; Does not do what we need it to do. I could do auto mu = *it; mu = 0; But that seems a bit complex and likely to be messed with by a maintainer unless I have a significant comment. So we are back to manually working out the type. – Martin York Dec 9 '18 at 0:56

I think it would be interesting to have such a collection of tools be standardized.

However, it will require a bit more thinking. I would not want to use these tools as they are now.

# mean()

I wrote this simple test program:

#include <vector_functionals.hpp>
#include <iostream>
#include <vector>

int main() {
std::vector<std::uint8_t> test1{ 5, 4, 4 };
std::vector<std::uint8_t> test2{ 4, 4, 5 };
std::cout << "mean of test1 = " << (int)boost::math::tools::mean( test1.begin(), test1.end() ) << '\n';
std::cout << "mean of test2 = " << (int)boost::math::tools::mean( test2.begin(), test2.end() ) << '\n';
}


This is the output:

mean of test1 = 5
mean of test2 = 4


The order of the elements in the array affects the mean? No!

Similarly, it computes the mean of integers {255, 0, 0, 0} not as 255.0/4=63.75, not as that value rounded up or down, but as 65.

mean_and_population_variance() will have similar issues, and it will overflow.

I would expect mean() to return a double if the input is integer. But determining the correct output type is not possible to do in a generic way. I'd suggest you require type traits to exist for the arithmetic type being used, and that one of these type traits be the type of the output for these functions. I believe other Boost modules work similarly.

# median()

This one is really only useful for arithmetic types, and that is a shame. It'd be interesting if it could return the median of a set of words, for example. But you can't do (word1 + word2) / 2.

• You're right; I think it should give a compilation error for integer types, or it should do "arg promotion" and return a double when integral types are provided. Why do you think mean_and_population_variance will overflow? – user14717 Dec 8 '18 at 16:37
• @user14717: again in the uint8 case: if the first element is 255, the second is 0, then tmp is 255 and tmp*tmp will overflow. If they are sorted the other way around, then tmp will be 0 instead of -255. Unsigned arithmetic cannot be used in this function, and variance computation just needs floating-point arithmetic. The mean can be computed correctly using a larger int type using the trivial sum/count method, but your implementation is less likely to overflow in the general case. I don’t know what to recommend there, if you need to keep these generic, sorry. – Cris Luengo Dec 8 '18 at 17:04
• @user14717: Also sorry for the short review, I got called away while writing it, and decided to just post what I had so far. Maybe I’ll have more time tonight to look at the other functions and your questions, but maybe you prefer to rethink the type logic and post a new question? – Cris Luengo Dec 8 '18 at 17:06
• Really useful: It's clear that this is almost totally useless for image processing guys who use 8/16 bit int logic. Unfortunately, algorithm design depends on whether the data is fixed or floating point, so I don't see any way I could generically support both. I'll think about this more. – user14717 Dec 8 '18 at 17:37
• @user14717: Yes, please don’t update your question. But feel free to post a new question with updates code. I will take another look at your code and the specific questions you had if I have time tonight. – Cris Luengo Dec 8 '18 at 20:18
Real i = 0;
for(auto it = first; it != last; ++it) {
mu = mu + (*it - mu)/i;
i += 1;
}


Unsure what your Real type is, technically, for sufficiently large containers, i in this loop, and similar ones, will soon stop changing.

l[02p]_norms would probably be better expressed via range-for or folds from <algorithm> (std::accumulate etc.)

• typedef typename std::remove_const<typename std::remove_reference<decltype(*std::declval<ForwardIterator>())>::type>::type Real; declared a few lines earlier – bruglesco Dec 9 '18 at 14:42
• Ah, yes, thank you. So could be, say, float thus growing no more past 2²⁴. – bipll Dec 9 '18 at 15:37