Python 运行速度比 MATLAB 慢答案

【问题标题】：Python running slower than MATLABPython 运行速度比 MATLAB 慢
【发布时间】：2012-12-25 10:56:53
【问题描述】：

我正在将我在 MATLAB 中编写的数值方法转换为 Python。由于某种原因，几乎相同的 Python 代码运行速度要慢得多。这里U 和V 是在每个时间步求解的未知数。 U[:,n] 和 V[:,n] 的大小为 700x1。其余变量（dt、A 和 denom）是常量。这是循环（numpy 已导入为*）：

for n in range(0, 400):
    UnVn2 = fft.fft(U[:, n] * V[:, n] ** 3)
    U[:, n +1 ] = fft.ifft((fft.fft(U[:, n]) / dt - UnVn2 + A) / denom)
    V[:, n + 1] = fft.ifft((fft.fft(V[:, n]) / dt + UnVn2) / denom)

有什么建议吗？非常感谢。

【问题讨论】：

什么python运行时？是否使用 JIT 编译器？ BLAS 和 FFT 例程是否与使用 MatLab 打包的例程一样快？ IIRC，MatLab 使用 Intel 和 AMD 二进制文件，这些二进制文件经过了非常严格的优化。
我对 Python 还是很陌生，但我只是通过将文件导入 IDLE 来运行代码。那叫JIT吗？ Python 运行时似乎至少是 MATLAB 的 3 倍。我曾假设 numpy 例程将被完全优化，因为它非常常用于科学计算。
numpy 能够使用这些优化的例程。请参阅software.intel.com/en-us/articles/numpyscipy-with-intel-mkl，但请使用最适合您的处理器的优化库。

标签： performance numpy python-3.x fft

【解决方案1】：

有关使 python 和 numpy 使用与 MATLAB 一起提供的相同加速 FFT 例程的说明，请参阅 this。

如果您有 AMD 处理器，请参阅 these 说明。

【讨论】：

【解决方案2】：

我不确定为什么 Python 比 Matlab 慢，但是...

作为傅立叶变换的 FFT 有许多 properties，这会产生大多数（所有）不必要的 FFT 操作：

def func1(U, V, dt, denom, A) :
    UnVn2 = np.fft.fft(U * V**3)
    U_ = np.fft.ifft((np.fft.fft(U) / dt - UnVn2 + A) / denom)
    V_ = np.fft.ifft((np.fft.fft(V) / dt + UnVn2) / denom)
    return np.vstack((U_, V_))

def func2(U, V, dt, denom, A) :
    UnVn2 = U * V**3
    U_ = (U / dt - UnVn2) / denom
    U_[0] += A / denom
    V_ = (V / dt + UnVn2) / denom
    return np.vstack((U_, V_))

U = np.random.rand(700)
V = np.random.rand(700)
dt, denom, A = tuple(np.random.rand(3))

>>> func1(U, V, dt, denom, A)
array([[ 2.35201751 -1.11022302e-16j,  0.81099082 -2.45463372e-16j,
         0.48451858 +2.15658782e-18j, ...,  2.23237712 -5.24753851e-16j,
         1.15264205 -2.31140087e-16j,  1.06670009 +1.28369537e-16j],
       [ 2.89314136 +8.67361738e-17j,  3.65612404 -7.80625564e-17j,
         3.31383830 +8.96916836e-17j, ...,  0.90415910 +6.27969898e-16j,
         3.03505664 +4.72358723e-16j,  0.64669863 +4.99600361e-16j]])
>>> func2(U, V, dt, denom, A)
array([[ 2.35201751,  0.81099082,  0.48451858, ...,  2.23237712,
         1.15264205,  1.06670009],
       [ 2.89314136,  3.65612404,  3.3138383 , ...,  0.9041591 ,
         3.03505664,  0.64669863]])
>>> np.max(np.abs(func1(U, V, dt, denom, A) - func2(U, V, dt, denom, A)))
1.5151595604785605e-15

当然还有：

>>> import timeit
>>> timeit.timeit('func1(U, V, dt, denom, A)', 'from __main__ import func1, U, V, dt, denom, A', number=400)
0.14169366197616284
>>> timeit.timeit('func2(U, V, dt, denom, A)', 'from __main__ import func2, U, V, dt, denom, A', number=400)
0.06098524703428154

我不得不承认这比我预期的要少，但它仍然快了将近 3 倍。

编辑不做 FFT 的速度似乎太小了，所以我修改了 func1 和 func2 以返回一个带有 (U_, V_) 的元组并运行以下代码：

from time import clock
U = np.zeros((700,400), dtype=np.float)
V = np.zeros((700,400), dtype=np.float)
U[:,0] = np.random.rand(700)
V[:,0] = np.random.rand(700)
dt, denom, A = tuple(np.random.rand(3))
t = clock()
for j in xrange(399) :
    U[:, j+1], V[:, j+1] = func1(U[:, j], V[:, j], dt, denom, A)
print clock() - t
t = clock()
for j in xrange(399) :
    U[:, j+1], V[:, j+1] = func2(U[:, j], V[:, j], dt, denom, A)
print clock() - t

打印输出是11.5148652438和0.321673111194，所以实际问题设置中的加速更像是x30。

我还计时了 pwuertz 的提案，没有明显改进，11.1805414552 和 0.297830755317 为以下代码：

U = np.zeros((400, 700), dtype=np.float)
V = np.zeros((400, 700), dtype=np.float)
U[0] = np.random.rand(700)
V[0] = np.random.rand(700)
dt, denom, A = tuple(np.random.rand(3))
t = clock()
for j in xrange(399) :
    U[j+1], V[j+1] = func1(U[j], V[j], dt, denom, A)
print clock() - t
t = clock()
for j in xrange(399) :
    U[j+1], V[j+1] = func2(U[j], V[j], dt, denom, A)
print clock() - t

不过，它看起来确实非常整洁。

【讨论】：

摆脱 fft 无疑是一个很大的进步。然而，内存对齐的影响取决于阵列的大小和硬件。在我的系统上，切换轴时我看到了大约 30% 的改进，对于较大的阵列 (2000x2000)，改进了 2 倍。记住这样的事情只是一般建议；）

【解决方案3】：

我不确定 MatLab 如何在多维数组中组织轴，但我很确定 numpy 使用类似 C 的 row-major order（编辑：维基百科甚至提到 MatLab 使用列优先顺序；）） .

由于您在单列上进行操作，因此您的所有操作都必须遍历行。对于行优先排序，这通常比遍历整行效率低。考虑转置二维数组的布局，您应该会获得显着的性能提升。

【讨论】：