2018年11月30日金曜日

Conditional Random Field

はじめに


 この文献をベースにしてCRFについてまとめる。

CRFモデル


 いま、入力系列を$x_{1:N}(=x_1,\cdots,x_N)$、出力系列を$z_{1:N}(=z_1,\cdots,z_N)$と置き、下のグラフィカルモデルを考える。
このとき、Linear Chain CRFは次の条件付き確率で定義される。 \begin{equation} p(z_{1:N}|x_{1:N})=\frac{1}{Z}\exp{\left( \sum_{n=1}^N\sum_{i=1}^F\lambda_i f_i\left(z_{n-1},z_n,x_{1:N},n\right) \right)} \end{equation} $f_i(\cdot)$は素性関数(feature function)と呼ばれ、あらかじめ定義しておく関数である。重み$\lambda_i$は学習により決定される。$Z$は分配関数と呼ばれ、次式で定義される。 \begin{equation} Z=\sum_{z_{1:N}}\exp{\left( \sum_{n=1}^N\sum_{i=1}^F\lambda_i f_i\left(z_{n-1},z_n,x_{1:N},n\right) \right)} \end{equation} 従って \begin{equation} \sum_{z_{1:N}}p(z_{1:N}|x_{1:N})=1 \end{equation} が成り立つ。

素性関数(Feature Function)


 指数関数の肩の部分を展開すると \begin{eqnarray} \sum_{n=1}^N\sum_{i=1}^F\lambda_i f_i\left(z_{n-1},z_n,x_{1:N},n\right) &=& \lambda_1\sum_{n=1}^N f_1\left(z_{n-1},z_n,x_{1:N},n\right)\\ &+& \lambda_2\sum_{n=1}^N f_2\left(z_{n-1},z_n,x_{1:N},n\right)\\ &+& \lambda_3\sum_{n=1}^N f_3\left(z_{n-1},z_n,x_{1:N},n\right)\\ &+& \cdots\nonumber\\ &+& \lambda_F\sum_{n=1}^N f_F\left(z_{n-1},z_n,x_{1:N},n\right)\nonumber\\ \end{eqnarray} となる。先に関数$f_i(\cdot)$はあらかじめ決めて置く量であると書いた。例えば以下のように与えることができる。 \begin{eqnarray} f_1\left(z_{n-1},z_n,x_{1:N},n\right)&=& \begin{cases} 1&\mbox{if $z_n$=PERSON and $x_n$=John}\\ 0&\mbox{otherwise} \end{cases}\\ f_2\left(z_{n-1},z_n,x_{1:N},n\right)&=& \begin{cases} 1&\mbox{if $z_n$=PERSON and $x_{n+1}$=said}\\ 0&\mbox{otherwise} \end{cases}\\ f_3\left(z_{n-1},z_n,x_{1:N},n\right)&=& \begin{cases} 1&\mbox{if $z_{n-1}$=OTHER and $z_{n}$=PERSON}\\ 0&\mbox{otherwise} \end{cases} \end{eqnarray} $f_1(\cdot)$は、同じ場所の入力値$x_n$と出力値$z_n$から決まる素性である。一方、$f_2(\cdot)$は、現在の出力値$z_n$と次の入力値$x_{n+1}$から決まる素性となっている。また、$f_3(\cdot)$は2つの隣接する出力値$z_{n-1},z_n$から決まる素性である。もちろん、$f_i(\cdot)$は2値である必要はなく、実数値を取る関数を与えることもできる。 いま \begin{eqnarray} F_1&=&\sum_{n=1}^N f_1\left(z_{n-1},z_n,x_{1:N},n\right) \\ F_2&=&\sum_{n=1}^N f_2\left(z_{n-1},z_n,x_{1:N},n\right) \\ F_3&=&\sum_{n=1}^N f_3\left(z_{n-1},z_n,x_{1:N},n\right) \end{eqnarray} などと置くと、指数型の部分は \begin{equation} \sum_{n=1}^N\sum_{i=1}^F\lambda_i f_i\left(z_{n-1},z_n,x_{1:N},n\right)=(\lambda_1,\cdots,\lambda_N) \cdot (F_1,\cdots,F_N)^T \end{equation} と書くことができる。$F=(F_1,\cdots,F_N)$を素性ベクトルと呼ぶ。

CRFの最適化


 先に書いたようにCRFのパラメータ$\lambda_i$は学習により決まる量である。最尤推定法を用いる。 \begin{equation} L=\left(\prod_{j=1}^M p(z_{1:N}^{(j)}|x_{1:N}^{(j)})\right)p(\lambda_{1:F}) \end{equation} ただし、$\lambda_i$についての事前確率 \begin{eqnarray} p(\lambda_{1:F})&=&\prod_{i=1}^F p(\lambda_{i})\\ p(\lambda_{i})&=&\mathcal{N}(\lambda_i|0,\sigma^2) \end{eqnarray} を導入した。両辺対数を取れば \begin{equation} \ln{L}=\sum_{j=1}^M \ln{p(z_{1:N}^{(j)}|x_{1:N}^{(j)})}-\frac{1}{2}\sum_i^F\frac{\lambda_i^2}{\sigma^2} \end{equation} $\lambda_k$で微分する。 \begin{equation} \frac{\partial L}{\partial\lambda_k} = \sum_j\sum_n f_k(z_{n-1}^{(j)},z_n^{(j)},x_{1:N}^{(j)},n)- \frac{\partial}{\partial\lambda_k}\sum_j\ln{Z^{(j)}}- \frac{\lambda_k}{\sigma^2} \end{equation} ここで \begin{eqnarray} \frac{\partial\ln{Z^{(j)}}}{\partial\lambda_k} &=& \frac{1}{Z^{(j)}}\frac{\partial Z^{(j)}}{\partial\lambda_k}\\ &=& \sum_n\sum_{z_{n-1},z_n}p(z_{n-1},z_{n}|x_{1:N})f_k(z_{n-1}^{(j)},z_n^{(j)},x_{1:N}^{(j)},n)\\ &=& \sum_n\mathbb{E}_{z_{n-1},z_n} \left[ f_k(z_{n-1}^{(j)},z_n^{(j)},x_{1:N}^{(j)},n) \right] \end{eqnarray} 従って \begin{equation} \frac{\partial L}{\partial\lambda_k} = \sum_j\sum_n f_k(z_{n-1}^{(j)},z_n^{(j)},x_{1:N}^{(j)},n)- \sum_j\sum_n\mathbb{E}_{z_{n-1},z_n} \left[ f_k(z_{n-1}^{(j)},z_n^{(j)},x_{1:N}^{(j)},n) \right] - \frac{\lambda_k}{\sigma^2} \end{equation} を得る。解析的には解けないので、例えば、勾配降下法を用いて近似解を計算することになる。 \begin{equation} \lambda_k\leftarrow\lambda_k-\eta \frac{\partial L}{\partial\lambda_k} \end{equation}

2018年11月18日日曜日

ディープニューラルネットワークとガウス過程の関係

はじめに


 ディープニューラルネットワークとガウス過程の関係をまとめる。先のページの続きである。

ディープニューラルネットワーク


 以下のネットワークを考える。
各量の定義は以下の通り。 \begin{eqnarray} z^{l-1},x^{l}&\in& \mathbb{R}^{N_l} \\ z^{l}&\in&\mathbb{R}^{N_{l+1}}\\ \end{eqnarray} \begin{eqnarray} x_j^l(x)&=&\phi\left(z_j^{l-1}(x)\right)\label{eq1}\\ z_i^l(x)&=&b_i^l+\sum_{j=1}^{N_l}W_{ij}^l x_j^l(x)\label{eq2} \end{eqnarray}

ガウス過程との関係


 $z^l_i$がガウス過程と等価となることを帰納法により証明する。

 $l=0,1$のとき、先のページでの議論から$z^0_i$と$z^1_i$はガウス過程に従う量である。 いま、$z^{l-1}_i$がガウス過程に従う量であるとする。このとき、式(\ref{eq1})から$x_j^l$と$x_{j^{\prime}}^l$は、$j\neq j^{\prime}$のとき、独立である。式(\ref{eq2})から$z_i^l$は独立同分布から生成される項の和となるから、$N_l\rightarrow\infty$の極限で、$z_i^l$はガウス分布に従う(中心極限定理)。以上により、$z^l_i$がガウス過程と等価となることが証明された。ただし、$N_{0}\rightarrow\infty,N_{1}\rightarrow\infty,\cdots,N_{l-1}\rightarrow\infty,N_{l}\rightarrow\infty$の順に極限を取る必要がある。

平均と共分散行列の導出


 $z^l_i$がガウス過程であることが示されたので、その平均と共分散行列を導く。平均は以下となる。 \begin{eqnarray} \mathbb{E}\left[z^l_i(x)\right] &=& \mathbb{E}\left[b_i^l\right]+\sum_{j=1}^{N_l}\mathbb{E}\left[W_{ij}^l x_j^l(x)\right]\\ &=& \mathbb{E}\left[b_i^l\right]+\sum_{j=1}^{N_l}\mathbb{E}_l\left[W_{ij}^l\right]\mathbb{E}_{l-1}\left[ x_j^l(x)\right] \end{eqnarray} ここで、$\mathbb{E}\left[b_i^l\right]=0,\mathbb{E}_l\left[W_{ij}^l\right]=0$であるから、$\mathbb{E}\left[z^l_i\right]=0$を得る。共分散行列は次式となる。 \begin{eqnarray} \mathbb{E}\left[z^l_i(x)z^l_i(x^{\prime})\right] &=& \mathbb{E}_l\left[b_i^l b_i^l\right]+ \sum_{j=1}^{N_l}\sum_{k=1}^{N_l} \mathbb{E}_l\left[W_{ij}^l W_{ik}^l\right] \mathbb{E}_{l-1}\left[x_j^l(x) x_k^l(x^{\prime})\right] \end{eqnarray} ここで \begin{eqnarray} \mathbb{E}_l\left[b_i^l b_i^l\right]&=&\sigma_b^2 \\ \mathbb{E}_l\left[W_{ij}^l W_{ik}^l\right]&=&\frac{\sigma_w^2}{N_l}\delta_{jk} \end{eqnarray} であるから \begin{eqnarray} \mathbb{E}\left[z^l_i(x)z^l_i(x^{\prime})\right] &=& \sigma_b^2+ \frac{\sigma_w^2}{N_l} \sum_{j=1}^{N_l} \mathbb{E}_{l-1}\left[x_j^l(x) x_j^l(x^{\prime})\right]\\ &=& \sigma_b^2+ \frac{\sigma_w^2}{N_l} \mathbb{E}_{l-1}\left[\sum_{j=1}^{N_l} x_j^l(x) x_j^l(x^{\prime})\right]\\ &=& \sigma_b^2+ \frac{\sigma_w^2}{N_l} \mathbb{E}_{l-1}\left[x^l(x)^T x^l(x^{\prime})\right]\\ &=& \sigma_b^2+ \frac{\sigma_w^2}{N_l} \mathbb{E}_{l-1}\left[ \phi\left(z^{l-1}\left(x\right)\right)^T \phi\left(z^{l-1}\left(x^{\prime}\right)\right) \right]\label{eq3}\\ &\equiv&K^l(x,x^{\prime}) \end{eqnarray} となる。$N$個の入力ベクトル$\left(x^{(1)},\cdots,x^{(N)}\right)$に対応する次の量 \begin{equation} {\bf z}_i^l=\left(z_i^l\left(x^{(1)}\right),\cdots,z_i^l\left(x^{(N)}\right)\right) \end{equation} は、ガウス分布$\mathcal{N}({\bf 0},K^l)$に従う。ここで、$K^l$の$(n,m)$成分は次式で与えられる。 \begin{equation} K_{nm}^l=K^l\left(x^{(n)},x^{(m)}\right) \end{equation} 式(\ref{eq3})の第2項にある期待値の計算を進める。表記を簡単にするため、$x^{(n)}$を$x_n$などと置くことにする。 \begin{eqnarray} A&\equiv& \mathbb{E}_{l-1}\left[ \phi\left(z^{l-1}\left(x_n\right)\right)^T \phi\left(z^{l-1}\left(x_m\right)\right) \right]\\ &=& \sum_{j=1}^{N_l}\mathbb{E}_{l-1}\left[ \phi\left(z^{l-1}_j\left(x_n\right)\right) \phi\left(z^{l-1}_j\left(x_m\right)\right) \right]\\ &=& \sum_{j=1}^{N_l} \int dz_j^{l-1}(x_1)\cdots dz_j^{l-1}(x_N) \phi\left(z^{l-1}_j\left(x_n\right)\right) \phi\left(z^{l-1}_j\left(x_m\right)\right)\\ &\times& \mathcal{N}(z_j^{l-1}(x_1)\cdots z_j^{l-1}(x_N)|{\bf 0},K^{l-1}) \end{eqnarray} 上の積分は、$z_j^{l-1}(x_n)$と$z_j^{l-1}(x_m)$以外の変数については実行できる。 \begin{eqnarray} A &=& \sum_{j=1}^{N_l} \int dz_j^{l-1}(x_n) dz_j^{l-1}(x_m) \phi\left(z^{l-1}_j\left(x_n\right)\right) \phi\left(z^{l-1}_j\left(x_m\right)\right)\\ &\times& \mathcal{N}(z_j^{l-1}(x_n),z_j^{l-1}(x_m)|{\bf 0},\Psi^{l-1}) \end{eqnarray} ここで \begin{equation} \Psi^{l-1}=\left( \begin{array}{cc} K^{l-1}_{nn} & K^{l-1}_{nm} \\ K^{l-1}_{mn} & K^{l-1}_{mm} \end{array} \right) \end{equation} とした。すなわち、$A$は3つの量、$K^{l-1}_{nn},K^{l-1}_{nm}(=K^{l-1}_{mn}),K^{l-1}_{mm}$と活性化関数$\phi$に依存する量であることがわかる。以上から \begin{equation} K^l_{nm}=\sigma_b^2+\sigma_w^2F_{\phi}(K^{l-1}_{nn},K^{l-1}_{nm},K^{l-1}_{mm}) \end{equation} と書くことができる。$F$は$\phi$の形が決まれば計算できる関数である。上の式から、$K^l$は$K^{l-1}$を用いて計算できることが分かる。

参考文献


Deep Neural Networks as Gaussian Processes

2018年11月10日土曜日

1層のニューラルネットワークとガウス過程の関係

はじめに


 1層のニューラルネットワークとガウス過程の関係をまとめる。

1層のニューラルネットワーク


 下図のニューラルネットワークを考える。
各量の定義は以下の通り。 \begin{eqnarray} x&\in&\mathbb{R}^{d_{in}}\\ z^0,x^1&\in&\mathbb{R}^{N_{1}}\\ z^1&\in&\mathbb{R}^{d_{out}}\\ W^0 &\in& \mathbb{R}^{N_1\times d_{in}} \\ b^0 &\in& \mathbb{R}^{N_1} \\ W^1 &\in& \mathbb{R}^{d_{out}\times N_1} \\ b^1 &\in& \mathbb{R}^{d_{out}} \end{eqnarray} \begin{eqnarray} z^0_j(x)&=&b_j^0+\sum_{k=1}^{d_{in}}W^0_{jk}x_k\\ x^1_j(x)&=&\phi\left(z^0_j(x)\right) \\ z^1_i(x)&=&b_i^1+\sum_{j=1}^{N_{1}}W^1_{ij}x_j^1(x) \end{eqnarray} パラメータ$W^l_{ij},b^l_i\;(l=0,1)$は次のガウス分布から独立に作られると仮定する。 \begin{eqnarray} W^l_{ij}&\sim&\mathcal{N}(W^l_{ij}|0, \sigma_w^2/N_l) \label{eq1}\\ b^l_{i}&\sim&\mathcal{N}(b^l_{i}|0, \sigma_b^2)\label{eq2} \end{eqnarray}

