7. アルゴリズム

7.1. 経路積分サンプリング

DSQSSでは, 分配関数を

\[Z \equiv \text{Tr} e^{-\beta H} =\sum_S W(S)\]

と経路積分表示したのち, マルコフ過程によって状態 \(S\) を確率的時系列的に発生し, これをサンプリングします(経路積分モンテカルロ法). ここで \(H\) は系を記述するハミルトニアンで, 「状態」とは \(d+1\) 次元時空上で定義された古典変数の場です (モデルが定義されている空間次元を \(d\) として, これに虚数時間軸を加えたものを \(d+1\) 次元時空と呼びます). また c-数 \(W(S)\) は状態 \(S\) のボルツマン重みです. 「状態」をひとつ定めることは, スピンや粒子の空間的な配置が虚数時間の増加に伴って変化する経路をひとつ定めることと等価であるため, 状態 \(S\) は「経路」とも呼ばれます. ボーズ粒子系のように局所的な粒子数保存則が成立する場合に, 状態を粒子の存在位置を実線でつないだ軌跡(=世界線)の集まりとして視覚化することが多いため, 経路積分モンテカルロ法は「世界線モンテカルロ法」とも呼ばれます. マルコフ過程の遷移確率は, 定常分布において状態の出現頻度が重み \(W(S)\) に比例するように定義されます. マルコフ過程で順次出現する状態を \(S_t (t=1, 2, 3....)\) とした時, 演算子 \(Q\) の期待値

\[\langle Q \rangle \equiv \mathrm{Tr}( Q e^{-\beta H} ) /\mathrm{Tr}( e^{-\beta H} )\]

\(Q\) に対応する観測量 \(Q(S)\) の統計的期待値

\[\langle Q \rangle_\text{MC} \equiv \frac{1}{N_\text{MCS}}\sum_{t = N_\text{therm}+1}^{N_\text{MCS}+N_\text{therm}} Q(S_t)\]

によって近似されます. この近似は初期条件の影響による系統誤差と, サンプリングによる統計誤差を含みます. 系統誤差が無視できるためには, 空回し数 \(N_\text{therm}\) が初期緩和時間よりも大きいことが必要で, その限りでは指数関数的に速やかに収束します. 一方統計誤差はサンプル数 \(N_\text{MCS}\) を大きくしたときに, この \(1/2\) 乗に逆比例して小さくなることが期待されます.

7.2. ワーム更新法

マルコフ過程の遷移確率の構成方法, すなわち状態更新方法にはいろいろな種類があり, それぞれ利点・欠点があります. スピン系, ボーズ系のシミュレーションで用いられる代表的な更新方法はループ更新とワーム更新です. ループ更新は全系をループと呼ばれるクラスターに分割してループごとに状態更新するもので, 高速な更新が可能だが, 一様磁場中での反強磁性体やボーズ系では効率が著しく低下してしまいます. 一方ワーム更新は, 全系に保存則を破る点(ワーム)を2つ導入して, これらの点を移動させることで状態を更新する方法です. ワームは非対角成分のみをもつ演算子に対応します. たとえば, 各点での粒子数が対角化されている表示を用いたボーズ系のシミュレーションでは, ワーム演算子として生成消滅演算子を用います. このとき, ワームの前後の虚数時刻では粒子数が1だけ異なっていて, このワームが移動すると, ワームが通過した部分の粒子数は通過前から1だけ変化し, ワームの移動によって状態が更新されていきます. ワーム更新法は, 上記のループ更新に比べると広い範囲のモデルで有効です.

オリジナルのワーム更新法では, ワームは虚時間方向および実空間方向にランダムウォークしていました. 一方, ワームの空間方向での移動が起きる候補であるバーテックスを, ループ更新と同様にしてあらかじめ作り, ワームはバーテックスにぶつかるまで虚時間方向で直進するようにしたものが向き付きループアルゴリズム(DLA)です. DSQSS/DLA はDLA を実装したプログラムです.

7.3. マルチワームアルゴリズム

DLAではただひとつのワームヘッドを動かすことで状態を更新しているため, 時空間分割などの非自明並列が行えません. この問題を解決するために提案されたアルゴリズムが, ワーム対の個数制限を取り除いたマルチワームアルゴリズム (MWA) および その時空間並列版アルゴリズムである PMWA です. DSQSS/PMWA はPMWA を実装したプログラムです.

7.4. on-the-fly バーテックス法

ワーム更新法では, ワームの生成, ワームの時間方向の移動, ワームの空間方向の移動(散乱), ワームの消滅などからなるサイクルが更新の時間的単位となっていて, このサイクルを繰り返すことで逐次的に状態が更新されます. このうちワームの空間方向の移動を実現する仕方として, 移動の起こる場所(バーテックス)をあらかじめ全系に配置しておくやり方(固定バーテックス法)と, ワームの進行方向に必要に応じて配置するやり方(on-the-fly バーテックス法)とがあります. 固定バーテックス法では, 状態更新の基本的単位(1モンテカルロステップ)は,

  1. バーテックスの配置

  2. ワームの生成から消滅までの(複数)サイクル

