Yin's Blog
0%

使用Numpy计算距离矩阵[翻译]

Posted on Edited on

这是一篇通过NumPy计算距离矩阵的文章，原作者是Jay Mody。文章介绍了距离矩阵的概念，并且通过层层递进的方式，特别清晰地讲解了使用NumPy进行向量计算的思路和技巧。

原文链接：Computing Distance Matrices with NumPy

超级推荐：）

(已获原文作者授权)

背景

一般来说，距离矩阵是一个方阵，它刻画了一组向量之间的成对距离。更形式化的定义如下：

给定一组向量$v_1, v_2, … v_n$和它的距离矩阵dist，矩阵元素$dist_{ij}$表示$v_i$和$v_j$之间的距离。注意：考虑到$dist_{ij} = dist_{ji}$，这意味着这个矩阵是对称的，并且矩阵的维度是$(n, n)$。

但是以上距离矩阵定义还是没有给出距离的定义。事实上，我们有很多的方式去定义和计算两个向量间的距离，但是通常而言，当我们提到向量间的距离时，我们指的是它们的欧几里得距离(欧氏距离)。 欧氏距离其实就是我们直觉上认为的距离（即：图上两点间的最短线段）。数学上，我们可以定义两向量$u$, $v$间的欧氏距离为：

$$||u - v||_2 = \sqrt{\sum_{k=1}^d (u_k - v_k)^2}$$

其中$d$是向量的维度。

距离矩阵在各类场景下被广泛应用，包括数学、计算机科学、图论、生物信息等领域。本文就来讨论一下距离矩阵在机器学习领域的应用。

引入的例子：k最近邻算法

k最近邻算法（kNN）是一类机器学习分类算法（也可以做回归），其底层计算充分利用了距离矩阵。它的想法很简单，我们可以通过查看其$k$个最近邻的标记数据点的类别来预测任何给定数据点的类别。

那我们如何确定哪个标注数据点是最近邻的呢？如果我们将每个数据点表示成向量的话，我们就可以计算他们的欧氏距离了。（有了距离，当然也就能判断近与远啦）

假如在这里，我们不仅仅想对单一数据点进行预测，而是对多个数据点进行预测。更形式化的说法是，给定$n$个标注数据点（训练数据），$m$个未标注数据点（测试数据，也即我们想要去分类的数据点）。这些数据点可以被表示为$d$维向量。为了实现kNN分类器，你需要计算所有标注-未标注数据对之间的距离。这些距离将被存储在一个维度为$(m，n)$的矩阵dist中，其中$dist_{ij}$表示第$i$个未标注数据点和第$j$个标注数据点之间的距离。如果我们将标注数据点表示为$(n, d)$维的矩阵$Y$，而未标注的的数据点表示为$(m, d)$维的矩阵$X$，那么距离矩阵就可以表示为：

$$dist_{ij} = \sqrt{\sum_{k=1}^d (X_{ik} - Y_{jk})^2}$$

这个距离计算便是算法的核心，也是本文的专注所在。那么下面我们就来实现它。

注意: 尽管矩阵不再是方阵，我在这里还是使用距离矩阵这个术语。（因为我们在这里是计算两组向量间的距离，而不是一组。）

使用三组循环

实现距离矩阵计算的最简单方法是对所有的对和元素遍历循环：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
X # 测试数据 (m, d)
X_train # 训练数据 (n, d)

m = X.shape[0]
n = X_train.shape[0]
d = X.shape[1]
dists = np.zeros((num_test, num_train)) # 距离矩阵 (m, n)

for i in range(m):
    for j in range(n):
        val = 0
        for k in range(d):
            val += (X[i][k] - X_train[j][k]) ** 2
        dist[i][j] = np.sqrt(val)

虽然上述代码是可行的，但是它却非常的低效，完全没有利用numpy的高效向量运算。那么让我们来改进一下：

使用两组循环

1
2
3
4
5

for i in range(m):
    for j in range(n):
        # 逐元素相减，逐元素平方操作，求和并计算平方根）
        dist[i][j] = np.sqrt(np.sum((X[i] - X_train[j]) ** 2))

这样是不是好点啦？而且代码也变得简洁易读了很多。但是我们还可以进一步改进。

使用单组循环

1
2
for i in range(m):
    dist[i, :] = np.sqrt(np.sum((X[i] - X_train) ** 2, axis=1))

哇，这里发生了什么？好吧，让我们来拆解分析一下。

首先，X[i] - X_train这句不应该会报错的吗？X[i]的维度是$(d)$，然而X_train的形状却是$(n,d)$啊。逐元素操作只有在两部分拥有相同维度的情况下才可以啊，这是啥情况？

