三维最近点对

基于分治法的三维最近点对问题的研究

摘要：最近点对问题被广泛应用于工程项目之中, 近年来, 由于大数据计算和人工智能等领域的兴起, 许多计算领域对三维最近点对问题的运用更加广泛. 分治法求解最近点对是计算几何领域的经典算法. 本论文通过对传统的二维最近点对问题进行分析并推广到三维空间下，利用分治法的思想和鸽巢原理，设计出两种高效求解三维最近点对的方法, 对比朴素做法，在时间复杂度上有显著的优化.

关键词：三维最近点对；计算几何；分治法；双指针; 鸽巢原理; 时间复杂度;

1. 引言

​ 分治法在计算机科学领域中占据着重要的地位.分治即“分而治之”，就是把一个复杂的问题分成两个或更多的相同或相似的子问题，再把子问题分成更小的子问题，直到最后子问题的规模足够小，不必再分割就可以直接得出它们的解，原问题的解即子问题的解的合并. 这个思想被应用于多种重要的算法设计中，如快速排序、归并排序、快速傅里叶算法等. 最近点对问题是计算几何领域的研究的重点问题，常常应用于大数据计算与空中交通的计算机自动控制系统中. 目前，研究方向已经往高维空间下发展.

​ 本论文所讨论的问题可以描述为：给定 个点，找出其中的一对点的距离，使得在这 个点的所有点对中，该距离为所有点对中最小的. 在实际情况中，可能会出现不止一个最近点对的情况，或者出现浮点误差. 为了规范化问题，我们只考虑得到在一定浮点误差内准确的答案即可.

2. 相关工作

2.1 二维最近点对朴素做法

2.1.1 算法描述

  对于二维最近点对的问题, 我们有朴素的做法: 维护一个最小值，直接对任意两个点，检验它们是否能够构成最小值，不断更新答案，最终得到的就是最近点对距离.

2.1.2 朴素做法的伪代码表示

step 1. input();
step 2. function solve()
            ans <- inf; /*设置为无穷大*/
            i <- 1;
            while i <= Point_num do 
                j <- i + 1;
                while j <= Point_num do
                       ans <- min(ans, distance(point[i], point[j])); /*对每两个点计算距离, 更新结果*/
                end while
            end while
        end function
step 3. output();

2.1.3 时间复杂度分析

  分析上述算法, 为了方便我们的研究，不妨视 函数和 函数为常数级别的运算, 则设总的点数为 总的计算次数是
$O(n^2)$级别的.

2.2 二维最近点对分治做法

2.2.1 算法描述

  利用分治法的主要思路是：对每一个点进行标号为 1 到 ，为了得到 区间里的最近点对，我们可以取一个中间值 , 把区间分成左右区间 . 分而治之求取 和 的最近点对，递归求解左右区间 , 最后再合并左右区间的结果。

  注意到我们是求 区间的平面最近点对, 那么我们从分治的思想去考虑, 能不能转化为求 和 的最近点对呢? 但是这样做有一个问题, 最近点对不一定单独出现在某个区间. 也就是说, 可能存在一个最近点对, 它的一个端点在, 另一个端点在.所以, 我们在分治的时候还需要对这种情况进行检验. 怎么检验呢？我们先对点集按 坐标升序排序, 当我们求得 和 的最近点对距离 和 后取， 在中找到跟中间点 的 距离小于 的点集 , 这样我们得到了候选点集 .

  实际上，我们得到的 是一个以 为中轴线，长为 , 宽为 的长方形， 中的点集具有稀疏的特点(这一点我们稍后会用鸽巢原理证明)，我们对 点集按 坐标排序, 可以在常数级别的时间内做到对点集进行检验，得到当前的最近点对距离.

2.2.2 分治法的实现

step1. input();
step2  merge_sort(point, x); /*对点集按x为关键字进行排序*/
step3. function solve(point, l, r)
            d <- inf;
            if l == r then
                return d;
            else if r - l == 1 then
                return distance(point[l], point[r]);
            end if;
            mid <- (l + r) / 2;
            d <- min(solve(l, mid), solve(mid + 1, r));
            num <- 0, i <- l;
            while i <= r do
                if fabs(point[i].y - point[mid].x) < d then /*得到点集P3*/
                    num <- num + 1;
                    tmp[num] <- point[i];
                end if;
            end while;
            merge_sort(tmp, y); /*对点集按y为关键字进行排序*/
            i <- 1;
            while i <= num do
                   j <- i + 1;
                while j <= num and fabs(tmp[j].y - tmp[i].y) < d do
                    d <- min(d, distance(tmp[j], tmp[i])); /*这里第二层循环运行次数为常数级别*/
                    j <- j + 1;
                end while;
            end while
       end function
