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

シミュレーション工学

Transcript

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

   View Slide

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

   View Slide

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

   View Slide

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

   View Slide

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

   View Slide

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

   View Slide

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

   View Slide

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

   View Slide

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

   View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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
    他のスレッドの結果を待たずに自分の担当のスピンをひっくり返してよい
    

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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)
    

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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
    分散が小さくなった
    

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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
    グラフを固定した場合のスピン状態の部分和が厳密に計算できた
    

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide