大規模並列化量子モンテカルロ法の開発と	 冷却ボーズ原子系のシミュレーション	 
正 木 晶 子	 
東京大学物性研究所	 	 	  本稿では筆者らが開発した大規模数値計算のための量子モンテカルロ法の新しいアルゴリズムである“並列化マルチワームアルゴリズム”の概要と，FX10 を用いて行ったアルゴリズムの性能評価の結果について報告する。	 	 
� � � � � � 	 	  量子モンテカルロ法には様々なアルゴリズムが存在するが，本研究で用いるファインマンの経路積分表示に基づいた世界線モンテカルロ法 1 は，有限温度下において，与えられた初期条件の影響を受けることなく，統計誤差の範囲内で厳密な計算が可能な優れた手法の一つである。
また，汎用性も高く，負符号問題2 がない範囲であれば様々な物理系に適応可能である。この手法は原子核物理の分野におけるハドロン模型や，物性物理における強相関模型など，幅広い分野で用いられているため，読者も一度は耳にしたことがあるかもしれない。ここ 20 年ほどで量子モンテカルロ法は大きく発展をとげ,様々な分野の量子系の理論的研究に大きな功績を残してきた。中でも世界線モンテカルロ法のための� ー � � ル � � ズ � 3 は最も汎用性が高く効率的なものとして知られている。	 	 	  しかし,現在残されている重要な問題には,既存の量子モンテカルロ法のアルゴリズムを用いてもなお到達できない様な,大規模な系での計算を必要とする問題が多く存在する。例えば対称性の破れを伴うような相転移の臨界現象やそれに類似する弱い一時転移である。系を構成する粒子やスピンが互いに影響を及ぼし合っている範囲（相関長）が転移点に近づくにつれ増大するが，これを覆う大きさの系で計算しなければ，相転移の性質を正確に知ることはできない。
また，液体ヘリウム４では超固体の存在の有無が古くから実験的にも理論的にも研究対象とされてきたが，量子モンテカルロ法でも多くの先行研究が存在するにもかかわらず，未だ説得力ある決定的な答えは示されていない。量子モンテカルロ法での第一原理的な実験の予測が重要視されているが，実験室系を再現できるような規模での精度の良い計算は難度の高い課題となっている。	 	  このような大規模系の問題解決の糸口として考えられるのが，アルゴリズムの並列化である。
TOP500 を見てわかるように，近年のハイパフォーマンス・コンピューティングの傾向は，多量のコアを搭載したスカラー計算機を用いた大規模並列計算である。もはやアルゴリズムの並列化は，数値計算をする者にとっては無視できない課題となっている。	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 1	 M.	 Suzuki,	 Prog.	 Theor.	 Phys.	 56,	 1454	 (1976);	 M.	 Suzuki,	 S.	 Miyashita,	 and	 A.	 Kuroda,	 2	 後述するような世界線の更新確率に負の符号が生じる問題。	 

3	 N.	 Prokof'ev,	 B.	 Svistunov,	 and	 I.	 Tupitsyn,	 Sov.	 Phys.	 JETP	 87,	 310	 (1998)など。	 	 

スーパーコンピューティングニュース- 	 15 -Vol.16,No.62014

	  残念なことに，ワームアルゴリズムは後述する理由から，並列化アルゴリズムは非自明であり実現されてこなかった。そこで我々は大規模並列計算が可能なワームアルゴリズムの開発に取り組んだ。完成したアルゴリズムを用いて，FX10 上でフラット MPI による 3200 コア並列のベンチマーク計算を行ったものを本稿にて紹介する。	 	 

� � � � � � � � � � � � � � � � � � � � � � � � � � � 	 	  本章では，格子模型のための世界線量子モンテカルロ法について説明する。我々が開発した並列化マルチワームアルゴリズムの概要と並列化の方法について紹介するために，まずは最初の節で，量子モンテカルロ法の基本的な考え方やアルゴリズムについて簡単に解説する。	 	 
� �� � � � � � � � � � � � � � 	 	  典型的な世界線モンテカルロ法では，ファインマンの経路積分に基づいて描かれた d+1 次元空間の世界線配置をサンプリングする。具体的には，逆温度βを Nτ 分割して Δτ = β Nτ とし，鈴木−トロッター分解を用いて，分配関数を，	 

