Slide 1

Slide 1 text

1 79 モンテカルロ法(3) 発展的なアルゴリズム シミュレーション工学 慶應義塾大学大学院理工学研究科基礎理工学専攻物理情報専修 渡辺宙志

Slide 2

Slide 2 text

2 79 • モンテカルロサンプリングと分配関数の関係 を知る • アルゴリズムを工夫することで、効率的なサ ンプリングができることを知る 本講義の目的 モンテカルロ法とは • 数値計算では、何かの和や積分の推定値を 計算することが多い • 和の形を変形することで、異なるアルゴリ ズムが生まれる

Slide 3

Slide 3 text

3 79 • 格子の各点にスピン(小さな磁石)がある • スピンは「上」と「下」の状態がある • 隣り合うスピンをつなぐ線をボンドと呼ぶ ボンドの両側の スピンの向き 同じ 逆 エネルギー −𝐽 𝐽 𝐽 > 0ならスピンは揃いたがる(強磁性的) 𝐽 < 0ならスピンは逆向きを好む(反強磁性)

Slide 4

Slide 4 text

4 79 格子の上の全てのボンドについてエネルギー の和を取り、この系の全エネルギーとする このような模型をイジング模型(Ising Model)と呼び、 磁性体の簡単なモデルになっている 以下、強磁性 (𝐽 > 0)の場合を考える 𝐸 = −2𝐽 + 2𝐽 = 0

Slide 5

Slide 5 text

5 79 𝑖番目のスピンの状態を𝜎𝑖 とする 𝜎𝑖 の値は1(上向き)か-1(下向き)のいずれか −𝐽 𝐽 𝜎𝑖 = 1 𝜎𝑗 = 1 𝜎𝑖 = 1 𝜎𝑗 = −1 両方まとめて−𝐽𝜎𝑖 𝜎𝑗 と書ける

Slide 6

Slide 6 text

6 79 全系のエネルギーは以下のように書ける 𝐻 = −𝐽 ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 系の全てのボンドについて和をとるという意味 全系のエネルギーを与える量をハミルトニアンとよぶ

Slide 7

Slide 7 text

7 79 系の状態に通し番号をつけ、𝑖番目の状態のエネルギーを𝐸𝑖 とする 状態𝑖 エネルギー 𝐸𝑖 = 4𝐽 exp(−𝛽𝐸𝑖 ) ボルツマン定数 𝛽 = 1/𝑘𝐵 𝑇 逆温度 𝑘𝐵 この状態の出現確率がボルツマン重みに比例する この時、エネルギーの温度依存性を知りたい

Slide 8

Slide 8 text

8 79 低温 高温 スピンがそろった状態 出現確率は高いが、総数が少ない →エネルギー重視 →低温で支配的 スピンがバラバラの状態 出現確率は低いが、総数が多い →エントロピー重視 →高温で支配的

Slide 9

Slide 9 text

9 79 エネルギーの期待値の温度依存性 𝑈 𝛽 = 𝑍−1 ෍ 𝑖 𝐸𝑖 exp −𝛽𝐸𝑖 𝑍 = ෍ 𝑖 exp −𝛽𝐸𝑖 エネルギー𝐸をとる状態の数を𝑔(𝐸)とすると 𝑈 𝛽 = 𝑍−1 න 𝐸𝑔(𝐸) exp −𝛽𝐸 𝑑𝐸 𝑍 = න 𝑔(𝐸) exp −𝛽𝐸 𝑑𝐸 状態からエネルギー の計算は簡単 𝐸 = 0 あるエネルギーをとる状態の 数の計算は大変 𝐸 = 0

Slide 10

Slide 10 text

10 79 𝐸 = 0 𝐸 = 4J 𝐸 = −4J 𝑔 −4𝐽 = 2 𝑔 0 = 12 𝑔 4𝐽 = 2 これを一般のNで求めるのは極めて大変 モンテカルロサンプリング

Slide 11

Slide 11 text

11 79 1. 全ての可能な状態から無作為に一つ選ぶ 2. その状態のエネルギー𝐸𝑖 と重み𝑤𝑖 を計算する 3. 1-2を繰り返し、重み付きで平均をとる スピンが𝑁個なら、状態は2𝑁個 ほとんどの状態は重みが小さいので サンプリング効率が悪い マルコフ連鎖モンテカルロ法 この中から無作為に一つ選ぶ

