1

我开发了一个图像处理软件,我需要对其进行数值分析,考虑到与其操作相关的错误传播以及浮点类型变量的不确定性,这些不确定性是由此类变量发生的固有舍入引起的。

考虑到 IEEE 754 标准machine epsilon,浮点类型变量是1.19e-07. 据我了解,这个值是到最近的可表示浮点的距离。

我做了一些测试,通过向这个 epsilon 添加一个浮点值来确定这是否是真的x + epsilon == x:这个概念不适用于浮点范围的每个值,这是可以理解的,因为浮点数的大值具有更多的不确定性,这是由舍入和用于表示它们的有限位数引起的。

我的问题是与浮点值相关的不确定性是什么,即(x + y) || (x - y) == x浮点值x和浮点不确定性y

可能是我对英语缺乏了解,但我似乎无法理解有关该主题的文献。

如果有人可以尽可能详细,您能否通过以下简单操作向我解释错误?

float result = valA * 0.587f + valB * 0.331f;

如果我知道浮点类型变量的不确定性,这个错误可以简单地用这个公式计算出来,对吧?

在此处输入图像描述

4

1 回答 1

3

介绍

该答案对以下错误进行了初步检查:

float result = valA * 0.587f + valB * 0.331f;

在这个答案中,浮点格式的值和使用浮点格式计算的表达式将用代码样式表示,如zor x * y。数学变量将使用斜体并且不会采用代码样式,如zxy

我假设所有算术都是使用 IEEE-754 基本 32 位二进制浮点数完成的。这种格式通常用于float类型,尽管一些编程语言实现混合精度,可能double在评估float类型表达式时使用或其他精度。我还假设所有算术都是使用四舍五入模式完成的,并与偶数低位的数字相关联。

这种格式为 24 位的有效位,因此最低精度单位 (ULP) 通常是最高有效位值的2 -23倍。这是可表示值之间的步长。例如,对于 [1, 2) 中的值,ULP 为 2 -23。对于 [128, 256) 中的值,ULP 为 2 7 •2 -23 = 2 -16。(对于次正规值,有效数字的位数较少。ULP 的最低值可以是 2 -149。超出最大的有限可表示值,到下一个可表示值的步长是无限的。但是,在这个问题中,只有适度的值涉及值,因此我们可以忽略无穷大。)

使用正确舍入计算任何运算的结果最多与正确答案相差 ½ ULP。也就是说,z = x + y例如,如果我们计算 ,计算结果z与精确的数学结果z = x+最多相差zy的 ½ ULP 。(虽然z是一个具有无限精度的精确数学结果,但我们使用它的大小来确定它在浮点格式中落入哪个范围,这决定了我们所说的z的 ULP 的含义。)误差最多为 ½ 的原因ULP 是,如果最接近z的两个可表示值是和,我们必须有≤ z ≤ ,并且如果 ½ ULP <z0z1z0z1z1z,然后zz0< ½ ULP(因为z1z0= 1 ULP,根据 ULP 的定义。)因此,在选择最接近的可表示值时,我们会选择 和 的更接近者z0z1因此误差永远不会超过 ½ ULP。

如评论中所述valAvalB、 和result位于 [0, 256) 中。

符号分析

当我们开始计算valA * 0.587f + valB * 0.331f时,之前valAvalB操作会出现一些错误。也就是说,理想情况下,使用精确的数学,我们会计算出一些数字AB,但计算机计算的是valAvalB,差异是eA = valA- AeB = valB- B

理想情况下,我们希望使用精确的数学计算R使得R为A • .587 + B • .331。当我们使用计算机算术时:

  • 0.587f将从 .587 转换为浮点格式,结果会有一些舍入误差e0,所以结果是0.587f= 0.587 + e0
  • 0.331f将从 .331 转换为浮点格式,结果会有一些舍入误差e1,所以结果是0.331f= .331 + e1
  • valA * 0.587f将计算出一些错误e2,因此结果将是valA * 0.587f= valA0.587f+ e2
  • valB * 0.331f将计算出一些错误e3,因此结果将是valB * 0.331f= valB0.331f+ e3
  • 这两个产品将被添加,有一些错误e4,所以结果将是valA * 0.587f + valB * 0.331f= valA * 0.587f++ valB * 0.331fe4