ガウス過程との関係


 パラメータ$W^0_{jk}$と$b^0_j$はそれぞれ独立同分布から生成されると仮定したので、$x_j^1$と$x_{j^{\prime}}^1$は、$j\neq j^{\prime}$のとき、独立である。さらに、$z^1_i(x)$は独立同分布から生成される項の和であるから、$N_1\rightarrow\infty$の極限でガウス分布に従う(中心極限定理)。その分布の平均値は以下のように計算される。 \begin{eqnarray} \mathbb{E}\left[z^1_i(x)\right]&=& \mathbb{E}\left[b_i^1+\sum_{j=1}^{N_{1}}W^1_{ij}x_j^1(x)\right]\\ &=& \mathbb{E}\left[b_i^1\right]+\mathbb{E}\left[\sum_{j=1}^{N_{1}}W^1_{ij}x_j^1(x)\right]\\ &=& \mathbb{E}_1\left[b_i^1\right]+\sum_{j=1}^{N_{1}}\mathbb{E}_1\left[W^1_{ij}\right]\mathbb{E}_0\left[x_j^1(x)\right] \end{eqnarray} ここで$\mathbb{E}_l[\cdot]$は$W^l,b^l$の量で期待値を取ることを意味する。式(\ref{eq1}),(\ref{eq2})から次式が成り立つ。 \begin{eqnarray} \mathbb{E}_1\left[b_i^1\right]&=&0\\ \mathbb{E}_1\left[W^1_{ij}\right]&=&0 \end{eqnarray} したがって \begin{equation} \mathbb{E}\left[z^1_i(x)\right]=0 \end{equation} である。分散は次式のようになる。 \begin{eqnarray} \mathbb{E}\left[z^1_i(x)z^1_i(x^{\prime})\right]&=&\mathbb{E}\left[\left(b_i^1+\sum_{j=1}^{N_{1}}W^1_{ij}x_j^1(x)\right)\left(b_i^1+\sum_{k=1}^{N_{1}}W^1_{ik}x_k^1(x^{\prime})\right)\right]\\ &=&\mathbb{E}_1\left[b^1_i b^1_i\right]+\sum_j^{N_1}\sum_k^{N_1}\mathbb{E}_1\left[W^1_{ij}W^1_{ik}\right]\mathbb{E}_0\left[x^1_j(x)x^1_k(x^{\prime})\right] \end{eqnarray} 式(\ref{eq1}),(\ref{eq2})から次式が成り立つ。 \begin{eqnarray} \mathbb{E}_1\left[b^1_i b^1_i\right]&=&\sigma_b^2\\ \mathbb{E}_1\left[W^1_{ij}W^1_{ik}\right]&=&\frac{\sigma^2_w}{N_1}\delta_{jk} \end{eqnarray} 故に \begin{eqnarray} \mathbb{E}\left[z^1_i(x)z^1_i(x^{\prime})\right] &=&\sigma^2_b+\sum_j^{N_1}\sum_k^{N_1}\delta_{jk}\frac{\sigma^2_w}{N_1}\mathbb{E}_0\left[x^1_j(x)x^1_k(x^{\prime})\right]\\ &=&\sigma^2_b+\frac{\sigma^2_w}{N_1}\sum_j^{N_1}\mathbb{E}_0\left[x^1_j(x)x^1_j(x^{\prime})\right]\\ &=&\sigma^2_b+\sigma^2_w\;C(x,x^{\prime}) \\ &\equiv&K^1(x,x^{\prime}) \end{eqnarray} を得る。ただし \begin{eqnarray} C(x,x^{\prime}) &=&\frac{1}{N_1}\sum_j^{N_1}\mathbb{E}_0\left[x^1_j(x)x^1_j(x^{\prime})\right]\\ &=&\frac{1}{N_1}\mathbb{E}_0\left[\sum_j^{N_1}x^1_j(x)x^1_j(x^{\prime})\right]\\ &=&\frac{1}{N_1}\mathbb{E}_0\left[x^1(x)^T\cdot x^1(x^{\prime})\right] \end{eqnarray} と置いた。いま入力ベクトル$x$が$N$個あるとする。これを${\bf x}=(x^{(1)},\cdots,x^{(N)})$と書くことにする。このとき$z_i^1(x)$も$N$個並べることができる。 \begin{equation} {\bf z}_i^1=(z_i^1(x^{(1)}),\cdots,z_i^1(x^{(N)})) \end{equation} この${\bf z}_i^1$が、$N_1\rightarrow\infty$のとき、ガウス分布$\mathcal{N}({\bf z}_i^1|{\bf 0},K^1)$に従う量と等価となる。そして、共分散行列$K^1$の成分$K^1_{nm}$が次式で与えられる。 \begin{equation} K^1_{nm}=K^1(x^{(n)},x^{(m)}) \end{equation} ${\bf 0}\in \mathbb{R}^N, K^1\in \mathbb{R}^{N\times N}$である。

 同じようにして$z_i^0(x)$の平均値を求めると \begin{eqnarray} \mathbb{E}\left[z^0_i(x)\right]&=& \mathbb{E}\left[b_i^0+\sum_{k=1}^{d_{in}}W^0_{ik}x_k\right]\\ &=& \mathbb{E}\left[b_i^0\right]+\mathbb{E}\left[\sum_{k=1}^{d_{in}}W^0_{ik}x_k\right]\\ &=& \mathbb{E}_0\left[b_i^0\right]+\sum_{k=1}^{d_{in}}\mathbb{E}_0\left[W^0_{ij}\right]x_k\\ &=&0 \end{eqnarray} となり、分散は \begin{eqnarray} \mathbb{E}\left[z^0_i(x)z^0_i(x^{\prime})\right] &=&\mathbb{E}\left[\left(b_i^0+\sum_{k=1}^{d_{in}}W^0_{ik}x_k\right)\left(b_i^0+\sum_{j=1}^{d_{in}}W^0_{ij}x_j^{\prime}\right)\right]\\ &=&\mathbb{E}\left[b_i^0 b_i^0\right]+ \sum_{k=1}^{d_{in}}\sum_{j=1}^{d_{in}}\mathbb{E}_0\left[ W^0_{ik}W^0_{ij} \right]x_k x_j^{\prime}\\ &=& \sigma_b^2+ \sum_{k=1}^{d_{in}}\sum_{j=1}^{d_{in}} \frac{\sigma_w^2}{d_{in}}\delta_{kj} x_k x_j^{\prime}\\ &=& \sigma_b^2+ \frac{\sigma_w^2}{d_{in}}\sum_{k=1}^{d_{in}} x_k x_k^{\prime}\\ &=& \sigma_b^2+ \frac{\sigma_w^2}{d_{in}} x^T\cdot x^{\prime}\\ &\equiv&K^0(x,x^{\prime}) \end{eqnarray} となる。 すなわち、${\bf z}_i^0$は、ガウス分布$\mathcal{N}({\bf z}_i^0|{\bf 0},K^0)$に従う量となる

参考文献


Deep Neural Networks as Gaussian Processes

2018年11月6日火曜日

ベイズ推論とガウス過程

はじめに


 ベイズ推論とガウス過程の関係を、回帰を例に取り確認する。

ベイズ線形回帰


 まず最初にベイズ推論を用いて線形回帰を考察する。いま観測値$X$と$Y$が与えられているとする。 \begin{eqnarray} X&=&\{x_1,\cdots,x_N\},x_n\in \mathbb{R}^M \\ Y&=&\{y_1,\cdots,y_N\},\;y_n\in\mathbb{R} \end{eqnarray} このとき次の尤度を考える。 \begin{equation} p(y_n|x_n,w)=\mathcal{N}(y_n|w^T\phi(x_n),\beta^{-1}) \end{equation} ここで、$\mathcal{N}(y_n|\mu,\sigma^2)$は平均$\mu$、分散$\sigma^2$の正規分布を表す。$w$はパラメータ($w\in\mathbb{R}^D$)、$\phi(\cdot)$は$M$次元空間内の点を$D$次元空間内の点に射影する関数である。次に、$w$の事前分布を次のように導入する。 \begin{equation} p(w)=\mathcal{N}(w|0,\alpha^{-1}I_D) \end{equation} $I_D$は$D\times D$の単位行列である。このとき事後分布$p(w|X,Y)$は厳密に計算でき \begin{eqnarray} p(w|X,Y)&=&\mathcal{N}(w|m,S) \\ m&=&S\beta\Phi^Ty \in \mathbb{R}^D\\ S&=&(\alpha I_D+\beta\Phi^T\Phi)^{-1} \in \mathbb{R}^{D\times D} \end{eqnarray} を得る。$y=(y_1,\dots,y_N)^T$、$\Phi$は$\Phi_{nd}=\phi_d(x_n)$を成分に持つ$N\times D$行列である。さらに、この事後分布から予測分布を求めることができる。 \begin{eqnarray} p(y_{*}|x_{*},X,Y)&=&\mathcal{N}(y_{*}|m^T\phi(x_*),\sigma^2(x_*)) \label{eq1}\\ \sigma^2(x_*)&=&\beta^{-1}+\phi^T(x_*)S\phi(x_*) \end{eqnarray} 以上が、ベイズ推論による線形回帰の結果である。

ガウス過程による回帰


 いま観測値$X$と$Y$が与えられているとする。 \begin{eqnarray} X&=&\{x_1,\cdots,x_N\},x_n\in \mathbb{R}^M \\ Y&=&\{y_1,\cdots,y_N\},\;y_n\in\mathbb{R} \end{eqnarray} 観測値$X$と$Y$の間には次の関係があるとする。 \begin{equation} y_n=z(x_n)+\epsilon_n \end{equation} $\epsilon_n$が正規分布に従うノイズであれば次の尤度を考えることができる。 \begin{equation} p(y_n|z_n)=\mathcal{N}(y_n|z_n,\beta^{-1}) \end{equation} ただし、$z_n=z(x_n)$とした。観測値は独立同分布から生成されると仮定すれば次式を得る。 \begin{eqnarray} p(y|z) &=&p(y_1,\dots,y_N|z_1,\cdots,z_N) \\ &=&\prod_{n=1}^N p(y_n|z_n) \\ &=&\prod_{n=1}^N \mathcal{N}(y_n|z_n,\beta^{-1}) \\ &=&\mathcal{N}(y|z,\beta^{-1}I_N) \end{eqnarray} ただし、$y=(y_1,\dots,y_N)^T$、$z=(z_1,\dots,z_N)^T$と置いた。$I_N$は$N\times N$の単位行列である。ここで、$z$がガウス過程から生成される量であると仮定すると、次の事前分布を使うことができる。 \begin{equation} p(z)=\mathcal{N}(z|0,K) \end{equation} 共分散行列$K$のサイズは$N\times N$であり、その成分は次式で定義される。 \begin{equation} K_{nm}=k(x_n,x_m) \end{equation} $k(x_n,x_m)$は2点から1つのスカラー量が決まる関数であり、カーネル関数と呼ばれる。これらを用いて、$p(y,z)$を$z$について周辺化する。 \begin{eqnarray} p(y)&=&\int dz\;p(y|z)\;p(z)\\ &=&\mathcal{N}(y|0,C) \end{eqnarray} $C$は$N\times N$の行列であり、その成分は次式で与えられる。 \begin{equation} C_{nm}=k(x_n,x_m)+\beta^{-1}\delta_{nm} \end{equation} 予測分布は次のようなる。 \begin{eqnarray} p(y_*|x_*,X,Y)&=&\mathcal{N}(y_*|m^{'}(x_*),\sigma^{'\;2}(x_*)) \label{eq2}\\ m^{'}(x_*)&=&k^TC^{-1}y \\ \sigma^{'\;2}(x_*)&=&c-k^TC^{-1}k \end{eqnarray} ただし \begin{eqnarray} k&=&(k(x_n,x_*),\cdots,k(x_N,x_*))^T \\ c&=&k(x_*,x_*)+\beta^{-1} \end{eqnarray} である。以上がガウス過程による回帰の結果である。

ベイズ推論とガウス過程の関係


 いま、カーネル関数が次式で定義されるとする。 \begin{equation} k(x_n,x_m)=\alpha^{-1}\phi^T(x_n)\phi(x_m) \end{equation} このとき、式(\ref{eq1})と(\ref{eq2})は等しくなる。

 ベイズ推論による解法では$D$次元空間を考えるが、ガウス過程による解法では$N$次元空間を考える。前者では$D\times D$行列の逆行列を、後者では$N\times N$行列の逆行列を計算する必要がある。それぞれの計算量は$\mathcal{O}(D^3)$と$\mathcal{O}(N^3)$である。$D\lt N$である場合、パラメータ空間で計算した方が計算量は少なくなる。ガウス過程を使う利点は、$K$として様々なカーネルを使うことができる点である。

ガウス過程のハイパーパラメータの決定


 いま共分散行列$K$がハイパーパラメータ$\theta$でモデル化されているとする。 \begin{equation} p(z|\theta)=\mathcal{N}(z|0,K_{\theta}) \end{equation} このとき \begin{eqnarray} p(y|\theta)&=&\int dz\;p(y|z)\;p(z|\theta)\\ &=&\mathcal{N}(y|0,C_{\theta}) \\ C_{\theta,nm}&=&K_{\theta,nm}+\beta^{-1}\delta_{nm} \end{eqnarray} が成り立つ。事前分布$p(\theta)$を考えると次式が成り立つ。 \begin{equation} p(\theta|y)=\frac{p(y|\theta)p(\theta)}{p(y)} \end{equation} 一般的に右辺分母は解析的に計算できないので、何らかの近似が必要になる。書籍「Pythonによるベイズ統計モデリング PyMCでのデータ分析実践ガイド」の第8章にPyMC3を用いた解法が掲載されている。

2018年10月8日月曜日

ベルマン最適方程式とQ-LearningとDQN

はじめに


  先の記事でベルマン最適方程式について書いた。今回は、この方程式とQ-Learningとを関連づけてみたい。最後にDQNについて少し言及する。

ベルマン最適方程式


 行動価値関数に対するベルマン最適方程式は次式で与えられた。 \begin{equation} Q^{*}(s,a)=\sum_{s^{\prime}} P(s^{\prime}|s,a)\left[r(s,a,s^{\prime})+\gamma \max_{a^{\prime}}Q^{*}(s^{\prime},a^{\prime})\right] \label{eq1} \end{equation} ここで、$Q^{*}(s,a)$は最適行動価値関数、$s$と$a$はそれぞれ状態と行動を表す観測値を表す。$P(s^{\prime}|s,a)$は$(s,a)$が実現しているとき$s^{\prime}$が起こる確率、$r(s,a,s^{\prime})$は$(s,a,s^{\prime})$が実現するときの報酬を表す。$\gamma$は収益を定義するときに導入した割引率である。

Q-Learningとの関係


 Q-Learningでは行動価値関数$Q(s,a)$を次式で更新していく。 \begin{equation} Q(s,a) \leftarrow (1-\alpha)Q(s,a)+\alpha\left[r(s,a,s^{\prime})+\gamma\max_{a^{\prime}}Q(s^{\prime},a^{\prime})\right] \label{eq3} \end{equation} ここで、$\alpha$は学習率である。現在の状態$s$とこの状態の下で選択した行動$a$から次の状態$s^{\prime}$ を決定論的に決めたあと、上式を用いて$Q(s,a)$を更新することに注意する。上式右辺を変形すると次式を得る。 \begin{equation} Q(s,a) \leftarrow Q(s,a)+\alpha\left[r(s,a,s^{\prime})+\gamma\max_{a^{\prime}}Q(s^{\prime},a^{\prime})-Q(s,a)\right] \end{equation} 更新を繰り返し$Q(s,a)$が収束したとき、上式右辺の第2項は0になるはずである。 \begin{equation} r(s,a,s^{\prime})+\gamma\max_{a^{\prime}}Q(s^{\prime},a^{\prime})-Q(s,a)=0 \end{equation} これより \begin{equation} Q(s,a)=r(s,a,s^{\prime})+\gamma\max_{a^{\prime}}Q(s^{\prime},a^{\prime}) \label{eq2} \end{equation} を得る。ここで、式(\ref{eq1})と式(\ref{eq2})を見比べて見ると、とても良く似ていることに気づく。違いは、式(\ref{eq1})では遷移確率$P(s^{\prime}|s,a)$が考慮されており、$s$の次に取り得る状態$s^{\prime}$について期待値が取られていることである。一方、式(\ref{eq2})では状態$s^{\prime}$は決定論的に決められる。

Deep Q Network(DQN)


 Q-Learningでは式(\ref{eq3})を用いて$Q$を更新する。アルゴリズムの手順は以下のようなる。
  1. 現在の状態$s$が与えられる。
  2. 現在までに確定している$Q$値の表を見て、$s$のとき$Q$値を最大にする行動$a$を求める。ただし、ランダムに行動を選択する仕組みも取り入れ、他にあるかしれないより良い行動が見つかるようにする。
  3. $(s,a)$の下で次の状態$s^{\prime}$を決める。
  4. 現在までに確定している$Q$値の表を見て、$s^{\prime}$のとき最大となる$Q$値を求め、式(\ref{eq3})を用いて$Q$値を更新する。
Q-Learningでは更新式(\ref{eq3})を用いて$Q$値の表を作成し、入力$s$に対する出力$a$を決定する。一方、DQNでは状態に対する行動をニューラルネットワークで決めることになる。

参考文献


2018年8月5日日曜日

Adversarial Variational Bayes 〜理論(その2)〜

はじめに


 前回の投稿では、論文Adversarial Variational Bayesの前半部分をまとめた。今回は後半部分をまとめる。

