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

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单独累计的尾数？

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可能会溢出,您可以将其定义为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;
}

Answer 2

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

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对产品中恒定数量的术语的想法可以避免分支.问题在于溢出/下溢可能仍然只在几个术语内发生,具体取决于值.