浮点精度和运算顺序

Question

我正在为 3D 矢量对象及其代数（点积、叉积等）的类编写单元测试，并且刚刚观察到我可以理解的行为，但还没有完全理解。

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

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

我在 python 中创建了一个 MWE（使用 numpy）只是为了更好地展示我的问题（但我实际上是在 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，您仍然需要将相当大的分量乘以一个数字
• 当取点积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

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

另外，这是否意味着，无论输入如何（即，无论两个向量的大小相似还是完全不同），最好总是首先执行所有涉及乘法（和除法？）的操作，那么那些涉及加法/减法的？

Answer 1

Miguel's answer 是正确的。就像附录一样，由于 OP 与 C++ 一起工作，我以我所知道的最准确的方式对计算进行了编码，尽可能地利用了融合乘加运算。此外，我尝试了补偿点积。可以将其视为将 Kahan 和扩展到点积的累积的想法。在这里没有显着差异。

当使用最严格的 IEEE-754 合规性编译器（对于我的 Intel 编译器，即 /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;
}

Answer 2

信息的大量丢失最有可能发生在点积中，而不是叉积中。在叉积中，你得到的结果仍然接近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 个非常接近的数字，我想最好最后还是这样做（不要把所有事情都搞砸），但不要引用我的话。