ゆるふわブログ

過去の私の墓場

ボロノイ図を描く！ ~Fortune's Algorithm~

こんにちは！ Ysmr-Ry です．
Twitter の知り合いにボロノイ図の論文を見せてもらって興味が湧いたのでボロノイ図の描き方を勉強してみました！
ボロノイ図とはなんなのか，そしてその描き方を見てみましょう．
「ビーチに波が押し寄せる」ようにボロノイ図ができてゆきます．
長くなってしまいましたが，なぜ放物線なのかがわかれば十分な気がします．

ボロノイ (Voronoi) 図とは

ボロノイ図とは，与えられたいくつかの母点のうちどれに一番近いかで領域を分けた図のことです．
実際にボロノイ図を見てみましょう．

f:id:Ysmr_Ry:20170906172456p:plain

これがボロノイ図です．黒い点が母点 (ボロノイ図を描くために与えられた点) で色分けされてるのがそれぞれのボロノイ領域です．
よく見てると気づくとは思いますが，隣接するボロノイ領域の境界線は対応する母点を結んだ線分の垂直二等分線です．

球面上で作ったり，いろいろ拡張はあるようです．
次のリンク先は地球上の空港を母点として球面上にボロノイ図を描いています．

World Airports Voronoi

Fortune's Algorithm

ボロノイ図の描き方は当然いろいろあるとは思いますが，今回紹介するのは Fortune's Algorithm という方法です．
Fortune's Algorithm は平面走査法の一種で，円錐曲線としての放物線に着想を得ています．
簡単言うと，ボロノイ図を描きたい平面上を横線が上から下に動いていってスキャンし，それに放物線を重ね合わせた "波" が続いて，波の継ぎ目がボロノイ領域の境界をなぞっていきます．
実際見た方が早いですね．こんな感じです．(点の配置によってはバグるのでそういう時はうまくいくまでページを更新してください)

assqr.github.io

水色の線は放物線の組み合わせでできている "波" で，それよりも上の領域が "海" と考えると良いです．
黒い横線が平面上を走査して行って，緑色のボロノイ図が描かれていくのがわかります．
赤い丸はなんなのか，右の木 (こういう形のデータ構造を木と言います) は何を表しているのか，はこの後説明します．

前提知識 1 (円錐曲線としての放物線)

何と言ってもこのアルゴリズムの面白いところは放物線が生きてくるということです．
放物線，2 次関数自体は中学生でもやると思いますが，ここで使う放物線の性質は現在の課程では数 III の 2 次曲線のところで出てきます．
それは，「ある定直線 d とその上にない定点 F があるとき， d と F から等距離にある点の軌跡は放物線である」です．
定直線 d を準線 (directrix)，定点 F を焦点 (focus) と言います．
アニメーションで見るとこんな感じです．

紫の定点と定直線がそれぞれ焦点と準線，赤いのがそれによってできる放物線，オレンジ色の円は焦点を通り定直線に常に接しているような円です．
ということはこの円の中心は焦点と準線から同じ距離 (半径) にあるので，放物線の上を動くというわけです．
これをより一般化して焦点と準線からの距離の比が \varepsilon となるような点の軌跡を離心率 \varepsilon の二次曲線といい，放物線・楕円・双曲線になります．

前提知識 2 (Binary Search Tree)

先ほどの Voronoi Visualizer の右側に表示されていた木を二分探索木 (Binary Search Tree; BST) といいます．
木というのは先ほど見ていただいた通りの形をしたデータ構造なのですが，中でも二分木はそれぞれの頂点 (丸のこと) から伸びる枝が高々二股以下に分かれているものを言います．
それぞれの頂点には例えば数字が入っています．二分探索木は，「ある頂点に注目した時にその中に入っている値に対して，そこから左に生えている枝の先の部分木に入っている値は全てそれよりも小さく，右は全て大きいという性質を保つもの」です．
図で見ていきましょう．次のような木は二分探索木と言えます．

f:id:Ysmr_Ry:20170906184935p:plain

一番上の 5 の入った頂点 (根と言います) に注目します．すると，左の部分木は下の四角で囲った部分です．