のふたつのフェーズからなりますが, on-the-flyバーテックス法ではこのうちの2のみを行います.

DSQSS/DLA では, on-the-fly バーテックス法を, DSQSS/PMWA では固定バーテックス法を採用します.

7.5. バーテックス配置

分配関数に現れる重み \(W(S)\) は格子点 \(i, j\) 間の相互作用を \(H_{ij}\) として,

\[W(S) = \prod_{(ij), \tau} \langle S_i (\tau + \Delta \tau) S_j(\tau + \Delta \tau)| e^{-\Delta \tau H_{ij}}| S_i (\tau) S_j(\tau)\rangle\]

と書けます. ここで, 虚数時間 \(\tau\) に関する積は \(\Delta \tau\) を単位として離散化されているとしました (これは説明のためであり, 実際の計算では連続虚数時間で行われます). \(H_{ij}\) で相互作用する実空間の2点 \(i, j\) に関して, 虚数時間区間 \([0, \beta)\) は状態 \(S_i (\tau)\)\(S_j(\tau)\) に不連続変化がある時刻によって有限個の区間に分割されます. 各区間内では, 2点の状態は \(S_i (\tau), S_j (\tau)\) は一定であり, 密度

\[\langle S_i S_j | E_0 - H_{ij} | S_i S_j \rangle\]

でバーテックスが一様ランダムに確率的に配置されます. ここで, エネルギーシフト \(E_0\) はある定数であり, 上記の密度が正である限り任意にとることができます. 状態一定の区間に配置されるこれらのバーテックスに加えて, 状態変化(キンク)のある虚数時刻に対しては, 確率1でその時刻にバーテックスが配置されます. バーテックスはファインマンダイアグラムを用いた分配関数の展開における相互作用グラフに対応します.

7.6. ワームの生成・消滅

ワームの生成にあたっては, まず状態遷移先の候補として, 時空から一様ランダムに時空点 \(i, \tau\) とワーム対を表す非対角演算子対 \(Q, Q'\) を選びます. 次に確率 \(p_\text{create}\) で実際にワーム対を生成します. 消滅はこの逆過程であり, ワーム対が同じ時空点に来た時に一定の確率 \(p_\text{annihilate}\) で消滅します. これらの確率は次の式で表される詳細つり合い条件

\[\frac{\Delta \tau}{\beta N_\text{site} N_Q} \times p_\text{create} = p_\text{annihilate} \times (\eta \Delta \tau)^2 \langle S_i (\tau) | Q_i | S_i'(\tau) \rangle \langle S_i' (\tau) | Q_i' | S_i(\tau) \rangle\]

が成り立つように定められます. ここで, \(N_Q\) は取りうる \(Q\) の数で, \(S'\)\(\langle S | Q | S' \rangle\) が非ゼロの値を取る(普通は \(S, Q\) の組に対して一意に定まる)状態です. また, \(\eta\) はワーム演算子に共役な場の量で, 例えばスピン系でワーム演算子を昇降演算子にとったときは横磁場(の半分)です. DLAでは, \(\eta\) の値は任意にとれるので, DSQSS/DLA の補助ツール dla_alg では \(\eta^{-2} = \beta N_\text{site} N_Q \Delta\tau \max_{S, Q} |\langle S | Q_i |S' \rangle|^2\) とすることで

\[p_{\rm create} = |\langle S_i (\tau) | Q_i | S_i'(\tau) \rangle|^2/\max_{S, Q} |\langle S | Q_i | S' \rangle|^2\]

かつ \(p_\text{annihilate}=1\) としています. algorithm.xml をユーザが独自に編集することで, ユーザの指定した生成消滅確率でシミュレーションを行うことも可能です.

7.7. ワームの散乱

バーテックスにおけるワームの散乱は, バーテックス自体の重みとワームの持つ重みとの間に詳細つり合いが成立するように決定されます. たとえば, 虚数時刻が増加する方向にサイト \(i\) 上を移動してきたワームが虚数時刻 \(\tau\) にあるバーテックスに当たり, サイト \(j\) を虚数時刻が減少する方向に出ていく過程を考えると, このような衝突が選ばれる確率 \(P\) と, 逆向きの衝突が選ばれる確率 \(P'\) との間には

\[\begin{split}& P \times \langle S_i (\tau+0) S_j (\tau+0) | E_0 -H_{ij} |S_i (\tau) S_j (\tau) \rangle \langle S_i (\tau) | Q_i | S_i'\rangle \\ & = P' \times \langle S_i (\tau+0) S_j (\tau+0) | E_0 -H_{ij} |S_i' (\tau) S_j' (\tau) \rangle \langle S_j' (\tau) | Q_j | S_j\rangle\end{split}\]

の関係が満たされる必要があります. ここで \(Q_i, Q_j\) はワームを表す非対角演算子, \(S'_i (\tau), S'_j (\tau)\) はワームが通過した後のそれぞれのサイト, 時刻の状態です. DSQSS/DLAでは, algorithm.xml ファイルをユーザが直接編集することによって, 任意の散乱確率を指定できるほか, ハミルトニアンを指定したときに, このような条件を満たす散乱確率を自動的に計算する補助ツール dla_alg を持っています.