K-D树学习笔记

K-dimension Tree，简称K-D树，适用于在$k$维欧几里得空间中维护信息的数据结构。
一般的数据结构只能维护线性或低维空间，而K-D树可以适用于$k$维空间，只是效率有所差异。
K-D树是从排序二叉树扩展而来，其本质还是一棵排序二叉树。

定义

比较维度

定义，一个结点的比较维度$d$指当前结点在$d$维上满足二叉树性质。
如，第一层的$x$维度满足排序二叉树性质，第二层的$y$维度满足排序二叉树性质。

分辨器

K-D树是一棵排序二叉树，它不能直接地维护$k$维空间，只能分层分别维护$k$维空间。
分辨器用于判断当前结点所比较的维度，第$i$层的分辨器所比较的维度是$i\bmod k$。
举个栗子：如图，第一层结点比较$x$维，第二层结点比较$y$维，第三层结点比较$x$维。

因此，我们可以将K-D树上的每一个非叶子节点视作一个超平面，按对应维度将空间分为两个半空间。
上图对应的2-d树，在平面上长这个样子：

不完全平衡性

和排序二叉树类似的，K-D树不是完全平衡的。
在某些情况下K-D树是平衡的，而某些情况K-D树可能不平衡。
在接下来的操作中我会说明其对平衡性的影响。

操作

构建

根据对分辨器的定义，K-D树的结构框架基本清晰了，接下来叙述如何构建一棵K-D树。

根据分辨器，选出对应维度靠中间的点$mid$。
将所有对应维度比$mid$小的点放到左子树递归处理，比$mid$大的点放到右子树递归处理。
这一步我们可以使用STL中的nth_element来完成。
nth_element(start,nth,end)表示按照小于符号的比较方式，将开始指针为start，结束指针后一位为end的区间，使得nth位置即为所有元素第$n$小的值，并将不大于$nth$的值放在nth位置前，不小于nth的放在nth位置后，但不一定有序，时间复杂度$O(end-start)$，完全满足需求。

一直递归直到仅有一个元素为止。
不难发现，构建一棵K-D树的时间复杂度是$O(n\log n)$的。
此时这棵K-D树是十分平衡的，叶子结点高度十分接近。

int D;

struct Point {
	int d[K];
	int ls,rs,val;
	int& operator [] (int x) {
		return d[x];
	}
	bool operator < (const Point& b) const {
		return d[D]<b.d[D];
	}
} p[maxn];

int build(int Left,int Right,int now) {
	int mid=(Left+Right)>>1,root=mid;
	D=now;
	nth_element(p+Left,p+mid,p+Right+1);
	if(Left<mid)p[root].ls=build(Left,mid-1,(now+1)%K);
	if(Right>mid)p[root].rs=build(mid+1,Right,(now+1)%K);
	return root;
}

插入

与二叉树的插入类似，根据分辨器决定比较维度，然后使用排序二叉树的方法在左子树或右子树插入。

#define lson p[index].ls
#define rson p[index].rs

int insert(int index,Point P,int now) {
	if(!index) {
		p[++size]=P;
		return size;
	}
	if(P[now]<=p[index][now])lson=insert(lson,P,(now+1)%K);
	else rson=insert(rson,P,(now+1)%K);
	return index;
}

在K-D树是平衡的时候，时间复杂度为$O(\log n)$。
在插入结点后，可能会影响K-D树的平衡性，导致K-D树失衡。

删除

删除？怎么做啊？不会啊？
这很简单嘛，打个标记就完了。

查询

查询精确点

和二叉查找树的查询类似，根据分辨器决定对应维度，根据排序二叉树性质选择子树递归。

bool check(Point x,Point y) {
	for(int i=0; i<K; i++)
		if(x[i]!=y[i])return false;
	return true;
}

int find(int index,Point P,int now) {
	if(check(p[index],P))return index;
	if(P[now]<=p[index][now])return find(lson,P,(now+1)%K);
	else return find(rson,P,(now+1)%K);
}

这是一个鸡肋的查询操作，基本上不会用到。
在K-D树是平衡的时候，时间复杂度为$O(\log n)$。

模糊查询

查询最远点

使用K-D树查询距离$P$最远点实际上是一个搜索剪枝的过程。

首先，我们要在构建与插入时维护两个信息：$Min[]$与$Max[]$，表示对应的维度的最小值与最大值。

定义乐观估价函数$get$-$max(index,P)$表示在$index$控制范围中包含的点与$P$最远距离的上界。（不一定能达到，但一定不大于它）
计算方法是对每个维度取最大差值计算距离。

LL sqr(LL x) {
	return x*x;
}

LL get_max(int index,Point P) { 
	if(!index)return 0;
	LL ans=0;
	for(int i=0; i<K; i++)ans+=max(sqr(p[index].Max[i]-P[i]),sqr(p[index].Min[i]-P[i]));
	return ans;
}

• 初始化$ans=+\infty$。
• 从根结点开始，先用根节点代表的点更新答案。左右儿子均代表一部分半空间，都可能更新答案。计算左右儿子的乐观估价$ldist,rdist$，若$ldist\gt rdist$，先递归左儿子，再递归右儿子，反之类似。在递归前需要判断乐观估价函数是否$\gt ans$，若不成立显然没有递归的意义。

LL dist(Point a,Point b) { 
	LL ans=0;
	for(int i=0; i<K; i++)ans+=sqr(a[i]-b[i]);
	return ans;
}

