レプリカ交換モンテカルロ法を用いた力学系の軌道・パラメータ探索

MCMC 力学系

最近,マルコフ連鎖モンテカルロ法の発展形である拡張アンサンブル法について勉強しています．ただ本を読むだけではなく,実際に手を動かして生きた知識としていくことが重要と考えています.そこで,まずは柳田先生と伊庭先生のレプリカ交換モンテカルロ法を用いた力学系の軌道・パラメータ探索に関する次の論文を追試することから始めたいと思います．

T.Yanagita, and Y.Iba, Exploration of order in chaos using the replica exchange Monte Carlo method, J.Stat.Mech., P02043, 2009.

この論文を理解するためには，レプリカ交換モンテカルロ法を含むマルコフ連鎖モンテカルロ法について，また力学系の基本的な知識が必要になります．それぞれについて順次簡単に説明し（7/5現在,概略しか書いていないので,余裕が出来次第,加筆します）,最後に柳田先生と伊庭先生の論文を追試します．なお,計算量も多いため本来はC++等を用いて実装すべきだと思われますが,効率を優先してRのパッケージを多用してプログラムを書きます（今後,C++への移植,またRでの並列計算の実装をしたいと思っています・・・）．

マルコフ連鎖モンテカルロ法

マルコフ連鎖モンテカルロ法(以下,MCMC)は,与えた確率分布から効率的にサンプリングを行うための一連の手法を指します．MCMCは高次元の確率分布に対しても有効で,非常に強力です.

MCMCの代表的な例として,メトロポリス法,メトロポリス・ヘイスティング法,ギブス・サンプラーなどが知られています.これらの手法の入門的な内容は@teramonagi氏のスライドが分かりやすくてお勧めです.

しかし,MCMCにも問題点があります．伊庭先生の計算統計IIでは,次の2つの問題点が指摘されています．

  1. 多峰性の分布の場合など,混合に時間がかかるケースでは計算に時間がかかる（あるいは事実上不可能な）こと
  2. 分布の正規化定数に相当する多重和や多重積分が直接計算できないこと

このうち,特に前者の問題点を解決するために提案されたのが拡張アンサンブルモンテカルロ法(extended ensemble Monte Carlo, generalized ensemble Monte Carlo)と呼ばれる一連の手法です．これらの方法では,パラメータの異なるいろいろな分布をまとめた分布を計算の対象として,効率よく混合することで問題点を回避します．今回,フォーカスを当てるレプリカ交換モンテカルロ法は,拡張アンサンブルモンテカルロ法の代表的な方法の一つです.

レプリカ交換モンテカルロ法

レプリカ交換モンテカルロ法(パラレル・テンパリング法,Metropolis-coupled MCMC)では,異なるパラメータの値 \theta_{k}, \, k=1, \, \dots, \, K を持つ分布 P(x|\theta_{k}) を複数まとめた同時分布
P(x_{1}, x_{2}, \dots, x_{K}) = \prod_{k=1}^{K} \, P(x_{k} \mid \theta_{k})
をマルコフ連鎖モンテカルロ法でサンプリングします．

上記の分布を定常分布にするような遷移として次の2種類を考えて,これらを交互に実行します．

  1. 各分布に対するMCMCの実行　各分布 P(x_{k} \mid \theta_{k}) について, P(\cdot \mid \theta_{k}) を不変にする x_{k} の遷移を考え,メトロポリス法やギブス・サンプラーを用いてサンプリングする.
  2. レプリカの確率的交換　適当なステップ数ごとに,ランダムに選んだ 1 \leq k < K-1 について状態 x_{k} と x_{k+1} を確率 \min(1, r) で交換する．

ここで,
r = \frac{P(x_{1}, \dots, x_{k+1}, x_{k}, \dots, x_{K})}{P(x_{1}, \dots, x_{k}, x_{k+1}, \dots, x_{K})} = \frac{P(x_{k+1} \mid \theta_{k}) P(x_{k} \mid \theta_{k+1})}{P(x_{k} \mid \theta_{k}) P(x_{k+1} \mid \theta_{k+1})}
です．

力学系とカオス

力学系は,初期値が与えられるとダイナミクスを支配する規則に従って時間発展するシステムを指します．典型的な例として,常微分方程式,写像(差分方程式)が挙げられます．

写像の例として,次式で表される1次元写像であるロジスティック写像の軌道を計算してみましょう．
x_{t+1}=ax_{t}(1-x_{t})
ここで,a はパラメータです．