Z = Tr e− β H = ∑ i e− β H iiNτ= lim ∏ ∑ ik e−Δτ H ik+1Nτ →∞k=1 {ik }	 と経路積分表示する。ただし H は系のハミルトニアンを， ik は完全直交系をそれぞれ表している。一行目から二行目の過程で， O(Δτ 2 ) 以上の項を無視し，また，Nτ 個の完全系を挿入した。しかし，このままでは各局所状態 ik e−Δτ H ik+1 において全ての相互作用を考慮する必要がある。
Nb

ここではハミルトニアンを， H = ∑ H と局所的なハミルトニアンの和で書ける場合に限定してbb=1考える。例えば b は相互作用のあるボンドやサイトを表しているが，本稿では b はボンドを表すこととする。すると，	 NτNb

Z = ∑ ∏ ∏ ik,b e{ik } k=1 b=1−Δτ H b= ∑W ({ik }){ik }ik+1,b	 	 	 	 (１)	 となり，各局所状態を２体問題として扱うことが可能となる。式（１）の W(…)の部分をグラフィカルに表現すると，図１の様になる。黒の太線が状態 ik の虚時間発展を追った世界線である。式（１）の和を実行するには，すべての世界線配置を再現すればよいが，そのまま計算するには膨大な状態を発生させなければならなくなる。そこで実際には，重み付きサンプリング（重要な寄与を持つ状態のみを発生させてサンプリングを行う）の精神に則り，エルゴード性と詳細釣り合いを満たすよう定義された確率過程のもとモンテカルロ更新を行うことで配置	 

スーパーコンピューティングニュース- 16 	 -

Vol.16,No.6

2014

をサンプリングする。シミュレーションの過程で，各配位空間での物理量をサンプルすることで，統計誤差の範囲内で厳密に物理量期待値を計算する。ただし，図１では，Nτ=4 の場合を図示したが,実際の計算では Nτ→∞の極限をとっている。	 	 

図１：鈴木−トロッター分解による世界線の配位空間の図。
分かりやすいように実空間次元が１次元の場合を示した。縦軸は虚時間で 0 とβで周期境界条件になっている。横軸は実空間軸であり，こちらも周期境界条件が課されている。世界線はグレーのボックスでのみ斜めに進める。

� �� ． � � � � � � � � � � � � � � � � 	 	  ワームアルゴリズムは，世界線モンテカルロ法のための状態更新アルゴリズムである。ワームアルゴリズムではΔτ→０の連続虚時間極限で配位空間を表現できる。配位空間の更新は，
� � � と呼ばれる非対角演算子の対を挿入した，「拡張配位空間」中で行われる。具体的にどのような演算子なのかは扱う基底によって異なるが，世界線の不連続点を生成するように定義される。一対のうちの片方のワームが釣り合い条件を満たす確率過程のもと配位空間中を動きながら世界線配置を更新していき，自身の片割れと再び出会ったときに切れていた世界線がつながり，新しい世界線配置が完成する仕組みとなっている（図２）。この演算子の対はグリーン関数に相当するものであり，よって，ワーム間には相関があり，この伸びた相関の範囲内での大域的な更新を行うのがワームアルゴリズムの特徴である。	 	  ワームアルゴリズムは，汎用性も高く最も優れたアルゴリズムの一つであるが，アルゴリズ
ムの並列化は困難な仕様となっており，これまで成し遂げられてこなかった。一つ注意しておきたいことは，ここで言う並列化とは，統計精度を上げるための単純な乱数並列計算（自明な並列化）のことではなく，１つの状態空間に対して施す�� � � � � � � のことである。我々は系をドメイン分割し,各コアに割り当てることを考える。	 	  これまでワームアルゴリズムの並列化が成し遂げられなかったのは，状態更新を行うワームが空間中の一点にのみ挿入されているためであると考えられる。このため，仮に配位空間をドメイン分割しても，ワームが挿入されているドメイン以外は更新されず並列化の意味をなさなくなる。	 	  それでは，ドメインごとにワームを配置できれば解決するのか，と言えば全くそうではない。	 

