计算多个数字的几何平均值的有效方法
我需要计算一大组数字的几何平均值，其值不是先验有限的。 天真的方式是
double geometric_mean(std::vector 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 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的指数和尾数（相当于）？
“分裂指数和尾数”解决方案：
double geometric_mean(std::vector 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 ::radix,ex * invN) * std::pow(m,invN); }
如果你担心 ex可能会溢出，你可以将它定义为double而不是long long ，并在每一步都乘以invN ，但是这种方法可能会失去很多精度。
编辑对于大型输入，我们可以将计算分成几个桶：
double geometric_mean(std::vector 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 ::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 ::radix,ex * invN ) * m; }
我想我找到了一种方法，它结合了问题中的两个例程，类似于彼得的想法。 这是一个示例代码。
double geometric_mean(std::vector 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对产品中恒定数量的术语的想法可以避免分支。 问题在于溢出/下溢可能仍然只在几个术语内发生，具体取决于值。
您可以通过将原始解决方案中的数字相乘来加速此操作，并且只能在每个特定数量的乘法中转换为对数（取决于初始数字的大小）。
与对数方法相比，提供更好的准确性和性能的不同方法是以固定量补偿超出范围的指数，保持取消的超额的精确对数。 像这样：
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 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 （一个超越的，因此特别不精确的函数）。
有简单的想法减少计算并防止溢出。 您可以将数字组合在一起，至少说两个，并计算他们的日志，然后评估他们的总和。
log(abcde) = 5*log(K) log(ab) + log(cde) = 5*log(k)
总结日志以稳定地计算产品是非常好的，并且非常有效（如果这还不够：有一些方法可以通过一些SSE操作获得矢量化对数 – 还有英特尔MKL的矢量运算）。
为了避免溢出，一种常见的技术是预先将每个数字除以最大或最小幅度条目（或将日志差异与log max或log min相加）。 如果数字变化很大，您也可以使用存储桶（例如，分别将小数字和大数字的对数相加）。 请注意，除了非常大的集合之外，通常都不需要这样做，因为double的日志永远不会很大（介于-700和700之间）。
此外，您需要单独跟踪标志。
计算log x通常保留与log x相同的有效位数，除非x接近1 ：如果需要计算带有小x_n prod(1 + x_n) ，则需要使用std::log1p 。
最后，如果在求和时遇到舍入误差问题，则可以使用Kahan求和或变量。
您可以直接使用2的幂来缩放结果，而不是使用非常昂贵的对数。 在伪代码中：
double huge = scalbn(1,512); double tiny = scalbn(1,512); int scale = 0; for(auto x:data) { if (x >= huge) { scalbn(x, -512); scale++; } else if (x <= tiny) { scalbn(x, 512); scale--; } product *= x; if (product >= huge) { scalbn(product, -512); scale++; } else if (product <= tiny) { scalbn(product, 512); scale--; } return exp2((512.0*scale + log2(product)) / data.size()); }