Union-Findアルゴリズムと計算科学の話
概要
Union-Findアルゴリズムは、グラフ系でよく使われるアルゴリズムですが、物性物理においても重要な役割を果たします。その実装と使われ方について少しメモしておきます。
Union-Findアルゴリズム
何か2つの物が与えられたとき、それらがある意味において「同じグループに所属するかどうか」が判定できるとき、その関係を同値関係と呼びます。例えば相似などが典型例で、図形Aと図形Bが相似、図形Bと図形Cが相似ならば、図形Aと図形Cも相似です。このように「友達の友達は友達」みたいな関係が同値関係です。とりあえず同値関係を
-
自分自身とは必ず同じグループ (反射律)a\sim a -
ならばa\sim b (対称律)b \sim a -
かつa\sim b ならばb \sim c (推移律)a \sim c
が成り立つとき、この二項関係
さて、ここからアルゴリズムの話です。入力として、二項関係がずらずらと与えられたとします。
ここから「ある要素
Union-Findアルゴリズムの実装はいくつかありますが、もっとも簡単なのは一次元配列を使った実装です。以下には概略のみ示しますが、実装のイメージについては例えばUnion-Find木のサンプルコードを参照してください。
全要素数がN
である場合、cluster[i] = i
であるようなサイズN
の配列を用意します。
cluster = [i for i in range(N)]
cluster[i]
は、要素i
の直接の親を指します。自分の親が自分自身、すなわちcluster[i]=i
である時、その要素は同じグループの代表ということにします。このグループのことをクラスターと呼びます。最初は全ての要素が、自分自身だけを含むクラスターの代表です。
さて、Union-Findアルゴリズムは、その名の通りunion
関数とfind
関数を使います。
同値関係union
です。こんな感じに実装できます。
def union(i, j, cluster):
if i < j:
cluster[j] = i
else:
cluster[i] = j
要するに、若い方の番号を親にする、という意味です。
これを全ての同値関係に適用すると、木構造ができます。ただし、cluster
には、例えばCの親はB、Bの親はAなどというバラバラの情報が格納されており、ここからCはAのグループに属す、という情報を引き出す必要があります。それがfind
関数です。
def find(i, cluster):
while i != cluster[i]:
i = cluster[i]
return i
これは、要素番号i
の親をたどっていって、グループの代表までたどり着いたらその番号を返す、というだけです。
これで、例えば任意のノードi
とj
が同じグループに所属するかどうかは
def isSameGroup(i, j, cluster):
i = find(i, cluster)
j = find(j, cluster)
return i==j
として実装できます。他にもクラスターのサイズ分布や、ネットワークの連結判定なんかも簡単にできます。
なお、ここで書いたUnion-Findアルゴリズムはもっとも基本的な実装で、より効率的な実装がいくつも提案されています。
物理との関わり
Union-Findアルゴリズムはグラフ操作でよく現れるアルゴリズムですが、物理現象の解析もよくでてきます。典型例はパーコレーションです。例えば正方格子に、ある確率
さて、物性物理において、Union-Findアルゴリズムはスピン系の更新アルゴリズムでよく使われます。単純な例は正方格子強磁性イジング模型です。この模型の自発磁化の温度依存性を知りたい場合、普通はマルコフ連鎖モンテカルロ法を使い、メトロポリス法などでスピンを一つ一つひっくり返してサンプリングをすることになります。
このようなスピン系のサンプリングをする際、効率的なアルゴリズムとしてSwendsen-WangアルゴリズムやWolffアルゴリズムがあります。
- R. H. Swendsen, J.-S. Wang, “Nonuniversal critical dynamics in Monte Carlo simulations”. Phys. Rev. Lett. vol. 58, pp. 86–88.
- U. Wolff, "Collective Monte Carlo Updating for Spin Systems", Phys. Rev. Lett. vol. 62, pp. 361 (1989)
これは、分配関数の和を、Fortuin-Kasteleyn表現と呼ばれるグラフの和に取り直し、そのグラフについてサンプリングするというものです。
その基礎アルゴリズムにUnion-Findアルゴリズムが使われます。詳細はSwendsen-Wangアルゴリズムのサンプルコードなどを参照してください。
Swendsen-Wangアルゴリズムは、簡単に書くと
- スピンを相互作用と温度から決まる、ある確率でつなぐ
- 繋がったスピンをクラスターに分け、クラスターごとに独立にひっくり返す
ということをします。この「スピンを繋いでクラスターに分ける」というところにUnion-Findが使われます。Swendsen-Wangは非常に強力なアルゴリズムで、たんに「一度に複数のスピンがひっくり返せる」から効率的なのではなく、本質的に臨界緩和の問題を改善します。したがって、Swendsen-Wangが使える場合、それを使わないでシングルスピンフリップでがんばるという選択はあり得ません。
このアルゴリズムをXYやハイゼンベルグスピンなど、連続スピンに適用したのがWolffアルゴリズムで、これも強力です。
これらは古典スピン系ですが、量子スピン系のサンプリング手法としてループアルゴリズムというものがあります。こちらも基礎はUnion-Findアルゴリズムになっており、一つ一つ状態を更新する方法に比べて非常に効率的なサンプリングができるようになります。
計算機科学から見たUnion-Findアルゴリズム
さて、Union-Findアルゴリズムは単純なわりに強力なアルゴリズムなのですが、その実装はあまり現代の計算機と相性がよくありません。
もう一度find
アルゴリズムを見てみましょう。
def find(i, cluster):
while i != cluster[i]:
i = cluster[i]
return i
この関数は
- 配列を用いた間接参照がある
- ループボディが小さく、かつ頻繁に条件分岐がある
という、CPUからすると性能がとても出しづらい形になっています。特に、間接参照が連続するため、メモリアクセスが本質的にランダムアクセスになるため、キャッシュ効率が悪く、メモリのレイテンシが直接見える形になってしまいます。
また、union
も、先程は
def union(i, j, cluster):
if i < j:
cluster[j] = i
else:
cluster[i] = j
と書きましたが、これでは出来上がるUnion-Find木の深さが深くなってしまうため、例えば
def union(i, j, cluster):
i = find(i, cluster)
j = find(j, cluster)
if i < j:
cluster[j] = i
else:
cluster[i] = j
などと、最初にクラスター代表番号を解決しておきたくなります。そうすると、2つをつなぐたびにfind
が走るため、メモリをランダムアクセスすることになります。他にもいろいろ工夫があるのですが、それがいずれもCPU的には負担が大きい形になります。
さらに問題なのは、このコードの並列化が難しい点です。例えば与えられた同値関係のデータをいくつかにわけて、それぞれ独立にUnion-Findアルゴリズムを走らせたとしても、どこかで同じ配列の同じ場所を触る可能性がでてきます。よく領域分割系の並列化と異なり「どのスレッドとどのスレッドがいつ同じ場所を触るか」は事前にわからないため、マルチスレッドで並列化をしようとすると、本質的に排他制御が必要になります。
Union-Findアルゴリズムの並列化
さて、様々な場所で使われるUnion-Findですが、計算機が本質的に並列化前提で進化しているため、当然Union-Findアルゴリズムも並列化したい、というニーズが出てきます。
GPGPU実装
例えばSwendsen-WangアルゴリズムをGPGPUで実装した例としては、小村、岡部によるものが有名です。
これはHawickらによる、グラフラベリングのGPGPU上での並列実装に基づいています。
グラフラベリングを並列化すると、どこかで衝突がおきてしまうため、それをアトミック演算(atomicMin()
)を使って排他制御をしているようです。
マルチスレッド実装
僕が知っているUnion-Findの並列実装の物理系への応用は、藤堂らによる量子スピン系で使われるループアルゴリズムのマルチスレッド実装です。
これは、大規模なUnion-Findアルゴリズムの並列化を、x86と「京」で実装したものです。やはり、ラベルの更新に排他制御が必要となるため、「compare-and-swap」と呼ばれる、アトミックな命令を使っています。これはマルチスレッドを実装する現代のCPUには必ず実装されている命令であり、x86ではcmpxchgl
を、京の石であったSPARKではcas
を、インラインアセンブリで呼び出して実装したようです。
これにより、スレッドセーフなUnion-Findアルゴリズムが実装されています。
まとめ
Union-Findアルゴリズムは広く使われているアルゴリズムですが、物性物理の分野でもSwendsen-Wang、Wolff、ループアルゴリズムのベースアルゴリズムとして使われています。これらはアルゴリズムとしては強力でありながら、メモリのレイテンシが直接見える、ランダムアクセスになるためキャッシュ効率が悪いなど、現代のアーキテクチャでは性能を出しづらい形になっており、どのように実装すればよいかは非自明です。また、大規模化しようとすると並列化が必要になりますが、そのためにはスレッドセーフな実装が必要であり、CPUでもGPGPUでもアトミックな命令が必要になります。
僕の知識が京コンピュータあたりで止まっているので、最近ではもっと進展があるのかもしれません。詳しい人は是非フォローの記事を書いてください。
Discussion