覚え書き:MCMCの収束速度を求めてみよう

 相変わらずMCMCについてぐちぐちと考えている。以下、Haggstrom(2017)「やさしいMCMC入門」(原著2002) からのメモ、第二弾。今回は7-8章からのメモで、MCMCの収束速度について。

 この本の特徴なんだけど、有限状態・離散時間のマルコフ連鎖に話が限定されている。ときどき「これ状態が可算無限だったり非可算無限だったらどうなるんだろうか」と疑問に思う箇所があるんだけど、千里の道も一歩から、である。まあそうやって時間切れで死んでいくんだけどな。

準備. マルコフ連鎖の定常分布
 任意の既約かつ非周期的なマルコフ連鎖は唯一の定常分布を持つ。初期分布がなんであれ、nが十分に大きければ、時点nにおける分布は定常分布に近づく。
 この話は前回のメモで学びました。

無作為q-彩色問題
 グラフ\(G = (V,E)\)の頂点に色を塗る。使ってよい色の数は\(q\)個。隣接している頂点は違う色にしなければならない。これがq-彩色問題である。
 \(q\)個の色の中から一様に色を選ぶことにしよう。解の確率分布(つまり、実行可能な解\(\xi \in \{1, \ldots, q\}^V\)上で一様な確率分布)を\(\rho_{G,q}\)とする。話を簡単にするために、解がひとつはあるということにしよう(実は解がないケースもある。三角形なのに色が2個だけとか)。
 この確率分布から、次のようにサンプリングを行う。

  1. 頂点を\(v_1\)から順繰りに選ぶ。これを\(v\)とする。
  2. \(v\)に隣接している頂点にない色の集合から、一様分布に従って色\(X_{n+1}(v)\)を選ぶ。他の任意の頂点\(w\)については色を変えない、つまり\(X_{n+1}(w) = X_n(w)\)とする。

こういうのを系統的巡回Gibbsサンプラーという。
 このマルコフ連鎖は非周期的だけど、既約かどうかは一概にいえない。でも\(\rho_{G,q}\)は定常分布である(その証明は演習問題になっている。答えをつけておいてほしい…)。

MCMCアルゴリズムの収束速度

定理8.1 \(G=(V,E)\)をグラフとする。\(k\)をグラフの頂点数とし、任意の頂点\(v \in V\)はたかだかd個の隣接頂点しか持たないとする。
 色の数\(q\)について\(q \gt 2d^2\)と仮定する。\(\epsilon \gt 0\)を任意に固定する。
 このとき、任意に固定したq-彩色から出発した、上記の系統的巡回Gibbsサンプラーが、目的の分布\(\rho_{G,q}\)から全変動距離で\(\epsilon\)以内に入ってくるために必要な繰り返し数の上界は次式となる。$$ k \left( \frac{\log k – \log \epsilon – \log d}{\log (\frac{q}{2d^2})} + 1 \right)$$