前半部分の結論


 従来のVAEの目的関数は \begin{equation} (\theta^{*},\phi^{*})=\arg\max_{\theta,\phi} {\rm E}_{p_{d}(x)}\left[ -{\rm KL} \left( q_{\phi}(z|x), p(z) \right) + {\rm E}_{q_{\phi}(z|x)} \left[\ln{p_{\theta}(x|z)}\right] \right] \label{eq1} \end{equation} である。本論文の前半部分では式(\ref{eq1})から次の2つの目的関数を導出した。 \begin{eqnarray} T^{*}(x,z)&=&\arg{\max_{T}}\; {\rm E}_{p_{d}(x)} \left[ {\rm E}_{q_{\phi}(z|x)} \left[ \ln{\sigma\left(T(x, z)\right)} \right] + {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right] \label{eq2}\\ (\theta^{*},\phi^{*})&=&\arg\max_{\theta,\phi} {\rm E}_{p_{d}(x)}\;{\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x|z)}-T^{*}(x,z) \right] \label{eq3} \end{eqnarray} ここで、$\sigma(\cdot)$はシグモイド関数である。実際に計算を行うと、式(\ref{eq2})から最適な値$T^*(x,z)$を求めることは難しいことが分かる。その原因は、式(\ref{eq1})の右辺にあるKullback Leibler divergence にある。この量は、事後分布$q_{\phi}(z|x)$の形状を事前分布$p(z)$の形状に近づけようとする。一般にこれら2つの分布は大きく異なる形状を持つため、$q_{\phi}(z|x)\approx p(z)$の実現は難しい。また、事後分布を事前分布に近づける操作は、事後分布から観測値$x$の依存性を失くすことに相当するため、Bayes推論の立場から考えても妥当ではない。本論文の後半部分ではこの操作を改善する。

後半部分の内容


 事後分布$q_{\phi}(z|x)$に良く似た形状を持ち、かつ、解析計算が可能な補助分布(auxiliary distribution)$r_{\alpha}(z|x)$を導入する。例えば、$r_{\alpha}(z|x)$として正規分布を考えることができる。この補助分布を用いて、式(\ref{eq1})の右辺を次のように書き換える。 \begin{equation} {\rm E}_{p_{d}(x)}\left[ -{\rm KL} \left( q_{\phi}(z|x), r_{\alpha}(z|x) \right) + {\rm E}_{q_{\phi}(z|x)} \left[ -\ln{r_{\alpha}(z|x)} +\ln{p_{\theta}(x,z)}\right] \right] \label{eq4} \end{equation} $r_{\alpha}(z|x)$を$q_{\phi}(z|x)$を近似する関数と置いたので、式(\ref{eq4})に現れるKullback Leibler divergenceの値は、式(\ref{eq1})のそれよりもずっと容易に小さな値とすることができる。

 いま、補助分布として次の正規分布を仮定する。 \begin{equation} r_{\alpha}(z|x)=\mathcal{N}(z|\mu(x),\sigma_s^2(x)) \end{equation} これは、変数変換 \begin{equation} \tilde{z}=\frac{z-\mu(x)}{\sigma_s(x)} \end{equation} の下で \begin{equation} r_{\alpha}(z|x)=\frac{1}{\sigma_s}\mathcal{N}(\tilde{z}|0,1)\equiv \frac{1}{\sigma_s}r_0(\tilde{z}) \end{equation} となるから \begin{eqnarray} {\rm KL} \left( q_{\phi}(z|x), r_{\alpha}(z|x) \right) &=& \int dz\;q_{\phi}(z|x)\ln{\frac{q_{\phi}(z|x)}{r_{\alpha}(z|x)}} \label{eq5} \\ &=& \int d\tilde{z}\;\sigma_s q_{\phi}(\sigma_s \tilde{z}+\mu|x) \ln{ \frac {\sigma_s q_{\phi}(\sigma_s\tilde{z}+\mu|x)} {r_0(\tilde{z})} }\\ &=& {\rm KL} \left( \tilde{q}_{\phi}(\tilde{z}|x), r_0(\tilde{z}) \right) \label{eq6} \end{eqnarray} となる。ただし、$\tilde{q}_{\phi}(\tilde{z}|x)\equiv \sigma_s q_{\phi}(\sigma_s\tilde{z}+\mu|x)$とした。$\tilde{q}_{\phi}(\tilde{z}|x)$は \begin{equation} \int d\tilde{z}\tilde{q}_{\phi}(\tilde{z}|x)=1 \end{equation} を満たす新たな分布である。さらに、次式が成り立つ。 \begin{eqnarray} {\rm E}_{q_{\phi}(z|x)} \left[ \ln{r_{\alpha}(z|x)} \right] &=& \int dz\;q_{\phi}(z|x)\ln{r_{\alpha}(z|x)} \\ &=& \int d\tilde{z}\;\tilde{q}_{\phi}(\tilde{z}|x)\ln{r_{0}(\tilde{z})}-\ln{\sigma_s} \\ &=& {\rm E}_{\tilde{q}_{\phi}(\tilde{z}|x)} \left[ \ln{r_{0}(\tilde{z})} \right]+{\rm const.} \end{eqnarray} 以上から、式(\ref{eq4})は、定数項を無視して、次式に変換される。 \begin{equation} {\rm E}_{p_{d}(x)}\left[ -{\rm KL} \left( \tilde{q}_{\phi}(\tilde{z}|x), r_{0}(\tilde{z}) \right) + {\rm E}_{\tilde{q}_{\phi}(\tilde{z}|x)} \left[ -\ln{r_{0}(\tilde{z})} \right]+ {\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x,z)}\right] \right] \label{eq7} \end{equation} ここで、$-{\rm KL} \left( \tilde{q}_{\phi}(\tilde{z}|x), r_{0}(\tilde{z}) \right)$に注目して、論文の前半部分で用いた議論を繰り返し、$T(x,\tilde{z})$を導入すると、式(\ref{eq7})は次の2式と等価となる。 \begin{eqnarray} T^{*}(x,\tilde{z})&=&\arg{\max_{T}} \;{\rm E}_{p_{d}(x)} \left[ {\rm E}_{\tilde{q}_{\phi}(\tilde{z}|x)} \left[ \ln{\sigma\left(T(x, \tilde{z})\right)} \right] + {\rm E}_{r_0(\tilde{z})} \left[ \ln{\left(1-\sigma\left(T(x, \tilde{z})\right)\right)} \right] \right] \label{eq8}\\ (\theta^{*},\phi^{*}) &=& \arg\max_{\theta,\phi} \;{\rm E}_{p_{d}(x)} \left[ {\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x,z)} \right] + {\rm E}_{\tilde{q}_{\phi}(\tilde{z}|x)} \left[ -T^*(x,\tilde{z})-\ln{r_0(\tilde{z})} \right] \right] \label{eq9} \end{eqnarray} これらに、再パラメータ化トリックを適用する。すなわち、 \begin{equation} z\sim q_{\phi}(z|x) \end{equation} を \begin{eqnarray} \epsilon&\sim&p(\epsilon)\\ z&=&z_{\phi}(x,\epsilon) \end{eqnarray} に置き換える。このとき \begin{equation} \tilde{z}\sim \tilde{q}_{\phi}(\tilde{z}|x) \end{equation} は \begin{eqnarray} \tilde{z}&=&\frac{z_{\phi}(x,\epsilon)-\mu}{\sigma_s}\equiv\tilde{z}_{\phi}(x,\epsilon) \end{eqnarray} に置き換わる。また、分布$r_0(\tilde{z})$は標準正規分布であるから \begin{equation} \tilde{z} \sim r_0(\tilde{z}) \end{equation} は \begin{eqnarray} \eta&\sim&\mathcal{N}(\eta|0,1)\equiv p(\eta) \\ \tilde{z}&=&\eta \end{eqnarray} と書くことができる。以上から式(\ref{eq8}),(\ref{eq9})は次式に変形される。 \begin{eqnarray} T^{*}(x,\tilde{z}) &=& \arg{\max_{T}} \;{ \rm E}_{p_{d}(x)} \Bigl[ {\rm E}_{p(\epsilon)} \left[ \ln{\sigma\left(T(x, \tilde{z}_{\phi}(x,\epsilon))\right)} \right] + {\rm E}_{p(\eta)} \left[ \ln{\left(1-\sigma\left(T(x, \eta)\right)\right)} \right] \Bigr] \label{eq10}\\ (\theta^{*},\phi^{*}) &=& \arg\max_{\theta,\phi} \;{\rm E}_{p_{d}(x)} \Bigl[ {\rm E}_{p(\epsilon)} \left[ \ln{p_{\theta}(x,z_{\phi}(x,\epsilon))} \right] + {\rm E}_{p(\epsilon)} \left[ -T^*(x,\tilde{z}_{\phi}(x,\epsilon))-\ln{r_0(\tilde{z}_{\phi}(x,\epsilon))} \right] \Bigr] \label{eq11} \end{eqnarray} $r_0$は標準正規分布であるから、式(\ref{eq11})の最後の項は \begin{equation} -\ln{r_0(\tilde{z}_{\phi}(x,\epsilon))}=\frac{1}{2}\|\tilde{z}_{\phi}(x,\epsilon)\|^2+{\rm const.} \end{equation} となる。最適化を行う際は定数項は無視すれば良い。また、$p_{\theta}(x,z)$の計算は \begin{equation} p_{\theta}(x,z)=p_{\theta}(x|z)p(z) \end{equation} を利用すれば良い。ここまでの議論を反映したアルゴリズムが論文に掲載されている。

2018年7月1日日曜日

Adversarial Variational Bayes 〜理論と実装〜

はじめに


 論文Adversarial Variational Bayesの前半部分の解説とChainerによる実装を行う。

論文概要


 深層学習において生成モデルを求める代表的な手法は、Generative Adversarial Network(GAN)とVariational AutoEncoder(VAE)である。最初にGANが最適化する目的関数を示す。 \begin{equation} \min_{\theta}\max_{\phi} \left\{ {\rm E}_{p_{d}(x)}\left[\log{D_{\phi}(x)}\right]+{\rm E}_{p(z)}\left[1-\log{D_{\phi}\left(G_{\theta}(z)\right)}\right] \right\} \end{equation} ここで、$p_d(x)$と$p(z)$はそれぞれ観測値$x$と潜在変数$z$の分布を表す。 ${\rm E}_{p}\left[\cdot\right]$は分布$p$による期待値である。$D_{\phi}$は識別器、$G_{\theta}$は生成器であり、それぞれパラメータ$\phi$と$\theta$で象徴される重みを持つネットワークである。上の式は次の2つの最適化を行うことを表す。
  1. $\theta$を固定して考える。観測値$x$と生成器$G_{\theta}(z)$の出力値を識別器$D_{\phi}$に与えたとき、$D_{\phi}(x)$の値を最大にし、$D_{\phi}(G_{\theta}(z))$の値を最小にするような$\phi$を見つける。
  2. $\phi$を固定して考える。観測値$x$と生成器$G_{\theta}(z)$の出力値を識別器$D_{\phi}$に与えたとき、$D_{\phi}(x)$の値を最小にし、$D_{\phi}(G_{\theta}(z))$の値を最大にするような$\theta$を見つける。
訓練の後、生成モデル$G_{\theta}(z)$を得ることになる。一方、VAEが最適化すべき式は次式である。 \begin{equation} \max_{\theta} \max_{\phi} \left\{ -{\rm KL} \left( q_{\phi}(z|x), p(z) \right) + {\rm E}_{q_{\phi}(z|x)} \left[\ln{p_{\theta}(x|z)}\right] \right\} \end{equation} ${\rm KL}$はKullback Leibler divergence、$q_{\phi}(z|x)$は観測値$x$が与えられた下での潜在変数$z$の事後確率を近似する分布である。上式第1項は正則化項、第2項は尤度を表す。すなわち、正則化項を考慮して尤度を最大化することになる。訓練の後、$q_{\phi}(z|x)$が推論モデル(inference model)、$p_{\theta}(x|z)$が生成モデル(generative model)になる。VAEの処理の流れを以下に示す。
GANとVAEを比べると、後者からは生成モデル以外に推論モデルも求めることができる。また、GANで決まる生成器は決定論的な関数であるが、VAEから求まる生成器は確率モデルである。さらに、GANには尤度に相当する量が存在せず、VAEには識別器に相当する仕組みが存在しない。

 GANを自然画像に適用すると、シャープな画像を生み出すが、VAEではボケた画像になることが知られている。これは、VAEの推論モデル$q_{\phi}(z|x)$が真の事後確率$p(z|x)$を再現できていないためである。今回紹介する文献は、任意の形状を持つ推論モデルを扱うことのできる手法(Adversarial Variational Bayes:AVB)を提案している。本文献の主張をまとめると次のようになる。
  1. 任意の複雑な推論モデルを扱うことができる。
  2. VAEの目的関数に数学的な変換を施すことによりGAN的な目的関数を導出する。すなわち、VAEに識別器に相当する量が導入される。
  3. 上の目的関数は、ある極限の下で、VAEの元の目的関数と厳密に一致する。

