浮点加法：精度损失问题答案

【问题标题】：Floating point addition: loss-of-precision issues浮点加法：精度损失问题
【发布时间】：2009-08-10 07:58:17
【问题描述】：

简而言之：我如何执行a+b 以使由于截断而导致的任何精度损失远离零而不是接近零？

长篇大论

我正在计算一长串浮点值的总和，以计算集合的样本均值和方差。由于 Var(X) = E(X²) - E(X)²，保持所有数字的连续计数就足够了，到目前为止所有数字的总和，以及迄今为止所有数字的平方和。

到目前为止一切顺利。

但是，E(X²) > E(X)² 是绝对必要的，因为浮点精度不是总是如此。在伪代码中，问题是这样的：

int count;
double sum, sumOfSquares;
...
double value = <current-value>;
double sqrVal = value*value; 

count++;
sum += value; //slightly rounded down since value is truncated to fit into sum
sumOfSquares += sqrVal; //rounded down MORE since the order-of-magnitude 
//difference between sqrVal and sumOfSquares is twice that between value and sum;

对于可变序列，这不是一个大问题 - 您最终会稍微低估方差，但这通常不是一个大问题。但是，对于具有非零均值的常数或几乎常数集，它可能意味着 E(X²) 2，导致计算出的方差为负，这违反了使用代码的预期。

现在，我知道了 Kahan Summation，这不是一个有吸引力的解决方案。首先，它使代码容易受到优化变幻莫测的影响（取决于优化标志，代码可能会或可能不会出现此问题），其次，由于精度，问题不是真的 - 这很好够了 - 这是因为加法将 systematic 错误引入零。如果我可以执行该行

sumOfSquares += sqrVal;

以确保 sqrVal 向上而不是向下舍入到 sumOfSquares 的精度的方式，我会有一个数值上合理的解决方案。但是我怎样才能做到这一点呢？

编辑：已完成的问题 - 为什么在标签字段的下拉列表中按 enter 无论如何都会提交问题？

【问题讨论】：

对集合进行排序并从较小的值开始进行相同的计算会改变这种情况吗？
对于当前解决方案，排序的计算成本要高得多，需要 O(n log(n)) 时间和 O(n) 存储，而不是线性时间（以及低得多的常数因子）和常数存储。我处理的数据集任意大（越大越好），所以这是个问题。
毫无疑问，但这首先会有所帮助吗？
这无济于事，因为即使在完全恒定的集合中也会出现问题 - 只要平方和和总和有差异截断错误（它们会），不幸的系列可能会导致负方差。

标签： c# c++ floating-point ieee-754

【解决方案1】：

还有另一种单程算法可以稍微重新安排计算。在伪代码：

n = 0
mean = 0
M2 = 0

for x in data:
    n = n + 1
    delta = x - mean
    mean = mean + delta/n
    M2 = M2 + delta*(x - mean)  # This expression uses the new value of mean

variance_n = M2/n         # Sample variance
variance = M2/(n - 1)     # Unbiased estimate of population variance

（来源：http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance）

就您指出的问题而言，这似乎表现得更好使用通常的算法。

【讨论】：

也就是说，我应该检查一下维基百科 ;-)。谢谢，看起来很有希望！
...而且维基百科甚至有一个加权版本，这是我真正追求的，但我认为我不会不必要地混淆水域。

【解决方案2】：

IEEE 提供四种舍入模式（朝向 -inf、朝向 +inf、朝向 0、最接近）。朝着 +inf 方向似乎是您想要的。 C90 或 C++ 中没有标准控件。 C99 添加了标题<fenv.h>，它在某些 C90 和 C++ 实现中也作为扩展存在。要遵守 C99 标准，您必须编写如下内容：

#include <fenv.h>
#pragma STDC FENV_ACCESS ON

int old_round_mode = fegetround();
int set_round_ok = fesetround(FE_UPWARD);
assert(set_round_ok == 0);
...
int set_round_ok = fesetround(old_round_mode);
assert(set_round_ok == 0);

众所周知，您使用的算法在数值上不稳定并且存在精度问题。对数据进行两次传递以提高精度。

【讨论】：

由于性能问题，使用两次传递确实很不幸（这也使 API 更丑陋）。据我所知，如果你只是四舍五入，算法应该是稳定的 - 对吧？
我想知道，像 "sumOfSquares += sqrVal + sumOfSquares/(1L
@Eamon，关于你的第一个问题，我没有时间做真正的稳定性分析。尤其是我不经常这样做，对结果充满信心。您的第二条评论中的代码似乎根本不等效（您是否打算改为划分 sqrVal？在这种情况下，缩放不会改变稳定性或精度）。
不，我想要 sumOfSquares。动机：double 的精度为 52 位，因此第 53 位是潜在的错误来源。为了确保估计值永远不会太高而不是太低，我也可以简单地添加第 53 位。大概 sqlVal 小到足以包含该位，然后我确信任何舍入误差都安全地低于 1/2^52 阈值。

【解决方案3】：

如果您不担心精度，而只担心负方差，您为什么不简单地做V(x) = Max(0, E(X^2) - E(X)^2)

【讨论】：

这是我最初的解决方法，但我希望利用 stackoverflow 的丰富智慧获得更好的解决方法。这是一个务实的解决方案 - 可能应该提到它;-)。