矢量化函数的numpy广播答案

【问题标题】：numpy broadcast of vectorized function矢量化函数的numpy广播
【发布时间】：2020-12-02 08:58:24
【问题描述】：

我有一个函数 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]，它可以工作，但为什么它可以工作？

我已经检查了this 和this，但仍然无法弄清楚。

【问题讨论】：

标签： python arrays numpy

【解决方案1】：

我们可以用添加说明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])

【讨论】：

我看到了你的回答，但已经写了很多文字，所以我已经发布了。

【解决方案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)--|

【讨论】：