浮点精度和操作顺序

如何解决浮点精度和操作顺序

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

我要做的实际上是生成两个伪随机向量b和c，以及伪随机标量s，然后检查这些向量上不同运算的结果。 / p>

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的乘积会导致如此大的差异吗？天真的想法是说，鉴于我的机器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

确实，在减法中，我失去了更多的信息。那是对的吗？怎么用浮点近似来解释？

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

解决方法

信息的最大损失很可能发生在点积而不是叉积中。在叉积中，您得到的结果仍然接近SELECT ... FROM bot b LEFT JOIN switch s ON s.botid = b.id AND s.date = ( SELECT MAX(s1.date) FROM switch s1 WHERE s1.User = 'SomeUsername' AND s1.botid = s.botid ) WHERE b.Active = 1中条目的数量级。这意味着您可能损失了大约一位数的精度，但是相对误差仍应在10 ^ -15左右。（减法c中的相对误差大约等于a-b）

点积是唯一涉及将两个非常接近的数字相减的运算。由于我们将先前的相对误差除以〜0，因此导致相对误差的极大增加。

现在来看您的示例，考虑到您拥有的数量：2*(|a|+|b|) / (a-b)，c和{{1} }的大小约为10 ^ 5，这意味着绝对误差最多为10 ^ -11。我不在乎e，因为它基本上等于1。

乘以d时的绝对误差大约为s（最坏的情况是误差不会抵消）。现在，在点积中，您将2个数量级乘以〜10 ^ 5，因此误差应该在a*b的范围内（并乘以3，因为您对每个分量执行了3次）。>

然后，如果10 ^ -6是预期的错误，我该如何解释您的结果？好吧，您很幸运：使用这些值（我更改了|a|*|err_b| + |b|*|err_a|和10^5*10^-11 + 10^5*10^-11 = 2*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

关于运算顺序，只要您不减去2个非常接近的数字，我认为这没什么大不了的。如果您仍然需要减去2个非常接近的数字，那么我想最好还是这样做（不要把所有东西搞砸了），但不要在那儿引用我。

Miguel's answer发现了。就像附录一样，由于OP与C ++一起使用，我以我所知道的最准确的方式对计算进行了编码，并尽可能地利用了融合的乘法加法运算。另外，我尝试了补偿点产品。可以认为这是因为Kahan sum的概念扩展到了点积的累积。在这里没有明显的区别。

下面的代码输出，如果使用最严格的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,a.x,b.z);
    r.z = diff_of_prod (a.x,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,-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.56005552104818945};
    double3 c = {-178151.26386435505,-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,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,compensated_dot (d,compensated_dot (e,s * compensated_dot (e,c));

    return EXIT_SUCCESS;
}

浮点精度和操作顺序

如何解决浮点精度和操作顺序

解决方法

相关推荐