この記事では、サンプリング手法のひとつであるメトロポリス法について、その基本原理と実装例を紹介します。最近注目されている画像生成AIなどの分野でサンプリングの話題を目にする機会が増えていますが、実際の実装となると意外と複雑です。そこで、基本となるマルコフ連鎖モンテカルロ法(MCMC)に焦点をあて、その１つの計算法であるメトロポリス法を通して理解を深めていきます。
なお、本記事は逆関数法の記事の続編という位置付けになります。
- 本記事の作成には以下の情報を大変参考にさせて頂きました
導入
逆関数法では確率分布関数から累積分布関数を算出し、逆変換により確率変数$x$の等式にすることで、一様乱数から確率分布関数に従う乱数を生成できることを確認しました。
実際には逆関数が得られないケースは多々あります。
例として１次元の混合ガウス分布を考えてみます。
p(x) = \sum_{k=1}^{K} \pi_k \cdot \mathcal{N}(x | \mu_k, \sigma_k^2)
K：混合成分の数
$\pi_i$：各成分の混合比$（\sum_{k} \pi_k = 1)$
$\mathcal{N}(x | \mu_k, \sigma_k^2)$：平均$\mu_i$、分散 $\sigma_i^2$のガウス分布
この時、正規分布の確率密度関数（PDF）は次の式になります。
\mathcal{N}(x | \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(x - \mu)^2}{2\sigma^2} \right)
これを0,1標準化すると以下のような累積分布関数となります。
\Phi(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{-t^2 / 2} dt
この逆関数は解析的に表せません。(近似計算はできようですが、誤差関数や級数展開などを利用する複雑な計算になる模様です。)
一方で、混合ガウスの累積分布関数は
F(x) = \sum_{i=1}^{K} \pi_i \Phi\left(\frac{x - \mu_i}{\sigma_i}\right)
となります。単一ガウスですら解析的には表せませんから、より難しいことが示唆されますね。
このように逆関数の形を得ることは一般的には容易ではありません。ただ、目的は確率分布に従うサンプリングです。
そのため、逆関数を使わない別の方法を考えます。その１つがMCMCになります。
※混合ガウス分布の近似分布を、確率分布の積の形で記述し、逐次的な形で解析的に近似解を得る方法があります。
これらは、EMアルゴリズムや変分ベイズ法と呼ばれ、これを利用してサンプリングすることも可能です。
筆者のこちらの記事ではEMアルゴリズムを利用して、モデルを比較することに取り組んでいます。
※ガウス分布は、ボックス＝ミュラー法と呼ばれる手法により一様分布を起点にサンプリングできることが知られています。こちらの説明が分かりやすいと思います。
MCMCの概要
この記事では、MCMC（Markov Chain Monte Carlo）が活躍する背景を説明した後、具体的な利用例としてのボルツマンマシンへの適用を通じて、MCMCの基本概念と代表的なアルゴリズム（メトロポリス法）について解説します。
MCMCが必要とされる背景
多くの応用分野で、分布の正規化定数（分母の計算）が解析的に求められない問題が存在します。ここでは主な例を2つ挙げます。
- ベイズ推定の場合
パラメータの事後分布を計算するケース。データ$x$に対するパラメータの事後分布$p(θ|x)$を知りたい時。
ベイズの定理より以下の式になります。
p(\theta|x) = \frac{p(x,\theta)}{\int p(x,\theta) d\theta}
ベイズ推定のケースでは、右辺の分子は自分たちが作ったモデルになるので、関数の形（関数の掛け算）を知ることができます。
一方、分母はその関数の積をとることになり、この関数の形を得ることは困難な場合が多く、その場合は事後分布も解析的な形式が得られません。
2. ボルツマンマシンの場合
ボルツマンマシンのモデルは以下のように表現されます。
p(x|\theta) = \frac{1}{Z(\theta)}\exp(-\Phi(x,\theta))
Z(\theta)=\sum_x\Phi(x,\theta)
右辺の分子はモデルとして構築できますが、分母に現れる正規化定数は、$x$のあらゆる状態の和を計算する必要があります。
例えば、$x=(x_1,x_2,..,x_N)$,$x_n=\{0,1\}$のケースを考えたとき、
$2^N$通りの和を計算する必要があり、計算コストが非常に大きくなります。
MCMCの役割と基本概念
上記のように、正規化定数の計算が困難なケースでは、直接解析的なサンプリングが難しいため、MCMC が非常に有効です。
MCMCは、マルコフ連鎖を用いて、対象の分布(たとえば$p(x|\theta)$,$p(\theta|x)$)に収束するサンプルを生成する手法です。具体的には、任意の初期状態から始め、遷移確率$W(x'|x)$に従って状態を逐次更新していくことで、最終的に目標とする定常分布に収束させます。
収束後のサンプルは、その定常分布に従うとみなすことができるため、正規化定数の計算が困難な分布（正規化されていない分布）のサンプリングが可能となります。
p_0(x|\theta)→p_1(x|\theta)→・・・→p_T(x|\theta)
釣り合い条件
概要
MCMCにより定常分布を見つける方法ですが、定常分布のマスター方程式に従う分布を見つけることに相当します。
まず非平衡定常マスター方程式は、以下の形になります。
\frac{dP(x,t)}{dt} = \sum_{x'} \left[ W(x|x')P(x',t) - W(x'|x)P(x,t) \right]
ここで$W(x'|x)$は以下の制約条件があります。
\sum_{x'}W(x'|x)=1
離散時間を考えると以下の式になります。
P(x, t+1) - P(x, t) = \sum_{x'} \left[ W(x|x') \, P(x', t) - W(x'|x) \, P(x, t) \right]
この式は、ある状態$x$の確率分布の変化は、ある状態$x$に対する「流入（gain）」と「流出（loss）」のバランスから成り立つことを明示しています。
さて、$P(x, t)$を右辺に移行し、$\sum_{x'}W(x'|x)=1$の積を取ると、左辺の第2項と打ち消しあい、以下の式になります。
P(x, t+1) = \sum_{x'} W(x|x')P(x', t)
この式は$x$の状態が遷移確率によって更新されていくことを示しています。
さらに平衡定常分布の条件に絞ると、方程式は$t$に依存しない形となります。
P(x) = \sum_{x'} W(x|x')P(x')
この式は、時間が経っても保存することから、この式が成り立つ条件は釣り合い条件と呼ばれます。
釣り合い条件の左辺を目標分布とすると、この関係式を満たす$W(x|x')$と$P(x')$の元でのサンプリングが課題となります。
これを実現する１つの方法が次に示す詳細釣り合い条件を満たすようなサンプリングになります。
より深い理解のための補足
上式は１つの状態$x$に着目していました。一般的に状態$x$はN個あると考えられます。すべての状態を$\textbf{x}=(x_1,x_2,...,x_N)$とすると、この時$P(\textbf{x})$はベクトルとなり、状態遷移は行列$W(\textbf{x})$となります。
この時、マスター方程式は以下の形になります。
P(\textbf{x}) = W(\textbf{x})P(\textbf{x})
となります。これは固有値１の固有ベクトル$P(\textbf{x})$を求める、固有値問題に相当します。つまり、上記のマスター方程式は固有値問題の理論によっても支えられていると考えることができます。詳細はこのサイトなどを参考にすると良いと思います。
詳細釣り合い条件
MCMCが定常分布に収束するための十分条件の1つとして詳細つり合い条件があります。（詳細釣り合い条件を満たすならばマスター方程式を満たすが逆は成り立たない、ということです。）
詳細つり合い条件は以下のような式になります。
W(x'|x)P(x) = W(x|x')P(x')
実際に、右辺の$W(x|x')P(x')$に関して、$\sum_{x'}$で和をとり、左辺に変換すると、釣り合いの条件を満たします。
よって詳細釣り合い条件を満たすような関数の形を作れば良いことになります。
ボルツマンマシンの詳細釣り合い条件
ボルツマンマシンを詳細釣り合い条件に適用します。この時以下の式が成り立ちます。
\frac{W(x'|x)}{W(x|x')} = \frac{e^{-E(x')}}{e^{-E(x)}} = \exp\left( -\left[ E(x') - E(x) \right] \right)
ここでは、規格化定数$Z$がキャンセルされ、計算の必要がなくなります。
この設計により、メトロポリス法やギブスサンプリングなどの具体的なアルゴリズムが実現されています。
メトロポリス法
ここでは、代表的なMCMCアルゴリズムであるメトロポリス法について解説します。
アルゴリズムの概要
メトロポリス法は、以下の2つの要素で遷移確率を構成します。
1.提案確率
現在の状態から次の候補状態を提案する確率。
2.受理確率
提案された状態を受け入れるかどうかを判断する確率。
受理確率は、エネルギー（あるいは対数尤度）の差に基づいて計算され、エネルギーが上がる場合でも、一定の確率で採択することで局所解からの脱出を可能にします。
ボルツマンマシンとメトロポリス法
メトロポリス法により、詳細釣り合い条件を満たすように遷移確率を設計することで、目標とするボルツマン分布からのサンプリングする方法を示します。アルゴリズムの手順は以下の通りです。
1.初期状態の設定
任意の状態$x$から開始する。
2.候補状態の提案
提案確率$q(x'|x)$に従って、現在の状態$x$から状態$x'$を提案する。
3.受理確率の計算
受理確率$A(x→x)$を次のように定義する。
A(x \to x') = \min\left\{ 1, \frac{q(x|x')}{q(x'|x)} \exp\left( -\left[ E(x') - E(x) \right] \right) \right\}
さて、上式ですが対称な提案確率（$𝑞(𝑥′∣𝑥)=𝑞(𝑥∣𝑥′)$）の場合、受理確率はより単純な式になります。
A(x \to x') = \min\left\{ 1, \exp\left( -\left[ E(x') - E(x) \right] \right) \right\}
この式について考察してみます。
$E(x')>>E(x)$では、$\exp\left( -\left[ E(x') - E(x) \right] \right)〜0$、$A(x \to x')〜0$となります。また$E(x') < E(x)$では、$\exp\left( -\left[ E(x') - E(x) \right] \right)>1$、$A(x \to x')=1$となります。
これより受託確率は、エネルギーの差により変化すること、非対称な性質を持つことが確認できました。
4.状態更新
提案された状態$x'$を、確率$𝐴(𝑥→𝑥^′)$で受理し、受理された場合は状態を$x'$に更新する。受理されなければ状態$𝑥$を維持する。
5.繰り返し
上記の手順を繰り返し、マルコフ連鎖が平衡状態に達したとみなし、その後のサンプルを目標分布からのサンプルとして利用する。
詳細釣り合い条件との関係
受理確率の設計は、詳細釣り合い条件を満たすために重要な役割を果たします。具体的には、前述の詳細釣り合い条件
\frac{W(x'|x)}{W(x|x')} = \frac{P(x')}{P(x)} = \exp\left( -\left[ E(x') - E(x) \right] \right)
となるように、受理確率$A(x→x')$がエネルギー差によって決定されます。つまり、エネルギーが低い方向への遷移が有利でありながら、エネルギーが高い方向への移行も一定の確率で許容することで、逆方向の遷移とのバランスが取られ、最終的に平衡状態（ボルツマン分布）に収束します。
ここで再度受理確率を確認します。
A(x \to x') = \min\left\{1, \exp\left(-\left[E(x') - E(x)\right]\right)\right\}
この式と詳細釣り合い条件の式を比べると、エネルギーが上がる方向を全て拒否すると、$A(x \to x')$=1のみの状態しか採択しない、つまり逆の方向の状態遷移が行われず、$W(x'∣x)$がゼロになり、詳細釣り合い条件の比が成立しなくなります。
逆に、エネルギーが低下する場合は、受理確率を1に制限することで確率としての範囲を守りながら、詳細釣り合い条件と整合する設計になっています。
おわりに
最後まで読んでいただきありがとうございました。
本記事を執筆する動機は、拡散モデルに向かって、一通りのサンプリング法の知識をまとめていくことにあります。MCMCは、正規化定数が求められない確率モデルからサンプリングする手法ですが、その手法を実現する際に、マスター方程式を解くという問題があり、またマスター方程式が固有値問題とも関連することから、理論を明らかにするには、より深い数学の知識が必要になると、執筆して感じました。
これらの理論を踏まえてアップデートすることや、図などを導入し、もう少し視覚的に読みやすくする工夫は今後の課題とします。
知識至らずで曖昧な点があると思います。ご指摘等ありましたら、編集いたしますので、コメントお待ちしております。