K-D Tree

369Pai · 2023-01-25 14:00:31 · 个人记录

• 数据结构
  • K-D Tree
    • 简介
    • 建树
    • 插入/删除
    • 重构
    • 邻域搜索
    • 范围搜索
    • 复杂度退化
    • 变体
    • 练习

简介

k-D Tree(KDT , k-Dimension Tree) 是一种可以 高效处理 k 维空间信息 的数据结构，是一种二叉搜索树，它的每个节点都为 k 维点。所有非叶子节点可以视作用一个超平面把空间分割成两个半空间。节点左边的子树代表在超平面左边的点，节点右边的子树代表在超平面右边的点。选择超平面的方法如下：每个节点都与 k 维中垂直于超平面的那一维有关。因此，如果选择按照 x 轴划分，所有 x 值小于指定值的节点都会出现在左子树，所有 x 值大于指定值的节点都会出现在右子树。这样，超平面可以用该 x 值来确定，其法线为 x 轴的单位向量。

在结点数 n 远大于 2^k 时，应用 k-D Tree 的时间效率很好。

在算法竞赛的题目中，一般有 k=2 。在本页面分析时间复杂度时，将认为 k 是常数。

const double ALPHA = 0.8 , DOUBLE_MAX = numeric_limits<double>::max();;
typedef array<double , K>Point;
Point mn[N] , mx[N] ,  val[N] , a[N];//以i为根的子树k维坐标最小/最大值；树上结点；原来的点
int rt , tot , lc[N] , rc[N] , cnt[N] , s[N] , sz[N] , sd[N];//树根，节点数，子树大小，左儿子，右儿子
int tc , tmp[N] , *sgt , zero; //重构用数组
double ansmin = DOUBLE_MAX , ansmax = 0;//最近/最远距离

建树

假设我们已经知道了 k 维空间内的 n 个不同的点的坐标，要将其构建成一棵 k-D Tree，步骤如下：

1. 若当前超长方体中只有一个点，返回这个点。

2. 选择一个维度，将当前超长方体按照这个维度分成两个超长方体。

3. 选择切割点：在选择的维度上选择一个点，这一维度上的值小于这个点的归入一个超长方体（左子树），其余的归入另一个超长方体（右子树）。

4. 将选择的点作为这棵子树的根节点，递归对分出的两个超长方体构建左右子树，维护子树的信息。

其构建出 k-D Tree 的形态可能是这样的：

其中树上每个结点上的坐标是选择的分割点的坐标，非叶子结点旁的 x 或 y 是选择的切割维度。

这样的复杂度无法保证。对于 2,3 两步，我们提出两个优化：

1. 选择的维度要满足其内部点的分布的差异度最大，即每次选择的切割维度是方差最大的维度。

2. 每次在维度上选择切割点时选择该维度上的 中位数，这样可以保证每次分成的左右子树大小尽量相等。

可以发现，使用优化 2 后，构建出的 k-D Tree 的树高最多为 O(\log n) 。

附：优化 1 可换成周期循环所选择的维度，实现更简单。

现在，构建 k-D Tree 时间复杂度的瓶颈在于快速选出一个维度上的中位数，并将在该维度上的值小于该中位数的置于中位数的左边，其余置于右边。如果每次都使用 sort 函数对该维度进行排序，时间复杂度是 O(n\log^2 n) 的。事实上，单次找出 n 个元素中的中位数并将中位数置于排序后正确的位置的复杂度可以达到 O(n) 。要找到 s[l] 和 s[r] 之间的值按照排序规则 cmp 排序后在 s[mid] 位置上的值，并保证 s[mid] 左边的值小于 s[mid]，右边的值大于 s[mid]，只需写 nth_element(s+l,s+mid,s+r+1,cmp)。

void Build(int l , int r , int &p , int d)
{
    if(l > r)return ;
    int mid = (l + r) >> 1;
    auto Cmp = [&](const int& x , const int& y)
    {
        if(val[x][d] != val[y][d])return val[x][d] < val[y][d];
        return val[x] < val[y];
    };
    nth_element(tmp + l , tmp + mid , tmp + r + 1 , Cmp);
    p = tmp[mid] , s[p] = sd[p] = 1 , sz[p] = cnt[p] , lc[p] = rc[p] = 0;
    mn[p] = mx[p] = val[p];
    if(l == r)return ;
    Build(l , mid - 1 , lc[p] , (d + 1) % K);
    Build(mid + 1 , r , rc[p] , (d + 1) % K);
    Pushup(p);
}

重构

如果维护的这个 k 维点集是可变的，即可能会插入或删除一些点，此时 k-D Tree 的平衡性无法保证。由于 k-D Tree 的构造，不能支持旋转，类似与 FHQ Treap 的随机优先级也不能保证其复杂度，可以保证平衡性的手段只有类似于替罪羊树的重构思想。

我们引入一个重构常数 \alpha ，对于 k-D Tree 上的一个结点 x ，若其有一个子树的结点数在以 x 为根的子树的结点数中的占比大于 \alpha (或未删除的结点数在以 x 为根的子树中的占比小于 \alpha) ，则认为以 x 为根的子树是不平衡的，需要重构。

重构时，先遍历子树求出一个序列，然后用以上描述的方法建出一棵 k-D Tree，代替原来不平衡的子树。