上の\(q \gt 2d^2\)というのはおおまかな評価で、実はもっと小さくても同じことがいえる。

 さあ、これを証明します。本では小さな字で6ページにわたっている。

 今回もカップリングを使います。
 \(\{1, \ldots, q\}^V\)に値をとる2つのマルコフ連鎖\((X_0, X_1, \ldots), (X’_0, X’_1, \ldots)\)を考える。いずれも上記の無作為q-彩色問題に対する系統的巡回Gibbsサンプラーで、推移行列も同じだが、出発点が異なる。1本目は\(X_0 = \xi\)から始まる。時点\(n\)での分布を\(\mu^{(n)}\)とする。2本目は定常分布\(\rho_{G,q}\)からドローした\(X_0\)から始まる。よってどの時点でも\(\rho_{G,q}\)に従う。
 これから全変動距離\(d_{TV}(\mu^{(n)},\rho_{G,q})\)を評価する。がんばるぞ。

 まず、実装の仕方を決めておく。
 頂点を選んだあとで、その頂点に隣接している頂点にない色の集合から、一様分布に従って色を選びますよね。このときこういう手順をとることにする。まず、あらかじめ、\(\{1, \ldots, q\}\)の異なる順列(\(q!\)個ある)のプールを作っておく。このなかからひとつをドローする。時点\(n\)にドローした奴を\(R_n = (R_n^1, \ldots, R_n^q)\)とする。1本目の連鎖では、\(R_n\)の要素を\(R_n^1\)から順にみていって、最初にみつかった「隣接頂点にない色」を取る。1本目の連鎖でも、同じ \(R_n\) (←ここがポイント) の要素を\(R_n^1\)から順にみていって、最初に見つかった「隣接頂点にない色」を取る。
 このやり方だと、ひとたび連鎖が一致したならば、以降は連鎖はずっと一致することになる。

 [ええええ? そうかなあ? たとえば時点\(n\)で2本の連鎖が頂点\(v\)を訪問し、どちらの連鎖も「周りにない色」として色2を選んだとしますわね。\(X_{n+1}(v) = X’_{n+1}(v) = 2\)である。どちらのグラフでも頂点1の色は2になる。でも連鎖1と連鎖2では、グラフ全体の頂点集合に対する色の割り当てはちがっているわけじゃないですか。だから次の頂点\(v+1\)では、連鎖1では色3が「周りにない色」だが連鎖2では色3が「周りにない色」になっている、ということが起こりうるのではないだろうか。なにか私が問題を誤解しているのだろうか? それとも「ひとたびグラフ全体の彩色が一致したならば以降は連鎖はずっと一致する」ということ? それとも「ひとたび\(\mu^{(n)}\)と\(\rho_{G,q}\)が一致したならば以降は連鎖はずっと一致する」ということ?
 → あー、わかった! \(X_n\)は「時点\(n\)において訪れた頂点の色」じゃなくて、「時点\(n\)における頂点集合の彩色」を表すベクトルなのだ。だからいちいち\(X_{n+1}(v)\)という書き方をしているのだ。「ひとたびグラフ全体の彩色が一致したならば以降は連鎖はずっと一致する」ってことですね。わかったよ先生]

 以降では、ある時点における\(q\)個の色を次の3つに分類する。

  1. タイプ0. どちらの連鎖でも隣接頂点にない色。\(B_0\)個とする。
  2. タイプ1. 一方の連鎖だけで隣接頂点にある色。\(B_1\)個とする。
  3. タイプ2. 両方の連鎖で隣接頂点にある色。\(B_2\)個とする。

 ある時点\(n\)で頂点\(v\)を更新したとき、\(X_{n+1}(v) = X’_{n+1}(v)\)になったら更新が「成功」したと呼び、そうでなかったら更新が「失敗」したと呼ぶことにする。

 まず、1巡目について。つまり時点\(n \leq k\)での更新について。
 失敗する確率は、「順列\(R_n\)をみたとき、タイプ0の色よりもタイプ1の色のほうが先に書いてある確率」、すなわち\(B_1/(B_0+B_1)\)である。これを書き換えていきます。
 $$ \frac{B_1}{B_0+B_1} = \frac{B_1}{q – B_2}$$ \(B_1\)について考えよう。隣接する頂点の数はたかだか\(d\)個だから、両方のグラフで隣接する頂点の延べ数はたかだか\(2d\)個。それらの色はタイプ1か2のどちらかで、そのうちタイプ2の色の頂点は延べ\(B_2\)個ある。だから\(B_1 \leq 2d – 2B_2\)である。$$ \leq \frac{2d – 2B_2}{q-B2} $$ これをさらに緩めて(え、もったいないな…) $$ \leq \frac{2d – B_2}{q-B_2} = \frac{2d(1-\frac{B_2}{2d})}{q(1-\frac{B_2}{q})} $$ 仮定\(q \geq 2d^2\)を、\(q \geq 2d\)と緩めれば、\(1-\frac{B_2}{2d} \geq 1-\frac{B_2}{q}\) だから$$ \leq \frac{2d}{q} $$ ということはですね。一巡した瞬間には、任意の頂点\(v\)について、$$ P(X_k(v) \neq X’_k(v)) \leq \frac{2d}{q}$$が成り立っていることになる。

 2巡目について。つまり時点\(k \lt n \leq 2k\)での更新について。
 今度は見方をちょっと変える。いま訪問した頂点から見た隣接頂点の色が、連鎖の間で異なっていることを「不一致」と呼ぼう。そもそも更新が失敗するのは不一致だからですよね。つまり、P(更新が失敗する) = P(更新が失敗する|不一致) P(不一致)といえる。
 [はいはい。2巡目からは不一致の確率を上から抑えることが出来るというわけね]
 P(更新が失敗する|不一致)は、1巡目と同じく\( \leq \frac{2d}{q}\)である。ではP(不一致)はどうなっているかというと、任意の頂点について色が異なる確率が\( \leq \frac{2d}{q}\)で、隣接頂点の数はたかだか\(d\)個だから、もっとも悲観的にみて、\(\frac{2d}{q}\)を\(d\)倍した値が上限となる。つまり\(\leq \frac{2d^2}{q}\)。
 よって、P(更新が失敗する) \(\leq \frac{2d}{q} \left( \frac{2d^2}{q} \right)\)。
 ということはですね、2巡した瞬間には、任意の頂点\(v\)について、$$ P(X_{2k}(v) \neq X’_{2k}(v)) \leq \frac{2d}{q} \left( \frac{2d^2}{q} \right)$$
 が成り立つわけです。

 これをさらに続けていくとこうなる。\(m)\)巡し終わった瞬間に $$ P(X_{mk}(v) \neq X’_{mk}(v)) \leq \frac{2d}{q} \left( \frac{2d^2}{q} \right)^{m-1} = \frac{1}{d} \left( \frac{2d^2}{q} \right)^m$$

 今度は、時点\(mk\)で2つの連鎖\(X_{mk}, X’_{mk}\)が正確に同じ配置を持つのに失敗する確率に注目する。つまり、\(k\)個の頂点のうち少なくともひとつの頂点\(v\)において\(X_{mk}(v) \neq X’_{mk}(v)\)である確率、である。$$ P(X_{mk} \neq X’_{mk}) \leq \sum_v P(X_{mk}(v) \neq X’_{mk}(v)) = \frac{k}{d}\left( \frac{2d^2}{q} \right)^m $$

 いよいよ全変動距離について。任意の部分集合\(A\)について考える (原文には任意の\(A \subset \{1, \ldots, k\}^V\)とあって首を捻ったが、\(A \subset \{1, \ldots, q\}^V\)の誤り、つまり任意の彩色の集合のことだろう)。
 全変動距離とは$$ d_{TV}(\mu^{(mk)}, \rho_{G,q}) = \max_A |\mu^{(mk)}(A) -\rho_{G,q}(A)| $$ のことであった。書き換えると$$ = \max_A |P(X_mk \in A) – P(X’_{mk} \in A)|$$
 絶対値記号の中味を取り出す。$$ P(X_{mk} \in A) – P(X’_{mk} \in A) $$ さらに場合分けする。$$ = P(X_{mk} \in A, X’_{mk} \in A) + P(X_{mk} \in A, X’_{mk} \notin A) – P(X_{mk} \in A, X’_{mk} \in A) – P(X_{mk} \notin A, X’_{mk} \in A)$$ 「どっちも入っている」項が2つ出てきたので消して $$ = P(X_{mk} \in A, X’_{mk} \notin A) – P(X_{mk} \notin A, X’_{mk} \in A)$$ 2項目を消して緩める。(ええええ、いいの? もったいないなあ) $$ \leq P(X_{mk} \in A, X’_{mk} \notin A) $$ ってことは、少なくとも同じじゃないって事だから、さらに緩めて $$ \leq P(X_{mk} \neq X’_{mk}) = \frac{k}{d} \left( \frac{2d^2}{q} \right)^m $$ 同様に $$ P(X’_{mk} \in A) – P(X_{mk} \in A) \leq \frac{k}{d} \left( \frac{2d^2}{q} \right)^m$$ ひっくるめて $$ |P(X’_{mk} \in A) – P(X_{mk} \in A)| \leq \frac{k}{d} \left( \frac{2d^2}{q} \right)^m$$ ここから全変動距離は $$ d_{TV}(\mu^{(mk)}, \rho_{G,q}) \leq \frac{k}{d} \left( \frac{2d^2}{q} \right)^m$$ となる。

 もう一息だぜ。
 右辺を\(\epsilon\)より小さくする\(m\)を求めたいわけだから、$$ \frac{k}{d} \left( \frac{2d^2}{q} \right)^m = \epsilon$$ を解けば良い。本には書いてないけどやってみます。両辺の対数をとって $$ \log k – \log d + m (\log 2d^2 – \log q) = \log \epsilon$$ $$ m(\log q – \log(2d^2)) = \log k – \log \epsilon – \log d $$ $$ m = \frac{\log k – \log \epsilon – \log d}{\log(\frac{q}{2d^2})} $$ これが、必要な巡回数\(m\)の上界を表している。マルコフ連鎖の実行数\(n\)の上界に換算するにあたっては、これを\(k\)倍し、さらに安全のために\(k\)回足してやる。よって答えは$$ k \left( \frac{\log k – \log \epsilon – \log d}{\log(\frac{q}{2d^2})} +1 \right)$$である。(証明終わり)

 いやー… 数学力がないもので、こういう話って、ほんとによくわかんない。正直、私にはついていけない。
 証明の途中でどんどん不等式をつないで、制約を緩めていくじゃないですか。そこがすごく恣意的に思えるのである。最終的に定理を証明できさえすればいいのだから、それでかまわないんだろうけれど… もっときちんと考えれば、収束に必要な実行数の上界はもっと小さいことがわかるんじゃないの??? もっと真面目にやろうよ??? と思えてしまう。ううむ。私の頭が悪いからなのかなあ。

 まあいいや、考え方がわかったので、本件に関してはこれで良しとしよう。この無作為q-彩色問題というのはレアケースであって、実務的には、マルコフ連鎖の収束速度なんてなかなかわかんないんじゃないかと思うし。
 えーっと、著者によれば、マルコフ連鎖の均衡状態への速い収束を証明するために、ここで扱ったカップリング法のほか、固有値限界、パスと流れの議論、異なる連鎖間の多様な比較、強定常双対性、などなど、いろんなテクニックが開発されているのだそうだ。研究者のみなさん、がんばってください。なにかわかったら結果だけ教えて下さいな。