Slide 12

Slide 12 text

12 79 現在の状態から遷移可能な状態を限定し、その中 から一つ無作為に選ぶ サンプリング候補の決め方 イジング模型への適用 1. 現在の状態を𝐴とし、スピンを一つ無作為に選ぶ 2. 選んだスピンを反転させた状態を提案状態𝐵とする 3. それぞれのエネルギーから遷移確率𝑃(𝐴 → 𝐵)を計算 し、遷移させるか決める 4. 遷移しなかった場合は状態はそのまま、遷移した場 合は提案状態を現在の状態にして1.へ

Slide 13

Slide 13 text

13 79 現在の状態 遷移可能な状態 スピンが𝑁個なら総状態数は2𝑁、遷移可能な状態は𝑁 この中から無作為に一つ選ぶ 遷移可能な状態が限定された

Slide 14

Slide 14 text

14 79 現在の状態𝐴 現在の状態𝐵 𝑃(𝐴 → 𝐵) 𝑃(𝐵 → 𝐴) 𝑃(𝐵 → 𝐵) 𝑃(𝐴 → 𝐴) 𝐸𝐴 = −4𝐽 𝐸𝐵 = 0 𝑤𝐴 = exp −𝛽𝐸𝐴 𝐸𝐴 < 𝐸𝐵 なので𝑤𝐴 > 𝑤𝐵 𝑤𝐵 = exp −𝛽𝐸𝐵

Slide 15

Slide 15 text

15 79 𝑃(𝐴 → 𝐵) 𝑃(𝐵 → 𝐴) 𝑃(𝐵 → 𝐵) 𝑃(𝐴 → 𝐴) 確率の保存条件 𝑃 𝐴 → 𝐴 + 𝑃 𝐴 → 𝐵 = 1 𝑃 𝐵 → 𝐵 + 𝑃 𝐵 → 𝐴 = 1 詳細釣り合い条件 𝑃 𝐴 → 𝐵 𝑃 𝐵 → 𝐴 = 𝑤𝐵 𝑤𝐴

Slide 16

Slide 16 text

16 79 メトロポリス法(Metropolis method) 重みが大きくなる場合は必ず遷移、そうでない場合は確率的に遷移させる 𝑃 𝐵 → 𝐴 = 1 𝑃 𝐴 → 𝐵 = 𝑤𝐵 𝑤𝐴 = exp −𝛽𝐸𝐵 exp(−𝛽𝐸𝐴 ) = exp(−𝛽Δ𝐸) Δ𝐸 = 𝐸𝐵 − 𝐸𝐴 > 0なので、 exp −𝛽Δ𝐸 < 1 提案状態のエネルギーが、現在の状態より低ければ必ず遷移 高ければ確率𝑝 = exp −𝛽Δ𝐸 < 1で遷移

Slide 17

Slide 17 text

17 79 イジング模型にマルコフ連鎖モンテカルロ法を適用した場 合のアルゴリズム(メトロポリス法を採用した場合) 1. スピンを一つ無作為に選ぶ 2. 選んだスピンを反転させた状態を提案状態とする 3. 現在の状態とのエネルギー差Δ𝐸を計算する 4. エネルギーが下がる場合(Δ𝐸 < 0)なら必ず遷移。そ うでなければ確率𝑝 = exp(−𝛽Δ𝐸)で遷移 5. 1-4を繰り返す 一度に一つのスピンだけひっくり返すので Single-Spin-Flip Algorithmと呼ぶ

Slide 18

Slide 18 text

18 79 状態に重みがある場合の期待値を厳密に計算することは難しい その理由は同じ重みを持つ状態の数がわからないから →サンプリングによる期待値の推定 単純サンプリング 全ての状態から無作為に状態をサンプリングする 重みの小さい状態ばかりサンプリングするため非効率 マルコフ連鎖モンテカルロ法 現在の状態をもとに、サンプリングする状態を限定する 現在の状態と提案状態の重みを比較し、適切な遷移確率を 作ることで重みに比例したサンプリングができる (重み付きサンプリング)→効率的

Slide 19

Slide 19 text

19 79 一度に一つしかスピンをひっくり返せないのは非効率的? 複数のスピンを同時にひっくり返せたら計算が加速 するのでは?