step4. output();            

2.2.3 算法复杂度分析

 在上面代码中有一步按 坐标排序后找最小值的过程， 可以看到它有两层循环, 但其实实际的运算次数很少, 可以看做是一个常数，我们可以做以下证明：

  如图2-1所示, 我们在某一递归过程中出现最近点对的两个端点分别处于不同区域, 即一个端点属于 , 另一个属于
那么我们沿着中线 画出一个长 , 宽 的小矩形区域, 对于这个矩形区域, 最多不超过八个点, 这个是可以严格证明的.
  因为我们知道当前的最小距离是 , 那么对于左边的边长为 的正方形, 任意两点的距离都不会小于 ，显然, 最多的情况不超过四个点, 怎么证明呢? 这里我们可以利用反证法和鸽巢原理.

  如图2-2所示， 在上述的 边长为 的小正方形中, 我们划分成四个 的小正方形, 若这个区域存在至少五个点, 由鸽巢原理得到至少有两个点位于同一个小正方形, 而两者最远的距离为对角线距离为
$$
​ 这与矛盾, 因此这个区域的点不超过四个.

  同理, 对于右边的小正方形, 最多也只有四个点那么总的这个矩形区域的, 点最多不会超过八个点. 因此, 我们在上述的二重循环中最多只会做每个点之后的七个点, 那么这个操作可以看做是在常数时间内完成的, 于是, 该算法的分割步骤和合并步骤总共耗时是的, 排序的耗时是, 而对 排序耗时为常数级别.

  根据我们的分析, 算法耗费时间 满足递归方程
$$

​ 因此，上文提到的分治法求解而为最近点对问题总体时间复杂度是 .

3. 基于分治法的三维最近点对算法

3.1 三维最近点对分治做法

3.1.1 算法描述

  我们都知道二维平面最近点对的分治做法, 那么我们能否将这个算法推广到三维的情况呢?还是那个思路, 在求 1 到 的最近点对问题, 可以转化为求 到 的最近点对.

  同样的思路, 我们先对点按 坐标排序, 此时可以截取一个平面 ，对于平面的左右部分分治, 找到左右两端点集中最近距离 , 此时同样有几率会出现最近点对在分界面两边的情况, 比较难处理的是合并部分. 设左平面为 , 右平面为 , 此时我们考虑做对点集按 坐标进行排序，并进行第二次分治——对 坐标进行分治。

  第二次分治的时候，找到左右区间的最小值 , 并更新当前最小值 ，此时已经转换成平面内的最近点对问题，我们对点集进行 坐标排序后对每个点进行检验，不断更新答案. 这里还有一个地方需要注意，由于理论分析上， 坐标的地位是一致的，但是在实际情况中，由于数据的随机性，并不一定是按 分治, 再按 坐标分治最优，极端情况下, 可能出现最小距离比另外两维都要大，所以对 的分治复杂度不能保证. 可以根据实际情况调换分治策略，找范围最小的一维作 分治.

3.1.2 算法的伪代码表示

  基于上文的分析，我们可以给出两次分治求解三维最近点对的算法如下：

step 1. input();
step 2. solveX(v, l, r);
step 3. output();

function solvey(*v, l, r)  /*v是点集，l是左端点，r是右端点 对y做分治*/
    if l == r then 
        return inf; /* 一个点的情况不合法 */
    end if
    if l + 1 == r then
        return distance(v[l], v[r]); /*返回两点距离*/
    end if
    mid <- (l + r) / 2; /*分割左右区间*/
    d = min(solveY(v, l, mid), solveY(v, mid + 1, r));
    merge_sort(v + l, v + r + 1, z); /*按z坐标归并排序*/
    i <- l;
    while i <= r do
        j <- i - 1;
        while j >= 0 do
            if fabs(v[i].z - v[j].z) <= d then /*z 坐标的绝对值小于d, 为可选方案*/
                d <- min(d, distance(v[i], v[j]));
            end if
            j <- j - 1;
        end while
    end while
    return d;
endfunction

