Numpy中类似于矩阵乘法的广播函数

Question

所以我正在为一些更大的数据集构建一个 KNN，并且我需要运行 Leave-one-out 交叉验证以选择正确的 K，因此速度很重要。

我正在尝试通过广播进行距离计算。

---

情况是：X 是我的训练矩阵，一个二维矩阵，行上有样本。 Q 是我的查询矩阵或测试数据，行上也有样本。

我需要运行类似于矩阵乘法的东西，其中我将 Q 的每一行与 X.T 的每一列（x 转置）匹配，并构建一个 sample x sample 矩阵，其中每个条目 [i,j] 是查询样本 i 与训练样本 j 的距离。然后我会从前 k 个样本中排序并选择类的模式。

---

无论如何，在矩阵乘法中，numpy 正是这样做的……但它不是进行距离计算，而是进行分段乘法和求和（点积）。如果我可以将距离函数插入到那个位置，我想我的 KNN 距离计算将与 numpy 矩阵乘法一样快。

有没有办法使用广播或其他一些 numpy 技术来做到这一点？

也许是一种并行化的方法？

---

示例代码：

import numpy as np

x1 = np.asarray([1.0,10.0,100.0])
x2 = np.asarray([40.0,60.0,80.0])
x3 = np.asarray([20.,30.,40.])
x = np.concatenate((x1.reshape(3,1),x2.reshape(3,1),x3.reshape(3,1)),axis=1)

y1 = np.asarray([4.0,88.0,35.0])
y2 = np.asarray([7.0,65.0,99.0])
y3 = np.asarray([40.0,13.0,27.0])
y = np.concatenate((y1.reshape(3,1),y2.reshape(3,1),y3.reshape(3,1)),axis=1)

def euclidean_distance(p1,p2): 
    return np.sqrt(np.sum((p1-p2)**2.0))

所以，我可以写：

distances = np.zeros((y.shape[0],x.shape[0]))
for i in range(y.shape[0]):
   for j in range(x.shape[0]):
       distances[i,j] = euclidean_distance(y[i,:],x[j,:])

这就是我接下来要排序的内容。在上面的当前 for 循环中，我将只选择我的 k 个最近邻居并在该内部循环中找到类......但它比在矢量化计算中计算所有距离要慢得多。

Answer 1

正如 Divakar 已经提到的，最简单的选择可能是scipy.spatial.distance.cdist：

from scipy.spatial.distance import cdist

distances = cdist(y, x)                 # Euclidean
distances = cdist(y, x, 'mahalanobis')  # Mahalanobis

这是单线程但速度很快。你也可以使用np.linalg.norm:

distances = np.linalg.norm(y[:, None, :] - x[None, :, :], axis=2)   # Euclidean

这会广播出对 x 和 y 中的行对的差异计算，以创建形状为 (3, 3, 3) 的中间数组，然后计算最后一个轴上的欧几里得范数。这是多线程的，但如果x 和y 有很多行（它也没有利用距离矩阵的对称性），则需要构造一个可能非常大的中间数组。

推广第二种计算马氏距离而不是欧几里得距离的方法是相当简单的（这部分我会留给你弄清楚......）。

Answer 2

我愿意：

• 重复 + 将两个数组重塑为 3D 形式 (3 x len(x) x len(y))
• 沿轴 = 0 取差、平方和和
• 现在您有了一个 2D 距离数组，并且可以沿相应的轴取最小值

这对你有帮助吗？或者我会尝试写得更明确...

对第 2 步的评论：您不必通过 sqrt 来找到最小值，您也可以将平方最小化

Answer 3

尝试广播采取交叉差异：

d = np.sqrt(np.sum((y[:,None,:]-x[None,:,:])**2,axis=-1))

我的测试脚本

import numpy as np

x1 = np.asarray([1.0,10.0,100.0])
x2 = np.asarray([40.0,60.0,80.0])
x3 = np.asarray([20.,30.,40.])
x = np.concatenate([i.reshape(-1,1) for i in [x1,x2,x3]], axis=1)
# see also column_stack

y1 = np.asarray([4.0,88.0,35.0])
y2 = np.asarray([7.0,65.0,99.0])
y3 = np.asarray([40.0,13.0,27.0])
"""
y1 = np.asarray([4.0,88.0])   # test 2d y
y2 = np.asarray([7.0,99.0])
y3 = np.asarray([13.0,27.0])
"""
y = np.concatenate([i.reshape(-1,1) for i in [y1,y2,y3]], axis=1)

def euclidean_distance(p1,p2):
    return np.sqrt(np.sum((p1-p2)**2.0))

distances = np.zeros((y.shape[0],x.shape[0]))
for i in range(y.shape[0]):
   for j in range(x.shape[0]):
       distances[i,j] = euclidean_distance(y[i,:],x[j,:])

print (distances)

d = np.sqrt(np.sum((y[:,None,:]-x[None,:,:])**2,axis=-1))
print(d)

生产

1230:~/mypy$ python2.7 stack35961972.py 
[[  38.70400496   54.2678542   120.60265337]
 [  90.79096871   79.98749902   33.13608305]
 [  68.45436436   46.42197755   68.95650803]]
[[  38.70400496   54.2678542   120.60265337]
 [  90.79096871   79.98749902   33.13608305]
 [  68.45436436   46.42197755   68.95650803]]