Slide 20

Slide 20 text

20 79 同じ向きのスピンの数 𝑛+ = 3 逆向きのスピンの数 𝑛− = 1 同じ向きのスピンの数 𝑛+ = 1 逆向きのスピンの数 𝑛− = 3 ひっくり返る 𝐸𝐴 = 𝑛− − 𝑛+ 𝐽 𝐸𝐵 = 𝑛+ − 𝑛− 𝐽 現在の状態𝐴 提案状態𝐵 Δ𝐸 = 𝐸𝐵 − 𝐸𝐴 = 2 𝑛+ − 𝑛− 𝐽 エネルギー差が同じ向きのスピンの数を数えるだけで計算できる

Slide 21

Slide 21 text

21 79 現在の状態𝐴 提案状態𝐵 この二つを同時に ひっくり返そうとする あるスピンの状態が確定していないと、隣のスピンを ひっくり返す試行ができない ? このスピンの状態が 確定しないと このスピンをひっくり返す ための重みが計算できない

Slide 22

Slide 22 text

22 79 ボンドで接続されていないスピンは影響を受けない あるスピンがひっくり返ると、ボンドで繋がったスピン の重みが変わる このスピンがひっくり返ると このスピンの重みに影響がでる このスピンがひっくり返っても このスピンの重みには影響しない

Slide 23

Slide 23 text

23 79 これらのスピンを同時にひっくり返そうとする ? このスピンをひっくり返す 試行の結果はまだ知らない わからないまま確率を計算して ひっくり返してしまう ?

Slide 24

Slide 24 text

24 79 もしスピンがひっくり返っていなければ? 一つスピンをひっくり 返しただけ もしスピンがひっくり返っていたら? 二つのステップを 逐次的に行ったのと等価 これらの状態間で詳細釣り合い条件が成り立っているので 正しいサンプリングができている

Slide 25

Slide 25 text

25 79 ボンドで繋がっていない二つのグループに分け、 Aサイト、Bサイトと名前をつける A A B B Bサイトのスピンを固定し、 Aサイトのスピンを全て独立 にひっくり返す試行を行う Aサイトのスピンを固定し、 Bサイトのスピンを全て独立 にひっくり返す試行を行う

Slide 26

Slide 26 text

26 79 Q. 何がうれしいの? → A. 並列化できる def flip(i, spin, beta): dE = calc_energy_diff(i) if check_flip(dE): spin[i] = -spin[i] flip(i) スレッド1 flip(j) スレッド2 他のスレッドの結果を待たずに自分の担当のスピンをひっくり返してよい

Slide 27

Slide 27 text

27 79 16個の状態を 4x4個の状態に分離した x =

Slide 28

Slide 28 text

28 79 𝑍 = ෍ {𝜎} exp(−𝛽𝐻) = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp(−𝛽𝐻) 全てのスピンの状態についての和 2𝑁通り A、Bサイトそれぞれのスピンの状態に ついての和 2𝑁/2通り 𝐻 = −𝐽 ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 ハミルトニアンを代入する 𝑍 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗

Slide 29

Slide 29 text

29 79 𝑍 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 外側を固定して 内側の和についてのみサンプリングする 𝜎𝑖 𝜎𝑗1 + 𝜎𝑖 𝜎𝑗2 + 𝜎𝑖 𝜎𝑗3 + 𝜎𝑖 𝜎𝑗4 = ℎ𝑖 𝜎𝑖 𝑖 𝑗3 𝑗2 𝑗4 𝑗1 外側のスピンが固定されているので、 それをまとめてしまう ℎ𝑖 ≡ 𝜎𝑗1 + 𝜎𝑗2 + 𝜎𝑗3 + 𝜎𝑗4

Slide 30

Slide 30 text

30 79 𝑍 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖 ℎ𝑖 𝜎𝑖 𝑍 = ෍ {𝜎} exp(−𝛽𝐻) = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp(−𝛽𝐻) だったから 𝐻 = −𝐽 ෍ 𝑖 ℎ𝑖 𝜎𝑖 というハミルトニアンを計算している気持ちになった ローカルな磁場下にある独立なスピンの集まり → 相互作用が消えてしまった → 相互作用がないので好き勝手できる(例えば並列化)

Slide 31

Slide 31 text

