计算多个数字的几何平均值的高效方法

Question

我需要计算一组数字的几何平均值，这些数字的值未事先限制。朴素的方法是

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());
}

但是，这种方法可能会因为在累积的 product 中出现下溢或上溢而失败（注意：即使使用 long double 也无法完全避免此问题）。因此，下一个选择是对对数进行求和：

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());
}

这个方法可以工作，但是对于每一个元素都调用了std::log()，可能会导致速度变慢。有没有办法避免这种情况？例如通过分别跟踪积累的product的指数和尾数（相当于）？

- Walter

---

数字的值是否有任何最大范围？ - Abhishek Bansal

@Walter：你尝试过使用长双精度吗？ - ahmedsafan86

@Ahmedsafan 如果有太多数字，就不会有帮助：请考虑范围在0.01到1之间的数字以及其中的10^6个... - Walter

你有没有检查过log()相对于其他替代方案的速度有多慢？我很想看到一个实际的性能比较... - comingstorm

对于double类型，可以非常高效地实现对数运算。你可能真的想在这里进行一次比较，结果会非常有趣 :) - filmor

Answer 1

“分离指数和尾数”的解决方案：

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);
}

如果您担心ex可能会溢出，可以将其定义为双精度而不是long long，并在每个步骤中乘以invN，但是这种方法可能会失去很多精度。

编辑对于大输入，我们可以将计算分成几个存储桶：

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

---

@Walter 我只进行了几项测试，没有尝试一些明显可能失败的边界情况。例如，如果你有超过约1022个数字，m就可能下溢。另外，经过再次考虑，ex实际上不会发生溢出（需要类似于10^16个输入才有可能发生溢出）。 - sbabbi

@stabbi确实，您可能仍会发生下溢，特别是因为始终f1<1。所以，经过再次考虑，我不应该接受这个...也许您可以解决这个问题。在实践中，我有超过1022个数字... - Walter

你仍然会对每个桶执行std::pow操作，而我认为这比使用std::log慢，所以肯定有提升的空间... - Walter

@AlexandreC。是这样吗？我没有计算各种选项的时间，但我认为frexp应该很快，因为它基本上不需要计算：它只是将参数分成尾数和指数。 - Walter

@AlexandreC。你是对的。使用GCC编译时，日志在使用-O2编译时更快。不确定为什么，但我尝试编写了一个定制的非便携式frexp，实际上非常快。我将感激任何关于这个事实的见解 :) - sbabbi

Answer 2

我认为我找到了一个方法来做到这一点，它结合了问题中的两个例程，类似于Peter的想法。以下是示例代码。

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());
}

坏消息是：这带来了一个分支。好消息是：分支预测器很可能几乎总是正确的（分支应该只在很少情况下被触发）。

可以使用彼得的想法，即产品中有固定数量的项，来避免分支。问题在于，在只有几个术语的情况下，根据其值可能仍会发生溢出/下溢。

- Walter

---

如果其中一个值非常大或小（>1e244或<1e-244），那么这仍然可能无法起作用。 - Jeffrey Sax

当然可以，但挑战在于即使在困难的情况下也要给出最佳结果。 - Jeffrey Sax

Answer 3

您可以通过像原始解决方案中那样乘以数字并仅在一定数量的乘法后（取决于您最初的数字的大小）才转换为对数来加速此过程。

- Peter de Rivaz

Answer 4

另一种比对数方法更准确和高效的方法是通过固定量抵消超出范围的指数，保持已取消多余部分的精确对数。具体实现如下：

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);

所有由BIG和SMALL进行的乘法都是精确的，且没有调用log（一个超越函数，因此特别不精确）。

- Sneftel

---

是的，我也考虑过这种方法。但是，我想所有这些算术运算可能会使它变得相当慢，所以“log”函数实际上可能更快。对数方法的准确性对我来说完全没问题。 - Walter

所有的算术运算？整数的增减和浮点数的乘法是处理器可以执行的最快速度之一。我保证这比对数方法更快。 - Sneftel

product *= A[i] 这里仍然可能会溢出。此外，while 循环将变成无限循环。 - ZachB

Answer 5

有一个简单的想法可以减少计算量，同时防止溢出。您可以将数字分组，至少每次两个，并计算它们的对数，然后计算它们的总和。

log(abcde) = 5*log(K)

log(ab) + log(cde)  = 5*log(k)

- Vikram Bhat

---

是的，这基本上与彼得的想法相同，但只涉及两个数字。 - Walter

Answer 6

不需要使用代价高昂的对数，可以直接通过二次幂来缩放结果。

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());
}

- Jeffrey Sax

---

看起来很熟悉。;-) 注意最终的对数并不是必要的，你只是在之后撤销它。 - Sneftel

是的，非常相似。不过请注意两件事情：首先，pow 隐式地调用了 log，所以你并没有节省任何东西。其次，如果任何一个值非常大或非常小（>DBL_MAX/BIG 或 <DBL_MIN/SMALL），你的结果将不正确。 - Jeffrey Sax

Answer 7

对于计算乘积，通过对数求和是完全可行的，而且相当高效（如果这还不够：有一些方法可以使用几个SSE操作获得向量化的对数 - 还有英特尔MKL的向量操作）。

为了避免溢出，常见的技术是事先将每个数字除以最大或最小幅度条目（或将对数差值相加到对数最大值或对数最小值）。如果数字变化很大，也可以使用桶（例如，分别对小数和大数的对数求和）。请注意，通常除非是非常大的集合，否则不需要这两种方法，因为double类型的对数永远不会很大（在-700到700之间）。

此外，您需要单独跟踪符号。

计算log x通常保留与x相同数量的有效数字，除非x接近1：如果需要计算prod（1 + x_n），其中x_n很小，则应使用std :: log1p。

最后，如果您在求和时遇到舍入误差问题，可以使用Kahan summation或其变体。

- Alexandre C.