的有效方式，加快一些numpy的操作

Question

我試圖找到一個有效的代碼，而不是下面的代碼（也就是我的代碼只是其中的一部分），以提高速度：

for pr in some_list: 
    Tp = T[partition[pr]].sum(0) 
    Tpx = np.dot(Tp, xhat) 
    hp = h[partition[[pr]].sum(0) 
    up = (uk[partition[pr][:]].sum(0))/len(partition[pr]) 
    hpu = hpu + np.dot(hp.T, up) 
    Tpu = Tpu + np.dot(Tp.T, up) 

我至少有兩個類似的代碼塊。正如你所看到的，我使用了三次幻想索引（真的找不到另一種方式）。在我的算法中，我需要很快完成這個部分，但現在不會發生。我會非常感謝任何建議。

謝謝大家。

最佳，

Answer 1

如果你的分區很少，而且有許多元件分別，你應該考慮在你的對象的索引交換。沿其第二維度總結形狀(30,1000)的陣列應比沿其第一尺寸求和形狀(1000,30)的陣列更快，因爲在前者的情況下，你總是求和對於每個剩餘的索引的存儲器連續的塊（即arr[k,:]每個k）。所以，如果你把求和指數放在最後（並且在你處理時去掉一些尾隨的單身維度），那麼你可能會加快速度。

由於hpaulj noted in a comment，目前尚不清楚你的循環如何可以量化。但是，由於它對性能至關重要，因此您仍然可以嘗試向量化某些的工作。

我建議您爲每個分區（在預分配之後）存儲hp,up和Tp，然後在單個矢量化步驟中執行標量/矩陣乘積。還要注意的是Tpx是你的榜樣未使用的，所以我在這裏省略它（無論你用它做，你同樣可以做到這一點其他的例子）：

part_len = len(some_list) # number of partitions, N 
Tpshape = (part_len,) + T.shape[1:] # (N,30,100) if T was (1000,30,100) 
hpshape = (part_len,) + h.shape[1:] # (N,30,1) if h was (1000,30,1) 
upshape = (part_len,) + uk.shape[1:] # (N,30,1) if uk was (1000,30,1) 
Tp = np.zeros(Tpshape) 
hp = np.zeros(hpshape) 
up = np.zeros(upshape) 

for ipr,pr in enumerate(some_list): 
    Tp[ipr,:,:] = T[partition[pr]].sum(0) 
    hp[ipr,:,:] = h[partition[[pr]].sum(0) 
    up[ipr,:,:] = uk[partition[pr]].sum(0)/len(partition[pr]) 

# compute vectorized dot products: 
#Tpx unclear in original, omitted 
# sum over second index (dot), sum over first index (sum in loop) 
hpu = np.einsum('abc,abd->cd',hp,up) # shape (1,1) 
Tpu = np.einsum('abc,abd->cd',Tp,up) # shape (100,1) 

顯然的關鍵球員是numpy.einsum。 ，當然，如果hpu和Tpu有循環之前的一些現有的價值觀，你有從上面einsum結果來增加這些值。

至於einsum，它執行求和和任意尺寸的陣列中的收縮。以上apearing，'abc,abd->cd'圖案，當施加到三維陣列A和B，將返回的2D陣列C，用下面的定義（數學僞代碼）：

C(c,d) = sum_a sum_b A(a,b,c)*B(a,b,d) 

對於給定的固定a求和索引，裏面有什麼是

sum_b A(a,b,c)*B(a,b,d) 

其中，如果c和d指數保持，將euqivalent到np.dot(A(a,:,:).T,B(a,:,:))。因爲我們這些矩陣相對於太總和爲a，我們應該做你的多圈版本究竟是幹什麼的，加起來的總金額中的每個np.dot()貢獻。