矢量化函数的numpy广播

Question

我有一个函数f，它计算2个向量的内积样值(该函数是对称的，所以f(x, y) = f(y, x)(。它采用2个相同大小的1d数组并输出单个值。现在我有一个形状为(n, d)的数据矩阵X(每行表示一个输入向量(，并且我想计算一个具有形状(n, n)的矩阵K，使得K[i][j] = f(X[i], X[j])。一种方法是使用以下代码。

import numpy as np
# ... some other code
K = np.zeros((n, n))
for i in range(n):
for j in range(n):
K[i][j] = f(X[i], X[j])

这很管用，但既不优雅也不表演。另一种方法是使用numpy 中的vectorize

import numpy as np
# ... some other code
vf = np.vectorize(f, signature="(n),(n)->()")
# this does not work (output shape (n,))
#K = vf(X, X)
# this works (output shape (n, n))
K = vf(X, X[: np.newaxis])
# this also works (output shape (n, n), result same as above)
#K = vf(X[: np.newaxis], X)

如果我们放入一个(n, d)的数组和另一个(n, d)的数组，它将输出一个一维的(n,)数组，本质上是将两个矩阵的n行中的每一行作为输入，从而产生n输出。然而，如果我们将其转换为形状为(n, 1, d)的X[: np.newaxis]，它是有效的，但为什么它有效呢？

我已经检查过这个和这个了，但还是搞不清楚。

Answer 1

我们可以用add:来说明这个broadcasting

In [111]: x = np.arange(12).reshape(3,4)
In [112]: x
Out[112]: 
array([[ 0,  1,  2,  3],
[ 4,  5,  6,  7],
[ 8,  9, 10, 11]])
In [113]: x + x[:,None]
Out[113]: 
array([[[ 0,  2,  4,  6],
[ 4,  6,  8, 10],
[ 8, 10, 12, 14]],
[[ 4,  6,  8, 10],
[ 8, 10, 12, 14],
[12, 14, 16, 18]],
[[ 8, 10, 12, 14],
[12, 14, 16, 18],
[16, 18, 20, 22]]])
In [114]: _.shape
Out[114]: (3, 3, 4)
In [115]: (x + x[:,None]).sum(axis=-1)
Out[115]: 
array([[12, 28, 44],
[28, 44, 60],
[44, 60, 76]])

广播进行到：(3，4(和(3，1，4(=>(1,3,4(和(3,1,4(=>(3,3,4(。然后我们在最后一个轴上求和，得到(3,3(

做和你一样的事情：

In [116]: vf = np.vectorize(lambda a,b:(a+b).sum(), signature="(n),(n)->()")
In [117]: vf(x, x[:,None])
Out[117]: 
array([[12, 28, 44],
[28, 44, 60],
[44, 60, 76]])

但如前所述，np.vectorize不是一个速度工具，尤其是在使用signature:时

In [118]: timeit (x + x[:,None]).sum(axis=-1)
13.8 µs ± 363 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [119]: timeit vf(x, x[:,None])
278 µs ± 7.68 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

不将维度添加到第二个x:

In [120]: vf(x, x)
Out[120]: array([12, 44, 76])
In [121]: (x + x).sum(axis=-1)
Out[121]: array([12, 44, 76])

Answer 2

除了@hpaulj的答案之外的另一个例子

广播可以用来有效地向量化代码。让我们举一个例子。我有一个数组(10,2(，我想取每个点之间的欧几里得距离。

X = np.random.uniform(-5,5,(10,2))

一种方法是使用支持广播的NumPy函数来广播我的欧几里得距离的完整公式。

np.sqrt(np.sum(np.square(np.subtract(X[:,None,:], X[None,:,:])), axis=-1))
#This returns a (10,10) matrix, with euclidean distance of each of the 10 points to the other

类似地，您可以做的另一件事是将函数矢量化，使其行为类似于numpy函数并支持广播。

def euclidean(point1, point2):
return ((point1[0]-point2[0])**2 + (point1[1]-point2[1])**2)**(1/2)
euclidean_v = np.vectorize(euclidean, signature="(n),(n)->()")
euclidean_v(X[:,None,:], X[None,:,:])
#Again returns a (10,10)

在这两种情况下-

|--- (10,1,2)--|
(10,2)---|              |---(10,10,2)--->(10,10)
|--- (1,10,2)--|