对单精度(浮点)值求和时的错误传播
Error Propagation upon Summing Single-Precision (float) Values
我正在学习单精度并想了解错误传播。根据this nice website,加法是危险的操作。
所以我写了一个小的 C 程序来测试错误累积的速度。我不完全确定这是否是一种有效的测试方法。如果是,我不确定如何解释结果,请参见下文。
#include <stdio.h>
#include <math.h>
#define TYPE float
#define NUM_IT 168600
void increment (TYPE base, const TYPE increment, const unsigned long num_iters) {
TYPE err;
unsigned long i;
const TYPE ref = base + increment * num_iters;
for (i=0; i < num_iters; i++ ) {
base += increment;
}
err = (base - ref)/ref;
printf("%lu\t%9f\t%9f\t%+1.9f\n", i, base, ref, err);
}
int
main()
{
int j;
printf("iters\tincVal\trefVal\trelErr\n");
for (j = 1; j < 20; j++ ) {
increment(1e-1, 1e-6, (unsigned long) (pow(2, (j-10))* NUM_IT));
}
return 0;
}
执行结果
gcc -pedantic -Wall -Wextra -Werror -lm errorPropagation.c && ./a.out | tee float.dat | column -t
是
iters incVal refVal relErr
329 0.100328 0.100329 -0.000005347
658 0.100657 0.100658 -0.000010585
1317 0.101315 0.101317 -0.000021105
2634 0.102630 0.102634 -0.000041596
5268 0.105259 0.105268 -0.000081182
10537 0.110520 0.110537 -0.000154624
21075 0.121041 0.121075 -0.000282393
42150 0.142082 0.142150 -0.000480946
84300 0.184163 0.184300 -0.000741986
168600 0.268600 0.268600 +0.000000222 <-- *
337200 0.439439 0.437200 +0.005120996
674400 0.781117 0.774400 +0.008673230
1348800 1.437150 1.448800 -0.008041115
2697600 2.723466 2.797600 -0.026499098
5395200 5.296098 5.495200 -0.036231972
10790400 10.441361 10.890400 -0.041232508
21580800 25.463778 21.680799 +0.174485177
43161600 32.000000 43.261597 -0.260313928 <-- **
86323200 32.000000 86.423195 -0.629729033
如果测试有效
- 为什么错误会改变符号?如果
0.1
表示为例如0.100000001
,无论求和次数如何,这不应该总是累积到相同的偏差吗?
168600
求和有什么特别之处(参见 *
)?误差变得很小。可能是巧合。
- 在
incVal = 32.00
处撞到了哪堵墙(参见 **
,最后两行)。我仍然远低于 unsigned long
限制。
在此先感谢您的努力。
回答您的问题...
1 - IEEE 浮点数舍入为 even 尾数。这样做是专门为了防止错误累积总是以一种或另一种方式产生偏差;如果它总是向下舍入或向上舍入,你的错误会大得多。
2 - 168600 本身并没有什么特别之处。我还没有计算出来,但它很可能最终在二进制表示中产生更清晰的值(即 rational/non-repeating 值)。查看二进制而不是十进制的值,看看该理论是否成立。
3 - 限制因素可能是由于浮点尾数为 23 位长。一旦 base
达到一定大小,increment
与 base
相比是如此之小以至于计算 base + increment
然后将尾数四舍五入到 23 位完全擦除变化.即 base
和 base + increment
之间的差异是舍入误差。
首先,重要的是要知道 0.1
不能准确表示,在二进制中它有周期性重复的数字。该值为 0.0001100110011...
。比较 1/3 和 1/7 是如何用十进制数字表示的。值得用增量 0.25
重复您的测试,它可以精确表示为 0.01
.
我会用十进制来说明错误,这是我们人类习惯的。让我们使用小数,并假设我们可以有 4 位精度。这些就是这里发生的事情。
除法:让我们计算1/11:
1/11 等于 0.090909...,这可能四舍五入为 0.09091。正如预期的那样,这正确到 4 位有效数字(粗体)。
量级差异:假设我们计算 10 + 1/11.
将1/11加到10时,我们必须做更多的舍入,因为10.09091是7位有效数字,而我们只有四位。我们要把1/11四舍五入到小数点后的两位数,计算出来的和是10.09。这是低估了。请注意如何仅保留 1/11 的一位有效数字。如果将很多小值加在一起,这将限制最终结果的精度。
现在计算 100 + 1/11。现在我们将 1/11 取整为 0.1,并将总和表示为 100.1。现在我们有轻微的高估而不是轻微的低估。
我猜你测试中的符号变化模式是系统性轻微低估与高估的影响,具体取决于 base
.
的大小
1000 + 1/11 呢?现在我们不能在点之后有任何数字,因为我们在点之前已经有 4 个有效数字。 1/11 现在四舍五入为 0,总和仍然是 1000。那就是你看到的墙。
您在测试中没有看到的另一件重要事情是:如果两个值的符号不同会发生什么。计算 1.234 – 1.243:两个数字都有 4 位有效数字。结果是 -0.009。现在结果只有一位正确的有效数字,而不是四位。
此处类似问题的答案:How does floating point error propagate when doing mathematical operations in C++?。它有一些指向更多信息的链接。
你打的"wall"与增量值无关,如果它通过加法是恒定的并且你从零开始。它必须与iters
。 2^23 = 800万,你做的是8600万次加法。所以一旦累加器比增量大 2^23,你就会碰壁。
尝试 运行 具有 86323200 次迭代的代码,但增量为 1 或 0.0000152587890625(或 2 的任何幂)。它应该与增量32具有相同的相对问题。
我正在学习单精度并想了解错误传播。根据this nice website,加法是危险的操作。
所以我写了一个小的 C 程序来测试错误累积的速度。我不完全确定这是否是一种有效的测试方法。如果是,我不确定如何解释结果,请参见下文。
#include <stdio.h>
#include <math.h>
#define TYPE float
#define NUM_IT 168600
void increment (TYPE base, const TYPE increment, const unsigned long num_iters) {
TYPE err;
unsigned long i;
const TYPE ref = base + increment * num_iters;
for (i=0; i < num_iters; i++ ) {
base += increment;
}
err = (base - ref)/ref;
printf("%lu\t%9f\t%9f\t%+1.9f\n", i, base, ref, err);
}
int
main()
{
int j;
printf("iters\tincVal\trefVal\trelErr\n");
for (j = 1; j < 20; j++ ) {
increment(1e-1, 1e-6, (unsigned long) (pow(2, (j-10))* NUM_IT));
}
return 0;
}
执行结果
gcc -pedantic -Wall -Wextra -Werror -lm errorPropagation.c && ./a.out | tee float.dat | column -t
是
iters incVal refVal relErr
329 0.100328 0.100329 -0.000005347
658 0.100657 0.100658 -0.000010585
1317 0.101315 0.101317 -0.000021105
2634 0.102630 0.102634 -0.000041596
5268 0.105259 0.105268 -0.000081182
10537 0.110520 0.110537 -0.000154624
21075 0.121041 0.121075 -0.000282393
42150 0.142082 0.142150 -0.000480946
84300 0.184163 0.184300 -0.000741986
168600 0.268600 0.268600 +0.000000222 <-- *
337200 0.439439 0.437200 +0.005120996
674400 0.781117 0.774400 +0.008673230
1348800 1.437150 1.448800 -0.008041115
2697600 2.723466 2.797600 -0.026499098
5395200 5.296098 5.495200 -0.036231972
10790400 10.441361 10.890400 -0.041232508
21580800 25.463778 21.680799 +0.174485177
43161600 32.000000 43.261597 -0.260313928 <-- **
86323200 32.000000 86.423195 -0.629729033
如果测试有效
- 为什么错误会改变符号?如果
0.1
表示为例如0.100000001
,无论求和次数如何,这不应该总是累积到相同的偏差吗? 168600
求和有什么特别之处(参见*
)?误差变得很小。可能是巧合。- 在
incVal = 32.00
处撞到了哪堵墙(参见**
,最后两行)。我仍然远低于unsigned long
限制。
在此先感谢您的努力。
回答您的问题...
1 - IEEE 浮点数舍入为 even 尾数。这样做是专门为了防止错误累积总是以一种或另一种方式产生偏差;如果它总是向下舍入或向上舍入,你的错误会大得多。
2 - 168600 本身并没有什么特别之处。我还没有计算出来,但它很可能最终在二进制表示中产生更清晰的值(即 rational/non-repeating 值)。查看二进制而不是十进制的值,看看该理论是否成立。
3 - 限制因素可能是由于浮点尾数为 23 位长。一旦 base
达到一定大小,increment
与 base
相比是如此之小以至于计算 base + increment
然后将尾数四舍五入到 23 位完全擦除变化.即 base
和 base + increment
之间的差异是舍入误差。
首先,重要的是要知道 0.1
不能准确表示,在二进制中它有周期性重复的数字。该值为 0.0001100110011...
。比较 1/3 和 1/7 是如何用十进制数字表示的。值得用增量 0.25
重复您的测试,它可以精确表示为 0.01
.
我会用十进制来说明错误,这是我们人类习惯的。让我们使用小数,并假设我们可以有 4 位精度。这些就是这里发生的事情。
除法:让我们计算1/11:
1/11 等于 0.090909...,这可能四舍五入为 0.09091。正如预期的那样,这正确到 4 位有效数字(粗体)。
量级差异:假设我们计算 10 + 1/11.
将1/11加到10时,我们必须做更多的舍入,因为10.09091是7位有效数字,而我们只有四位。我们要把1/11四舍五入到小数点后的两位数,计算出来的和是10.09。这是低估了。请注意如何仅保留 1/11 的一位有效数字。如果将很多小值加在一起,这将限制最终结果的精度。
现在计算 100 + 1/11。现在我们将 1/11 取整为 0.1,并将总和表示为 100.1。现在我们有轻微的高估而不是轻微的低估。
我猜你测试中的符号变化模式是系统性轻微低估与高估的影响,具体取决于
base
. 的大小
1000 + 1/11 呢?现在我们不能在点之后有任何数字,因为我们在点之前已经有 4 个有效数字。 1/11 现在四舍五入为 0,总和仍然是 1000。那就是你看到的墙。
您在测试中没有看到的另一件重要事情是:如果两个值的符号不同会发生什么。计算 1.234 – 1.243:两个数字都有 4 位有效数字。结果是 -0.009。现在结果只有一位正确的有效数字,而不是四位。
此处类似问题的答案:How does floating point error propagate when doing mathematical operations in C++?。它有一些指向更多信息的链接。
你打的"wall"与增量值无关,如果它通过加法是恒定的并且你从零开始。它必须与iters
。 2^23 = 800万,你做的是8600万次加法。所以一旦累加器比增量大 2^23,你就会碰壁。
尝试 运行 具有 86323200 次迭代的代码,但增量为 1 或 0.0000152587890625(或 2 的任何幂)。它应该与增量32具有相同的相对问题。