代码详解：用Python构建邻域矩阵

全文共1522字，预计学习时长3分钟

如何构建邻域矩阵？只掌握理论知识，很难让人深入地了解这一问题。

本文我们以这样一个模型为例：

\tilde = W \cdot YY~=W⋅Y

其中W是一个矩阵，可以定义为（例如）：

w_ = 1wij=1

当j在最临近i的K内，Wij=1，否则结果为0。

下面，我们将详细介绍通过使用numpy、scipy和matplotlib进行可视化来创建W矩阵的方法。

样本数据

为达到演示目的，我们创建了一个虚拟数据集，大小为N = 12个训练样本，M = 3个测试样本：

import numpy as np

XY_test = np.array([[1, 1], [-1, 1], [-1, -1], [1, -1]])

让我们来看看这些点的分布情况：红点是训练数据，而绿点是测试数据。

寻找邻域

用现代工具寻找邻域是非常简单的。在这里，我们选择使用scipy，因为稍后将使用这个软件包中的其他工具，但sklearn或其他软件包也可以完成这项工作。使用scipy时，首先使用训练数据集创建一个cKDTree：

from scipy.spatial import cKDTree

tree = cKDTree(XY_train)

可再次查询该tree：

K = 3

result = tree.query(XY_test, k=K)

在这里，我们需要测试样本元素中训练样本的三个最近邻域。默认情况下，tree.query返回邻域索引和相关距离。我们将保留两者。

distances, indices = result

让我们集中讨论索引数组。

array([[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]])

索引数组是一个numpy数组，包含M（测试样本数）行和K（邻域数）列。那么如何将其转换为我们需要的矩阵呢？示例如下：

在一个情节中看到所选邻域会很有趣：

import matplotlib.pyplot as plt

n = 0 # first element in the test dataset

xy_test = XY_test[n]

index = indices[n]

neighbours = XY_train[index]

plt.clf()

plt.scatter(xy_test[0], xy_test[1], color="red")

plt.scatter(neighbours[:,0], neighbours[:,1], color="blue")

plt.xlabel("x")

plt.ylabel("y")

plt.xlim(-2, 2)

plt.ylim(-2, 2)

plt.show()

好吧，所以邻域的寻找似乎和预期的一样有效！让我们了解如何将索引数组转换为完全可用的矩阵，我们的目的应该是：

1 1 1 0 0 0 0 0 0 0 0 0

0 0 0 1 1 1 0 0 0 0 0 0

0 0 0 0 0 0 1 1 1 0 0 0

0 0 0 0 0 0 0 0 0 1 1 1

因为测试观察0（第一行）的邻域是训练观察0,1和2，所以测试观察1（第二行）的邻域是训练观测3,4和5，依此类推。

创建矩阵

首先，我们将使用numpy索引创建这样的矩阵：

import numpy as np

a = np.array([1, 2, 3, 4, 5, 6])

i = [0, 0, 1, 1, 2, 2]

a[i]

# array([1, 1, 2, 2, 3, 3])

但你会发现它不适用于多维度数组。

我们选择的解决方案是使用scipy稀疏矩阵，其可在索引列表中创建。例如，使用稀疏矩阵创建大小为N = 4的对角矩阵，可以表示为：

from scipy import sparse

i_index = [0, 1, 2, 3]

j_index = [0, 1, 2, 3]

values = [1, 1, 1, 1]

matrix = sparse.coo_matrix((values, (i_index, j_index)), shape=(4, 4))

print(matrix)

# (0, 0)1

# (1, 1)1

# (2, 2)1

# (3, 3)1

因此scipy获取i_index和j_index数组的第一个元素i和j，并将values数组的第一个元素放在最终矩阵中的[i，j]位置。或者换句话说，元素（0,0）的值是1，元素（1,1）的值也是1 ...... 其他未特殊说明的元素的值都是零。

如果你更喜欢数组表示，可以输入以下代码查看结果：

matrix.toarray() # transforms sparse matrix into numpy array just for visualization

#array([[1, 0, 0, 0],

# [0, 1, 0, 0],

# [0, 0, 1, 0],

# [0, 0, 0, 1]])

这里你可以看到对角线矩阵。

让我们用第二个例子来更清楚地说明一切。现在要创建的是逆对角矩阵：

array([[0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]])

此次的编码是：

i_index = [3, 2, 1, 0] #

j_index = [0, 1, 2, 3]

values = [1, 1, 1, 1]

matrix = sparse.coo_matrix((values, (i_index, j_index)), shape=(4, 4))

注意：只有当矩阵规模相对较小时才能从稀疏表示切换到密集表示，否则会出现内存问题（稀疏矩阵存在的原因！）

如何创建W矩阵？

对于w矩阵，j_index（即“列”）对应相邻索引：

j_index = indices.flatten()

#array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])

然后行索引i_index 对应测试样本中的索引，但需要重复K次以匹配j_index 排序：

i_index = np.repeat(np.array(range(M), dtype=int), repeats=K, axis=0).ravel()

#array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3])

这意味着在第一行（行索引0）的0、1和2列将有一个索引。在第二行（1）的第3、4、5列中会有一个索引……如果你再次查看测试/训练样本位置（第一个图），结果是一致的！

我们假设所有值为“1”：

values = np.ones(M * K) # M = number of test sample, K = number of neighbours

或者用取决于距离的函数表示，例如：

values = 1. / distances.flatten()**2

最后，我们的矩阵看起来像（值为“1”）：

matrix = sparse.coo_matrix((values, (i_index, j_index)), shape=(M, N))

# array([[1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],

# [0., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0.],

# [0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0.],

# [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1.]])

回到我们最初的问题

现在，我们可以计算点积（使用稀疏或密集矩阵）：

y_tilde = matrix.dot(y) # where y has shape (N, )

最后，问题解决啦！

留言 点赞 发个朋友圈

我们一起分享AI学习与发展的干货

编译组：王玲、狄思云