これは Competitive Programming Advent Calendar Div2013 の 20 日目の記事です．最近点対問題の話をします．
最近点対問題は，空間上に点の集合が与えられた時に，その中で最も距離が近いペアを探す問題です．
応用としては，何らかのオブジェクトを特徴ベクトルに写した後で，最も類似したペアを探したいときなどに役に立つのではないかと思います．競技プログラミング界では，今年のICPC会津大会で3次元上の最近点対問題が出題されました．
空間が平面のときは分割統治による $O(n\log{n})$ 時間アルゴリズムが有名かと思います．今回は一般の次元で，(ハッシュマップを用いて)線形時間で動く実装が楽そうな乱択アルゴリズムの紹介をします．参考にした文献は以下のサーベイです．
Smid, Michiel. Closest point problems in computational geometry. MPI Informatik, Bibliothek & Dokumentation, 1995.
定義
問題の定義は以下のとおりです．
入力 : $D$ 次元空間上の $n$ 個の点集合 $S=\{p_1,p_2,\cdots,p_n\}$
出力 : $\min_{i\not= j} d(p_i, p_j)$ (ここで $d(\cdot, \cdot)$ は2点の距離)
距離は色々種類がありますが，とりあえずユークリッド距離だとします．(が，後の証明を見て分かる通り，一般の $L_\tau$-距離 ($1 \le \tau \le \infty$) で成立します．$\tau=1$ のときはマンハッタン距離です．) あと，次元 $D$ は十分小さく，定数であるとみなします．
以降，記号を以下で定義します．
- 点集合 $T$，点 $p$ に対し $d(p, T) := \min_{q \in T, p \not= q} d(p, q)$ とします．
- 点集合 $T$ に対し $d(T) := \min_{p,q \in T, p \not= q} d(p,q)$ とします．
また，$\delta := d(S)$ とします．$\delta$ がこの問題で求めたい値です．
グリッド分割
低次元の最近点対問題を解く際，以下の補題が便利です．
補題 $c>0$ とする．任意の幅 $c\delta$ の超立方体の中に含まれる $S$ の点の個数は高々 $(cD+c)^D$ である．
証明 背理法による．点が $(cD+c)^D$ 個より多くあるとする．
元の超立方体をさらに小さい超立方体で分割し，小さい超立方体の個数が $(cD+c)^D$ になるように分割したとする．鳩ノ巣原理より，少なくとも1つの小さい超立方体は2つの点を含むが，このとき対角線の長さは高々 $c\delta/(cD+c) \cdot D \lt \delta$ なので矛盾．□
これより，最近点対の長さに比例する幅で空間をグリッド状に切ると，グリッドの各セルには高々定数個の点しか含まれないことがわかります．
また，もしこのグリッドの幅が $\delta$ より大きい (すなわち $c \ge 1$) であるとすると，任意の点 $p \in S$ について，$d(p,q) = d(p,S)$ となるような点 $q$ は，$p$ と同じセルか，あるいは $p$ が属するセルに隣接するセルに存在するはずなので，高々 $3^D \cdot (cD+c)^D$ 個の点しか調べなくて良いことになります．(ここで，2つのセルは少なくとも1点を共有する時に隣接しているとみなします．)
ハッシュマップを使うことで任意のセルにある点集合に $O(1)$ でアクセスできると仮定すると，このような幅がわかっていれば，最初にまずグリッドを構成し，各点 $p$ について，$p$ が属するセルと隣接するセル内にある点の中で $p$ と最も近いものの距離を計算することで，$O(n)$ 時間で $S$ の最近点対を計算できることがわかります．
というわけで，適切な幅を前もって見積もることができれば，問題を高速に解けることになります．が，以下の点に注意が必要です．
- 幅を小さくしすぎると (つまり $c < 1$ のとき)，最近点対は隣り合うセルに無いかもしれない(飛び越えてしまうかもしれない)ので良くない．
- 幅を大きくしすぎると1つのセルに大量の点が存在するかもしれない．
そのため，幅を適切に設定しなければ効率の良いアルゴリズムを設計することができません．
乱択アルゴリズム
以下ではこの補題を用いた2つの乱択アルゴリズムを紹介します．どちらのアルゴリズムも，最初に幅を十分大きな値に設定した後，何回か反復して少しずつ幅を小さな値に改良していきます．
篩(ふるい)法では，各反復で次のことをします：「現在の幅を使ってグリッドを構築し，グリッド上で孤立した頂点(隣接するセルに自分以外いないような頂点)の集合を取り除き，その上で幅を更新する．」こうすると，反復が進んで幅が小さくなっていくにつれて孤立した頂点の個数が増えていくので，反復が進む度に候補となる点集合が減っていきます．
逐次構成法では，まず最初に点の番号をランダムにシャッフルし，$i$ 回目の反復では点 $\{p_1, p_2, p_3, \cdots, p_i\}$ の最近点対の長さを計算するとともに，その長さを幅とするグリッドを構成します．こうすると毎回の反復でグリッドを構成しないといけないのでとても効率が悪い用に見えます．が，ここで点を追加した時に最近点対の長さが小さくなったときだけグリッドを更新すれば良いということを用います．最近点対の長さが次の反復で小さくなることは，最初に点をランダムにシャッフルしているおかげで滅多に発生しないため，実際には効率の良いアルゴリズムになります．
篩法
アルゴリズムは以下の通りです．
- $S_1:=S, i:=1$ とする．
- 以下を $S_i = \emptyset$ になるまで反復する．
- ランダムに点 $p_i \in S_i$ を選ぶ．$p_i$ と $S_i$ の最近点の距離を $d_i := d(p_i, S_i)$ とする．
- 幅 $d_i / (4D)$ でグリッドを作る．
- $S_i$ で孤立した点(グリッド上でその点が属するセル及び隣接するセルに，自分以外存在しないような点)を探す．そのような点の集合を $A$ とする．
- $S_{i+1} := S_i \setminus A$ とし，$i := i+1$ とする．
- $d_{i-1}$ を幅とするグリッドを作って最近点対の長さを求める．
このアルゴリズムによる解の正当性と，計算時間について見て行きましょう．以下，アルゴリズムが終了するタイミングの $i$ の値を $l$ と記します．
解の正当性
最終的には $d_{l-1}$ を幅とするグリッドを用いる手法に帰結しているので，$\delta \le d_{l-1}$ であることが言えれば良いです．が，これについては明らかです．($\delta$ が全点ペアの距離の最小値だったのに対し，$d_{l-1}$ は一部のペアの距離の最小値であるため．)
計算時間
計算時間の解析は，反復パートにかかる期待時間と，最後のグリッドに帰結したあとにかかる時間(幅 $d_{l-1}$ が大きすぎると計算時間がかかりすぎることに注意)の2つの部分に分かれます．
まず反復パートについて．
$i$ 回目の反復でかかる計算時間は，$O(|S_i|)$ です．よって全体の期待計算時間を見積もるには $i$ 回目の反復で残っている頂点の個数の期待値 $E[|S_i|]$ を見積もればいいことになります．
まず，$i$ 回目の反復とき，グリッドの幅を $d_i/(4D)$ と少し小さい値に設定していることから，以下のことが成り立ちます．
(♣) 点 $p\in S_i$ について，$d(p, S_i) > d_i/2$ ならば $p$ に隣接するセルに点は存在しない．そのため $p$ は次の反復では削除される：$p \not\in S_{i+1}$．
いま，$E[|S_{i+1}|] \le E[|S_i|]/2$ を示します．これが言えれば，等比級数の和の公式から，全体の期待計算時間を $\sum_i E[|S_i|] \le 2n$ で抑える事ができます．
ある数 $k$ を固定し，$S_i = k$ であるときを考えます．$S_i = \{r_1, r_2, \cdots, r_k\}$ を，$d(r_1, S_i) \le d(r_2, S_i) \le \cdots \le d(r_k, S_i)$ となるように番号付けします．このとき，性質 (♣) から，もし $p_i$ として $r_j$ が選ばれた場合，次の反復では $r_{j}, r_{j+1}, \cdots, r_{k}$ は削除されることが分かります．$p_i$ はこれら $k$ 個の中からランダムに選ばれることを考えると，(条件付き)期待値について，$E[|S_{i+1}| \mid |S_i|=k] \le k/2$ が成り立ちます．よって，
$E[|S_{i+1}|] = \sum_{k} E[|S_{i+1}| \mid |S_i|=k] \cdot Pr[|S_i|=k] \le E[|S_i|]/2$.
となり，反復パートの期待計算量が線形であることが言えました．
グリッドに帰結した後のパートについて．
これについては，$d_{l-1}/(4D) \le \delta$ が成立します．(∵2点 $p,q$ が最近点対で $d(p,q)=\delta$ だった場合，アルゴルズムが終了するまでのタイミングで $p,q$ は隣接するセルに入っていてはいけないが，そのためには刻み幅について $d_i / (4D) \le \delta$ でなければならない．)
よって刻み幅が区間 $[\delta, 4D\delta]$ に入っているので，効率よく計算することができます．
逐次構成法
今度は逐次構成法についてです．アルゴリズムは以下のとおりです．
- 入力の点集合をランダムにシャッフルし，$S = \{p_1, p_2, \cdots, p_n\}$ とする．以降，$S_i=\{p_1,p_2,\cdots,p_i\}$ と表記する．
- $\delta=d(S_2)$ とする．また，幅 $\delta$ のグリッドを構築する．
- 以下を $i=3,4,5,\cdots,n$ について反復する．
- $d_i = d(p_i, S_i)$ を計算する．これには前の反復で構築したグリッドを利用する．
- もし $d_i < d_{i-1}$ なら $\delta := d_i$ とし，幅が $d_i$ のグリッドを新しく構築する．
- $\delta$ を出力して終了．
このアルゴリズムが正しい解を出力するのは明らかです．なぜなら $i$ 回目の反復が終わった時点で $\delta = d(S_i)$ になっているからです．
というわけで計算時間について見ましょう．各 $d_i$ を求めるプロセスでは前のグリッドを利用できるので $O(1)$ で可能です．
よって，問題はグリッドを新しく構築するにどれだけ時間がかかるか，ということになります．
新しいグリッドの構築にかかる期待計算時間
今回のように，「ランダムに頂点を並べて順番に追加していく」というアルゴリズムでは，後向き解析という手法がよく取られるようです(参考：コンピュータ・ジオメトリの色んな章)． 後ろ向き解析では，頂点を順番に追加していく過程を逆順に考えます．つまり，最初にもう完成したものがあって，そこから頂点を1つずつ取り除く(このとき取り除く頂点は毎反復ごとにランダムに選ぶものと見なす)ときに発生する期待計算量を考える，というものです．期待値の線形性から，後向き解析が正しい期待計算量を与える事に注意します．
以降，後向き解析を考えます．いま，$i=n,n-1,n-2,\cdots$ の順で，頂点集合 $S_i$ からランダムに点 $p$ を選んで削除するというプロセスを考えます．(これは，通常のプロセスなら $S_{i-1} = S_i \setminus \{p\}$ のところに点 $p$ を追加することに相当します．) このとき，3つのケースが考えられます．
- どの点を削除しても最近点対の長さが変わらない．i.e. $\forall q \in S_i, d(S_i \setminus \{q\}) = d(S_i)$．
- $S_i$ の中に最近点ペアが複数個あり，すべてのペアがある点 $q$ を共有する．
- $S_i$ の中に最近点ペアが唯一つのみ存在し，$p$ を削除すると最近点対の長さが変わる．
ケース1の場合，$p$ の選択によらず再構築が発生しないので，計算コストは発生しません．
ケース2の場合，$p=q$ となった場合に限り再構築が発生します．これは確率 $O(\frac{1}{i})$ で起きます．ところで再構築に要する計算時間は $O(i)$ ですので，期待計算時間は $O(i \cdot \frac{1}{i}) = O(1)$ ということになります．
ケース3も同様です．
よって，全体で期待計算時間 $O(n)$ で最近点対の長さを計算することができます．
おわりに
最近点対問題の乱択アルゴリズムを紹介しました．篩法，逐次構成法について書かれた元論文は以下です．
- Khuller, Samir, and Yossi Matias. "A simple randomized sieve algorithm for the closest-pair problem." Information and Computation 118.1 (1995): 34-37.
- Golin, Mordecai J., et al. "Simple randomized algorithms for closest pair problems." Nord. J. Comput. 2.1 (1995): 3-27.
元のサーベイ論文では最近点対問題の色んな変種(動的な最近点対とか)も扱ってるっぽくてちょっと気になります(まだ読めてない)．
次は @skyaozora さんです．
