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のような行優先順位を使用していることは確かです（編集：Wikipediaは、MatLabが列優先順序を使用しているとさえ述べています;））。

単一の列で操作しているため、すべての操作のみが行を反復処理する必要があります。行優先の順序付けでは、これは通常、行全体を反復処理するよりも効率的ではありません。2Dアレイのレイアウトを入れ替えることを検討してください。そうすれば、パフォーマンスが著しく向上するはずです。

于 2013-01-10T23:16:20.997 に答える