logo标签列表

Matlab/Octave/Numpy数值差异

matlabnumpyoctave
6

我正在将一个在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为什么会这样吗？

3个回答
6
他们在内部使用双精度浮点数，其精度仅约为15位数字。您的数学运算很可能超过了这个范围，导致错误的数学结果。
值得一读：http://floating-point-gui.de/ 编辑：从文档中我发现Numpy没有更高的精度可用。看起来SymPy虽然可能会给你所需的精度——如果该库也适用于您。
那不太对。虽然有float128数据类型，但它的精度可能并不总是被很好地定义。 - seberg
@Sebastian，我找不到任何关于float128类型的参考资料 - 只有complex128（因为这是两个float64视为一个具有实部和虚部的数字）。http://docs.scipy.org/doc/numpy/user/basics.types.html - Lucero
是的...那是因为float128只在运行它的计算机上可用。但在通常的个人电脑上是可以使用的。 - seberg
4

不管价值多少，我在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），则会遇到精度问题，并得到不同的答案，如下所述。您可以尝试在这种情况下明确指定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)
1
我并没有真正看到32位与64位的区别（Matlab + Numpy默认使用双精度）。然而，Blas库可能是关键点，当我使用ATLAS时，我得到了与@user1655812相同的结果，否则我得到了完全不同的结果。另外，np.einsum可以产生相同的结果，同时避免使用ATLAS。 - seberg
32位和64位之间的区别在这种情况下足以显现。无论如何，我得到了一个明显不同的答案，但这并不能完全解释OP的结果。我同意，这可能是由于BLAS库，但我没有想到去测试它。（很高兴你这样做！）我的结果没有使用ATLAS。（关于einsum的好观点！） - Joe Kington
1
抱歉，是的，但它始终是64位浮点数，系统并不重要。 - seberg
我原本以为在32位系统上，numpy默认使用32位精度数组。但是我可能完全错了，这就是我所依据的。编辑：撤销刚才的话，我记错了。 - Joe Kington
0

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]]

第二行和第二列的非对角线项不同，但其他位置都是0。

使用双点符号，P.dot(A.T)在添加乘积时会产生舍入误差。这些误差会传播到下一个dot中。einsum生成所有乘积，并进行一次求和。我怀疑MATLAB解释器会识别这种情况，并执行特殊计算以最小化舍入误差。

Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix - 是一个最近的问题，其解释相同。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接