function solveX(*v, l, r) /*v是点集，l是左端点，r是右端点 对x做分治*/
    if l == r then 
        return inf;
    end if
    if l + 1 == r then
        return distance(v[l], v[r]); /*返回两点距离*/
    end if
    mid <- (l + r) / 2;
    midx <- v[mid].x;
    d <- min(solveX(v, l, mid), solveX(v, mid + 1, r));
    merge_sort(v + l, v + r + 1, y); /*按y坐标归并排序*/ 
    tot <- 0;
    i <- l;
    while i <= r do 
        if fabs(v[i].x - midx) then
            tmp[++tot] <- v[i];
        end if
        i <- i + 1;
    end while
    if tot > 0 then
        d <- min(d, solveY(tmp, 1, tot)); /*对y坐标分治*/
    end if
    return d;
end function

3.1.3 算法复杂度分析

  对于 中每个点在至多在 个区间内被纳入 ，所以他至多参与 次 ， 复杂度 ，所以总复杂度是 .

3.2 三维最近点对分治做法

3.2.1 算法描述

  我们考虑位于平面 左右两端的点， 找出左端距离平面小于 的点加入到点集 , 找出右端距离平面小于的点加入到点集 , 那么点集 到 的任意两点的距离都需要我们去验证是否为最近点对.

  基于上文的分析，当我们截取一个平面 时，它分隔了左右两个区间 和 , 最坏情况下, 和 的点集规模都是 的, 直接找复杂度是 . 但是, 假设当前我们得到的最近点对值为 , 我们只需要到平面 的距离小于 的点集 和 即可.

  对 的每一个点, 将它投影到平面 , 都可以对它建立一个 的长方体, 所有可选方案都将包含在长方体中, 如图3-1所示:

  可以通过鸽巢原理证明, 最多不超过 个点存在于这个长方体中. 那么我们先对 和 的点进行 为第一关键字, 为第二关键字的排序. 由于点集的 值具有上升的特点，可以通过双指针的方法维护一个满足条件的 坐标区间，其中左指针代表符合条件的最小 值，右指针代表符合条件的最大 值，那么在这个区间中，我们只需验证, 其中 为可以维护区间的点集规模.

3.2.2 复杂度证明

  证明思路是利用鸽巢原理进行构造, 只需进行适当划分，构造出矛盾关系. 把长方体分割成 个大小相等的 的小长方体, 若长方体中有多于 个点, 那么至少有一个小长方体里有两个点,.

   它们最远的距离是对角线的距离. 这与 是最近点对距离的事实矛盾, 因此大长方体里最多有 个点, 于是我们只需要通过双指针法维护一个可行的长方体区域即可，利用尺取法的思想, 每一次取一段满足条件的区域, 如果遍历超过 个点就跳出.

  根据我们的分析, 算法耗费时间 满足递归方程
可得, 故总体时间复杂度是.

3.2.3 伪代码实现

1. 
merge_sort(point, l, r, x); /*以x为关键字进行排序*/
2. 
function solve(l, r)
    if l == r then
        return inf;
    else if r - l == 1 then
        return distance(point[l], point[r]);
    end if;
    mid <- (l + r) / 2;
    d <- min(solve(l, mid), solve(mid + 1, r));
    num1 <- 0, num2 <- 0;
    i <- l;
    while i <= mid do 
        if fabs(point[i].x - point[mid].x) < d then
            tmp1[++num1] = point[i];
        end if
    end while
      while i <= r do 
        if fabs(point[i].x - point[mid].x) < d then
            tmp2[++num2] = point[i];
        end if
    end while
    merge_sort(tmp1, y, z); /*以y为第一关键字，z为第二关键字排序*/
    merge_sort(tmp2, y, z);
    pl <- 1, pr <- 1, i <- 1;
    while i <= num1 do
        if tmp1[i].y < tmp2[pl].y then 
            cur <- pl;
            while cur <= num2 and fabs(tmp2[cur].y - tmp1[cur].y) < d then 
                d <- min(d, distance(tmp1[i], tmp2[cur]));
                cur <- cur + 1;
            end while
        end if
    end while
    while pr <= num2 and pr < tmp1[i].y do
        pr <- pr + 1;
    end while
    while pl <= num1 and tmp2[pl].y < tmp1[i].y && fabs(tmp2[pl].y - tmp1[i].y) > d do
        pl <- pl + 1;
    end while
    while pr <= num2 and fabs(tmp2[pr].y - tmp1[i].y) do
        pr <- pr + 1;
    end while
    if(pr > pl) then
        cur <- 0;
        j <- pl;
        while j < pr do
            d <- min(d, distance(tmp1[i], tmp2[j]));
            ++cur;
            if cur == 24 then
                break;
            end if
           end while
     end if
     return d;