f:id:Ysmr_Ry:20170906185322p:plain

見てみるとこの中のどの値も 5 より小さくなっていることがわかります．(1, 2, 3 < 5)
また右の部分技についてはどの値も 5 より大きくなっていることがわかります．(5 < 6, 8, 11)
今回は根に注目しましたがこれの性質はどの頂点についても成り立っています．(例えば 2 の頂点なら 1 < 2 で 2 < 3 です)
こういう風に数値のデータを持っていて何がうれしいかというと，探索木ですから，この中からある値を探す時に効率がいいのです．
例えばこの二分探索木から 6 を探したいとしましょう．
根から順にたどっていきます．根に入っている値は 5 です．5 < 6 なので，先ほど見た性質から，6 は少なくとも右の部分木にあることがわかるので，右の枝を辿ります．
これを繰り返せば，最悪の場合でも木の深さ分だけ比較をすれば見つかることがわかります．

f:id:Ysmr_Ry:20170906190834p:plain

木の深さ d に対してその深さにある頂点はだいたい 2^d 個あるので計算量は (平衡状態であれば) O(\log n) となります．
Fortune's Algorithm では二分探索木を放物線の波の最前線を管理するのに使います．

Fortune's Algorithm (なぜ放物線を考えるのか)

では本題の Fortune's Algorithm について考えていきましょう．
Fortune's Algorithm は平面走査法なので，走査線 (Sweep Line) が平面上を上から下に動きます．線が掃いていくから Sweep です．
ただこのアルゴリズムの場合は特殊で，もう一つ放物線の組み合わせの "波" である汀線 (Beach Line) を考えます．
今回はボロノイ図を構成するためのアルゴリズムなので，走査線よりも上の掃き終わった領域においてはボロノイ図の一部が組みあがっていて欲しいのです．
しかしながら，走査線よりも上ならば下にある母点たちの影響を全く受けないとは言えません．
そこで，それよりも上ならば下の母点たちから一切影響を受けないというもう一つの境界線，汀線を考えるのです．
先ほど言ったように汀線の上側の領域を "海" と呼称すると，ボロノイ図の出来上がった領域である "海" がじわじわと平面全体を覆っていくということです．
ではなぜその汀線が放物線となるのか．それは前提知識 1 で見た放物線の性質によるものです．
汀線は走査線が掃き終わった，上にある全ての母点をそれぞれ焦点とし，現在の走査線を準線とした放物線よりも上の領域のうち，もっとも下にある境界線をつぎはぎしてできた最前線として決めます．
例えば，先ほどの Visualizer で放物線を全て表示すると

f:id:Ysmr_Ry:20170906193730p:plain

となるのですが，実際に考えるのは一番下にある，最前線の放物線の組み合わせです．
そうするとそれらの放物線は「走査線と対応する母点から等距離にある点の軌跡」なので，放物線の上の領域は「走査線よりも母点の方が近い点の集まり」になります．走査線よりも下の点は当然走査線よりも遠いですから，"海" は走査線よりも下の点から影響を受けないことがわかります．
また，隣接する放物線の交点を Breakpoint と呼ぶのですが，その点は準線と二つの母点からの距離が等しいのでボロノイ図の辺の上に乗ります．

Fortune's Algorithm (2 つのイベント)

平面走査法はアニメーションとしては連続的に走査線が動いていくのですが，実際のプログラムでは特別なイベントが起こる時点をとびとびにたどっていきます．
今回の場合そのイベントは 2 つあります．

Circle Event

f:id:Ysmr_Ry:20170906201010p:plain

図のように 3 つの母点が同一円周上にあり，その円の下側に走査線が接する時点を言います．(したがって上の図は Circle Event が起こる少し前)
このまま走査線が下に行くと，赤い円の中心に集まってきている 3 つの放物線の一部のうち，真ん中が完全に汀線からなくなります．
また，円の中心は Circle Event のまさにその時点では 3 つの放物線に同時に乗っていると言うことができ，それに対応する 3 つの母点 (と走査線) から等距離にあるため，ボロノイ図における点 (ボロノイ点) になります．
したがってこのイベントが起きた時は真ん中の放物線を汀線の放物線のリストの中から消去し，また，ボロノイ点で終わっている辺の終わりにボロノイ点をセットし，ボロノイ点から始まる辺の始点をセットします．

