Vectorise Python代碼

Question

我編寫了一個克里格算法，但我覺得它很慢。特別是，你對我怎麼能vectorise的一段代碼在利弊的想法如下功能：

import time 
import numpy as np 

B = np.zeros((200, 6)) 
P = np.zeros((len(B), len(B))) 

def cons(): 
    time1=time.time() 
    for i in range(len(B)): 
    for j in range(len(B)): 
     P[i,j] = corr(B[i], B[j]) 
    time2=time.time() 
    return time2-time1 

def corr(x,x_i): 
    return np.exp(-np.sum(np.abs(np.array(x) - np.array(x_i))))  

time_av = 0. 
for i in range(30): 
    time_av+=cons() 
print "Average=", time_av/100. 

---

編輯：獎金問題

1. 發生了廣播解決方案如果我想要corr(B[i], C[j])與C相同的維度比B

2. 如果我的p範數秩是一個數組，它會發生什麼scipy解決方案：

   p=np.array([1.,2.,1.,2.,1.,2.]) 
   def corr(x, x_i): 
       return np.exp(-np.sum(np.abs(np.array(x) - np.array(x_i))**p)) 
   

   對於2.，我試過P = np.exp(-cdist(B, C,'minkowski', p))，但scipy期待一個標量。

Answer 1

你的問題似乎很簡單矢量化。對於每對B行要計算

P[i,j] = np.exp(-np.sum(np.abs(B[i,:] - B[j,:]))) 

你可以利用陣列廣播，並介紹了第三個維度，沿着過去的一個總結：

P2 = np.exp(-np.sum(np.abs(B[:,None,:] - B),axis=-1)) 

的想法是重塑第一次出現B以塑造(N,1,M)，而第二個B留下形狀(N,M)。與陣列廣播，後者相當於(1,N,M)，所以

B[:,None,:] - B 

是形狀(N,N,M)的。沿最後一個索引進行求和將導致您正在查找的(N,N)-形相關數組。

---

請注意，如果你使用scipy，你將能夠做到這一點使用scipy.spatial.distance.cdist（或等價的，scipy.spatial.distance.pdist和scipy.spatial.distance.squareform的組合），而不必計算這個的Symmetrix矩陣的下三角一半。在註釋中使用@Divakar的建議爲最簡單的辦法是這樣的：

from scipy.spatial.distance import cdist 
P3 = 1/np.exp(cdist(B, B, 'minkowski',1)) 

cdist將計算在1範數，這正是座標差的絕對值之和閔可夫斯基距離。