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

小编典典

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

algorithm

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

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？

阅读 667

2020-07-28

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

2020-07-28