需要幫助的矢量化代碼或優化

Question

我想通過第一次內插數據來做一個雙積分做一個表面。我正在使用numba來加速這個過程，但這只是需要很長時間。

Here is my code,帶有運行位於here和here的代碼所需的圖像。

Answer 1

注意到你的代碼有一個for循環的四重嵌套集，我專注於優化內對。這裏的舊代碼：

for i in xrange(K.shape[0]): 
    for j in xrange(K.shape[1]): 

     print(i,j) 
     '''create an r vector ''' 
     r=(i*distX,j*distY,z) 

     for x in xrange(img.shape[0]): 
      for y in xrange(img.shape[1]): 
       '''create an ksi vector, then calculate 
        it's norm, and the dot product of r and ksi''' 
       ksi=(x*distX,y*distY,z) 
       ksiNorm=np.linalg.norm(ksi) 
       ksiDotR=float(np.dot(ksi,r)) 

       '''calculate the integrand''' 
       temp[x,y]=img[x,y]*np.exp(1j*k*ksiDotR/ksiNorm) 

     '''interpolate so that we can do the integral and take the integral''' 
     temp2=rbs(a,b,temp.real) 
     K[i,j]=temp2.integral(0,n,0,m) 

由於K和IMG每個約爲2000×2000，最裏面的語句需要執行160000億次。這對於使用Python來說並不實用，但我們可以使用NumPy將工作轉換爲C和/或Fortran以進行矢量化。我一次只做了一步，試圖確保結果一致;這裏是我結束了：

'''create all r vectors''' 
R = np.empty((K.shape[0], K.shape[1], 3)) 
R[:,:,0] = np.repeat(np.arange(K.shape[0]), K.shape[1]).reshape(K.shape) * distX 
R[:,:,1] = np.arange(K.shape[1]) * distY 
R[:,:,2] = z 

'''create all ksi vectors''' 
KSI = np.empty((img.shape[0], img.shape[1], 3)) 
KSI[:,:,0] = np.repeat(np.arange(img.shape[0]), img.shape[1]).reshape(img.shape) * distX 
KSI[:,:,1] = np.arange(img.shape[1]) * distY 
KSI[:,:,2] = z 

# vectorized 2-norm; see http://stackoverflow.com/a/7741976/4323              
KSInorm = np.sum(np.abs(KSI)**2,axis=-1)**(1./2) 

# loop over entire K, which is same shape as img, rows first               
# this loop populates K, one pixel at a time (so can be parallelized)            
for i in xrange(K.shape[0]):                      
    for j in xrange(K.shape[1]):                     

     print(i, j) 

     KSIdotR = np.dot(KSI, R[i,j]) 
     temp = img * np.exp(1j * k * KSIdotR/KSInorm) 

     '''interpolate so that we can do the integral and take the integral''' 
     temp2 = rbs(a, b, temp.real) 
     K[i,j] = temp2.integral(0, n, 0, m) 

內部對循環現在已經完全消失了，提前完成矢量操作替代（在輸入端的大小的空間成本直線）。

這樣可以在不使用Numba的情況下，將我的Macbook Air 1.6 GHz i5上的外部兩個循環的每次迭代的時間從340秒減少到1.3秒。在每次迭代1.3秒中，0.68秒用於rbs函數，即scipy.interpolate.RectBivariateSpline。有可能有進一步優化的空間 - 這裏有一些想法：

1. Reenable Numba。我沒有在我的系統上。它在這一點上可能沒有太大的區別，但很容易讓你測試。
2. 做更多特定領域的優化，例如試圖簡化正在完成的基本計算。我的優化旨在無損，並且我不知道您的問題域，因此無法儘可能深入地進行優化。
3. 嘗試矢量化剩餘的循環。這可能很難，除非你願意用每次調用多次計算的東西來替換scipy RBS函數。
4. 獲得更快的CPU。我的速度很慢;通過使用比我的微型筆記本電腦更好的計算機，你可以獲得至少2倍的加速比。
5. 對您的數據進行降採樣。您的測試圖像是2000x2000像素，但包含的細節很少。如果你將他們的線性尺寸減少2-10倍，你會得到巨大的加速。

這就是我現在的情況。這在哪裏離開你？假設一臺稍微好一點的計算機並且沒有進一步的優化工作，即使優化的代碼也需要大約一個月的時間來處理你的測試圖像。如果你只需要做一次，也許沒關係。如果您需要更頻繁地執行此操作，或者在嘗試不同的操作時需要迭代代碼，那麼您可能需要繼續優化 - 從現在消耗一半以上時間的RBS函數開始。

特別提示：您的代碼會更容易處理，如果它不具有幾乎相同的變量名稱，如k和K，也沒用過j作爲變量名，也可以作爲一個複雜的數字後綴很多（0j） 。