浮点型变量不确定性
Float type variable uncertainty
我开发了一个图像处理软件，我需要对其进行数值分析，考虑与其操作相关的错误传播以及由此类变量发生的固有舍入引起的浮点型变量的不确定性.
考虑到 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;
如果我知道float型变量的不确定性，这个误差就可以简单的用这个公式计算出来吧？
简介
此答案对以下错误进行了初步检查：
float result = valA * 0.587f + valB * 0.331f;
在此答案中，floating-point 格式的值和使用 floating-point 格式计算的表达式将使用代码样式表示，如 z 或 x * y。数学变量将使用斜体，并且不会采用代码样式，如 z 或 x • y .
我假设所有算法都是用 IEEE-754 基本 32 位二进制 floating-point 完成的。这种格式通常用于 float 类型，尽管一些编程语言实现混合精度，在计算 float 类型的表达式时可能使用 double 或其他精度。我还假设所有算术都是使用 round-to-nearest 模式完成的，与偶数低位的数字有关。
这种格式在有效数中为 24 位，因此最低精度单位 (ULP) 通常是最高有效位值的 2−23 倍。这是可表示值之间的步长。例如，对于 [1, 2) 中的值，ULP 为 2−23。对于 [128, 256) 中的值，ULP 为 27•2−23 = 2−16. （对于次正规值，有效位数的位数较少。最低 ULP 可以是 2−149。超出最大的有限可表示值，到下一个可表示值的步长是无限的。不过这道题只涉及到中等值的值，可以忽略无穷大。)
使用正确舍入计算任何运算的结果最多与正确答案相差 ½ ULP。也就是说，如果我们计算 z = x + y，例如，计算结果 z 不同于精确的数学结果 z = x + y最多为z的½ ULP。 （虽然z是一个精确的无限精度的数学结果，但是我们用它的大小来确定它落在floating-point格式中的哪个范围，这就决定了我们所说的ULP z。）误差最多为 ½ ULP 的原因是，如果最接近 z 的两个可表示值是 z0 和z1，我们必须有 z0 ≤ z ≤ z1，并且如果 ½ ULP < z1 − z ，然后 z − z0 < ½ ULP（因为 z1 − z0 = 1 ULP，根据 ULP 的定义。 ) 因此，在选择最接近的可表示值时，我们会选择 z0 和 z1 中较接近的一个，因此误差永远不会超过 ½ ULP。
如评论中所述，valA、valB 和 result 在 [0, 256] 中。
符号分析
当我们开始计算 valA * 0.587f + valB * 0.331f 时，valA 和 valB 有一些先前操作的错误。也就是说，理想情况下，使用精确的数学，我们会计算一些数字 A 和 B，但计算机计算的是 valA 和valB，相差eA = valA − A 和eB = valB - B.
理想情况下，我们想计算数字 R 使得 R 使用精确的数学计算 A • .587 + B • .331。当我们使用计算机算术时：
0.587f会从.587转换成floating-point格式，结果会有一些舍入误差e0，所以结果是0.587f = 0.587 + e0.0.331f会从.331转换成floating-point格式，结果会有一些舍入误差e1，所以结果是0.331f = .331 + e1.valA * 0.587f 将计算出一些错误 e2，因此结果将是 valA * 0.587f = valA • 0.587f + e2.valB * 0.331f 将计算出一些错误 e3，因此结果将是 valB * 0.331f = valB • 0.331f + e3.- 两个产品会相加，有一些错误e4，所以结果会是
valA * 0.587f + valB * 0.331f = valA * 0.587f + valB * 0.331f + e4.
现在我们可以替换表达式了，所以：
valA * 0.587f + valB * 0.331f = (valA • 0.587f + e2) + (valB • 0.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) + e.
数值分析
接下来，我们可以对错误进行一些限制。 e0 和 e1 很容易，它们的量级最多分别为 .587 和 .331 的 ½ ULP。 .587 在 [½, 1), 所以它的 ULP 是 2−24, .331 在 [¼, ½), 所以它的 ULP 是 2− 25。所以 |e0| ≤= 2−25，且|e1| ≤= 2−26.
e2 和 e3 的界限取决于 valA * 0.587f 和 valB * 0.331f 的大小。由于val < 256，valA * 0.587f < 256，所以其ULP至多为2−16，且|e2 | ≤ 2−17。有了valB，我们可以看到valB * 0.331f < 128，所以valB * 0.331f的ULP最多为2−17，|e3| ≤ 2−18.
最后，我们在最后添加 valA * 0.587f + valB * 0.331f 时出现错误 e4。我们假设它小于 256，所以它的 ULP 最多为 2−16，并且 |e4| ≤ 2−17.
看计算结果的数学表达式，((A + eA) • (0.587 + e0) + e2) + ((B + eB) • (. 331 + e1) + e3) + e4, 可以看出最大的可能误差发生在e0、e1、e2、e3 和 e4 具有最大值（除非 eA 或 eB 是巨大的负数，我们假设这不是真的）。所以我们可以替换我们为这些错误准备的上限：
((A + eA) • (0.587 + 2−25) + 2−17) + ((B + eB) • (.331 + 2−26) + 2−18) + 2−17.
为了赶时间，我用Maple测评了这个。 （手动扩展表达式并保留一些因子而不是将系数合并为单个数字可能更有启发性，但我将其留给 reader。）结果是：
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 代替 A对于 B，产生：
1/32768 + 2462056573/4194304000 • eA + 2776629373/8388608000 • eB.
逆向一点Maple的算法，即：
2−15 + (.587 + 2−25) • eA + (.331 + 2−26) • eB.
因此，这是 valA * 0.587f + valB * 0.331f 中误差的上限。使用有关 valA 和 valB 之间关系的附加信息，它可能会减少更多。此外，将 .587 和 .331 转换为 float 的错误是完全已知的，因此应该使用这些错误而不是我在此答案中用作说明的界限。
还需要确定误差的下限。舍入误差可能是负数，我们不得不问 ((A + eA) • (0.587 + e0) + e2) + ((B + eB) • (.331 + e1) + e3) + e4 是。由于我现在没时间，所以留给 reader.
附录
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.
我开发了一个图像处理软件，我需要对其进行数值分析，考虑与其操作相关的错误传播以及由此类变量发生的固有舍入引起的浮点型变量的不确定性.
考虑到 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;
如果我知道float型变量的不确定性，这个误差就可以简单的用这个公式计算出来吧？
简介
此答案对以下错误进行了初步检查：
float result = valA * 0.587f + valB * 0.331f;
在此答案中，floating-point 格式的值和使用 floating-point 格式计算的表达式将使用代码样式表示，如 z 或 x * y。数学变量将使用斜体，并且不会采用代码样式，如 z 或 x • y .
我假设所有算法都是用 IEEE-754 基本 32 位二进制 floating-point 完成的。这种格式通常用于 float 类型，尽管一些编程语言实现混合精度，在计算 float 类型的表达式时可能使用 double 或其他精度。我还假设所有算术都是使用 round-to-nearest 模式完成的，与偶数低位的数字有关。
这种格式在有效数中为 24 位，因此最低精度单位 (ULP) 通常是最高有效位值的 2−23 倍。这是可表示值之间的步长。例如，对于 [1, 2) 中的值，ULP 为 2−23。对于 [128, 256) 中的值，ULP 为 27•2−23 = 2−16. （对于次正规值，有效位数的位数较少。最低 ULP 可以是 2−149。超出最大的有限可表示值，到下一个可表示值的步长是无限的。不过这道题只涉及到中等值的值，可以忽略无穷大。)
使用正确舍入计算任何运算的结果最多与正确答案相差 ½ ULP。也就是说，如果我们计算 z = x + y，例如，计算结果 z 不同于精确的数学结果 z = x + y最多为z的½ ULP。 （虽然z是一个精确的无限精度的数学结果，但是我们用它的大小来确定它落在floating-point格式中的哪个范围，这就决定了我们所说的ULP z。）误差最多为 ½ ULP 的原因是，如果最接近 z 的两个可表示值是 z0 和z1，我们必须有 z0 ≤ z ≤ z1，并且如果 ½ ULP < z1 − z ，然后 z − z0 < ½ ULP（因为 z1 − z0 = 1 ULP，根据 ULP 的定义。 ) 因此，在选择最接近的可表示值时，我们会选择 z0 和 z1 中较接近的一个，因此误差永远不会超过 ½ ULP。
如评论中所述，valA、valB 和 result 在 [0, 256] 中。
符号分析
当我们开始计算 valA * 0.587f + valB * 0.331f 时，valA 和 valB 有一些先前操作的错误。也就是说，理想情况下，使用精确的数学，我们会计算一些数字 A 和 B，但计算机计算的是 valA 和valB，相差eA = valA − A 和eB = valB - B.
理想情况下，我们想计算数字 R 使得 R 使用精确的数学计算 A • .587 + B • .331。当我们使用计算机算术时：
0.587f会从.587转换成floating-point格式，结果会有一些舍入误差e0，所以结果是0.587f= 0.587 + e0.0.331f会从.331转换成floating-point格式，结果会有一些舍入误差e1，所以结果是0.331f= .331 + e1.valA * 0.587f将计算出一些错误 e2，因此结果将是valA * 0.587f=valA•0.587f+ e2.valB * 0.331f将计算出一些错误 e3，因此结果将是valB * 0.331f=valB•0.331f+ e3.- 两个产品会相加，有一些错误e4，所以结果会是
valA * 0.587f + valB * 0.331f=valA * 0.587f+valB * 0.331f+ e4.
现在我们可以替换表达式了，所以：
valA * 0.587f + valB * 0.331f= (valA•0.587f+ e2) + (valB•0.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) + e.
数值分析
接下来，我们可以对错误进行一些限制。 e0 和 e1 很容易，它们的量级最多分别为 .587 和 .331 的 ½ ULP。 .587 在 [½, 1), 所以它的 ULP 是 2−24, .331 在 [¼, ½), 所以它的 ULP 是 2− 25。所以 |e0| ≤= 2−25，且|e1| ≤= 2−26.
e2 和 e3 的界限取决于 valA * 0.587f 和 valB * 0.331f 的大小。由于val < 256，valA * 0.587f < 256，所以其ULP至多为2−16，且|e2 | ≤ 2−17。有了valB，我们可以看到valB * 0.331f < 128，所以valB * 0.331f的ULP最多为2−17，|e3| ≤ 2−18.
最后，我们在最后添加 valA * 0.587f + valB * 0.331f 时出现错误 e4。我们假设它小于 256，所以它的 ULP 最多为 2−16，并且 |e4| ≤ 2−17.
看计算结果的数学表达式，((A + eA) • (0.587 + e0) + e2) + ((B + eB) • (. 331 + e1) + e3) + e4, 可以看出最大的可能误差发生在e0、e1、e2、e3 和 e4 具有最大值（除非 eA 或 eB 是巨大的负数，我们假设这不是真的）。所以我们可以替换我们为这些错误准备的上限：
((A + eA) • (0.587 + 2−25) + 2−17) + ((B + eB) • (.331 + 2−26) + 2−18) + 2−17.
为了赶时间，我用Maple测评了这个。 （手动扩展表达式并保留一些因子而不是将系数合并为单个数字可能更有启发性，但我将其留给 reader。）结果是：
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 代替 A对于 B，产生：
1/32768 + 2462056573/4194304000 • eA + 2776629373/8388608000 • eB.
逆向一点Maple的算法，即：
2−15 + (.587 + 2−25) • eA + (.331 + 2−26) • eB.
因此，这是 valA * 0.587f + valB * 0.331f 中误差的上限。使用有关 valA 和 valB 之间关系的附加信息，它可能会减少更多。此外，将 .587 和 .331 转换为 float 的错误是完全已知的，因此应该使用这些错误而不是我在此答案中用作说明的界限。
还需要确定误差的下限。舍入误差可能是负数，我们不得不问 ((A + eA) • (0.587 + e0) + e2) + ((B + eB) • (.331 + e1) + e3) + e4 是。由于我现在没时间，所以留给 reader.
附录
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.