なぜB = numpy.dot（A、x）はB [i、：、：] = numpy.dot（A [i、：、：]、x））を実行するとループが遅くなるのですか？

Question

私は説明できない効率テストの結果を得ています。

i番目のエントリB [i、：、：] = A [i、：、：] .dot（x）を持つ行列Bをアセンブルしたいと考えており、各A [i、：、：]は2次元行列、xも同様です。

パフォーマンスをテストするために3つの方法があります。ランダム（numpy.random.randn）の行列A =（10,1000,1000）、x =（1000,1200）を作成します。

（2）Iをループと2D内積を実行

# initialize B = np.zeros([dim1, dim2, dim3]) 
    for i in range(A.shape[0]): 
     B[i,:,:] = A[i,:,:].dot(x) 

total time: 0.826 s 

（1）単一の多次元内積

B = A.dot(x) 

total time: 102.361 s 

（3）numpy.einsum

：私は、次の時間結果を得ます

B3 = np.einsum("ijk, kl -> ijl", A, x) 

total time: 8.289 s 

したがって、オプション（2）は、最も速いです。しかし、ちょうど（1）と（2）を考えると、私はそれらの間に大きな違いは見当たりません。どのようにして二次元ドット製品をループしたり、やったりするのが〜124倍速くなるのでしょうか？彼らは両方ともnumpy.dotを使用します。どんな洞察？

私はちょうど下の上記の結果のために使用されるコードが含ま：

import numpy as np 
import numpy.random as npr 
import time 

dim1, dim2, dim3 = 10, 1000, 1200 
A = npr.randn(dim1, dim2, dim2) 
x = npr.randn(dim2, dim3) 

# consider three ways of assembling the same matrix B: B1, B2, B3 

t = time.time() 
B1 = np.dot(A,x) 
td1 = time.time() - t 
print "a single dot product of A [shape = (%d, %d, %d)] with x [shape = (%d, %d)] completes in %.3f s" \ 
    % (A.shape[0], A.shape[1], A.shape[2], x.shape[0], x.shape[1], td1) 


B2 = np.zeros([A.shape[0], x.shape[0], x.shape[1]]) 
t = time.time() 
for i in range(A.shape[0]): 
    B2[i,:,:] = np.dot(A[i,:,:], x) 
td2 = time.time() - t 
print "taking %d dot products of 2D dot products A[i,:,:] [shape = (%d, %d)] with x [shape = (%d, %d)] completes in %.3f s" \ 
    % (A.shape[0], A.shape[1], A.shape[2], x.shape[0], x.shape[1], td2) 

t = time.time() 
B3 = np.einsum("ijk, kl -> ijl", A, x) 
td3 = time.time() - t 
print "using np.einsum, it completes in %.3f s" % td3 

Answer 1

は、私はあなたの場合と同じくらいではないが、同様の順位

In [355]: %%timeit 
    .....: B=np.zeros((N,M,L)) 
    .....: for i in range(N): 
       B[i,:,:]=np.dot(A[i,:,:],x) 
    .....: 
10 loops, best of 3: 22.5 ms per loop 
In [356]: timeit np.dot(A,x) 
10 loops, best of 3: 44.2 ms per loop 
In [357]: timeit np.einsum('ijk,km->ijm',A,x) 
10 loops, best of 3: 29 ms per loop 

In [367]: timeit np.dot(A.reshape(-1,M),x).reshape(N,M,L) 
10 loops, best of 3: 22.1 ms per loop 

In [375]: timeit np.tensordot(A,x,(2,0)) 
10 loops, best of 3: 22.2 ms per loop 

itererativeが速くなるに取得します。

これは、反復ディメンションが他のディメンションに比べて小さい限り、おそらく真です。この場合、反復（関数呼び出しなど）のオーバーヘッドは、計算時間に比べて小さい。一度にすべての値を実行すると、より多くのメモリが使用されます。

dotバリエーションを試しましたが、Aを2dに再構成しました。dotは内部的にそのような形に整形しています。私はそれが実際に最速であることに驚いています。 tensordotは、おそらく同じ再構成（Pythonが読める場合はそのコード）を行っています。

---

einsum 4つの変数、i,j,k,mを伴う「和製品の」反復を設定する - それはCレベルnditerとdim1*dim2*dim2*dim3手順です。したがって、より多くのインデックスを持つほど、反復スペースは大きくなります。

Answer 2

私はnumpyののC-APIとあまりにも精通していないよ、とnumpy.dotは以前に_dotblas下に使用されるものな組み込み関数ですがバージョン。

しかし、ここに私の考えがあります。

1）numpy.dotは、2次元配列とn次元配列に対して異なるパスをとります。 numpy.dotのonline documentation：それは行列の乗算に相当し、ベクトルの内積に1-Dアレイ（複素共役せず）のための2次元アレイで

。 N次元の場合は、aの最後の軸およびbの2番目の最後の軸上の合計積である。

ドット（a、b）[i、j、k、m] = sum（a [i、j 、：] * B [K、：、M]）

だから2次元配列のためにあなたが常にBLASのdgemmへの呼び出しを1つ持つことが保証されている、しかし、NDアレイのnumpyのは、アレイ用の乗算軸を選択する場合があります（私が掲載した抜粋から分かるように）最も速い変化軸に対応していない可能性があり、その結果、dgemmのフルパワーを逃してしまう可能性があります。

2）Aの配列が大きすぎてCPUキャッシュにロードできません。あなたの例では、あなたは、ほとんど80MBあなたのキャッシュのサイズよりもはるかに大きい

In [1]: A.nbytes 
80000000 
In [2]: 80000000/1024 
78125 

を与える寸法(10,1000,1000)でAを使用しています。だからもう一度dgemmの力の大部分を失う。

3）また、関数のタイミングもやや間違っています。 Pythonの関数は正確ではないことが知られています。代わりにtimeitを使用してください。

：それでは、

dim1, dim2, dim3 = 20, 20, 20 
A = np.random.rand(dim1, dim2, dim2) 
x = np.random.rand(dim2, dim3) 

def for_dot1(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[i,:,:], x) 

def for_dot2(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[:,i,:], x)  

def for_dot3(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[:,:,i], x) 

、ここでは、私は（OpenBLAS 0.2.14に対して構築numpy 1.9.2を使用して）取得のタイミングですキャッシュにロードすることができます配列を実験してみましょう、心の中で上記のすべての点を有する

In [3]: %timeit np.dot(A,x) 
10000 loops, best of 3: 174 µs per loop 
In [4]: %timeit np.einsum("ijk, kl -> ijl", A, x) 
10000 loops, best of 3: 108 µs per loop 
In [5]: %timeit np.einsum("ijk, lk -> ijl", A, x) 
10000 loops, best of 3: 97.1 µs per loop 
In [6]: %timeit np.einsum("ikj, kl -> ijl", A, x) 
1000 loops, best of 3: 238 µs per loop 
In [7]: %timeit np.einsum("kij, kl -> ijl", A, x) 
10000 loops, best of 3: 113 µs per loop 
In [8]: %timeit for_dot1(A,x) 
10000 loops, best of 3: 101 µs per loop 
In [9]: %timeit for_dot2(A,x) 
10000 loops, best of 3: 131 µs per loop 
In [10]: %timeit for_dot3(A,x) 
10000 loops, best of 3: 133 µs per loop 

時間差はまだありますが、桁違いにはないことに注意してください。 choosing the axis of multiplicationの重要性にも注意してください。おそらく、numpyの開発者は、N-Dアレイのフードの下で実際に何を行うのかについていくつかの光を当てることができます。小さい暗くなる10,100,200で