现在我们可以替换表达式,所以:

  • valA * 0.587f + valB * 0.331f= ( valA0.587f+ e2 ) + ( valB0.331f+ e3 + e4
  • valA * 0.587f + valB * 0.331f= ( valA• (0.587 + e0 ) + e2 ) + ( valB• (.331 + e1 ) + e3 ) + e4
  • valA * 0.587f + valB * 0.331f= (( A + eA ) • (0.587 + e0 ) + e2 ) + (( B + eB ) • (.331 + e1 ) + e3 ) + e4

这样,我们将计算结果 表示valA * 0.587f + valB * 0.331f为精确的数学表达式(具有不完全已知值的变量), (( A + eA ) • (0.587 + e0 ) + e2 ) + (( B + eB ) • (. 331 + e1 ) + e3 ) + e4

数值分析

接下来,我们可以对错误设置一些界限。e0e1很容易,它们的幅度分别为 0.587 和 0.331 的 ½ ULP。.587 在 [½, 1) 中,所以它的 ULP 是 2 -24,而 .331 在 [¼, ½) 中,所以它的 ULP 是 2 -25。所以| e0 | ≤= 2 -25和 | e1 | ≤= 2 -26

e2e3的边界取决于 和 的valA * 0.587f大小valB * 0.331f。因为val< 256, valA * 0.587f< 256,所以它的 ULP 最多为 2 -16,并且 | e2 | ≤ 2 -17。有了valB,我们可以看到valB * 0.331f< 128,因此 ULPvalB * 0.331f最多为 2 -17,并且 | e3 | ≤ 2 -18

最后,我们在最后添加valA * 0.587f + valB * 0.331f. 我们假设它小于 256,所以它的 ULP 最多为 2 -16,并且 | e4 | ≤ 2 -17

查看计算结果的数学表达式 (( A + eA ) • (0.587 + e0 ) + e2 ) + (( B + eB ) • (.331 + e1 ) + e3 ) + e4,我们可以看到当e0e1e2e3e4具有最大值时,会发生最大可能的错误(除非eAeB是巨大的负数,我们假设这不是真的)。因此,我们可以替换为这些错误准备的上限:

(( A + eA ) • (0.587 + 2 -25 ) + 2 -17 ) + (( B + eB ) • (.331 + 2 -26 ) + 2 -18 ) + 2 -17

为了时间的利益,我用 Maple 对此进行了评估。(手动扩展表达式并保留一些因子而不是将系数合并为单个数字可能更有启发性,但我将其留给读者。)结果是:

2462056573/4194304000 • A + 2462056573/4194304000 • eA + 5/262144 + 2776629373/8388608000 • B + 2776629373/8388608000 • eB

理想的结果是A • .587 + B • .331。当我们从上面减去它时,结果是计算误差的界限:

1/33554432 • A + 2462056573/4194304000 • eA + 5/262144 + 1/67108864 * B + 2776629373/8388608000 • eB

由于A < 256 和B < 256,我们可以用 256 代替AB,得到:

1/32768 + 2462056573/4194304000 • eA + 2776629373/8388608000 • eB

反转一点 Maple 的算术,即:

2 -15 + (.587 + 2 -25 ) • eA + (.331 + 2 -26 ) • eB

因此,这是 中错误的上限valA * 0.587f + valB * 0.331fvalA通过有关和之间关系的附加信息,它可能会进一步减少valB。此外,将 .587 和 .331 转换为的错误float是完全已知的,因此应该使用这些错误,而不是我在此答案中用作说明的边界。

还需要建立误差的下限。舍入误差可能是负数,我们必须问 (( A + eA ) • (0.587 + e0 ) + e2 ) + (( B + eB ) • (.331 + e1 ) + e3 )的最低可能值是多少+ e4是。由于我现在没有时间,这留给读者。

附录

e0是 13/1048576000。e1是 1/4194304000。那么误差的上限可以减少到 731/32768000 + 4924113/8388608 • eA + 11106517/33554432 • eB,即:

.731•2 -15 + (.587 + .013•2 -20 ) • eA + (.331 + .001•2 -22 ) • eB

于 2018-03-26T14:40:25.393 回答