答案是：Numpy自动对X[i]进行了广播操作，从而使其匹配了X_train的维度。你可以认为是X[i]堆叠了$n$次，进而产生了一个维度为$(n,d)$的矩阵。其中该矩阵的每一行都是X[i]的一个拷贝。这样的话，当进行减法操作的时候，X_train的每一行就可以和X[i]进行相减了（这里的减与被减没所谓的，反正我们后续要对结果取平方）。如果你想的话，在numpy中，你甚至可以通过np.tile自行创建”堆叠”的矩阵。不过这样操作会比让numpy自动通过广播方式进行处理慢。现在我们就获得了一个$(n,d)$维度的矩阵，其中每行是X[i] - X_train[j]。

那么下一步就容易很多啦，我们进行逐元素的平方操作。然后我们需要对每一行求和，因此我们使用带着参数axis=1的np.sum操作。该操作让numpy通过第一轴向进行求和（即通过行）。注意，如果不带上轴向参数的话，np.sum将会对矩阵中的每个元素进行求和，并输出一个标量值。而通过带着参数axis=1的np.sum操作，我们将获得一个大小为$n$的向量。

最后，我们对该向量进行逐元素的平方根操作，并将其存储到dist[i]中。

好了，下面是一个更容易理解的注释版本：

1
2
3
4
5
6
7
8
9
10
11
12
13
for i in range(m):
    # X[i] 被广播化 (d) -> (n, d)
    # (每行都是X[i]的一份拷贝)
    diffs = X[i] - X_train

    # 逐元素平方操作
    squared = diffs ** 2

    # 对每行进行求和 (n, d) -> (n)
    sums = np.sum(squared, axis=1)

    # 逐元素的平方根操作，并存储到dist中
    dist[i, :] = np.sqrt(sums)

无循环版本？！

其实我们还可以做的更好，只使用向量/矩阵操作，连循环都不用。这是啥情况？我们再仔细看看我们的方程：

$$dist_{ij} = \sqrt{\sum_{k=1}^d (x_{ik} - y_{jk})^2}$$

想想看，如果我们把求和号里的表达式展开呢？

$$dist_{ij} = \sqrt{\sum_{k=1}^d x^2_{ik} - 2x_{ik}y_{jk} + y^2_{jk}}$$

有趣吧？让我们把这个求和分配开：

$$dist_{ij} = \sqrt{\sum_{k=1}^d x^2_{ik} - 2 \sum_{k=1}^d x_{ik}y_{jk} + \sum_{k=1}^dy^2_{jk}}$$

你会发现，上式的每一个求和其实都是点乘。现在我们就替换掉原来的符号，使用更干净的表达：

$$dist_{ij} = \sqrt{x_i \cdot x_i - 2x_i \cdot y_j + y_j \cdot y_j}$$

我们注意到，对于所有$i,j$的组合，中间项是唯一的，而左右两项都是重复的。考虑固定$i$或者$j$并对其它变量进行迭代，你会发现$x_i \cdot x_i$出现$j$次而$y_j \cdot y_j$出现$i$次。因此，这里的挑战是搞清楚如何计算所有可能的$x_i \cdot x_i$，$x_i \cdot y_j$和$y_j \cdot y_j$，然后将它们正确地加和。所有这些都不涉及循环。我们来试一下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# this has the same affect as taking the dot product of each row with itself
x2 = np.sum(X**2, axis=1) # shape of (m)
y2 = np.sum(X_train**2, axis=1) # shape of (n)

# we can compute all x_i * y_j and store it in a matrix at xy[i][j] by
# taking the matrix multiplication between X and X_train transpose
# if you're stuggling to understand this, draw out the matrices and
# do the matrix multiplication by hand
# (m, d) x (d, n) -> (m, n)
xy = np.matmul(X, X_train.T)

# each row in xy needs to be added with x2[i]
# each column of xy needs to be added with y2[j]
# to get everything to play well, we'll need to reshape
# x2 from (m) -> (m, 1), numpy will handle the rest of the broadcasting for us
# see: https://numpy.org/doc/stable/user/basics.broadcasting.html
x2 = x2.reshape(-1, 1)
dists = np.sqrt(x2 - 2*xy + y2) # (m, 1) repeat columnwise + (m, n) + (n) repeat rowwise -> (m, n)

-1组循环！！？！

哈哈哈哈哈

1
from sklearn.neighbors import KNeighborsClassifier

速度对比

为了比较每种实现的速度，我们可以使用cs231n课程作业1-knn中的cifar-10数据集的子集进行测试。

1
2
3
Two loop version took 39.707250 seconds
One loop version took 28.705156 seconds
No loop version took 0.218127 seconds

很明显，无循环版本是最后的胜利者，其以数量级的优势击败两组循环和一组循环版本。注意噢，我们在这里并没有使用三组循环的版本，因为这会跑好几个小时啊！在10个训练和10个测试例子的情况下，三组循环版本都花了0.5秒。作为参考，上述测试是基于5000个训练数据和500个测试数据的结果。好耶！给向量操作点赞 ：）