k近邻法的实现：kd树

参考文献：【KD树简介】，【统计学习方法】k近邻 kd树的python实现，【kd树求k近邻 python 代码】

1. 前言

实现kkk近邻法时，主要考虑的问题是如何对训练数据进行快速kkk近邻搜索。
kkk近邻法最简单的实现方法是线性扫描（linear scan）。这时要计算输入实例与每一个训练实例的距离。当训练数据容量很大时，计算非常耗时。
为了提高kkk近邻搜索的效率，可以考虑使用特殊的结构存储训练数据。以减少计算距离的次数。

2. 构造kd树

kdkdkd树是一种对kkk维空间中的实例点进行存储以便对其快读检索的树形数据结构。kdkdkd树是二叉树，表示对kkk维空间的一个划分（partition）。构造kdkdkd树相当于不断地用垂直于坐标轴的超平面将kkk维空间切分，构成一系列的kkk维超矩形区域。kdkdkd树的每个结点对应于一个kkk维超矩形区域。
方法：
构造根结点，使根结点对应于kkk维空间中包含所有实例点的超矩形区域；通过下面的递归方法，不断地对kkk维空间进行切分，生成子结点。在超矩形区域（结点）上选择一个坐标轴和在此坐标轴上的一个切分点，确定一个超平面，这个超平面通过选定的切分点并且垂直于选定的坐标轴，将当前超矩形区域切分为两个子区域（子结点）；这时，实例被分到两个子区域。这个过程直到子区域内没有实例时终止（终止时的结点为叶结点）。在此过程中，将实例保存在相应的结点上。

一般情况下，依次选择坐标轴对空间切分，选择训练实例点在选定坐标轴上的中位数为切分点，这样得到的kdkdkd树是平衡的。注意，平衡的kdkdkd树搜索时的效率未必是最优的。

3. 算法（构造平衡kd树）

输入：kkk维空间数据集T={x1,x2,...,xN}T=\{x_1,x_2,...,x_N\}T={x1,x2,...,xN}，其中xi=(xi(1),xi(2),...,xi(k))Tx_i=(x_i^{(1)},x_i^{(2)},...,x_i^{(k)})^Txi=(xi(1),xi(2),...,xi(k))T，i=1,2,...,Ni=1,2,...,Ni=1,2,...,N；
输出：kdkdkd树
第一步：构造根结点，根结点对应于包含TTT的kkk维空间的超矩形区域。
选择x(1)x^{(1)}x(1)为坐标轴，以TTT中所有实例的x(1)x^{(1)}x(1)坐标的中位数为切分点，将根结点对应的超矩形区域切分为两个子区域。切分由通过切分点并与坐标轴x(1)x^{(1)}x(1)垂直的超平面实现。
由根结点生成深度为1的左、右子结点对应坐标x(1)x^{(1)}x(1)小于切分点的子区域，右子结点对应于坐标x(1)x^{(1)}x(1)大于切分点的子区域。
将落在切分平面上的实例点保存在根结点中。
第二步：重复：对深度为jjj的结点，选择x(l)x^{(l)}x(l)为切分的坐标轴
l=j(modk)+1l=j \ (mod \ k) + 1 l=j (mod k)+1
以该结点的区域中所有实例的x(l)x^{(l)}x(l)坐标的中位数为切分点，将该结点对应的超矩形区域切分为两个子区域。切分由通过切分点并与坐标轴x(l)x^{(l)}x(l)垂直的超平面实现。
由该结点生成深度为j+1j+1j+1的左、右子结点：左子结点对应坐标x(l)x^{(l)}x(l)小于切分点的子区域，右子结点对应坐标x(l)x^{(l)}x(l)大于切分点的子区域。
将落在切分超平面上的实例点保存在该结点。
第三步：
直到两个子区域没有实例存在时停止。从而形成kdkdkd树的区域划分。

例题：给定一个二维空间的数据集：
T={(2,3)T,(5,4)T,(9,6)T,(4,7)T,(8,1)T,(7,2)T}T=\{(2,3)^T,(5,4)^T,(9,6)^T,(4,7)^T,(8,1)^T,(7,2)^T\} T={(2,3)T,(5,4)T,(9,6)T,(4,7)T,(8,1)T,(7,2)T}
构造一个平衡kdkdkd树。
解：根结点对应包含数据集TTT的矩形，选择x(1)x^{(1)}x(1)轴，6个数据点的x(1)x^{(1)}x(1)坐标的中位数是7，以平面x(1)=7x^{(1)}=7x(1)=7将空间分为左、右两个子矩形（子结点）；接着，左矩形以x(2)=4x^{(2)}=4x(2)=4分为两个子矩形，右矩形以x(2)=6x^{(2)}=6x(2)=6分为两个子矩形，如此递归，最后得到如下划分。

树形图：

4. 搜索kd树