論文詳細


 観測値$x$の分布$p(x)$を考え、潜在変数$z$を導入する。 \begin{eqnarray} \ln{p(x)} &=&\ln\int dz\;p(x|z)p(z) \\ &=&\ln\int dz\;q(z|x) \frac{p(x|z)p(z)}{q(z|x)} \end{eqnarray} Jensenの不等式を用いて \begin{eqnarray} \ln{p(x)} &\geqq&\int dz\;q(z|x) \ln\frac{p(x|z)p(z)}{q(z|x)} \\ &=&\int dz\;q(z|x)\ln{p(x|z)}-\int dz\;q(z|x)\ln{\frac{q(z|x)}{p(z)}} \\ &=&{\rm E}_{q(z|x)}\left[\ln{p(x|z)}\right]-{\rm KL}\left[q(z|x),p(z)\right] \\ &\equiv&\mathcal{L} \end{eqnarray} を得る。上式右辺はEvidence Lower Bound(ELOB)と呼ばれる量である。ここで、次の量を考える。 \begin{eqnarray} {\rm KL}\left(q(z|x),p(z|x)\right) &=&\int dz\;q(z|x)\ln{\frac{q(z|x)}{p(z|x)}} \\ &=&\int dz\;q(z|x)\ln{\frac{p(x)q(z|x)}{p(x,z)}} \\ &=&\int dz\;q(z|x)\ln{\frac{p(x)q(z|x)}{p(x|z)p(z)}} \\ &=&\int dz\;q(z|x)\ln{p(x)}-\int dz\;q(z|x)\ln{p(x|z)}+{\rm KL}(q(z|x),p(z)) \\ &=&\ln{p(x)}-\int dz\;q(z|x)\ln{p(x|z)}+{\rm KL}(q(z|x),p(z)) \\ \end{eqnarray} 上式は先に定義したELBO($\mathcal{L}$)を用いると \begin{equation} {\rm KL}\left(q(z|x),p(z|x)\right)=\ln p(x)-\mathcal{L} \end{equation} と書くことができる。すなわち次式が成り立つ。 \begin{equation} \ln p(x)=\mathcal{L}+{\rm KL}\left(q(z|x),p(z|x)\right) \end{equation} $q(z|x)$は分布$p(z|x)$を近似するために導入された分布であり、パラメータ$\phi$を持つとする。一方、分布$p(x|z)$を表現するモデルはパラメータ$\theta$を持つとする。これらパラメータを顕に書くと \begin{eqnarray} \ln p(x)&=&\mathcal{L}+{\rm KL}\left(q_{\phi}(z|x),p(z|x)\right) \label{eq2}\\ \ln p(x)&\geqq&\mathcal{L} \label{eq1}\\ \mathcal{L}&=&{\rm E}_{q_{\phi}(z|x)}\left[\ln{p_{\theta}(x|z)}\right]-{\rm KL}\left[q_{\phi}(z|x),p(z)\right] \end{eqnarray} となる。VAEでは式(\ref{eq1})の右辺を$\phi$と$\theta$について最大化する。式(\ref{eq2})から、$\ln p(x)=\mathcal{L}$が成り立つのは$q_{\phi}(z|x)=p(z|x)$となる時であることが分かるが、一般にこの等式が成り立つことはない。通常のVAEでは実際に計算を進める際、$q_{\phi}(z|x)$を正規分布で近似する。本文献では$q_{\phi}(z|x)$に対してそのような近似を行わない。

 $\mathcal{L}$は以下のように書き換えることができる。 \begin{eqnarray} \mathcal{L} &=&{\rm E}_{q_{\phi}(z|x)}\left[\ln{p_{\theta}(x|z)}\right]-{\rm KL}\left[q_{\phi}(z|x),p(z)\right] \\ &=&{\rm E}_{q_{\phi}(z|x)}\left[\ln{p_{\theta}(x|z)}\right]-\int dz\;q_{\phi}(z|x)\ln{\frac{q_{\phi}(z|x)}{p(z)}} \\ &=&{\rm E}_{q_{\phi}(z|x)}\left[\ln{p_{\theta}(x|z)}\right]-\int dz\;q_{\phi}(z|x) \left[ \ln{q_{\phi}(z|x)}-\ln{p(z)} \right] \\ &=&{\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x|z)}-\ln{q_{\phi}(z|x)}+\ln{p(z)} \right]\\ \end{eqnarray} 従って、VAEの目的関数は \begin{equation} \max_{\theta}\max_{\phi} {\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x|z)}-\ln{q_{\phi}(z|x)}+\ln{p(z)} \right] \label{eq4} \end{equation} となる。ここまでは通常のVAEである。本文献では、ここで、関数$T^*(x,z)$を次式で導入する。 \begin{equation} T^*(x,z)=\arg\max_{T} \left\{ {\rm E}_{q_{\phi}(z|x)} \left[ \ln{\sigma\left(T(x, z)\right)} \right] + {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right\} \label{eq3} \end{equation} この式は、$q_{\phi}(z|x)$からサンプリングされた$z$のとき$T(x,z)$を大きくし、$p(z)$からサンプリングされた$z$のとき$T(x,z)$を小さくすることを意味する。すなわち、$T(x,z)$は$q_{\phi}(x|z)$と$p(z)$を識別する識別器である。実際に$T^*(x,z)$を求めるため、式(\ref{eq3})の期待値を積分に書き換える。 \begin{equation} \int dz \left[ q_{\phi}(z|x) \ln{\sigma\left(T(x, z)\right)} + p(z) \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \end{equation} 上式の被積分関数の最大値を求めるため、$a=q_{\phi}(z|x)$、$b=p(z)$、$t=\sigma\left(T(x,z)\right)$として次式を考える。 \begin{equation} y=a\ln{t}+b\ln{(1-t)} \end{equation} 両辺$t$で微分して \begin{equation} \frac{dy}{dt}=\frac{a}{t}-\frac{b}{1-t} \end{equation} 右辺を0とおいて \begin{equation} t=\frac{a}{a+b} \end{equation} を得る。各変数を元に戻して計算すると \begin{equation} T(x,z)=\ln{q_{\phi}(z|x)}-\ln{p(z)}\equiv T^{*}(x,z) \end{equation} を得る。$T^{*}$を用いると、式(\ref{eq4})は以下のように書くことができる。 \begin{equation} \max_{\theta}\max_{\phi} {\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x|z)}-T^{*}(x,z) \right] \label{eq5} \end{equation} ただし \begin{equation} T^{*}(x,z)=\arg{\max_{T}} \left\{ {\rm E}_{q_{\phi}(z|x)} \left[ \ln{\sigma\left(T(x, z)\right)} \right] + {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right\} \label{eq6} \end{equation} である。式(\ref{eq5})と(\ref{eq6})に再パレメータ化トリックを用いると \begin{eqnarray} T^{*}(x,z)&=&\arg{\max_{T}} \left\{ {\rm E}_{p(\epsilon)} \left[ \ln{\sigma\left(T(x, z_{\phi}(x,\epsilon)\right)} \right] + {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right\} \\ (\theta^*,\phi^*)&=&\arg\max_{\theta,\phi} {\rm E}_{p(\epsilon)} \left[ \ln{p_{\theta}\left(x|z_{\phi}\left(x,\epsilon\right)\right)}-T^{*}\left(x,z_{\phi}\left(x,\epsilon\right)\right) \right] \label{eq7} \end{eqnarray} を得る。$p(\epsilon)$は標準正規分布とすれば良い。観測値の分布$p_d(x)$による期待値も考慮して最終的に次の2つの目的関数を考えることになる。 \begin{eqnarray} T^{*}(x,z)&=&\arg{\max_{T}} \left\{ {\rm E}_{p_d(x)} {\rm E}_{p(\epsilon)} \left[ \ln{\sigma\left(T(x, z_{\phi}(x,\epsilon)\right)} \right] + {\rm E}_{p_d(x)} {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right\} \label{eq9}\\ (\theta^*,\phi^*)&=&\arg\max_{\theta,\phi} {\rm E}_{p_d(x)}{\rm E}_{p(\epsilon)} \left[ \ln{p_{\theta}\left(x|z_{\phi}\left(x,\epsilon\right)\right)}-T^{*}\left(x,z_{\phi}\left(x,\epsilon\right)\right) \right] \label{eq8} \end{eqnarray} 上の2式が本文献の最初の手法である。$T$を最適化したあと$(\theta,\phi)$を最適化する。アルゴリズムが文献に掲載されている。
ニューラルネットワークで表現される量は次の3つである。
  1. $z_{\phi}(x,\epsilon)$
  2. $p_{\theta}(x|z)$のパラメータ
  3. $T_{\psi}(x,z)$
文献の図5は、本アルゴリズムを評価するために用いられた実験データである。4次元データであり1つの成分だけが1、残りが0のデータである。
図6は$z_{\phi}(z|x)$から求まる潜在変数$z$の分布である。実験データは4種類のベクトルからなる。それぞれのデータごとに対応する$z$が色分けされている。同じデータに対する従来のVAEの結果も示されている。AVBの方が密に分布していることが分かる。 VAEの方はそれぞれの分布がガウス分布しているように見える。
下の表1は、各種評価値である。
表中のlog-likelihoodは$\ln{p_{\theta}(x|z)}$を、reconstruction errorは観測値$x$と再構築した$x$の間のクロスエントロピーを表す。いずれの値も従来法より精度が上がっていることが分かる。ここまでが本文献の前半の内容である。

Chainerによる実装


 ここからは、Chainerを用いて実際に実装を行い、上の図5に示された実験データから図6(b)の結果を再現するまでの様子を記載する。ソースコードはここにある。既存の実装として参考にしたのは以下の3つである。
  1. https://gist.github.com/poolio/b71eb943d6537d01f46e7b20e9225149
  2. https://github.com/gdikov/adversarial-variational-bayes
  3. https://github.com/LMescheder/AdversarialVariationalBayes
1はtensorflow、2はkeras、3は本文献の著者によるtensorflowによる実装である。

Encoder


 まず最初にEncoder($z_{\phi}(x,\epsilon)$)の実装を示す(encoder.py内のコードである)。
class Encoder_2(chainer.Chain):

    def __init__(self, x_dim, eps_dim, h_dim=512):
        super(Encoder_2, self).__init__()
        with self.init_scope():
            self.l1 = L.Linear(x_dim + eps_dim, h_dim, initialW=xavier.Xavier(x_dim + eps_dim, h_dim))
            self.l2 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.l3 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.l4 = L.Linear(h_dim, eps_dim, initialW=xavier.Xavier(h_dim, eps_dim))

    def update(self, updates):
        update_links(self, updates)

    def __call__(self, xs, es, activation=F.relu):
        xs = 2 * xs - 1 
        h = F.concat((xs, es), axis=1)

        h = self.l1(h)
        h = activation(h)

        h = self.l2(h)
        h = activation(h)

        h = self.l3(h)
        h = activation(h)

        h = self.l4(h)
        return h
  1. 入力値は観測値xsと標準正規分布からのサンプル値esである。
  2. 15行目:xsは0と1である。これを-1と1に置き換える。
  3. 16行目:xsとesを連結する。
  4. あとは、全結合層と活性化関数の繰り返しである。
  5. 出力層には活性化関数を適用しない。
  6. 活性化関数としてreluを採用した。

Decoder


 次はDecoderである(decoder.py内のコードである)。実験データは0と1から構成されるので$p_{\theta}(x|z)$としてBernoulli分布を用いる。 \begin{eqnarray} p_{\theta}(x|z) &=&{\rm Bern}(x|\mu(z))\\ &=&\mu(z)^x(1-\mu(z))^{1-x} \end{eqnarray} Decoderはパラメータ$\mu(z)$を計算する。
class Decoder_1(chainer.Chain):

    def __init__(self, z_dim, x_dim=1, h_dim=512):
        super(Decoder_1, self).__init__()
        with self.init_scope():
            self.l1 = L.Linear(z_dim, h_dim, initialW=xavier.Xavier(z_dim, h_dim))
            self.l2 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.l3 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.l4 = L.Linear(h_dim, x_dim, initialW=xavier.Xavier(h_dim, x_dim))

    def update(self, updates):
        update_links(self, updates)

    def __call__(self, zs, activation=F.tanh, is_sigmoid=False):
        h = self.l1(zs)
        h = activation(h)

        h = self.l2(h)
        h = activation(h)

        h = self.l3(h)
        h = activation(h)

        h = self.l4(h)
        if is_sigmoid:
            h = F.sigmoid(h)
        return h
  1. 入力値は潜在変数zsである。
  2. 最終層にだけ仕掛けがしてある。訓練時はis_sigmoid=False、テスト時にはis_sigmoid=Trueとする。

Discriminator


 次はDiscriminator($T_{\psi}(x,z)$)である(discriminator.py)。
class Discriminator_1(chainer.Chain):

    def __init__(self, x_dim, z_dim, h_dim=512):
        super(Discriminator_1, self).__init__()
        self.h_dim = h_dim
        with self.init_scope():
            self.xl1 = L.Linear(x_dim, h_dim, initialW=xavier.Xavier(x_dim, h_dim))
            self.xl2 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.xl3 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.zl1 = L.Linear(z_dim, h_dim, initialW=xavier.Xavier(z_dim, h_dim))
            self.zl2 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.zl3 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))

    def update(self, updates):
        update_links(self, updates)

    def __call__(self, xs, zs, activation=F.relu):
        xs = 2 * xs - 1 
        hx = self.xl1(xs)
        hx = activation(hx)
        hx = self.xl2(hx)
        hx = activation(hx)
        hx = self.xl3(hx)
        hx = activation(hx)

        hz = self.zl1(zs)
        hz = activation(hz)
        hz = self.zl2(hz)
        hz = activation(hz)
        hz = self.zl3(hz)
        hz = activation(hz)
        h = F.sum(hx * hz, axis=1) / self.h_dim
        return h
  1. 入力値は観測値xsと潜在変数zsである。
  2. 18行目:xsは0と1である。これを-1と1に置き換える。
  3. 19行目から24行目:xsについての計算である。
  4. 26行目から31行目:zsについての計算である。
  5. 32行目:それぞれの結果の内積を取る。

$(\theta,\phi)$についての目的関数(式(\ref{eq8}))


 式(\ref{eq8})を実装したものが以下である(phi_loss_calculator.py)。
class PhiLossCalculator_2(chainer.Chain):
    
    def __init__(self, encoder, decoder, discriminator):
        super(PhiLossCalculator_2, self).__init__()
        with self.init_scope():
            self.encoder = encoder
            self.decoder = decoder
            self.discriminator = discriminator

    def __call__(self, xs, zs, es):
        batch_size = xs.shape[0]
        encoded_zs = self.encoder(xs, es)
        ys = self.decoder(encoded_zs)
        d_loss = F.bernoulli_nll(xs, ys) / batch_size
        t_loss = F.sum(self.discriminator(xs, encoded_zs)) / batch_size
        return t_loss + d_loss, encoded_zs
  1. 式(\ref{eq8})では最大化しているが、実装では符号を反転させて最小化させる。

$\psi$についての目的関数(式(\ref{eq9}))


 式(\ref{eq9})を実装したものが以下である(psi_loss_calculator.py)。
class PsiLossCalculator_3(chainer.Chain):

    def __init__(self, encoder, discriminator):
        super(PsiLossCalculator_3, self).__init__()
        with self.init_scope():
            self.encoder = encoder
            self.discriminator = discriminator

    def __call__(self, xs, zs, es):
        # batch_size = xs.shape[0]
        encoded_zs = self.encoder(xs, es) 
        posterior = self.discriminator(xs, encoded_zs)
        prior = self.discriminator(xs, zs) 
        a = F.sigmoid_cross_entropy(posterior, np.ones_like(posterior).astype(np.int32))
        b = F.sigmoid_cross_entropy(prior, np.zeros_like(prior).astype(np.int32))
        c = F.sum(a + b)
        return c
  1. こちらの目的関数も符号を反転させて最小化させる。

訓練コード


 訓練コードの一部を掲載する。
    for epoch in range(args.epochs):
        with chainer.using_config('train', True):
            # shuffle dataset
            sampler.shuffle_xs()

            epoch_phi_loss = 0 
            epoch_psi_loss = 0 

            for i in range(batches):
                xs = sampler.sample_xs()
                zs = sampler.sample_zs()
                es = sampler.sample_es()

                # compute psi-gradient(eq.3.3)
                update_switch.update_models(enc_updates=False, dec_updates=False,
                                            dis_updates=True)
                psi_loss = psi_loss_calculator(xs, zs, es) 
                update(psi_loss, psi_loss_calculator, psi_optimizer)
                epoch_psi_loss += psi_loss

                # compute phi-gradient(eq.3.7)
                update_switch.update_models(enc_updates=True, dec_updates=True,
                                            dis_updates=False)
                phi_loss, _ = phi_loss_calculator(xs, zs, es) 
                update(phi_loss, phi_loss_calculator, phi_optimizer)
                epoch_phi_loss += phi_loss

            # end for ...
            # see loss per epoch
            epoch_phi_loss /= batches
            epoch_psi_loss /= batches

        # end with ...
        print('epoch:{}, phi_loss:{}, psi_loss:{}'.format(epoch, epoch_phi_loss.data,
                                                          epoch_psi_loss.data))
        epoch_phi_losses.append(epoch_phi_loss.data)
        epoch_psi_losses.append(epoch_psi_loss.data)

    # end for ...
  1. 文献に掲載されているアルゴリズム図をほぼ踏襲している。
  2. 15行目から19行目:目的関数(\ref{eq9})を更新する。その際、EncoderとDecoderの重みを固定する(15行目)。
  3. 22行目から26行目:目的関数(\ref{eq8})を更新する。その際、Discriminatorの重みを固定する(22行目)。

予測のためのコード


 訓練後に実行するコートはpredict.pyである。掲載は略。

実行


 次を実行する。
$> ./train.py
$> ./predict.py


結果


 結果の表示はvisualize.ipynbで行った。最初に式(\ref{eq9})と(\ref{eq8})の変化を示す。
 以下が学習後の$z_{\phi}(x,\epsilon)$の結果である。
Lossの振る舞いは良くわからないが、潜在変数$z$の分布は$(0,0)$を中心に綺麗に4分割されている。文献掲載の結果をそれなりに再現できているように見えるが。。。

まとめ


 論文Adversarial Variational Bayesの前半部分の解説と実装を示した。正直なところ、上のような形に到達するまでかなりの時間を費やして試行錯誤を繰り返した。2つの目的関数を交互に最適化するのは大変難しい。ご批判いただければ幸いである。

2018年5月2日水曜日

Variational Auto Encoder 〜その3〜

はじめに


 先のページで、Variational Auto Encoder(VAE)を実装したChainerのサンプルコードをそのまま動かし、生成される画像を見た。符号化・復号化して得られる画像は、入力画像をそれなりに再現していたが、乱数から生成される画像は精度が良くないことを示した。今回は、サンプルコードに手を加えて実験を行う。

ソースコードの場所


 今回のソースコードはここにある。

評価基準


 このページで示したように、勾配降下法で最適化すべき式は次式である。 \begin{equation} \min_{\phi} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}|X) \right] = \min_{\phi} {\left[ D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]-E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right] \right] } \end{equation} この右辺第1項は次のように書けた。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]= \frac{1}{2} \sum_{d=1}^{D}\left\{ -\ln{\sigma^2_{\phi,d}(X)}-1+\sigma^2_{\phi,d}(X)+\mu_{\phi,d}^2(X) \right\} \label{eq1} \end{equation} また、潜在変数$\vec{z}$の成分は次式で与えられた。 \begin{equation} z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d \end{equation} 式(\ref{eq1})で理想的な最適化を実現できれば、$\mu_{\phi,d}(X)\rightarrow 0,\sigma_{\phi,d}(X)\rightarrow 1$とできる。このとき、$\vec{z}$は標準正規分布から生成される値$\vec{\epsilon}$に置き換えることができる。net.py内にある関数decodeは本来ならば$\vec{z}$を与えて画像を生成するコードであるが、実際に使われる際は標準正規分布から作られる値を与えている。従って、上記の収束が不十分であると、意図した振る舞いをしないことになる。最適化を行う際に、$\mu_{\phi,d}(X)$と$\sigma_{\phi,d}(X)$の値の変化も追跡すべきである。以上の考察に基づきnet.pyに実装されている関数get_loss_funcを以下のように変更した。
    def get_loss_func(self, C=1.0, k=1):
        """Get loss function of VAE.
        The loss value is equal to ELBO (Evidence Lower Bound)
        multiplied by -1.
        Args:
            C (int): Usually this is 1.0. Can be changed to control the
                second term of ELBO bound, which works as regularization.
            k (int): Number of Monte Carlo samples used in encoded vector.
        """
        def lf(x):
            mu, ln_var = self.encode(x)
            mean_mu, mean_sigma = calculate_means(mu, ln_var)
            batchsize = len(mu.data)
            # reconstruction loss
            rec_loss = 0
            for l in six.moves.range(k):
                z = F.gaussian(mu, ln_var)
                rec_loss += F.bernoulli_nll(x, self.decode(z, sigmoid=False)) \
                    / (k * batchsize)
            self.rec_loss = rec_loss
            kl = gaussian_kl_divergence(mu, ln_var) / batchsize
            self.loss = self.rec_loss + C * kl
            chainer.report(
                {
                    'rec_loss': rec_loss,
                    'loss': self.loss,
                    'kl': kl,
                    'mu': mean_mu,
                    'sigma': mean_sigma,
                },
                observer=self)
            return self.loss
        return lf