end function

4. 性能测试

  通过在随机生成的数据下考查算法的时间效率，并采用的在线判题系统的评测机作为速度评判标准，测试结果如表4-1所示：

输入点集规模 n	做法平均运行时间/	做法平均运行时间/
	25	19
	371	350
	1433	1140
	2951	2285

  由于普通个人计算机在1000 内的计算次数级别约为 , 大数据测试难以在个人计算机中实现，我们可以对上述结果进行拟合，画出图像进行预测,如图3-1所示:

  在实际的测试中，我们的算法实际表现的时间复杂度与理论分析差别不大，分别属于 和 级别.

5. 结论

  本论文通过对二维最近点对问题的分析，引入了三维空间下的最近点对问题，推广并提出了两种利用分治法求三维最近点对问题的高效算法. 在分析三维空间下的最近点对问题时，从朴素的做法开始入手，一步一步优化，先通过分治的思想将时间复杂度从 优化到了 ，再利用鸽巢原理证明了点集具有稀疏的特点，进一步将时间复杂度优化到了 . 在现今大数据和高性能计算领域里，对算法的时间复杂度提出了更高的要求，论文中提出的优化在实际计算中能够大大提升解决问题的效率. 例如，当问题的规模 达到 级别时，在普通个人计算机下朴素的 算法运行时间约为 100 秒，而使用 和 的算法可以在 1 秒内得到答案. 因此，本文提出的两种算法对比朴素做法性能提升较为显著.

参考文献：

[1] 王晓东. 计算机算法设计与分析[M]. 电子工业出版社, 2018.

[2] 严蔚敏. 数据结构[M]. 清华大学出版社, 2007.

[3] 张晓红，胡金初. 分治法实现最接近点对问题的三维推广算法. 山西师范大学学报(自然科学报)，2006.

[4] Donald E.Knuth. Sorting and Searching, volume 3 of The Art of Computer Programming. Addison-Wesley, 1973

[5] Robert E.Tarjan. Data Structures and Network Algorithms. Society for Industrial and Applied Mathematics, 1983.

[6] Alfred V.Aho, John E.Hopcroft, and Jeffrey D, Ullman. The Design and Analysis of Computer Algorithms. Addison- Wesley, 1974.

[7] Alfred V.Aho, John E.Hopcroft, and Jeffrey D, Ullman. Data Structures and Algorithms. Addison-Wesley, 1983.

[8] J.L.Bently. Writing Efficient Programs. Prentice-Hall, 1982.

[9] J.L.Bently. Programming Pearls. Addison-Weslet, Reading 1982.

[10] Jon Kleinberg, Va Tardos. Algorithm Design. Pearson Education, 2013.

[11] Steven S Skiena. The Algorithm Design Manual, 2nd Edition. Springer, 2011.

[12] E.Horowitz, S.Sahni. and S.Rajasekeran. Computer Algorithms/C++. Computer Science Press, 1996.

[13] Michael Ben-Or. Lower bounds for algebraic computation trees. In Proceedings of the Fifteenth Annual ACM Symposium on Theory of Computing, 80-86, 1983.

[14] Michael R.Garey and David S. Johnson. Computers and Intractability : A guide to the Theory of NP-Completeness. W.H.Freeman, 1979.

[15] Michael T.Goodrich and Roberto Tamassia. Algorithm Design: Foundations, Analysis, and Internet Examples. Wiley, 2001.

[16] Sara Baa***puter Algorithms Introduction to Design and Analysis, 3 rd edition[M]. Pearson Education, 2000.

[17] Shyu J, Yin P, Lin B. An ant colony optimization algorithm for the minimun weight vertex cover problem [J]. Annals of Operations Research, 2004, 131(2): 283-304.

[18] Raka Jovanovic, Milan Tubab. An ant colony optimization algorithm with improved pheromone correction strategy for the minimum weight vertex cover proble [J]. Applied Soft Computing, 2011, 11: 5360-5366.

[19] Satoshi Taoka, Toshimasa Watanabe. Performance comparison of approximation algorithms for the minimum weight vertex cover problem [J]. IEEE, 2012, 209(2): 632-636.

[20] Mletu R R, Plaxton C G. The online median problem[J]. SIAM J Comput, 2003, 32(3): 816-832.

[21] Bentely J L, Shamos M I. Divide and conquer in multidimensional space [C}]//Eighth ACM Symposium on Theory of Computing. ACM, 1976:220-230.