3

我正在将我在 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)

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

4

3 回答 3

3

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

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

于 2013-01-10T22:36:55.920 回答
2

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

作为傅立叶变换的 FFT 具有许多属性，它们会产生大多数（全部）不必要的 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

不过，它看起来确实要整洁得多。

于 2013-01-11T00:07:23.453 回答
1

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

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

于 2013-01-10T23:16:20.997 回答