浮点精度和运算顺序

Question

我正在为 3D 矢量对象及其代数（点积、叉积等）编写 class 的单元测试，并且刚刚观察到一种我能以某种方式理解但不完全理解的行为范围。

我所做的实际上是生成 2 个伪随机向量，b 和 c，以及一个伪随机标量，s，然后检查对这些向量进行不同操作的结果.

b 的组件在 [-1, 1] 范围内生成，而 c 的组件在 [-1e6, 1e6] 范围内生成，因为在我的用例中我将遇到类似的情况，可能会导致尾数中的信息大量丢失。 s 也在 [-1, 1] 范围内生成。

我在 python（使用 numpy）中创建了一个 MWE，只是为了更好地揭示我的问题（但我实际上是在用 C++ 编码，问题本身与语言无关）：

b = np.array([0.4383006177615909, -0.017762134447941058, 0.56005552104818945])
c = np.array([-178151.26386435505, 159388.59511391702, -720098.47337336652])
s = -0.19796489160874975

我再定义

d = s*np.cross(b,c)
e = np.cross(b,c)

最后计算

In [7]: np.dot(d,c)
Out[7]: -1.9073486328125e-06

In [8]: np.dot(e,c)
Out[8]: 0.0

In [9]: s*np.dot(e,c)
Out[9]: -0.0

因为d和e都垂直于b和c，上面计算的标量积应该都是0（代数）。

现在，我很清楚，在真正的计算机中，这只能在浮点运算的限制下实现。然而，我想更好地理解这个错误是如何产生的。

让我有点吃惊的是三个结果中第一个的准确性很差。

我将尝试在以下方面表达我的想法：

np.cross(b, c) 基本上是 [b[1]*c[2]-b[2]*c[1], b[2]*c[0]-b[0]*c[2], ...]，它涉及大数和小数的乘法以及随后的减法。 e（叉积 b x c）本身保留了相对较大的分量，即 array([-76475.97678585, 215845.00681978, 66695.77300175])
因此，要获得 d，您仍然需要将相当大的分量乘以 <1 的数字。这当然会导致一些截断错误。
点积e . c的结果是正确的，而d . c的结果差了差不多2e-6。最后乘以 s 会导致这么大的差异吗？一个天真的想法是，考虑到我的机器 epsilon 2.22045e-16 和分量 d 的大小，误差应该在 4e-11.
叉积减法中尾数的信息丢失了吗？

为了检查最后的想法，我做了：

In [10]: d = np.cross(s*b,c)                                                    

In [11]: np.dot(d,c)                                                            
Out[11]: 0.0

In [12]: d = np.cross(b,s*c)                                                    

In [13]: np.dot(d,c)                                                            
Out[13]: 0.0

而且看起来确实在减法中我丢失了更多信息。那是对的吗？如何用浮点数近似来解释？

此外，这是否意味着，无论输入如何（即，无论两个向量大小相似还是完全不同），最好始终先执行所有涉及乘法（和除法）的运算？），然后涉及 addition/subtraction?

Answer 1

大量信息丢失很可能发生在点积中，而不是叉积中。在叉积中，您得到的结果仍然接近 c 中条目的数量级。这意味着您可能损失了大约一位数的精度，但相对误差仍应在 10^-15 左右。（减法的相对误差a-b约等于2*(|a|+|b|) / (a-b)）

点积是唯一涉及两个彼此非常接近的数字相减的运算。这会导致相对误差的大幅增加，因为我们将之前的相对误差除以~0。

现在来看你的例子，你得到的错误 (~10^-6) 实际上是你所期望的，考虑到你拥有的数量：c、e 和 d 的大小约为 10^5，这意味着绝对误差最多在 10^-11 左右。我不关心 s 因为它基本上等于 1.

乘以 a*b 时的绝对误差约为 |a|*|err_b| + |b|*|err_a|（误差未抵消的最坏情况）。现在在点积中，你将 2 个数量乘以 ~10^5，所以误差应该在 10^5*10^-11 + 10^5*10^-11 = 2*10^-6 的范围内（并乘以 3，因为你对每个分量执行了 3 次）。

那么如果 10^-6 是预期的误差，我该如何解释你的结果？嗯，你很幸运：使用这些值（我更改了 b[0] 和 c[0]）

b = np.array([0.4231830061776159, -0.017762134447941058, 0.56005552104818945])
c = np.array([-178151.28386435505, 159388.59511391702, -720098.47337336652])
s = -0.19796489160874975

我得到了（按顺序）

-1.9073486328125e-06
7.62939453125e-06
-1.5103522614192943e-06

-1.9073486328125e-06
-1.9073486328125e-06

此外，当您查看相对误差时，它做得很好：

In [10]: np.dot(d,c)
Out[11]: -1.9073486328125e-06

In [11]: np.dot(d,c) / (np.linalg.norm(e)*np.linalg.norm(c))
Out[11]: -1.1025045691772927e-17

