NumPy einsumで上三角要素のみを処理する

Question

私はnumpy einsumを使用して、shape（3、N）の列ベクトルptsの配列のドット積を計算します。 N、N）をすべてのドット積と置き換えます。これは私が使用するコードです：

dotps = np.einsum('ij,ik->jk', pts, pts) 

これは動作しますが、私は主対角線の上の値が必要です。すなわち、結果の上三角部分は対角線なし。 einsumでこれらの値だけを計算することは可能ですか？または行列全体を計算するためにeinsumを使うよりも速い他の方法でもよいでしょうか？

私のpts配列はかなり大きいので、必要な値だけを計算すると計算速度が倍になります。

Answer 1

あなたは、関連する列をスライスして、np.einsum使用することができます - さらに最適化

In [109]: N = 5 
    ...: pts = np.random.rand(3,N) 
    ...: dotps = np.einsum('ij,ik->jk', pts, pts) 
    ...: 

In [110]: dotps 
Out[110]: 
array([[ 0.26529103, 0.30626052, 0.18373867, 0.13602931, 0.51162729], 
     [ 0.30626052, 0.56132272, 0.5938057 , 0.28750708, 0.9876753 ], 
     [ 0.18373867, 0.5938057 , 0.84699103, 0.35788749, 1.04483158], 
     [ 0.13602931, 0.28750708, 0.35788749, 0.18274288, 0.4612556 ], 
     [ 0.51162729, 0.9876753 , 1.04483158, 0.4612556 , 1.82723949]]) 

In [111]: R,C = np.triu_indices(N,1) 
    ...: out = np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 
    ...: 

In [112]: out 
Out[112]: 
array([ 0.30626052, 0.18373867, 0.13602931, 0.51162729, 0.5938057 , 
     0.28750708, 0.9876753 , 0.35788749, 1.04483158, 0.4612556 ]) 

---

- -

R,C = np.triu_indices(N,1) 
out = np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 

のサンプル実行を

レッツ・時間我々のアプローチをし、いずれかがありますかどうかを確認しますパフォーマンスの向上のための範囲。メモリの制約内に留まる

In [126]: N = 5000 

In [127]: pts = np.random.rand(3,N) 

In [128]: %timeit np.triu_indices(N,1) 
1 loops, best of 3: 413 ms per loop 

In [129]: R,C = np.triu_indices(N,1) 

In [130]: %timeit np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 
1 loops, best of 3: 1.47 s per loop 

我々はnp.einsumの最適化について多くを行うことができますように、それは見ていません。だから、焦点をnp.triu_indicesに移してみましょう。 N = 4については

、我々は持っている：

In [131]: N = 4 

In [132]: np.triu_indices(N,1) 
Out[132]: (array([0, 0, 0, 1, 1, 2]), array([1, 2, 3, 2, 3, 3])) 

それはしかし、ソートのシフト1のように、規則的なパターンを作成しているようです。これは、3と5の位置にシフトした累積合計で書くことができます。それは様々なN'sため

def triu_indices_cumsum(N): 

    # Length of R and C index arrays 
    L = (N*(N-1))/2 

    # Positions along the R and C arrays that indicate 
    # shifting to the next row of the full array 
    shifts_idx = np.arange(2,N)[::-1].cumsum() 

    # Initialize "shift" arrays for finally leading to R and C 
    shifts1_arr = np.zeros(L,dtype=int) 
    shifts2_arr = np.ones(L,dtype=int) 

    # At shift positions along the shifts array set appropriate values, 
    # such that when cumulative summed would lead to desired R and C arrays. 
    shifts1_arr[shifts_idx] = 1 
    shifts2_arr[shifts_idx] = -np.arange(N-2)[::-1] 

    # Finall cumsum to give R, C 
    R_arr = shifts1_arr.cumsum() 
    C_arr = shifts2_arr.cumsum() 
    return R_arr, C_arr 

レッツ・時間 - 一般的に考えると、我々はこのような何か、それをコーディング終わるだろう！

In [133]: N = 100 

In [134]: %timeit np.triu_indices(N,1) 
10000 loops, best of 3: 122 µs per loop 

In [135]: %timeit triu_indices_cumsum(N) 
10000 loops, best of 3: 61.7 µs per loop 

In [136]: N = 1000 

In [137]: %timeit np.triu_indices(N,1) 
100 loops, best of 3: 17 ms per loop 

In [138]: %timeit triu_indices_cumsum(N) 
100 loops, best of 3: 16.3 ms per loop 

したがって、それはまともなN'sのためにのように見える、triu_indices基づいてカスタマイズCUMSUMは一見の価値があるかもしれません！