格子場の理論における精度や乱数
あまりまとまりのない話です石川健一

2019/3/18（広島大学理学研究科）

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

1

もくじ
1. 格子場の理論での計算精度
2. 格子場の理論での乱数
3. まとめ2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

2

１．格子場の理論での計算精度• ざっくりいうと𝑥• 格子場の理論＝多重積分• 力学変数＝４次元時空中の場• 連続時空を離散化⇒格子時空• 場⇒配列𝑛𝜙 𝑥, 𝑦, 𝑧, 𝜏𝑥 = 𝑎𝑛𝑥 , 𝑥 ∈ 𝑅, 𝑛𝑥 ∈ 𝑍
𝜙 𝑥, 𝑦, 𝑧, 𝜏 ⇒ 𝜙(𝑛𝑥 , 𝑛𝑦 , 𝑛𝑧 , 𝑛𝜏 )• 力学：運動方程式（作用積分 𝑆[𝑔, 𝜙]）⇒差分方程式（作用積分 𝑆𝐿 [𝑔, 𝜙]）• 経路積分⇒多重積分

𝑂 =1𝑍𝑥𝑑𝜙(𝑥) 𝑂[𝜙]exp[−𝑆[𝑔, 𝜙]]

𝑂𝐿 =1𝑍𝑛𝑑𝜙 𝑛 𝑂𝐿 [𝜙]exp[−𝑆𝐿 [𝑔, 𝜙]]• 連続極限＝理論パラメータ空間における２次相転移点＋くりこみ群の解析
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

3

１．格子場の理論での計算精度• 計算精度絡みでの困難の原因
一．多重積分の自由度（次元）がとても大きいことに起因する精度上の困難
一．運動方程式を解く際に生じる困難。特に差分方程式＝大規模連立方程式の係数行列の条件数がとても大きいときの困難が考えられます。

2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

4

１．格子場の理論での計算精度
一．多重積分の自由度（次元）がとても大きいことに起因する精度上の困難• 示量変数（大きさの絶対値は自由度に比例）の差を計算するときに桁落ち• 2564 格子での各格子点からの何らかの量の総和から 大
きさ１程度の量を差から計算する（例）メトロポリス法に使う重みとなるエネルギー変化の計算
2564 = 4294967296 : すでに10桁の数値。倍精度で１０桁同士の差から大きさ１程度の量を計算すると、５桁から６桁程度しか有効数字がない。
時空の領域ごとに差を計算するなどの工夫が必要。時空の領域が計算機を移動していく場合など領域トレースが困難な場合は４倍精度で足し上げて引く。
ほかの分野（流体、構造解析）ではいかがでしょうか？
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

5

１．格子場の理論での計算精度
一．多重積分の自由度（次元）がとても大きいことに起因する精度上の困難• 手前味噌な話• 数値確率過程摂動論を用いた摂動計算• 結合定数による展開∞∞

𝑔 𝑗 𝜙 𝑛 (𝑗)𝜙(𝑛) =𝑔 𝑗 𝜓 𝑛 (𝑗)𝜓(𝑛) =

𝑗=0∞

𝑗=0

𝑗=0

𝜒(𝑛) = 𝜙 𝑛 𝜓(𝑛)

𝜒 𝑛(𝑗)= 𝜙 𝑛 ⋆𝜓 𝑛(𝑗)

多項式積

𝑗=𝜙 𝑛

𝑔 𝑗 𝜒 𝑛 (𝑗)𝜒(𝑛) =(𝑘)𝜒 𝑛(𝑗−𝑘)

高次打切

畳み込み積
FFTで高速化

𝑘=0• 展開係数が確率変数

𝜙 𝑛 (𝑘)• 物理量は確率変数の関数。量子的な期待値は統計平均で計算する。
𝑁𝑚𝑎𝑥∞

𝑔 𝑗 𝑂 𝜙 (𝑗)𝑂=𝑂 𝜙 =𝑗=0

𝑔 𝑗 〈𝑂 𝜙 𝑗 〉

𝑂 ≃
𝑗=0

𝑂[𝜙]に𝜙(𝑛)の展開式を代入してgでテイラー展開。 𝑂 𝜙 (𝑗) は畳み込み計算となる。
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

6

１．格子場の理論での計算精度
一．多重積分の自由度（次元）がとても大きいことに起因する精度上の困難• 手前味噌な話