给定一个目标点，搜索其最近邻。首先找到包含目标点的叶结点；然后从该叶结点出发，依次回退到父结点；不断查找与目标点最邻近的结点，当确定不可能存在更近的结点时终止。这样搜索就被限制在空间的局部区域上，效率大为提高。
包含目标点的叶结点对应包含目标点的最小超矩形区域。以此叶结点的实例点作为当前最近点。目标点的最近邻一定在以目标点为中心并通过当前最近点的超球体的内部。然后返回当前结点的父结点，如果父结点的另一子结点的超矩形区域与超球体相交，那么在相交的区域内寻找与目标点更近的实例点。如果存在这样的点，将此点作为新的当前最近点。算法转到更上一级的父结点，继续上述过程。如果父结点的另一子结点的超矩形区域与超球体不相交，或不存在比当前最近点更近的点，则停止搜索。

算法：（用kdkdkd树的最近邻搜索）
输入：已构造的kdkdkd树，目标点xxx；
输出：xxx的最近邻
第一步：在kdkdkd树中找出包含目标点xxx的叶结点：从根结点出发，递归地向下访问kdkdkd树。若目标点xxx当前维的坐标小于切分点的坐标，则移动到左子结点，否则移动到右子结点。直到子结点为叶结点为止。
第二步：以此叶结点为“当前最近点”。
第三步：递归地向上回退，在每个结点进行以下操作：

如果该结点保存的实例点比当前最近点距离目标点更近，则以该实例点为“当前最近点”。
当前最近点一定存在于该结点一个子结点对应的区域。检查该子结点的父结点的另一子结点对应的区域是否有更近的点。具体地，检查另一结点对应的区域是否与==以目标点为球心、以目标点与“当前最近点”间的距离为半径的超球体相交。、
–如果相交，可能在另一子结点对应的区域内存在距目标点更近的点，移动到另一个子结点。接着，递归地进行最近邻搜索；
–如果不相交，向上回退。

第四步：当回退到根结点时，搜索结束。最后的”当前最近点“即为xxx的最近邻点。

如果实例点是随机分布的，kdkdkd树搜索的平均计算复杂度是O(log⁡N)O(\log \ N)O(log N)，这里NNN是训练实例数。kdkdkd树更适用于训练实例数远大于空间维数时的kkk近邻搜索。当空间维数接近训练实例数时，它的效率会迅速下降，几乎接近线性扫描。

5. 完整代码

5.1 最近邻

import numpy as np
class Node:  # 结点def __init__(self, data, lchild=None, rchild=None):self.data = dataself.lchild = lchildself.rchild = rchildclass KdTree:  # kd树def __init__(self):self.kdTree = Nonedef create(self, dataSet, depth):  # 创建kd树，返回根结点  参数1：数据集   参数2：树的深度（同时决定：排序依据的维度）if (len(dataSet) > 0):m, n = np.shape(dataSet)  # 求出样本行，列midIndex = int(m / 2)  # 中间数的索引位置axis = depth % n  # 判断以哪个轴划分数据sortedDataSet = self.sort(dataSet, axis)  # 进行排序node = Node(sortedDataSet[midIndex])  # 将节点数据域设置为中位数，具体参考下书本# 将两边的数切割成两个集合，然后递归调用create方法leftDataSet = sortedDataSet[: midIndex]  # 将中位数的左边创建2改副本rightDataSet = sortedDataSet[midIndex + 1:]print(leftDataSet)print(rightDataSet)node.lchild = self.create(leftDataSet, depth + 1)  # 将中位数左边样本传入来递归创建树node.rchild = self.create(rightDataSet, depth + 1)return nodeelse:return Nonedef sort(self, dataSet, axis):  # 采用冒泡排序，利用aixs作为轴进行划分sortDataSet = dataSet[:]  # 由于不能破坏原样本，此处建立一个副本m, n = np.shape(sortDataSet)for i in range(m):for j in range(0, m - i - 1):if (sortDataSet[j][axis] > sortDataSet[j + 1][axis]):  #把数值大的一个一个向下沉temp = sortDataSet[j]sortDataSet[j] = sortDataSet[j + 1]sortDataSet[j + 1] = tempprint(sortDataSet)return sortDataSedef preOrder(self, node):  # 前序遍历if node != None:print("tttt->%s" % node.data)self.preOrder(node.lchild)self.preOrder(node.rchild)def search(self, tree, x):  # 搜索self.nearestPoint = None  # 保存最近的点self.nearestValue = 0  # 保存最近的值def travel(node, depth=0):  # 递归搜索if node != None:  # 递归终止条件n = len(x)      # 特征数axis = depth % n  # 计算轴if x[axis] < node.data[axis]:  # 如果数据小于结点，则往左结点找travel(node.lchild, depth + 1)else:travel(node.rchild, depth + 1)# 以下是递归完毕后，往父结点方向回朔，对应算法3.3(3)distNodeAndX = self.dist(x, node.data)  # 目标和节点的距离判断if (self.nearestPoint == None):  # 确定当前点，更新最近的点和最近的值，对应算法3.3(3)(a)self.nearestPoint = node.dataself.nearestValue = distNodeAndXelif (self.nearestValue > distNodeAndX):self.nearestPoint = node.dataself.nearestValue = distNodeAndXprint(node.data,'\t', depth,'\t',self.nearestValue,distNodeAndX,'\t',node.data[axis], x[axis])if (abs(x[axis] - node.data[axis]) <= self.nearestValue):  # 确定是否需要去子节点的区域去找（圆的判断），对应算法3.3(3)(b)if x[axis] < node.data[axis]:travel(node.rchild, depth + 1)else:travel(node.lchild, depth + 1)print('位置\t树深度\t当前最近距离\t当前点与目标点距离\t')travel(tree)        #tree的属性与根结点是一样的，所以可以带入定义方法参数时，里面的node属性return self.nearestPointdef dist(self, x1, x2):  # 欧式距离的计算return ((np.array(x1) - np.array(x2)) ** 2).sum() ** 0.5if __name__ == '__main__':dataSet = [[2, 3],[5, 4],[9, 6],[4, 7],[8, 1],[7, 2]]x = [3, 6]kdtree = KdTree()tree = kdtree.create(dataSet, 0)   #创建KD树print('打印先根遍历的输出顺序')kdtree.preOrder(tree)              #先序遍历测试print("给定的数据的最近邻",kdtree.search(tree, x))      #在树中搜索给定的 X的最近邻

