Zenn
🤖

Union-Findアルゴリズムと計算科学の話

2024/02/10に公開
Python
#
数値計算
tech

概要

Union-Findアルゴリズムは、グラフ系でよく使われるアルゴリズムですが、物性物理においても重要な役割を果たします。その実装と使われ方について少しメモしておきます。

Union-Findアルゴリズム

何か2つの物が与えられたとき、それらがある意味において「同じグループに所属するかどうか」が判定できるとき、その関係を同値関係と呼びます。例えば相似などが典型例で、図形Aと図形Bが相似、図形Bと図形Cが相似ならば、図形Aと図形Cも相似です。このように「友達の友達は友達」みたいな関係が同値関係です。とりあえず同値関係を\simで表現することにしましょう。ある集合の要素a,b,cについて、

  • a\sim a 自分自身とは必ず同じグループ (反射律)
  • a\sim b ならば b \sim a (対称律)
  • a\sim b かつ b \sim c ならば a \sim c (推移律)

が成り立つとき、この二項関係\simを同値関係と呼びます。

さて、ここからアルゴリズムの話です。入力として、二項関係がずらずらと与えられたとします。

\begin{aligned} a \sim b \\ c \sim d \\ b \sim c \\ \cdots \end{aligned}

ここから「ある要素xとyは同じグループか？」「最大のグループのサイズはどれくらいか？」などが知りたいとします。このときに使われるアルゴリズムがUnion-Findアルゴリズムです。

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関数を使います。

同値関係i\sim jが与えられたとき、それを反映させる関数が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アルゴリズムはグラフ操作でよく現れるアルゴリズムですが、物理現象の解析もよくでてきます。典型例はパーコレーションです。例えば正方格子に、ある確率pで粒子があり、1-pで粒子が存在しないとき、上下左右に粒子があるときだけ移動できるとするとき、この系全体にわたって移動可能かどうか、という問題を考えることができます。これは例えば不純物を含む金属のモデルになっており、系全体が移動可能であるときには電気が流れ、そうでない場合は絶縁体になっている、と考えることができます。パーコレーションについてはこちらも参考にしてください。

さて、物性物理において、Union-Findアルゴリズムはスピン系の更新アルゴリズムでよく使われます。単純な例は正方格子強磁性イジング模型です。この模型の自発磁化の温度依存性を知りたい場合、普通はマルコフ連鎖モンテカルロ法を使い、メトロポリス法などでスピンを一つ一つひっくり返してサンプリングをすることになります。

このようなスピン系のサンプリングをする際、効率的なアルゴリズムとしてSwendsen-WangアルゴリズムやWolffアルゴリズムがあります。

これは、分配関数の和を、Fortuin-Kasteleyn表現と呼ばれるグラフの和に取り直し、そのグラフについてサンプリングするというものです。

その基礎アルゴリズムにUnion-Findアルゴリズムが使われます。詳細はSwendsen-Wangアルゴリズムのサンプルコードなどを参照してください。

Swendsen-Wangアルゴリズムは、簡単に書くと

  1. スピンを相互作用と温度から決まる、ある確率でつなぐ
  2. 繋がったスピンをクラスターに分け、クラスターごとに独立にひっくり返す

ということをします。この「スピンを繋いでクラスターに分ける」というところに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

この関数は

  1. 配列を用いた間接参照がある
  2. ループボディが小さく、かつ頻繁に条件分岐がある

という、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で実装した例としては、小村、岡部によるものが有名です。

Y. Komura and Y. Okabe, "CUDA programs for the GPU computing of the Swendsen–Wang multi-cluster spin flip algorithm: 2D and 3D Ising, Potts, and XY models", Comput. Phys. Commun. vol. 185, pp 1038-1034 (2014)

これはHawickらによる、グラフラベリングのGPGPU上での並列実装に基づいています。

K. A. Hawick, A. Leist, and D. P. Playne, "Parallel graph component labelling with GPUs and CUDA", Parallel Comput., 36 (2010), pp. 655-678

グラフラベリングを並列化すると、どこかで衝突がおきてしまうため、それをアトミック演算(atomicMin())を使って排他制御をしているようです。

マルチスレッド実装

僕が知っているUnion-Findの並列実装の物理系への応用は、藤堂らによる量子スピン系で使われるループアルゴリズムのマルチスレッド実装です。

S. Todo, H. Matsuo, and H. Shitara, "Parallel loop cluster quantum Monte Carlo simulation of quantum magnets based on global union-find graph algorithm", Comput. Phys. Commun. vol. 239, pp. 84-93 (2019)

これは、大規模なUnion-Findアルゴリズムの並列化を、x86と「京」で実装したものです。やはり、ラベルの更新に排他制御が必要となるため、「compare-and-swap」と呼ばれる、アトミックな命令を使っています。これはマルチスレッドを実装する現代のCPUには必ず実装されている命令であり、x86ではcmpxchglを、京の石であったSPARKではcasを、インラインアセンブリで呼び出して実装したようです。

これにより、スレッドセーフなUnion-Findアルゴリズムが実装されています。

まとめ

Union-Findアルゴリズムは広く使われているアルゴリズムですが、物性物理の分野でもSwendsen-Wang、Wolff、ループアルゴリズムのベースアルゴリズムとして使われています。これらはアルゴリズムとしては強力でありながら、メモリのレイテンシが直接見える、ランダムアクセスになるためキャッシュ効率が悪いなど、現代のアーキテクチャでは性能を出しづらい形になっており、どのように実装すればよいかは非自明です。また、大規模化しようとすると並列化が必要になりますが、そのためにはスレッドセーフな実装が必要であり、CPUでもGPGPUでもアトミックな命令が必要になります。

僕の知識が京コンピュータあたりで止まっているので、最近ではもっと進展があるのかもしれません。詳しい人は是非フォローの記事を書いてください。

GitHubで編集を提案

Discussion