Site Event (NewPoint Event)

f:id:Ysmr_Ry:20170906210749p:plain

Site Event は母点の一つに走査線が触れた時点で起きます．上の図は走査線が 3 の母点を掃いた直後です．
母点 3 を焦点とする放物線が生じるため，それを汀線に組み込まねばなりません．

Fortune's Algorithm (実装の詳細)

先ほどのイベントは上から順番に並べて順に処理します．この時 y 座標が最小のものを効率よく取り出すデータ構造が必要で，Priority Queue を使うみたいです．最初は母点から Site Event だけを作って並べて置いて，それを処理ながら Circle Event を発生させてイベントリストに追加します．
このアルゴリズムで使うもう一つのデータ構造は先ほども紹介したように二分探索木です．
頂点の持つ値として何を入れているかというと，Breakpoint の x 座標です．
Visualizer で (1,2) と頂点に表示されていたら，母点 1 と母点 2 を焦点とする放物線が左からこの順で交わっている時の交点の x 座標が入っています．
Site Event のところの図がわかりやすいと思います．(再掲)
f:id:Ysmr_Ry:20170906210749p:plain

この場合は，(0,2) (2,3) (3,1) (1,0) の順で左から Breakpoint が並んでいるのがわかります．(このうち画面内に見えているのは (2,3) (3,1) だけですが…)
ちなみにこの順序は二分探索木を深さ優先探索 (Depth-First Search) すると得られます．

Circle Event においては真ん中の放物線を消し (つまり真ん中と左，真ん中と右の作る Breakpoint を消す)，新たな Breakpoint として円の中心を追加します．
Site Event では新たに放物線ができるので，既存の汀線とどこでぶつかるかを調べて，新たに Breakpoint を追加しなければなりません．
関係する汀線の部分を探索するのに効率が良いのです．

また実装の際には母点の座標と準線の位置さえあれば放物線も Breakpoint も計算できるのでそれさえ持っておけば OK です．
辺は始点だろうが終点だろうがどのみち Circle Event の時の円の中心になるのでその時決定すれば良いです．

詳しくは Haskell で実装してる人がいるのでそちらでどうぞ．

github.com

また，他の言語での実装もあります．(英語だけど)

Fortune’s algorithm and implementation

おわり

日本語の資料少ないので苦労しました…

http://www.cs.sfu.ca/~binay/813.2011/Fortune.pdf
http://cs.nyu.edu/yap/bks/egc/09/7Vor.pdf

ちなみに僕の Visualizer のソースコードは以下にあります．Haskell の Haste というライブラリを使って JavaScript のコードを生成しています．
したがって描画は HTML5 の Canvas です．
僕が確認しただけでも以下のようなバグがあるので，Pull-Request をいただけると助かります…

・最初の母点 0, 1 だけの放物線の場合，片方が表示されずかくかくする (治った？)
・点の配置による原因不明のフリーズ
・放物線の継ぎ目がはみだしてる (CircleEvent が発生していない) ← CircleEvent がおきその先は画面によって切り取られているときにおかしくなる

これアルゴリズムの性質上 Circle Event が起きないと辺が終わらないので画面に切り取られているような辺は自分で計算しなきゃいけないんですよね…
僕が改造したんですが，ボロノイ辺がありえないようなものになる場合があるみたいで…

だいぶ下調べはしたはずなんですけど，間が空いちゃって今は Haskell のソースちらちら見ながら書いてるだけなので，間違い等あったらご指摘お願いします…

転送行列の意味 (1 次元 Ising モデルの厳密解)

こんにちは！Ysmr-Ry です！

最近僕は統計力学を勉強しているのですが，その中の有名な話について紹介します．転送行列についてです．

僕の読んでいる本では，行間が空いており飛躍があるように感じましたので，そのあたりを補完できればと思います．

物体の磁性と Ising モデル

