モンテカルロ法(3) 発展的アルゴリズム / Simulation 04

シミュレーション工学

Transcript

1. 1 79 モンテカルロ法(3) 発展的なアルゴリズム シミュレーション工学 慶應義塾大学大学院理工学研究科基礎理工学専攻物理情報専修 渡辺宙志
2. 2 79 • モンテカルロサンプリングと分配関数の関係 を知る • アルゴリズムを工夫することで、効率的なサ ンプリングができることを知る 本講義の目的 モンテカルロ法とは • 数値計算では、何かの和や積分の推定値を 計算することが多い • 和の形を変形することで、異なるアルゴリ ズムが生まれる
3. 3 79 • 格子の各点にスピン(小さな磁石)がある • スピンは「上」と「下」の状態がある • 隣り合うスピンをつなぐ線をボンドと呼ぶ ボンドの両側の スピンの向き 同じ 逆 エネルギー −𝐽 𝐽 𝐽 > 0ならスピンは揃いたがる(強磁性的) 𝐽 < 0ならスピンは逆向きを好む(反強磁性)
4. 4 79 格子の上の全てのボンドについてエネルギー の和を取り、この系の全エネルギーとする このような模型をイジング模型(Ising Model)と呼び、 磁性体の簡単なモデルになっている 以下、強磁性 (𝐽 > 0)の場合を考える 𝐸 = −2𝐽 + 2𝐽 = 0
5. 5 79 𝑖番目のスピンの状態を𝜎𝑖 とする 𝜎𝑖 の値は1(上向き)か-1(下向き)のいずれか −𝐽 𝐽 𝜎𝑖 = 1 𝜎𝑗 = 1 𝜎𝑖 = 1 𝜎𝑗 = −1 両方まとめて−𝐽𝜎𝑖 𝜎𝑗 と書ける
6. 6 79 全系のエネルギーは以下のように書ける 𝐻 = −𝐽 ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 系の全てのボンドについて和をとるという意味 全系のエネルギーを与える量をハミルトニアンとよぶ
7. 7 79 系の状態に通し番号をつけ、𝑖番目の状態のエネルギーを𝐸𝑖 とする 状態𝑖 エネルギー 𝐸𝑖 = 4𝐽 exp(−𝛽𝐸𝑖 ) ボルツマン定数 𝛽 = 1/𝑘𝐵 𝑇 逆温度 𝑘𝐵 この状態の出現確率がボルツマン重みに比例する この時、エネルギーの温度依存性を知りたい
8. 8 79 低温 高温 スピンがそろった状態 出現確率は高いが、総数が少ない →エネルギー重視 →低温で支配的 スピンがバラバラの状態 出現確率は低いが、総数が多い →エントロピー重視 →高温で支配的
9. 9 79 エネルギーの期待値の温度依存性 𝑈 𝛽 = 𝑍−1 ෍ 𝑖 𝐸𝑖 exp −𝛽𝐸𝑖 𝑍 = ෍ 𝑖 exp −𝛽𝐸𝑖 エネルギー𝐸をとる状態の数を𝑔(𝐸)とすると 𝑈 𝛽 = 𝑍−1 න 𝐸𝑔(𝐸) exp −𝛽𝐸 𝑑𝐸 𝑍 = න 𝑔(𝐸) exp −𝛽𝐸 𝑑𝐸 状態からエネルギー の計算は簡単 𝐸 = 0 あるエネルギーをとる状態の 数の計算は大変 𝐸 = 0
10. 10 79 𝐸 = 0 𝐸 = 4J 𝐸 = −4J 𝑔 −4𝐽 = 2 𝑔 0 = 12 𝑔 4𝐽 = 2 これを一般のNで求めるのは極めて大変 モンテカルロサンプリング
11. 11 79 1. 全ての可能な状態から無作為に一つ選ぶ 2. その状態のエネルギー𝐸𝑖 と重み𝑤𝑖 を計算する 3. 1-2を繰り返し、重み付きで平均をとる スピンが𝑁個なら、状態は2𝑁個 ほとんどの状態は重みが小さいので サンプリング効率が悪い マルコフ連鎖モンテカルロ法 この中から無作為に一つ選ぶ
12. 12 79 現在の状態から遷移可能な状態を限定し、その中 から一つ無作為に選ぶ サンプリング候補の決め方 イジング模型への適用 1. 現在の状態を𝐴とし、スピンを一つ無作為に選ぶ 2. 選んだスピンを反転させた状態を提案状態𝐵とする 3. それぞれのエネルギーから遷移確率𝑃(𝐴 → 𝐵)を計算 し、遷移させるか決める 4. 遷移しなかった場合は状態はそのまま、遷移した場 合は提案状態を現在の状態にして1.へ
13. 13 79 現在の状態 遷移可能な状態 スピンが𝑁個なら総状態数は2𝑁、遷移可能な状態は𝑁 この中から無作為に一つ選ぶ 遷移可能な状態が限定された
14. 14 79 現在の状態𝐴 現在の状態𝐵 𝑃(𝐴 → 𝐵) 𝑃(𝐵 → 𝐴) 𝑃(𝐵 → 𝐵) 𝑃(𝐴 → 𝐴) 𝐸𝐴 = −4𝐽 𝐸𝐵 = 0 𝑤𝐴 = exp −𝛽𝐸𝐴 𝐸𝐴 < 𝐸𝐵 なので𝑤𝐴 > 𝑤𝐵 𝑤𝐵 = exp −𝛽𝐸𝐵
15. 15 79 𝑃(𝐴 → 𝐵) 𝑃(𝐵 → 𝐴) 𝑃(𝐵 → 𝐵) 𝑃(𝐴 → 𝐴) 確率の保存条件 𝑃 𝐴 → 𝐴 + 𝑃 𝐴 → 𝐵 = 1 𝑃 𝐵 → 𝐵 + 𝑃 𝐵 → 𝐴 = 1 詳細釣り合い条件 𝑃 𝐴 → 𝐵 𝑃 𝐵 → 𝐴 = 𝑤𝐵 𝑤𝐴
16. 16 79 メトロポリス法(Metropolis method) 重みが大きくなる場合は必ず遷移、そうでない場合は確率的に遷移させる 𝑃 𝐵 → 𝐴 = 1 𝑃 𝐴 → 𝐵 = 𝑤𝐵 𝑤𝐴 = exp −𝛽𝐸𝐵 exp(−𝛽𝐸𝐴 ) = exp(−𝛽Δ𝐸) Δ𝐸 = 𝐸𝐵 − 𝐸𝐴 > 0なので、 exp −𝛽Δ𝐸 < 1 提案状態のエネルギーが、現在の状態より低ければ必ず遷移 高ければ確率𝑝 = exp −𝛽Δ𝐸 < 1で遷移
17. 17 79 イジング模型にマルコフ連鎖モンテカルロ法を適用した場 合のアルゴリズム(メトロポリス法を採用した場合) 1. スピンを一つ無作為に選ぶ 2. 選んだスピンを反転させた状態を提案状態とする 3. 現在の状態とのエネルギー差Δ𝐸を計算する 4. エネルギーが下がる場合(Δ𝐸 < 0)なら必ず遷移。そ うでなければ確率𝑝 = exp(−𝛽Δ𝐸)で遷移 5. 1-4を繰り返す 一度に一つのスピンだけひっくり返すので Single-Spin-Flip Algorithmと呼ぶ
18. 18 79 状態に重みがある場合の期待値を厳密に計算することは難しい その理由は同じ重みを持つ状態の数がわからないから →サンプリングによる期待値の推定 単純サンプリング 全ての状態から無作為に状態をサンプリングする 重みの小さい状態ばかりサンプリングするため非効率 マルコフ連鎖モンテカルロ法 現在の状態をもとに、サンプリングする状態を限定する 現在の状態と提案状態の重みを比較し、適切な遷移確率を 作ることで重みに比例したサンプリングができる (重み付きサンプリング)→効率的
19. 19 79 一度に一つしかスピンをひっくり返せないのは非効率的？ 複数のスピンを同時にひっくり返せたら計算が加速 するのでは？
20. 20 79 同じ向きのスピンの数 𝑛+ = 3 逆向きのスピンの数 𝑛− = 1 同じ向きのスピンの数 𝑛+ = 1 逆向きのスピンの数 𝑛− = 3 ひっくり返る 𝐸𝐴 = 𝑛− − 𝑛+ 𝐽 𝐸𝐵 = 𝑛+ − 𝑛− 𝐽 現在の状態𝐴 提案状態𝐵 Δ𝐸 = 𝐸𝐵 − 𝐸𝐴 = 2 𝑛+ − 𝑛− 𝐽 エネルギー差が同じ向きのスピンの数を数えるだけで計算できる
21. 21 79 現在の状態𝐴 提案状態𝐵 この二つを同時に ひっくり返そうとする あるスピンの状態が確定していないと、隣のスピンを ひっくり返す試行ができない ？ このスピンの状態が 確定しないと このスピンをひっくり返す ための重みが計算できない
22. 22 79 ボンドで接続されていないスピンは影響を受けない あるスピンがひっくり返ると、ボンドで繋がったスピン の重みが変わる このスピンがひっくり返ると このスピンの重みに影響がでる このスピンがひっくり返っても このスピンの重みには影響しない
23. 23 79 これらのスピンを同時にひっくり返そうとする ？ このスピンをひっくり返す 試行の結果はまだ知らない わからないまま確率を計算して ひっくり返してしまう ？
24. 24 79 もしスピンがひっくり返っていなければ？ 一つスピンをひっくり 返しただけ もしスピンがひっくり返っていたら？ 二つのステップを 逐次的に行ったのと等価 これらの状態間で詳細釣り合い条件が成り立っているので 正しいサンプリングができている
25. 25 79 ボンドで繋がっていない二つのグループに分け、 Aサイト、Bサイトと名前をつける A A B B Bサイトのスピンを固定し、 Aサイトのスピンを全て独立 にひっくり返す試行を行う Aサイトのスピンを固定し、 Bサイトのスピンを全て独立 にひっくり返す試行を行う
26. 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 他のスレッドの結果を待たずに自分の担当のスピンをひっくり返してよい
27. 27 79 16個の状態を 4x4個の状態に分離した x =
28. 28 79 𝑍 = ෍ {𝜎} exp(−𝛽𝐻) = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp(−𝛽𝐻) 全てのスピンの状態についての和 2𝑁通り A、Bサイトそれぞれのスピンの状態に ついての和 2𝑁/2通り 𝐻 = −𝐽 ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 ハミルトニアンを代入する 𝑍 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗
29. 29 79 𝑍 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 外側を固定して 内側の和についてのみサンプリングする 𝜎𝑖 𝜎𝑗1 + 𝜎𝑖 𝜎𝑗2 + 𝜎𝑖 𝜎𝑗3 + 𝜎𝑖 𝜎𝑗4 = ℎ𝑖 𝜎𝑖 𝑖 𝑗3 𝑗2 𝑗4 𝑗1 外側のスピンが固定されているので、 それをまとめてしまう ℎ𝑖 ≡ 𝜎𝑗1 + 𝜎𝑗2 + 𝜎𝑗3 + 𝜎𝑗4
30. 30 79 𝑍 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖 ℎ𝑖 𝜎𝑖 𝑍 = ෍ {𝜎} exp(−𝛽𝐻) = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp(−𝛽𝐻) だったから 𝐻 = −𝐽 ෍ 𝑖 ℎ𝑖 𝜎𝑖 というハミルトニアンを計算している気持ちになった ローカルな磁場下にある独立なスピンの集まり → 相互作用が消えてしまった → 相互作用がないので好き勝手できる(例えば並列化)
31. 31 79 𝑍 = ෍ {𝜎𝐴} ෍ {𝜎𝐵} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 = ෍ {𝜎𝐵} ෍ {𝜎𝐴} exp ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 和の順序を入れ替えれば、Aサイトについて同じことができる アルゴリズム 1. Aサイトのスピンを固定し、Bサイトのスピンだけ更新する 2. Bサイトのスピンを固定し、Aサイトのスピンだけ更新する 3. 以下繰り返し 計算の手間は変わらないが、更新しようとするスピンが全て独立で あるために並列化が可能になる
32. 32 79 • 分配関数の和の取り方を工夫すると、もとの系とは異な る形のハミルトニアンが現れる→異なる物理系をシミュ レーションしているかのように見える • 多重和の外側を固定し、内側だけサンプリングすると、 スピンが独立になるために計算が楽になることがある もっともうまく分配関数を変形すれば、より効率 的なサンプリングが可能になるのでは？ クラスターアルゴリズム
33. 33 79 𝑍 = ෍ {𝜎} exp(−𝛽𝐻) 𝐻 = −𝐽 ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 イジング模型のハミルトニアン その分配関数 𝑍 = ෍ {𝜎} exp ෍ 𝑖,𝑗 𝐾𝜎𝑖 𝜎𝑗 = ෍ {𝜎} ෑ 𝑖,𝑗 exp 𝐾𝜎𝑖 𝜎𝑗 ハミルトニアンを代入した 𝐾 ≡ 𝛽𝐽 指数関数の肩の和を外に出して 積に書き換えた これを書き直す
34. 34 79 exp 𝐾𝜎𝑖 𝜎𝑗 = ൝ 𝑒𝐾 𝜎𝑖 = 𝜎𝑗 𝑒−𝐾 𝜎𝑖 ≠ 𝜎𝑗 であるから exp 𝐾𝜎𝑖 𝜎𝑗 = 𝑒−𝐾 + 𝛿𝜎𝑖,𝜎𝑗 𝑒𝐾 − 𝑒−𝐾 = 𝑒𝐾 𝑒−2𝐾 + 𝛿𝜎𝑖,𝜎𝑗 1 − 𝑒−2𝐾 = 𝑒𝐾 (1 − 𝑝) + 𝑝𝛿𝜎𝑖,𝜎𝑗 ≡ 𝑝 ここまではただの等式変形
35. 35 79 以下のデルタ関数を使った恒等式を使って書き直す 𝑥 + 𝑦 = ෍ 𝑏=0 1 𝑥𝛿𝑏,0 + 𝑦𝛿𝑏,1 𝑒𝐾 (1 − 𝑝) + 𝑝𝛿𝜎𝑖,𝜎𝑗 = ෍ 𝑏𝑖𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 新たに導入された𝑏𝑖𝑗 はスピン𝑖とスピン𝑗を結ぶボンド上で 定義されており、ボンド自由度と呼ぶ
36. 36 79 𝑏𝑖𝑗 はスピンiとjをつなぐボンド上で定義される 値として1か0をとる 1の時にactive、0の時にinactiveであると定義する ボンド変数 スピン変数 𝜎𝑖 はサイトi上で定義される 値として1か-1をとる 1の時に上向き、-1の時に下向きと定義する 𝜎𝑖 𝜎𝑗 𝑏𝑖𝑗
37. 37 79 exp(𝐾𝜎𝑖 𝜎𝑗 ) = ෍ 𝑏𝑖𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 = + 新たにボンド自由度を導入し、スピンの相互作用の 重みをactiveなボンドとinactiveなボンドの和に分けた inactive active
38. 38 79 以上の式変形をまとめると 𝑍 = ෍ {𝜎} exp(−𝛽𝐻) = ෍ {𝜎𝑖} ෍ {𝑏𝑖𝑗} ෑ 𝑖,𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 隣り合うスピンに関する積 ボンドに関する和 スピンに関する和 等価変形なのにボンドに関する自由度が増えた
39. 39 79 𝑍 = ෍ {𝜎𝑖} ෍ {𝑏𝑖𝑗} ෑ 𝑖,𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 ≡ 𝑆𝑖𝑗 • 新たにボンド自由度がうまれ𝑏𝑖𝑗 = 0か𝑏𝑖𝑗 = 1の値をとる • 分配関数が多重和で書かれているので、スピンに関するサ ンプリングとボンドに関するサンプリングを交互にとれる • ボンド𝑏𝑖𝑗 がactive(𝑏𝑖𝑗 = 1)であるか、inactiveであるかの確 率は、スピン状態に依存する重みで決まる スピン状態を固定した状態で ボンド状態のサンプリングを考える
40. 40 79 両側のスピンが逆向き(𝜎𝑖 ≠ 𝜎𝑗)の時 𝑆𝑖𝑗 = 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 𝑆𝑖𝑗 = 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 𝑏𝑖𝑗 = 0の場合しか重みを持たない → 必ずinactive 両側のスピンが同じ向き(𝜎𝑖 ≠ 𝜎𝑗)の時 𝑆𝑖𝑗 = 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝑏𝑖𝑗,1 確率1 − 𝑝でinactive、確率𝑝でactive スピン変数を固定
41. 41 79 現在の状態 ボンド状態を忘れる 新しいボンド状態を与える 隣り合うスピンが同じ向きの場合 • 確率𝑝でactive • 確率1 − 𝑝でinactive 隣り合うスピンが逆向きの場合 • 必ずinactive 𝑝 ≡ 1 − 𝑒−2𝐾
42. 42 79 𝑆𝑖𝑗 = 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 ボンド変数を固定 ボンドがactive(𝑏𝑖𝑗 = 1)の時 𝑆𝑖𝑗 = 𝑒𝐾𝑝𝛿𝜎𝑖,𝜎𝑗 ボンドがinactive(𝑏𝑖𝑗 = 0)の時 𝑆𝑖𝑗 = 𝑒𝐾(1 − 𝑝) スピンの相互作用が含まれない → スピンの相互作用が切れている
43. 43 79 activeなボンドでつながっているスピンをクラスターと呼ぶ
44. 44 79 inactiveなボンドでつながるスピンは相互作用を持たない → activeなボンドでつながるスピンを一斉にひっくり返しても 状態の重みが変わらない 状態𝐴 状態𝐵 𝑤𝐴 = 𝑤𝐵
45. 45 79 𝑃(𝐴 → 𝐵) 𝑃(𝐵 → 𝐴) 熱浴法を採用すると 𝑃 𝐴 → 𝐵 = 𝑤𝐵 𝑤𝐴 + 𝑤𝐵 = 1 2 𝑃 𝐵 → 𝐴 = 𝑤𝐴 𝑤𝐴 + 𝑤𝐵 = 1 2 状態𝐴 状態𝐵
46. 46 79 𝑍 = ෍ {𝑏𝑖𝑗} ෍ {𝜎𝑖} ෑ 𝑖,𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 ボンド状態を固定した状態で スピン状態のサンプリングを考える activeなボンドでつながるスピンを一斉にひっくり返しても重みが等しい →同じクラスターに所属するスピンを全てランダムにひっくり返してよい クラスターごとに、元のスピン配置を忘れて ランダムにスピンを割り当ててよい ※同じクラスターに属すスピンは必ず同じ向き
47. 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)
48. 48 79 R. H. Swendsen N. Ito H. Watanabe
49. 49 79 𝑍 = ෍ {𝑏𝑖𝑗} ෍ {𝜎𝑖} ෑ 𝑖,𝑗 𝑒𝐾 1 − 𝑝 𝛿𝑏𝑖𝑗,0 + 𝑝𝛿𝜎𝑖,𝜎𝑗 𝛿𝑏𝑖𝑗,1 外側にボンド状態に関する和がある ボンド状態は、たとえば左図のよう なグラフを表している →あらゆるグラフに関する和 • このグラフ表現をFortuin-Kasteleyn表現と呼ぶ • グラフ表現を利用してスピンを更新するアルゴリズムを クラスターアルゴリズムと呼ぶ • Swendsen-Wangはクラスターアルゴリズムの一種 C.M.Fortuin, P.W.Kasteleyn, Physica, vol. 57, pp 536 (1972).
50. 50 79 クラスターアルゴリズムは非常に早い しかし、早い理由はスピンをまとめて更新するからではない 低温である時、以下のようなドメインウォールができたとする 低温だから全てのスピンがそろった状態になりたい
51. 51 79 Single-Spin-Flipでドメインウォール解消を試みる 途中で必ずエネルギーが高い状態を経由する必要がある 低温であるほど、高エネルギー状態は出現しづらい →極めて緩和が遅くなる エネルギーの損 4J エネルギーの損 6J エネルギーの損 0
52. 52 79 Swendsen-Wangの場合 現在の状態 ボンド状態を更新して スピンを忘れる 確率1/2で全ての スピンがそろう クラスターアルゴリズムでは、ドメインウォールを非常 に早く解消できる→緩和が本質的に早い
53. 53 79 • 分配関数を式変形し、新たにボンドの自由度を追加すると グラフの和の形に書ける(Fortuin-Kasteleyn表現) • グラフ表現を利用したサンプリングアルゴリズムをクラス ターアルゴリズムと呼ぶ • Inactiveなボンドでつながったスピンは重みがスピン状態に 依存しない→実質的に相互作用が切れている • スピンをお互いに相互作用しないクラスターに分けること ができる • クラスターはそれぞれ独立にスピン状態を与えられる • ボンド変数とスピンを交互に更新することで、効率的なサ ンプリングができる • Swendseng-Wangはクラスターアルゴリズムの一種で、他 にもWolffのアルゴリズムなどがある
54. 54 79 • 状態𝑖の出現確率が𝑤𝑖 に比例する時、𝑤𝑖 を状態𝑖の重みと呼ぶ • 重みの総和𝑍 = σ𝑖 𝑤𝑖 を分配関数と呼ぶ • マルコフ連鎖モンテカルロ法とは分配関数をサンプリングに より推定する手法 • 分配関数を変形すると、異なる物理系のように見える • 同じ分配関数であるにも関わらず、変形をすることで並列化 できたり緩和を加速できたりする Fortuin-Kasteleyn表現は1972年には知られていた。 しかしSwendsen-Wangアルゴリズムの提案は1987年 まだ見ぬ新たなブレイクスルーがあるかも・・・？
55. 55 79 確率変数 ෠ 𝑋の母集団のパラメータを知りたい 𝐸 ෠ 𝑋 = ෍ 𝑘 𝑌𝑘 𝑝𝑘 ෠ 𝑋の取りうる値𝑌𝑘 とその実現確率𝑝𝑘 が既知なら 母集団の期待値の厳密な表記は 一般にはこの和を厳密に計算できない →サンプリングで推定
56. 56 79 確率変数 ෠ 𝑋の母集団のパラメータを知りたい 確率変数 ෠ 𝑋を𝑁回観測して標本 {𝑋𝑖 }を得た ここから ෠ 𝑋の母集団のパラメータを推定したい 𝐸 ෠ 𝑋 ∼ ത 𝑋 = 1 𝑁 ෍ 𝑖 𝑋𝑖 母集団の期待値 期待値のestimator 知りたい量を標本で表す関数をestimatorと呼ぶ
57. 57 79 𝑁が大きい時、中心極限定理により、 ത 𝑋の分布はガウス分布に近づく ത 𝑋の分散のestimator 𝑉 ത 𝑋 = 1 𝑁 𝑉 ෠ 𝑋 = 1 𝑁(𝑁 − 1) ෍ 𝑖 𝑋𝑖 − ത 𝑋 2 ത 𝑋の分散 ෠ 𝑋の分散(母分散) これが統計誤差だから、なるべく小さくしたい
58. 58 79 𝐸 ෠ 𝑋 = ෍ 𝑘 𝑌𝑘 𝑝𝑘 全体では厳密な和を取れないが、状態がいくつか のグループに分かれ、その中の部分和は厳密に取 れるとする ここは計算できないが ここは厳密に計算できる = ෍ 𝑔 ෍ 𝑘∈𝑔 𝑌𝑘 𝑝𝑘 = ෍ 𝑔 ത 𝑌 𝑔 𝑝𝑔 すると、対応するestimatorの分散が小さくなる これをimproved estimatorと呼ぶ
59. 59 79 𝐸 ෠ 𝑋 = ෍ 𝑘 𝑌𝑘 𝑝𝑘 = ෍ 𝑔 ෍ 𝑘∈𝑔 𝑌𝑘 𝑝𝑘 = ෍ 𝑔 ത 𝑌 𝑔 𝑝𝑔 確率𝑝𝑔 で状態𝑔となる確率変数 ෠ 𝑋I を考えると もとの確率変数 ෠ 𝑋よりも分散が小さくなる 𝑉 ෠ 𝑋 > 𝑉 ෠ 𝑋I 対応するestimatorの分散も小さくなる これをimproved estimatorと呼ぶ
60. 60 79 サイコロの場合 母集団 𝐸 ෠ 𝑋 = 7 2 𝑉 ෠ 𝑋 = 35 12 期待値 分散
61. 61 79 偶数と奇数でグループ分けする 奇数グループ 偶数グループ 平均3 平均4 全体の和は取れないが、グループ内では期待値が 厳密に求まったとする 6つの事象があるサイコロ→2つの事象があるコイン
62. 62 79 3 4 確率変数 ෠ 𝑋I は、確率1/2で3、確率1/2で4 コインを弾いて、表ならサイコロの奇数、裏なら 偶数が出たことにする
63. 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 分散が小さくなった
64. 64 79 もとの確率変数 ෠ 𝑋 = グループ分けして部分和をとった確率変数 ෠ 𝑋I = 部分和により「グループ内の期待値が十分に求められるだ けのサンプル数」を得している →実質的にサンプル数が増えている→分散が減る
65. 65 79 • 全体では厳密な和を取れないが • 状態がいくつかのグループに分かれ • その中の部分和は厳密に取れるなら • グループに関するサンプリングにより • エラーバーが小さくなる そんな都合の良い話があるか？ →グラフ表現によるImproved Estimator
66. 66 79 𝑚 = 1 𝑁 ෍ 𝑖 𝜎𝑖 自発磁化 𝑁はスピン数 低温 高温 𝑚2 > 0 𝑚2 ∼ 0 自発磁化の二乗𝑚2はこの系の秩序変数となる
67. 67 79 𝑚2 = 𝑍−1 ෍ 𝜎𝑖 1 𝑁 ෍ 𝑖 𝜎𝑖 2 𝑤({𝜎𝑖 }) 自発磁化の期待値の表式 これのグラフ表現を考える
68. 68 79 クラスターを「グループ」だと思う クラスターに通し番号𝑘をつける クラスター内のスピンは全て同じ向き → ෤ 𝜎𝑘 とする クラスター内のスピン数を𝑛𝑘 とする 𝑚 = 1 𝑁 ෍ 𝑖 𝜎𝑖 = 1 𝑁 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 スピン個別の和を クラスターごとの和にまとめる
69. 69 79 𝑚2 = 𝑍−1 𝑁2 ෍ 𝑔 𝑤 𝑔 ෍ ෥ 𝜎𝑘 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 2 この部分和が厳密に計算できる 𝑚 = 1 𝑁 ෍ 𝑖 𝜎𝑖 = 1 𝑁 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 を 𝑚^2 の式へ代入する
70. 70 79 𝑚2 = 𝑍−1 𝑁2 ෍ 𝑔 𝑤 𝑔 ෍ ෥ 𝜎𝑘 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 2 ここでボンド状況を固定 各クラスターのスピン状態は 全て同じ重み
71. 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 グラフを固定した場合のスピン状態の部分和が厳密に計算できた
72. 72 79 𝑚2 = 𝑍−1 𝑁2 ෍ 𝑔 𝑤 𝑔 ෍ ෥ 𝜎𝑘 ෍ 𝑘 𝑛𝑘 ෤ 𝜎𝑘 2 グラフ表現した物理量 部分和をとった物理量 𝑚2 = 𝑍−1 𝑁2 ෍ 𝑔 𝑤 𝑔 ෍ 𝑘 𝑛𝑘 2 系をグラフ表現した時、クラスターサイズ の二乗の平均が自発磁化の二乗になる
73. 73 79 グラフの更新 スピンの更新 グラフの更新 これらのグラフは重み𝑤(𝑔)で出現する 𝑚2 = 1 𝑁2 ෍ 𝑘 𝑛𝑘 2 グラフ更新のタイミングで以下の量を計算すると、これが improved estimatorになっている
74. 74 79 スピン系で秩序変数としてよく用いられるBinder比を計算する 𝑈 = 𝑚4 𝑚2 2 𝑚2 = 1 𝑁2 ෍ 𝑘 𝑛𝑘 2 高温で1、低温で3となる それぞれのモーメントのimproved estimator 𝑚4 = 1 𝑁4 ෍ 𝑘 𝑛𝑘 2 + 6 ෍ 𝑘<𝑘′ 𝑛𝑘 2𝑛 𝑘′ 2
75. 75 79 低温ではどちらも高精度だが、高温で通常のestimatorの精度が悪くなる 𝑈 = 𝑚4 𝑚2 2
76. 76 79 グラフ状態 クラスター数𝑁𝑐 対応するスピンの状態数= 2𝑁𝑐 ボンドがつながる確率 𝑝 = 1 − 𝑒2𝐾 高温(𝐾 → 0)で𝑝 → 0 高温ではほとんどのボンドが切れている →多数の小さいクラスターができる →部分和が厳密に取れるスピンの状態数が増える →精度がよくなる
77. 77 79 𝑚2 𝑚4 2次と4次のモーメントそれぞれの推定の精度はどちらも問題ない 𝑈 = 𝑚4 𝑚2 2 しかし、高温で微小な量の比を計算する 必要があり、精度が極端に落ちる
78. 78 79 • 全体の和は計算できないが、状態をグループにわ けた時に、グループ内の和は厳密に取れる場合が ある • 部分和を厳密に取った場合に対応するestimatorを improved estimatorと呼ぶ • Improved estimatorの分散は必ず小さくなる
79. 79 79 • 状態の重みの総和を分配関数とよぶ • 分配関数を書き直す＝見方を変える • 見方を変えると、効率がよくなる場合がある • 同じものを異なる数式で表現すると、異なる世界が 見えてくる • 何十年も前の発見を活用してブレイクスルーが生ま れることがある ※ Fortuin-Kasteleyn表現が1972年、Swendsen-Wangが1987年