Matlab/Octave/Numpy 数值差答案

【问题标题】：Matlab/Octave/Numpy numeric differenceMatlab/Octave/Numpy 数值差
【发布时间】：2012-09-01 19:20:58
【问题描述】：

我正在将一个在 Matlab 中工作的算法移植到 numpy 并观察到一个奇怪的行为。相关的代码段是

P = eye(4)*1e20;
A = [1 -0.015 -0.025 -0.035; 0.015 1 0.035 -0.025; 0.025 -0.035 1 0.015; 0.035 0.025 -0.015 1];
V1 = A*(P*A')
V2 = (A*P)*A'

当我使用 Matlab 运行时，这段代码提供了以下矩阵：

V1 = 1.0021e+20            0  -8.0000e+00            0
              0   1.0021e+20            0            0
    -8.0000e+00            0   1.0021e+20            0
              0            0            0   1.0021e+20


V2 = 1.0021e+20            0  -8.0000e+00            0
              0   1.0021e+20            0            0
    -8.0000e+00            0   1.0021e+20            0
              0            0            0   1.0021e+20

请注意，正如预期的那样，V1 和 V2 是相同的。

当相同的代码在 Octave 中运行时，它提供：

V1 = 1.0021e+20   4.6172e+01  -1.3800e+02   1.8250e+02
    -4.6172e+01   1.0021e+20  -1.8258e+02  -1.3800e+02
     1.3801e+02   1.8239e+02   1.0021e+20  -4.6125e+01
    -1.8250e+02   1.3800e+02   4.6125e+01   1.0021e+20

V2 = 1.0021e+20  -4.6172e+01   1.3801e+02  -1.8250e+02
     4.6172e+01   1.0021e+20   1.8239e+02   1.3800e+02
    -1.3800e+02  -1.8258e+02   1.0021e+20   4.6125e+01
     1.8250e+02  -1.3800e+02  -4.6125e+01   1.0021e+20

在numpy中，段变为

from numpy import array, dot, eye
A = numpy.array([[1, -0.015, -0.025, -0.035],[0.015, 1, 0.035, -0.025],[0.025, -0.035, 1, 0.015],[0.035, 0.025, -0.015, 1]])
P = numpy.eye(4)*1e20
print numpy.dot(A,numpy.dot(P,A.transpose()))
print numpy.dot(numpy.dot(A,P),A.transpose())

哪个输出

[[  1.00207500e+20   4.61718750e+01  -1.37996094e+02   1.82500000e+02]
 [ -4.61718750e+01   1.00207500e+20  -1.82582031e+02  -1.38000000e+02]
 [  1.38011719e+02   1.82386719e+02   1.00207500e+20  -4.61250000e+01]
 [ -1.82500000e+02   1.38000000e+02   4.61250000e+01   1.00207500e+20]]
[[  1.00207500e+20  -4.61718750e+01   1.38011719e+02  -1.82500000e+02]
 [  4.61718750e+01   1.00207500e+20   1.82386719e+02   1.38000000e+02]
 [ -1.37996094e+02  -1.82582031e+02   1.00207500e+20   4.61250000e+01]
 [  1.82500000e+02  -1.38000000e+02  -4.61250000e+01   1.00207500e+20]]

因此，Octave 和 numpy 都提供了相同的答案，但与 Matlab 的答案大不相同。第一点是 V1 != V2，这似乎不对。另一点是，虽然非对角元素比对角元素小很多数量级，但这似乎在我的算法中造成了一些问题。

有人知道 numpy 和 Octave 的行为方式吗？

【问题讨论】：

标签： matlab numpy octave

【解决方案1】：

他们在内部使用双精度，只有大约 15 位精度。您的数学运算可能会超出此范围，从而导致数学上错误的结果。

值得一读：http://floating-point-gui.de/

编辑：我从docs 得知 Numpy 没有更高的精度。似乎 SymPy 虽然 may give you the needed precision - 如果该库也适合您。

【讨论】：

这不太对。有一个 float128 数据类型，但我认为它的精度可能并不总是很好定义。
@Sebastian，我没有发现对 float128 类型的任何引用 - 只有 complex128 （因为它们是两个 float64 被视为具有实部和虚部的一个数字）。 docs.scipy.org/doc/numpy/user/basics.types.html
是的......那是因为 float128 仅根据其运行的计算机可用。但在通常的 PC 上是这样。

【解决方案2】：

无论如何，我在 64 位系统上得到的结果与 matlab 相同：

[[  1.00207500e+20   0.00000000e+00  -8.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   1.00207500e+20   0.00000000e+00   0.00000000e+00]
 [ -8.00000000e+00   0.00000000e+00   1.00207500e+20   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.00207500e+20]]
[[  1.00207500e+20   0.00000000e+00  -8.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   1.00207500e+20   0.00000000e+00   0.00000000e+00]
 [ -8.00000000e+00   0.00000000e+00   1.00207500e+20   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.00207500e+20]]
[[  1.00207500e+20   0.00000000e+00  -8.00000000e+00   0.00000000e+00]

如果您使用的是 32 位系统（或者您在 64 位系统上安装了 32 位版本的 python 和 numpy），您将遇到精度问题，并得到不同的答案，如前所述由@Lucero 下面。在这种情况下，您可以尝试显式指定 64 位浮点数（但操作会更慢）。例如，尝试使用np.array(..., dtype=np.float64)。

如果您认为需要更高的精度，可以使用np.longdouble（在 64 位系统上与 np.float128 相同，在 32 位系统上与 np.float96 相同），但这可能并非在所有平台上都支持线性代数函数会将事物截断回原始精度。

另外，您使用的是什么 BLAS 库？ numpy 和 octave 的结果可能相同，因为它们使用的是相同的 BLAS 库。

最后，您可以将 numpy 代码简化为：

import numpy as np
A = np.array([[1,     -0.015, -0.025, -0.035],
              [0.015, 1,      0.035,  -0.025],
              [0.025, -0.035, 1,      0.015],
              [0.035, 0.025,  -0.015, 1]])
P = np.eye(4)*1e20
print A.dot(P.dot(A.T))
print A.dot(P).dot(A.T)

【讨论】：

我并没有真正看到 32 位和 64 位的点（matlab+numpy 都默认使用双精度）。然而，Blas 库可能是重点，当使用 ATLAS 时，我得到了与 @user1655812 相同的结果，但没有得到确切的结果。加入其他东西 np.einsum 可以做同样的结果，同时在这里避免 ATLAS。
32 位和 64 位之间的差异足以在这种情况下显示出来。无论如何，我得到了一个截然不同的答案，但这还不足以完全解释 OP 的结果。我同意，这可能是由于 BLAS 库，但我没想过要测试它。（很高兴你做到了！）我的结果没有使用 ATLAS。（关于 einsum 的好处！）
对不起，是的，但它始终是 64 位浮点数，系统无关紧要。
我的印象是 numpy 在 32 位系统上默认为 32 位精度数组，但是？（我可能完全错了，但这就是我要说的。）编辑：从头开始，我记错了。

【解决方案3】：

np.einsum 更接近 MATLAB

In [1299]: print(np.einsum('ij,jk,lk',A,P,A))
[[  1.00207500e+20   0.00000000e+00  -5.07812500e-02   0.00000000e+00]
 [  0.00000000e+00   1.00207500e+20   5.46875000e-02   0.00000000e+00]
 [ -5.46875000e-02   5.46875000e-02   1.00207500e+20   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   1.00207500e+20]]

第 2 行和第 2 列的非对角项不同，但其他地方的 0 相同。

使用双点，P.dot(A.T) 在添加产品时会产生舍入错误。这些传播到下一个dot。 einsum 生成所有产品，只有一个总和。我怀疑 MATLAB 解释器会识别这种情况，并执行旨在最小化舍入误差的特殊计算。

Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix - 是一个更新的问题，具有相同的解释。

【讨论】：