統計力学の応用の中で魅力的なものに，物体の磁性があります．物体の磁性とは物体が磁石にくっつくかということです．
物体の磁性の起源は電子のスピンと言われるものです．
スピンは簡単に言えば電子が上向きか下向きかということで，スピン変数 \sigma = \pm1 で表し， \sigma=+1 ならば上向き， \sigma=-1 ならば下向きということにします．
簡単のため，ここでは 1 次元，つまり L 個の電子が一列に並んでいる状況を考えます．
例えばこんな感じです．

f:id:Ysmr_Ry:20170705204451p:plain

上向き・下向きの矢印が電子を表しています．この場合は L=5 です．
以下周期的境界条件を考え，この場合ならば L=5 の列が無限に繰り返されているという状況を考えます．
すると，左から i 番目のスピン変数を \sigma_i と表すと， \sigma_{L+1} = \sigma_1 となります．
さて，Ising モデルは，これらのスピンのうち，隣り合っているもの同士が相互作用を持つ時を考える簡単なモデルです．
スピン配位 (\sigma_1,\cdots,\sigma_L) に対し，そのエネルギー E_{(\sigma_1,\cdots,\sigma_L)} は以下のように定義されます．

f:id:Ysmr_Ry:20170705205232p:plain

J は交換相互作用定数という物質固有の正の定数です．
第 2 項は磁場 \vec{H} = (0,0,H) が物質にかかった時の各スピンに対応するエネルギーです．
第 1 項が隣り合うスピン同士の相互作用に対応しています．
今，スピン変数は \pm1 の値しか取らないため，その積 \sigma_1\sigma_2 は \sigma_1, \sigma_2 が同じ値ならば 1，違うなら -1 になります．よって， -J\sigma_1\sigma_2 は同じ値ならば -J，違うならば +J になります．
エネルギーは低いほど安定*1ですから，隣り合うスピンが同じ方が安定になるようになっています．このように相互作用を組み込むのが Ising モデルです．

転送行列が出てくるまで

先のエネルギーの式を考えやすいように少し変形しましょう．周期的境界条件が成り立っていますから，

f:id:Ysmr_Ry:20170705210600p:plain

です．
ここで，カノニカル分布を考えます．
カノニカル分布では，エネルギー E_i に対する確率を \frac{e^{-\beta E_i}}{\sum_i e^{-\beta E_i}} とします．
Z(\beta) = \sum_i e^{-\beta E_i} は確率の総和を 1 にするためのつじつま合わせの定数 (規格化定数) で，分配関数と呼ばれています．
カノニカル分布は， U = -\frac{\partial}{\partial\beta} \log Z(\beta),\ F = -\frac{1}{\beta} \log Z(\beta) のように分配関数 Z(\beta) からマクロな熱力学関数をひねり出すので，分配関数が解析的に閉じた形で表せれば勝ちなのです．
この場合の分配関数は，

f:id:Ysmr_Ry:20170705211844p:plain

と変形できます．あとはこの積を解析的に表せば良いのです．
( \prod は \sum の積に対応するものです．)
ここで， \sigma_i, \sigma_{i+1} = \pm1 しかとりませんから， f(\sigma_i,\sigma_{i+1}) のとり得る値は f(\pm1,\pm1) の 4 通りです．
これを 2 次の正方行列で表すと，

f:id:Ysmr_Ry:20170705212503p:plain

となります．この行列 T を転送行列といいます．

転送行列の意味

さて，この行列を使って先ほどの積を計算するのですが，まず， \prod 以降の式に注目してみましょう．

f:id:Ysmr_Ry:20170705213013p:plain

