点を円で囲む方法のあれこれ

2013/5/8 不正確な点もあるので、そのうち書き直します
2014/10/15 警告文追加しました

MM56の問題でこんなものがありました。

平面上のすべての点をM個以下の円で覆い、それらの半径の2乗和を最小化せよ

細かい条件は省略していますが、これの厳密解を求めるのは困難ですので、クラスタリングを行い、それぞれのクラスタについてすべてのの点を囲む最小円(最小包含円)を求めるのがたぶん定石だと思います。ここで問題になるのは、どうやって最小包含円を求めるかです。そのとき調べたところいろいろあったのでまとめてみました。計算量のところの、nは点の個数、epsは半径に対する相対誤差です。

予備知識

まず、知っておくべきことは三角形の最小包含円です。これは、鋭角三角形と鈍角三角形で場合分けできて、円の中心は

鋭角三角形のとき、三角形の外心
鈍角三角形のとき、鈍角の対辺の中点

となります。

n^3の方法

一般の多数の点について最小包含円を求めるには、局所改善を考えると、すべての3点の組み合わせについて最小包含円を求めれば、これらの円のどれかはすべての点を囲むことができるはずなので、それぞれ包含できているかを確認していけばよいことになります。各組み合わせに関して、最小包含円は定数時間で求められるので、全体ではO(n^3)です。

n^2log(1/eps)の方法

ご存じのように円を指定する方法には、2点と半径を決めるというものもあります。半径rを固定すれば、2点を通る円は高々2個なので、これらの円の中に包含円があるかはO(n^2)でできます。最小半径よりもrが大きいときには必ず包含円が存在するので、これを使って2分探索などで最小値を求めることができて、必要な反復回数はO(log(1/eps))であり、全体でO(n^2log(1/eps))となります。

誤差があると思うかもしれませんが、十分高い精度で中心を求めることができれば、最小包含円を決める3角形をO(n)で正しく求めることができるので、そこから誤差なしで求めることもできます。

nlog(n)の方法

計算幾何の本などによく載っている方法で、最遠点ボロノイ図を用いる方法があります。それを使うとO(nlog(n))でできるそうですが、詳細は省略します。

そのほか、確率的にO(nlog(n))になる方法があります。まず点を定数個選んでその最小包含円を求め、それが全体の最小包含円になっていたら終了、それ以外なら、円に入らなかった点の選ばれる確率を2倍にするという方法です。この方法では、包含円上の点を含む外側の点は、他の点の組み合わせに対する包含円に入りにくいはず、ということを利用して高速化を図っています。だいたいlog(n)程度繰り返すと外側の点の重みが十分大きくなって、最小包含円が求まります。

同じような考え方として、前の方にある点の組み合わせで包含円を作り、入らなかった点を前の方に持ってきて優先的に調べるようにするという方法も考えることもできて、たぶんO(nlog(n))くらいになると思いますが、正確なことはちょっと分かりません。

枝狩りの方法

実際のところ、凸包上の点だけで最小包含円を決めることができます。一応これは最遠点ボロノイ図の性質とも関係しています。凸包の計算はO(n^2)、あるいはO(nlog(n))でできるので、これによって候補点を減らして計算量を節約できる可能性があります。

nlog(1/eps)の方法

プロムナードの記事の後ろの方に書いてある方法ですが、nステップ目に一番遠い点との差分ベクトルの1/2^n倍を移動することによって、円の中心を求めることができます。円内の点、すなわち囲まれる点のどれかを初期点に選ぶと、1回の反復で真の中心からの距離の最大値が少なくとも1/sqrt(2)倍になるので、2log(1/eps)回反復すれば、所望の相対誤差で中心を求めることができます。1番遠い点を求めるのにO(n)かかるので、全体でO(nlog(1/eps))で求まります。

2014/10/15 追記: 上の方法は距離をおそらく縮めすぎです。1/sqrt(2)^n倍くらいに抑えないといけないはずです。(未検証)

その他Spaghttiに載っている方法ではO(n)でできるそうですが、あまり理解できていないので省略します。

結論

最後の点を移動させていくやり方が1番速くて、しかも式が単純で分かりやすいのでそれを使うのがいいと思います。工夫すると計算量が劇的に落ちるといういい例でしょう。

気が向いたら、多次元とか円が複数とかについて補足の記事も書こうかと思います。