𝜒 𝑛 (𝑗) = 𝜙 𝑛 ⋆ 𝜓 𝑛(𝑗)𝑗𝜙 𝑛 (𝑘) 𝜒 𝑛 (𝑗−𝑘)=𝑘=0• 多項式積の係数値の大きさが組み合わせ爆発で大きくなる例）𝜙 𝑛 (𝑗) ≃ 𝑂(1)を模して

𝜙 = ∑𝑔 𝑗 𝜙 𝑗 ≃ 1 + 𝑔 + 𝑔2 + ⋯ + 𝑔63

𝜙11 ≃ 1 + 𝑔 + 𝑔2 + ⋯ + 𝑔63 11
1 + 𝑔 + 𝑔2 + ⋯ +𝑔63 11 = 1 + 11𝑔 + 66𝑔2 + …
… + 536211932256𝑔62 + 621324937376𝑔63• FFTを用いた畳み込みの繰り返しで計算すると最後の方で桁落ちが生じて精度が出ない。

• 𝑔 のスケーリングで指数的爆発は吸収できるが、組み合わせ爆発は追いつかない。高次摂動計算で複雑な物理量（ベキが高くなるような物理量）は４倍精度FFTが必要。
• 多項式行列（複素数行列が多項式係数）の場合∞∞

𝑔 𝑗 𝐴(𝑗) , 𝐵 =𝐴=
𝑗=0

𝑔 𝑗 𝐵 (𝑗) , 𝐶 =𝑗=0

𝐶 (𝑗) = 𝐴 ⋆ 𝐵 (𝑗) =𝑔 𝑗 𝐶 (𝑗)𝑗=0(𝑤)

𝐴

2019/3/18𝑗∞= 𝐹𝐹𝑇 𝐴(𝑗),𝑤 ,𝐶

第3回HPC-Phys 勉強会(2019/3/18)@広島大学(𝑗)= 𝐼𝐹𝐹𝑇 𝐶𝑤

𝐴(𝑘) 𝐵(𝑗−𝑘)𝑘=0, 𝑗 , 𝐶 (𝑤) = 𝐴(𝑤) 𝐵(𝑤) ,7

１．格子場の理論での計算精度
一．運動方程式を解く際に生じる困難。特に差分方程式＝大規模連立方程式の係数行列の条件数がとて𝐴𝑥𝑏も大きいときの困難𝐴𝑥 = 𝑏

行列𝐴の大きさは格子点数×格子点数に比例。
𝛾𝜇 𝜕𝜇 + 𝑖𝐴𝜇 𝑥+ 𝑚 𝑆(𝑥, 𝑦) = 𝛿(𝑥 − 𝑦)𝐷 𝑛, ℓ 𝑆(ℓ, 𝑚) = 𝛿(𝑛, 𝑚)ℓ𝑆(𝑥, 𝑦) : ディラック方程式のグリーン関数（ディラック方程式：電子とかクォークとかスピン1/2を持つ相対論的粒子の運動方程式）

𝑆 (𝑛, 𝑚) : 格子上のディラック方程式（＝ウィルソン・ディラック方程式）のグリーン関数

41 − 𝛾𝜇 𝑈𝜇 𝑛 𝛿(𝑛 + 𝜇, 𝑚) + 1 + 𝛾𝜇 𝑈𝜇† (𝑚)𝛿(𝑛 − 𝜇, 𝑚)

𝐷 𝑛, 𝑚 = 𝛿 𝑛, 𝑚 − 𝜅
𝜇=1
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

8

１．格子場の理論での計算精度

一．運動方程式を解く際に生じる困難。特に差分方
程式＝大規模連立方程式の係数行列の条件数がとて𝐴𝑥𝑏も大きいときの困難
𝐴𝑥 = 𝑏

行列𝐴の大きさは格子点数×格子点数に比例。
2564 だと4294967296ｘ4294967296• 一般に反復法で解く。
• 共役勾配法（CG法）、安定化双共役勾配法（BiCGStab法）など
条件数

𝐾 𝐴 ≡ | 𝐴 |/| 𝐴−1 |• 条件数（最大固有値と最小固有値の比、ダイナミックレンジ）が大きくなると、反復法が収束しなくなる• 自由度（次元）が大きいと反復法内で計算しているベクトルの内積（示量性をもつ）の桁あふれとか桁落ちも気になる。
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

9１．格子場の理論での計算精度

一．運動方程式を解く際に生じる困難。特に差分方
程式＝大規模連立方程式の係数行列の条件数がとても大きいときの困難𝐴𝑥𝑏

𝐴𝑥 = 𝑏• 格子量子色力学の場合は？

