30

I need to compute the geometric mean of a large set of numbers, whose values are not a priori limited. The naive way would be

double geometric_mean(std::vector<double> const&data) // failure
{
  auto product = 1.0;
  for(auto x:data) product *= x;
  return std::pow(product,1.0/data.size());
}

However, this may well fail because of underflow or overflow in the accumulated product (note: long double doesn't really avoid this problem). So, the next option is to sum-up the logarithms:

double geometric_mean(std::vector<double> const&data)
{
  auto sumlog = 0.0;
  for(auto x:data) sum_log += std::log(x);
  return std::exp(sum_log/data.size());
}

This works, but calls std::log() for every element, which is potentially slow. Can I avoid that? For example by keeping track of (the equivalent of) the exponent and the mantissa of the accumulated product separately?

Walter
  • 44,150
  • 20
  • 113
  • 196

7 Answers7

13

The "split exponent and mantissa" solution:

double geometric_mean(std::vector<double> const & data)
{
    double m = 1.0;
    long long ex = 0;
    double invN = 1.0 / data.size();

    for (double x : data)
    {
        int i;
        double f1 = std::frexp(x,&i);
        m*=f1;
        ex+=i;
    }

    return std::pow( std::numeric_limits<double>::radix,ex * invN) * std::pow(m,invN);
}

If you are concerned that ex might overflow you can define it as a double instead of a long long, and multiply by invN at every step, but you might lose a lot of precision with this approach.

EDIT For large inputs, we can split the computation in several buckets:

double geometric_mean(std::vector<double> const & data)
{
    long long ex = 0;
    auto do_bucket = [&data,&ex](int first,int last) -> double
    {
        double ans = 1.0;
        for ( ;first != last;++first)
        {
            int i;
            ans *= std::frexp(data[first],&i);
            ex+=i;
        }
        return ans;
    };

    const int bucket_size = -std::log2( std::numeric_limits<double>::min() );
    std::size_t buckets = data.size() / bucket_size;

    double invN = 1.0 / data.size();
    double m = 1.0;

    for (std::size_t i = 0;i < buckets;++i)
        m *= std::pow( do_bucket(i * bucket_size,(i+1) * bucket_size),invN );

    m*= std::pow( do_bucket( buckets * bucket_size, data.size() ),invN );

    return std::pow( std::numeric_limits<double>::radix,ex * invN ) * m;
}
sbabbi
  • 11,070
  • 2
  • 29
  • 57
  • @Walter I made only a few tests, I did not try some obvious border-cases that can fail. For example `m` can underflow if you have more than ~1022 numbers. Also, on second thought, `ex` won't realistically overflow (need something like 10^16 inputs to have a chance of overflowing). – sbabbi Nov 14 '13 at 16:50
  • 1
    @stabbi Indeed, you may still underflow, in particular since always `f1`<1. So, on second thought, I should not have accepted this ... Perhaps you can fix this problem. In practice, I have many more than 1022 numbers ... – Walter Nov 14 '13 at 18:33
  • @Walter fixed and tested with input sizes up to 10^8. – sbabbi Nov 14 '13 at 18:58
  • you still do a `std::pow` per bucket, which I think it slower than `std::log`, so there must be potential for improvement ... – Walter Nov 14 '13 at 19:01
  • @Walter you can do it in the "adaptive" way like in your solution. Also consider that a bucket is quite large (~1000 elements). Moreover, this solution can be easily parallelized (either with openmp or by replacing the first for loop with a `parallel_for`), you just have to protect `ex` with std::atomic. – sbabbi Nov 14 '13 at 19:06
  • On MSVC `frexp` and `ldexp` are typically slower than logarithms. – Alexandre C. Nov 20 '13 at 00:44
  • @sbabbi Isn't `std::numeric_limits::max_exponent` the same as `-std::log2( std::numeric_limits::min()` (if the radix is 2)? – Walter Nov 20 '13 at 11:28
  • @AlexandreC. Is this so? I havn't timed the various options, but I would think that `frexp` should be fast, as it essentially needs no computation: it just splits the argument into mantissa and exponent. – Walter Nov 20 '13 at 11:29
  • @AlexandreC. You are right. With GCC log is faster when compiling with -O2. Not sure why, but I tried to write a custom-made nonportable frexp which is actually very fast. I'll appreciate any insight on this fact :) – sbabbi Nov 20 '13 at 15:38
11

I think I figured out a way to do it, it combined the two routines in the question, similar to Peter's idea. Here is an example code.

double geometric_mean(std::vector<double> const&data)
{
    const double too_large = 1.e64;
    const double too_small = 1.e-64;
    double sum_log = 0.0;
    double product = 1.0;
    for(auto x:data) {
        product *= x;
        if(product > too_large || product < too_small) {
            sum_log+= std::log(product);
            product = 1;      
        }
    }
    return std::exp((sum_log + std::log(product))/data.size());
}

The bad news is: this comes with a branch. The good news: the branch predictor is likely to get this almost always right (the branch should only rarely be triggered).

The branch could be avoided using Peter's idea of a constant number of terms in the product. The problem with that is that overflow/underflow may still occur within only a few terms, depending on the values.

Walter
  • 44,150
  • 20
  • 113
  • 196
4

You may be able to accelerate this by multiplying numbers as in your original solution and only converting to logarithms every certain number of multiplications (depending on the size of your initial numbers).

Peter de Rivaz
  • 33,126
  • 4
  • 46
  • 75
3

A different approach which would give better accuracy and performance than the logarithm method would be to compensate out-of-range exponents by a fixed amount, maintaining an exact logarithm of the cancelled excess. Like so:

const int EXP = 64; // maximal/minimal exponent
const double BIG = pow(2, EXP); // overflow threshold
const double SMALL = pow(2, -EXP); // underflow threshold

double product = 1;
int excess = 0; // number of times BIG has been divided out of product

for(int i=0; i<n; i++)
{
    product *= A[i];
    while(product > BIG)
    {
        product *= SMALL;
        excess++;
    }
    while(product < SMALL)
    {
        product *= BIG;
        excess--;
    }
}

double mean = pow(product, 1.0/n) * pow(BIG, double(excess)/n);

All multiplications by BIG and SMALL are exact, and there's no calls to log (a transcendental, and therefore particularly imprecise, function).

Sneftel
  • 40,271
  • 12
  • 71
  • 104
  • Yes, I thought about this approach too. However, I would imagine that all this arithmetic makes is rather slow, so that the `log` function may actually be faster. The accuracy of the logarithm method is perfectly okay for me. – Walter Nov 14 '13 at 16:35
  • 2
    All the arithmetic? Integer increment/decrement and floating-point multiplication are two of the fastest things your processor can do. I guarantee this is faster than the logarithm method. – Sneftel Nov 14 '13 at 17:07
  • The `product *= A[i]` mul can still overflow here. Furthermore, the `while`s then become infinite loops. – ZachB Dec 26 '18 at 06:05
1

There is simple idea to reduce computation and also to prevent overflow. You can group together numbers say atleast two at time and calculate their log and then evaluate their sum.

log(abcde) = 5*log(K)

log(ab) + log(cde)  = 5*log(k)
Vikram Bhat
  • 6,106
  • 3
  • 20
  • 19
1

Instead of using logarithms, which are very expensive, you can directly scale the results by powers of two.

double geometric_mean(std::vector<double> const&data) {
  double huge = scalbn(1,512);
  double tiny = scalbn(1,-512);
  int scale = 0;
  double product = 1.0;
  for(auto x:data) {
    if (x >= huge) {
      x = scalbn(x, -512);
      scale++;
    } else if (x <= tiny) {
      x = scalbn(x, 512);
      scale--;
    }
    product *= x;
    if (product >= huge) {
      product = scalbn(product, -512);
      scale++;
    } else if (product <= tiny) {
      product = scalbn(product, 512);
      scale--;
    }
  }
  return exp2((512.0*scale + log2(product)) / data.size());
}
ZachB
  • 13,051
  • 4
  • 61
  • 89
Jeffrey Sax
  • 10,253
  • 3
  • 29
  • 40
  • Looks familiar. ;-) note that the final logarithm isn't necessary, you're just undoing it afterward. – Sneftel Nov 14 '13 at 15:10
  • 1
    Yes, very similar. Note 2 things though: First, `pow` implicity uses a call to `log` so you're not saving anything. Second, if any of the values are very large or very small (>DBL_MAX/BIG or – Jeffrey Sax Nov 14 '13 at 15:15
1

Summing logs to compute products stably is perfectly fine, and rather efficient (if this is not enough: there are ways to get vectorized logarithms with a few SSE operations -- there are also Intel MKL's vector operations).

To avoid overflow, a common technique is to divide every number by the maximum or minimum magnitude entry beforehand (or sum log differences to the log max or log min). You can also use buckets if the numbers vary a lot (eg. sum the log of small numbers and large numbers separately). Note that typically neither of this is needed except for very large sets since the log of a double is never huge (between say -700 and 700).

Also, you need to keep track of the signs separately.

Computing log x keeps typically the same number of significant digits as x, except when x is close to 1: you want to use std::log1p if you need to compute prod(1 + x_n) with small x_n.

Finally, if you have roundoff error problems when summing, you can use Kahan summation or variants.

Alexandre C.
  • 55,948
  • 11
  • 128
  • 197