加快矩陣計算（在子陣列上循環）[numpy]

我正在嘗試爲我的學生作業編寫一個算法，它運行良好。但是，計算需要很長時間，尤其是對於大數組。這部分代碼會減慢所有程序。加快矩陣計算（在子陣列上循環）[numpy]

Shapes: X.shape = mask.shape = logBN.shape = (500,500,1000), 
     F.shape = (20,20), 
     A.shape = (481,481), 
     s2 -- scalar.

我應該如何改變這種代碼，使其更快？

h = F.shape[0] 
w = F.shape[1] 
q = np.zeros((A.shape[0], A.shape[1], X.shape[2])) 
for i in range(A.shape[0]): 
    for j in range(A.shape[1]): 
     mask[:,:,:] = 0 
     mask[i:i + h,j:j + w,:] = 1 
     q[i,j,:] = ((logBN*(1 - mask)).sum(axis=(0,1)) + 
        (np.log(norm._pdf((X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])/s2)/s2)).sum(axis=(0,1))

來源

2016-11-01 Acapello

這是不完整的 - 將所有變量（F，A，X）放在一起，以便人們可以使用某些東西。如果在數組上迭代，通常最好轉換爲python列表，因爲它非常慢 - 使用向量操作最快。 – kabanus

@kabanus我不能在程序工作期間生成它們。 – Acapello

我建議打印一次，並在開始時粘貼結果。 – kabanus

只是試圖讓你的內環感

mask[:,:,:] = 0 
    mask[i:i + h,j:j + w,:] = 1 
    q[i,j,:] = ((logBN*(1 - mask)).sum(axis=(0,1)) + 
       (np.log(norm._pdf((X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])/s2)/s2)).sum(axis=(0,1))

看起來像

idx = (slice(i,i+h), slice(j,j_w), slice(None)) 
mask = np.zeros(X.shape) 
mask(idx) = 1 
mask = 1 - mask 
# alt mask=np.ones(X.shape);mask[idx]=0 
term1 = (logBN*mask).sum(axis=(0,1)) 
term2 = np.log(norm._pdf((X[idx] - F[...,None])/s2)/s2).sum(axis=(0,1)) 
q[i,j,:] = term1 + term2

所以idx和mask在A定義一個子陣列。您正在陣列外部使用logBN;裏面還有term。您正在對第1個2個dim進行求和，因此term1和term2的形狀爲X.shape[2]，您可以將其保存在q中。

該面具/窗口是20x20。

作爲第一次切割，我試圖一次計算所有i,j的term2。這看起來像一個典型的滑動窗口問題。我也嘗試將term1表示爲減法 - 整個logBN減去此窗口。

來源

2016-11-01 17:42:48 hpaulj

後通過log，exp，power代數運算重雜耍，這一切來到這個 -

# Params 
m,n = F.shape[:2] 
k1 = 1.0/(s2*np.sqrt(2*np.pi)) 
k2 = -0.5/s2**2 
k3 = np.log(k1)*m*n 

out = np.zeros((A.shape[0], A.shape[1], X.shape[2])) 
for i in range(A.shape[0]): 
    for j in range(A.shape[1]): 
     mask[:] = 1 
     mask[i:i + h,j:j + w,:] = 0 
     XF = (X[i:i + h,j:j + w,:]-F[:,:,np.newaxis])   
     p1 = np.einsum('ijk,ijk->k',logBN,mask) 
     p2 = k2*np.einsum('ijk,ijk->k',XF,XF) 
     out[i,j,:] = p1 + p2 
out += k3

使用幾件事情是 -

1] norm._pdf基本上是：norm.pdf(x) = exp(-x**2/2)/sqrt(2*pi)。所以，我們可以內嵌實現並在腳本級別對其進行優化。

2]由標量劃分將不會有效，所以這些被它們的倒數乘法代替。所以，在進入循環之前，作爲一個預處理存儲他們的倒數。

來源

2016-11-01 18:25:42 Divakar

這就是爐排！謝謝！我有想法改變一個算法，但是你用'einsum'的解決方案工作速度快8倍。雖然我不明白爲什麼它運行速度更快，普通的numpy功能。 – Acapello

@Acapello相信我，我不得不通過代數函數，因此在帖子頂部的評論。是一個艱苦的會議！ :)這更快，因爲我們正在做更小的操作。 – Divakar

加快矩陣計算（在子陣列上循環）[numpy]

回答

相關問題