2020-08-12

ウェーブレット変換

　ウェーブレット変換は、フーリエ変換のように、時間座標を持つ関数 $f(t)$ を基底関数の和で表現する方法です。しかし、フーリエ変換とは異なり、ウェーブレット変換では時間方向に局在した基底関数を用います。それによって、時間領域に応じて変化する関数の振る舞いを解析できます。

　フーリエ変換と同様に、連続ウェーブレット変換を離散化すると離散ウェーブレット変換になりますが、離散ウェーブレット変換では、正規直交基底を選ぶことができるほか、周波数に応じて基底展開の「解像度」が上がっていくといった便利な性質を持ちます。それによって、関数 $f(t)$ を、振る舞いの周波数に応じて解像度を変えながら分解していく「多重解像度解析」を行うことができます。

　そのため、今回は先に離散ウェーブレット変換について説明します。また、ウェーブレット変換は再生核ヒルベルト空間を構成する表現になっているので、その点にも簡単に触れます。

フーリエ変換とウェーブレット変換
離散ウェーブレット変換
多重解像度解析
- 関数空間と多重解像度表現
- 2スケール方程式
連続ウェーブレット変換
再生核ヒルベルト空間（RKHS）
リースの表現定理

フーリエ変換とウェーブレット変換

　フーリエ変換の場合、基底関数は周波数が異なる多数の波でした。基底関数の一つを $\psi (t) = e^{2\pi it}$ としたとき、基底関数全体は、 $a\in \mathbb{R} \backslash \{0\}$ に対して、 $\psi(\frac{t}{a})$ に定数項を加えたものです。各基底関数の周波数は $\frac{1}{a}$ です。

　フーリエ変換では、基底関数が時間方向に対して一様に広がっているので、各時間領域に応じて振る舞いが変化するような関数の特徴を抽出することが困難でした。 $t$ を一定区間ごとに窓関数で区切って、各々でフーリエ変換を行うこともできますが、そうすると窓関数のスケールより小さいスケールの振る舞いは平均化されてしまいますし、窓関数より大きいスケールの振る舞いは分断されてしまいます。

　ウェーブレット変換では、基底関数として、 $\psi (t) = (1-t^2)e^{-\frac{t^2}{2}}$ のように、時間に対して一様ではなく、中心付近に局在しているようなものをベースに選びます。そして、基底関数全体を $b\in \mathbb{R} \backslash \{0\}, a\in \mathbb{R}$ に対して $\psi \left( \frac{t-b}{a} \right)$ のように取ります。

　フーリエ変換の場合と比較して、 $t$ 方向のシフト量 $b$ が追加されています。フーリエ変換では基底関数が一様だったので、時間方向にシフトさせた基底関数を別に考える意味はありませんでした。実際、フーリエ変換における時間方向のシフトは、展開の係数の位相成分にしか現れません。

　ウェーブレット変換では基底関数が局在しているので、時間方向のある領域に局在した特徴を捉えることができるようになります。スケール因子 $a$ と時間シフト因子 $b$ の組み合わせによって、関数 $f(t)$ が持つ様々なスケールと時間領域の振る舞いを解析することができます。

離散ウェーブレット変換

　関数 $f(t)$ を、整数 $j,k \in \mathbb{Z}$ でラベル付けた正規直交基底 $\psi _{j,k} (t)$ によって次のように級数展開することを考えます。

$\begin{align} f(t) = \sum_{j,k} {d_{j,k}\psi_{j,k}(t)} \end{align} \tag{1}$

　基底の正規直交条件は次のように表せます。

$\begin{align} \int_{-\infty}^{\infty} {\psi^*_{j,k}(t) \psi_{j',k'}(t) dt} = \delta _{j,j'} \delta_{k,k'} \end{align} \tag{2}$

　基底関数 $\psi _{j,k} (t)$ は、ある関数 $\psi (t)$ を使って次のように構成します。

$\begin{align} \psi_{j,k}(t) = 2^\frac{j}{2} \psi (2^j t - k) \end{align} \tag{3}$

　基底関数系 $\psi_{j,k}(t)$ を生成する元になる関数 $\psi(t)$ は、マザーウェーブレットと呼ばれます。 $2^\frac{j}{2}$ は正規化係数です。マザーウェーブレットは、直流成分を含まない、次のような条件を満たすべきとされています。その理由は、連続ウェーブレット変換の項目で説明します。

$\begin{align} \int_{-\infty}^{\infty} {\psi(t)dt} = 0 \end{align} \tag{4}$

　 $f(t)$ の級数展開の係数 $d_{j,k}$ は、式(1)、式(2) より、次のように表すことができます。これが、離散ウェーブレット変換です。

$\begin{align} d_{j,k} = \int_{-\infty}^{\infty} {f(t)\psi^*_{j,k}(t) dt} \end{align} \tag{5}$

　式(1) は、離散ウェーブレット変換に対する逆変換になっています。

　マザーウェーブレットの最も簡単な例として、ハール・ウェーブレット $\psi_H(t)$ があります。

$\begin{align} \psi_H(t) = \begin{cases} 1 & (0 < t < \frac{1}{2}) \\ -1 & (\frac{1}{2} \leq t < 1) \\ 0 & (\mbox{otherwise}) \end{cases} \end{align} \tag{6}$

これが式(2) を満たしていることはすぐに分かります。

多重解像度解析

　関数 $f(t)$ を級数展開するにあたって、先程はマザーウェーブレット $\psi (t)$ を元にして基底関数系を構成しました。これでも級数展開は可能なのですが、マザーウェーブレットは式(4) の条件に拘束されているため、非効率な級数表現になっている可能性があります。また、周波数に対応する値 $2^j$ が大きくなるほどに、解像度が高くなっていく（ $f(t)$ をより細かく表現できるようになる）べきですが、そのことが明示されていません。そのような性質を持つ表現は、多重解像度表現と呼ばれます。

　そこで、上記の $2^j$ に応じて解像度が高くなっていくという条件、すなわち多重解像度表現を満たしつつ、マザーウェーブレットを補助して効率的な表現を得られるような追加の基底関数系を考えます。追加の基底関数系を生成するベースの関数を、ファーザーフェーブレット（またはスケーリング関数）と呼びます。

　ファーザーウェーブレット $\phi (t)$ を用いて、マザーウェーブレットの場合と同じように基底関数系を生成します。

$\begin{align} \phi_{j,k}(t) = 2^\frac{j}{2} \phi (2^j t - k) \end{align} \tag{7}$

　正規直交条件もマザーウェーブレットと同じように定義します。さらに、マザーウェーブレットが生成する基底関数系とフェーザーウェーブレットが生成する基底関数系の間にも、正規直交条件を課します。すなわち、

$\begin{align} \int_{-\infty}^{\infty} {\psi^*_{j,k}(t) \phi_{j',k'}(t) dt} = \delta _{j,j'} \delta_{k,k'} \end{align} \tag{8}$

となります。

関数空間と多重解像度表現

　多重解像度表現を定義するために、ファーザーウェーブレットによる基底関数系の $j$ をある値に固定した場合に表現可能な関数 $f_{j}(t)$ 全体の集合を $V_j$ とします。 $V_j$ は次のように表せます。 $L^2(\mathbb{R})$ は2乗可積分関数全体です。

$\begin{align} V_j = \left\{ f_j (t)\in L^2(\mathbb{R}) \middle| \exists c_{j,k}\in \mathbb{C}, f_j (t) = \sum_k c_{j,k} \phi_{j,k} (t)\right\} \end{align} \tag{9}$

　これから、多重解像度表現は次のように定義できます。

$\begin{align} \cdots \subset V_{-2} \subset V_{-1} \subset V_0 \subset V_1 \subset V_2 \subset \cdots \end{align} \tag{10}$

　ある $J_0$ に対して、基底 $\phi_{J_0,k}$ の $k$ に関する和で表せる関数は全て、 $\phi_{J_0+1,k}$ の和でも表せるということを意味しています。逆は成り立ちません。 $j$ が大きくなるほどよりたくさんの関数を表現できて、 $j$ が小さくなるとほとんど表現できなくなってしまうということです。その極限として、 $V_{\infty} = L^2(\mathbb{R})$ 、 $V_{-\infty} = \{0\}$ というものも多重解像度表現の条件です。

　フェーザーウェーブレットのもう一つの条件は、マザーウェーブレットに関係します。上記と同様に、 $j$ を固定した場合にマザーウェーブレットが生成する基底 $\psi_{j,k}$ の張る空間を $W_j$ とするとき、次の関係が成り立つようにします。

