最近點對問題Closest Pair of Points是我當(dāng)年學(xué)算法時印象最深的一道題平面上給一堆點找出距離最近的兩個點。暴力解法分分鐘想到但數(shù)據(jù)量一上去就卡死。分治法的處理思路極其經(jīng)典配合畫圖理解后整個證明和實現(xiàn)都會變得異常清晰。這篇文章我不僅會把分治法講透還會用一套具體數(shù)據(jù)把每一步的圖都畫出來包括劃分、遞歸、合并、候選區(qū)域這些關(guān)鍵環(huán)節(jié)看完你就能完整復(fù)現(xiàn)整個算法。1. 最近點對問題從暴力到分治的演進1.1 問題定義平面上找距離最近的兩個點題目本身不復(fù)雜。假設(shè)平面上有 n 個點每個點用坐標(biāo) (x, y) 表示要求找出歐幾里得距離最近的一對點。注意是點到點之間的直線距離d sqrt((x1-x2)2 (y1-y2)2)這個問題的應(yīng)用場景比我最初想的多。比如地理位置聚類里要找出距離最近的兩個基站碰撞檢測里要找出空間中距離最近的兩個對象還有圖像處理中的特征點匹配第一步往往都是找最近鄰居。很多更大規(guī)模的問題會把最近點對當(dāng)作一個子過程反復(fù)調(diào)用所以這個基礎(chǔ)算法能不能跑得快直接影響整個系統(tǒng)的性能。我剛學(xué)的時候第一反應(yīng)就是暴力法把所有點兩兩組合都算一遍距離記錄最小值。兩層循環(huán)搞定代碼簡單到不用思考。但問題在于復(fù)雜度。n 個點兩兩組合有 C(n,2) n(n-1)/2 種情況時間復(fù)雜度是 O(n2)。n 等于 1000 的時候大約要算 50 萬次距離n 等于 10 萬的時候就是 50 億次。這個增長趨勢會迅速擊穿計算資源。實際項目中如果點集規(guī)模上了萬O(n2) 的暴力解法基本就是不可用的必須找更聰明的辦法。1.2 為什么暴力解會卡死規(guī)模增長的速度我拿一組實際數(shù)據(jù)感受一下。假設(shè)一秒鐘能算 1000 萬次距離已經(jīng)是相當(dāng)樂觀的估算不同規(guī)模的點集會呈現(xiàn)這樣的耗時點數(shù)量 n兩兩組合數(shù)量估算耗時1,000約 50 萬0.05 秒10,000約 5000 萬5 秒100,000約 50 億500 秒1,000,000約 5000 億50000 秒n 從 1 千漲到 100 萬點數(shù)量只多了 1000 倍耗時卻多了約 100 萬倍。這就是 O(n2) 的可怕之處。而且這還只是距離計算還沒算上距離比較、臨時變量維護這些開銷。真實場景里點動輒幾十萬個暴力法基本不可行。所以問題的關(guān)鍵變成了能不能把復(fù)雜度壓到 O(n log n) 甚至 O(n)分治法正好能做到 O(n log n)而且思路非常符合直覺就是大事化小、小事化了。2. 分治法的整體設(shè)計與畫圖拆解2.1 三步走劃分、遞歸、合并分治法處理最近點對的標(biāo)準(zhǔn)流程可以總結(jié)成三步。第一步劃分。把所有點按照 x 坐標(biāo)排序然后取中位數(shù)用一條垂直于 x 軸的線把點集分成左半和右半兩部分兩邊點數(shù)盡量均衡。第二步遞歸。分別對左半部分和右半部分遞歸地求解最近點對得到左半內(nèi)部的最近距離 dL 和右半內(nèi)部的最近距離 dR。遞歸出口是區(qū)間里只有一個點或兩個點一個點沒有距離直接返回?zé)o窮大兩個點直接算它們之間的距離。第三步合并。這是最關(guān)鍵的步驟。全局最近點對除了可能出現(xiàn)在左半內(nèi)部或右半內(nèi)部還有一種特殊情況一個點在左半一個點在右半它們跨越了分割線。這類跨線點對需要在合并時專門處理。合并的目標(biāo)就是找出所有跨線點對中的最近距離然后跟 dL、dR 比較取最小值。這個三步走看起來平淡無奇但合并步驟里藏著一個幾乎所有初學(xué)者都會忽略的優(yōu)化點跨線點對到底要檢查哪些點如果無腦把所有左半點和所有右半點兩兩比較那復(fù)雜度又回到 O(n2)分治就白做了。2.2 用一張圖畫清分治過程理解分治最好的方式就是畫圖。我用一組具體數(shù)據(jù)演示點集一共 8 個點隨便編一組坐標(biāo)但特意設(shè)計成跨線點對會成為最終答案這樣合并步驟的意義才能體現(xiàn)出來。點集分布示意 y 31 L3 30 L1 28 L2 ... 16 R3 15 R1 13 R2 8 6 4 B 3 L4 ----------------------------------------------------------- x 1 2 3 4 5 6 7 8 9 ... 20 21 228 個點按 x 坐標(biāo)排序后是編號坐標(biāo)左/右半?yún)^(qū)L1(1, 30)左半L2(2, 28)左半L3(3, 31)左半L4(5, 3)左半B(6, 4)右半R1(20, 15)右半R2(21, 13)右半R3(22, 16)右半初始狀態(tài)按 x 坐標(biāo)排序后取中間位置把點集切成左右各 4 個。分割線畫在 x 5.5 左右也就是 L4 和 B 之間。分割示意 y 31 L3 30 L1 28 L2 ... 16 R3 15 R1 13 R2 4 B 3 L4 |← 左半4個點 →|← 右半4個點 →| x5.5后面整個算法的推演我都會基于這張圖來講。你看到的是先豎著切一刀左右各自找最短再檢查跨線區(qū)域這套邏輯配上圖就完全不會亂。3. 合并步驟真正的大 Boss3.1 strip 區(qū)域哪些點值得檢查合并時的樸素做法是把左半所有點和右半所有點都兩兩算一遍。但這顯然不行。仔細(xì)想想就會發(fā)現(xiàn)大部分跨線點對根本不可能成為最優(yōu)解。已知當(dāng)前最佳的左右內(nèi)部距離是 d min(dL, dR)。如果跨線點對的距離要小于 d那么這兩個點的 x 坐標(biāo)都不能離分割線太遠(yuǎn)。具體來說任何一個點如果它的 x 坐標(biāo)與分割線的水平距離已經(jīng)大于等于 d那它跟對面任何一個點的距離都不可能小于 d。因為歐幾里得距離必然大于等于水平方向的差距。所以真正值得檢查的候選區(qū)域是以分割線為中心、寬度為 2d 的一條豎直帶狀區(qū)域通常叫 strip。分割線左邊 d 寬右邊 d 寬。strip 區(qū)域示意 y 31 L3 30 L1 28 L2 ... 16 R3 15 R1 13 R2 4 B 3 L4 | | | |← d →| | strip左邊界 分割線 strip右邊界 x5.5-d x5.5 x5.5dstrip 區(qū)域之外的點直接忽略因為它們在水平方向上就不夠格。這個剪枝讓候選點數(shù)量大幅減少。但問題還沒完全解決如果點分布很密strip 里可能依然有大量點全兩兩比較還是災(zāi)難。3.2 在 strip 內(nèi)按 y 坐標(biāo)排序滑動窗口的思想strip 內(nèi)的點不能繼續(xù)用 x 方向來剪枝了因為它們 x 方向都在分割線附近。這里換一個思路把 strip 內(nèi)的所有點按 y 坐標(biāo)排序。排序之后固定當(dāng)前點 p只需要檢查排在 p 后面的那些點 q。如果 p 和 q 的 y 坐標(biāo)差已經(jīng)大于等于 d那兩點距離必然大于等于 d直接停止檢查后面的點。因為 y 坐標(biāo)遞增后面的點只會離得更遠(yuǎn)。這個思路非常像滑動窗口對每個點只在一個高度為 d 的窗口里向后看。strip 內(nèi)按 y 排序后的檢查示意 y 坐標(biāo)升序排列 (5,3) ← p1 (6,4) ← p2 → 與 p1 的 y 差 1 d需要比較 如果某點的 y 與 p1 的差超過 d就跳過。這樣一來每個點不會跟 strip 里所有點比較最多只會跟 y 坐標(biāo)差距小于 d 的那些點比較。那么接下來還有一個關(guān)鍵問題在一個高度為 d 的窗口里最多可能有多少個點如果這個數(shù)量是常數(shù)整個合并步驟就能做到 O(n)。3.3 幾何證明的畫圖講解為什么最多只能檢查常數(shù)個點這個證明是整個分治法最近點對里最精彩的部分。想象一個寬 2d、高 d 的矩形p 點在這個矩形的左下角附近矩形內(nèi)部包含了所有可能與 p 距離小于 d 的候選點。2d × d 矩形內(nèi)的劃分示意 --------------------------- | | | | | 格 A | 格 B | 格 C | | | | | --------------------------- | p所在 | 格 E | 格 F | | 格 D | | | ---------------------------把這個矩形切成 2 行 3 列共 6 個小矩形。每個小矩形的寬是 2d/3高是 d/2。算一下對角線長度sqrt((2d/3)2 (d/2)2) d × sqrt(4/9 1/4) d × sqrt(16/36 9/36) d × sqrt(25/36) 5d/6對角線約等于 0.833d嚴(yán)格小于 d。這意味著什么如果一個格子里面放了兩個點它們之間的距離會被限制在這個格子的對角線以內(nèi)也就是必須小于 d。但我們知道 d 已經(jīng)是左右兩個子區(qū)域內(nèi)找出來的最短距離了整個平面上不可能再有距離小于 d 的點對存在否則左右遞歸早就發(fā)現(xiàn)了。所以每個小矩形里最多只能有 1 個點。6 個小矩形加起來最多 6 個點再加上 p 點自己也不超過 7 個點。就算不同教材畫法有差異有的切成 8 個小正方形核心結(jié)論都一樣這個區(qū)域內(nèi)點的數(shù)量是常數(shù)。所以每個點只需要向后檢查固定數(shù)量的幾個點就夠了不需要看完整條 strip。4. 手把手畫圖演示完整流程4.1 遞歸過程的畫圖推演回到我設(shè)計的那組數(shù)據(jù)。8 個點按 x 排序后先分成左半和右半然后兩邊各自遞歸。遞歸樹如下全部 8 個點 / \ 左半 4 個點 右半 4 個點 / \ / \ (1,30),(2,28) (3,31),(5,3) (6,4),(20,15) (21,13),(22,16) d√5 d√788 d√317 d√10 \ / \ / 左半合并 d√5 右半合并 d√5 \ / 全局合并跨線候選 d√2逐個看。左半 4 個點是 (1,30)、(2,28)、(3,31)、(5,3)。繼續(xù)劃分(1,30) 和 (2,28) 一組距離是 sqrt((1-2)2(30-28)2) sqrt(14) √5(3,31) 和 (5,3) 一組距離是 sqrt(4784) √788。左半合并時d 取 min(√5, √788) √5strip 區(qū)域?qū)挾?2√5里面包含 (1,30)、(2,28)、(3,31)進一步檢查后仍然沒有更小值所以左半內(nèi)部最近距離就是 √5。右半 4 個點是 (6,4)、(20,15)、(21,13)、(22,16)。(6,4) 和 (20,15) 距離 sqrt(196121) √317(21,13) 和 (22,16) 距離 sqrt(19) √10。右半合并時 d √10strip 內(nèi)按 y 排序后有 (21,13)、(20,15)、(22,16)其中 (20,15) 和 (21,13) 的距離是 sqrt(14) √5比 √10 更小于是右半最近距離更新成 √5。此時全局上界 d min(√5, √5) √5。4.2 合并過程的畫圖推演跨線點對登場合并左半和右半時分割線取 x 5.5當(dāng)前最優(yōu)距離 d √5 ≈ 2.236。strip 區(qū)域的范圍就是[5.5 - 2.236, 5.5 2.236] [3.264, 7.736]所有點里x 坐標(biāo)落在這個區(qū)間的只有兩個(5,3) 和 (6,4)。其他點要么在 strip 左邊要么在右邊直接排除。全局合并時的 strip y 4 B 3 L4 |← d2.236 →|← d2.236 →| x3.264 x5.5 x7.736 strip 內(nèi)只有兩個點 (5,3) 和 (6,4)strip 內(nèi)按 y 排序后依然是 (5,3)、(6,4)。兩個點之間 y 坐標(biāo)差 1小于 d所以必須計算距離sqrt((5-6)2 (3-4)2) sqrt(11) √2 ≈ 1.414√2 小于當(dāng)前全局上界 √5于是答案更新為 √2。最終最近點對就是跨線的 (5,3) 和 (6,4)。這組數(shù)據(jù)特別適合演示分治的合并威力左半內(nèi)部和右半內(nèi)部各自的最短距離都是 √5但全局最短卻是跨線的 √2。如果只看左右兩半的結(jié)果永遠(yuǎn)發(fā)現(xiàn)不了這個更優(yōu)答案。這也是為什么合并步驟絕不能省略的原因。4.3 畫圖理解為什么合并不會漏掉答案很多第一次學(xué)的人都會擔(dān)心strip 里只檢查這么點候選會不會把正確答案漏掉我當(dāng)時也有這個懷疑。畫一遍圖之后其實就清楚了。假設(shè)真正的最優(yōu)跨線點對是 a 和 b距離小于當(dāng)前 d。a 和 b 分別在分割線兩側(cè)它們每個點到分割線的水平距離都必然小于它們之間的距離也就必然小于 d。所以 a 和 b 都一定在 strip 區(qū)域內(nèi)不會被剪枝剪掉。這是合并步驟正確性的第一重保障。第二重保障是 y 坐標(biāo)窗口。a 和 b 的 y 坐標(biāo)差不會超過它們的距離也就小于 d。所以在按 y 排序后的序列里a 和 b 一定落在同一個高度為 d 的窗口內(nèi)檢查 a 時一定會看到 b。兩重保障疊加答案不可能漏只是需要確保 check 的時候用嚴(yán)格的小于當(dāng)前最優(yōu)來判斷。這不是玄學(xué)是幾何性質(zhì)給出的硬保證。5. 復(fù)雜度與正確性分析5.1 遞推公式與 O(n log n)分治最近點對的時間復(fù)雜度需要認(rèn)真算一下。設(shè) T(n) 是處理 n 個點的時間遞歸時切成兩個規(guī)模為 n/2 的子問題所以有 2T(n/2)。合并部分構(gòu)造 strip 需要遍歷 n 個點判斷是否在區(qū)域內(nèi)一次線性掃描O(n)。strip 內(nèi)按 y 排序時如果每一層都重新排序排序是 O(m log m)其中 m 是 strip 內(nèi)點數(shù)最壞情況下 m 接近 n合并就是 O(n log n)。于是遞推公式有兩種情況如果合并前對 strip 重新排序T(n) 2T(n/2) O(n log n)解出來 T(n) O(n log2 n)。如果做預(yù)排序優(yōu)化遞歸過程中始終保持 y 序列有序合并時直接線性歸并T(n) 2T(n/2) O(n)解出來 T(n) O(n log n)。實際工程里我大多數(shù)時候用第一種簡化實現(xiàn)因為數(shù)據(jù)量在幾萬以內(nèi)兩者差距不大。但如果是競賽題或超大數(shù)據(jù)量必須用預(yù)排序優(yōu)化。預(yù)排序的做法是所有點一開始按 x 排序進入遞歸同時在主函數(shù)里維護一個按 y 排序的全局序列遞歸時把序列按 x 中位數(shù)分成兩份合并階段線性地把兩個有序 y 序列 merge 成一個。這樣每一層都是 O(n)總復(fù)雜度穩(wěn)定 O(n log n)。5.2 正確性論證的三個關(guān)鍵點分治算法能保證正確性靠的是三個環(huán)環(huán)相扣的結(jié)論。第一最優(yōu)解只有三種可能位置全在左半、全在右半、跨分割線。遞歸處理前兩種合并處理第三種覆蓋了所有情況。第二跨線最優(yōu)解的兩個點必然都在 strip 內(nèi)。這點前面已經(jīng)證明過水平距離小于等于兩點實際距離實際距離小于 d 就意味著水平距離小于 d。第三strip 內(nèi)按 y 排序后每個點只需要檢查后面有限個點。因為在一個高為 d 的窗口內(nèi)點的數(shù)量有常數(shù)上界。只要窗口內(nèi)的候選點沒有被漏掉檢查就不會錯。這三個結(jié)論放在一起算法在邏輯上就是完備的。你可以用隨機數(shù)據(jù)反復(fù)驗證我實測過大量隨機點集和暴力法對拍結(jié)果完全一致。5.3 排序優(yōu)化的細(xì)節(jié)為什么預(yù)處理后合并可以做到線性經(jīng)常有人問既然 strip 要按 y 排序能不能每層直接 sort可以但不優(yōu)雅。每層 sort 會讓復(fù)雜度多一個 log。預(yù)排序的核心思路是用空間換時間。一開始就把點按 y 坐標(biāo)排好序然后在遞歸過程中保持這個順序不被打亂。具體實現(xiàn)是每層遞歸拿到的是已經(jīng)按 x 分好的左右兩個子數(shù)組以及一個按 y 排好序的全局?jǐn)?shù)組。按 y 排序的數(shù)組里逐個判斷點屬于左半還是右半分別放進兩個臨時數(shù)組合并階段再把兩個臨時數(shù)組歸并回 y 有序數(shù)組。這樣每個點只在每一層被常數(shù)次掃描總時間 O(n log n)。這個優(yōu)化對理解很重要。我在一開始學(xué)的時候只滿足于能算出正確答案沒想過排序帶來的額外 log后來做大規(guī)模數(shù)據(jù)處理時才意識到這個差異。如果你的數(shù)據(jù)量超過幾十萬每層 sort 的 O(n log2 n) 會明顯拖慢速度預(yù)排序幾乎是必須的。6. 工程實現(xiàn)與踩坑筆記6.1 一份可以直接跑的 Python 參考實現(xiàn)下面是我實際常用的分支限界風(fēng)格寫法加了不少注釋。我用平方距離比較來避免開根號的浮點誤差最后一步再開根。import math def sqdist(p1, p2): dx p1[0] - p2[0] dy p1[1] - p2[1] return dx * dx dy * dy def closest_pair(points): pts sorted(points) # 按 x 排序相同 x 按 y n len(pts) if n 2: return float(inf), None, None return _solve(pts, 0, n - 1) def _solve(pts, l, r): if l r: return float(inf), None, None if l 1 r: return sqdist(pts[l], pts[r]), pts[l], pts[r] mid (l r) // 2 mid_x pts[mid][0] dl, al, bl _solve(pts, l, mid) dr, ar, br _solve(pts, mid 1, r) if dl dr: best, pa, pb dl, al, bl else: best, pa, pb dr, ar, br # 收集 strip 內(nèi)的候選點 strip [] for i in range(l, r 1): if abs(pts[i][0] - mid_x) best: strip.append(pts[i]) # 按 y 排序 strip.sort(keylambda p: p[1]) # 滑動窗口檢查 m len(strip) for i in range(m): for j in range(i 1, m): if strip[j][1] - strip[i][1] best: break d sqdist(strip[i], strip[j]) if d best: best, pa, pb d, strip[i], strip[j] return best, pa, pb調(diào)用示例points [(1,30), (2,28), (3,31), (5,3), (6,4), (20,15), (21,13), (22,16)] d2, a, b closest_pair(points) print(a, b, math.sqrt(d2))輸出結(jié)果是(5, 3)(6, 4)1.4142135623730951和前面手工推演完全一致。6.2 邊界條件與遞歸出口遞歸出口是三件容易出問題的事。首先是單點情況區(qū)間里只有一個點沒有點對直接返回?zé)o窮大。我的代碼里用l r處理。其次是兩個點的情況直接計算距離返回避免繼續(xù)劃分。不寫這個出口的話兩個點會被切成一個點加空區(qū)間空區(qū)間返回?zé)o窮大合并時 strip 掃一遍也能得到正確結(jié)果但多一層無意義遞歸代碼也不清晰。第三是重復(fù)點。如果兩個點坐標(biāo)完全相同理論上它們的距離是 0是全局最優(yōu)算法會立刻發(fā)現(xiàn)。但要注意 strip 收集條件里的浮點比較。如果 best 已經(jīng)變成 0那 strip 里只有 x 完全相同的點滑動窗口依然能處理。最穩(wěn)妥的做法還是在主函數(shù)一開始先做一次去重檢測存在重復(fù)點直接返回 0。6.3 浮點精度與比較策略距離比較有一條小經(jīng)驗不要在遞歸內(nèi)部用math.sqrt。開根號是浮點運算不但慢而且累積誤差會讓比較變得不干凈。我全程用平方距離最后在返回時才開一次根。這樣既快又準(zhǔn)。另一個容易踩的坑是 strip 邊界的等號判斷。我建議用而不是。因為如果兩個點的距離恰好等于當(dāng)前最優(yōu)值用把它們排除在 strip 外也不影響更新答案反正不會變得更優(yōu)但用可以避免一些邊界爭議實現(xiàn)上更安全。不過要注意如果整個算法允許返回等距的多個解那必須用并且在比較時用更新。通常題目只要求任意一個最近點對所以代碼里不必過度糾結(jié)。6.4 常見問題速查表現(xiàn)象可能原因解決方法答案偏大合并時漏了 strip 區(qū)域內(nèi)的候選點檢查 strip 收集條件確認(rèn)用的是當(dāng)前最優(yōu)值 d答案偏大滑動窗口內(nèi)提前 break確認(rèn) break 條件用的是 y 差值大于等于 d而不是大于復(fù)雜度退化到 O(n2)strip 內(nèi)所有點兩兩都比較必須按 y 排序并用 y 差值剪枝重復(fù)點返回異常沒有處理坐標(biāo)完全相同的點主函數(shù)開頭特判存在重復(fù)點直接返回 0浮點誤差導(dǎo)致邊界判斷失誤遞歸內(nèi)部用開根后的距離做比較全程用平方距離比較最后再開根輸出點對與暴力法不一致多個最優(yōu)解選擇了不同的一對只要距離相等任何一對都算正確不用糾結(jié)7. 寫在最后畫圖理解與個人體會我個人在實際學(xué)習(xí)中的體會是像最近點對這種幾何味很重的算法光看文字推導(dǎo)效率很低畫圖幾乎是必須的。把點標(biāo)在坐標(biāo)系上畫出分割線、畫出 strip、畫出按 y 排序后的滑動窗口整個算法的每一步都有了直觀載體。我當(dāng)年就是親手在紙上把 8 個點的例子推演了一遍才徹底搞懂為什么合并階段只需要檢查常數(shù)個點。如果你也在學(xué)這個算法建議拿筆照著文中的 8 個點自己畫一遍遞歸樹再畫一遍 strip 區(qū)域。畫到合并那一步看到 (5,3) 和 (6,4) 這兩個跨線的點在 strip 里相遇時那種原來如此的感覺是看多少篇文章都換不來的。最后再分享一個小技巧寫完分治實現(xiàn)后一定寫一個暴力解法做對拍隨機生成幾千組點集逐一比對結(jié)果。這個習(xí)慣幫我排查掉了大量邊界 bug也是我寫任何算法題都會做的第一件事。算法寫出來是給人看的但正確性最終要靠數(shù)據(jù)說話。