12行目でmuとsigmaの平均値を計算している。これに伴い、実行中にコマンドラインに表示する項目を増やしている( 25行目から29行目)。関数calculate_meansの中身は以下の通りである。
def calculate_means(mu, ln_var):
    xp = chainer.cuda.get_array_module(mu)
    mean_mu = xp.mean(mu.data)
    sigma = xp.exp(ln_var.data / 2)
    mean_sigma = xp.mean(sigma)
    return mean_mu, mean_sigma

実験-1


 先のページと同じパラメータ(epoch=100,dimz=20)で実行したときのmuとsigmaの変化は以下の通りである。
mu

sigma

訓練データに対するmuの値は0近傍を推移しているがepochの後半部分で振動しており、テストデータに対するそれは終始振動していることが分かる。sigmaの方は1からは程遠い値である。このときのトータルのlossは以下のようになる。
loss


実験-2


 次に、こちらで紹介されている実装で実験を行った。ソースコード内にnet_2.pyとして保存してある。epoch=100、dimz=100とした結果を以下に示す。
mu

sigma

muもsigmaもかなり改善されていることが分かる。muは0に、sigmaは1に近づいている。このときのトータルのlossは以下のようになる。

こちらも少しだけ改善されている。

生成画像の比較


 実験-1と2で、標準正規分布から生成した同じ乱数値から復号された画像を示す。
実験-1

実験-2

後者の方が数字であることを判別できるので精度は良くなっている(と思う)。

まとめ


 今回は、Chainerのサンプルコードに手を加えて実験を行った。muは0に、sigmaは1に近い方が復号化される画像の精度も上がることを示した。VAEを使用する場合、lossだけでなく、muやsigmaの値の推移も確認した方が良い。

2018年4月27日金曜日

Conditional Variational Auto Encoder

はじめに


 先の2回の投稿(ここここ)では、Variational Auto Encoder(VAE)をBayes推論の枠組みで解説した。今回は、Conditional Variational Auto Encoder(CAVE)をBayes推論の枠組みで説明する。

問題設定


 回帰問題を考え、$N$個のペア$(\vec{x}_n, \vec{y}_n)$が観測されているとする。$X=\{\vec{x}_1,\cdots,\vec{x}_N\}, Y=\{\vec{y}_1,\cdots,\vec{y}_N\}$と置いたとき、未観測データ$\vec{x}_{\alpha}$に対応する$\vec{y}_{\alpha}$を生成する確率分布$p(\vec{y}_{\alpha}|\vec{x}_{\alpha},X,Y)$を求めたい。潜在変数$\vec{z}$を導入し、$X$、$Y$、$\vec{z}$の同時確率分布$p(X,Y,\vec{z})$を考え、Bayesの定理を適用すると次式を得る。 \begin{equation} p(\vec{z}|X,Y) = \frac{p(Y|X,\vec{z})p(\vec{z})}{p(Y|X)} \label{eq9} \end{equation} ただし、式変形の途中で$p(X|\vec{z})=p(X)$を用いた。事後確率$p(\vec{z}|X,Y)$が求まれば、次式により$\vec{y}_{\alpha}$を生成する確率分布を求めることができる。 \begin{equation} p(\vec{y}_{\alpha}|\vec{x}_{\alpha},X,Y)=\int d\vec{z}\;p(\vec{y}_{\alpha}|\vec{x}_{\alpha},\vec{z})p(\vec{z}|X,Y) \label{eq3} \end{equation} 事後確率$p(\vec{z}|X,Y)$を求めることが目的である。

最適化すべき量


 $p(\vec{z}|X,Y)$を直接求めることはせず、パラメータ$\phi$を持つ関数$q_{\phi}(\vec{z}|X,Y)$を導入し、次のKullback Leibler divergenceを最小にすることを考える。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X,Y)||p(\vec{z}|X,Y) \right]=\int d\vec{z}\;q_{\phi}(\vec{z}|X,Y) \ln{ \frac{ q_{\phi}(\vec{z}|X,Y) } { p(\vec{z}|X,Y) } } \end{equation} これを変形すると次式を得る。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X,Y)||p(\vec{z}|X,Y) \right] = D_{KL} \left[ q_{\phi}(\vec{z}|X,Y)||p(\vec{z}) \right]-E_{q_{\phi}(\vec{z}|X,Y)}\left[\ln{p(Y|X,\vec{z})}\right]+\ln{p(Y|X)} \label{eq1} \end{equation} ただし、式変形の途中で式(\ref{eq9})を用いた。式(\ref{eq1})右辺にある$\ln{p(Y|X)}$は$\phi$に依存せず、観測値だけから決まる定数である。従って、次式が成り立つ。 \begin{equation} \min_{\phi} D_{KL} \left[ q_{\phi}(\vec{z}|X,Y)||p(\vec{z}|X,Y) \right] = \min_{\phi} {\left[ D_{KL} \left[ q_{\phi}(\vec{z}|X,Y)||p(\vec{z}) \right]-E_{q_{\phi}(\vec{z}|X,Y)}\left[\ln{p(Y|X,\vec{z})}\right] \right] } \label{eq2} \end{equation} 式(\ref{eq2})の右辺第1項を小さく、第2項の期待値を大きくすれば良い。第1項は$q_{\phi}(\vec{z}|X,Y)$をできるだけ$p(\vec{z})$に近い形の分布にすることを要求し、この分布の下で対数尤度$\ln{p(Y|X,\vec{z})}$の期待値を大きくすることを第2項は要求する。第1項は正則化項に相当する。

KL divergenceの計算


 式(\ref{eq2})の右辺第1項を考える。いま次の仮定をおく。 \begin{eqnarray} q_{\phi}(\vec{z}|X,Y)&=&\mathcal{N}(\vec{z}|\vec{\mu}_{\phi}(X,Y),\Sigma_{\phi}(X,Y)) \\ p(\vec{z})&=&\mathcal{N}(\vec{z}|\vec{0},I_D) \end{eqnarray} ここで、$\vec{z}$の次元を$D$とした。$I_D$は$D\times D$の単位行列である。どちらの分布も正規分布とし、$q_{\phi}(\vec{z}|X,Y)$の平均と共分散行列は$\phi,X,Y$から決まる量とする。これらは、入力$X,Y$、パラメータ$\phi$のニューラルネットワークを用いて計算される。一方、$p(\vec{z})$の方は平均0、分散1の標準正規分布である。このとき、$D_{KL} \left[ q_{\phi}(\vec{z}|X,Y)||p(\vec{z}) \right]$は解析的に計算することができる。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X,Y)||p(\vec{z}) \right]=\frac{1}{2}\left[ -\ln{|\Sigma_{\phi}(X,Y)|} -D +\mathrm{Tr}\left(\Sigma_{\phi}(X,Y)\right)+\vec{\mu}_{\phi}(X,Y)^T\vec{\mu}_{\phi}(X,Y) \right] \label{eq4} \end{equation}

ここまでの処理の流れ


 式(\ref{eq2})を計算する際の手順は以下のようになる。
分布$q_{\phi}(\vec{z}|X,Y)$は$X$と$Y$から$\vec{z}$を生成するEncoder、$p(Y|X,\vec{z})$は$X$と$\vec{z}$から$Y$を生成するDecoderとみなすことができる。青色で示した部分は最小化すべき量であり、赤字はサンプリングするステップである。青色の式の和を勾配降下法により最小にするが、その際、誤差逆伝播ができなければならない。$q_{\phi}(\vec{z}|X,Y)$はその$\phi$依存性のため誤差逆伝播時の微分鎖の中に組み込まれるが、サンプリングという処理の勾配を定義することができない。対数尤度の期待値の計算に工夫が必要である。

対数尤度の期待値の計算


 計算したい式は次式である。 \begin{equation} E_{q_{\phi}(\vec{z}|X,Y)}\left[\ln{p(Y|X,\vec{z})}\right]=\int d\vec{z}\;q_{\phi}(\vec{z}|X,Y)\ln{p(Y|X,\vec{z})} \end{equation} この式に再パラメータ化トリック(re-parametrization trick)を適用する。すなわち \begin{equation} \vec{z}\sim\mathcal{N}(\vec{z}|\vec{\mu}_{\phi}(X,Y),\Sigma_{\phi}(X,Y)) \end{equation} の代わりに \begin{eqnarray} \vec{\epsilon}&\sim&\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\\ \vec{z}&=& \vec{\mu}_{\phi}(X,Y)+\Sigma_{\phi}^{1/2}(X,Y)\vec{\epsilon} \label{eq7} \end{eqnarray} を用いてサンプリングを行う。これを用いて期待値を書き直すと次式を得る。 \begin{equation} E_{q_{\phi}(\vec{z}|X,Y)}\left[\ln{p(Y|X,\vec{z})}\right]=\int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\ln{p(Y|X, \vec{z}=\vec{\mu}_{\phi}(X,Y)+\Sigma_{\phi}^{1/2}(X,Y)\vec{\epsilon})} \label{eq11} \end{equation} 処理の流れは以下のように変更される。
上図であれば誤差逆伝播が可能となる。

未観測データの生成


 未観測データ$\vec{y}_{\alpha}$を生成する確率分布は次式で与えられた。 \begin{equation} p(\vec{y}_{\alpha}|\vec{x}_{\alpha},X,Y)=\int d\vec{z}\;p(\vec{y}_{\alpha}|\vec{x}_{\alpha},\vec{z})p(\vec{z}|X,Y) \end{equation} 事後確率$p(\vec{z}|X,Y)$の近似解$q_{\phi}(\vec{z}|X,Y)$を用いると \begin{equation} p(\vec{y}_{\alpha}|\vec{x}_{\alpha},X,Y)\approx\int d\vec{z}q_{\phi}(\vec{z}|X,Y)p(\vec{y}_{\alpha}|\vec{x}_{\alpha},\vec{z}) \end{equation} を得る。先と同様に再パラメータ化トリックを適用すると \begin{equation} p(\vec{y}_{\alpha}|\vec{x}_{\alpha},X,Y)\approx\int d\vec{\epsilon}\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)p(\vec{y}_{\alpha}|\vec{x}_{\alpha}, \vec{z}=\vec{\mu}_{\phi}(X,Y)+\Sigma_{\phi}^{1/2}(X,Y)\vec{\epsilon}) \label{eq10} \end{equation} となる。

実装に向けた詳細な計算


 最初に$\vec{\mu}_{\phi}(X,Y)$と$\Sigma_{\phi}(X,Y)$を次のように置く。 \begin{eqnarray} \vec{\mu}_{\phi}(X,Y)&=&(\mu_{\phi,1}(X,Y),\cdots,\mu_{\phi,D}(X,Y))^T \label{eq5}\\ \Sigma_{\phi}(X,Y)&=&\mathrm{diag}(\sigma^2_{\phi,1}(X,Y),\cdots,\sigma^2_{\phi,D}(X,Y)) \label{eq6} \end{eqnarray} このとき式(\ref{eq4})は次式となる。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X,Y)||p(\vec{z}) \right]= \frac{1}{2} \sum_{d=1}^{D}\left\{ -\ln{\sigma^2_{\phi,d}(X,Y)}-1+\sigma^2_{\phi,d}(X,Y)+\mu_{\phi,d}^2(X,Y) \right\} \label{eq8} \end{equation} また、$\vec{z}$の成分は次式で与えられる。 \begin{equation} z_d=\mu_{\phi,d}(X,Y)+\sigma_{\phi,d}(X,Y)\epsilon_d \end{equation} 観測値が独立同分布に従うと仮定すると、式(\ref{eq11})の中にある対数尤度は以下のように変形される。 \begin{equation} \ln{p(Y|X,\vec{z})}= \sum^{N}_{n=1}\ln{p(\vec{y}_n|\vec{x}_n,\vec{z})} \end{equation} さらに計算を進めるには、具体的に$X,Y$として、何を与えるか決定しなければならない。 ここでは、$X$として0から9までのラベルを、$Y$としてMNISTの画像(2値画像)を与えることにする。$X$の各観測値$\vec{x}_n$は$9$次元のone-hotベクトルで表現される。各画素が独立同分布に従うと仮定すると、$\vec{y}_n$の次元数を$M$として \begin{equation} \ln{p(\vec{y}_n|\vec{x}_n,\vec{z})}=\sum_{m=1}^{M}\ln{p(y_{n,m}|\vec{x}_n,\vec{z})} \end{equation} と書くことができる。いま考える画像は0と1から構成されるから、$p(y_{n,m}|\vec{x}_n,\vec{z})$として0と1を生成するBernoulli分布を仮定する。 \begin{eqnarray} p(y_{n,m}|\vec{x}_n,\vec{z})&=&\mathrm{Bern}\left(y_{n,m}|\eta_{\theta,m}\left(\vec{x}_n,\vec{z}\right)\right) \\ \mathrm{Bern}(x|\eta)&=&\eta^{x}(1-\eta)^{1-x} \end{eqnarray} $\eta_{\theta,m}\left(\vec{x}_n,\vec{z}\right)$は、入力が$\vec{x}_n$と$\vec{z}$、パラメータとして$\theta$を持つニューラルネットワークで学習される。以上を踏まえて処理の流れを書き直すと下図となる。

次に式(\ref{eq10})を考える。これは、観測値$X,Y$とラベル$\vec{x}_{\alpha}$が与えられときの$\vec{y}_{\alpha}$の実現確率である。 \begin{equation} p(\vec{y}_{\alpha}|\vec{x}_{\alpha},X,Y)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;p(\vec{y}_{\alpha}|\vec{x}_{\alpha},\vec{z}) \end{equation} ここで、$z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d$である。上式は以下のように書くことができる。 \begin{equation} \prod_{m=1}^M p(y_{\alpha,m}|\vec{x}_\alpha,X,Y)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\prod_{m=1}^M p(y_{\alpha,m}|\vec{x}_\alpha,\vec{z}) \end{equation} すなわち、要素$y_{\alpha,m}$ごとに次式が成り立つ。 \begin{equation} p(y_{\alpha,m}|\vec{x}_\alpha,X,Y)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;p(y_{\alpha,m}|\vec{x}_\alpha,\vec{z}) \end{equation} $p(y_{\alpha,m}|\vec{x}_\alpha,\vec{z})$としてBernoulli分布を仮定したから \begin{equation} p(y_{\alpha,m}|\vec{x}_\alpha,X,Y)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\mathrm{Bern}(y_{\alpha,m}|\eta_{\theta,m}(\vec{x}_\alpha,\vec{z})) \end{equation} となる。確率分布$p(y_{\alpha,m}|\vec{x}_\alpha,X,Y)$の下での$y_{\alpha,m}$の期待値は \begin{eqnarray} <y_{\alpha,m}>&=&\sum_{y_{\alpha,m}=0,1} y_{\alpha,m}\;p(y_{\alpha,m}|\vec{x}_\alpha,X,Y) \\ &\approx& \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\sum_{y_{\alpha,m}=0,1} y_{\alpha,m} \mathrm{Bern}(y_{\alpha,m}|\eta_{\theta,m}(\vec{x}_\alpha,\vec{z})) \\ &=& \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\eta_{\theta,m}(\vec{x}_\alpha,\vec{z}) \end{eqnarray} となる。$\vec{z}$は$\vec{\epsilon}$に依存する項であることに注意する。$\eta_{\theta,m}(\vec{x}_\alpha,\vec{z})$はDecoderの出力である。上式から、復号化した結果を得るには、$\eta_{\theta,m}(\vec{x}_\alpha,\vec{z})$を標準正規分布に従ってサンプリングすれば良いことが分かる。さらに、式(\ref{eq8})のKullback Leibler divergenceを十分小さくできれば、すなわち、$\sigma_{\phi,d}(X,Y)\rightarrow 1,\mu_{\phi,d}(X,Y)\rightarrow 0$とできれば、$\vec{z}=\vec{\epsilon}$とすることができるので、標準正規分布から生成した値$\vec{\epsilon}$と$\vec{x}_\alpha$からDecoderの出力を直接得ることができる。

まとめ


 今回は、CVAEをBayes推定の枠組みで説明した。前回のVAEの論法とほとんど同じである。VAEでは未観測データ$\vec{x}$が従う確率分布$p(\vec{x}|X)$を求める過程でVAEの構造を見出した。一方、CVAEでは未観測データ$\vec{x}$に対応する$\vec{y}$が従う確率分布$p(\vec{y}|\vec{x},X,Y)$を求める過程でCVAEの構造が現れることを見た。その構造は、VAEに少し手を加えれば実現できる程度のものである。ChainerのVAEのサンプルコードをベースにすればすぐに実装できそうである。

2018年4月22日日曜日

Variational Auto Encoder 〜その2〜

はじめに


 先のページで、Variational Auto Encoder(VAE)をBayes推論の枠組みで解説し、Chainerのサンプルコードをみた。今回は、サンプルコードを実際に動かし、その動作を確認する。