• Dの最小固有値の大きさは格子間隔aを単位としたクォーク質量 am• Dの最大固有値は無次元量の 1 程度
𝐾 𝐷 ≡𝐷
𝐷−1=1𝑎𝑚• 世の中の最も軽いクォークの質量は大体 2MeV ぐらい• 格子間隔 a は

256 a = 10 fm とすると、 a = 1/(5 GeV), 1/a = 5 000MeV

𝐾 𝐷 =15000== 2500𝑎𝑚2• まだ大丈夫？
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学10

2. 格子場の理論での乱数• 経路積分＝多重積分の計算にモンテカルロ法を用います。
• モンテカルロ法：何らかの形で乱数を使ってなにか計算する方法を一般的にモンテカルロ法と言うそうです。

• 格子場の理論では多重積分の計算にマルコフチェーンモンテカルロ法を用います。
1
𝑂𝐿 =𝑍𝑛重み exp[−𝑆𝐿 [𝑔, 𝜙]] 付き積分
重みが非負実数ならモンテカルロ法が使える

𝑑𝜙 𝑛 𝑂𝐿 [𝜙]exp[−𝑆𝐿 [𝜙]]• 場の変数がある値𝜙にあるとき、次の場の値𝜙′を何らかの確率過程で生成します。次の場の値は一個前の値の場の値と遷移確率行列𝑃(𝜙′, 𝜙)で結びつきます。
{… , 𝜙 −1 , 𝜙 0 , 𝜙 1 , 𝜙 2 , 𝜙 3 , … , 𝜙 𝑁𝑠𝑎𝑚𝑝 −1 , 𝜙 𝑁𝑠𝑎𝑚𝑝 }• この遷移確率行列がある性質を満たすと、定常分布が多重積分に現れている（測度付）重みexp −𝑆𝐿 𝑔, 𝜙 𝐷𝜙になり、統計平均で多重積分が評価できま𝑁𝑠𝑎𝑚𝑝
す。
1𝑂𝐿 =

2019/3/18

lim

𝑁𝑠𝑎𝑚𝑝 →∞ 𝑁𝑠𝑎𝑚𝑝

𝑂𝐿 [𝜙(𝑖) ]𝑘=1

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

11

2. 格子場の理論での乱数• 経路積分＝多重積分の計算にモンテカルロ法を用います。
• 場の変数がある値𝜙にあるとき、次の場の値𝜙′を何らかの確率過程で生成します。次の場の値は一個前の値の場の値と遷移確率行列𝑃(𝜙′, 𝜙)で結びつきます。
{… , 𝜙 −1 , 𝜙 0 , 𝜙 1 , 𝜙 2 , 𝜙 3 , … , 𝜙 𝑁𝑠𝑎𝑚𝑝 −1 , 𝜙 𝑁𝑠𝑎𝑚𝑝 }• 格子上に乗っている場の変数をMCMCで変化させて状態遷移させていくのですが、ここに乱数（擬似乱数）を用います。
• 基本的に各格子点上で乱数を振ります。
• 格子点ごとに独立な乱数が望ましいです。

• 格子点以外に自由度を保つ場合もありますが、並列化を考えると格子の少領域ごと（Open MP のスレッドレベル）に乱数生成器があると良いです。
• ガウス分布、指数分布、変な分布を実数一様乱数から作って場にセットしてターゲットの場を変化させます。
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

12

2. 格子場の理論での乱数• 経路積分＝多重積分の計算にモンテカルロ法を用います。
• 場の変数がある値𝜙にあるとき、次の場の値𝜙′を何らかの確率過程で生成します。次の場の値は一個前の値の場の値と遷移確率行列𝑃(𝜙′, 𝜙)で結びつきます。
{… , 𝜙 −1 , 𝜙 0 , 𝜙 1 , 𝜙 2 , 𝜙 3 , … , 𝜙 𝑁𝑠𝑎𝑚𝑝 −1 , 𝜙 𝑁𝑠𝑎𝑚𝑝 }• ガウス分布、指数分布、変な分布を実数一様乱数から作って場にセットしてターゲットの場を変化させます。
• 例）確率過程量子化、分子動力学法、ハイブリドモンテカルロ法、確率微分方程式⇒ランジュバン方程式𝑑𝜙(𝑡)𝜕𝑆𝑛 =−𝜙 𝑡 + 𝜂 𝑡 (𝑛)𝑑𝑡𝜕𝜙 𝑛
オイラー法によるランジュバン時間積分：
各時刻ステップでガウス分布する数値を𝑊(𝑛)にセット

𝜂𝑡 𝑛

𝜂 𝑡 𝑛 𝜂 𝑡′ 𝑚