スーパーコンピューティングニュース- 	 17 -Vol.16,No.6

2014従来のアルゴリズムと同様な仕様では，ワームがすべて消えるまで待つことになり，並列化効率が実行途中から著しく低下し，やはり並列化の意味をなさなくなる。しかし，それ以前に，アルゴリズムの根本にあるイベント駆動性が悪さをする。すなわち，ワームがドメイン境界に達するたびに通信を行えば，かなり不規則かつ頻繁に通信が発生してしまい，非常に非効率的である。そもそも，別の領域で発生する予定の更新の終状態を予測することは，ワームの動きをトレースしない限り困難であるため，移動するタイミングや移動のための更新確率を決めるのも困難である。	 	  並列化の際には，これらの問題を解決するようにアルゴリズムを設計することはもちろんのこと，重み付きサンプリングが正しく実行されるように注意する必要もある。	 

図２：一対のワームが世界線の配置を更新する過程の模式図。
太線が世界線。三角と丸がワーム対であり，シミュレーション上では三角が動き回り（ワームのあたまと呼ばれる），丸が止まっている（ワームのしっぽと呼ばれる）とする。最左図でワームが挿入され，中２つの図であたまが移動しながら状態を更新する。最右図でしっぽと出会い対消滅する。

� �３ ． � � � � � � � � � � � � � � � � � � 	 	  我々が開発した並列化マルチワームアルゴリズムの特徴は，複数のワームを配位空間中に挿入でき，ドメイン分割による並列計算が可能であるところである。図３−１に模式図を示す。複数ワームの導入のために次のような，求めたい系のハミルトニアン H にソース場 Q を付加した
拡張ハミルトニアン Hη を定義する。	 

Hη = H + ηQ 	 	  ηは任意のパラメータで，これを操作することでワームの数を調節できる。従来のアルゴリ
ズムとの大きな違いは，ワームを対消滅させず，η→0 への外挿によって求めたいハミルトニアンでの物理量を求める点である。	 	  更新過程については，ワームの挿入方法が新しく，隣接する世界線同士の相互作用項（��� � � � ＝図３−１の説明文参照）で区切られたセグメントごとに「� � � � � � 」を	 定義し（図３−２の説明文参照），セグメントごとに独立に挿入されるワーム数とその位	 置が決定される。この「ワームの重み」の概念は，ドメイン並列の際に威力を発揮する。各ドメインはフラット MPI で各コアに割り当てられ，ワームのドメイン間移動は，ワームの重みをMPI 通信によりやり取りすることで間接的に実行できる。これは２−２節で挙げた問題点を克服できている。これについては次の節で詳しく説明する。もしワームのドメイン間移動が行われないと，ドメイン境界上の状態が更新されなくなる（エルゴード性が破れる）。また,正しく移	 

スーパーコンピューティングニュース- 18 	 -

Vol.16,No.62014

動できてないと，ワーム間の相関が一つのドメイン内にとどまってしまい，本来の系の大きさで計算ができなくなってしまうので注意が必要である。	 	  また,ここでは詳しく述べないが，複数のワームを挿入した場合，ワームが従来のアルゴリズムと同じ動きをするのでは詳細釣り合いが満足されなくなる。我々のアルゴリズムでは，この点も克服するようなワームの動きを考えた。	 	 

図３−１：並列化マルチワームアルゴリズムの配位空間とドメイン分割の模式図。
左の図は配位空間をドメイン分割したもの。右の図（吹き出し内）は一つのドメインに着目して詳しく示したもので，横線はバーテックス，三角はワームを表しており，×は隣のドメインに所属する格子との相互作用のバーテックスを表している。	 