前回の補足


 前回の式(16)は、観測値$X$の下での未知変数$\vec{x}$の実現確率であった。 \begin{equation} p(\vec{x}|X)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;p(\vec{x}|\vec{z}) \end{equation} ここで、$z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d$である。$\vec{x}$の各成分が独立同分布で生成されると仮定すると、以下のように書き換えることができる。 \begin{equation} \prod_{m=1}^M p(x_{m}|X)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\prod_{m=1}^M p(x_m|\vec{z}) \end{equation} すなわち、要素$x_m$ごとに次式が成り立つ。 \begin{equation} p(x_{m}|X)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;p(x_m|\vec{z}) \end{equation} $p(x_m|\vec{z})$としてBernoulli分布を仮定したから \begin{equation} p(x_{m}|X)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\mathrm{Bern}(x_m|\eta_{\theta,m}(\vec{z})) \end{equation} となる。確率分布$p(x_{m}|X)$の下での$x_m$の期待値は \begin{eqnarray} <x_m>&=&\sum_{x_m=0,1} x_m\;p(x_{m}|X) \\ &\approx& \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\sum_{x_m=0,1} x_m \mathrm{Bern}(x_m|\eta_{\theta,m}(\vec{z})) \\ &=& \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\eta_{\theta,m}(\vec{z}) \end{eqnarray} となる。$\vec{z}$は$\vec{\epsilon}$に依存する項であることに注意する。ここまでの議論で言えることは、復号化した結果を得るには、$\eta_{\theta,m}(\vec{z})$を標準正規分布に従ってサンプリングすれば良いということである。$\eta_{\theta,m}(\vec{z})$はDecoderの出力である。

Chainerの実装の確認


 前回は、net.pyを見たので、今回はtrain_vae.pyを見る。trainerを用いた実装部分は特に指摘することはないので、結果を描画している箇所を解説する。最初は訓練データに関わる描画部分である。 適当に画像を9枚選択し、これを関数__call__で符号化・復号化している。前回指摘したように、復号化の際、$\sigma_{\phi,d}(X)$は無視されている。計算のあと元画像と復号化画像を保存している。epochを100とした結果は以下の通りである。

訓練画像


復号化した訓練画像

次はテスト画像に関わる部分である。 ここでも適当に9枚の画像を選択し、元画像と復号化画像を保存している。epochを100とした結果は以下の通りである。

テスト画像


復号化したテスト画像

最後に、標準正規分布に従う値$\vec{z}$から復号化する部分である。 9個の乱数$\vec{z}$を作り、関数decodeを呼び出している。epochを100とした結果は以下の通りである。
訓練・テスト画像を符号化・復号化した結果と比べるとかなり精度の悪い結果である。関数decodeは$\eta_{\theta,m}(\vec{z})$を出力する。上で見たように本来の$\vec{z}$は、$z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d$として与えらるべきものである。精度が悪いのは、$\vec{z}$を標準正規分布に置き換えたことが原因である。ところで、勾配降下法で最小にすべき式の1つが次式であった(前回の式(19))。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]= \frac{1}{2} \sum_{d=1}^{D}\left\{ -\ln{\sigma^2_{\phi,d}(X)}-1+\sigma^2_{\phi,d}(X)+\mu_{\phi,d}^2(X) \right\} \end{equation} これを十分小さくできれば、すなわち、$\sigma_{\phi,d}\rightarrow 1, \mu_{\phi,d}\rightarrow 0$とできれば、標準正規分布による置き換えは意味があるものになる。残念ながらepochを1000としても大して精度は良くならない。Chainerのサンプル実装を変更する必要があるかもしれない。

まとめ


 今回は、$<x_m>$が$\eta_{\theta,m}(\vec{z})$を求めることに帰着すること示し、Chainerのサンプルコードの計算結果を考察した。 さらに、標準正規分布による$\vec{z}$から計算した画像の精度が悪い理由についても述べた。改善するには、epoch数を増やすだけではなくコードの見直し(多層化、初期化関数の変更)も必要であろう。次回はこの辺りのことについてまとめたい。

2018年4月15日日曜日

Variational Auto Encoder

はじめに


 Variational Auto Encoder(VAE)をBayes推論の枠組みで解説し、Chainerのサンプルコードを読解する。

問題設定


 観測値$X=\{\vec{x}_1,\cdots,\vec{x}_N\}$が与えられたとき、未知の値$\vec{x}_*$を生成する確率分布$p(\vec{x}_*|X)$を求めたい。潜在変数$\vec{z}$を導入し、$X$と$\vec{z}$の同時確率分布$p(X,\vec{z})$を考え、Bayesの定理を適用すると次式を得る。 \begin{equation} p(\vec{z}|X) = \frac{p(X|\vec{z})p(\vec{z})}{p(X)} \label{eq9} \end{equation} 事後確率$p(\vec{z}|X)$が求まれば、次式により$\vec{x}_*$を生成する確率分布を求めることができる。 \begin{equation} p(\vec{x}_*|X)=\int d\vec{z}\;p(\vec{x}_*|\vec{z})p(\vec{z}|X) \label{eq3} \end{equation} 事後確率$p(\vec{z}|X)$を求めることが目的である。

最適化すべき量


 $p(\vec{z}|X)$を直接求めることはせず、パラメータ$\phi$を持つ関数$q_{\phi}(\vec{z}|X)$を導入し、次のKullback Leibler divergenceを最小にすることを考える。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}|X) \right]=\int d\vec{z}\;q_{\phi}(\vec{z}|X) \ln{ \frac{ q_{\phi}(\vec{z}|X) } { p(\vec{z}|X) } } \end{equation} これを変形すると次式を得る。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}|X) \right] = D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]-E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right]+\ln{p(X)} \label{eq1} \end{equation} ただし、式変形の途中で式(\ref{eq9})を用いた。式(\ref{eq1})右辺にある$\ln{p(X)}$は$\phi$に依存せず、観測値だけから決まる定数である。従って、次式が成り立つ。 \begin{equation} \min_{\phi} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}|X) \right] = \min_{\phi} {\left[ D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]-E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right] \right] } \label{eq2} \end{equation} 式(\ref{eq2})の右辺第1項を小さく、第2項の期待値を大きくすれば良い。第1項は$q_{\phi}(\vec{z}|X)$をできるだけ$p(\vec{z})$に近い形の分布にすることを要求し、この分布の下で対数尤度$\ln{p(X|\vec{z})}$の期待値を大きくすることを第2項は要求する。第1項は正則化項に相当する。

KL divergenceの計算


 式(\ref{eq2})の右辺第1項を考える。いま次の仮定をおく。 \begin{eqnarray} q_{\phi}(\vec{z}|X)&=&\mathcal{N}(\vec{z}|\vec{\mu}_{\phi}(X),\Sigma_{\phi}(X)) \\ p(\vec{z})&=&\mathcal{N}(\vec{z}|\vec{0},I_D) \end{eqnarray} ここで、$\vec{z}$の次元を$D$とした。$I_D$は$D\times D$の単位行列である。どちらの分布も正規分布とし、$q_{\phi}(\vec{z}|X)$の平均と共分散行列は$\phi$と$X$から決まる量とする。これらは、入力$X$、パラメータ$\phi$のニューラルネットワークを用いて計算される。一方、$p(\vec{z})$の方は平均0、分散1の標準正規分布である。このとき、$D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]$は解析的に計算することができる。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]=\frac{1}{2}\left[ -\ln{|\Sigma_{\phi}(X)|} -D +\mathrm{Tr}\left(\Sigma_{\phi}(X)\right)+\vec{\mu}_{\phi}(X)^T\vec{\mu}_{\phi}(X) \right] \label{eq4} \end{equation}

ここまでの処理の流れ


 式(\ref{eq2})の最適化を行う際に勾配降下法を用いる。処理の流れは以下のようになる(下図参照)。
分布$q_{\phi}(\vec{z}|X)$は$X$から$\vec{z}$を生成するEncoder、$p(X|\vec{z})$は$\vec{z}$から$X$を生成するDecoderとみなすことができる。青色で示した部分は最小化すべき量であり、赤字はサンプリングするステップである。勾配降下法を実現するには、誤差逆伝播ができなければならない。$q_{\phi}(\vec{z}|X)$はその$\phi$依存性のため誤差逆伝播時の微分鎖の中に組み込まれるが、サンプリグという処理の勾配を定義することができない。対数尤度の期待値の計算に工夫が必要である。

対数尤度の期待値の計算


 計算したい式は次式である。 \begin{equation} E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right]=\int d\vec{z}\;q_{\phi}(\vec{z}|X)\ln{p(X|\vec{z})} \end{equation} この式に再パラメータ化トリック(re-parametrization trick)を適用する。すなわち \begin{equation} \vec{z}\sim\mathcal{N}(\vec{z}|\vec{\mu}_{\phi}(X),\Sigma_{\phi}(X)) \end{equation} の代わりに \begin{eqnarray} \vec{\epsilon}&\sim&\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\\ \vec{z}&=& \vec{\mu}_{\phi}(X)+\Sigma_{\phi}^{1/2}(X)\vec{\epsilon} \label{eq7} \end{eqnarray} を用いてサンプリングを行う。これを用いて期待値を書き直すと次式を得る。 \begin{equation} E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right]=\int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\ln{p(X| \vec{z}=\vec{\mu}_{\phi}(X)+\Sigma_{\phi}^{1/2}(X)\vec{\epsilon})} \end{equation} このときの処理の流れは以下のようになる(下図参照)。
上図であれば誤差逆伝播が可能となる。

未知変数の生成


 未知変数を生成する確率分布は次式で与えられた。 \begin{equation} p(\vec{x}_*|X)=\int d\vec{z}\;p(\vec{x}_*|\vec{z})p(\vec{z}|X) \end{equation} 事後確率$p(\vec{z}|X)$の近似解$q_{\phi}(\vec{z}|X)$を用いると \begin{equation} p(\vec{x}_*|X)\approx\int d\vec{z}q_{\phi}(\vec{z}|X)p(\vec{x}_*|\vec{z}) \end{equation} を得る。先と同様に再パラメータ化トリックを適用すると \begin{equation} p(\vec{x}_*|X)\approx\int d\vec{\epsilon}\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)p(\vec{x}_*| \vec{z}=\vec{\mu}_{\phi}(X)+\Sigma_{\phi}^{1/2}(X)\vec{\epsilon}) \end{equation} となる。

Chainer実装の確認


 ChainerのサンプルコードにVAEがある。これはMNISTデータセットにVAEを適用したものである。MNISTは2値画像であるから$\vec{x}_n$として0と1が784(=28$\times$28)個並んだベクトルを考えることになる。実際にコードを見て行く前に先に導出した式をもう少し詳細に計算しておく。

 最初に$\vec{\mu}_{\phi}(X)$と$\Sigma_{\phi}(X)$を次のように置く。 \begin{eqnarray} \vec{\mu}_{\phi}(X)&=&(\mu_{\phi,1}(X),\cdots,\mu_{\phi,D}(X))^T \label{eq5}\\ \Sigma_{\phi}(X)&=&\mathrm{diag}(\sigma^2_{\phi,1}(X),\cdots,\sigma^2_{\phi,D}(X)) \label{eq6} \end{eqnarray} このとき式(\ref{eq4})は次式となる。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]= \frac{1}{2} \sum_{d=1}^{D}\left\{ -\ln{\sigma^2_{\phi,d}(X)}-1+\sigma^2_{\phi,d}(X)+\mu_{\phi,d}^2(X) \right\} \label{eq8} \end{equation} また、$\vec{z}$の成分は次式で与えられる。 \begin{equation} z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d \end{equation} 観測値が独立同分布に従うと仮定すると対数尤度は以下のように変形される。 \begin{equation} \ln{p(X|\vec{z})}= \sum^{N}_{n=1}\ln{p(\vec{x}_n|\vec{z})} \end{equation} $\vec{x}_n$の次元数を$M$(=784)とすると \begin{equation} \ln{p(\vec{x}_n|\vec{z})}=\sum_{m=1}^{M}\ln{p(x_{n,m}|\vec{z})} \end{equation} となる。いま考える画像は0と1から構成される。従って、$p(x_{n,m}|\vec{z})$として0と1を生成する確率分布であるBernoulli分布を仮定する。 \begin{eqnarray} p(x_{n,m}|\vec{z})&=&\mathrm{Bern}\left(x_{n,m}|\eta_{\theta,m}\left(\vec{z}\right)\right) \\ \mathrm{Bern}(x|\eta)&=&\eta^{x}(1-\eta)^{1-x} \end{eqnarray} $\eta_{\theta,m}\left(\vec{z}\right)$は、入力$\vec{z}$、パラメータ$\theta$のネットワークで学習される量である。 以上を踏まえて処理の流れを書き直すと下図となる。

