Python NumPy与Octave/MATLAB精度比较

Question

Python NumPy与Octave/MATLAB精度比较

pythonmatlabnumpyscipyoctave

6

这个问题是关于使用NumPy与Octave/MATLAB进行计算精度的比较（下面的MATLAB代码仅在Octave上进行过测试）。我知道Stackoverflow上有一个类似的问题，即这个，但那似乎与我下面要问的有些远离。

设置

所有东西都在Ubuntu 14.04上运行。

Python版本为3.4.0。

NumPy版本为1.8.1，编译时使用了OpenBLAS。

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

我们现在运行以下内容。

np.min(evolve(2, build_laplacian(500)))

这将会得到大约为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测试过代码，但我希望看到类似的结果）。

当然，也可以通过首先将矩阵对角化来计算矩阵指数。我已经这样做了，并得到了类似或更糟糕的结果（使用NumPy）。

编辑

我的scipy版本是0.14.0。我知道Octave / MATLAB使用Pade逼近方案，并熟悉此算法。我不确定scipy做了什么，但我们可以尝试以下操作。

使用numpy的eig或eigh（在我们的情况下，后者很好用，因为矩阵是Hermitian）对矩阵进行对角化。结果我们得到两个矩阵：一个对角线矩阵D，和矩阵U，其中D由原始矩阵的特征值组成，而U由相应的特征向量组成为列；因此原始矩阵由U.T.dot(D).dot(U)给出。
对D进行指数化（现在易于做到，因为D是对角线）。
现在，如果M是原始矩阵，d是原始向量d=[1]+[0]*n，我们得到scipy.linalg.expm(-1j*s*M).dot(d)=U.T.dot(numpy.exp(-1j*s*D).dot(U.dot(d))。

不幸的是，这产生了与以前相同的结果。因此，这可能与numpy.linalg.eig和numpy.linalg.eigh的工作方式有关，或者与numpy内部执行算术的方式有关。

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

- user2321808

你正在使用scipy中的expm方法。你使用的是哪个版本的scipy? - Warren Weckesser

1

expm是一个棘手的函数。http://www.mathworks.com/help/matlab/ref/expm.html引用了Moler关于“19个可疑方法”的文章。Pade逼近法强烈依赖于预处理策略。Octave文档中的`expm`也讨论了这个逼近问题。 - hpaulj

@hpaulj：是的，我知道这个问题，并且也熟悉Pade近似算法。我的观点是Octave和MATLAB产生的结果非常好，而scipy/numpy则相差几个数量级。我们可以使用numpy的eig或eigh先对矩阵进行对角化，然后取对角矩阵的指数。这会产生相同的结果。我不知道scipy使用哪种算法，但似乎这要么与numpy对角化矩阵的方式有关，要么与内部算术工作有关。问题是，我们如何增加numpy的精度？ - user2321808

@William：expm在scipy中已经实现。开发版本的源代码可以在https://github.com/scipy/scipy/blob/master/scipy/linalg/matfuncs.py找到；那里的`expm`函数只是一个薄包装，真正的实现在https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/matfuncs.py中。 - Warren Weckesser

请在 https://github.com/scipy/scipy/issues 上尝试问题搜索。 - hpaulj

显示剩余2条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- pv. · Accepted Answer

以下是代码：

import numpy as np
from scipy import linalg as la
import scipy

print np.__version__
print scipy.__version__
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]=la.norm(wave[i])**2
  return wave

r = evolve(2, build_laplacian(500))
print np.min(abs(r))
print r[59]

输出结果为：

1.8.1
0.14.0
0
(2.77560227344e-101+0j)

使用OpenBLAS 0.2.8-6ubuntu1，对我而言，似乎没有立即重现您的问题。您上面的代码示例无法直接运行（有错别字）。

如scipy.linalg.expm文档所述，该算法来自Al-Mohy和Higham（2009年），与Octave中的简单比例和平方Pade不同。

因此，我从Octave获得的结果略有不同，尽管结果在矩阵范数（1,2,inf）中非常接近eps。 MATLAB使用Higham（2005）的Pade方法，似乎会给出与上面的Scipy相同的结果。