� �� ． � � � � � � � � � � � � � � � � � ��� � � 	 	  ワームのドメイン間移動は，ワームが境界に達した時に行うのではなく，一定プロセスごとにワームの重みをドメイン間通信することで実現されるため，通信コストはそれほど大きくならない。図３−２はワームの重みのドメイン間通信を模式的に示した図である。	 	 ド メ イ ン 境 界 の 状 態 と ワ ー ム を 更 新 す る に は ， 最 初 に 図 ３ − ２ の 左 図 の 矢 印 の 方 向 へ
MPI_Sendrecv	 を実行し，Receive した情報からドメイン境界を更新し，先ほどの矢印と逆方向へ MPI_Sendrecv を再度実行することで更新された状態を隣のドメインと共有する。	 	  ドメイン表面に隣接するセグメントの情報とドメインをまたぐ相互作用の情報のみ通信すればよいため，比較的少ない通信量で済む。我々のコードは C++で設計されており，セグメントの情報とは，そのセグメントに所属するバーテックスのクラスの型とワームの重みであるdouble 型の情報のことである。バーテックスクラスは，クラスの型のまま送信するのではなく，必要なメンバ変数だけをデータ型ごとに配列に格納し直し MPI_Sendrecv を行う方が，全体の計算時間が速くなった。	 	 

スーパーコンピューティングニュース- 	 19 -Vol.16,No.62014図３−２：図３−１の配位空間の並列計算における通信の模式図。
左の図で，出て行く矢印の先に send し，入ってくる矢印のもとから receive するように MPI_Sendrecv を実行する。矢印の右の図（吹き出し）は一つのドメインに着目して詳しく示したもので，四角は横線で区切られたセグメント上に割り当てられたワームの重みを表している。白抜きの四角は通信が必要な箇所で，塗りつぶしの四角はそうでないものである。矢印の配列に情報を格納して	 １回目の MPI_Sendrecv する。矢印がついてない白抜きの四角は，２回目の MPI_Sendrecv で通信を行う。	 	 

３ ． � � � � � � � � � � � � � � � � � � � � 	 
３ �� ． � � � � � ・ � � � � � � � � � � � � ・ � � � � � � � � � � � � � � 	 	  ベンチマーク計算では，東大情報基盤センターFX10 を用いて，次のような２重占有を禁止した２次元正方格子上のハードコア・ボーズ粒子模型の計算を行った。	 

"%µH = ∑$t ( bi†b j + bi b†j ) +Vbi†bi b†j b j + ∑ ( bi†bi + b†j b j )' 	 
2 i&(i, j ) #	  ここで， bi† (bi ) は格子点 i でのハードコア・ボーズ粒子の生成（消滅）演算子である。(i,j)
は隣接サイトの対（ボンド）を意味する。t はホッピング，V は隣接粒子間相互作用，μは化学ポテンシャルである。このモデルでもワームのソース場は	 Q = ∑ ( bi† + bi ) 	 iである。	 	 スーパーコンピューティングニュース- 20 	 -

Vol.16,No.62014

図４：ボーズ・アインシュタイン凝縮秩序変数のη依存性。
パラメータはβ=16,	 V=3.0,	 μ=1.2 の超流動相を選んだ。L=8〜10240 のサイズでの値をそれぞれlog-log プロットした。ただしインセットは線形プロットである。L＝∞の外挿値は L＝8〜1024 の値から導いた。点線は L=∞をべきでフィッティングしたものである。参考文献にある筆者らの論文参照。	 	  この模型では低温で，超流動相4 が安定化されるが,この相内で有限温度下でのボーズ・アインシュタイン凝縮秩序変数<Q>/L2 のη依存性を調べたのが図４である。なおすべての計算結果
では，ホッピングを t=1.0 にとっている。	 	  各点はηについての外挿から求めたものである。系の一方向の格子点数を L とした。我々は3200 のドメイン分割及び MPI 並列により，従来のアルゴリズムでは到達不可能なサイズである10240×10240 格子点の計算を達成した。結果を無限系（L による外挿値）と比較したところ，誤差の範囲内で一致することが確認できた。これは熱力学極限に対応するような系のサイズを直接計算できたと言える。<Q>は，従来のアルゴリズムでは直接計算できないものであり，２次元系で絶対零度ではηに関して定数となる。また，有限温度相転移がおこる臨界領域では<Q>のηに関するべきの指数から相転移現象の臨界指数を知ることができる。図４はこのどちらでもなく，ηはべき的振る舞いをするが，η〜0.01 と非常に小さい指数をとる結果となった。	 	 
� �� � � � � � � � � � � � � � � 	 	  アルゴリズムの緩和効率は非常に重要である。いくら大きな計算ができても，緩和に途方もない時間がかかるようでは意味がない。前節で述べたように，このアルゴリズムでは，並列化効率を上げる目的から，ドメイン間通信の回数を制限している。通信がない間，ワームは境界にはじかれて動き続けるため，ドメイン境界が増えるにしたがってその影響は増すことになる。
もちろん通信時には境界は更新されるためエルゴード性と詳細釣り合いは満たしているが，更	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 有限温度の２次元系では真の意味での超流動長距離秩序は存在しないため，正確に言うと，準長距離秩序を持つコスタリッツ−サウレス相である。そのため，図４で示した秩序変数がη→0 でべき減衰している。	 
4	 

