Python NumPy与Octave/MATLAB精度对比

Question

这个问题是关于使用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对矩阵进行对角化（在我们的例子中，后者效果很好，因为矩阵是厄米矩阵）。结果我们得到两个矩阵：一个对角矩阵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在这种情况下似乎做得更好。

numpy matlab 算法优化特征值 Octave 计算差异精度对比对角化

Python NumPy与Octave/MATLAB精度对比

1 个回答

撰写回答