31 79 𝑍 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 = ෍ {𝜎𝐵} ෍ {𝜎𝐴} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 和の順序を入れ替えれば、Aサイトについて同じことができる アルゴリズム 1. Aサイトのスピンを固定し、Bサイトのスピンだけ更新する 2. Bサイトのスピンを固定し、Aサイトのスピンだけ更新する 3. 以下繰り返し 計算の手間は変わらないが、更新しようとするスピンが全て独立で あるために並列化が可能になる

Slide 32

Slide 32 text

32 79 • 分配関数の和の取り方を工夫すると、もとの系とは異な る形のハミルトニアンが現れる→異なる物理系をシミュ レーションしているかのように見える • 多重和の外側を固定し、内側だけサンプリングすると、 スピンが独立になるために計算が楽になることがある もっともうまく分配関数を変形すれば、より効率 的なサンプリングが可能になるのでは? クラスターアルゴリズム

Slide 33

Slide 33 text

33 79 𝑍 = ෍ {𝜎} exp(−𝛽𝐻) 𝐻 = −𝐽 ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 イジング模型のハミルトニアン その分配関数 𝑍 = ෍ {𝜎} exp ෍ 𝑖,𝑗 𝐾𝜎𝑖 𝜎𝑗 = ෍ {𝜎} ෑ 𝑖,𝑗 exp 𝐾𝜎𝑖 𝜎𝑗 ハミルトニアンを代入した 𝐾 ≡ 𝛽𝐽 指数関数の肩の和を外に出して 積に書き換えた これを書き直す

Slide 34

Slide 34 text

34 79 exp 𝐾𝜎𝑖 𝜎𝑗 = ൝ 𝑒𝐾 𝜎𝑖 = 𝜎𝑗 𝑒−𝐾 𝜎𝑖 ≠ 𝜎𝑗 であるから exp 𝐾𝜎𝑖 𝜎𝑗 = 𝑒−𝐾 + 𝛿𝜎𝑖,𝜎𝑗 𝑒𝐾 − 𝑒−𝐾 = 𝑒𝐾 𝑒−2𝐾 + 𝛿𝜎𝑖,𝜎𝑗 1 − 𝑒−2𝐾 = 𝑒𝐾 (1 − 𝑝) + 𝑝𝛿𝜎𝑖,𝜎𝑗 ≡ 𝑝 ここまではただの等式変形

Slide 35

Slide 35 text

35 79 以下のデルタ関数を使った恒等式を使って書き直す 𝑥 + 𝑦 = ෍ 𝑏=0 1 𝑥𝛿𝑏,0 + 𝑦𝛿𝑏,1 𝑒𝐾 (1 − 𝑝) + 𝑝𝛿𝜎𝑖,𝜎𝑗 = ෍ 𝑏𝑖𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 新たに導入された𝑏𝑖𝑗 はスピン𝑖とスピン𝑗を結ぶボンド上で 定義されており、ボンド自由度と呼ぶ

Slide 36

Slide 36 text

36 79 𝑏𝑖𝑗 はスピンiとjをつなぐボンド上で定義される 値として1か0をとる 1の時にactive、0の時にinactiveであると定義する ボンド変数 スピン変数 𝜎𝑖 はサイトi上で定義される 値として1か-1をとる 1の時に上向き、-1の時に下向きと定義する 𝜎𝑖 𝜎𝑗 𝑏𝑖𝑗

Slide 37

Slide 37 text

37 79 exp(𝐾𝜎𝑖 𝜎𝑗 ) = ෍ 𝑏𝑖𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 = + 新たにボンド自由度を導入し、スピンの相互作用の 重みをactiveなボンドとinactiveなボンドの和に分けた inactive active

Slide 38

Slide 38 text

38 79 以上の式変形をまとめると 𝑍 = ෍ {𝜎} exp(−𝛽𝐻) = ෍ {𝜎𝑖} ෍ {𝑏𝑖𝑗} ෑ 𝑖,𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 隣り合うスピンに関する積 ボンドに関する和 スピンに関する和 等価変形なのにボンドに関する自由度が増えた

Slide 39

Slide 39 text