となっています． f の引数に着目すると， (\sigma_1,\sigma_2) \to (\sigma_2,\sigma_3) \to \cdots \to (\sigma_{L-1},\sigma_L) \to (\sigma_L,\sigma_1) と， \sigma_1 \to \sigma_2 \to \cdots \to \sigma_L \to \sigma_1 のような値の列，道順のようなものが決まれば \prod で積をとっている各々の項は決まり，全ての道順， 2^L 通り全てについて総和を取れば良いことがわかります．
例えば， L=5 の時， +1\to-1\to+1\to+1\to-1\to+1 という道順がありえます．
(上の図に対応しています)
この時，対応する項は f(+1,-1)\,f(-1,+1)\,f(+1,+1)\,f(+1,-1)\,f(-1,+1)\,f(+1,+1) です．
また，道順の総数は， \sigma_i\,(i=1,2,\cdots,5) の各々が \pm1 の 2 通りの値を取りうるので，重複順列より， 2^5=32 通りです．
それら各々に対応する項をこの要領で求め，総和を求めるということです．
これをグラフに帰着します．
下図のように +1, -1 の 2 つの頂点を持つグラフを考えます．

f:id:Ysmr_Ry:20170705215510p:plain

有限オートマトンのような感じです．
辺を通ったら，辺の上に書いてある式を掛けると約束します．
すると， +1\to-1\to+1\to+1\to-1\to+1 とこのグラフ上を移動することが， f(+1,-1)\,f(-1,+1)\,f(+1,+1)\,f(+1,-1)\,f(-1,+1)\,f(+1,+1) の項を作ることに対応します．
さらに，転送行列 T がこのグラフの隣接行列になります．
グラフに対する隣接行列 A=(a_{ij})とは，上のグラフで言えば，頂点 1 を 1 番，頂点 -1 を 2 とすると， a_{ij} = (\mbox{頂点 i から 頂点 j に行くコスト (ここでは通ったら掛ける式)}) と定義される行列のことです．
隣接行列には便利な性質があって， A^n を考えると，その (i,j) 成分は頂点 i から頂点 j まで，合計辺を n 本通って行く全ての道順に対応する項の和となります*2．
よって， \sigma_1 \to \sigma_2 \to \cdots \to \sigma_L \to \sigma_1 の道順は出発した頂点と同じ頂点に戻る長さ L の道順 (これを閉路といいます) なので， T^L の (1,1) 成分と (2,2) 成分を足したもの，つまり， \mathrm{Tr}(T^L) を求めればよいのです．
ここまでが転送行列の意味 (この記事のメインテーマ) です．あとは，線形代数の初等的知識をもとに計算をするだけです．

トレースを対角化によって求める

行列 A に対する対角化とは， D=P^{-1}AP が対角行列*3となるような正則な行列 P を見つけることです．
よって， A=PDP^{-1} と表せます．今，転送行列 T は対称行列なので，直交行列 P を用いて T=PDP^{-1} と対角化できます．
また， \mathrm{Tr}(AB) = \mathrm{Tr}(BA) が成り立つので，

f:id:Ysmr_Ry:20170705221848p:plain

となります．ただし，対角行列の表す線形変換は拡大変換であるため， D^L は単に D の対角成分 (固有値といいます) を L 乗するだけで OK です．
これが分配関数になるので，あとは具体的に転送行列の固有値を求め，分配関数から熱力学関数をひねり出せばよいです．
やってみると，

f:id:Ysmr_Ry:20170705224550p:plain

と，磁化率などが不連続点を持たないため，1 次元 Ising モデルでは相転移は起こらないことがわかります．

おわり

どうだったでしょうか…
この話は極めて有名なように思うのですが，そこまで簡単な話ではないと思います．現に Ising の学位論文ではこんなスマートな方法は使ってなくて，分配関数の求め方がかなり複雑になっているそうです．

例によって，コメント・ブクマ・スターとくれると泣きながら飛び上がって喜ぶのでよろしければレスポンスをください．

*1:エネルギーは変化する能力なので，それが低ければ変化しにくく安定．

*2:行列の積の定義を考えるとわかります． A^2 の (i,j) 成分は \sum_{k=1}^N a_{ik}a_{kj} なので，これで総和をとっている各項は頂点 i → 頂点 k → 頂点 j と行く通り数であり，k を全ての頂点について動かしているので，頂点 i → 頂点 j へと行く長さ 2 の全ての道順を試していることになります．以降帰納的に示せます．

*3:対角成分のみがあるような行列， \begin{pmatrix}\lambda_1 & 0 \\ 0 & \lambda_2\end{pmatrix} のような行列のことです．