numpy中叉积的数值精度

Question

我遇到了我认为在 numpy 中使用叉积且值有点大的奇怪行为。

例如，以下似乎是正确的：

r = 1e15
a = array([1, 2, 3]) * r;
b = array([-1, 2, 1]) * r;

c = cross(a / norm(a), b / norm(b));

print(dot(c, a))  # outputs 0.0

但是如果我们将指数增加 1，我们得到：

r = 1e16
a = array([1, 2, 3]) * r;
b = array([-1, 2, 1]) * r;

c = cross(a / norm(a), b / norm(b));

print(dot(c, a))  # outputs 2.0

对于更大的指数值，数字会变得更加奇怪。有谁知道这里发生了什么？谢谢！

Answer 1

您看到舍入错误。默认情况下，array() returns 具有 dtype=float64 的对象。随着 r 越来越大，您运行的尾数空间无法准确表示数组乘积。这是一种测试方法：

def testcross(r, dt):
    a = array([1, 2, 3], dtype=dt)*r
    b = array([-1, 2, 1], dtype=dt)*r
    c = cross(a/norm(a), b/norm(b))
    return dot(c, a)

for rr in logspace(4, 15, 10):
    print "%10.2f %10.2f %g" % (testcross(rr, float32), testcross(rr, float64)

结果：

     -0.00       0.00 10000
      0.00      -0.00 166810
      0.00       0.00 2.78256e+06
     -4.00       0.00 4.64159e+07
    -64.00       0.00 7.74264e+08
   1024.00       0.00 1.29155e+10
      0.00       0.00 2.15443e+11
-524288.00       0.00 3.59381e+12
      0.00      -0.02 5.99484e+13
-134217728.00       0.00 1e+15

请注意，即使 float64 和 r=5.99484e13 也不是 "perfect"。这表明在您达到 r=1e15 之前很久，精度就开始下降，即使是 float64。正如预期的那样，使用不太精确的 float32.

情况会更糟

跟进OP的建议：32位和64位浮点表示的尾数字段分别为24位和53位（包括隐含位）。取 log10([2**24, 2**53])，我们看到这分别对应大约 7 和 16 个数量级。这与最初指出的 float32 和 r=1e16 的 r=4.6e7 附近出现的表格错误相对应。当点积导致基础矩阵计算减去大数时，会发生舍入，并且差异不能表示一个或另一个大数的尾数。

numpy中叉积的数值精度

numerical precision of cross product in numpy

python

precision

numpy

cross-product