39 79 𝑍 = ෍ {𝜎𝑖} ෍ {𝑏𝑖𝑗} ෑ 𝑖,𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 ≡ 𝑆𝑖𝑗 • 新たにボンド自由度がうまれ𝑏𝑖𝑗 = 0か𝑏𝑖𝑗 = 1の値をとる • 分配関数が多重和で書かれているので、スピンに関するサ ンプリングとボンドに関するサンプリングを交互にとれる • ボンド𝑏𝑖𝑗 がactive(𝑏𝑖𝑗 = 1)であるか、inactiveであるかの確 率は、スピン状態に依存する重みで決まる スピン状態を固定した状態で ボンド状態のサンプリングを考える

Slide 40

Slide 40 text

40 79 両側のスピンが逆向き(𝜎𝑖 ≠ 𝜎𝑗)の時 𝑆𝑖𝑗 = 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 𝑆𝑖𝑗 = 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 𝑏𝑖𝑗 = 0の場合しか重みを持たない → 必ずinactive 両側のスピンが同じ向き(𝜎𝑖 ≠ 𝜎𝑗)の時 𝑆𝑖𝑗 = 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝑏𝑖𝑗,1 確率1 − 𝑝でinactive、確率𝑝でactive スピン変数を固定

Slide 41

Slide 41 text

41 79 現在の状態 ボンド状態を忘れる 新しいボンド状態を与える 隣り合うスピンが同じ向きの場合 • 確率𝑝でactive • 確率1 − 𝑝でinactive 隣り合うスピンが逆向きの場合 • 必ずinactive 𝑝 ≡ 1 − 𝑒−2𝐾

Slide 42

Slide 42 text

42 79 𝑆𝑖𝑗 = 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 ボンド変数を固定 ボンドがactive(𝑏𝑖𝑗 = 1)の時 𝑆𝑖𝑗 = 𝑒𝐾𝑝𝛿𝜎𝑖,𝜎𝑗 ボンドがinactive(𝑏𝑖𝑗 = 0)の時 𝑆𝑖𝑗 = 𝑒𝐾(1 − 𝑝) スピンの相互作用が含まれない → スピンの相互作用が切れている

Slide 43

Slide 43 text

43 79 activeなボンドでつながっているスピンをクラスターと呼ぶ

Slide 44

Slide 44 text

44 79 inactiveなボンドでつながるスピンは相互作用を持たない → activeなボンドでつながるスピンを一斉にひっくり返しても 状態の重みが変わらない 状態𝐴 状態𝐵 𝑤𝐴 = 𝑤𝐵

Slide 45

Slide 45 text

45 79 𝑃(𝐴 → 𝐵) 𝑃(𝐵 → 𝐴) 熱浴法を採用すると 𝑃 𝐴 → 𝐵 = 𝑤𝐵 𝑤𝐴 + 𝑤𝐵 = 1 2 𝑃 𝐵 → 𝐴 = 𝑤𝐴 𝑤𝐴 + 𝑤𝐵 = 1 2 状態𝐴 状態𝐵

Slide 46

Slide 46 text

46 79 𝑍 = ෍ {𝑏𝑖𝑗} ෍ {𝜎𝑖} ෑ 𝑖,𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 ボンド状態を固定した状態で スピン状態のサンプリングを考える activeなボンドでつながるスピンを一斉にひっくり返しても重みが等しい →同じクラスターに所属するスピンを全てランダムにひっくり返してよい クラスターごとに、元のスピン配置を忘れて ランダムにスピンを割り当ててよい ※同じクラスターに属すスピンは必ず同じ向き

Slide 47

Slide 47 text

47 79 1. スピンを固定してボンド状態を更新する • 隣り合うスピンが同じ向きなら確率pでactive • 逆向きならinactive 2. ボンド状態を固定してスピン状態を更新する 1. activeなボンドでつながったスピンをクラスターに 分ける 2. クラスターごとに新しいスピン状態を与える 3. 1-2のステップを繰り返す 以上の方法をSwendsen-Wangのアルゴリズムと呼ぶ R. H. Swendsen, J.-S. Wang, Phys. Rev. Lett. vol. 58, pp. 86–88 (1987)

Slide 48

Slide 48 text

48 79 R. H. Swendsen N. Ito H. Watanabe

Slide 49

Slide 49 text

