[目次][前ページ][次ページ][ＯＨＰ][質疑応答]
(3/4)

３．量子モンテカルロ法の並列化


３．１ 量子モンテカルロ法

前節で述べた厳密対角化法では、必要となる計算時間、および主記憶 容量が系のサイズに対して指数的に増加するため、扱うことのできる系のサイズ は非常に限られる。前述の並列化手法を用いたとしても、L=34程度が、現在の スーパーコンピュータでの限界である。さらに大きな系を扱う方法として、広く 用いられている方法の一つに「モンテカルロ法」がある。この方法では、広大は 位相空間の中から、代表的な少数の点をある「重み」に従って確率的にサンプリ ングすることにより、系の物理量を求める。

ここでは、量子反強磁性ハイゼンベルグ模型と呼ばれる量子スピン系を考える。 前節で扱ったイジング模型はスピンのもつ量子性を無視した模型であり、スピン は上向き、下向きの2つの自由度しか取り得なかったが、量子スピン系では、不 確定性原理に起因する量子ゆらぎ(あるいは演算子の非可換性と呼んでもよい)の ため、スピンを「単なる方向をもった矢印」として扱うことはできない。現実の 系においても、低温になればなるほどこの量子ゆらぎの効果が強くなり、シミュ レーションにおいてもこの効果を正しく取り入れた模型を考える必要がある。

この量子的にゆらいだスピン状態を計算機の中で表現するために、一般には「経 路積分(path integral)」と呼ばれる手法が用いられる。すなわち、d次元の量 子系を(d+1)の次元をもつ古典系へとマップする手法である。スピン状態はこ の余分な1次元方向(虚時間方向と呼ぶ)にゆらいだ「世界線(world line)」とし て表わされる(図4a)。この世界線の配置は、計算機の中では、 ダブルリンクリストを用いて、時刻0での状態の情報に加え、世界線がジャンプ した時刻と場所を覚えておくことにより表現される。

量子モンテカルロ法では、この世界線の配置を「重み」に従って確率的に変化さ せていくことにより位相空間のサンプリングを行なう。以前は、世界線を局所的 に変形させて量子モンテカルロを行なっていたが、この方法では、非常に低温、 もしくは相転移の臨界点近傍になると、状態の変化が非常に遅くなり、効率的な サンプリングが不可能になる。この問題は「臨界減速(critical slowing down)」 と呼ばれる。最近、この問題を回避する方法として「ループアルゴリズム」が提 案された。この方法では、世界線の配置の更新は、局所的に行なわれるのではな く、グローバルなオブジェクト(ループと呼ばれる)を反転させることにより行な われる。1モンテカルロステップは、具体的には以下のような手順となる。

(1)
ラベル付け(labeling)
世界線の配置(図4a)の中から、世界線が通っている部 分と通っていない部分が隣りあっている領域を探し出し、そこにある 「密度」で確率的に「リンク」を配置する。また、世界線がジャンプし ている箇所には確率1で「リンク」を配置する(図4b)。
(2)
ループ認識 (loop identification)
上記のラベル付け操作により作成されたグラフの上をある点から辿って いく。途中で「リンク」に出会ったら、隣りにジャンプし、進む方向を 反転する。この操作をループが閉じる(出発点に戻る)まで続ける。次に、 まだループの通っていない点から出発し、同じ操作を繰り返す。全ての 点が必ずどれか一つのループに属するまで繰り返す(図 4c)。
(3)
ループ反転 (loop flip)
得られたそれぞれのループ毎に、そのループ上のスピン状態をある確率 (一般的には 1/2)で反転する。これにより次のステップの世界線の配置 が得られる(図4d)。
この操作を何万回、もしくはそれ以上の回数繰すことにより物理量の測定を行なう。


３．２ ベクトル化

ループアルゴリズムでは、世界線の配置はダブルリンクリストを用い て表現され、このリストに対してリンクの挿入・削除や、リンクリストを辿ると いった操作が主体となっており、残念ながらベクトル化は全く不可能である。実 際、VPP500では、Alpha チップなどのRISC CPU と比較すると数十分の一から数 百分の一程度の性能しか出ない。ループアルゴリズムを用いた実際のシミュレー ションには、Compaq Alphaワークステーション、SGI Origin 2000、および Hitachi SR-2201 等のスカラ(並列)機を使用している。


３．３ 並列化

一般の局所的な状態更新に基づくモンテカルロ法は、(d+1)次元空 間を各プロセッサに分割し、各モンテカルロステップで境界上のスピン状態を交 換することにより、比較的自明に並列化可能である。一方、ループアルゴリズム では、グローバルな操作が必要なので、並列化は非自明となる。我々は以下に挙 げる方針に基づき、並列化を行なった。

(d+1)次元空間のうち、(d次元の)実空間方向は離散的な格子であるが、虚時 間方向に関しては、連続変数となっている。したがって、プロセッサ数によらず 均等な分割が可能となる。また、アルゴリズムが空間次元によらない、というメ リットもある。


３．４ 並列性能評価

上記の並列化アルゴリズムの理論的な並列化効率は以下のように見積 もることができる。

実空間方向の系のサイズをL、虚時間方向のサイズ(これは系の温度の逆数に比 例する)を$\beta$とすると、逐次実行の場合の実行時間は、ステップ(1)、(2)、 (3)共に $\beta L^d$ に比例する。$N_{\rm p}$プロセッサで並列実行した場合、 ステップ(1)、(2-1)、および(3)の所要時間は単純に $1/N_{\rm p}$倍となる。一 方、ステップ(2-2)については、 $O(L^d\log N_{\rm p})$の時間が必要であるの で、並列化によるスピードアップは、cをある定数とすると、

\begin{displaymath}\frac{\beta L^d}{\beta L^d/N_{\rm p} + c L^d\log N_{\rm p}} \... ... \left\{ 1 - c \frac{N_{\rm p} \log N_{\rm p}}{\beta }\right\} \end{displaymath}

すなわち、 $N_{\rm p} \log N_{\rm p}$に比べて$\beta$が十分に大きければ(温度が十分に低ければ)、 高い並列化効率が期待される。

図5に、SGI Origin 2000 および Hitachi SR-2201 での実 測例を示す。特にSGI Origin 2000では、非常に高い並列化効率が得られている。 実際のシミュレーションでは、この並列化したアルゴリズムを用いて、256プロ セッサまでの並列計算を行なっている。


  
図 4: 量子モンテカルロ法(ループアルゴリズム)における1モンテカルロス テップ。水平方向が実空間方向、垂直方向が虚時間方向を表す。虚時間方向に は、周期的境界条件が課されている。
\resizebox{.70\textwidth}{!}{\includegraphics{loop.eps}}


  
図 5: 量子モンテカルロ法(ループアルゴリズム)の並列化効率。 横軸はプロセッサ数、縦軸は速度の向上の度合を示す。
\resizebox{.6\textwidth}{!}{\includegraphics{loop-para.eps}}

[目次][前ページ][次ページ][ＯＨＰ][質疑応答]