KD树-预排序加速以及动态化分割

背景

kd树可认为是二叉树的高维衍生, 如果二叉树是对一维的数据二分来建立, 那KD树就是对$K$维的点集不断切换维度二分. 比如一个二维点集有x,y两个维度, 则在建立树的第一层节点时, 便先对x轴分割, 到了第二层, 再对Y轴分割. 如此循环.

分割时需要用到取中位数算法(假设中位数是第(len-1)/2大的数), 有两种方案:

1, 排序, 并取中间值. ( 平均时间开销: $O(nlogn)$)
2, 通过topK查找第(len-1)/2大的数. ( 平均时间开销: $O(n)$)

『1』的平均开销虽大, 但是它能返回有序序列增加划分质量. 且可以通过一次预排序, 让后来每次分割都加速到稳定的$O(n)$. 但这部分处理, 需要一定的问题意识. 因此, 这个周末研究了一下如何预排序, 来达成这个的目的. 以下是遇到的两个问题以及解决:

---

两个问题

1. 本层的分割输入, 依赖上层的分割输出.——如下图, 红线在x方向上进行了一次分割后, 划分结果{E, D}和{A, B, C}两个集合, 对于下次划分来说是输入, 但它依赖上次划分…
   
2. 更新麻烦. 每当你插入或者移动一个新点时, 在树上行进的路径遇到的点, 你打算更新还是不更新?数组的新建需要增加开销, 何况会破坏有序性, 还引来代码膨胀.

问题一的本质在于, 每次的分割虽然是在单一维度上进行的, 但本质都是用对"点"的分割, 而非对某个坐标孤立的分割. 因此, 每次对坐标去分割时, 需要考虑这个操作映射到其他维度的变化情况.
如下例子(满秩)中, 假设点集有x和y两个维度(为了方便扩展到高维, 可以理解成dim0维和dim1维).

• 在两个方向上分别排序, 得到的数组分别为:
  x| {A, B, C, D, E, F, G, H, I, J, K ,L};
y | {A, F, B,I, G, D, H, K, C, J, L, E };
  
  之后针对x方向进行了一次分割, 得到结果(本例中只抓小于中位数的左集合):
  x| {A, B, C, D, E, -, - , - ,- ,-,- };
y | {A, F, B, -, -, D, -, -, C, -, -, E };
  如下所示(红线是分割的后的LeftChild点集). 可以看到, 对x方向取中位数的同时, 等于对y方向上用曲线进行了一次分割

x维的分割线(红色虚线)映射到y维上.

• 第2轮分割: 在这一轮中, 我们在水平方向取一次中点. 因此, 是对上一轮有序分割后的结果先拼接, 再取中点:
  下图是本轮的输入. y | {A F B D C E}

右下: 拼接后的元素.

结果如下:

如此不断分割、拼接、拼接, 便完成最终的划分. 整个过程刻画如下:
研究每一轮的横向的排列, 可以看到一组取中点对应一组过滤.

我们可以做如下归纳:

• 每个维度的有序序列, 都是同一组点的不同维度视角而已.
• 本维度序列"按某维度中位数分割"的同时还需要映射到其他的维度上
• "中位数分割"这项操作的本质, 就是用一个正交于自己维度的超平面去分类.
• 这个超平面在其他维度的视角下, 是平行于自己维度的,

总之, 我们每次二分, 都要同时对所有序列上的值进行下滤.
如果将2维的有序序列用一个arr[2,len]的数组来表达, 那么伪代码可以表示如下(双指针法):

SPLIT(arr,dim):                       %输入arr, 维度为2 X N的点集. 
                                      %以及本轮分割的维度dim, 
-------             
left = new Point[2,len/2];                    %左数组大小为 len/2, 
right= new Point[2,len-len/2];                %右数组大小为 len-len/2
for d=1:2:
	ptrL=1;         					    %左右指针
	ptrR=1;
	for i=1:len:
		if arr[d,i]<= arr[dim,len/2]           
			left[d,ptrL]=arr[d,i];      %小于本轮中位数的部分.  
			ptrL++;
		else
		    right[d,prtR]=arr[d,i];      %大于本轮中位数的部分.
		    prtR++;
return (left,right);

这个过程有bug吗? 我们考虑循环体维持了一个什么样的循环不变式. dim维最大值始终小于指定值的有序数组.

其他都好, 但发现这并不能明确指向中位值. 因为它的分割条件是arr[d,i]<= arr[dim,len/2];
因此在有重复的点时就分割不下去了(秩较低). 或者有大量重复元素时(比如大量相同Z值的点.), 可能会做无效分割. 抗退化性不是很好.
于是可以这么改写: if arr[d,i]<= arr[dim,len/2] && ptrL<= len/2
实际写的代码我附在文末, 没debug, 也懒得管了, 毕竟我就一刚学编程的新手. 还可以干脆就在Build前, 在维度分割前做一些筛选(根据序列的最大最小值). 输入上也做一些去重的限制

.
时间复杂度: 不含$O(nlogn)$的排序, 单次取中位数的时间复杂度是$O(n)$, 整树的build的耗时$O(nlogn)$ (有$logn$层, 每一层有一次$O(n)$)
而且缓存友好哦…

