Python NumPy与Octave/MATLAB精度

Question

这个问题是关于使用NumPy与Octave/MATLAB进行计算的精度(下面的MATLAB代码仅用Octave进行了测试).我知道Stackoverflow上有一个类似的问题,即这个问题,但这似乎与我在下面提出的问题有些不同.

建立

一切都在Ubuntu 14.04上运行.

Python版本3.4.0.

针对OpenBLAS编译的NumPy版本1.8.1.

针对OpenBLAS编译的Octave版本3.8.1.

示例代码

示例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).

EDITS

我的scipy版本是0.14.0.我知道Octave/MATLAB使用Pade近似方案,并熟悉这个算法.我不确定是什么scipy,但我们可以尝试以下方法.

1. 用numpy's eig或eigh(在我们的例子中后者工作正常,因为矩阵是Hermitian)对矩阵进行对角化.结果我们得到两个矩阵:一个对角矩阵D和一个矩阵U,D由对角线上原始矩阵的特征值组成,U并由相应的特征向量作为列组成; 这样原始矩阵就是由U.T.dot(D).dot(U).

2. Exponentiate D(现在这很简单,因为D是对角线).

3. 现在,如果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似乎做得更好.

Answer 1

以下代码

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中更简单的scale-and-square-Pade不同.

因此,我从Octave获得的结果略有不同,尽管结果在矩阵范数(1,2,inf)中接近eps.MATLAB使用Higham(2005)的Pade方法,它似乎与上面的Scipy给出了相同的结果.