$\begin{align} V_{j+1} = V_{j} \oplus W_{j} \end{align} \tag{11}$

　ファーザーウェーブレットが張る空間とマザーウェーブレットが張る空間が直交し、その和空間が、一つ上の $j$ においてファーザーウェーブレットが張る空間に等しくなります。つまり、マザーウェーブレットは、解像度が一段階変化したときの表現能力の差分とみなせます。式(11) と $V_{\infty} = L^2(\mathbb{R})$ から、

$\begin{align} L^2(\mathbb{R}) = V_{J_0} \oplus \bigoplus_{j=J_0}^{\infty} {W_{j}} \end{align} \tag{12}$

が成り立ちます。これから、マザーウェーブレットとファーザーウェーブレットを用いた多重解像度表現が得られます。

$\begin{align} f(t) = \sum_{k} c_{J_0, k}\phi_{J_0, k}(t) + \sum_{j\geq J_0,k} d_{j,k}\psi_{j,k}(t) \end{align} \tag{13}$

　関数 $f(t)$ の概形が $\phi_{J_0, k}(t)$ で抑えられるので、式(1) のマザーウェーブレットによる級数展開と比べて少数の係数で表現ができ、 $j$ の値を大きくするたびに、解像度が高くなっていく表現になっています。

2スケール方程式

　最後に、ファーザーウェーブレットを構成する条件をまとめます。 $V_{j},W_{j}$ が直交するので、式(8) 等の直交条件が必要です。また、式(10) と式(11) を合わせて次のように書けます。

$\begin{align} \phi(t) = \sum_{n}{g_0(n)\sqrt{2}\phi(2t-n)} \\ \psi(t) = \sum_{n}{g_1(n)\sqrt{2}\phi(2t-n)} \end{align} \tag{14}$

　式(14) は、2スケール方程式と呼ばれます。係数は次の式で決まります。

$\begin{align} g_0(n) = \int_{-\infty}^{\infty}{\phi(t)\sqrt{2}\phi^*(2t-n)dt} \\ g_1(n) = \int_{-\infty}^{\infty}{\psi(t)\sqrt{2}\phi^*(2t-n)dt} \end{align} \tag{15}$

　また、式(11) を満たすためには、例えば $\phi(2t) = \frac{\phi(t)+\psi(t)}{2}$ のように、 $\phi(2t)$ が $\phi(t-k)$ と $\psi(t-k)$ の和で表せることが必要です。ハール・ウェーブレット（ $\psi_H(t)$ は式(6) で与えられ、 $\phi_H(t)$ は、 $0 < t < 1$ でのみ値 1 を取り、他では 0 とする）やシャノン・ウェーブレット $\psi_S(t) = \frac{\sin{2\pi t} - \sin{\pi t}}{\pi t}, \phi_S(t) = \frac{sin{\pi t}}{\pi t}$ は、これらの条件を全て満たしています。

連続ウェーブレット変換

　連続ウェーブレット変換は、マザーウェーブレット $\psi(t)$ によって、次のように定義します。 $W_{\psi}(a,b)$ の $a$ がスケール方向、 $b$ が時間方向の成分を表します。

$\begin{align} W_{\psi}(a,b) = \frac{1}{\sqrt{|a|}} \int_{-\infty}^{\infty}{f(t) \psi^*\left( \frac{t-b}{a} \right) dt} \end{align} \tag{16}$

　 $\frac{1}{\sqrt{|a|}}$ は正規化係数です。マザーウェーブレットが $\int_{\infty}^{\infty} {|\psi(t)|^2 dt} = 1$ であることによります。逆変換は、下記のようになります。

$\begin{align} f(t) = \frac{1}{C_{\psi}} \int_{-\infty}^{\infty} {W_{\psi}(a,b) \frac{1}{\sqrt{|a|}} \psi\left( \frac{t-b}{a}\right) \frac{da db}{a^2}} \end{align} \tag{17}$

ただし、

$\begin{align} C_{\psi} = \int_{-\infty}^{\infty} {\frac{|\tilde{\psi}(\omega)|^2}{|\omega |}} \end{align} \tag{18}$

で、 $\tilde{\psi}(\omega)$ は $\psi(t)$ のフーリエ変換です（変換基底を $e^{-2\pi i\omega t}$ とした場合の表現）。式(17) の逆変換が存在するためには、 $C_{\psi}$ が有限の値を持たなければなりません。このことを、アドミッシブル条件と言います。 $\psi(t)$ 自体は2乗可積分なので、アドミッシブル条件は、 $\tilde{\psi}(0)=0$ とほとんど同じことです。これは、式(4) でみた条件 $\int_{-\infty}^{\infty} {\psi(t)dt} = 0$ になっています。

　式(17) が連続ウェーブレット変換の逆変換になっていることを確認します。式(16) を式(17) の右辺に代入し、右辺を直接計算した結果が $f(t)$ に戻ってくることを見ます。

$\begin{align*} \frac{1}{C_{\psi}} \int{f(t')\psi^* \left( \frac{t'-b}{a}\right) \psi \left( \frac{t-b}{a}\right) dt' db \frac{da}{|a|^3}} \\ = \frac{1}{C_{\psi}} \int{\tilde{f}(\omega)e^{2\pi i\omega t}\psi^* \left( \frac{t'-b}{a}\right) \psi \left( \frac{t-b}{a}\right) dt' db \frac{da}{|a|^3} d\omega} \\ \\ \mbox{ここで, }t'=t+(u-v)a, b=t-av \mbox{ と変数変換.} \\ \\ \frac{1}{C_{\psi}} \int{\tilde{f}(\omega)e^{2\pi i\omega t + 2\pi i\omega au - 2\pi i\omega av}\psi^* \left( u \right) \psi \left( v \right) du dv \frac{da}{|a|} d\omega} \\ = \frac{1}{C_{\psi}} \int{\tilde{f}(\omega)e^{2\pi i\omega t} \tilde{\psi}^* \left( a\omega \right) \tilde{\psi} \left( a\omega \right) \frac{da}{|a|} d\omega} \\ = \frac{1}{C_{\psi}} \int{\tilde{f}(\omega)e^{2\pi i\omega t} C_{\psi} d\omega} \\ =f(t) \end{align*} \tag{19}$

計算の途中で、フーリエ変換と逆フーリエ変換を何度か行っています。積分の交換をはじめとして、このような式変形ができるのは、フビニの定理によります。従って、積分の中の関数を多変数関数とみたときに、(絶対)可積分である必要があります。

再生核ヒルベルト空間（RKHS）

ヒルベルト空間

　再生核ヒルベルト空間は、ヒルベルト空間の一つです。ヒルベルト空間 $\mathcal{H}$ は、内積が定義された完備距離空間でした。関数空間の場合、 $f,g\in \mathcal{H}$ の内積は

$\begin{align} \langle f, g \rangle = \int f(t) g^*(t) dt \end{align} \tag{20}$

で、距離は $||f-g||^2 = \langle f-g, f-g \rangle$ となります。距離空間なので、任意のコーシー列が $\mathcal{H}$ 自身に収束すれば完備です。

再生核ヒルベルト空間の概要

　さて、集合 $\Omega$ でラベル付けたヒルベルト空間上の線形汎関数 $L_t(f)=f(t)$ を考えます。 $t\in \Omega$ です。 $L_t: \mathcal{H}\rightarrow \mathbb{C}$ となっています。

　この $L_t$ が任意の $t\in \Omega$ に対して $\mathcal{H}$ 上で連続写像になっているとき、 $\mathcal{H}$ を再生核ヒルベルト空間（RKHS）と呼びます。*1

　関数 $f, g \in \mathcal{H}$ が近ければ、各点での値 $f(t), g(t)$ も近くならなければなりません。例えば、ある測度0の点 $t_0 \in \Omega$ でだけ $f(t_0) \neq g(t_0)$ で、他の点では $f(t) = g(t)$ となっているような場合は、 $L_{t_0}$ は連続写像になっておらず、再生核ヒルベルト空間にはなりません。

　再生核ヒルベルト空間では、 $L_t$ は連続線形汎関数になっています。すると、リースの表現定理（後述）より、任意の $f\in \mathcal{H}$ に対して次の式が成り立つような $K_t \in \mathcal{H}$ が、各々の $t$ に対してただ一つずつ存在します。