.
.实际代码我附在最后(C#), 仅作概念说明, 非bug free.

---

.

问题2-动态性.

第二个问题, 是可以专门划出一个专题去讲的
我进行一些学习后, 在这里采用了二进制分解的方式来进行, 正好可以将更新问题解决——二进制分解的特点是, 只有重建, 没有更新. 因此能回避任何更新问题.
简单减少一下二元分解(ref:《计算机图形学--几何体数据结构>):
假设有11个元素需要考虑, 111的二进制是1011, 一共四位,
那我们可以将kd树分解4棵. 每一棵上挂载$2^i$个元素($i$是树的编号, 从0开始). 由于这个二进制表达是唯一的, 因此任何数字都能表示为不同二进制基底的线性组合, 而查询的结果也等于四颗之和.

1. :[1]:1个元素,
2. :[1]:2个元素,
3. :[0]:0个元素;
4. :[1]:8个元素.

而当插入一颗元素时, 数量变更为1100, 树也重构为:

1. :[0]:0个元素,
2. :[0]:0个元素,
3. :[1]:4个元素;
4. :[1]:8个元素.

其中3th树是由低位重构而来的.

从这个角度来说, 树的更新只需重建, 而无需删除了.
关于删除, 这里有个细节:

经典KD树中会通过弱删除去加速——即在删除时只打标记, 直到超过阈值再重构(一般是一半的时候).

这个原则对二进制分解后的kd树也同样有效.
不过我在这里做了一些小创新, 重构变为局部重构, 重构线变为"局部容量"的一半 —— 当第$0-i$棵树上的弱删除总数突破第$0-i$上容量的一半时, 进行重构.
例如1111个元素中被删了4个, 而且只让低位变得不平衡, 那就仅对$0-2$号树进行重构(下图打底色部分).

我把Insert()部分的代码挂在了后头. 思路很简单, 但细节要正确表达出来还是可能有点困难的. (顺便说下个人试验性质的东西不保证真的更有效. ) 反正没啥人看, 就不解释了.

这里还有一点小心得, 如果你的任务局部性较强, 那么按有序列表建立的kd树可能更容易失去平衡. 比如你删除了xy值都比较小的元素, 那么树的左子树可能被一次性清空(如下图).
因此如果可能的话, 在预排序阶段对每个元素的索引进行排序, 而非元素本身.

---

XYZ是自定义的类, 一个double型的三维点.
class XYZ
{
   public double[] Coords = new double[3];
}
    
public (XYZ[,] L, XYZ[,] R) Split(XYZ[,] D_order, int dim)
{   //Contract:: D==D.Distinct
    int len = D_order.GetLength(0);
    if (len == 3 || len == 0) {
        return (D_order, new XYZ[3, 0]);
    }
    double med_dim = D_order[dim, (len - 1) / 2].Coords[dim];                       //本次中位数.
                                                                                    //ptrs
    int ptrL=0;
    int ptrR=0;
    //returnee
    XYZ[,] LSet_ord = new XYZ[3, (len - 1) / 2];
    XYZ[,] RSet_ord = new XYZ[3, len - (len - 1) / 2];
    
    for (int d = 0; d < 3; d++) {
        if (d != dim) {
            for (int i = 0; i < len; i++) {
                //dim0
                if (D_order[d, i].Coords[dim] <= med_dim) {
                    LSet_ord[d, ptrL++] = D_order[d, i];
                }
                else {
                    RSet_ord[d, ptrR++] = D_order[d, i];
                }
                ptrL = ptrR = 0;
            }
        }//直接取中位数, 没有这个if也行.
        else {
            for (int i = 0; i < LSet_ord.GetLength(d); i++) 
            { LSet_ord[d, i] = D_order[d, i]; }
            for (int i = (len - 1) / 2 + 1; i < len; i++) 
            { RSet_ord[d, ptrR] = D_order[d, i]; }
        }
    }
    return (LSet_ord, RSet_ord);
}

/// <summary>
/// 从高位开始逆序寻找重构位置. 
/// 这是一张『融合牌』. 每当你打出一组Insert, 它会在重建前计算最高受影响的位数; (你可以在ConsoleApp9的KDINSERT中找到它的简单测试.)
public void Insert(T[] D)
{
    //模型:
    //对w数组上的累积有效元素(acc_active)自高位向低位进行判断.存在2种state:
    // state ①: <进位态>, 此时假定元素超出低位容纳的极限, 或者低于降位的极限, 自己的低位是 (0000)1  or (1111)1, 存在两个临界点导致进降位.cap_carry /cap_down
    // state ②: <本位态>. 此时假定前方(1010), 而自身被进位时, 不会导致进位

    int acc_active = D.Length + roots.TotalActive; //
    int lim_carry; //carry  critical );   //进位极限
    int lim_Down; //[down critical  ;     //降位极限
    int rebuilt = roots.W.Count - 1;      //待重建的树的编号
    var capt = roots.Capacity;            //roots是树的集合. capacity是容量
    var binr = roots.BinaryRep;           //11001,,,二进制表达. 
    var actc = roots.ActiveCounts;        //激活的数量(对应删除的数量). 

    for (int i = roots.W.Count - 1; i >= 0; i--) {
        //state ①: 
        lim_carry = ((capt[i] << binr[i]) - 1);
        lim_Down = ((capt[i]) - 1) * binr[i];
        if (acc_active > lim_carry) { //e.g.: {1,2,4}  >7
            rebuilt = 1 + i; //point to the high 1.
            break;
        }
        else if (binr[i] == 1 && acc_active <= lim_Down) {
            rebuilt = i;
            break;
        }
        if (i == 0 && acc_active > 0) { //i=1 case otherwhise never change the first position.
            rebuilt = -1 + 1;
            break;
        }
        //state ②: 
        rebuilt--;
        acc_active -= actc[i];
        bool state2 = (acc_active - (capt[i] - 1)) > 0 && (acc_active - (capt[i] - 1)) <= roots.DeletedCounts[i] && binr[i] != 0;
        if (state2) {
            rebuilt = i;
            break;
        }
    }
    //
    this.Rebuild(rebuilt > roots.W.Count - 1 ? roots.W.Count - 1 : rebuilt, D);
}