Pythonのループがnumpyの配列操作よりも速くなる

Question

いくつかのフーリエ変換に基づいたモデルのコードを作成しています。現在、私はその一部を最適化しようとしているので、いくらかの大量のデータでも使用できるようになっています。それをしようとしている間、私は奇妙な振る舞いを見つけました、主に私のコードのループバージョンは、numpyで書かれた同じコードより速くなっています。次のようにテストコードは次のとおりです。

# -*- coding: utf-8 -*- 
import numpy as np 
import timeit 


def fourier_coef_loop(ts, times, k_max): 
    coefs = np.zeros(k_max, dtype=float) 
    t = 2.0 * np.pi * (times - times[0])/times[-1] 
    x = np.dot(np.arange(1,k_max+1)[np.newaxis].T,t[np.newaxis]) 
    for k in xrange(1, k_max + 1): 
     cos_k = np.cos(x[k-1]) 
     coefs[k - 1] = (ts[-1] - ts[0]) + (ts[:-1] * (cos_k[:-1] - cos_k[1:])).sum() 
    return(coefs) 


def fourier_coef_np(ts, times, k_max): 
    coefs = np.zeros(k_max, dtype=float) 
    t = 2.0 * np.pi * (times - times[0])/times[-1] 
    x = np.dot(np.arange(1,k_max+1)[np.newaxis].T,t[np.newaxis]) 
    coefs = np.add(np.einsum('ij,j->i',np.diff(np.cos(x)), -ts[:-1]), (ts[-1] - ts[0])) 
    return(coefs) 


if __name__ == '__main__': 
    iterations = 10 
    size = 20000 
    setup = "from __main__ import fourier_coef_loop, fourier_coef_np, size\n" \ 
      "import numpy as np" 

# arg = np.random.normal(size=size) 
# print(np.all(fourier_coef_np(arg, np.arange(size,dtype=float), size/2) == fourier_coef_loop(arg, np.arange(size,dtype=float), size/2))) 

    time_loop = timeit.timeit("fourier_coef_loop(np.random.normal(size=size), np.arange(size,dtype=float), size/2)", 
           setup=setup, number=iterations) 
    print("With loop: {} s".format(time_loop)) 

    time_np = timeit.timeit("fourier_coef_np(np.random.normal(size=size), np.arange(size,dtype=float), size/2)", 
          setup=setup, number=iterations) 
    print("With numpy: {} s".format(time_np)) 

それは次のような結果を与える：ループバージョンはより速く、純粋にnumpyのバージョンよりも、なぜ

With loop: 60.8385488987 s 
With numpy: 64.9192998409 s 

誰かが私に教えていただけますか？私は完全にアイデアがなくなった。この特定の機能をより速くする方法に関する提案もありがとう。

Answer 1

私がコメントしたところでは、ループとベクトル化されたバージョンの違いに関係しているとは思っていません。20000の正規擬似乱数の生成も含まれています。正確な画像を取得したい場合は、タイミングの外にできるだけ多くの設定を行うようにしてください。

とにかく、あなたのコードは、いくつかの奇妙なポイントを持っているので、ここでは私自身の提案です：

def fourier_coef_new(ts,times,k_max): 
    # no need to pre-allocate coefs, you're rebinding later 
    t = 2.0 * np.pi * (times - times[0])/times[-1] 
    x = np.arange(1,k_max+1)[:,None] * t # make use of array broadcasting 
    coefs = -np.dot(np.diff(np.cos(x)),ts[:-1]) + ts[-1]-ts[0] # einsum was dot 
    return(coefs) 

は私が入力のセットに対してこの機能をテストし、それはあなたの関数と同じ結果を与えました。配列にシングルトンディメンションを導入するには、[:,None]（またはそれに相当する[:,np.newaxis]）の方法に注意してください。 (N,1)の形状の配列と形状(M,)（後者は(1,M)と互換性があります）のいずれかを取得すると、それらの製品はnumpyの配列ブロードキャスト規則に従って(N,M)になります。

あらかじめ生成されたデータセットを使用して、3つの機能のすばやいタイミングをとりましたが、1ではPython 3,2になりました。size = 2000で、3.IPythonの%timeitビルトインを使用しています。私は、これらの結果は、任意のより信頼性の高いあなたよりであることを主張していないんだけど、私は上記のバージョンは最速であることを疑う：

In [37]: %timeit fourier_coef_loop(ts,times,k_max) 
1000 loops, best of 3: 1.09 ms per loop 

In [38]: %timeit fourier_coef_np(ts,times,k_max) 
1000 loops, best of 3: 1.06 ms per loop 

In [39]: %timeit fourier_coef_new(ts,times,k_max) 
1000 loops, best of 3: 1.05 ms per loop 

あなたが見ることができるように、numpyのバージョンはわずかに速くなるように見えます。 2つの時限の間でコードの小さな部分だけが異なるので（両方の場合に重い三角関数が含まれます）、これは妥当と思われます。