以上で準備が整ったので順にコードを見ていく。まず最初にネットワークを定義したクラスVAEを見る。コンストラクタは以下の通り。 Encoderとして$\vec{\mu}_{\phi}(X)$と$\Sigma_{\phi}(X)$を生成する層がそれぞれ1層ずつ定義されている。Decoderとして$\vec{\eta}_{\theta}(\vec{z})$を生成する2層が定義されている。次に関数encodeを見る。 $\mu_{\phi,d}$と$\ln{\sigma_{\phi,d}^2}$を生成する処理が記述されている。次は関数decodeである。 $\vec{z}$を受け取り$\vec{\eta}_{\theta}(\vec{z})$を返す処理が記述されている。次は__call__である。 Encoderで計算した平均を入力としてDecoderを呼び出している。分散を無視して復号化している。次はget_loss_funcである。
  • 13行目:Encoderで平均と対数分散を計算する。
  • 17行目:ここから始まるループは$\vec{z}$のサンプリングのためのものである。
  • 18行目:正規分布からサンプリングする。関数F.gaussianの中で再パラメータ化トリックが実行されていることに注意する。
  • 19行目:Bernoulli分布の対数にマイナスを付けたものが計算される。
  • 23行目:式(\ref{eq8})が計算される。

  • まとめ


     今回は、VAEをBayes推論の枠組みで解説し、Chainerのサンプルコードを見た。ニューラルネットワークで計算されるのは確率分布のパラメータであることを明確に示した。言い換えると他の手法でパラメータを計算できるのであればそれでも構わないということである。計算の仮定で確率分布にいくつかの仮定(ガウス分布やBernoulli分布)をしていることに注意しなければならない。今回Bernoulli分布を使用したのはターゲットとした観測値が0と1から構成される2値画像であるためである。2値でない観測値を対象とするならそれに見合った確率分布を導入する必要がある。次回はChainerのコードを動かして得られる結果を考察したい。

    参考文献


  • Tutorial on Variational Autoencoders
  • Pattern Recognition and Machine Learning
  • ベイズ推論による機械学習入門
  • Variational Autoencoder徹底解説
  • 2018年3月11日日曜日

    強化学習 〜Fisher情報行列と自然勾配法〜

    はじめに


     先のページで方策勾配定理を導いた。すなわち、時間ステップ$t=0$から始める状態価値関数$J(\theta;s_0)$に対して次式が成り立つことを示した。 \begin{equation} \frac{\partial J(\theta;s_0)}{\partial \theta} = {\bf E}\left[\frac{\partial \ln{\pi_{\theta}(a|s)}}{\partial \theta}\;Q(s,a)\right] \end{equation} 今回は、$\theta$をベクトル$\vec{\theta}$に拡張した式 \begin{equation} \vec{\nabla}_{\theta}J(\vec{\theta};s_0) = {\bf E} \left[ \vec{\nabla}_{\theta} \ln{\pi_{\theta}(a|s)}\;Q(s,a) \right] \label{eq2} \end{equation} を考え、$J(\vec{\theta};s_0)$を勾配法で最大化する過程で用いられる手法をまとめる。以降、$J(\vec{\theta};s_0)$を$J(\vec{\theta})$と書く。

    Fisher情報行列の導出


     2つの確率分布$P(w),Q(w)$を考え、これらの間のKullback-Leibler divergenceを考える。 \begin{equation} D(P,Q)=\int dw\; P(w) \ln{\frac{P(w)}{Q(w)}} \end{equation} この量を対称化して、2つの確率分布間の「距離」を表すことができるようにする。 \begin{eqnarray} D_s(P,Q) &=&\int dw\; P(w) \ln{\frac{P(w)}{Q(w)}}+\int dw\; Q(w) \ln{\frac{Q(w)}{P(w)}}\\ &=&\int dw\; \left(P(w)-Q(w)\right) \ln{\frac{P(w)}{Q(w)}} \end{eqnarray} いま、パラメータ$\vec{\theta}$を持つモデルで表現された確率分布$P(w|\vec{\theta})$を考え、パラメータをわずかに変えた後の分布$P(w|\vec{\theta}+d\vec{\theta})$との間の距離$D_s$を計算する。 \begin{eqnarray} D_s(P(w|\vec{\theta}+d\vec{\theta}),P(w|\vec{\theta})) &=&\int dw\; \left(P(w|\vec{\theta}+d\vec{\theta})-P(w|\vec{\theta})\right) \ln{\frac{P(w|\vec{\theta}+d\vec{\theta})}{P(w|\vec{\theta})}} \\ &=&\int dw\; \delta(w|\vec{\theta}) \ln{\frac{\delta(w|\vec{\theta})+P(w|\vec{\theta})}{P(w|\vec{\theta})}} \label{eq0} \end{eqnarray} ここで、 \begin{equation} \delta(w|\vec{\theta})\equiv P(w|\vec{\theta}+d\vec{\theta})-P(w|\vec{\theta}) \end{equation} と置いた。式(\ref{eq0})は$\delta(w|\vec{\theta})$が微小量であることに注意すると以下のように変形することができる。 \begin{eqnarray} D_s(P(w|\vec{\theta}+d\vec{\theta}),P(w|\vec{\theta})) &=&\int dw\; \delta(w|\vec{\theta}) \ln{\frac{\delta(w|\vec{\theta})+P(w|\vec{\theta})}{P(w|\vec{\theta})}} \\ &=&\int dw\; \delta(w|\vec{\theta}) \ln{\left(1+\frac{\delta(w|\vec{\theta})}{P(w|\vec{\theta})}\right)} \\ &\simeq&\int dw\; \delta(w|\vec{\theta}) \frac{\delta(w|\vec{\theta})}{P(w|\vec{\theta})} \label{eq1} \end{eqnarray} ところで、$\delta(w|\vec{\theta})$は次のように計算される。 \begin{eqnarray} \delta(w|\vec{\theta}) &=&P(w|\vec{\theta}+d\vec{\theta})-P(w|\vec{\theta}) \\ &\simeq&\sum_i\;\frac{\partial P(w|\vec{\theta})}{\partial \theta_i} d\theta_i \\ &=&P(w|\vec{\theta})\sum_i\;\frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_i}d\theta_i \end{eqnarray} これを式(\ref{eq1})に代入する。 \begin{eqnarray} D_s(P(w|\vec{\theta}+d\vec{\theta}),P(w|\vec{\theta})) &\simeq& \int dw\; \left[P(w|\vec{\theta})\sum_i\;\frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_i}d\theta_i\right] \frac{P(w|\vec{\theta})\sum_j\;\frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_j}d\theta_j}{P(w|\vec{\theta})} \\ &=& \int dw\; \left[P(w|\vec{\theta})\sum_i\;\frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_i}d\theta_i\right] \sum_j\;\frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_j}d\theta_j \\ &=& \sum_{i,j}\;d\theta_i \left[\int dw\; P(w|\vec{\theta}) \frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_i} \frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_j}\right] d\theta_j \\ &=& \sum_{i,j}\;d\theta_i\; {\bf E}\left[ \frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_i} \frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_j}\right] d\theta_j \end{eqnarray} ここで、${\bf E}[\cdot]$は確率$P(w|\vec{\theta})$の下での期待値を表す。いま \begin{equation} F_{i,j}(\vec{\theta})\equiv {\bf E}\left[ \frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_i} \frac{\partial\ln{ P(w|\vec{\theta})} }{\partial\theta_j}\right] \label{eq4} \end{equation} と置くと \begin{eqnarray} D_s(P(w|\vec{\theta}+d\vec{\theta}),P(w|\vec{\theta})) &\simeq& \sum_{i,j}\;d\theta_i\;F_{i,j}(\vec{\theta})\;d\theta_j\\ &=&d\vec{\theta}^{T}F(\vec{\theta})\;d\vec{\theta} \end{eqnarray} を得る。$F(\vec{\theta})$をFisher情報行列と呼ぶ。$D_s$は微小距離の2乗に相当する量である。従って、$F(\vec{\theta})$は今考えているパラメータ空間(リーマン空間)における計量テンソルと見ることができる。

    自然勾配法


     パラメータ$\vec{\theta}=(\theta_1,\cdots,\theta_n)$の張る空間上の関数$f(\vec{\theta})$を考える。$f(\vec{\theta})$の微小変化は次式で与えられる。 \begin{equation} \delta f\equiv f(\vec{\theta}+d\vec{\theta})-f(\vec{\theta})=\vec{\nabla}f^T d\vec{\theta} \end{equation} $\delta f$が最大となる向き$d\vec{\theta}$を、$\|d\vec{\theta}\|^2=\epsilon^2$の条件の下で考える。ここで、$\epsilon$は微小な定数である。ところで、リーマン空間において 、$\|d\vec{\theta}\|^2$は計量テンソル$G(\vec{\theta})$を用いて、次のように書くことができる。 \begin{equation} \|d\vec{\theta}\|^2=d\vec{\theta}^{T}G(\vec{\theta})\;d\vec{\theta} \end{equation} 従って、$\delta f$が最大となる向きは、Lagrangeの未定乗数法を用いて、次を最大化することで求めることができる。 \begin{equation} L=\vec{\nabla}f^T d\vec{\theta} -\lambda\left(d\vec{\theta}^{T}G(\vec{\theta})\;d\vec{\theta}-\epsilon^2 \right) \end{equation} $d\theta_i$で偏微分して0と置く。 \begin{equation} \frac{\partial L}{\partial d\theta_i}=\frac{\partial f}{\partial \theta_i}-2\lambda\sum_j\;G_{i,j}(\vec{\theta})d\theta_j=0 \end{equation} これより \begin{equation} \vec{\nabla}f = 2\lambda\;G(\vec{\theta})d\vec{\theta} \end{equation} を得る。$d\vec{\theta}$について解くと \begin{equation} d\vec{\theta}\propto G^{-1}(\vec{\theta})\vec{\nabla}f \end{equation} を得る。これが微小変化$\delta f$が最大となる向きである。この向きを利用する勾配法を自然勾配法と呼ぶ。ユークリッド空間の場合は$G$を単位行列と置けば良い。

    方策勾配への適用


     $J(\vec{\theta})$の勾配を最大にする向きは、先の議論より \begin{equation} d\vec{\theta}\propto F^{-1}(\vec{\theta})\vec{\nabla}_{\theta}J(\vec{\theta}) \label{eq5} \end{equation} で与えられる。ただし、ここでの情報行列$F(\vec{\theta})$は式(\ref{eq4})の$P(w|\vec{\theta})$を$\pi_{\theta}(a|s)$に置き換えた次式で定義される。 \begin{equation} F(\vec{\theta})\equiv{\bf E} \left[ \vec{\nabla}_{\theta}\ln{\pi_{\theta}(a|s)} \vec{\nabla}_{\theta}^T\ln{\pi_{\theta}(a|s)} \right] \end{equation} 式(\ref{eq5})に式(\ref{eq2})を代入する。 \begin{equation} d\vec{\theta}\propto F^{-1}(\vec{\theta})\;{\bf E} \left[ \vec{\nabla}_{\theta} \ln{\pi_{\theta}(a|s)}\;Q(s,a) \right] \label{eq3} \end{equation} $Q(s,a)$をパラメータ$\vec{w}$を用いた次のモデルで近似することがよく行われる。 \begin{equation} Q(s,a)=\vec{w}^T \vec{\nabla}_{\theta}\ln{\pi_{\theta}(a|s)} \end{equation} これを、式(\ref{eq3})に代入すると \begin{eqnarray} d\vec{\theta} &\propto& F^{-1}(\vec{\theta})\;{\bf E} \left[ \vec{\nabla}_{\theta} \ln{\pi_{\theta}(a|s)}\;\left[\vec{w}^T \vec{\nabla}_{\theta}\ln{\pi_{\theta}(a|s)}\right] \right] \\ &=& F^{-1}(\vec{\theta})\;{\bf E} \left[ \vec{\nabla}_{\theta}\ln{\pi_{\theta}(a|s)} \vec{\nabla}_{\theta}^T\ln{\pi_{\theta}(a|s)} \right] \vec{w}\\ &=& F^{-1}(\vec{\theta})\;F(\vec{\theta})\;\vec{w}\\ &=& \vec{w} \end{eqnarray} を得る。つまり、$J(\vec{\theta})$を更新する際の$\vec{\theta}$の向きは、$\vec{w}$と一致することになる(というより、こうなるように$Q(s,a)$を近似するのだろう)。

    参考文献


    2018年2月18日日曜日

    強化学習 〜ベルマン最適方程式〜

    はじめに


     先のページでベルマン方程式の導出を行った。この議論に追記を行い、ベルマン最適方程式を示す。

    不動点方程式


     状態価値関数$V(s)$に対し、次のベルマン方程式が成り立つことを見た。 \begin{equation} V(s)=\sum_{s^{\prime},a} P(s^{\prime}|s,a)\pi(a|s)\left[r(s,a,s^{\prime})+\gamma V(s^{\prime})\right] \end{equation} ここで、右辺第1項,第2項に対し \begin{eqnarray} \eta^{V}(s)&\equiv&\sum_{s^{\prime},a}P(s^{\prime}|s,a)\pi(a|s)r(s,a,s^{\prime}) \\ M^{V}(s,s^{\prime})&\equiv&\sum_{a}P(s^{\prime}|s,a)\pi(a|s) \end{eqnarray} を定義すると \begin{equation} V(s)=\eta^{V}(s)+\gamma \sum_{s^{\prime}}M^{V}(s,s^{\prime})V(s^{\prime}) \end{equation} を得る。いま、状態が$D$個に離散化されている場合を考えると、$V(s)$は$D$次元ベクトル$\vec{V}$の第$s$番目の成分と見ることができる。$M^{V}(s,s^{\prime})$は行列$M^{V}$の$(s,s^{\prime})$成分である。従って、次式を得る。 \begin{equation} \vec{V}=\vec{\eta}^{\;V}+\gamma M^{V}\vec{V} \end{equation} この式の右辺を眺めると、平行移動と線形変換の組み合わせ、すなわち、アフィン変換であることが分かる。このアフィン変換をベルマン作用素と呼ぶ。ベルマン作用素を$T^{V}$で表すと \begin{equation} \vec{V}=T^{V}\vec{V} \end{equation} となる。これは、変換後の自身が自身に等しくなることから不動点方程式であり、$0<\gamma<1$のとき唯一解$\vec{V}$を持つことが知られている。

     次に、行動価値関数$Q(s,a)$について考える。$Q(s,a)$のベルマン方程式は次式であった。 \begin{equation} Q(s,a)=\sum_{s^{\prime}} P(s^{\prime}|s,a)\left[r(s,a,s^{\prime})+\gamma \sum_{a^{\prime}}\pi(a^{\prime}|s^{\prime})Q(s^{\prime},a^{\prime})\right] \end{equation} ここで \begin{eqnarray} \eta^{Q}(s,a)&\equiv&\sum_{s^{\prime}}P(s^{\prime}|s,a)r(s,a,s^{\prime}) \\ M^{Q}(s,a,s^{\prime},a^{\prime})&\equiv&P(s^{\prime}|s,a)\pi(a^{\prime}|s^{\prime}) \end{eqnarray} を導入すると \begin{equation} Q(s,a)=\eta^{Q}(s,a)+\gamma\sum_{s^{\prime},a^{\prime}}M^{Q}(s,a,s^{\prime},a^{\prime})Q(s^{\prime},a^{\prime}) \end{equation} を得る。今度は状態だけでなく行動も$D^{\prime}$個に離散化されている場合を考える。このとき$Q(s,a)$は$D\times D^{\prime}$次元のベクトルになる。$(s,a)$を1つの添え字で表せば、状態価値関数のときと同じ議論を繰り返すことができ、次式を得る。 \begin{equation} \vec{Q}=\vec{\eta}^{\;Q}+\gamma M^{Q}\vec{Q} \end{equation} 右辺はアフィン変換であり、ベルマン作用素$T^{Q}$を使うと \begin{equation} \vec{Q}=T^{Q}\vec{Q} \end{equation} と書くことができる。これも不動点方程式であり、$0<\gamma<1$のとき唯一解$\vec{Q}$を持つことが知られている。

    ベルマン最適方程式


     $V(s)$と$Q(s,a)$のベルマン方程式を再度示す。 \begin{eqnarray} V(s)&=&\sum_{s^{\prime},a} P(s^{\prime}|s,a)\pi(a|s)\left[r(s,a,s^{\prime})+\gamma V(s^{\prime})\right]\\ Q(s,a)&=&\sum_{s^{\prime}} P(s^{\prime}|s,a)\left[r(s,a,s^{\prime})+\gamma \sum_{a^{\prime}}\pi(a^{\prime}|s^{\prime})Q(s^{\prime},a^{\prime})\right] \end{eqnarray} これらは、方策$\pi$を用いて、ある状態に対する行動を選択している。いま、方策$\pi$を介さずに、最適な行動を直接選択することを考える。上式の$\sum_{a}\pi(a|s)$を形式的に$\max_a$で置き換えると次式を得る。 \begin{eqnarray} V^{*}(s)&=&\max_{a}\sum_{s^{\prime}} P(s^{\prime}|s,a)\left[r(s,a,s^{\prime})+\gamma V^{*}(s^{\prime})\right]\\ Q^{*}(s,a)&=&\sum_{s^{\prime}} P(s^{\prime}|s,a)\left[r(s,a,s^{\prime})+\gamma \max_{a^{\prime}}Q^{*}(s^{\prime},a^{\prime})\right] \end{eqnarray} これらは、最適状態価値関数$V^{*}(s)$と最適行動価値関数$Q^{*}(s,a)$に対し、厳密に成り立つことが知られており、ベルマン最適方程式と呼ばれる。これらについては不動点方程式の形に変形することはできないが、やはり、唯一解が存在することが知られている。

    参考文献


    2018年2月13日火曜日

    強化学習 〜方策勾配定理の導出〜

    はじめに


     前回に引き続き、強化学習のテキスト「これからの強化学習」に出てくる方策勾配定理を導出する。自身のための覚書である。

    最大化すべき目的関数


     状態価値関数$V(s)$は次式で定義された。 \begin{equation} V(s)={\bf E}\left[G_t|S_t=s\right] \end{equation} いま、 方策を表す確率$\pi(a|s)$を、パラメータ$\theta$に依存する微分可能な関数でモデル化し、 時間ステップ$t=0$から始める状態価値関数を考える。 \begin{eqnarray} V(s_0)&=&{\bf E}\left[G_0|S_0=s_0\right] \\ &=&{\bf E}\left[\sum_{t=1}^{\infty}\gamma^{t-1}R_t|S_0=s_0\right] \\ &\equiv& J(\theta;s_0) \end{eqnarray} これを$\theta$に関して最大化する。

    方策勾配定理の導出


     先に見たように、状態価値関数$V(s)$と行動価値関数$Q(s,a)$の間には次式が成り立つ。 \begin{equation} V(s)=\sum_a \pi(a|s)Q(s,a) \end{equation} $\pi(a|s)$が$\theta$に依存するとき、$V(s)$と$Q(s,a)$も$\theta$に依存する。従って次式が成り立つ。 \begin{equation} \frac{\partial V(s)}{\partial \theta}=\sum_a \left[\frac{\partial \pi(a|s)}{\partial \theta}Q(s,a)+\pi(a|s)\frac{\partial Q(s,a)}{\partial \theta}\right] \label{eq0} \end{equation} ここで、先に示した$Q(s,a)$についてのベルマン方程式 \begin{equation} Q(s,a)=\sum_{s^{\prime}} P(s^{\prime}|s,a)\left[r(s,a,s^{\prime})+\gamma\;V(s^{\prime})\right] \end{equation} の両辺を$\theta$で微分する。上式の右辺第1項は$\theta$に依存しないことに注意すると \begin{equation} \frac{\partial Q(s,a)}{\partial \theta}=\gamma\;\sum_{s^{\prime}} P(s^{\prime}|s,a)\;\frac{\partial V(s^{\prime})}{\partial \theta} \end{equation} を得る。これを、式($\ref{eq0}$)に代入する。 \begin{eqnarray} \frac{\partial V(s)}{\partial \theta} &=& \sum_a \left\{\frac{\partial \pi(a|s)}{\partial \theta}Q(s,a) + \pi(a|s) \left[ \gamma\;\sum_{s^{\prime}} P(s^{\prime}|s,a)\;\frac{\partial V(s^{\prime})}{\partial \theta} \right] \right\}\\ &=& f(s) + \sum_a \pi(a|s) \gamma\;\sum_{s^{\prime}} P(s^{\prime}|s,a)\;\frac{\partial V(s^{\prime})}{\partial \theta} \end{eqnarray} ここで、$f(s)\equiv\sum_a \frac{\partial \pi(a|s)}{\partial \theta}Q(s,a)$とした。再帰的に代入を繰り返す。 \begin{eqnarray} \frac{\partial V(s)}{\partial \theta} &=& f(s) + \sum_a \pi(a|s) \gamma\;\sum_{s^{\prime}} P(s^{\prime}|s,a) \left[ f(s^{\prime}) + \sum_{a^{\prime}} \pi(a^{\prime}|s^{\prime}) \gamma\;\sum_{s^{\prime\prime}} P(s^{\prime\prime}|s^{\prime},a^{\prime}) \frac{\partial V(s^{\prime\prime})}{\partial \theta} \right]\\ &=& f(s) + \sum_a \pi(a|s) \gamma\;\sum_{s^{\prime}} P(s^{\prime}|s,a) f(s^{\prime}) \\ &&+ \sum_a \pi(a|s) \gamma\;\sum_{s^{\prime}} P(s^{\prime}|s,a) \sum_{a^{\prime}} \pi(a^{\prime}|s^{\prime}) \gamma\;\sum_{s^{\prime\prime}} P(s^{\prime\prime}|s^{\prime},a^{\prime}) \frac{\partial V(s^{\prime\prime})}{\partial \theta} \\ &=& f(s) + \sum_a \pi(a|s) \gamma\;\sum_{s^{\prime}} P(s^{\prime}|s,a) f(s^{\prime}) \\ &&+ \sum_a \pi(a|s) \gamma\;\sum_{s^{\prime}} P(s^{\prime}|s,a) \sum_{a^{\prime}} \pi(a^{\prime}|s^{\prime}) \gamma\;\sum_{s^{\prime\prime}} P(s^{\prime\prime}|s^{\prime},a^{\prime}) f(s^{\prime\prime})+\cdots \end{eqnarray} ここで、 \begin{equation} \sum_a \pi(a|s)P(s^{\prime}|s,a)=\sum_a P(s^{\prime}|s,a)\pi(a|s)=P(s^{\prime}|s) \end{equation} が成り立つから次式を得る。 \begin{equation} \frac{\partial V(s)}{\partial \theta} =f(s)+\gamma \sum_{s^{\prime}} P(s^{\prime}|s)f(s^{\prime}) +\gamma^2 \sum_{s^{\prime},s^{\prime\prime}} P(s^{\prime\prime}|s^{\prime}) P(s^{\prime}|s)f(s^{\prime\prime})+\cdots \end{equation} 右辺第2項の$P(s^{\prime}|s)$は1ステップで状態$s$から$s^{\prime}$へ遷移する確率、第3項の$\sum_{s^{\prime}}P(s^{\prime\prime}|s^{\prime}) P(s^{\prime}|s)$は2ステップで状態$s$から$s^{\prime\prime}$へ遷移する確率を表す。これを一般化し、$k$ステップで状態$s$から$x$へ遷移する確率を$P(s\rightarrow x,k)$と書くことにすると \begin{eqnarray} \frac{\partial V(s)}{\partial \theta} &=&f(s)+\gamma \sum_{x} P(s\rightarrow x,1)f(x) +\gamma^2 \sum_{x} P(s\rightarrow x,2)f(x)+\cdots \\ &=& \sum_{k=0}^{\infty}\gamma^{k}\sum_x P(s\rightarrow x,k)f(x) \end{eqnarray} を得る。ただし、$k=0$のとき状態は変化しないので次式が成り立つことを用いた。 \begin{equation} \sum_x P(s\rightarrow x,0)=1 \end{equation} ところで、$J(\theta;s_0)$は$V(s_0)$であったから \begin{eqnarray} \frac{\partial J(\theta;s_0)}{\partial \theta} &=& \frac{\partial V(s_0)}{\partial \theta} \\ &=& \sum_{k=0}^{\infty}\gamma^{k}\sum_x P(s_0\rightarrow x,k)f(x) \end{eqnarray} が成り立つ。$f(x)$を元の式に戻して \begin{eqnarray} \frac{\partial J(\theta;s_0)}{\partial \theta} &=& \sum_s \left[\sum_{k=0}^{\infty}\gamma^{k}P(s_0\rightarrow s,k)\right]\sum_a \frac{\partial \pi(a|s)}{\partial \theta}Q(s,a)\\ &=& \sum_s d(s)\sum_a \frac{\partial \pi(a|s)}{\partial \theta}Q(s,a) \end{eqnarray} を得る。ここで、$d(s)\equiv \sum_{k=0}^{\infty}\gamma^{k}P(s_0\rightarrow s,k)$と置いた。上式をさらに変形すると \begin{eqnarray} \frac{\partial J(\theta;s_0)}{\partial \theta} &=& \sum_{s,a} d(s)\pi(a|s)\frac{1}{\pi(a|s)} \frac{\partial \pi(a|s)}{\partial \theta}Q(s,a) \label{eq1} \end{eqnarray} を得る。ここで、右辺の$d(s)\pi(a|s)$は以下のように変形できる。 \begin{eqnarray} d(s)\pi(a|s) &=& \sum_{k=0}^{\infty}\gamma^{k}P(s_0\rightarrow s,k)\pi(a|s)\\ &=& \sum_{k=0}^{\infty}\gamma^{k}P(S_k=s|S_0=s_0)\pi(a|s)\\ &=& \sum_{k=0}^{\infty}\gamma^{k}P(S_k=s|S_0=s_0)\;P(A_k=a|S_k=s)\\ &=& \sum_{k=0}^{\infty}\gamma^{k}P(S_k=s, A_k=a|S_0=s_0) \end{eqnarray} 上式は、時間ステップ$t=0$において$s_0$であった状態が、最終的に状態$s$・行動$a$に遷移する全てのステップを足し合わせた確率を表している。割引率$\gamma$により、ステップ数が多いほど確率が低くなることが考慮されている。以上の考察から、式($\ref{eq1}$)は期待値の記号を用いて表すことができる。 \begin{eqnarray} \frac{\partial J(\theta;s_0)}{\partial \theta} &=& {\bf E}\left[\frac{1}{\pi(a|s)} \frac{\partial \pi(a|s)}{\partial \theta}Q(s,a)\right]\\ &=& {\bf E}\left[\frac{\partial \ln{\pi(a|s)}}{\partial \theta}Q(s,a)\right] \end{eqnarray} 上式を方策勾配定理と呼ぶ。

    参考文献


    2018年2月12日月曜日

    強化学習 〜ベルマン方程式の導出〜

    はじめに


     強化学習のテキスト「これからの強化学習」に出てくるベルマン方程式を導出する。自身のための覚書である。

    収益


     時間ステップ$t$における収益$G_t$として、次式で定義される割引報酬和を考える。 \begin{equation} G_t=\sum_{\tau=0}^{\infty}\gamma^{\tau}R_{t+1+\tau} \end{equation} ここで、$\gamma$は割引率と呼ばれる量であり、$0\leqq\gamma\leqq 1$を満たす。$R_t$は時間ステップ$t$における報酬を表す。上式を展開すると \begin{equation} G_t=R_{t+1}+\gamma R_{t+2}+\gamma^2 R_{t+3}+\cdots \end{equation} となる。つまり収益$G_t$とは、次の時間ステップ$t+1$から未来に向かって得られる報酬の和である。

    状態価値関数


     状態価値関数$V(s)$は、時間ステップ$t$における状態を表す確率変数$S_t$の観測値が$s$である条件の下で計算される収益$G_t$の期待値である。 \begin{equation} V(s)={\bf E}\left[G_t|S_t=s \right] \end{equation}

    行動価値関数


     行動価値関数$Q(s,a)$は、確率変数$S_t$の観測値が$s$、行動を表す確率変数$A_t$の観測値がaである条件の下で計算される収益$G_t$の期待値である。 \begin{equation} Q(s,a)={\bf E}\left[G_t|S_t=s,A_t=a \right] \end{equation}  ところで、$S_t=s$のとき収益$G_t$が実現する確率$P(G_t|S_t=s)$と$S_t=s, A_t=a$のとき収益$G_t$が実現する確率$P(G_t|S_t=s,A_t=a)$の間にはベイズの定理から次式が成り立つ。 \begin{equation} P(G_t|S_t=s)=\sum_a P(A_t=a|S_t=s)P(G_t|S_t=s,A_t=a) \end{equation} ここで、$P(A_t=a|S_t=s)$は、状態$s$のとき行動$a$が選択される確率であり、特に$\pi(a|s)$と表記される。これはエージェントの方策を決める確率である。上の関係式から、状態価値関数と行動価値関数の間には次式が成り立つ。 \begin{equation} V(s)=\sum_{a}\pi(a|s)Q(s,a) \end{equation}

    ベルマン方程式


     状態価値関数$V(s)$に式(2)を代入して \begin{eqnarray} V(s)&=&{\bf E}\left[G_t|S_t=s \right] \\ &=&{\bf E}\left[R_{t+1}+\gamma R_{t+2}+\gamma^2 R_{t+3}+\cdots|S_t=s \right] \\ &=&{\bf E}\left[R_{t+1}|S_t=s \right] +\gamma{\bf E}\left[R_{t+2}+\gamma R_{t+3}+\cdots|S_t=s \right] \end{eqnarray} を得る。右辺第1項は \begin{eqnarray} {\bf E}\left[R_{t+1}|S_t=s \right] &=&\sum_{s^{\prime},a}P(S_{t+1}=s^{\prime},A_t=a|S_t=s)\;R_{t+1} \\ &=&\sum_{s^{\prime},a}P(S_{t+1}=s^{\prime}|S_t=s,A_t=a)P(A_t=a|S_t=s)\;R_{t+1} \\ &=&\sum_{s^{\prime},a}P(s^{\prime}|s,a)\pi(a|s)\;r(s,a,s^{\prime}) \end{eqnarray} となる。ここで、報酬関数$r(S_t,A_t,S_{t+1})$を導入した。これは現在の状態と行動、および次の状態から報酬が決定されることを表している。また、確率$P$は、確率変数ではなくその実現値に依存するものである。従って、確率変数を省略することができる。
    次に、式(9)の第2項を考える。 \begin{equation} {\bf E}\left[R_{t+2}+\gamma R_{t+3}+\cdots|S_t=s \right] = {\bf E}\left[R_{t+2}|S_t=s \right]+ \gamma{\bf E}\left[R_{t+3}|S_t=s \right]+\cdots \end{equation} これの右辺第1項は \begin{eqnarray} {\bf E}\left[R_{t+2}|S_t=s \right] &=& \sum_{s^{\prime\prime},s^{\prime},a^{\prime},a} P(S_{t+2}=s^{\prime\prime},A_{t+1}=a^{\prime},S_{t+1}=s^{\prime},A_t=a|S_t=s)\;R_{t+2} \\ &=& \sum_{s^{\prime\prime},s^{\prime},a^{\prime},a} P(S_{t+2}=s^{\prime\prime},A_{t+1}=a^{\prime}|S_{t+1}=s^{\prime})P(S_{t+1}=s^{\prime},A_t=a|S_t=s)\;R_{t+2} \\ &=&\sum_{s^{\prime},a} P(S_{t+1}=s^{\prime},A_t=a|S_t=s)\sum_{s^{\prime\prime},a^{\prime}}P(S_{t+2}=s^{\prime\prime},A_{t+1}=a^{\prime}|S_{t+1}=s^{\prime})\;R_{t+2} \\ &=&\sum_{s^{\prime},a} P(s^{\prime},a|s)\;{\bf E}\left[R_{t+2}|S_{t+1}=s^{\prime} \right] \\ &=&\sum_{s^{\prime},a} P(s^{\prime}|s,a)\pi(a|s)\;{\bf E}\left[R_{t+2}|S_{t+1}=s^{\prime} \right] \end{eqnarray} となる。式(16)からの(17) への変形では、確率変数を省略し、式(10)を用いた。式(13)の他の項についても同様に計算することにより次式を得る。 \begin{eqnarray} {\bf E}\left[R_{t+2}+\gamma R_{t+3}+\cdots|S_t=s \right]&=&\sum_{s^{\prime},a} P(s^{\prime}|s,a)\pi(a|s)\;{\bf E}\left[R_{t+2}+\gamma R_{t+3}+\cdots|S_{t+1}=s^{\prime} \right] \\ &=& \sum_{s^{\prime},a} P(s^{\prime}|s,a)\pi(a|s)\;V(s^{\prime}) \end{eqnarray} 式(19)から(20)への変形では$V(s)$の定義を用いた。最後に式(12)と(20)から次式を得る。 \begin{equation} V(s)=\sum_{s^{\prime},a} P(s^{\prime}|s,a)\pi(a|s)\left[r(s,a,s^{\prime})+\gamma V(s^{\prime})\right] \end{equation} これを状態価値関数に関するベルマンの方程式と呼ぶ。

     次に行動価値関数についてのベルマン方程式を求める。 \begin{eqnarray} Q(s,a)&=&{\bf E}\left[G_t|S_t=s,A_t=a \right] \\ &=&{\bf E}\left[R_{t+1}+\gamma R_{t+2}+\gamma^2 R_{t+3}+\cdots|S_t=s,A_t=a \right] \\ &=&{\bf E}\left[R_{t+1}|S_t=s,A_t=a \right] +\gamma{\bf E}\left[R_{t+2}+\gamma R_{t+3}+\cdots|S_t=s,A_t=a \right] \end{eqnarray} 右辺第1項は \begin{eqnarray} {\bf E}\left[R_{t+1}|S_t=s,A_t=a \right] &=&\sum_{s^{\prime}}P(S_{t+1}=s^{\prime}|S_t=s,A_t=a)\;R_{t+1} \\ &=&\sum_{s^{\prime}}P(s^{\prime}|s,a)\;r(s,a,s^{\prime}) \end{eqnarray} 右辺第2項は \begin{equation} {\bf E}\left[R_{t+2}+\gamma R_{t+3}+\cdots|S_t=s,A_t=a \right] = {\bf E}\left[R_{t+2}|S_t=s,A_t=a \right]+ \gamma{\bf E}\left[R_{t+3}|S_t=s,A_t=a \right]+\cdots \end{equation} これの右辺第1項は \begin{eqnarray} &&{\bf E}\left[R_{t+2}|S_t=s,A_t=a \right] \\ &=& \sum_{s^{\prime\prime},s^{\prime},a^{\prime}} P(S_{t+2}=s^{\prime\prime},A_{t+1}=a^{\prime},S_{t+1}=s^{\prime}|S_t=s,A_t=a)\;R_{t+2} \\ &=& \sum_{s^{\prime\prime},s^{\prime},a^{\prime}} P(S_{t+2}=s^{\prime\prime},A_{t+1}=a^{\prime}|S_{t+1}=s^{\prime})P(S_{t+1}=s^{\prime}|S_t=s,A_t=a)\;R_{t+2} \\ &=&\sum_{s^{\prime}} \sum_{s^{\prime\prime},s^{\prime},a^{\prime}} P(S_{t+2}=s^{\prime\prime}|S_{t+1}=s^{\prime},A_{t+1}=a^{\prime}) P(A_{t+1}=a^{\prime}|S_{t+1}=s^{\prime}) \times\\ && \hspace{200pt}P(S_{t+1}=s^{\prime}|S_t=s,A_t=a)\;R_{t+2} \\ &=&\sum_{s^{\prime}} P(S_{t+1}=s^{\prime}|S_t=s,A_t=a) \sum_{a^{\prime}} P(A_{t+1}=a^{\prime}|S_{t+1}=s^{\prime})\times \\ && \hspace{200pt}\sum_{s^{\prime\prime}} P(S_{t+2}=s^{\prime\prime}|S_{t+1}=s^{\prime},A_{t+1}=a^{\prime})\;R_{t+2} \\ &=&\sum_{s^{\prime}} P(S_{t+1}=s^{\prime}|S_t=s,A_t=a) \sum_{a^{\prime}} P(A_{t+1}=a^{\prime}|S_{t+1}=s^{\prime}) \;{\bf E}\left[R_{t+2}|S_{t+1}=s^{\prime},A_{t+1}=a^{\prime} \right] \\ &=&\sum_{s^{\prime}} P(s^{\prime}|s,a) \sum_{a^{\prime}} \pi(a^{\prime}|s^{\prime}) \;{\bf E}\left[R_{t+2}|S_{t+1}=s^{\prime},A_{t+1}=a^{\prime} \right] \end{eqnarray} 他の項についても同様な計算を行えば次式を得る。 \begin{eqnarray} &&{\bf E}\left[R_{t+2}+\gamma R_{t+3}+\cdots|S_t=s,A_t=a \right]\\ &=& \sum_{s^{\prime}} P(s^{\prime}|s,a)\sum_{a^{\prime}} \pi(a^{\prime}|s^{\prime}) \;{\bf E}\left[R_{t+2}+\gamma R_{t+3}+\cdots|S_{t+1}=s^{\prime},A_{t+1}=a^{\prime} \right] \\ &=& \sum_{s^{\prime}} P(s^{\prime}|s,a)\sum_{a^{\prime}} \pi(a^{\prime}|s^{\prime}) \;Q(s^{\prime},a^{\prime}) \\ \end{eqnarray} 式(26)と(37)から最終的に次式を得る。 \begin{equation} Q(s,a)=\sum_{s^{\prime}} P(s^{\prime}|s,a)\left[r(s,a,s^{\prime})+\gamma \sum_{a^{\prime}}\pi(a^{\prime}|s^{\prime})Q(s^{\prime},a^{\prime})\right] \end{equation} 式(6)を使えば \begin{equation} Q(s,a)=\sum_{s^{\prime}} P(s^{\prime}|s,a)\left[r(s,a,s^{\prime})+\gamma\;V(s^{\prime})\right] \end{equation} とも書ける。これらが、行動価値関数についてのベルマン方程式である。