floating-point - 几何平均数的安全计算

Question

我正在寻找一个理由来选择以下方法之一来计算一长串 float x 的几何平均值:

1. 取每个x的n次方根，然后将它们相乘
2. 将它们全部相乘，然后取n次根

我听说对于 float ，乘法和除法丢失的信息少于加法和减法。因此，我不考虑求和技巧。

我应该通过 1 还是 2 来计算几何平均值，为什么？

---

更新 1，回应评论:

所有 x 均小于 1 且为 double 。它们的数量级在 10^-1 到 10^-6 之间。请假设计算 n 次根的最常用方法，因为我使用的是编程语言的内置函数。我担心的不是溢出，而是下溢(？)，因为所有值都小于 1。您可以假设 x 系列的长度为 10^8

Answer 1

一般而言，在还涉及平方根或立方根等收缩运算的浮点运算序列中，从精度的角度来看，最后执行收缩运算是有利的。例如sqrt(1.0/x)比1.0/sqrt(x)更准确，sqrt(a*b)更准确比 sqrt(a)*sqrt(b)，cbrt(a*b*c) 比 cbrt(a)*cbrt(b) 更准确*cbrt(c).

因此，除非存在上溢或下溢所选浮点格式的危险，例如 IEEE-754 binary64(例如 C/C++ 中的 double ), 在中间计算时，应选择方法[2]。与准确性相关的其他方面:如果第 n 个根是通过求幂计算的，例如 C/C++ 中的 pow()，则每个计算的根都会引入额外的错误，正如我对 this question 的回答中关于立方根的解释.最后，第 n 个根的计算将比乘法慢，因此最后只进行最终根计算的乘法也是一种性能优越的方法。

通过使用补偿乘积(类似于 Kahan summation 提供的补偿加法)，方法 [2] 可以获得非常准确的结果。有关详细信息，请参阅以下论文:

Stef Graillat，“精确的浮点乘积和求幂”，IEEE 计算机交易，卷。 58, No. 7, July 2009, pp. 994-1000 ( online )

在硬件中提供 FMA(融合乘加)运算的系统上，可以特别有效地计算此补偿乘积。所有常见的现代处理器架构(包括 CPU 和 GPU)都是这种情况。 C/C++ 通过标准数学函数 fma()、fmaf() 提供对此的便捷访问。

更新:提问者在评论中澄清说，下溢的风险迫在眉睫，因为 [10^-6, 10^-1]. @Yves Daoust 在评论中提到的一种可能的解决方法是将因子分成尾数和指数并分别累加。这是否可行取决于浮点环境。虽然 C 和 C++ 提供了执行此拆分的标准函数 frexp()，但此函数可能不是很快。