49 79 𝑍 = ෍ {𝑏𝑖𝑗} ෍ {𝜎𝑖} ෑ 𝑖,𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 外側にボンド状態に関する和がある ボンド状態は、たとえば左図のよう なグラフを表している →あらゆるグラフに関する和 • このグラフ表現をFortuin-Kasteleyn表現と呼ぶ • グラフ表現を利用してスピンを更新するアルゴリズムを クラスターアルゴリズムと呼ぶ • Swendsen-Wangはクラスターアルゴリズムの一種 C.M.Fortuin, P.W.Kasteleyn, Physica, vol. 57, pp 536 (1972).

Slide 50

Slide 50 text

50 79 クラスターアルゴリズムは非常に早い しかし、早い理由はスピンをまとめて更新するからではない 低温である時、以下のようなドメインウォールができたとする 低温だから全てのスピンがそろった状態になりたい

Slide 51

Slide 51 text

51 79 Single-Spin-Flipでドメインウォール解消を試みる 途中で必ずエネルギーが高い状態を経由する必要がある 低温であるほど、高エネルギー状態は出現しづらい →極めて緩和が遅くなる エネルギーの損 4J エネルギーの損 6J エネルギーの損 0

Slide 52

Slide 52 text

52 79 Swendsen-Wangの場合 現在の状態 ボンド状態を更新して スピンを忘れる 確率1/2で全ての スピンがそろう クラスターアルゴリズムでは、ドメインウォールを非常 に早く解消できる→緩和が本質的に早い

Slide 53

Slide 53 text

53 79 • 分配関数を式変形し、新たにボンドの自由度を追加すると グラフの和の形に書ける(Fortuin-Kasteleyn表現) • グラフ表現を利用したサンプリングアルゴリズムをクラス ターアルゴリズムと呼ぶ • Inactiveなボンドでつながったスピンは重みがスピン状態に 依存しない→実質的に相互作用が切れている • スピンをお互いに相互作用しないクラスターに分けること ができる • クラスターはそれぞれ独立にスピン状態を与えられる • ボンド変数とスピンを交互に更新することで、効率的なサ ンプリングができる • Swendseng-Wangはクラスターアルゴリズムの一種で、他 にもWolffのアルゴリズムなどがある

Slide 54

Slide 54 text

54 79 • 状態𝑖の出現確率が𝑤𝑖 に比例する時、𝑤𝑖 を状態𝑖の重みと呼ぶ • 重みの総和𝑍 = σ𝑖 𝑤𝑖 を分配関数と呼ぶ • マルコフ連鎖モンテカルロ法とは分配関数をサンプリングに より推定する手法 • 分配関数を変形すると、異なる物理系のように見える • 同じ分配関数であるにも関わらず、変形をすることで並列化 できたり緩和を加速できたりする Fortuin-Kasteleyn表現は1972年には知られていた。 しかしSwendsen-Wangアルゴリズムの提案は1987年 まだ見ぬ新たなブレイクスルーがあるかも・・・?

Slide 55

Slide 55 text

55 79 確率変数 ෠ 𝑋の母集団のパラメータを知りたい 𝐸 ෠ 𝑋 = ෍ 𝑘 𝑌𝑘 𝑝𝑘 ෠ 𝑋の取りうる値𝑌𝑘 とその実現確率𝑝𝑘 が既知なら 母集団の期待値の厳密な表記は 一般にはこの和を厳密に計算できない →サンプリングで推定

Slide 56

Slide 56 text

56 79 確率変数 ෠ 𝑋の母集団のパラメータを知りたい 確率変数 ෠ 𝑋を𝑁回観測して標本 {𝑋𝑖 }を得た ここから ෠ 𝑋の母集団のパラメータを推定したい 𝐸 ෠ 𝑋 ∼ ത 𝑋 = 1 𝑁 ෍ 𝑖 𝑋𝑖 母集団の期待値 期待値のestimator 知りたい量を標本で表す関数をestimatorと呼ぶ

Slide 57

Slide 57 text

57 79 𝑁が大きい時、中心極限定理により、 ത 𝑋の分布はガウス分布に近づく ത 𝑋の分散のestimator 𝑉 ത 𝑋 = 1 𝑁 𝑉 ෠ 𝑋 = 1 𝑁(𝑁 − 1) ෍ 𝑖 𝑋𝑖 − ത 𝑋 2 ത 𝑋の分散 ෠ 𝑋の分散(母分散) これが統計誤差だから、なるべく小さくしたい

