29

我需要计算大量数字的几何平均值，其值不受先验限制。天真的方法是

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？

4

7 回答 7

12

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

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可能会溢出，您可以将其定义为 double 而不是 a 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;
}
于 2013-11-14T15:54:54.223 回答
11

我想我想出了一种方法，它结合了问题中的两个例程，类似于彼得的想法。这是一个示例代码。

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

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

使用 Peter 的产品中的常数数量的概念可以避免该分支。问题在于，根据值的不同，上溢/下溢仍可能仅在几个术语内发生。

于 2013-11-14T14:47:39.363 回答
4

您可以通过在原始解决方案中将数字相乘并仅在每一定次数的乘法中转换为对数（取决于初始数字的大小）来加速这一过程。

于 2013-11-14T14:46:03.157 回答
3

与对数方法相比，可以提供更好的精度和性能的另一种方法是将超出范围的指数补偿一个固定的量，同时保持消除的过量的精确对数。像这样：

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（一个超越的，因此特别不精确的函数）。

于 2013-11-14T15:00:46.040 回答
1

对日志求和以稳定地计算产品非常好，而且相当有效（如果这还不够：有一些方法可以通过一些 SSE 操作获得矢量化对数 - 还有英特尔 MKL 的矢量操作）。

为了避免溢出，一种常见的技术是预先将每个数字除以最大或最小幅度条目（或将对数差异与对数最大值或对数最小值相加）。如果数字变化很大，您也可以使用存储桶（例如，分别对小数和大数的对数求和）。请注意，除了非常大的集合之外，通常都不需要这些，因为 a 的日志double永远不会很大（比如在 -700 和 700 之间）。

此外，您需要单独跟踪标志。

计算log x通常保持与 相同数量的有效数字x，除非x接近1:std::log1p如果您需要prod(1 + x_n)使用 small计算，您希望使用x_n。

最后，如果求和时出现舍入误差问题，可以使用Kahan 求和或变体。

于 2013-11-20T00:58:18.547 回答
1

您可以直接通过 2 的幂来缩放结果，而不是使用非常昂贵的对数。

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());
}
于 2013-11-14T15:07:48.827 回答
1

有一个简单的想法可以减少计算并防止溢出。您可以将至少两个数字组合在一起，计算它们的对数，然后评估它们的总和。

log(abcde) = 5*log(K)

log(ab) + log(cde)  = 5*log(k)
于 2013-11-14T14:48:18.683 回答