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

シミュレーション工学

Transcript

1. 1 82 モンテカルロ法(3) マルコフ遷移行列 シミュレーション工学 慶應義塾大学大学院理工学研究科基礎理工学専攻物理情報専修 渡辺宙志
2. 2 82 はじめに • 状態遷移図を理解し、描けるようになる • マルコフ遷移行列の固有値、固有ベクトルの 意味を理解する 本講義の目的 モンテカルロ法とは • 数値計算では、何かの和や積分の推定値を 計算することが多い • 和の形を変形することで、異なるアルゴリ ズムが生まれる
3. 3 82 状態遷移図 • 3つのマスがある「すごろく」を考える • マスをそれぞれ「マス1」「マス2」「マス3」と名前をつける サイコロを振って・・・ • マス1にいるとき • 1が出たらそのまま • 2,3が出たらマス2へ • 4,5,6が出たらマス3へ • マス2にいるとき • 1,2が出たらそのまま • 3が出たらマス1へ • 4,5,6が出たらマス3へ • マス3にいるとき • 1,2,3が出たらそのまま • 4が出たらマス1へ • 5,6が出たらマス2へ 1 2 3
4. 4 82 状態遷移図 • 3つのマスがある「すごろく」を考える • マスをそれぞれ「マス1」「マス2」「マス3」と名前をつける サイコロを振って・・・ • マス1にいるとき • 1/6の確率でそのまま • 1/3の確率でマス2へ • 1/2の確率でマス3へ • マス2にいるとき • 1/3の確率でそのまま • 1/6の確率でマス1へ • 1/2の確率でマス3へ • マス3にいるとき • 1/2の確率でそのまま • 1/6の確率でマス1へ • 1/3の確率でマス2へ 1 2 3 1/6 1/3 1/2
5. 5 82 状態遷移図 ある状態から、次のステップにどんな確率でどこに行くかを 表した図を状態遷移図と呼ぶ 1 2 3 1/6 1/3 1/2 1/3 1/2 1/6 1/2 1/6 1/3
6. 6 82 マルコフ連鎖 サイコロを振って・・・ • マス1にいるとき • 1/6の確率でそのまま • 1/3の確率でマス2へ • 1/2の確率でマス3へ 2 3 1 2 2 2 1 1 ・・・ ・・・ どのような履歴をたどっていたとしても、次の状態が現在の状態だけで決まる この性質をマルコフ性(Markov property)と呼ぶ この性質を持つ状態遷移の連鎖をマルコフ連鎖(Markov chain)と呼ぶ
7. 7 82 マルコフ行列 以下の状態遷移図において、以下のようなことを知りたい • 最初に状態1にいたとき、tステップ後に各状態にいる確率 • 十分に時間が経過したとき、各状態にいる確率 1 2 3 1/6 1/3 1/2 1/3 1/2 1/6 1/2 1/6 1/3
8. 8 82 マルコフ行列 𝑝𝑖 (𝑡) 時刻𝑡において、状態𝑖にいる確率 𝑝1 𝑡 + 1 = 1 6 𝑝1 𝑡 + 1 6 𝑝2 𝑡 + 1 6 𝑝3 (𝑡) 𝑝2 𝑡 + 1 = 1 3 𝑝1 𝑡 + 1 3 𝑝2 𝑡 + 1 3 𝑝3 (𝑡) 𝑝3 𝑡 + 1 = 1 2 𝑝1 𝑡 + 1 2 𝑝2 𝑡 + 1 2 𝑝3 (𝑡) 行列とベクトルの形にかけそう 1 2 3 1/6 1/3 1/2 1/3 1/2 1/6 1/2 1/6 1/3
9. 9 82 マルコフ行列 Ԧ 𝑝 𝑡 = 𝑝1 𝑡 𝑝2 (𝑡) 𝑝3 (𝑡) 時刻tにおける状態ベクトル Ԧ 𝑝 𝑡 + 1 = 𝑀 Ԧ 𝑝(𝑡) 𝑀 = 1/6 1/6 1/6 1/3 1/3 1/3 1/2 1/2 1/2 ただし この行列をマルコフ行列(もしくは遷移確率行列)と呼ぶ
10. 10 82 マルコフ行列の性質 • すべての要素が0から1の間の実数(確率だから) • 各列の値の総和は1(確率の保存則) 𝑀 = 1/6 1/6 1/6 1/3 1/3 1/3 1/2 1/2 1/2 • 状態ベクトルにかけると次の状態ベクトルが得られる Ԧ 𝑝 𝑡 + 1 = 𝑀 Ԧ 𝑝(𝑡) • n回かけるとnステップ後の状態ベクトルが得られる Ԧ 𝑝 𝑡 + 𝑛 = 𝑀𝑛 Ԧ 𝑝(𝑡) • 最大固有値は1であり、対応する固有ベクトルが定常状態 ※マルコフ過程が非周期的、正再帰的、既約である場合
11. 11 82 マルコフ行列 最初に状態1にいたとき、tステップ後に各状態にいる確率 Ԧ 𝑝 0 = 1 0 0 として Ԧ 𝑝 𝑡 = 𝑀𝑡 Ԧ 𝑝(0) を計算すれば良い 行列のベキ乗を求めるには対角化する 𝑀 = 𝑃−1𝐴𝑃 𝑀𝑡 = 𝑃−1𝐴𝑃 𝑡 = 𝑃−1𝐴𝑃𝑃−1𝐴𝑃 … 𝐴𝑃 = 𝑃−1𝐴𝑡𝑃
12. 12 82 マルコフ行列 十分に時間が経過したとき、各状態にいる確率 Ԧ 𝑝(∞) = 𝑀∞ Ԧ 𝑝(0) を求めたい マルコフ行列の固有値を𝜆𝑖 、対応する固有ベクトルを Ԧ 𝑣𝑖 とする マルコフ行列の最大固有値の絶対値は1であるから、 1 = |𝜆1 | > 𝜆2 > |𝜆3 | とすると 𝑀 Ԧ 𝑣1 = Ԧ 𝑣1 𝑀 Ԧ 𝑣𝑖 = 𝜆𝑖 Ԧ 𝑣𝑖 定常状態が存在するなら Ԧ 𝑝(∞) = 𝑀 Ԧ 𝑝(∞) であるから Ԧ 𝑝(∞) = Ԧ 𝑣1 定常状態とは、マルコフ行列の最大固有値に対応する固有ベクトル
13. 13 82 詳細釣り合い条件 定常状態はマルコフ行列の最大固有状態だが、行列の対角化 をせずに定常状態を求めたい Ԧ 𝑝(∞) = 𝜋1 𝜋2 𝜋3 𝜋𝑖 定常状態において状態𝑖にいる確率 ここで、マルコフ遷移図のある2状態の遷移に注目する 1 2 3 1/6 1/3 1/2 1/3 1/2 1/6 1/2 1/6 1/3
14. 14 82 詳細釣り合い条件 1 2 1/3 1/6 • 状態1から2へは確率1/3で遷移する • 状態2から1へは確率1/6で遷移する 定常状態なら、上記のやりとりで確率が変わらない 1 3 𝜋1 = 1 6 𝜋2 状態1から2に行く流れ 状態2から1に行く流れ 𝜋1 : 𝜋2 = 1: 2
15. 15 82 詳細釣り合い条件 2 3 1/2 1/3 • 状態2から3へは確率1/2で遷移する • 状態3から2へは確率1/3で遷移する 1 2 𝜋2 = 1 3 𝜋3 状態2から3に行く流れ 状態3から2に行く流れ 𝜋2 : 𝜋3 = 2: 3
16. 16 82 詳細釣り合い条件 以上から 𝜋1 : 𝜋2 : 𝜋3 = 1: 2: 3 Ԧ 𝑝(∞) = 𝜋1 𝜋2 𝜋3 であり 𝜋1 + 𝜋2 + 𝜋3 = 1 だから Ԧ 𝑝(∞) = 1/6 2/6 3/6 定常状態が、マルコフ行列の固有値、固有ベクトルを求めずに決まった
17. 17 82 詳細釣り合い条件 𝑖 𝑗 𝑃(𝑖 → 𝑗) 𝑃(𝑗 → 𝑖) 状態𝑖の定常状態の確率を𝜋𝑖 として、任意の2状態𝑖, 𝑗に ついて以下が成立することを詳細釣り合い条件と呼ぶ 𝜋𝑖 𝑃 𝑖 → 𝑗 = 𝜋𝑗 𝑃(𝑗 → 𝑖) 遷移確率が決まると、定常分布の比が決まる 𝜋𝑖 𝜋𝑗
18. 18 82 ここまでのまとめ • 履歴によらず、次の状態が現在の状態だけ決まること マルコフ性 • マルコフ性を持つ離散的な確率過程のこと マルコフ連鎖 詳細釣り合い条件 • 任意の二状態において、定常分布の比(重みの比) と遷移確率の比が満たす性質 • 遷移確率が決まると、定常分布の比が決まる
19. 19 82 重みと遷移確率 1 2 3 1/6 1/3 1/2 1/3 1/2 1/6 1/2 1/6 1/3 遷移確率が決まると、定常状態の分布の比が決まる 1 2 3 𝜋1 𝜋2 𝜋3 𝜋1 : 𝜋2 : 𝜋3 = 1: 2: 3 逆に、希望する定常状態が実現するように遷移確率を決められないだろうか？ 1 2 3 𝜋1 𝜋2 𝜋3 𝜋1 : 𝜋2 : 𝜋3 = 1: 2: 3 1 2 3 ??? ??? ??? ??? ??? ??? ??? ??? ???
20. 20 82 詳細釣り合い条件 𝑖 𝑗 𝑃(𝑖 → 𝑗) 𝑃(𝑗 → 𝑖) 一般に、定常分布の確率𝜋𝑖 は求められないが、それに 比例する重み𝑤𝑖 はわかる 𝜋𝑖 𝑃 𝑖 → 𝑗 = 𝜋𝑗 𝑃(𝑗 → 𝑖) 𝑤𝑖 ∝ 𝜋𝑖 𝑤𝑗 ∝ 𝜋𝑗 𝑃 𝑗 → 𝑖 𝑃 𝑖 → 𝑗 = 𝜋𝑖 𝜋𝑗 = 𝑤𝑖 𝑤𝑗 上記を満たすように遷移確率を決めれば良い
21. 21 82 詳細釣り合い条件 𝑖 𝑗 𝑃(𝑖 → 𝑗) 𝑃(𝑗 → 𝑖) 𝑤𝑖 ∝ 𝜋𝑖 𝑤𝑗 ∝ 𝜋𝑗 𝑃(𝑖 → 𝑖) 𝑃(𝑗 → 𝑗) 決めるべき確率は以下の4つ 𝑃(𝑖 → 𝑖), 𝑃(𝑖 → 𝑗), 𝑃(𝑗 → 𝑖), 𝑃(𝑗 → 𝑗) 確率の保存則 𝑃 𝑖 → 𝑖 + 𝑃 𝑖 → 𝑗 = 1 𝑃 𝑗 → 𝑗 + 𝑃 𝑗 → 𝑖 = 1 𝑤𝑖 𝑃 𝑖 → 𝑗 = 𝑤𝑗 𝑃(𝑗 → 𝑖) 詳細釣り合い条件 式が一本足りない
22. 22 82 確率の決め方 𝑤𝑖 𝑃 𝑖 → 𝑗 = 𝑤𝑗 𝑃(𝑗 → 𝑖) を満たすように遷移確率を決めたい 熱浴法 (heat bath method) 𝑃 𝑖 → 𝑗 = 𝑤𝑗 𝑤𝑖 + 𝑤𝑗 𝑃 𝑗 → 𝑖 = 𝑤𝑖 𝑤𝑖 + 𝑤𝑗 メトロポリス法 (Metropolis method) 𝑤𝑖 < 𝑤𝑗 である時 𝑃 𝑖 → 𝑗 = 1 𝑃 𝑗 → 𝑖 = 𝑤𝑖 𝑤𝑗 重みが増える向きには 必ず遷移する 重みが減る向きには 確率的に遷移する
23. 23 82 イジング模型 • 格子の各点にスピン(小さな磁石)がある • スピンは「上」と「下」の状態がある • 隣り合うスピンをつなぐ線をボンドと呼ぶ ボンドの両側の スピンの向き 同じ 逆 エネルギー −𝐽 𝐽 𝐽 > 0ならスピンは揃いたがる(強磁性的) 𝐽 < 0ならスピンは逆向きを好む(反強磁性)
24. 24 82 イジング模型 格子の上の全てのボンドについてエネルギー の和を取り、この系の全エネルギーとする このような模型をイジング模型(Ising Model)と呼び、 磁性体の簡単なモデルになっている 以下、強磁性 (𝐽 > 0)の場合を考える 𝐸 = −2𝐽 + 2𝐽 = 0
25. 25 82 イジング模型 𝑖番目のスピンの状態を𝜎𝑖 とする 𝜎𝑖 の値は1(上向き)か-1(下向き)のいずれか −𝐽 𝐽 𝜎𝑖 = 1 𝜎𝑗 = 1 𝜎𝑖 = 1 𝜎𝑗 = −1 両方まとめて−𝐽𝜎𝑖 𝜎𝑗 と書ける
26. 26 82 イジング模型 全系のエネルギーは以下のように書ける 𝐻 = −𝐽 ෍ 𝑖,𝑗 𝜎𝑖 𝜎𝑗 系の全てのボンドについて和をとるという意味 全系のエネルギーを与える量をハミルトニアンとよぶ
27. 27 82 ボルツマン重み 系の状態に通し番号をつけ、𝑖番目の状態のエネルギーを𝐸𝑖 とする 状態𝑖 エネルギー 𝐸𝑖 = 4𝐽 𝑤𝑖 = exp(−𝛽𝐸𝑖 ) ボルツマン定数 𝛽 = 1/𝑘𝐵 𝑇 逆温度 𝑘𝐵 この状態の出現確率がボルツマン重みに比例する この時、エネルギーの温度依存性を知りたい
28. 28 82 エネルギーの振る舞い 低温 高温 スピンがそろった状態 出現確率は高いが、総数が少ない →エネルギー重視 →低温で支配的 スピンがバラバラの状態 出現確率は低いが、総数が多い →エントロピー重視 →高温で支配的
29. 29 82 エネルギーの振る舞い エネルギーの期待値の温度依存性 𝑈 𝛽 = 𝑍−1 ෍ 𝑖 𝐸𝑖 exp −𝛽𝐸𝑖 𝑍 = ෍ 𝑖 exp −𝛽𝐸𝑖 エネルギー𝐸をとる状態の数を𝑔(𝐸)とすると 𝑈 𝛽 = 𝑍−1 න 𝐸𝑔(𝐸) exp −𝛽𝐸 𝑑𝐸 𝑍 = න 𝑔(𝐸) exp −𝛽𝐸 𝑑𝐸 状態からエネルギー の計算は簡単 𝐸 = 0 あるエネルギーをとる状態の 数の計算は大変 𝐸 = 0
30. 30 82 N=4の場合の状態数 𝐸 = 0 𝐸 = 4J 𝐸 = −4J 𝑔 −4𝐽 = 2 𝑔 0 = 12 𝑔 4𝐽 = 2 これを一般のNで求めるのは極めて大変 モンテカルロサンプリング
31. 31 82 単純サンプリング 1. 全ての可能な状態から無作為に一つ選ぶ 2. その状態のエネルギー𝐸𝑖 と重み𝑤𝑖 を計算する 3. 1-2を繰り返し、重み付きで平均をとる スピンが𝑁個なら、状態は2𝑁個 ほとんどの状態は重みが小さいので サンプリング効率が悪い マルコフ連鎖モンテカルロ法 この中から無作為に一つ選ぶ
32. 32 82 マルコフ連鎖モンテカルロ法 現在の状態から遷移可能な状態を限定し、その中 から一つ無作為に選ぶ サンプリング候補の決め方 イジング模型への適用 1. 現在の状態を𝐴とし、スピンを一つ無作為に選ぶ 2. 選んだスピンを反転させた状態を提案状態𝐵とする 3. それぞれのエネルギーから遷移確率𝑃(𝐴 → 𝐵)を計算 し、遷移させるか決める 4. 遷移しなかった場合は状態はそのまま、遷移した場 合は提案状態を現在の状態にして1.へ
33. 33 82 マルコフ連鎖モンテカルロ法 現在の状態 遷移可能な状態 スピンが𝑁個なら総状態数は2𝑁、遷移可能な状態は𝑁 この中から無作為に一つ選ぶ 遷移可能な状態が限定された
34. 34 82 マルコフ連鎖モンテカルロ法 現在の状態𝐴 現在の状態𝐵 𝑃(𝐴 → 𝐵) 𝑃(𝐵 → 𝐴) 𝑃(𝐵 → 𝐵) 𝑃(𝐴 → 𝐴) 𝐸𝐴 = −4𝐽 𝐸𝐵 = 0 𝑤𝐴 = exp −𝛽𝐸𝐴 𝐸𝐴 < 𝐸𝐵 なので𝑤𝐴 > 𝑤𝐵 𝑤𝐵 = exp −𝛽𝐸𝐵
35. 35 82 マルコフ連鎖モンテカルロ法 𝑃(𝐴 → 𝐵) 𝑃(𝐵 → 𝐴) 𝑃(𝐵 → 𝐵) 𝑃(𝐴 → 𝐴) 確率の保存条件 𝑃 𝐴 → 𝐴 + 𝑃 𝐴 → 𝐵 = 1 𝑃 𝐵 → 𝐵 + 𝑃 𝐵 → 𝐴 = 1 詳細釣り合い条件 𝑃 𝐴 → 𝐵 𝑃 𝐵 → 𝐴 = 𝑤𝐵 𝑤𝐴
36. 36 82 マルコフ連鎖モンテカルロ法 メトロポリス法(Metropolis method) 重みが大きくなる場合は必ず遷移、そうでない場合は確率的に遷移させる 𝑃 𝐵 → 𝐴 = 1 𝑃 𝐴 → 𝐵 = 𝑤𝐵 𝑤𝐴 = exp −𝛽𝐸𝐵 exp(−𝛽𝐸𝐴 ) = exp(−𝛽Δ𝐸) Δ𝐸 = 𝐸𝐵 − 𝐸𝐴 > 0なので、 exp −𝛽Δ𝐸 < 1 提案状態のエネルギーが、現在の状態より低ければ必ず遷移 高ければ確率𝑝 = exp −𝛽Δ𝐸 < 1で遷移
37. 37 82 マルコフ連鎖モンテカルロ法 イジング模型にマルコフ連鎖モンテカルロ法を適用した場 合のアルゴリズム(メトロポリス法を採用した場合) 1. スピンを一つ無作為に選ぶ 2. 選んだスピンを反転させた状態を提案状態とする 3. 現在の状態とのエネルギー差Δ𝐸を計算する 4. エネルギーが下がる場合(Δ𝐸 < 0)なら必ず遷移。そ うでなければ確率𝑝 = exp(−𝛽Δ𝐸)で遷移 5. 1-4を繰り返す 一度に一つのスピンだけひっくり返すので Single-Spin-Flip Algorithmと呼ぶ
38. 38 82 ここまでのまとめ 状態に重みがある場合の期待値を厳密に計算することは難しい その理由は同じ重みを持つ状態の数がわからないから →サンプリングによる期待値の推定 単純サンプリング 全ての状態から無作為に状態をサンプリングする 重みの小さい状態ばかりサンプリングするため非効率 マルコフ連鎖モンテカルロ法 現在の状態をもとに、サンプリングする状態を限定する 現在の状態と提案状態の重みを比較し、適切な遷移確率を 作ることで重みに比例したサンプリングができる (重み付きサンプリング)→効率的
39. 39 82 モンテカルロ法の応用 アニーリング • 数値計算では、系の基底状態に興味がある ことが多い。しかし、基底状態を探すこと は難しい→アニーリングによる探索 レプリカ交換モンテカルロ法 • 複数の安定な状態がある場合、どちらかし かサンプリングできなくなることがある。 しかし、両方を効率よくサンプリングした い→レプリカ交換モンテカルロ法
40. 40 82 基底状態探索 タンパク質は、アミノ酸がつながってできている それが水の中で折りたたまれて、なんらかの形を持つ この形は、水の中でエネルギーが低い状態であると考えられる どのようにその構造を予測すればよいだろうか？
41. 41 82 基底状態探索 一次元タンパク質を二次元格子上で折りたたませる ただし、同じ色はなるべく上下左右に隣り合うようにしたい 最終的にこんな形状を探したい
42. 42 82 基底状態探索 状態 Ԧ 𝑥に対して、自由エネルギー𝐹( Ԧ 𝑥)が定義される時、 𝐹( Ԧ 𝑥)を最小にする Ԧ 𝑥を探したい ※系の状態空間は一般に超高次元だが、以後は一次元的に表現する 𝑥 𝐹(𝑥) この場所を知りたい この関数を自由エネルギー ランドスケープと呼ぶ
43. 43 82 基底状態探索 自由エネルギーランドスケープがすり鉢型の場合 𝑥 𝐹(𝑥) 現在地 目的地 𝑥𝑡 𝑥𝑡+1 = 𝑥𝑡 + Δ𝑥もしくは𝑥𝑡+1 = 𝑥𝑡 − Δ𝑥として 𝐹(𝑥𝑡+1 )が低い方へ進んでいけば良い→(確率的)最急降下法
44. 44 82 基底状態探索 𝑥 𝐹(𝑥) 目的地 現在地 必ずエネルギーが低い方向に進ませると、局所的にエネルギーが低い 状態(ローカルミニマム)にトラップされ、正しい基底状態(グローバル ミニマム)に到達できない→エネルギーが高い方向にも進ませる
45. 45 82 基底状態探索 ෤ 𝑥𝑡+1 = 𝑥𝑡 + Δ𝑥 2 Ƹ 𝑟 − 1 Ƹ 𝑟: 0から1までの一様乱数 𝐸(෤ 𝑥𝑡+1 ) < 𝐸(𝑥𝑡 ) まず、現在の状態𝑥𝑡 から提案状態෤ 𝑥𝑡+1 を作る エネルギーが下がったら採用 𝑥𝑡+1 ← ෤ 𝑥𝑡+1 エネルギーが上がっても、確率的に採用 Δ𝐸 ≡ 𝐸(෤ 𝑥𝑡+1 ) − 𝐸(𝑥𝑡 ) ቊ 𝑥𝑡+1 ← ෤ 𝑥𝑡+1 Ƹ 𝑟 < exp 𝛽Δ𝐸 𝑥𝑡+1 ← 𝑥𝑡+1 otherwise
46. 46 82 基底状態探索 𝑥 𝐹(𝑥) 目的地 現在地 提案状態 低確率だが、エネルギーが高い方向にも遷移する 長時間極限において必ず一番エネルギーが低い状態に到達する
47. 47 82 ローカルミニマムの突破 長時間極限において必ず一番エネルギーが低い状態に到達する どのくらい待てばよいか？ Δ𝐸 𝜏 ∝ exp 𝛽Δ𝐸 エネルギーバリア エネルギーバリアに対して指数関数的に脱出時間が長くなる
48. 48 82 温度の効果 𝜏 ∝ exp 𝛽Δ𝐸 𝛽 ≡ 1/𝑘𝑇 脱出時間： 逆温度： 温度が低い(𝛽が大きい)ほど待ち時間が長い 温度が高い(𝛽が小さい)ほど待ち時間が短い 温度をめちゃくちゃ高くすれば、あっという間に 基底状態を見つけられるのでは？
49. 49 82 温度の効果 温度が低すぎる 温度が高すぎる ローカルミニマムにひっかかりまくる 細かい構造が全て消えてしまい、正確 な情報が得られない
50. 50 82 アニーリング 高温 中温 低温 まず高温である程度 の目処をつける 徐々に温度を下げて いくと、構造が見え てくる 最終的に温度ゼロに 持っていくことで、 グローバルミニマム を見つける 高温から始めて、徐々に温度を下げていく手法をアニーリングと呼ぶ 十分遅く温度を下げれば、必ずグローバルミニマムに到達 アニーリングの温度の下げ方には経験と工夫が必要
51. 51 82 レプリカ交換モンテカルロ法 アニーリングは温度の下げ方に対して経験が必要 長時間計算するのに最終的に基底状態しか得られない 温度の下げ方を気にしたくない 途中の温度の状態もちゃんと調べたい レプリカ交換モンテカルロ法 (Replica Exchange Monte Carlo Method, REMC) R. H. Swendsen and J-S. Wang, Phys. Rev. Lett. 57 2607 (1986) K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65 1604 (1996)
52. 52 82 レプリカ交換モンテカルロ法 全く同じ系を複数用意する(レプリカ) 異なる温度で計算を実行。たまに温度を交換する。 高温 低温 𝛽1 𝛽2 𝛽3 𝛽4 ・・・
53. 53 82 レプリカ交換モンテカルロ法 𝐸𝑖 , 𝐸𝑗 ：レプリカiとjのエネルギー 𝛽𝑖 , 𝛽𝑗 ：レプリカiとjの逆温度 上記２つで温度交換をする確率𝑝を以下のように定める 𝑝 = min 1, exp 𝐸𝑖 − 𝐸𝑗 (𝛽𝑖 − 𝛽𝑗 ) 温度が高い方がエネルギーが低く、温度が低い方がエネルギーが高い場合は必ず交換 そうでない場合はたまに交換 • 十分緩和したら、それぞれの温度の平衡状態が一度に得られる • 系の温度が勝手に上下するため、自動的にアニーリングの効果が得られる • 並列計算と相性が良い • 分子動力学法にも応用可能 • 温度間隔の設定には経験が必要
54. 54 82 モンテカルロ法の応用のまとめ • 基底状態を調べたい時、最急降下法やモンテカルロ法ではローカル ミニマムにトラップされてしまい、効果的にサンプリングできない • 温度が高いほどローカルミニマムの脱出時間が短くなるが、その分 細かい構造が見えなくなる 問題点 アニーリング • 高温の状態から徐々に温度を下げていくことで、ローカルミニマム にはトラップされず、最終的に基底状態を得ることができる • 温度の下げ方には経験が必要 レプリカ交換モンテカルロ法 • 様々な温度の「レプリカ」を作り、同時にシミュレーションする。 たまに温度を交換することで、自動的にアニーリングする • 並列計算と相性が良いが、計算コストは高め
55. 55 82 詳細釣り合い条件の応用 時間の矢の起源 貧富の差の起源 今日は話しません 興味ある人は 「リンゴゲームと貧富の差」で検索
56. 56 82 時間の矢とは？ 我々の感じる「時間」は一方向に流れている 過去の記憶はあるが、未来の記憶はない これはなぜだろう？
57. 57 82 時間反転対称性 ミクロな支配方程式は時間反転対称性を持つ 𝑚 𝑑2𝑥 𝑑𝑡2 = 𝐹(𝑥) 例：ニュートンの運動方程式 運動を録画したビデオを逆再生してもどちらが 正方向か区別がつかない 時間に関して二階微分 →もし𝑥(𝑡)が解なら、𝑥(−𝑡)も解
58. 58 82 時間反転対称性 マクロな支配方程式は時間反転非対称性 𝜕𝜌 𝜕𝑡 = 𝐷 𝜕2𝜌 𝜕𝑥2 例：拡散方程式 時間に関して一階微分 →もし𝜌(𝑥, 𝑡)が解でも、𝜌(𝑥, −𝑡)は解にならない 拡散現象を録画したビデオを逆再生したら逆再 生とわかる
59. 59 82 ミクロからマクロへ 水にインクを垂らすと拡散していく 水原子の動き 𝑚 𝑑2𝑥 𝑑𝑡2 = 𝐹(𝑥) 𝜕𝜌 𝜕𝑡 = 𝐷 𝜕2𝜌 𝜕𝑥2 マクロには時間反転非対称 ミクロには時間反転対称
60. 60 82 エーレンフェストの壺 1 2 3 4 5 6 • 2つ壺を用意する • 数字が書かれた玉をN個用意する • 一つの玉をランダムに選んで、その玉をもう一方に移す • 最初は片方にすべての玉を入れておく • 片方の壺の玉の数の時間発展を追う
61. 61 82 エーレンフェストの壺 N=10の場合 ステップ 壺 の 中 の 玉 の 数
62. 62 82 エーレンフェストの壺 N=100の場合 ステップ 壺 の 中 の 玉 の 数
63. 63 82 エーレンフェストの壺 N=1000の場合 ステップ 壺 の 中 の 玉 の 数 どんな状態からスタートしても 玉が半分ずつの状態に収束する
64. 64 82 エーレンフェストの壺 1 2 3 4 5 6 ミクロな操作は可逆 マクロな観測事実は不可逆 時間の矢 どこで時間反転対称性が破れたのか？ ※逆過程が等確率で起きる ※初期条件を忘れる
65. 65 82 玉が1個の場合 玉が一つの場合、右側の壺の状態は二通り 1 玉がない 玉がある
66. 66 82 玉が1個の場合 状態遷移図 1 1 状態がくるくる回ってしまって定常状態にならない (偶数回と奇数回でそれぞれ必ず異なる状態になる) 1
67. 67 82 エーレンフェストの壺(改) 1 2 3 4 5 6 • 2つ壺を用意する • 数字が書かれた玉をN個用意する • 一つの玉をランダムに選んでその玉をもう一方に移すが、 確率ε(0< ε <1)でなにもしない • 最初は片方にすべての玉を入れておく • 片方の壺の玉の数の時間発展を追う
68. 68 82 玉が1個の場合 状態遷移図 1 − 𝜀 同じ状態にとどまる可能性があるため、定常状態が存在する 1 − 𝜀 𝜀 𝜀 1
69. 69 82 確率の時間発展と定常状態 1 tステップ目に玉がない確率 𝑝𝜙 𝑡 𝑝1 𝑡 tステップ目に玉がある確率
70. 70 82 確率の時間発展と定常状態 𝑝𝜙 𝑡+1= 𝜀𝑝𝜙 𝑡 + (1 − 𝜀)𝑝1 𝑡 𝑝1 𝑡+1= (1 − 𝜀)𝑝𝜙 𝑡 + 𝜀𝑝1 𝑡 確率の時間発展 Ԧ 𝑝𝑡 = 𝑝𝜙 𝑡 𝑝1 𝑡 と書くと Ԧ 𝑝𝑡+1 = 𝑀 Ԧ 𝑝𝑡 𝑀 = 𝜀 1 − 𝜀 1 − 𝜀 𝜀 ただし
71. 71 82 確率の時間発展と定常状態 𝑀 = 𝜀 1 − 𝜀 1 − 𝜀 𝜀 マルコフ行列の最大固有値は1 最大固有値に対応する固有ベクトルが定常状態 Ԧ 𝑝∞ = 𝑀 Ԧ 𝑝∞ もし Ԧ 𝑝∞ が定常状態なら、 𝑀をかけても状態がかわらない に対応する固有ベクトルは Ԧ 𝑝∞ = 1/2 1/2 定常状態は2つの状態が等確率で現れる 1 =
72. 72 82 玉が2個の場合 1 tステップ目に玉がない確率 𝑝𝜙 𝑡 𝑝1 𝑡 tステップ目に玉1がある確率 2 𝑝2 𝑡 tステップ目に玉2がある確率 𝑝12 𝑡 tステップ目に玉1,2がある確率 2 1
73. 73 82 状態遷移図(N=2) 1 2 2 1 𝜀 1 2 (1 − 𝜀) 1 2 (1 − 𝜀) 1 2 (1 − 𝜀) 1 2 (1 − 𝜀) 1 2 (1 − 𝜀) 1 2 (1 − 𝜀) 1 2 (1 − 𝜀) 1 2 (1 − 𝜀) 𝜀 𝜀 𝜀
74. 74 82 ミクロな対称性 1 確率1/2で玉1が選ばれ、かつ確率(1 − 𝜀)で玉を移す 1 2 (1 − 𝜀) 1 1 2 (1 − 𝜀) すべての2つの状態間の遷移確率は等しい 全ての過程と逆過程は等確率で起きる →可逆過程
75. 75 82 遷移行列 𝑀 = 𝜀 𝑐 𝑐 0 𝑐 𝜀 0 𝑐 𝑐 0 𝜀 𝑐 0 𝑐 𝑐 𝜀 𝑐 ≡ 1 2 1 − 𝜀 Ԧ 𝑝∞ = 1/4 1/4 1/4 1/4 最大固有値に対応する固有ベクトル 十分時間が経つと、全てのミクロな状態は等確率で出現する →等重率の原理 1 2 2 1 = = =
76. 76 82 粗視化 1 2 = = 玉の数字を見ないことにする 同一視
77. 77 82 粗視化 1 2 = = 玉の数字を見ないことにする 同一視
78. 78 82 粗視化 十分時間がたった後に片方の壺を観察すると 𝑝0 = 1 4 𝑝1 = 1 2 𝑝2 = 1 4 玉がない 玉が1個 玉が2個 1 2 玉が1つ(=N/2)ある状態を観測する確率が最も高くなった
79. 79 82 粗視化：玉がN個の場合 十分時間がたった後に片方の壺を観察すると 玉がない 玉が1個 𝑝0 = 1 2𝑁 𝑝1 = 𝑁 2𝑁 … 玉がn個 𝑝𝑛 = 𝐶𝑛 2𝑁 𝑁 N個の玉がある→ミクロな状態は2𝑁個 n個の玉がある状態→ 𝐶𝑛 個 𝑁
80. 80 82 粗視化：玉がN個の場合 𝑁 = 1000 Nが大きい時にN/2を中心とするガウス分布に収束 玉の数がほぼN/2である状態が観測される 𝑛 𝑝𝑛
81. 81 82 時間の矢のまとめ • エーレンフェストの壺はミクロには可逆、マクロには不可逆 • ミクロとは「全ての玉の番号を知っている状態」 • マクロとは「玉の区別をなくした状態」 • ミクロにはすべての状態が等確率で出現する →等重率の原理 • マクロには玉が半分ずつに分かれる状態に収束する →時間の矢 • 古典的には「粗視化」が時間反転対称性を破る https://www.gakushuin.ac.jp/~881791/materials/Irreversiblity09.pdf 参考資料：田崎晴明 なぜ時間に向きがあるのだろう？ 現実のこの世界は？
82. 82 82 まとめ • 履歴によらず次の状態が決まることをマルコフ性 と呼ぶ • 状態が離散的であり、マルコフ的に次の状態が決 まる連鎖をマルコフ連鎖と呼ぶ • 状態と遷移確率を描いた図を状態遷移図と呼ぶ • 遷移確率の行列表現をマルコフ行列と呼ぶ • 任意の二状態間で、定常分布の確率と遷移確率の 間に成り立つ条件を詳細釣り合い条件と呼ぶ • 遷移確率を決めると定常分布が決まるし、逆に定 常分布から遷移確率を決めることもできる • 遷移確率の決め方は一意に決まらない • 熱浴法：分布の比で決める • メトロポリス法：片方の遷移確率を1にする