= 2𝛿 𝑡 − 𝑡 ′ 𝛿𝑛,𝑚

𝜂 𝑡 𝑛 ：外場、ノイズ場、ランジュバン時刻と格子点座標に関して独立𝜕𝑆
𝜙(𝑗+1) 𝑛 = 𝜙(𝑗) 𝑛 − Δ𝑡
𝜙 𝑗 + 𝑊 𝑗 (𝑛)𝜕𝜙 𝑛
2019/3/18=0

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

𝑊𝑗 𝑛=0𝑊 𝑗 𝑛 𝑊 𝑘′ 𝑚

= 2Δ𝑡𝛿𝑗,𝑘 𝛿𝑛,𝑚13

2. 格子場の理論での乱数• 例）確率過程量子化、分子動力学法、ハイブリドモンテカルロ法、確率微分方程式⇒ランジュバン方程式
オイラー法によるランジュバン時間積分：
各時刻ステップでガウス分布する数値を𝑊(𝑛)にセット

𝜕𝑆
𝜙(𝑗+1) 𝑛 = 𝜙(𝑗) 𝑛 − Δ𝑡
𝜙 𝑗 + 𝑊 𝑗 (𝑛)𝜕𝜙 𝑛𝑊𝑗 𝑛=0

𝑊 𝑗 𝑛 𝑊 𝑘′ 𝑚

= 2Δ𝑡𝛿𝑗,𝑘 𝛿𝑛,𝑚

• 計算機でこれを計算する。

• 格子サイズ（自由度）が大きくなると並列化が必要となる。
領域毎に独立な乱数系
列乱数生成器

乱数生成器

乱数生成器
2019/3/18

並列度（MPI+OpenMP）は近年では1000～10000 にもなる乱数生成器

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

14

2. 格子場の理論での乱数• 例）確率過程量子化、分子動力学法、ハイブリドモンテカルロ法、確率微分方程式⇒ランジュバン方程式
オイラー法によるランジュバン時間積分：
各時刻ステップでガウス分布する数値を𝑊(𝑛)にセット

𝜕𝑆
𝜙(𝑗+1) 𝑛 = 𝜙(𝑗) 𝑛 − Δ𝑡
𝜙 𝑗 + 𝑊 𝑗 (𝑛)𝜕𝜙 𝑛

𝑊𝑗 𝑛=0𝑊 𝑗 𝑛 𝑊 𝑘′ 𝑚= 2Δ𝑡𝛿𝑗,𝑘 𝛿𝑛,𝑚• 計算機でこれを計算する。

• 格子サイズ（自由度）が大きくなると並列化が必要となる。
• 独立性（相関がない）ように• 乱数の種を変えるのはあまり望ましくない。重なりがあるかも。

• 乱数の話 ：松本さん斎藤さん（メルセンヌツイスタ）• 1系列を

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html• 飛ばし飛ばしに使う• うんと飛ばしたものを分割して使う( jump ahead Mersenne Twister)• 複数ストリーム： Intel MKL のライブラリにはあった。CUDA(CURAND)?
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

15

2. 格子場の理論での乱数• 乱数の話 ：松本さん斎藤さん（メルセンヌツイスタ）• 1系列を

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html• 飛ばし飛ばしに使う• うんと飛ばしたものを分割して使う(Jump ahead Mersenne Twister[H. Haramoto, M. Matsumoto, T. Nishimura, F. Panneton, and P. L‘Ecuyer, ``Efficient Jump Ahead for
F_2-Linear Random Number Generators’‘ (2008)]• 石川：

Fortranner用：Multiple stream Mersenne Twister PRNG

http://theo.phys.sci.hiroshima-u.ac.jp/~ishikawa/PRNG/mt_stream.html• 格子場の理論では乱数生成の計算コストの重みはまだ低い（ランジュバン法、クラマース法はやや重め、分子動力学法、HMC法では低い）• 陥りやすいミス：• 初心者：初めて書いた並列プログラムの計算で各プロセスが同じ乱数を使っていて失敗。変な相関がある…
• われわれ：パラメータスキャンで複数の異なるジョブを管理するときにうっかり乱数の種の同じジョブを動かしてしまった。変な相関がある…
本当は乱数の重なりがないことも保証しないといけないけど…
2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

16

3．まとめ
格子場の理論での計算精度と乱数の利用について紹介しました。
自由度が大きくなった時や、大規模並列計算で陥りそうな問題を紹介しました。
• 以上です。ありがとうございました。

2019/3/18

第3回HPC-Phys 勉強会(2019/3/18)@広島大学

17