用两个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例如故意制造过于简单化：在现实中，你可能要添加一些形状/步幅检查，你的函数等