scipy vs matlab,Python NumPy vs Octave/MATLAB精度

这个问题是关于使用NumPy与Octave/MATLAB的计算精度(但是下面的MATLAB代码只在Octave下进行了测试)。我知道Stackoverflow上有一个类似的问题，即this，但这似乎与我下面要问的有点远。在

设置

一切都在ubuntu14.04上运行。在

Python 3.4.0版。在

根据OpenBLAS编译的NumPy版本1.8.1。在

Octave版本3.8.1是根据OpenBLAS编译的。在

示例代码

示例Python代码。在import numpy as np

from scipy import linalg as la

def build_laplacian(n):

lap=np.zeros([n,n])

for j in range(n-1):

lap[j+1][j]=1

lap[j][j+1]=1

lap[n-1][n-2]=1

lap[n-2][n-1]=1

return lap

def evolve(s, lap):

wave=la.expm(-1j*s*lap).dot([1]+[0]*(lap.shape[0]-1))

for i in range(len(wave)):

wave[i]=np.linalg.norm(wave[i])**2

return wave

我们现在运行以下内容。在

^{pr2}$

它给出了e-34的顺序。在

我们可以在Octave/MATLAB中生成类似的代码：function lap=build_laplacian(n)

lap=zeros(n,n);

for i=1:(n-1)

lap(i+1,i)=1;

lap(i,i+1)=1;

end

lap(n,n-1)=1;

lap(n-1,n)=1;

end

function result=evolve(s, lap)

d=zeros(length(lap(:,1)),1); d(1)=1;

result=expm(-1i*s*lap)*d;

for i=1:length(result)

result(i)=norm(result(i))^2;

end

end

然后我们跑min(evolve(2, build_laplacian(500)))

得到0。事实上，evolve(2, build_laplacian(500)))(60)给出了大约e-100或更小的值(如预期)。在

问题

有人知道是什么导致了NumPy和Octave之间如此大的差异吗(同样，我还没有用MATLAB测试过代码，但我希望看到类似的结果)。在

当然，我们也可以通过首先对角化矩阵来计算矩阵的指数。我已经做了这个，并得到了类似或更差的结果(与纽比)。在

编辑

我的scipy版本是0.14.0。我知道Octave/MATLAB使用的是Pade近似格式，并且对这个算法很熟悉。我不确定scipy是做什么的，但是我们可以试试下面的方法。在用numpy的eig或{}对角化矩阵(在我们的例子中，后者工作得很好，因为矩阵是Hermitian)。得到了两个矩阵：对角矩阵D，矩阵U，其中{}由原矩阵在对角线上的特征值组成，U由相应的特征向量列组成，因此原始矩阵由U.T.dot(D).dot(U)给出。

指数D(这现在很容易，因为D是对角线的)。

现在，如果M是原始矩阵，d是原始向量{}，我们得到{}。

不幸的是，这会产生和以前一样的结果。因此，这可能与numpy.linalg.eig和{}的工作方式有关，或者与{}内部进行算术的方式有关。在

所以问题是：我们如何提高numpy的精度？实际上，如上所述，Octave在这种情况下似乎做得更好。在