关于运算顺序，我认为这并不重要，只要您不减去 2 个非常接近的数字即可。如果您仍然需要减去 2 个非常接近的数字，我想最好在最后做（不要搞砸一切）但不要引用我的话。

Answer 2

是正确的。正如附录一样，由于 OP 使用 C++，我以我所知道的最准确的方式对计算进行编码，尽可能利用融合 multiply-add 操作。此外，我尝试了补偿点积。人们可以认为这是将 Kahan 和的概念扩展到点积的累加。在这里没有显着差异。

我下面的代码的输出，当使用最严格的 IEEE-754 合规编译器编译时（对于我的英特尔编译器，即 /fp:strict），应该看起来类似于：

Using FMA-based dot product:
dot(d,c)   = -1.0326118360251935e-006
dot(e,c)   =  4.3370577648224470e-006
s*dot(e,c) = -8.5858517031396220e-007
Using FMA-based compensated dot product:
dot(d,c)   = -1.1393800219802703e-006
dot(e,c)   =  3.0970281801622503e-006
s*dot(e,c) = -6.1310284799506335e-007

#include <cstdio>
#include <cstdlib>
#include <cmath>

typedef struct {
    double x;
    double y;
} double2;

typedef struct {
    double x;
    double y;
    double z;
} double3;

/*
  diff_of_prod() computes a*b-c*d with a maximum error < 1.5 ulp

  Claude-Pierre Jeannerod, Nicolas Louvet, and Jean-Michel Muller, 
  "Further Analysis of Kahan's Algorithm for the Accurate Computation 
  of 2x2 Determinants". Mathematics of Computation, Vol. 82, No. 284, 
  Oct. 2013, pp. 2245-2264
*/
double diff_of_prod (double a, double b, double c, double d)
{
    double w = d * c;
    double e = fma (-d, c, w);
    double f = fma (a, b, -w);
    return f + e;
}

double3 scale (double3 a, double s)
{
    double3 r;
    r.x = s * a.x;
    r.y = s * a.y;
    r.z = s * a.z;
    return r;
} 

double dot (double3 a, double3 b)
{
    return fma (a.x, b.x, fma (a.y, b.y, a.z * b.z));
}

double3 cross (double3 a, double3 b)
{
    double3 r;
    r.x = diff_of_prod (a.y, b.z, a.z, b.y);
    r.y = diff_of_prod (a.z, b.x, a.x, b.z);
    r.z = diff_of_prod (a.x, b.y, a.y, b.x);
    return r;
}

/* returns the sum of a and b as a double-double */
double2 TwoProdFMA (double a, double b)
{
    double2 r;
    r.x = a * b;
    r.y = fma (a, b, -r.x);
    return r;
}

/* returns the product of a and b as a double-double. Knuth TAOCP */
double2 TwoSum (double a, double b)
{
    double2 res;
    double s, r, t;
    s = a + b;
    t = s - a;
    r = (a - (s - t)) + (b - t);
    res.x = s;
    res.y = r;
    return res;
}

/*
  S. Graillat, Ph. Langlois and N. Louvet, "Accurate dot products with FMA",
  In: RNC-7, Real Numbers and Computer Conference, Nancy, France, July 2006,
  pp. 141-142
*/
double compensated_dot (double3 x, double3 y)
{
    double2 t1, t2, t3;
    double sb, cb, pb, pi, sg;

    t1 = TwoProdFMA (x.x, y.x);
    sb = t1.x;
    cb = t1.y;

    t2 = TwoProdFMA (x.y, y.y);
    pb = t2.x;
    pi = t2.y;
    t3 = TwoSum (sb, pb);
    sb = t3.x;
    sg = t3.y;
    cb = (pi + sg) + cb;

    t2 = TwoProdFMA (x.z, y.z);
    pb = t2.x;
    pi = t2.y;
    t3 = TwoSum (sb, pb);
    sb = t3.x;
    sg = t3.y;
    cb = (pi + sg) + cb;

    return sb + cb;
}

int main (void)
{
    double3 b = {0.4383006177615909, -0.017762134447941058, 0.56005552104818945};
    double3 c = {-178151.26386435505, 159388.59511391702, -720098.47337336652};
    double s = -0.19796489160874975;
    double3 d = scale (cross (b, c), s);
    double3 e = cross (b, c);

    printf ("Using FMA-based dot product:\n");
    printf ("dot(d,c)   = % 23.16e\n", dot (d, c));
    printf ("dot(e,c)   = % 23.16e\n", dot (e, c));
    printf ("s*dot(e,c) = % 23.16e\n", s * dot (e, c));

    printf ("Using FMA-based compensated dot product:\n");
    printf ("dot(d,c)   = % 23.16e\n", compensated_dot (d, c));
    printf ("dot(e,c)   = % 23.16e\n", compensated_dot (e, c));
    printf ("s*dot(e,c) = % 23.16e\n", s * compensated_dot (e, c));

    return EXIT_SUCCESS;
}

浮点精度和运算顺序

Floating point accuracy and order of operations

floating-point

language-agnostic

precision

linear-algebra

floating-accuracy