用兩個numpy向量中的元素對來填充矩陣的最快方法？

Question

我有兩個1維numpy向量va和vb，它們被用來通過將所有對組合傳遞給函數來填充矩陣。

na = len(va) 
nb = len(vb) 
D = np.zeros((na, nb)) 
for i in range(na): 
    for j in range(nb): 
     D[i, j] = foo(va[i], vb[j]) 

既然這樣，這段代碼需要很長的時間來運行由於以下事實：VA和VB是相對大的（4626和737）。不過，我希望這可以得到改善，因爲使用來自scipy的cdist方法執行類似的過程，並且性能非常好。

D = cdist(va, vb, metric) 

我明明知道SciPy的具有運行這段代碼用C而不是蟒蛇的好處 - 但我希望有一些numpy的功能即時不知道，可以快速執行此。

Answer 1

之一中分配內存文檔稱functional programming routines爲np.frompyfunc的最不爲人知的numpy功能。這從Python函數創建了一個numpy ufunc。沒有其他的東西可以模擬出一種numpy的ufunc，但是它是一個適當的ufunc，它具有所有的花裏胡哨的功能。雖然行爲是在很多方面非常相似，np.vectorize，它有一些獨特的優勢，有希望下面的代碼應突出：

In [2]: def f(a, b): 
    ...:  return a + b 
    ...: 

In [3]: f_vec = np.vectorize(f) 

In [4]: f_ufunc = np.frompyfunc(f, 2, 1) # 2 inputs, 1 output 

In [5]: a = np.random.rand(1000) 

In [6]: b = np.random.rand(2000) 

In [7]: %timeit np.add.outer(a, b) # a baseline for comparison 
100 loops, best of 3: 9.89 ms per loop 

In [8]: %timeit f_vec(a[:, None], b) # 50x slower than np.add 
1 loops, best of 3: 488 ms per loop 

In [9]: %timeit f_ufunc(a[:, None], b) # ~20% faster than np.vectorize... 
1 loops, best of 3: 425 ms per loop 

In [10]: %timeit f_ufunc.outer(a, b) # ...and you get to use ufunc methods 
1 loops, best of 3: 427 ms per loop 

因此，儘管它仍然是明顯不如正確矢量實現，這是一個有點更快（循環是在C中，但你仍然有Python函數調用開銷）。

Answer 2

cdist是快的，因爲它是寫在高度優化的C代碼（如你已經指出的那樣），和它僅支持小組預定義的metric秒。

由於您希望將該操作一般地應用於任何給定的foo函數，因此別無選擇，只能調用該函數na -times- nb次。這部分不太可能進一步優化。

剩下的要優化的是循環和索引。一些建議嘗試：

1. 使用xrange而不是range（如果python2.x在python3，範圍已經是一個發電機等）
2. 使用enumerate，而不是範圍+明確索引
3. 使用蟒蛇速度「magic」，例如cython或numba，來加速循環過程。

如果您可以對foo作進一步的假設，則可以進一步加速。

Answer 3

像@ shx2說的，這一切都取決於什麼是foo。如果你能在numpy的ufuncs方面表達出來，然後用outer方法：

In [11]: N = 400 

In [12]: B = np.empty((N, N)) 

In [13]: x = np.random.random(N) 

In [14]: y = np.random.random(N) 

In [15]: %%timeit 
for i in range(N): 
    for j in range(N): 
    B[i, j] = x[i] - y[j] 
    ....: 
10 loops, best of 3: 87.2 ms per loop 

In [16]: %timeit A = np.subtract.outer(x, y) # <--- np.subtract is a ufunc 
1000 loops, best of 3: 294 µs per loop 

否則，您可以推動循環下來用Cython水平。繼續以上一個簡單的例子：

In [45]: %%cython 
cimport cython 
@cython.boundscheck(False) 
@cython.wraparound(False) 
def foo(double[::1] x, double[::1] y, double[:, ::1] out): 
    cdef int i, j 
    for i in xrange(x.shape[0]): 
     for j in xrange(y.shape[0]): 
      out[i, j] = x[i] - y[j] 
    ....: 

In [46]: foo(x, y, B) 

In [47]: np.allclose(B, np.subtract.outer(x, y)) 
Out[47]: True 

In [48]: %timeit foo(x, y, B) 
10000 loops, best of 3: 149 µs per loop 

的用Cython例如故意製造過於簡單化：在現實中，你可能要添加一些形狀/步幅檢查，你的函數等