スーパーコンピューティングニュース- 	 21 -Vol.16,No.62014新効率に影響を及ぼすことが懸念される。この様子を Strong	 scaling で確認したのが図５である。まずサンプル数を固定し，標準誤差の並列数 N に対する変化を調べた（図５−１）。結果，標準誤差は分割数に応じて大きくなっているが，それは非常に緩やか（〜N0.09）であるため，緩和に大きな影響はないことが分かった。N ごとに適当な虚時間方向の分割数 Nβと実空間方向
Nr の分割数を割り当て，N=Nβ×Nr となるようにしている。このアルゴリズムでは虚時間方向と実空間方向とではドメイン境界の更新方法が，本質的には同じであるが，詳細は異なるため，結果の図では若干のがたつきが生じているが，誤差の範囲内である。	 	  さらに，計算のウォール・クロック時間を固定したときの標準誤差の N 依存性も調べ，結果を従来のワームアルゴリズムの一種である向き付きループアルゴリズム（DLA）5 と比較した(図５−２)。我々のアルゴリズムは DLA の考え方をベースに取り入れているため比較対象として選んだ。これにより，N が 8 以上で DLA の精度を上回ることを確認できた。また,標準誤差の大きさの変化は実質 100 パーセントの効率で小さくなっていくのが確認できた。これは Strong	 scaling が非常に効率よく働いているためであり，並列数を稼ぐことにより，従来のアルゴリズムを超える精度の向上を期待できる。	 

図５−１：サンプル数を固定したときの標準

図５−２：実時間を固定したときの標準誤差σ

誤差σの分割数（並列数）依存性。の分割数（並列数）依存性。

L=128,	 β=128,	 V=3.0,	 μ=1.2,	 η=0.002 とした。

L=256,	 β=16,	 V=3.0,	 μ=1.2,	 η=0.004 とした。

2

<Q>/L とエネルギーの標準誤差で，統計平均を求観測実時間を 10000 秒に固定したときのエネルギーめる際のサンプル数 10000 に固定した。参考文献の標準誤差で，向き付きループアルゴリズム（DLA）にある筆者らの論文参照。の結果と比較している。	 
４ ． � � � 価 	 ４ -１ ． ����	 ������� � � � � � � � � 	 	  Weak	 scaling による並列化効率についても検証した。加速効率を図６に示す。Weak	 scalingの加速効率を以下のように定義した。	 

A(N node ) =T1

TNnode× N node 	 

ここで Nnode はノード数である。TNnode は計算に Nnode 用いた時のウォール・クロック時間である。
１ノード用いたときのウォール・クロック時間を基準に取ってある。μは超流動相と固体相を	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 5	 O.	 F.	 Syljuasen	 and	 A.	 W.	 Sandvik,	 Phys.	 Rev.	 E	 66,	 046701	 (2002).	 	 

スーパーコンピューティングニュース- 22 	 -Vol.16,No.62014