パラメータa を変化させて,初期値からスタートして最終的に落ち着く軌道をプロットすると,リンク先にあるような図が得られます.この図は「分岐図」と呼ばれます．パラメータaが3.6の手前までは周期解が得られていますが,それよりもaが大きいと複雑な軌道となっています(aが3.8より少し大きい値などで周期解が得られていますが,こうしたパラメータaの領域を「窓」と呼びます).この複雑な軌道がいわゆるカオス軌道です.

常微分方程式の例として，次式で表されるLorenz方程式の軌道を計算してみましょう．
\left\{\begin{array}{l} \frac{\mathrm{d}x}{\mathrm{d}t}= -\sigma(x - y) \\ \frac{\mathrm{d}y}{\mathrm{d}t}= (r - z)x -y \\ \frac{\mathrm{d}z}{\mathrm{d}t}= xy -bz \end{split}\right.
ここで,\sigma, r, bはそれぞれパラメータです．
Rの常微分方程式を解くodesolveパッケージ，OpenGLを使用して3次元のプロットを行うrglパッケージを用いて軌道を計算し表示します．plot3d関数でプロットしたLorenzアトラクタは，回転させて眺めることができます．

> install.packages(odesolve); library(odesolve)
> install.packages(rgl); library(rgl)
> # Lorenz方程式のパラメータ
> params <- list(s=10, b=8/3, r=28)
> # 軌道計算を行う時刻
> times <- seq(0, 1000,  by=0.01)
> # 軌道計算
> y <- lsoda(c(1,0,0), times, lorenz, params)
> plot3d(y[, 2:4], type="l", color="red", xlab="x", ylab="y", zlab="z")

なお,rglパッケージをインストールするためには，あらかじめOpenGLがインストールされている必要があります．例えば，ubuntuの場合は次のようなコマンドを打ってインストールを行います．

$ sudo apt-get install freeglut3 freeglut3-dev freeglut3-dbg

力学系の不安定周期軌道

分岐図では安定な軌道,すなわち長時間経過後に与えた初期値から出発する軌道が吸引される軌道がプロットされます．

しかし，実は分岐図で描かれる安定的な軌道以外にも，状態空間の中には不安定な軌道が存在しています．この不安定周期軌道は，カオスアトラクタの骨格を担っています．ところが，不安定周期軌道は一般的に吸引領域(時間発展させた結果，不安定周期軌道に引き込まれていく初期値の集合のこと)が非常に狭く，検出しづらいことが多いです(以下，省略．余裕ができたら加筆します)．

以下で説明する柳田先生と伊庭先生の論文では，レプリカ交換モンテカルロ法を用いて吸引領域が非常に狭い不安定周期軌道や系のパラメータなどの「まれな」現象を捉えます．

力学系の軌道探索へのレプリカ交換モンテカルロ法の適用

さて,以上で準備は終わりです．ここからは柳田先生と伊庭先生の論文の内容紹介と追試を行います．なお，論文では適用例として，Lorenz方程式と二重振り子の2つの例が載っていますが,ここでは前者のみを扱います.

提案された手法の概要
本論文では,力学系の「まれな」初期条件や系のパラメータを求めるために，レプリカ交換モンテカルロ法を用います.

軌道が周期的であるほど値が小さくなるエネルギー関数 E(X) を定義します．
エネルギー関数が定義できれば，以下のようにこのエネルギーに関してのギブス分布が定義されます.
p(X) = \frac{\exp{(-\beta E(X))}}{Z(\beta)}, \, \, \, Z=\int \exp{(-\beta E(X)) \, \mathrm{d}X.

p(X) は, \beta が小さくすると一様分布に近づき， \beta が大きくするとある領域に集中します．すなわち， \beta の値を制御することにより，非常に不安定な周期軌道等の望んだ性質を持つ軌道や系のパラメータのサンプリングを行えるようになります．本論文では,様々な \betaの値からサンプリングを行うとともに,局所領域に陥ることを回避するためにレプリカ交換を行います.

アルゴリズム
アルゴリズムは以下のとおりです。 xが初期時刻での状態空間上での位置，Tが周期，\alphaが系のパラメータを表します．なお,モンテカルロシミュレーションのステップ数を K, レプリカ数を N とします.

k = 1, \, \dots, K に対して以下を実行する.

レプリカ i=1, \, \dots, \, N に対して,
(i) 現在の状態に摂動 (\Delta x_{i}, \Delta T_{i}, \Delta \alpha_{i}) を加える．\Delta x_{i}は,s\times d \times 10^{-e}の形で与える．ここで，d\in \{1, \, \dots, \, 9\}, s\in\{-1, 1\}, e\in\{e_{\min}, e_{\max}\}であり，各区間上の一様分布から1つずつサンプリングする.
現在の値が (x^{(n)}_{i}, T^{(n)}_{i}, \alpha^{(n)}_{i}) ならば，次の状態の候補は (x'_{i}, T'_{i}, \alpha'_{i}) = (x^{(n)}_{i} + \Delta x_{i}, T^{(n)}_{i} + \Delta T_{i}, \alpha^{(n)}_{i} + \Delta \alpha_{i}) となる．
(ii) 初期状態 x'_{i} からスタートして，パラメータ \alpha'_{i} の力学系の長さ T'_{i} の軌道を計算する．その結果から， E(\alpha'_{i}, T'_{i}, x'_{i}) を計算する．
(iii) 一様乱数 R \in [0, 1] を生成し，
R < \frac{\exp{(-\beta E(x'_{i}, T'_{i}, \alpha'_{i}))}}{\exp{(-\beta E(x^{(n)}_{i}, T^{(n)}_{i}, \alpha^{(n)}_{i}))}}
となる場合のみ提案を受理して，
\left\{\begin{array}{c} x^{(n+1)}_{i} = x'_{i} \\ T^{(n+1)}_{i} = T'_{i} \\ \alpha^{(n+1)}_{i} = \alpha' \end{array} \right.
とする．それ以外の場合は，
\left\{\begin{array}{c} x^{(n+1)}_{i} = x^{(n)}_{i} \\ T^{(n+1)}_{i} = T^{(n)}_{i} \\ \alpha^{(n+1)}_{i} = \alpha^{(n)}_{i} \end{array} \right.
とする．

kがあらかじめ決めたレプリカ交換を行うステップ数である場合は，添字の組の集合\{(1,2), (2,3), \dots, (N-1, N)\}からそれぞれ同一の確率で1組(m,m+1)を選択し， x_{m}と x_{m+1}を確率\min\left(1, (\beta_{m+1} - \beta_{m})(E_{m+1} - E_{m})\right)で交換する．

Lorenz方程式に対しては,エネルギー関数を次式で定義します．
E(\vec{x}_{0}, T) = \log{(|\phi^{T}(\vec{x}_{0}) - \vec{x}_{0}| + P(\vec{x}_{0}, T) + \epsilon),
P(\vec{x}_{0}< T) = g\left(\frac{1}{T}\int_{0}^{T} \, \phi^{\tau}(\vec{x}_{0})^2 \, \mathrm{d}\tau, v_{c}\right) + g(T, T_{c})
ここで関数 g=g(s,s_{c}) は，
g=g(s, s_{c})=\Theta(s_{c} - s)\left(\frac{1}{s} - \frac{1}{s_{c}}\right)
で定義される． \Theta(s) はHeaviside関数であり,
\Theta(s) = \left\{\begin{array}{cc} 1 & (s\geq 0) \\ 0 & (s < 0)\end{array}\right.
です．

初期条件のサンプリング
方程式の対称性 (x, y) \leftrightarrow (-x, -y) を考慮して,論文ではPoincare断面 \vec{x}_{0} = (x_{0}, y_{0}, r-1) から初期値をサンプリングしています．ここで, x_{0} \in [0, 20], y_{0} \in [0, 20]です．

軌道の周期 T は,区間 [0, T_{\max}]に存在すると仮定しています． T_{\max}=6.0 とします．

実装

コードは次のようになります．少し長くなるので,コンソール上に打ち込むのではなくファイルに保存してコンソールから実行します．

rpm_with_ode.r

replica.exchange.montecarlo.with.ode <-
  function(step.replica, step.exchange, step.trans, n.replica, func.energy,
            inv.temp,
            e.min, e.max,
            ode, params,
            period.min, period.max,
            init.state, t.start, t.delta,
            exec.which,
            outputdir, fn, result.inv.temp, pch, col
  )
{
  # レプリカ交換モンテカルロ法を用いた力学系の軌道探索

  # step.replica: レプリカ交換モンテカルロ法の全ステップ数
  # step.exchange: レプリカ交換の時間間隔
  # step.trans: 過渡ステップ数
  # n.replica: レプリカ数
  # func.energy: エントロピー関数
  # inv.temp: 逆温度
  # e.min: eの最小値
  # e.max: eの最大値
  # ode: 常微分方程式
  # params: 常微分方程式のパラメータ
  # period.min: 最小周期
  # period.max: 最大周期
  # init.state: 初期時刻における状態(位置,周期)
  # t.start: 常微分方程式の軌道を計算する際の初期時刻
  # t.delta: 常微分方程式の軌道を計算する際の時刻の刻み幅
  # exec.which: レプリカモンテカルロの対象
  # (順に軌道の初期位置, 周期, 常微分方程式のパラメータ)
  # outputdir: 出力ディレクトリ
  # fn: 出力ファイル名
  # result.inv.temp: 結果表示に用いる逆温度

  pair.all <- rbind(1:(n.replica-1), 2:n.replica) #隣り合うペアのインデクス

  r.adoption <- runif(step.replica * n.replica) # 採択判定用一様乱数
  r.exchange <- runif(floor(step.replica/step.exchange)) # レプリカ交換判定用一>様乱数
  energy <- rep(NA, n.replica) # あるステップにおける各レプリカのエネルギー

  # 結果の格納配列
  pos <- match(result.inv.temp, inv.temp)
  vars <- c("T", "E")
  result <- array(NA, dim=c(step.replica-step.trans, length(vars), length(pos)),
                  dimnames=list((step.trans+1):step.replica, vars,
                                result.inv.temp))

cat("モンテカルロシミュレーション\n")
  n <- 1 # レプリカ番号
  sapply(1:step.replica,
          function(k) {
cat("k=", k, "\n")
            N <- length(init.state) * n.replica # パラメータ数×レプリカ数
            perturb <- sample(c(-1, 1), N, replace=T) *
                        sample(1:9, N, replace=T) *
                        10^(-sample(e.min:e.max, N, replace=T)) # 摂動
            perturb <- matrix(perturb, ncol=ncol(init.state)) # 摂動行列

            y <- apply(init.state, 1,
                      function(x) {
                        ans <- metropolis.with.ode(
                                  x, ode, params,
                                  period.min, period.max,
                                  t.start, t.delta,
                                  perturb[n, ], inv.temp[n],
                                  r.adoption[n], n)
                        energy[n] <<- ans$energy
                        n <<- n+1
                        return(ans$x)
                      }
                  )
            init.state <<- t(y)
            # レプリカ交換
            if (k %% step.exchange == 0) {
cat("レプリカ交換\n")
              idx <- sample(1:ncol(pair.all), 1)
              pair.idx <- pair.all[, idx]
              r <- exp(diff(inv.temp[pair.idx]) * diff(energy[pair.idx]))
              if (r > r.exchange[as.integer(k/step.exchange)]) {
                init.state[pair.idx, ] <<- init.state[rev(pair.idx), ]
                energy[pair.idx] <<- energy[rev(pair.idx)]
              }
            }
            if (k > step.trans) {
              result[k - step.trans, ,] <<-
                t(cbind(init.state[pos, "period"], energy[pos]))
            }
            n <<- 1
          }
  )
  # 結果を配列から行列に変換
  result <- do.call(rbind, result)

  jpeg(paste(outputdir, fn, sep=""))
  plot(result, pch=pch, col=col)
  dev.off()
}

metropolis.with.ode <-
  function(x, ode, params, period.min, period.max, t.start, t.delta,
            perturb, inv.temp, r, n)
{
  # メトロポリス法

  # 現在の初期条件のエネルギー算出
  period <- x["period"]
  if (period > period.max) x["period"] <- period <- period.max
  if (period < period.min) x["period"] <- period <- period.min
  xx <- solve.ode(x, period, ode, params, t.start, t.delta)
  energy.present <- energy.lorenz(xx, period)

  # 候補初期条件のエネルギー算出
  x.cand <- x + perturb
  period.cand <- x.cand["period"]
  if (period.cand > period.max) x.cand["period"] <- period.cand <- period.max
  if (period.cand < period.min) x.cand["period"] <- period.cand <- period.min
  xx.cand <- solve.ode(x.cand, period.cand, ode, params, t.start, t.delta)
  energy.cand <- energy.lorenz(xx.cand, period)

  # 採択の判断
  z <- exp(-inv.temp * (energy.cand - energy.present))
  if (z > r) return(list(x=x.cand, energy=energy.cand))
    else return(list(x=x, energy=energy.present))
}

solve.ode <- function(x, period, ode, params, t.start, t.delta)
{
  # x: 初期値
  # ode: 常微分方程式
  # params: 常微分方程式のパラメータ
  # t.start: 軌道計算の初期時刻
  # t.delta: 軌道計算の時刻刻み幅

  # 軌道計算
  library(odesolve)
  xx <- x[1:3]
  period <- x[4]

  if (t.start <= period) {
    times <- seq(t.start, period, by=t.delta)
    y <- lsoda(xx, times, ode, params)
    return(y)
  } else {
    return(matrix(x, nrow=1, byrow=T))
  }
}

lorenz <- function(t, x, params)
{
  # Lorenz方程式
  s <- params[["s"]]
  r <- params[["r"]]
  b <- params[["b"]]
  list(c(s * (x[2] - x[1]),
      x[1] * (r - x[3]) - x[2],
      x[1] * x[2] - b * x[3]))
}

energy.lorenz <- function(y, period, eps=0.001, Vc=4, Tc=4)
{
  # Lorenz方程式のエネルギー関数
  # y: 計算された軌道
  # period: 軌道の周期
  # eps:
  # Vc: 速度の閾値
  # Tc: 閾値

  if (identical(period, 0.0)) {
    return(log(eps))
  } else { 
    y0 <- y[1, 2:4]
    yT <- y[nrow(y), 2:4]

    period.measure <- log(sqrt(sum((yT - y0)^2))) # 周期性を表す項
    penalty.velocity <- auxiliary.func(sum(y[, 2:4]^2)/period, Vc)
    penalty.state <- auxiliary.func(period, Tc)
    penalty <- penalty.velocity + penalty.state

    return(log(sqrt(sum((yT - y0)^2)) + penalty + eps))
  }
}

auxiliary.func <- function(s, sc)
{
  return(heaviside(sc - s) * (1/s - 1/sc))
}

heaviside <- function(s)
{
  # Heaviside関数
  ifelse(s>=0, 1, 0)
}

最後に,実行関数です．Rだと計算速度が遅いため,モンテカルロシミュレーションのステップ数K=6000, レプリカ交換は50ステップに1回,最初の4000ステップは過渡状態として捨てることにします.

exec.r

exec <- function()
{
  # パラメータ設定
  step.replica <- 6000
  step.exchange <- 50
  step.trans <- 4000
  n.replica <- 31
  func.energy <- energy.lorenz
  e.min <- 1
  e.max <- 5
  ode <- lorenz
  params <- list(s=10, b=8/3, r=28)

  inv.temp <- seq(2, by=0.1, length.out=n.replica) # 逆温度
  period.min <- 0 # 周期の最小値
  period.max <- 6 # 周期の最大値
  t.start <- 0 #
  x0 <- runif(n.replica, 0, 20) # 初期位置(x)
  y0 <- runif(n.replica, 0, 20) # 初期位置(y)
  z0 <- params[["r"]] - 1 # 初期位置(z)
  period <- runif(n.replica, period.min, period.max) # 周期
  init.state <- cbind(x0, y0, z0, period) # 初期位置(x, y, z, T)
  t.delta <- 0.01
  which.exec <- c(rep(T, 2), F)
  outputdir <- "../output/"
  fn <- "test.jpg"
  result.inv.temp <- seq(2, 5, by=1)
  each <- rep(step.replica-step.trans, length(result.inv.temp))
  pch <- rep(c(4, 2, 20, 1), each=each)
  col <- rep(c("blue4", "magenta4", "chocolate4", "darkgreen"), each=each)

  # レプリカモンテカルロ法実行
  replica.exchange.montecarlo.with.ode(
    step.replica, step.exchange, step.trans, n.replica, func.energy,
    inv.temp,
    e.min, e.max,
    ode, params,
    period.min, period.max,
    init.state, t.start, t.delta,
    exec.which,
    outputdir, fn, result.inv.temp, pch, col
  )
}

Rのコンソール上で以上のように実行します．

> source("rpm_with_ode.r"); source("exec.r"); exec()

次の図が得られます.ステップ数が少ないためか,乱数生成器や乱数種の違いを考慮しても再現できているとは言い難いです.とりあえず,今後はステップ数を増やしてみようかと思っていますが,上記のステップ数K=6000の場合でも1時間強かかっているので,C++に移植するなりして対応したいと思います.

参考文献
レプリカ交換モンテカルロ法を含めてMCMCを勉強する際は,「計算統計II」は必読です.「ベイズ統計と統計物理」も分かりやすく書かれています．

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

岩波講座 物理の世界 物理と情報〈3〉ベイズ統計と統計物理

岩波講座 物理の世界 物理と情報〈3〉ベイズ統計と統計物理

力学系については,非数学系の方は"Chaos in dynamical systems", 数学系の方は「力学系」で勉強すると良いのではないかと思います.
Chaos in Dynamical Systems

Chaos in Dynamical Systems

力学系〈下〉

力学系〈下〉