bool IsUnbalanced(int p){return cnt[p] && (s[p] * ALPHA <= max(s[lc[p]] , s[rc[p]]) || sd[p] <= s[p] * ALPHA);}
void Travel(int p)
{
    if(!p)return ;
    Travel(lc[p]);
    if(cnt[p])tmp[++tc] = p;
    Travel(rc[p]);
}
void Rebuild(int &p)
{
    tc = 0;
    Travel(p);
    Build(1 , tc , p , 0);
}

插入/删除

插入一个 k 维点时，先根据记录的分割维度和分割点判断应该继续插入到左子树还是右子树，如果到达了空结点，新建一个结点代替这个空结点。成功插入结点后回溯插入的过程，维护结点的信息，如果发现当前的子树不平衡，则重构当前子树。

删除操作，则使用 惰性删除，即删除一个结点时打上删除标记，而保留其在 k-D Tree 上的位置。如果这样写，当未删除的结点数在以 x 为根的子树中的占比小于 \alpha 时，同样认为这个子树是不平衡的，需要重构。

void Insert(int &p , const Point& v , int d)
{
    if(!p)
    {
        p = New(v);
        return ;
    }
    if(v[d] < val[p][d])Insert(lc[p] , v , (d + 1) % K);
    else if(v == val[p])cnt[p]++;
    else Insert(rc[p] , v , (d + 1) % K);
    Pushup(p);
    if(IsUnbalanced(p))sgt = &p;
}
void Delete(int &p , const Point& v , int d)
{
    if(!p)return ;
    if(v[d] < val[p][d])Insert(lc[p] , v , (d + 1) % K);
    else if(v == val[p])cnt[p]--;
    else Insert(rc[p] , v , (d + 1) % K);
    Pushup(p);
    if(IsUnbalanced(p))sgt = &p;
}

邻域查询

问题：给定 N 个点的坐标，求 欧氏距离/曼哈顿距离 下的 最近/最远/第 K 远点对。

枚举每个结点，对于每个结点找到不等于该结点且距离最小的点，即可求出答案。每次暴力遍历 2-D Tree 上的每个结点的时间复杂度是 O(n) 的，需要剪枝。我们可以维护一个子树中的所有结点在每一维上的坐标的最小值和最大值。假设当前已经找到的最近点对的距离是 ans ，如果查询点到子树内所有点都包含在内的长方形的 最近 距离大于等于 ans ，则在这个子树内一定没有答案，搜索时不进入这个子树。

此外，还可以使用一种启发式搜索的方法，即若一个结点的两个子树都有可能包含答案，先在与查询点距离最近的一个子树中搜索答案。可以认为，查询点到子树对应的长方形的最近距离就是此题的估价函数。

注意：虽然以上使用的种种优化，但是使用 k-D Tree 单次查询最近点的时间复杂度最坏还是 O(n) 的，但不失为一种优秀的骗分算法，使用时请注意。在这里对邻域查询的讲解仅限于加强对 k-D Tree 结构的认识。

对于 K 远点对，可以将前 K 远的距离存在堆中，每次更新答案时保留前 K 远的距离，估价函数也用第 K 远距离即可。

double Dis(const Point& a , const Point& b)
{
    double res = 0;
    for(int i = 0 ; i < K ; i++)
        res += pow(a[i] - b[i] , 2);
    return res;
}
double Fmin(int p , const Point& v)
{
    double res = 0;
    for(int i = 0 ; i < K ; i++)
        res += pow(max(mn[p][i] - v[i] , 0.0) , 2) + pow(max((v[i] - mx[p][i]) , 0.0) , 2);
    return res;
}
double Fmax(int p , const Point& v)
{
    double res = 0;
    for(int i = 0 ; i < K ; i++)
        res += pow(max(abs(mx[p][i] - v[i]) , abs(mn[p][i] - v[i])) , 2);
    return res; 
}
void Qmin(int &p , const Point& v)
{
    if(!p)return ;
    ansmin = min(ansmin , Dis(val[p] , v));
    double disl = Fmin(lc[p] , v) , disr = Fmin(rc[p] , v);
    if(disl < ansmin)
    {
        if(disr < disl)Qmin(rc[p] , v) , disr = ansmin + 1;
        if(disl < ansmin)Qmin(lc[p] , v);
    }
    if(disr < ansmin)Qmin(rc[p] , v);
}
void Qmax(int &p , const Point& v)
{
    if(!p)return ;
    ansmax = max(ansmax , Dis(val[p] , v));
    double disl = Fmax(lc[p] , v) , disr = Fmax(rc[p] , v);
    if(disl > ansmax)
    {
        if(disr > disl)Qmax(rc[p] , v) , disr = ansmax - 1;
        if(disl > ansmax)Qmax(lc[p] , v);
    }
    if(disr > ansmax)Qmax(rc[p] , v);
}

范围搜索

复杂度退化

变体

//TODO K-D Tree 的线段树实现

练习

K-D Tree 基础练习 - Wallace

参考资料

K-D Tree - OI Wiki k-d tree - Wikipedia KDT 小记 - command_block 浅谈偏序问题与K-D Tree - 光与影的彼岸 K-D Tree - Tenshi