Slide 58

Slide 58 text

58 79 𝐸 ෠ 𝑋 = ෍ 𝑘 𝑌𝑘 𝑝𝑘 全体では厳密な和を取れないが、状態がいくつか のグループに分かれ、その中の部分和は厳密に取 れるとする ここは計算できないが ここは厳密に計算できる = ෍ 𝑔 ෍ 𝑘∈𝑔 𝑌𝑘 𝑝𝑘 = ෍ 𝑔 ത 𝑌 𝑔 𝑝𝑔 すると、対応するestimatorの分散が小さくなる これをimproved estimatorと呼ぶ

Slide 59

Slide 59 text

59 79 𝐸 ෠ 𝑋 = ෍ 𝑘 𝑌𝑘 𝑝𝑘 = ෍ 𝑔 ෍ 𝑘∈𝑔 𝑌𝑘 𝑝𝑘 = ෍ 𝑔 ത 𝑌 𝑔 𝑝𝑔 確率𝑝𝑔 で状態𝑔となる確率変数 ෠ 𝑋I を考えると もとの確率変数 ෠ 𝑋よりも分散が小さくなる 𝑉 ෠ 𝑋 > 𝑉 ෠ 𝑋I 対応するestimatorの分散も小さくなる これをimproved estimatorと呼ぶ

Slide 60

Slide 60 text

60 79 サイコロの場合 母集団 𝐸 ෠ 𝑋 = 7 2 𝑉 ෠ 𝑋 = 35 12 期待値 分散

Slide 61

Slide 61 text

61 79 偶数と奇数でグループ分けする 奇数グループ 偶数グループ 平均3 平均4 全体の和は取れないが、グループ内では期待値が 厳密に求まったとする 6つの事象があるサイコロ→2つの事象があるコイン

Slide 62

Slide 62 text

62 79 3 4 確率変数 ෠ 𝑋I は、確率1/2で3、確率1/2で4 コインを弾いて、表ならサイコロの奇数、裏なら 偶数が出たことにする

Slide 63

Slide 63 text

63 79 確率変数 ෠ 𝑋I は、確率1/2で3、確率1/2で4 𝐸 ෠ 𝑋I = 3 2 + 4 2 = 7 2 𝑉 ෠ 𝑋I = 1 4 × 1 2 + 1 4 × 1 2 = 1 4 期待値 分散 𝑉 ෠ 𝑋I = 1 4 < 𝑉 ෠ 𝑋 = 35 12 分散が小さくなった

Slide 64

Slide 64 text

64 79 もとの確率変数 ෠ 𝑋 = グループ分けして部分和をとった確率変数 ෠ 𝑋I = 部分和により「グループ内の期待値が十分に求められるだ けのサンプル数」を得している →実質的にサンプル数が増えている→分散が減る

Slide 65

Slide 65 text

65 79 • 全体では厳密な和を取れないが • 状態がいくつかのグループに分かれ • その中の部分和は厳密に取れるなら • グループに関するサンプリングにより • エラーバーが小さくなる そんな都合の良い話があるか? →グラフ表現によるImproved Estimator

Slide 66

Slide 66 text

66 79 𝑚 = 1 𝑁 ෍ 𝑖 𝜎𝑖 自発磁化 𝑁はスピン数 低温 高温 𝑚2 > 0 𝑚2 ∼ 0 自発磁化の二乗𝑚2はこの系の秩序変数となる

Slide 67

Slide 67 text

67 79 𝑚2 = 𝑍−1 ෍ 𝜎𝑖 1 𝑁 ෍ 𝑖 𝜎𝑖 2 𝑤({𝜎𝑖 }) 自発磁化の期待値の表式 これのグラフ表現を考える

Slide 68

Slide 68 text

68 79 クラスターを「グループ」だと思う クラスターに通し番号𝑘をつける クラスター内のスピンは全て同じ向き → ෤ 𝜎𝑘 とする クラスター内のスピン数を𝑛𝑘 とする 𝑚 = 1 𝑁 ෍ 𝑖 𝜎𝑖 = 1 𝑁 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 スピン個別の和を クラスターごとの和にまとめる

Slide 69

Slide 69 text