5.2 k近邻

1 import numpy as np2 from numpy import array3 class decisionnode:4     def __init__(self,value=None,col=None,rb=None,lb=None):5         self.value=value6         self.col=col7         self.rb=rb8         self.lb=lb9
10 #读取数据并将数据转换为矩阵形式
11 def readdata(filename):
12     data=open(filename).readlines()
13     x=[]
14     for line in data:
15         line=line.strip().split('\t')
16         x_i=[]
17         for num in line:
18             num=float(num)
19             x_i.append(num)
20         x.append(x_i)
21     x=array(x)
22     return x
23
24 #求序列的中值
25 def median(x):
26     n=len(x)
27     x=list(x)
28     x_order=sorted(x)
29     return x_order[n//2],x.index(x_order[n//2])
30
31 #以j列的中值划分数据，左小右大，j=节点深度%列数
32 def buildtree(x,j=0):
33     rb=[]
34     lb=[]
35     m,n=x.shape
36     if m==0: return None
37     edge,row=median(x[:,j].copy())
38     for i in range(m):
39         if x[i][j]>edge:
40             rb.append(i)
41         if x[i][j]<edge:
42             lb.append(i)
43     rb_x=x[rb,:]
44     lb_x=x[lb,:]
45     rightBranch=buildtree(rb_x,(j+1)%n)
46     leftBranch=buildtree(lb_x,(j+1)%n)
47     return decisionnode(x[row,:],j,rightBranch,leftBranch)
48
49 #搜索树：输出目标点的近邻点
50 def traveltree(node,aim):
51     global pointlist  #存储排序后的k近邻点和对应距离
52     if node==None: return
53     col=node.col
54     if aim[col]>node.value[col]:
55         traveltree(node.rb,aim)
56     if aim[col]<node.value[col]:
57         traveltree(node.lb,aim)
58     dis=dist(node.value,aim)
59     if len(knears)<k:
60         knears.setdefault(tuple(node.value.tolist()),dis)#列表不能作为字典的键
61         pointlist=sorted(knears.items(),key=lambda item: item[1],reverse=True)
62     elif dis<=pointlist[0][1]:
63         knears.setdefault(tuple(node.value.tolist()),dis)
64         pointlist=sorted(knears.items(),key=lambda item: item[1],reverse=True)
65     if node.rb!=None or node.lb!=None:
66         if abs(aim[node.col] - node.value[node.col]) < pointlist[0][1]:
67             if aim[node.col]<node.value[node.col]:
68                 traveltree(node.rb,aim)
69             if aim[node.col]>node.value[node.col]:
70                 traveltree(node.lb,aim)
71     return pointlist
72
73 def dist(x1, x2): #欧式距离的计算
74     return ((np.array(x1) - np.array(x2)) ** 2).sum() ** 0.5
75
76 knears={}
77 k=int(input('请输入k的值'))
78 if k<2: print('k不能是1')
79 global pointlist
80 pointlist=[]
81 file=input('请输入数据文件地址')
82 data=readdata(file)
83 tree=buildtree(data)
84 tmp=input('请输入目标点')
85 tmp=tmp.split(',')
86 aim=[]
87 for num in tmp:
88     num=float(num)
89     aim.append(num)
90 aim=tuple(aim)
91 pointlist=traveltree(tree,aim)
92 for point in pointlist[-k:]:
93     print(point)