void find_max(int index,Point P) {
	if(!index)return;
	LL Dist=dist(p[index],P);
	ans=max(ans,Dist);
	LL ldist=get_max(lson,P),rdist=get_max(rson,P);
	if(ldist>rdist) {
		if(ldist>ans)find_max(lson,P);
		if(rdist>ans)find_max(rson,P);
	} else {
		if(rdist>ans)find_max(rson,P);
		if(ldist>ans)find_max(lson,P);
	}
}

查询最近点

使用K-D树查询距离$P$最近点也是一个搜索剪枝的过程。

同样维护出$Min[],Max[]$。

因为查询最近点，因此不能再使用上述乐观估价函数。
我们重新定义乐观估价函数。
若一个矩阵在某一维度上的最小值比$P$的该维度还要大，或最大值比它还小才能更新估价。

LL get_min(int index,Point P) {
	if(!index)return LLONG_MAX/2;
	LL ans=0;
	for(int i=0; i<K; i++) {
		if(p[index].Min[i]-P[i]>0)ans+=sqr(p[index].Min[i]-P[i]);
		if(P[i]-p[index].Max[i]>0)ans+=sqr(P[i]-p[index].Max[i]);
	}
	return ans;
}

除了乐观估价函数，最近点还可以使用边界相交判断。

先递归其中一个儿子，若以目标点为圆心，以答案为半径作圆，不与当前边界相交，则无递归另一个儿子的必要。

接下来寻找过程与最大值类似。

void find_min(int index,Point P,int now) {
	if(!index)return;
	LL Dist=dist(p[index],P);
	ans=min(ans,Dist);
	LL ldist=get_min(lson,P),rdist=get_min(rson,P);
	if(ldist<rdist) {
		if(ldist<ans)find_min(lson,P,(now+1)%K);
		if(rdist<ans&&P[now]+ans>=p[index][now])find_min(rson,P,(now+1)%K);
	} else {
		if(rdist<ans)find_min(rson,P,(now+1)%K);
		if(ldist<ans&&P[now]-ans<=p[index][now])find_min(lson,P,(now+1)%K);
	}
}

查询$k$远/近点

将询问过程中的更新答案部分使用堆维护。
如在查询第$k$远点时，开一个小根堆（初始化$k$个$0$），当$Dist$比堆顶大时，弹出堆顶，$Dist$入堆，最后堆顶即为第$k$远点。

void find_max(int index,Point P) {
	if(!index)return;
	LL Dist=dist(p[index],P);
	if(Dist>Q.top()) {
		Q.pop();
		Q.push(Dist);
	}
	LL ldist=get_max(lson,P),rdist=get_max(rson,P);
	if(ldist>rdist) {
		if(ldist>Q.top())find_max(lson,P);
		if(rdist>Q.top())find_max(rson,P);
	} else {
		if(rdist>Q.top())find_max(rson,P);
		if(ldist>Q.top())find_max(lson,P);
	}
}

第$k$近点与之类似。
之前提到的所有剪枝均适用。

范围查询信息

除了查询点对距离问题，K-D树还能用于定位一部分高维空间，可以相当于高维空间的平衡树维护信息。

如查询一部分$k$维空间中的最大权值。
使用$Left[],Right[]$表示询问的$k$维空间边界。
使用$p[index].max$记录$index$代表的空间中的最大权值。

bool inspace(Point p) { //判断是否完全覆盖
	for(int i=0; i<K; i++)
		if(p[i]<Left[i]||p[i]>Right[i])return false;
	return true;
}
bool check(int index) { //判断是否相交
	if(!index)return 0;
	for(int i=0; i<K; i++)
		if(p[index].Max[i]<Left[i]||p[index].Min[i]>Right[i])return false;
	return true;
}
void query(int index) {
	if(!index)return;
	if(inspace(p[index]))ans=max(ans,p[index].val);
	bool lable=check(lson),rable=check(rson);
	if(p[lson].max>p[rson].max) {
		if(lable&&p[lson].max>ans)query(lson);
		if(rable&&p[rson].max>ans)query(rson);
	} else {
		if(rable&&p[rson].max>ans)query(rson);
		if(lable&&p[lson].max>ans)query(lson);
	}
}

模板

见这里

时间复杂度

在K-D树是平衡的时候，$n$个结点的$k$维K-D树模糊查询搜索过程时间复杂度上界为：

$$O(kn^{1-\frac 1k})$$

在$k=2$时，时间复杂度为$O(\sqrt n)$，在$k$很大的时候，和暴力几乎无区别。

优化

我们发现K-D树在多次插入/删除后失衡，导致时间复杂度严重退化。
对于本问题较好的解决方案是使得K-D树平衡化。
因为K-D树不支持旋转操作，因此加入非旋转机制是较好的选择。

• 定期重构
  每经过$O(\sqrt q)$的时间就对K-D树重新build一次，可以有效防止退化为$O(n)$的复杂度。
• 替罪羊化
  除了定期重构，还可以对K-D树进行替罪羊化，发现依然满足优美的均摊性质，因此在不计算模糊查询的复杂度时，可以保持$O(\log n)$的优秀复杂度。

参考资料

• k-d tree - Wikepedia
• k-d树.ppt - 巴蜀中学讲义
• k-d树学习笔记 - sengxian
• 从K近邻算法谈到KD树、SIFT+BBF算法 - _harvey