とる値をそれぞれ用いた。超流動相と固体相の結果はほぼ一致しており，特に違いは見られなかった。２つの相内ではワームの数が根本的に異なるが，本アルゴリズムでは，ワームの数を積分した重みを用いて通信を行うため，並列化効率には影響がなかったものと推測できる。肝
心のコア数依存性であるが，1024 コア（16 ノード）までほぼ１００パーセントの効率を実現し,その後徐々に衰退している様子が見られたが，約３パーセント程度の減速にとどまった。今回調べた 214 ノード並列の範囲では並列化効率の飽和の振る舞いは見られなかった。	 

図６：Weak	 scaling での加速効率。	 
丸は超流動相でμ=-0.2,	 四角は固体相でμ=5.2。ドメインのサイズは実空間方向 LD=16,虚時間方向βD=8とし,	 V=3.0,η=0.006 とした。黒の点線は 100 パーセントの効率のときの加速効率である。	 
４ -２ ． � � � � � � � � � � 	 	  最後に，fapp によるプロファイリングを行ったので，結果を図７に示す。パラメータは
L=256,β=16,μ=1.2,V=3.0,η=0.008 とし，Strong	 scaling による計算を行った。Performance	 monitor で命令処理回数（MIPS）に対するピーク性能を確認したところ，約 12〜18 パーセントあり,並列数を上げるにつれてピーク性能も高くなると分かり，比較的よい結果であることが分かった。	 

図７：命令処理(MIPS)に対するピーク性能とロードインバランスの分割数（並列数）依存性 。	 	 

スーパーコンピューティングニュース- 	 23 -Vol.16,No.62014	 	 一方，ロードバランスは並列数を増やすに従いインバランスが大きくなっていることが分かった。これはアルゴリズムの性質上，Strong	 scaling で同じηを用いると，系全体のワームの数は同程度でも，ドメインあたりのワーム数は減っていき，ドメインごとのワームの数にばらつきが目立つようになってくる。そうすると待ちのノードも増えることになり，インバランスが増すと推測される。これを解決するにはプログラムのチューニングよりも，ηのチューニング（ワーム数の調整）が必要になると言える。ワームの数が一桁にならず，相が変化しない様なηを選ぶことが好ましいだろう。	 	 

５ ． � � � � � � � � � 	 	  複数ワームの導入によって，並列化可能なワームアルゴリズムを開発し，熱力学極限に対応する規模のハードコア・ボーズ粒子系のシミュレーションに成功した。今回のベンチマーク計
算では単純な模型を用いたが，このアルゴリズムは原理的には負符号問題が出ない範囲において一般的なスピン模型やソフトコア・ボーズ粒子模型に適応することが可能である。我々のアルゴリズムは Weak	 scaling,	 Strong	 scaling ともに非常に高い並列化効率を持つが，物理量測定のプロセス中で統計平均をとる際などで，特定のランクのみで処理を行う箇所が存在しているため，そこでは Open	 MP とのハイブリッド並列を取り入れるなどすれば更なる効率向上が期待できる。	 	 

� � 	 	  若手・女性利用支援推薦制度を通じてお世話になりました，東京大学情報基盤センターの関係者の皆様に深く感謝致します。また,本研究の共同研究者ならびに計算物質科学イニシアティ
ブ第一部会・川島・藤堂グループの皆様には，日頃から有益な議論とご教授を頂きましたこと，感謝申し上げます。本研究は，文部科学省の HPCI 戦略プログラム，および，計算物質科学イニシアティブの助成を受けており,計算の一部は，計算物質科学イニシアティブが確保した FX１０の利用枠を使用しています。また，計算物質科学イニシアティブのアプリ高速化支援サービスを通じてお世話になりました富士通の石附さんに感謝申し上げます。	 	 

参 	  考 	  文 	  献 	 
N.	 Kawashima	 and	 K.	 Harada,	 J.	 Phys.	 Soc.	 Jpn.	 73,	 1379(2004).	 
A.	 Masaki-Kato,	 T.	 Suzuki,	 K.	 Harada,	 S.	 Todo,	 N.	 Kawashima,	 Phys.	 Rev.	 Lett.,	 1 12,	 140603	 
(2014).	 	 

スーパーコンピューティングニュース- 24 	 -

Vol.16,No.6

2014