69 79 𝑚2 = 𝑍−1 𝑁2 ෍ 𝑔 𝑤 𝑔 ෍ ෥ 𝜎𝑘 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 2 この部分和が厳密に計算できる 𝑚 = 1 𝑁 ෍ 𝑖 𝜎𝑖 = 1 𝑁 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 を 𝑚^2 の式へ代入する

Slide 70

Slide 70 text

70 79 𝑚2 = 𝑍−1 𝑁2 ෍ 𝑔 𝑤 𝑔 ෍ ෥ 𝜎𝑘 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 2 ここでボンド状況を固定 各クラスターのスピン状態は 全て同じ重み

Slide 71

Slide 71 text

71 79 クラスター1 𝑛1 = 5 クラスター2 𝑛2 = 4 ෍ ෥ 𝜎𝑘 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 2 = ෍ {෥ 𝜎𝑘} 𝑛1 2 ෤ 𝜎1 2 + 𝑛1 𝑛2 ෤ 𝜎1 ෤ 𝜎2 + 𝑛2 2 ෤ 𝜎2 2 = ෍ 𝑘 𝑛𝑘 2 ෤ 𝜎𝑘 = ±1だから ෤ 𝜎𝑘 2 = 1 ෤ 𝜎𝑘 = ±1だから0 グラフを固定した場合のスピン状態の部分和が厳密に計算できた

Slide 72

Slide 72 text

72 79 𝑚2 = 𝑍−1 𝑁2 ෍ 𝑔 𝑤 𝑔 ෍ ෥ 𝜎𝑘 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 2 グラフ表現した物理量 部分和をとった物理量 𝑚2 = 𝑍−1 𝑁2 ෍ 𝑔 𝑤 𝑔 ෍ 𝑘 𝑛𝑘 2 系をグラフ表現した時、クラスターサイズ の二乗の平均が自発磁化の二乗になる

Slide 73

Slide 73 text

73 79 グラフの更新 スピンの更新 グラフの更新 これらのグラフは重み𝑤(𝑔)で出現する 𝑚2 = 1 𝑁2 ෍ 𝑘 𝑛𝑘 2 グラフ更新のタイミングで以下の量を計算すると、これが improved estimatorになっている

Slide 74

Slide 74 text

74 79 スピン系で秩序変数としてよく用いられるBinder比を計算する 𝑈 = 𝑚4 𝑚2 2 𝑚2 = 1 𝑁2 ෍ 𝑘 𝑛𝑘 2 高温で1、低温で3となる それぞれのモーメントのimproved estimator 𝑚4 = 1 𝑁4 ෍ 𝑘 𝑛𝑘 2 + 6 ෍ 𝑘<𝑘′ 𝑛𝑘 2𝑛 𝑘′ 2

Slide 75

Slide 75 text

75 79 低温ではどちらも高精度だが、高温で通常のestimatorの精度が悪くなる 𝑈 = 𝑚4 𝑚2 2

Slide 76

Slide 76 text

76 79 グラフ状態 クラスター数𝑁𝑐 対応するスピンの状態数= 2𝑁𝑐 ボンドがつながる確率 𝑝 = 1 − 𝑒2𝐾 高温(𝐾 → 0)で𝑝 → 0 高温ではほとんどのボンドが切れている →多数の小さいクラスターができる →部分和が厳密に取れるスピンの状態数が増える →精度がよくなる

Slide 77

Slide 77 text

77 79 𝑚2 𝑚4 2次と4次のモーメントそれぞれの推定の精度はどちらも問題ない 𝑈 = 𝑚4 𝑚2 2 しかし、高温で微小な量の比を計算する 必要があり、精度が極端に落ちる

Slide 78

Slide 78 text

78 79 • 全体の和は計算できないが、状態をグループにわ けた時に、グループ内の和は厳密に取れる場合が ある • 部分和を厳密に取った場合に対応するestimatorを improved estimatorと呼ぶ • Improved estimatorの分散は必ず小さくなる

Slide 79

Slide 79 text

79 79 • 状態の重みの総和を分配関数とよぶ • 分配関数を書き直す=見方を変える • 見方を変えると、効率がよくなる場合がある • 同じものを異なる数式で表現すると、異なる世界が 見えてくる • 何十年も前の発見を活用してブレイクスルーが生ま れることがある ※ Fortuin-Kasteleyn表現が1972年、Swendsen-Wangが1987年