$\begin{align} L_t(f) = \langle f | K_t \rangle \end{align} \tag{21}$

　明示的に書き直すと、次のようになります。

$\begin{align} f(t) = \int f(t') K^*_t(t') dt' \end{align} \tag{22}$

　 $k(t',t) = K_t(t')$ が再生核で、カーネル関数とも呼ばれます。 $k: \Omega \times \Omega \rightarrow \mathbb{C}$ です。`

　式(22) から、次のように再生核の一意性を直接確認することもできます。
　再生核を $K_t(t'), K'_t(t')$ とおいて、 $f(t)$ として再生核自身 $K_s(t)$ を与えれば、 $K_s(t) = \int K_s(t') K^*_t(t') dt' = \int K_s(t') K'^*_t(t') dt'$ となりますが、二項目は $K^*_t(s)$ に等しく、三項目は $K'^*_t(s)$ に等しいので、 $K_t = K'_t$ となります。

　また、式(22) から、任意の $f\in \mathcal{H}$ に対して $\int f(t) k(t', t) f^*(t') dt dt' = \langle f | f \rangle \geq 0$ なので、カーネル関数 $k(t', t) = K_t(t')$ は半正定値を取ることが分かります。

ウェーブレットにおける再生核ヒルベルト空間

　
　正規直交基底を用いる離散ウェーブレット変換では、式(1) と式(5) から、例えば下記がカーネル関数になっており、その時の再生核ヒルベルト空間は（理想的には）二乗可積分関数全体です。

$\begin{align} k(t', t) = \sum_{j,k} \psi_{j,k}(t') \psi^*_{j,k}(t) \end{align} \tag{23}$

　それだけでなく、式(10) でみたような部分空間 $V_j$ も再生核ヒルベルト空間です。式(11) の $W_j$ もまた、別の再生核ヒルベルト空間です。これらの再生核は、基底関数 $\phi_{j,k}(t)$ や $\psi_{j,k}(t)$ の $j$ を固定して、式(23) のように和を取ったものになります。

　連続ウェーブレット変換では、例えば式(19) から、下記がカーネル関数になっていることが分かります。

$\begin{align} k(t', t) = \frac{1}{C_{\psi}} \int \psi \left( \frac{t'-b}{a}\right) \psi^* \left( \frac{t-b}{a}\right) db \frac{da}{|a|^3} \end{align} \tag{24}$

リースの表現定理

　再生核ヒルベルト空間の構成に必要な、リースの表現定理を証明します。リースの表現定理は、関数を実数に移す汎関数による写像は、実は特別な一つの関数との内積による写像として表現できることを示しています。これは、汎関数の具体的な表現方法を与える、非常に有用な定理です。*2

　リースの表現定理は次のようなものです。簡単のため、実数の空間を考えます。

【リースの表現定理】
　ヒルベルト空間 $\mathcal{H}$ 上の連続な線形汎関数を $L$ とする。 $L$ に対し、 $h_L \in \mathcal{H}$ がただ一つ定まり、任意の $f\in \mathcal{H}$ に対し、下記の式が成り立つ。
$\begin{align} L(f) = \langle f, h_L \rangle \end{align} \tag{24}$

　【証明】

　全ての $f \in \mathcal{H}$ に対して恒等的に $L(f)=0$ である場合、 $h_L = 0$ とすれば良いです。以下、そうでない場合を考えます。

　 $\mathcal{K} = \left\{ f \in \mathcal{H} \middle| L(f)=0 \right\}$ とします。

　 $\mathcal{K} = L^{-1}(0)$ で $\{0 \}$ は閉部分空間、 $L$ は連続写像なので、 $\mathcal{K}$ も閉部分空間です。直交分解定理より、 $g \in \mathcal{H} \backslash \mathcal{K}$ であって、 $\mathcal{K}$ の全ての元と直交する（内積が0になる）ような $g$ が存在します。

　ここで、次のような $u \in \mathcal{H}$ を考えます。

$\begin{align} u = L(f) \frac{g}{||g||} - L\left( \frac{g}{||g||} \right) f \end{align} \tag{25}$

　 $L$ が線形汎関数であることから、 $L(u) = 0$ となります。従って、 $u \in \mathcal{K}$ です。すると、 $\langle u, g \rangle = 0$ となっています。

$\begin{align} \langle u, g \rangle = L(f) ||g|| - L\left( \frac{g}{||g||} \right) \langle f, g \rangle \end{align} \tag{26}$

これが0になるので、式(24) の $h_L$ が次の形で得られます。

$\begin{align} h_L = L\left( \frac{g}{||g||} \right) \frac{g}{||g||} \end{align} \tag{27}$

　一意性については、もし $L(f) = \langle f, h_L \rangle = \langle f, h'_L \rangle$ なら、全ての $f \in \mathcal{H}$ に対して $\langle f, h_L - h'_L \rangle=0$ が成立します。このような場合、最初の仮定より $h_L - h'_L = 0$ となります。

*1:参考文献：L. Rosasco , “Reproducing Kernel Hilbert Spaces,” http://www.mit.edu/~9.520/fall14/slides/class03/class03_rkhsPart1.pdf , 2014.

*2:参考文献：田中勝，「エントロピーの幾何学」，コロナ社，2019年5月（初版）

2020-08-09

シンプレクティック形式と数値積分

　ハミルトン力学の位置座標と運動量座標からは、シンプレクティック形式を自然に定義することができます。そうすると、例えば、シンプレクティック形式を保存する写像によって、正準変換やシンプレクティック積分による数値積分を記述することができます。

準備
ハミルトンベクトル場
シンプレクティック形式
正準変換
シンプレクティック数値積分法
高次のシンプレクティック数値積分法と近似精度

準備

　n次元ベクトル場 $X$ を次のように書きます。 $x^i$ は $i$ 番めの座標です。 $X^i$ は $i$ 番めの成分です。一般には $\frac{\partial}{\partial x^i}$ は基底を表す単なる記号であり、偏微分ではありません。しかし、変換則が分かりやすいので、この形で書かれます。

$\begin{align} X=\sum_{i=1}^{n} X^i \frac{\partial}{\partial x^i} \end{align} \tag{1}$

　特に、多様体 $M$ 上では、各点ごとに接ベクトル空間を定義できます。点 $p \in M$ での接ベクトル空間を、 $T_p (M)$ と書きます。

　また、微分 $k$ 形式は一般に次のように書けます。 $f_{l_1, \cdots, l_k}(x_1, \cdots, x_n)$ は関数です。

$\begin{align} u=\sum_{l_1< \cdots < l_k} {f_{l_1, \cdots, l_k}(x_1, \cdots, x_n) dx^{l_1} \wedge \cdots \wedge dx^{l_k}} \end{align} \tag{2}$

　関数 $f$ の微分 $df$ は、ベクトル場 $X$ を実数に移す写像 $df: T_p(X)\rightarrow \mathbb{R}$ とみなせます。微分形式の微分が次のようになるので、外積の反対称性から $d(du)=0$ が成り立ちます。

$\begin{align} du=\sum_{i=1}^n \sum_{l_1< \cdots < l_k} {\frac{\partial f_{l_1, \cdots, l_k}}{\partial x^i} dx^i \wedge dx^{l_1} \wedge \cdots \wedge dx^{l_k}} \end{align} \tag{3}$

　次に、開集合 $U$ から $V$ への可微分同相写像を $\Phi$ とします。 $(y^1, \cdots , y^n)=\Phi (x^1 ,\cdots , x^n)$ です。座標変換もこの写像で表せます。 $\Phi$ によるベクトル場の変換を次のように表します。

$\begin{align} \Phi_{*} X = \Phi_{*} \left( \sum_{i=1}^{n} X^i \frac{\partial}{\partial x^i} \right) = \sum_{j=1}^{n} \left( \sum_{i=1}^n {X^I \frac{\partial y^j}{\partial x^i}} \right) \frac{\partial}{\partial y^j} \end{align} \tag{4}$

　逆に、 $V$ から $U$ への微分形式の引き戻しを、次のように表します。ここでの $f$ は先程とは異なり、 $V$ 上の関数です。

$\begin{align} \Phi^* (u) = \Phi^* \left( \sum_{l_1< \cdots < l_k} {f_{l_1, \cdots, l_k} dy^{l_1} \wedge \cdots \wedge dy^{l_k}} \right) \\ = \sum_{l_1< \cdots < l_k} {(f_{l_1, \cdots, l_k} \circ \Phi) d\Phi^{l_1} \wedge \cdots \wedge d\Phi^{l_k}} \end{align} \tag{5}$

ハミルトンベクトル場

　ハミルトン力学系の話に移ります。ハミルトニアン $H$ によるハミルトンベクトル場は次のように定義されます。 $q^i$ は位置座標、 $p^i$ は運動量座標です。

$\begin{align} X_H = \sum_{i=1}^n {\frac{\partial H}{\partial p^i} \frac{\partial}{\partial q^i} - \frac{\partial H}{\partial q^i} \frac{\partial}{\partial p^i}} \end{align} \tag{6}$

　これは単なるベクトル場ですが、形式的にリウヴィル演算子 $i\mathcal{L}_H$ と同じ形になっています。ポアソン括弧を使うと $i\mathcal{L}_H=\{\cdot, H\}$ です。相空間上で定義される関数 $\gamma(q,p,t)$ について、次の式が成り立ちます。

$\begin{align} \frac{d\gamma}{dt} = i\mathcal{L}_H \gamma + \frac{\partial \gamma}{\partial t} \end{align} \tag{7}$

　例えば、 $\gamma =p,q$ の場合は、ハミルトン方程式になります。

　
　また、ハミルトン方程式の解の各々 $q(t), p(t)$ を積分曲線と言い、関数 $G(p, q)$ のうち、任意の積分曲線上で $\frac{dG}{dt}=0$ となるものを第一積分または保存量と言います。式(7) から、 $G(p,q)$ が第一積分である必要十分条件は $i\mathcal{L}_HG=0$ またはポアッソン括弧を使って $\{G,H \}=0$ です。

シンプレクティック形式

　次の微分2形式がシンプレクティック形式です。これは閉形式（ $d\omega = 0$ ）です。

$\begin{align} \omega = \sum_{I=1}^n {dp^i \wedge dq^i} \end{align} \tag{8}$

　写像 $\Phi: U \rightarrow V$ による変換 $(P, Q)=\Phi (p, q)$ を考えます。シンプレクティック形式が変換によって不変であることと、ハミルトンベクトル場が変換によって不変であることが等価であることを確認します。このことは、次のように表せます。 $\Omega$ は $(P,Q)$ 空間のシンプレクティック形式です。

$\begin{align} \omega = \Phi^* \Omega \longleftrightarrow \Phi_* X_{H(p,q)} = X_{H(P.Q)} \end{align} \tag{9}$

　道具として、微分形式に対する次のような線形作用素 $i_X$ を考えます。

$\begin{align} i_{\frac{\partial}{\partial x^{l_j}}} (dx^{l_1} \wedge \cdots \wedge dx^{l_n}) = (-1)^{j-1} dx^{l_1} \wedge \cdots \wedge dx^{l_{j-1}}\wedge dx^{l_{j+1}} \wedge \cdots \wedge dx^{l_n} \end{align} \tag{10}$
$\begin{align} i_{fX+gX'}(u)=f i_{X}(u) + g i_{X'}(u) \\ i_{X}(fu+gu') = f i_{X}(u) + g i_{X}(u') \end{align} \tag{11}$

$i_{\frac{\partial}{\partial x^i}}$ は、微分形式から $dx^i$ の項を先頭に移した後に、この項を除く操作を表しています。微分形式 $u$ が $dx^i$ を含まない場合には、 $i_{\frac{\partial}{\partial x^i}}(u)=0$ となります。

　一般に、 $\Phi : U \rightarrow V$ が定義された時、式(4)、(5)、(10)、(11) から、次の関係が成り立つことが分かります。 $u$ は $V$ 上の微分形式です。

$\begin{align} \Phi^* i_{\Phi_* X}(u) = i_{X}(\Phi^* u) \end{align} \tag{12}$

　 $\Phi$ が相空間 $U: (p,q)$ から $V: (P,Q)$ への写像であったことから、式(6)、式(8) を参照すると、

$\begin{align} i_{X_H}(\omega) = dH \end{align} \tag{13}$

が成り立ちます。式(12) によって、 $\omega = \Phi^* \Omega$ から、 $\Phi^* i_{\Phi_* X_{H(p,q)}}(\Omega) = i_{H(p,q)}(\omega)$ となります。 $\Phi^* dH(P,Q)$ を $dH(p,q)$ と同一視することで、式(13) と一意性から、 $\Phi_* X_{H(p,q)}=X_{H(P,Q)}$ となっていることが分かります。

正準変換

　変換 $\Phi :U \rightarrow V$ によってシンプレクティック形式が不変であるとき、この変換は正準変換です。シンプレクティック形式が不変ならハミルトンベクトル場も不変であり、変換後の座標でもハミルトニアンが存在します。

　正準変換によって結びつく二つのシンプレクティック形式が等しいことから、

$\begin{align} \sum_i ( dp^i \wedge dq^i -dP^i \wedge dQ^i)=0 \end{align} \tag{14}$

ですが、これを

$\begin{align} d\left( \sum_i ( p^i \wedge dq^i -P^i \wedge dQ^i)\right) =0 \end{align} \tag{15}$

と書くことで、次のような関数 $S(q,Q)$ が存在することが分かります（ $U$ が単連結な空間の場合）。

$\begin{align} dS = \sum_i ( p^i \wedge dq^i -P^i \wedge dQ^i) \end{align} \tag{16}$

　 $S(q,Q)$ を生成関数と呼びます。生成関数を用いて下記の式が得られます。

$\begin{align} p^i = \frac{\partial S}{\partial q^i} \\ P^i = -\frac{\partial S}{\partial Q^i} \end{align} \tag{17}$

　式(15) の書き方に応じて、4通りの生成関数 $S(q,Q), S(q,P), S(p,Q), S(p,P)$ を定義できます。生成関数を用いて、正準変換 $(P,Q)=\Phi(p,q)$ を構成することができます。

シンプレクティック数値積分法

　ハミルトン方程式の解を数値積分によって近似的に求めることを考えます。リウヴィル演算子を使って

$\begin{align} \frac{d}{dt}\left( \begin{array}{c} p \\ q \end{array} \right) = i\mathcal{L}\left( \begin{array}{c} p \\ q \end{array} \right) \end{align} \tag{18}$

ですから、 $\exp = \sum_n {\frac{1}{n!} \frac{\partial ^n}{\partial t^n}}$ と定義してテイラー展開と比較することで、時刻 $t$ から $t+\tau$ までの積分を次のような時間並進演算子として書くことができます。

$\begin{align} q(t+\tau) = \exp {(i\mathcal{L}_H\tau)}q(t) \end{align} \tag{19}$

$p$ についても同様です。 $\exp {(i\mathcal{L}_H\tau)}$ を $\tau$ が小さい場合の近似として具体的に書き下すことで、数値積分を行うことができます。

　ハミルトニアンを分解し、 $H(p,q)=K(p)+U(q)$ とします。リウヴィル演算子の定義から、 $i\mathcal{L}_H = i\mathcal{L}_K + i\mathcal{L}_U$ です。一粒子系なら $i\mathcal{L}_f = \frac{\partial f}{\partial p}\frac{d}{dq} - \frac{\partial f}{\partial q}\frac{d}{dp}$ ですから、

$\begin{align} i\mathcal{L}_K q = \frac{\partial K}{\partial p} \\ i\mathcal{L}_K p = 0 \\ i\mathcal{L}_U q = 0 \\ i\mathcal{L}_U p = -\frac{\partial U}{\partial q} \end{align} \tag{20}$

これらを使って、次の式が得られます。

$\begin{align} \exp{(i\mathcal{L}_K\tau)} q = q + \frac{\partial K}{\partial p} \tau \\ \exp{(i\mathcal{L}_K\tau)} p = p \\ \exp{(i\mathcal{L}_U\tau)} q = q \\ \exp{(i\mathcal{L}_U\tau)} p = p - \frac{\partial U}{\partial q} \tau \end{align} \tag{21}$

式(21) から、 $\exp{(i\mathcal{L}_K\tau)}$ や $\exp{(i\mathcal{L}_U\tau)}$ による変換がシンプレクティック形式 $dp \wedge dq$ を不変に保つことはすぐに分かります。これらをシンプレクティック写像と言います。シンプレクティック写像の積（合成写像）もシンプレクティック写像です。そのため、次のように時間並進演算子を近似することでシンプレクティック写像が得られます。

$\begin{align} \exp{(i\mathcal{L}_H\tau)} \approx \exp{(i\mathcal{L}_K\tau)}\exp{(i\mathcal{L}_U\tau)} \end{align} \tag{22}$

　式(22) に従って位置と運動量を更新するのが、一次のシンプレクティック数値積分法です。式(22) 右辺の二つの項は逆に並べても構いません。具体的なアルゴリズムは、式(21) を用いて次のように簡単に実行できます。

$\begin{align} q(t+\tau)=q(t) + \frac{\partial K(p(t))}{\partial p} \tau \\ p(t+\tau)=p(t) - \frac{\partial U(q(t+\tau))}{\partial q} \tau \end{align} \tag{23}$

このアルゴリズムは、最も単純な数値積分法であるオイラー法と同じ程度に簡単ですが、実際のハミルトニアン $H(p,q)$ に近い近似的なハミルトニアン $\tilde{H}(p,q)$ に厳密に従う時間発展を記述しています。

　式(22) の右辺から決まる下記の $\tilde{H}$ が、一次のシンプレクティック数値積分法において厳密に保存される近似的なハミルトニアンです。

$\begin{align} \exp{(i\mathcal{L}_\tilde{H}\tau)} = \exp{(i\mathcal{L}_K\tau)}\exp{(i\mathcal{L}_U\tau)} \end{align} \tag{23}$

　真のハミルトニアン $H$ と近似ハミルトニアン $\tilde{H}$ との誤差は、時間幅 $\tau$ に対して、シンプレクティック法の次数程度までしか拡大しません。従って、シンプレクティック法で数値積分を行う限り、エネルギー誤差が発散することはありません。これが、シンプレクティック形式を不変に保たないアルゴリズムと比較した場合の大きな利点です。

高次のシンプレクティック数値積分法と近似精度

　一次のシンプレクティック数値積分法は式(22) の近似に依拠しますが、式(22) の両辺は一致しません。ベーカー・キャンベル・ハウスドルフの公式を使って、 $\tau$ の次数に対する両辺の差を評価することができます。

$\begin{align} \exp{(A)}\exp{(B)}=\exp{\left( A+B+\frac{1}{2}[A,B] + \cdots \right) } \end{align} \tag{24}$

例えば、式(22) は

$\begin{align} \exp{(i\mathcal{L}_K\tau)}\exp{(i\mathcal{L}_U\tau)}=\exp{\left( i\mathcal{L}_H\tau+ O(\tau^2) \right) } \end{align} \tag{25}$

なので、一次精度です。

　より一般に、シンプレクティック形式を不変に保つように $\exp{( i\mathcal{L}_H\tau)}$ を近似する方法として、次の形が考えられます。

$\begin{align} \exp{(i\mathcal{L}_H\tau)} \approx \prod_j {\exp{(c_j (i\mathcal{L}_K\tau))}\exp{(d_j (i\mathcal{L}_U\tau))}} \end{align} \tag{26}$

ただし、一次近似を成立させるために $\sum_j {c_j} = 1, \sum_j {d_j} = 1$ である必要があります。

　例えば、次の形で近似すると、二次のシンプレクティック数値積分法（リープ・フロッグ法）が得られます。

$\begin{align} \exp{(i\mathcal{L}_H\tau)} \approx \exp{(\frac{i\mathcal{L}_K\tau}{2})}\exp{(i\mathcal{L}_U\tau)}\exp{(\frac{i\mathcal{L}_K\tau}{2})} \end{align} \tag{27}$

　式(24) で評価してみると、ちょうど $[i\mathcal{L}_K, i\mathcal{L}_U] \tau^2$ の項が打ち消しあうことで、 $\exp{(i\mathcal{L}_K\tau)}\exp{(i\mathcal{L}_U\tau)} = \exp{(i\mathcal{L}_H\tau) + O(\tau^3)}$ となり、二次精度の近似になっていることが分かります。さらに、式(27) は演算子として右から作用させても左から作用させても同じ、対称な形をしており、時間反転対称性を満たしています。

　さらに高次のシンプレクティック数値積分法も発見されています。例えば四次のシンプレクティック法は下記のように与えられています（Forest, 1989）。 $\alpha = 1 - 2^{\frac{1}{3}}$ として、

$\begin{align} \exp{(i\mathcal{L}_H\tau)} \approx \exp{\left( \frac{\tau}{2(1+\alpha)}i\mathcal{L}_K \right) } \exp{\left( \frac{\tau}{1+\alpha}i\mathcal{L}_U\right) } \exp{\left( \frac{\alpha \tau}{2(1+\alpha)}i\mathcal{L}_K\right) } \\ \exp{\left( \frac{(\alpha -1)\tau}{1+\alpha}i\mathcal{L}_U\right) } \exp{\left( \frac{\alpha \tau}{2(1+\alpha)}i\mathcal{L}_K\tau\right) } \exp{\left( \frac{\tau}{1+\alpha}i\mathcal{L}_U\right) } \exp{\left( \frac{\tau}{2(1+\alpha)}i\mathcal{L}_K\right) } \end{align} \tag{27}$

　その他、シンプレクティック積分については、散逸系のような非ハミルトン系へ適用するための Split 法や、Partitioned-ルンゲクッタ法としての表現といった話題があります。

2020-08-01

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

　ある確率密度関数 $f(z)$ が得られたけれど、その平均値や標準偏差、最大値などを解析的に計算できないことがよくあります。数値積分すれば良さそうですが、今度は数値積分の方法が問題になります。統計量や積分を直接計算する代わりに、もし、確率分布 $f(z)$ に従うサンプル $z_1 , z_2 , \cdots , z_N$ が得られれば、ある量 $g(z)$ の統計平均は $\langle g(z) \rangle = \int{g(z)f(z)dz}=\frac{1}{N} \sum_n {g(z_n)}$ と近似計算できます。

　最も単純なサンプリング法は、モンテカルロ法の一種である棄却サンプリングです。棄却サンプリングのうち、提案分布を一様分布としたものが、最も簡単なモンテカルロ法です。これは、確率変数 $z$ を一様分布からサンプルし、得られた $z$ を確率 $f(z)/M$ で採用し、確率 $1-f(z)/M$ で棄却するものです（M は $f(z)$ の最大値以上の値）。しかし、この方法では、棄却が非常に多くなってしまう問題や、適切な提案分布を選ぶのが難しいという問題があります。

　マルコフ連鎖モンテカルロ法は、提案分布の細かい構造を決めることで、上記の問題を解決します。そのマルコフ連鎖モンテカルロ法のために適切な提案分布を作る方法の一つが、メトロポリス・ヘイスティングス法です。高次元分布からサンプルする場合は、ギブスサンプリングが便利です。また、メトロポリス・ヘイスティングス法を効率的に構成する方法の一つが、ハミルトニアン・モンテカルロ法です。

マルコフ連鎖モンテカルロ法
詳細釣り合い条件
メトロポリス・ヘイスティングス法（MHアルゴリズム）
- MHアルゴリズムと詳細釣り合い条件との関係
- ギブスサンプリング
もう少し詳しい議論
ハミルトニアン・モンテカルロ法

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

　以降、確率分布を離散分布とし、 $p(z)$ と表します。 $z$ は $1〜k$ までの値をとる $k$ 状態変数とします。 $\int{g(z)f(z)dz} = \sum^{k}_{z=1} {g(z)p(z)}$ という対応関係があります。

　ステップ $t$ に対して、変数が確率的に変化することを考えます。ステップ $t$ における変数の値が、それ以前の $m$ ステップ、すなわちステップ $t-m〜t-1$ における変数の値に依存するとき、これを $m$ 階の離散時間マルコフ過程と呼びます。状態変数が離散的な場合はマルコフ連鎖とも言います。

　特に1階マルコフ過程では、ステップ $t$ で変数が $z'$ であったときに、ステップ $t+1$ で $z$ に遷移する確率 $q(z|z')$ が定義されます。 $q(z|z')$ は、マルコフ核とも呼ばれます。ステップ $t$ における確率分布を $q_t(z)$ と書くと、次の関係が成り立ちます。

$\begin{align} q_{t+1}(z)=\sum^{K}_{z'=1} q(z|z')q_t(z') \end{align} \tag{1}$

　一般に、ステップ $t$ での分布 $q_t(z)$ とステップ $t+1$ での分布 $q_{t+1}(z)$ とは異なる分布になりますが、ステップを進めることで特定の分布に収束し、それ以上変化しなくなることがあります（詳細は後述）。この分布を定常分布と呼び、ここでは $q_{\infty}(z)$ と書きます。すなわち、

$\begin{align} q_{\infty}(z)=\sum^{K}_{z'=1} q(z|z')q_{\infty}(z') \end{align} \tag{2}$

となります。

　もし、当初の目的としていた確率分布 $p(z)$ が定常分布 $q_{\infty}(z)$ に一致するならば、適当な初期値 $z_1$ からはじめて、定常分布に収束するまで $z_1 , z_2 , \cdots$ とサンプリングを繰り返すことで、最終的に $p(z)$ に従う多数のサンプルを得ることができます。もちろん、定常分布に収束する前の前半のサンプルは捨てなければなりません。また、サンプル間の相関が小さくなるように、採用するサンプル間のステップ幅を離すことが必要です。

　これが、マルコフ連鎖モンテカルロ法（MCMC）です。問題は、どのようにして $p(z)$ を定常分布に持つマルコフ過程を準備するかです。

詳細釣り合い条件

　次の関係が成立しているとき、マルコフ核 $q(z|z')$ は分布 $p(z)$ に対して詳細釣り合い条件を満たしていると言います。

$\begin{align} q(z'|z)p(z) = q(z|z')p(z') \end{align} \tag{3}$

　両辺それぞれで $z'$ に関して和をとると、

$\begin{align} p(z) = \sum^{K}_{z'=1}q(z|z')p(z') \end{align} \tag{4}$

となります。これは、（詳細でない）釣り合いの式です。式(2) と比較すると、 $p(z)$ はこのマルコフ過程の定常分布になっていることがわかります。定常分布に対して、（詳細でない）釣り合いの式(2),(4) は必ず成り立ちますが、詳細釣り合いの式(3) は、可逆過程に対応する特殊な条件であり、常に成り立つとは限りません。

　詳細釣り合い条件が満たされていれば、確率分布 $q_t(z)$ は $t \rightarrow \infty$ で目的の分布 $p(z)$ に収束します（詳細は後述）。残る問題は、詳細釣り合い条件を満たすマルコフ核をどうやって求めるかです。

メトロポリス・ヘイスティングス法（MHアルゴリズム）

　MHアルゴリズムは、適当な提案分布 $\tilde{q}(z|z')$ からはじめて、詳細釣り合い条件を満たすマルコフ核 $q(z|z')$ を作ることで、目的の分布 $p(z)$ からのサンプリングを行うアルゴリズムです。

　アルゴリズムは次のようなものです。

適当な事前分布に従って、初期値 $z_1$ をサンプルする。

ステップ $t$ において値 $z'$ をサンプルした時、次のステップでのサンプル候補 $z$ を、適当な提案分布 $\tilde{q}(z|z')$ に従ってサンプルする。

次の値を計算する。 $\rho(z, z') = \frac{\tilde{q}(z'|z) p(z)}{\tilde{q}(z|z') p(z')} \tag{5}$

$\rho(z,z') \geq 1$ なら、サンプル候補 $z$ を次のステップ $t+1$ でのサンプル値とする。 $\rho(z,z') < 1$ なら、確率 $\rho$ で $z$ を次のステップ $t+1$ でのサンプル値とし、確率 $1 - \rho(z,z')$ で $z'$ （遷移前と同じ値）を次のステップ $t+1$ でのサンプル値とする。

　このアルゴリズムによって生成されたマルコフ連鎖は、目的の分布 $p(x)$ に収束し、安定分布となります。提案分布 $\tilde{q}(z|z')$ は何でも良いのですが、 $z'$ を中心としたガウス分布等の単峰性の分布や矩形分布が用いられます。中心からの幅が狭すぎると、隣接ステップの相関が大きくなり、無駄なサンプルが増えます。一方で中心からの幅が広すぎると、マルコフ性の利点が失われてしまいます。提案分布の形状（単峰性か多峰性か等）が目的の分布に類似しているほど、マルコフ性の利点が増します。

MHアルゴリズムと詳細釣り合い条件との関係

　それでは、MHアルゴリズムが詳細釣り合い条件を満たすマルコフ核を作ることを見てみましょう。アルゴリズムをマルコフ過程として書き下すと、 $z \neq z'$ の場合は

$\begin{align} q_{t+1}(z)=\sum^K_{z'=1} \min{(\rho(z, z'), 1)}\tilde{q}(z|z')q_{t}(z) \end{align} \tag{6}$

となります。式(1) と比較することで、 $z \neq z'$ の場合のマルコフ核は

$\begin{align} q(z|z') = \min{(\rho(z, z'), 1)}\tilde{q}(z|z') \end{align} \tag{7}$

であることが分かります。 $q(z|z')$ は確率分布ですから、 $z = z'$ の場合は

$\begin{align} q(z'|z') = 1 - \sum_{z \neq z'} q(z|z') \end{align} \tag{8}$

となります。

　式(7) のマルコフ核の性質を調べます。式(5) から、 $\rho(z, z')=\frac{1}{\rho(z', z)}$ です。従って、 $\rho(z, z')$ と $\rho(z, z')$ の一方は 1 以上で、他方は 1 以下です。例えば $\rho(z, z') \leq 1$ ならば

$\begin{align} q(z|z') &=& \rho(z, z') \tilde{q}(z|z') \\ q(z'|z) &=& \tilde{q}(z'|z) \end{align} \tag{9}$

という関係が成り立っており、 $\rho(z,z')$ に式(5) を代入すると、確かに詳細釣り合い条件 $q(z'|z)p(z) = q(z|z')p(z')$ が成立しています。

ギブスサンプリング

　詳細釣り合いに依拠するMHアルゴリズムは汎用的で、目的の分布に確実に収束する方法です。しかし、変数 $z$ の次元が大きければ、マルコフ核の底空間も巨大になり、サンプリングが追いつかなくなります。ギブスサンプリングはMHアルゴリズムほど汎用的ではありませんが、高次元空間においては非常に強力なサンプリング法です。

　ギブスサンプリングは、目的の分布から各変数ごとの条件付き分布を求め、各変数に対して逐次的にサンプリングを繰り返し、収束を狙うものです。まず、 $z=(x_1, x_2, \cdots, x_d)$ として、目的の確率分布 $p(x_1, x_2, \cdots, x_d)$ から、次のように変数ごとの条件付き分布を求めます。

$\begin{align} p(x_1 | x_2, x_3, \cdots, x_d), \ \ \ p(x_2 | x_1, x_3, \cdots, x_d), \ \ \ \cdots , \ \ \ p(x_d| x_1, x_2, \cdots, x_{d-1}) \end{align} \tag{10}$

　そして、各変数を順番にサンプリングしながら更新し、収束するまで繰り返します。

もう少し詳しい議論

　これまでの説明の中には、いくつかごまかしがあります。下記の項目について、もう少し詳しく見てみましょう。

安定分布は、遷移確率（マルコフ核）だけによって表現できるのか
詳細釣り合い条件が成り立てば、なぜ定常分布に収束するのか
詳細釣り合い条件が成り立たないのはどのような場合か

マルコフ過程と行列方程式

　マルコフ核は、 $K_{zz'}=q(z|z')$ という $k\times k$ 行列 $K$ で表せます。確率変数の分布 $q_t(z)$ も同様に $k$ 次元ベクトル $q_t$ で表せます。式(1) は

$\begin{align} q_{t+1}=Kq_t \end{align} \tag{11}$

と表せます。定常分布を $K$ 次元ベクトル $p$ と書くと、定常分布を表す行列方程式

$\begin{align} p=Kp \end{align} \tag{12}$

が得られます。これが自明でない解を持つので、 $I-K$ の各行（または列）は線形独立ではありません。実際、下記の拘束条件が含まれています。

$\begin{align} \sum_z K_{zz'} &=& 1 \\ \sum_z p_z &=& 1 \end{align} \tag{13}$

　そこで、 $p_k = 1 - p_1 - p_2 - \cdots - p_{k-1}$ を用いて、式(12) から $p_k$ を消去すると、次の行列方程式が得られます。

$\begin{align} p'=K'p'+b' \end{align} \tag{14}$

　ただし、 $K'$ は $(k-1)\times (k-1)$ 行列、 $p', b'$ はいずれも $k-1$ 次元ベクトルで、

$\begin{align} p'_z &=& p_z \\ K'_{zz'} &=& K_{zz'} - K_{zk} \\ b'_z &=& K_{zk} \end{align} \tag{15}$

を満たします。式(14) の各行は独立に与えられるので、 $I-K'$ が逆行列を持てば、

$\begin{align} p'=(I-K')^{-1}b' \end{align} \tag{16}$

によって、遷移確率から安定分布を求めることができます。

詳細釣り合い条件と収束性

　式(1) で表されるマルコフ過程において、式(3) の詳細釣り合い条件が満たされていれば、 $q_t(z)$ は $p(z)$ に収束することを示します。そのために、 $q_t(z)$ と $p(z)$ の間の、 $p(z)$ を基準としたカルバック・ライブラー・ダイバージェンス $D( p||q_t )$ が、ステップ毎に単調に減少することを確認します。

$\begin{align} D(p||q_{t+1}) &=& -H(p) - \sum_z p(z) \log{q_{t+1}(z)} \\ &=& -H(p) - \sum_z p(z) \log{\sum_{z'}{q(z|z')q_t(z')}} \\ &=& -H(p) - \sum_z p(z) \log{\sum_{z'}{\frac{q(z'|z)p(z)}{p(z')}q_t(z')}} \\ &=& - \sum_z p(z) \log{\sum_{z'}{\frac{q(z'|z)q_t(z')}{p(z')}}} \\ &\leq & - \sum_z p(z) \sum_{z'}{q(z'|z)\log{\frac{q_t(z')}{p(z')}}} \\ &=& - \sum_{z, z'} p(z')q(z|z')\log{\frac{q_t(z')}{p(z')}} \\ &=& - \sum_{z'} p(z')\log{\frac{q_t(z')}{p(z')}} \\ &=& D(p||q_t) \end{align} \tag{17}$

となり、確かに $D(p||q_{t+1}) \leq D(p||q_t)$ となっています。 $q_t(z)$ が $p(z)$ に一致するまで、 $D(p||q_t)$ が減少し続けることが分かります。ただし、より厳密には、減少速度を評価する必要があります。 $H$ はエントロピーを表します。なお、上記の式変形において、1行目から2行目へは式(1) を用い、2行目から3行目へは詳細釣り合いの式(3) を用いました。4行目と5行目の間の不等号は、イェンゼンの不等式によります。5行目から6行目へは再び詳細釣り合いの式(3) を用いました。

詳細釣り合い条件と自由度

　マルコフ核は、式(13) の拘束条件の範囲で自由な正値を取ることができます。式(13) の拘束条件は式(8) に反映されていますから、一般にマルコフ核の自由度は、 $z=z'$ を除き、 $k^2 - k$ 個です。

　一方で、詳細釣り合い条件が満たされている場合、式(3) より $q(z,z')=\frac{q(z',z)p(z)}{p(z')}$ なので、 $q(z',z)$ と $p(z)$ が決まれば $q(z', z)$ も決まります。従って、この場合の自由度はマルコフ核のうち $z < z'$ であるものと $p(z)$ のうち一つを除いたものの個数の和であり、 $\frac{k(k+1)}{2}-1$ 個です。

　自由度が詳細釣り合い条件の有無によらないのは、 $k=1, 2$ の場合のみです。 $k=2$ の場合は、 $q(2|1)p(1)=q(1|2)p(2)$ が必ず成り立ちます。 $k \geq 3$ では、詳細釣り合い条件を満たさないマルコフ過程が存在します。

ハミルトニアン・モンテカルロ法

　メトロポリス・ヘイスティングス法（MHアルゴリズム）では、遷移確率（マルコフ核）の提案分布 $\tilde{q}(z|z')$ をどのように決めるかが問題です。 $\tilde{q}(z|z')$ の決め方の一つが、ハミルトニアンモンテカルロ法です。

　そもそも、MHアルゴリズムの目的は、サンプリングを効率よく行うために、モンテカルロ法での棄却を抑制することでした。MHアルゴリズムにおいて、式(5) で表される $\rho(z, z')$ を1 に近づけることができれば、棄却率は低くなります。

　ハミルトニアンモンテカルロ法では、ハミルトン方程式に従って、時刻 $t$ での位置座標 $z'$ から時刻 $t+1$ での位置座標 $z$ への変化をシミュレートします。通常、ハミルトン方程式による時間発展は確率的ではなく、確定的な変化です。そのため、運動量の初期値をサンプリングすることによって、確率的な遷移を表現します。このことから分かるように、マルコフ核 $\tilde{q}(z|z')$ は明示的な関数形では与えられません。運動量の初期値の分布が決まっていて、そこから得られた運動量のサンプルと位置の初期値 $z'$ を、ハミルトン方程式に従う時間発展という複雑な関数に通すことで、遷移後のサンプル $z$ が確率的に得られます。

　ハミルトニアンモンテカルロ法のポイントは、ハミルトニアンが次の二つの性質を持っていることです。一つめは、ハミルトニアンによる時間発展は可逆であるという性質です。二つめは、ハミルトニアンに従う時間発展において、エネルギーの値（ハミルトニアンに位置と運動量を代入した値）は不変であるという性質です。特に、二つめの性質によって、式(5) の $\rho (z,z')$ の値が 1 に近くなります。

　ハミルトニアンモンテカルロ法のアルゴリズムを説明します。以下、 $z$ を位置座標として扱い、運動量を $P$ とします（確率分布の記号ではないことに注意してください）。まず、ハミルトニアン $H(k,z)=K(P)+U(z)$ は次のように定義します。

$\begin{align} K(P)=\frac{P^2}{2} \\ U(z)=-\log{p(z)} \tag{18} \end{align}$

　 $p(z)$ は、目的の確率分布です。運動エネルギーの項は質量 1 の粒子に対応し、位置エネルギーの項は逆温度 1 に対応しています。

　運動量 $P$ の初期値は、平均 0、標準偏差 1 のガウス分布からサンプリングします。なぜなら、 $\exp{(-H(P,z))}=p(z)\exp{(-\frac{P^2}{2})}$ だからです。従って、運動エネルギーの項が質量 $m$ の粒子を表すように決めることは、標準偏差 $\sqrt{m^{-1}}$ のガウス分布を選ぶことに対応します。

　シミュレーションは、ハミルトン方程式を数値積分することで行います。ハミルトン方程式は以下の通りです。

$\begin{align} \frac{dz}{dt} = \frac{\partial H}{\partial P} \\ \frac{dP}{dt} = -\frac{\partial H}{\partial z} \end{align} \tag{18}$

これから、エネルギーが時間不変であることは、次のようにして分かります。

$\begin{align} \frac{d H(P, z)}{dt}=\frac{dP}{dt} \frac{\partial H}{\partial P} + \frac{dz}{dt} \frac{\partial H}{\partial z}=0 \end{align} \tag{19}$

　可逆性については、 $t^*=-t, P^*=-P, z^*=z$ と置き換えても式(18) が成立することから分かります。

　式(11) を数値積分する方法としては、ルンゲ・クッタ法が最も汎用的ですが、ハミルトン力学系ではシンプレクティック法が良い性質を示します。ここでは、可逆性を持つ二次のシンプレクティック法を記載します。 $\tau$ を時間ラベルとして、十分小さい値 $\epsilon$ に対し、以下のように時間発展を行います。

$\begin{align} z(\tau + \frac{\epsilon}{2})=z(\tau)+\frac{\partial K (P(\tau))}{\partial P} \frac{\epsilon}{2}\\ P(\tau + \epsilon)=P(\tau)-\frac{\partial U(z(\tau + \frac{\epsilon}{2}))}{\partial z}\epsilon \\ z(\tau + \epsilon)=z(\tau + \frac{\epsilon}{2})+\frac{\partial K (P(\tau + \epsilon))}{\partial P} \frac{\epsilon}{2} \end{align} \tag{20}$

　この計算を $L$ 回繰り返すことで、時刻 $\epsilon L$ だけ先の $z$ が得られます。この値を、提案分布による遷移先のサンプルとします。 $L=1$ でも構いません。

　結局、ハミルトニアンモンテカルロ法のアルゴリズムは次のようになります。

サンプルの初期値 $z_1$ を決める。

サンプル $z'$ が得られたとき、次の値をサンプリングするための運動量の初期値 $P(0)$ を、平均 0、標準偏差 1 のガウス分布からサンプルする。

位置座標の初期値 $z(0)=z'$ と運動量座標の初期値 $P(0)$ から、式(20) の更新式を $L$ 回繰り返し適用する。

得られた $z(\epsilon L)$ を、サンプル $z$ の候補とし、 $\rho(z,z')=\exp{[-\{H(P(\epsilon L), z(\epsilon L))-H(P(0), z(0))\}]}$ を計算する。

確率 $\min{(\rho (z,z'), 1)}$ でサンプルを採択し、 $z=z(\epsilon L)$ とする。確率 $1-\min{(\rho (z,z'), 1)}$ でサンプルを棄却し、 $z=z'$ とする。

　アルゴリズムの中で、 $\rho(z,z')=\exp{[-\{H(P(\epsilon L), z(\epsilon L))-H(P(0), z(0))\}]}$ が非自明に見えるかもしれません。 $\rho(z,z')$ の定義は、式(5) で与えられていました。そして、最終的な遷移確率が運動量の初期値の確率分布に依存することから、 $\tilde{q}(z|z')=\frac{1}{\sqrt{2\pi}}\exp{(-\frac{P(0)^2}{2})}$ です。一方で、ハミルトン方程式の可逆性より、 $\tilde{q}(z'|z)=\frac{1}{\sqrt{2\pi}}\exp{(-\frac{P(\epsilon L)^2}{2})}$ も成り立ちます。

2019-07-24

揺らぎの定理

揺らぎの定理は，ある任意の過程（あるいは軌道） $\Gamma \colon \{ x(t_1)=x_1, x(t_2)=x_2 \cdots x(t_L)=x_L \}$ とその逆過程 $\Gamma^* \colon \{ x(t_1)=x_L^*, x(t_2)=x_{L-1}^* \cdots x(t_L)=x_1^* \}$ に対して，それらの生起確率の間の関係が $\begin{align} \frac{p(\Gamma)}{p(\Gamma^*)}=e^{\frac{\Delta S}{k_B}} \end{align} \tag{1}$ であるということを主張するものです．ここで， $\Delta S$ は過程 $\Gamma$ におけるエントロピーの増加量です．*1

揺らぎの定理から，例えば $e^{-\frac{\Delta S}{k}}$ の平均値 $\langle e^{-\frac{\Delta S}{k}}\rangle$ が $\begin{align} \langle e^{-\frac{\Delta S}{k}}\rangle = \int e^{-\frac{\Delta S}{k}} p(\Gamma) d\Gamma = \int p(\Gamma) d\Gamma = 1 \end{align} \tag{2}$ であることが分かります．ここで， $d\Gamma$ は全軌道空間上に取った確率測度で，全空間での積分を $1$ としています．Jensen の不等式より $\langle e^{-\frac{\Delta S}{k_B}} \rangle \geqq e^{\langle -\frac{\Delta S}{k_B}\rangle }$ ですから，(2)は $\langle \Delta S \rangle \geqq 0$ となり，熱力学第二法則が自然に導かれます．
(2)から分かるように，熱力学第二法則は無数の過程の統計平均について述べています．それに対して(1)は一つ一つの微細な過程間の定量的な関係について述べる，極めて根本的な表現です．

揺らぎの定理は以下のように考えることができます．

時刻 $t_1,t_2\cdots$ の間隔を限りなく小さくできるものとします．
$\begin{align} \frac{p(\Gamma)}{p(\Gamma^*)}=\frac{p(x_1,x_2\cdots x_L )}{p( x_L^*,x_{L-1}^*\cdots x_1^* )}=\frac{p(x_L | x_{L-1})\cdots p(x_2 |x_1)p(x_1)}{p(x_L^*)p(x_{L-1}^* | x_L^*)\cdots p(x_1^* | x_2^*)}\end{align}\tag{3}$
と表しておきます． $p$ の中の変数は，時間順に並べています．

微小な状態変化過程とその逆過程との間に詳細釣り合いが成立すると仮定すると， $p(x_n, x_{n+1})=p(x_{n+1}^*, x_n^*)$ より $\begin{align} \frac{p(x_{n+1} | x_n)}{p(x_n^* | x_{n+1}^*)}=\frac{p(x_{n+1}^*)}{p(x_n)}=e^{-\frac{\Delta E_n}{kT}}=e^{\frac{\Delta S_n^{e}}{k}} \end{align} \tag{4}$ と表せます．ここで， $T$ は系が接している熱浴の温度で， $\Delta E_n$ は $x_n \rightarrow x_{n+1}$ における系のエネルギー変化です．従って $-\Delta E_n$ は逆に熱浴が受け取る熱量となり，熱浴のエントロピー変化 $\Delta S_n^e$ を導きます．ただし，ここでの議論が，局所平衡の前提に基づいていることには注意する必要があります．

一方， $\begin{align} \frac{p(x_1)}{p(x_L^*)}=e^{\frac{s_L - s_1}{k_B}}=e^{\frac{\Delta S^i}{k_B}} \end{align}\tag{5}$ が成り立ちます． $s_n$ は状態 $x_n$ の微視的エントロピーです。また， $\Delta S^i$ は系の微視的エントロピーの増分を表しています．

(4)，(5) を (3) に代入し，最終的な熱浴のエントロピー変化を $\Delta S^e = S_1^e+S_2^e+\cdots +S_L^e$ ，系と熱浴を合わせた全体系のエントロピー変化を $\Delta S = \Delta S^i + \Delta S^e$ と書くと，揺らぎの定理(1)が妥当であることが分かります．

Jarzynski 等式

揺らぎの定理(1)からは，系のヘルムホルツ自由エネルギー $F$ と系に加えられた仕事 $W$ との間で成り立つ，Jarzynski等式と呼ばれる次の関係を導くことができます．
\begin{align}
\langle e^{-\beta W} \rangle = e^{-\beta \Delta F} \tag{6}
\end{align}
Jensen の不等式から $\langle e^{-\beta W} \rangle \geqq e^{-\langle \beta W \rangle}$ が成り立ちますから，(6)は $\Delta F \leqq \langle W \rangle$ となります．系の自由エネルギー増加は加えられた仕事以下になり，逆に $\Delta F$ ， $W$ が負である場合を考えると，系から取り出せる仕事が自由エネルギー以下になることが分かります．これらも，熱力学第二法則を与えています．

さて，(6)を導出するために，系の最初の状態を $x_1$ ，最後の状態を $x_L$ とします．これらの間を結ぶ（拘束された）軌道を $\tilde{\Gamma}$ と書きます．
系の最初と最後の状態が確定していますから，この場合揺らぎの定理は
\begin{align}
\frac{p(x_L, \tilde{\Gamma} | x_1)}{p(x_1^*, \tilde{\Gamma}^* | x_L^*)}=e^{\frac{S^e}{k_B}} \tag{7}
\end{align}
と書けます．
状態 $x$ における系の内部エネルギー $U(x)$ ，系に与えられた仕事 $W$ とすると， $W=U(x_L)-U(x_1)+T\Delta S^e$ です．
これらから， $e^{-\beta W}$ の統計平均 $\langle e^{-\beta W} \rangle$ を計算してみます．
\begin{align}
\langle e^{-\beta W} \rangle &=&\int dx_1 dx_L d\tilde{\Gamma} e^{-\beta (U(x_L)-U(x_1)+T\Delta S^e)}p(x_L, \tilde{\Gamma} | x_1)p(x_1) \\
&=& \int dx_1 dx_L d\tilde{\Gamma} e^{-\beta (U(x_L)-U(x_1)+T\Delta S^e)}p(x_L, \tilde{\Gamma} | x_1) \cdot \frac{e^{-\beta U(x_1)}}{Z(x_1)} \\ &=& \int dx_1 dx_L d\tilde{\Gamma} \frac{e^{-\beta U(x_L)} }{Z(x_1)}p(x_1^*, \tilde{\Gamma}^* | x_L^*) \\
&=& \int dx_L \frac{e^{-\beta U(x_L)} }{Z(x_1)}=\frac{Z(x_L)}{Z(x_1)}=e^{-\beta \Delta F} \tag{8}
\end{align}
1行目から2行目に移るのに，最初の状態 $x_1$ が平衡状態であることを仮定しました．これは非自明な仮定です．2行目から3行目では(7)を用いています．最後に，自由エネルギーと分配関数との関係 $-\beta F = \log{Z}$ を用いて自由エネルギーの形に書き換えました．

※追記予定

参考文献

https://www.jps.or.jp/books/gakkaishi/2014/10/69-10trends1.pdf

*1: