論文読み「Your Classifier Is Secretly An Energy Based Model And You Should Treat It Like One」

はじめに

　「Your Classifier Is Secretly An Energy Based Model And You Should Treat It Like One」をまとめる。

概要

同時分布$p({\bf x},y)$をEnergy Based Modelで表すことにより、従来の識別器（discriminative classifier：$p(y|{\bf x})$）から$p({\bf x},y)$を導出することができる。

$p({\bf x})$も容易に求めることができる。

従来の手順で求められる識別器が出発点となる。ラベル付けされていないデータを用いて訓練される。

本手法を用いて同時分布$p({\bf x},y)$を訓練することにより、較正度合い、敵対的サンプルに対する頑健さ、分布外検知を改善することができる。それと同時に、最近のGANに匹敵する高画質画像を生成することができる。

Energy Based Modelの訓練を大規模に行うテクニックを改善し、通常の識別器の訓練とほとんど同じオーバヘッドで行う手法を与える。

Energy Based Model

　パラメータ$\theta$を持つ同時確率分布$p_{\theta}({\bf x},y)$をボルツマン分布を用いて表す。 \begin{equation} p_{\theta}({\bf x},y)=\frac{\exp{ \left( -E_{\theta} ({\bf x},y) \right) } }{Z(\theta)} \label{eq1} \end{equation} $E_{\theta}({\bf x},y)$はエネルギーであり、その定義については後述する。$Z(\theta)$は次式で定義される分配関数である。 \begin{equation} Z(\theta)=\int dx\int dy \exp{ \left( -E_{\theta} ({\bf x},y) \right) } \label{eq2} \end{equation} 式(\ref{eq1})のように確率分布をボルツマン分布を用いて表すモデルを、Energy Based Modelと呼ぶ。式(\ref{eq1})を$y$について周辺化して \begin{equation} p_{\theta}({\bf x}) = \int dy\;p_{\theta}({\bf x},y) \label{eq3} \end{equation} を得る。いま、$p_{\theta}({\bf x})$に対してもボルツマン分布を仮定する。 \begin{equation} p_{\theta}({\bf x})=\frac{\exp{\left(-E_{\theta}({\bf x})\right)}}{Z(\theta)} \label{eq5} \end{equation} これを変形すると \begin{eqnarray} \exp{\left(-E_{\theta}({\bf x})\right)} &=& Z(\theta)p_{\theta}({\bf x})\nonumber\\ &=& Z(\theta) \int dy\;p_{\theta}({\bf x},y)\nonumber\\ &=& Z(\theta) \int dy \frac{\exp{ \left( -E_{\theta} ({\bf x},y) \right) }}{Z(\theta)}\nonumber\\ &=& \int dy \exp{ \left( -E_{\theta} ({\bf x},y) \right) } \label{eq6} \end{eqnarray} 両辺の対数を取り \begin{equation} E_{\theta}({\bf x})=-\ln{\left( \int dy \exp{ \left( -E_{\theta} ({\bf x},y) \right) } \right) } \label{eq7} \end{equation} を得る。ここまでの手順をまとめると以下のようになる。

1. エネルギー$E_{\theta}({\bf x},y)$を与える。
2. 式(\ref{eq1})から$p_{\theta}({\bf x},y)$を計算する。
3. 式(\ref{eq7})から$E_{\theta}({\bf x})$を計算する。
4. 式(\ref{eq5})から$p_{\theta}({\bf x})$を計算する。

つまり、エネルギー$E_{\theta}({\bf x},y)$が決まれば、さまざまな量を求めることができる。

エネルギーの定義

　$E_{\theta}({\bf x},y)$について考えるため、従来の識別器の構成を示す（下の左図）。

Softmaxに入る前の出力を$f_{\theta}({\bf x})$と書くと、Softmaxによる出力が従来の識別器$p_{\theta}(y|{\bf x})$に相当する。 \begin{equation} p_{\theta}(y|{\bf x}) = \frac { \exp{\left(f_{\theta}({\bf x})[y] \right)} } { \sum_{y^{\prime}}\exp{\left(f_{\theta}({\bf x})[y^{\prime}] \right)} } \end{equation} ここで、$f_{\theta}({\bf x})[y]$は$f_{\theta}({\bf x})$の$y$番目の成分である。本文献では次式によりエネルギーを定義する。 \begin{equation} E_{\theta}({\bf x},y)=-f_{\theta}({\bf x})[y] \end{equation} 上の右側の図は、出力$f_{\theta}({\bf x})$から$p_{\theta}({\bf x},y)$や$p_{\theta}({\bf x})$が導出される様子を示している。本文献の手法をJoint Energy Model（JEM）と呼ぶ。

最適化

　学習時に最大化する目的関数は次式である。 \begin{equation} \ln{p_{\theta}({\bf x},y)}=\ln{p_{\theta}({\bf x})}+\ln{p_{\theta}(y|{\bf x})} \label{eq9} \end{equation} 第2項の最大化は従来とおりクロスエントロピーを用いれば良い。第1項の最大化については以下のようにする。最初に、式(\ref{eq5})の両辺の対数をとる。 \begin{equation} \ln{p_{\theta}({\bf x})}=-E_{\theta}({\bf x})-\ln{Z(\theta)} \end{equation} $\theta$で微分すると \begin{equation} \frac{\partial\ln{p_{\theta}({\bf x})}}{\partial\theta}=-\frac{\partial E_{\theta}({\bf x})}{\partial\theta} -\frac{1}{Z(\theta)}\frac{\partial Z(\theta)}{\partial\theta} \end{equation} ここで \begin{equation} Z(\theta)=\int dx \exp{\left(-E_{\theta}({\bf x})\right)} \end{equation} であるから、これを$\theta$で微分すると \begin{equation} \frac{\partial Z(\theta)}{\partial\theta}=\int dx \exp{\left(-E_{\theta}({\bf x})\right)}\left(-\frac{\partial E_{\theta}({\bf x})}{\partial\theta}\right) \end{equation} 故に \begin{eqnarray} \frac{1}{Z(\theta)}\frac{\partial Z(\theta)}{\partial\theta} &=& -\int dx \frac{\exp{\left(-E_{\theta}({\bf x})\right)}}{Z(\theta)} \left(\frac{\partial E_{\theta}({\bf x})}{\partial\theta}\right)\nonumber\\ &=& -\int dx p_{\theta}({\bf x}) \left(\frac{\partial E_{\theta}({\bf x})}{\partial\theta}\right)\nonumber\\ &=&-\mathbb{E}_{ p_{\theta}({\bf x}) } \left[ \frac{\partial E_{\theta}({\bf x})}{\partial\theta} \right] \end{eqnarray} となる。最後の等式は$p_{\theta}({\bf x})$による期待値であることを表す。以上から \begin{equation} \frac{\partial\ln{p_{\theta}({\bf x})}}{\partial\theta}=-\frac{\partial E_{\theta}({\bf x})}{\partial\theta} +\mathbb{E}_{ p_{\theta}({\bf x}) } \left[ \frac{\partial E_{\theta}({\bf x})}{\partial\theta} \right] \end{equation} を得る。$E_{\theta}({\bf x})$は式(\ref{eq7})により計算することができる。上式の第2項は分布$p_{\theta}({\bf x})$からサンプリングする必要がある。 本文献で使用するサンプリング手法は以下の通りである。 \begin{eqnarray} {\bf x}_{0}&\sim&p_0({\bf x}) \nonumber\\ {\bf x}_{i+1}&=&{\bf x}_i-\frac{\alpha}{2}\frac{\partial E_{\theta}({\bf x}_i)}{\partial {\bf x}_i}+\epsilon\nonumber\\ \epsilon&\sim&\mathcal{N}(0,\alpha) \end{eqnarray} ここで、$p_0({\bf x})$は一様分布である。本手法は、Stochastic Gradient Langevin Dynamics（SGLD）をベースにしたものである（詳細は略）。式(\ref{eq9})の第1項は教師なし学習、第2項は教師あり学習を行う。

実験

識別モデル、生成モデルのどちらに対しても、現在SOTAとされている手法に匹敵する性能を実現した。

識別モデルにおいて、不確かさ評価（uncertainty quantification）、分布外検知、敵対的サンプル（adversarial example）に対する頑健さ、が改善された。　

本文献で用いるネットワークは、Wide Residaul Networkである。このネットワークから$f_{\theta}({\bf x})$を取り出す。

データセットCIFAR10を用いてハイパーパラメータを調節した。本文献で検証するモデルは全て同じハイパーパラメータを用いた。

SGLDサンプラーで使うパラメータは、本文献で扱うデータセットやモデルアーキテクチャの全てで共有できた。

訓練時間は、1つのGPUで36時間ほどであった。

Hybrid Modelint

　従来手法と比較したのが下図である。使用したデータセットは、CIFAR10である。

Class欄に「Hybrid」とあるのが、識別モデルと生成モデルを同時に構築する手法（本手法JEMも含まれる）、「Disc.」欄、「Gen.」欄にあるのが　それぞれ識別モデルだけ、生成モデルだけの手法である。IS（Inception Score）は大きいほど、FID（Frechet Inception Distance）は小さいほど良い指標である。FID以外はJEMが優位であることが分かる。 上の表の4行目にある「JEM$p({\bf x}|y)$ factored」は式(\ref{eq9})の代わりに次式を最大化した結果である。 \begin{equation} \ln{p_{\theta}({\bf x},y)}=\ln{p_{\theta}(y)}+\ln{p_{\theta}({\bf x}|y)} \end{equation} このとき、精度が大きく落ちることが分かる。　

Calibration

　較正済み識別器（calibrated classifier）とは、$\max_y{p(y|{\bf x})}$（確信度）が正解率（accuracy）と比例するような識別器のことである。高い確信度を持って予測されたラベルは高い確率で正解となり、低い確信度で予測されたラベルは、誤認識である確率が高くなる。較正されているか否かを可視化する手順は以下の通り。

1. データセット内の各サンプル${\bf x}_i$ごとに$\max_y{p(y|{\bf x}_i)}$を算出し、これらを等区間のバケットに振り分ける（ヒストグラムを作る）。例えば0から1までの区間を20等分した場合、20個のバケットに振り分ける。
2. 各バケットに属するサンプルから平均精度を算出し、このバケットの高さとする。

上の手順を経ると以下のような図を得る。

上の右図のように確信度と精度が比例するとき、識別器は良く較正されているとみなすことできる。左図は従来手法で訓練した識別器による結果、右図はJEMで訓練した識別器による結果である。全サンプルを用いて算出した精度は、前者が74.2%、後者が72.2%であり、較正された識別器の方が悪い。従って、較正の良し悪しと精度の良さは一致するとは限らないことに注意しなければならない。実環境で有用とされる識別器は、較正されたものである。
　較正度合いを定量化する指標が、次式で定義されるExpected Calibration Error（ECE）である。 \begin{equation} {\rm ECE}=\sum_{m=1}^{M}\frac{|B_m|}{n}|{\rm acc}(B_m)-{\rm conf}(B_m)| \end{equation} ここで、$M$はバケットの数、$n$はデータセット内のサンプル総数、$B_m$は１つののバケット中に格納されるサンプルの集合、$|B_m|$はサンプル数、${\rm acc}(B_m)$は集合$B_m$を用いて計算された平均精度、${\rm conf}(B_m)$は集合$B_m$を用いて計算された平均確信度である。${\rm ECE}$が小さくなるほど較正度合いは良くなり、完全に較正されている場合は0になる。上図に見るように、従来手法による識別器より、JEMで訓練した識別器の方が${\rm ECE}$は小さい。

Out-Of-Distribution Detection

　分布外検知とは例えば以下のような識別器を作ることである（引用元）。

猫の画像が与えられる。この画像で訓練する（教師なし学習）。出来上がった識別器に猫の画像を与えると正常（in-distribution）を、フクロウを撮影した画像を与えると異常（out-of-distribution）を返す。

犬・鳥の2種類の画像で訓練し識別器を作る（教師あり学習）。出来上がった識別器に犬・鳥の画像を与えると正常を、馬の画像を与えると異常を返す。

　データセットCIFAR10で訓練した識別器に別の3つのデータセットSVHN、CIFAR100、CelebAを与えときのサンプルの対数分布（$\ln{p(x)}$）を示したのが下図である。

緑の分布がCIFAR10の分布、赤の分布が他のデータセットの分布を表す。「Glow $\log{p(x)}$」の列はGlowと呼ばれる手法による結果、「JEM $\log{p(x)}$」の列にある結果がJEMによるものである。赤い分布と緑の分布が良く分離できているのは後者の方である。

　分布外検知に成功している多くの手法が採用している指標が次式である。 \begin{equation} s_{\theta}({\bf x})=\max_y{p_{\theta}(y|{\bf x})} \end{equation}

表の真ん中の行は上の指標で各種手法を比較したものである。CIFAR10で訓練した場合、データセットSVHNとInterpのサンプルを分布外であると一番良く認識できているのはWide-ResNetである。一方、CIFAR100とCelebAをサンプルとして与えたときは、JEMの認識が一番良い。

Robustness

　JEMは敵対的サンプルに対しても頑健であると主張されている（詳細は略）。

まとめ

同時分布$p({\bf x},y)$をEnergy Based Modelで表すことにより、従来の識別器（discriminative classifier：$p(y|{\bf x})$）から$p({\bf x},y)$を導出することができる。

$p({\bf x})$も容易に求めることができる。

従来の手順で求められる識別器が出発点となる。ラベル付けされていないデータを用いて訓練される。

本手法を用いて同時分布$p({\bf x},y)$を訓練することにより、較正度合い、敵対的サンプルに対する頑健さ、分布外検知を改善することができる。それと同時に、最近のGANに匹敵する高画質画像を生成することができる。

Energy Based Modelの訓練を大規模に行うテクニックを改善し、通常の識別器の訓練とほとんど同じオーバヘッドで行う手法を与える。

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

はじめに

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

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

　以下のネットワークを考える。

各量の定義は以下の通り。 \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

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

ベルマン最適方程式と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では状態に対する行動をニューラルネットワークで決めることになる。

参考文献

• これからの強化学習
• Pythonによる深層強化学習入門 ChainerとOpenAI Gymではじめる強化学習

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} を利用すれば良い。ここまでの議論を反映したアルゴリズムが論文に掲載されている。

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の値の推移も確認した方が良い。

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の画像（２値画像）を与えることにする。$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のサンプルコードをベースにすればすぐに実装できそうである。

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数を増やすだけではなくコードの見直し（多層化、初期化関数の変更）も必要であろう。次回はこの辺りのことについてまとめたい。

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徹底解説

強化学習 〜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情報行列の導出

　２つの確率分布$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)$を近似するのだろう）。

参考文献

• 西尾さんのスライド。とても勉強になります。
• これからの強化学習
• 自然勾配学習法

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

はじめに

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

不動点方程式

　状態価値関数$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)$に対し、厳密に成り立つことが知られており、ベルマン最適方程式と呼ばれる。これらについては不動点方程式の形に変形することはできないが、やはり、唯一解が存在することが知られている。

参考文献

• 速習強化学習
• これからの強化学習

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

はじめに

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

最大化すべき目的関数

　状態価値関数$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)f(x)=f(s) \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} 上式を方策勾配定理と呼ぶ。

参考文献

• https://www.slideshare.net/nishio/3-71708970
• これからの強化学習

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

はじめに

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

収益

　時間ステップ$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} とも書ける。これらが、行動価値関数についてのベルマン方程式である。

word2vec その2

はじめに

　このページの続きです。

やりたいこと

　単語列 \begin{equation} W = \{w_1,\cdots, w_T\} \end{equation} が与えられているとき、この単語列に含まれる任意の単語$w_t$を考える。$w_t$を含む文からその周辺の単語の集合$\Xi_{w_t}$を作ることができる。 \begin{equation} \Xi_{w_t} = \{\xi_1^t, \cdots, \xi_C^t\} \end{equation} このとき、以下の対数尤度を考える。 \begin{eqnarray} L&=&\sum_{t=1}^T \sum_{c=1}^{C} \biggl\{ \log{\sigma(s(\xi_c^t,w_t))}+\sum_{w \in {\rm Ng}}\log{\sigma(-s(w,w_t))} \biggr\} \label{obj}\\ \sigma(x) &=& \frac{1}{1+\exp{(-x)}}\\ s(a,b) &=& \vec{u}^T_a \cdot \vec{v}_b \end{eqnarray} ここで、${\rm Ng}$は分布 \begin{equation} p(w) = \frac{U(w)^{0.75}}{\sum_{t=1}^{T}U(w_t)^{0.75}} \end{equation} に従ってサンプリングした$k$個の単語の集合である。$U(w)$は単語の出現頻度である。本ページでは、上記の対数尤度の最大化問題を、Neural Networkを用いて解くことができることを示す。

Neural Networkによる表現

　単語列$W$内の$t$番目の単語$w_t$を以下のone-hot vector $\vec{x}_t$で表現する。 \begin{equation} \vec{x}_t = \left( \begin{array}{c} 0 \\ \vdots \\ 1 \\ 0 \\ \vdots \\ 0 \end{array} \right) \end{equation} $t$番目の要素だけが1の$T$次元ベクトルである。次に、行列$W_I$と$W_O$を次式で定義する。 \begin{eqnarray} W_I&=&\left[\vec{v}_1, \cdots, \vec{v}_T\right] \\ W_O&=&\left( \begin{array}{c} \vec{u}_1^T \\ \vdots \\ \vec{u}_T^T \end{array} \right) \end{eqnarray} ここで、$\vec{v}_t$と$\vec{u}_t$は$M$次元ベクトルとする。従って、$W_I$は$M\times T$行列、$W_O$は$T\times M$行列である。これらを用いると \begin{eqnarray} (W_O\; W_I\; \vec{x}_t)_i &=& (W_O)_{ik}\;(W_I\;\vec{x}_t)_k \\ &=&(W_O)_{ik}\;(W_I)_{km} (\vec{x}_t)_{m}\\ &=& u_{ik}\;v_{mk}\;x_{tm} \\ &=& u_{ik}\;v_{mk}\;\delta_{tm} \\ &=& u_{ik}\;v_{tk} \\ &=& \vec{u}_{i}^T \cdot \vec{v}_{t} \\ &=& s(i, t) \end{eqnarray} となる。これを図にすると以下のようなる。

以上から以下のことが分る。

1. 式(\ref{obj})の第1項に含まれる$s(\xi_c^t, w_t)$を計算するには、$w_t$に相当するone-hot vector $\vec{x}_t$をネットワークに入力し、その出力ベクトル（$T$次元ペクトル）の成分のうち、$\xi_c^t$に相当するものを取り出せば良い。
2. 式(\ref{obj})の第2項に含まれる$s(w, w_t)$の計算も同様である。

これらの計算のあと活性化関数としてsigmoid関数を作用させ、従来の手法（確率的最急降下法など）を使って$L$を最適化すれば良い。この計算で求まる行列$W_I$を使って単語$w_t$を表現するベクトル$\vec{v}_t$が次式で決まる。 \begin{equation} \vec{v}_t = W_I\;\vec{x}_t \end{equation} $T \gg M$と取るので低次元で効率良く単語を表現できるベクトルを得ることができる。

Chainerによる実装

　Chainerのソースには、word2vecのサンプルプログラムが含まれている。このサンプルプログラムには、Skip Gram以外にContinuous BoWも実装されており、さらに、Negative Sampling以外の損失関数も選択できるようになっている。以下に示すのは、サンプルプログラムを参考にして、Skip GramかつNegative Samplingの場合のクラスを実装したものである。

1. 5行目で定義されるL.EmbedIDは、$W_I$に相当する。
2. 6行目で定義されるL.NegativeSamplingは、式(\ref{obj})を計算する。この関数の中に$W_O$に相当するものがある。
3. n_vocabは$T$、n_unitsは$M$に相当する。
4. countsは$U(w)$に相当する。
5. sample_sizeは$k$に相当する。
6. 13行目の引数xはいま対象にする単語$w_t$、contextは$\Xi_{w_t}$に相当する。ただし、batch処理が考慮されている。

ここまでは式と対応付けることができたが、関数__call__の中で行っている処理で分らなくなる。最初にself.embed(context)を実行している。$\Xi_{w_t}$の中の1つの単語を$\vec{x}_t$とおくと、この計算は$W_I\;\vec{x}_t$に相当してしまう。定式化では$\Xi_{w_t}$に含まれる単語には$W_O$が作用するはずである。さてどこでおかしくなったのか？

word2vec その1

はじめに

　word2vecについてまとめます（その2を書きました）。

やりたいこと

　ある程度まとまった量の文章$A$を考える。これに出現する単語を重複を許さずに取り出すと、長さ$T$の単語列が出来上がる。 \begin{equation} W = \{w_1, \cdots, w_T\} \end{equation} いま、この単語列に含まれる単語$w_t$を含む文を$A$から取り出し、$w_t$の前後$b$個の範囲にある単語の集合$\Xi_{w_t}$を考える。 \begin{equation} \Xi_{w_t} = \{\xi_1^t, \cdots, \xi_C^t\} \end{equation} これは文脈を反映した集合である。たとえば、

Linux was originally developed as a free operating system for personal computers based on the Intel x86 architecture.

という文からは、$w_t={\rm developed}$、$b=1$のとき \begin{equation} \Xi_{{\rm developed}} = \{{\rm originally}, {\rm as}\} \end{equation} を、$w_t={\rm free}$のとき \begin{equation} \Xi_{{\rm free}} = \{{\rm a}, {\rm operating}\} \end{equation} を得る。$w_t$を指定したとき$\Xi_{w_t}$が実現する確率 \begin{equation} p(\Xi_{w_t}|w_t) \end{equation} を最大化するような単語の表現（ベクトル）を求めることがword2vec（Skip-Gram model）の目的である。

対数尤度の最大化

\begin{eqnarray} p(\Xi_t|w_t) &=& p(\xi_1^t, \cdots, \xi_C^t|w_t) \\ &=& \prod_{c=1}^{C}p(\xi_c^t|w_t) \end{eqnarray} 対数をとり、全ての単語について足し合わせたもの \begin{equation} L=\sum_{t=1}^T \sum_{c=1}^{C}\log{p(\xi_c^t|w_t)} \label{obj} \end{equation} を最大化する。

Softmax関数による計算

　$p(\xi_c^t|w_t)$としてSoftmax関数を仮定する。 \begin{equation} p(\xi_c^t|w_t)=\frac{\exp{\vec{u}^T_{\xi_c^{\;t}} \cdot \vec{v}_{w_t}}}{\sum_{w \in W} \exp{\vec{u}^T_w \cdot \vec{v}_{w_t}}} \end{equation} ここで、$\vec{v}$は入力単語を表現するベクトル、$\vec{u}$は出力単語を表現するベクトルである。$\sum_{w \in W}$は全単語について足し合わすことを意味する。モデルに入力する単語とモデルが出力する単語は区別されることに注意する。式(\ref{obj})に代入すると \begin{equation} L = \sum_{t=1}^T \sum_{c=1}^{C} \left\{ s(\xi_c^t, w_t) - \log{\sum_{w \in W} \exp{s(w, w_t)}} \right\} \end{equation} を得る。ただし、$s(a,b) = \vec{u}^T_a \cdot \vec{v}_b$とした。$W$の要素数が少なければ総和（$\sum_{w \in W}$）を取ることは可能であるが、現実の問題では$10^5\sim10^7$のオーダーであるため計算は困難である。

1次元ロジスティック回帰による計算（2分類問題への置き換え）

　単語$a$と$b$が与えられたとき、これらが同じ文脈上に存在する確率を$p(D=1|a, b)$、同じ文脈上に存在しない確率を$p(D=0|a, b)$とする。 このとき、次式のように書くことができる。 \begin{equation} p(D|a,b)=p(D=1|a, b)^{D}\; p(D=0|a, b)^{1-D} \end{equation} ここで以下を考える。

1. $w_t$が与えられた時の集合$\Xi_{w_t}$の各要素は、$w_t$と同じ文脈上に存在する。
2. $w_t$が与えられた時の集合$\Xi_{w_t}$に属さない要素は、$w_t$と同じ文脈上に存在しない。

これらを実現するには、次の対数尤度を最大化すればよい。 \begin{equation} L = \sum_{t=1}^T \sum_{c=1}^{C} \biggl\{ D_t \log{p(D_t=1|\xi^t_c, w_t)}+\sum_{w \not\in \Xi_{w_t}}(1-D_t)\log{p(D_t=0|w, w_t)} \biggr\} \label{f1} \end{equation} さらに \begin{eqnarray} p(D_t=1|a, b) &=& \sigma(s(a,b)) \\ p(D_t=0|a, b) &=& 1- \sigma(s(a,b)) = \sigma(-s(a,b)) \end{eqnarray} を仮定すれば（$\sigma(x)$はシグモイド関数）、1次元ロジスティック回帰の問題に帰着する。 式(\ref{f1})を最大にすることにより

1. 同じ文脈上にある2語の確率を高くし（第1項からの寄与）
2. 同じ文脈上にない2語の確率を高くする（第2項からの寄与）

ことが実現される。上記とほとんど同じことであるが次式を最大化する手法が存在する。 \begin{equation} L = \sum_{t=1}^T \sum_{c=1}^{C} \biggl\{ \log{p(D_t=1|\xi^t_c, w_t)}+\sum_{w \in {\rm Ng}}\log{p(D_t=0|w, w_t)} \biggr\} \end{equation} ${\rm Ng}$は、分布 \begin{equation} p(w) = \frac{U(w)^{0.75}}{\sum_{t=1}^{T}U(w_t)^{0.75}} \end{equation} に従ってサンプリングした$k$個の単語の集合である。$U(w)$は単語の頻度である。人為的に外れ(negative)データを導入することから、本手法をNegative Samplingと呼ぶ（その2を書きました）。

OpenGLで影を付ける方法

はじめに

　物体と平面を考える。OpenGLにて物体の影を平面上に描画する最も簡単な方法は影行列なるものを使う方法である。ここでは、この影行列を求め、実装してみる。

影行列の導出

　同次座標系で考える。すなわち、ベクトルは全て4次元ベクトル、行列は$4\times4$行列である。物体上の任意の点を $\vec{p}$、光源の位置を $\vec{l}$、平面の法線ベクトルを $\vec{n}$、点 $\vec{p}$ を平面上に投影した点を $\vec{q}$ とおく（上図参照）。点 $\vec{q}$ は平面上にあるから \begin{equation} <\vec{n},\vec{q}> = 0 \label{eq_1} \end{equation} が成り立つ。ここで、$<\vec{n},\vec{q}>$ は内積 $\vec{n}^T \cdot \vec{q}$ を表す。これはスカラー量である。また、次式が成り立つ。 \begin{equation} \vec{q}-\vec{l} = t\;(\vec{p}-\vec{l}) \label{eq_2} \end{equation} $t$ は求めるべき実数である。式(\ref{eq_2})の両辺に$\vec{n}$をかけて \begin{equation} <\vec{n},\vec{q}>-<\vec{n},\vec{l}>=t<\vec{n},\vec{p}-\vec{l}> \end{equation} 式(\ref{eq_1})より \begin{equation} -<\vec{n},\vec{l}>=t<\vec{n},\vec{p}-\vec{l}> \end{equation} よって \begin{equation} t=\frac{-<\vec{n},\vec{l}>}{<\vec{n},\vec{p}-\vec{l}>} \label{eq_3} \end{equation} を得る。式(\ref{eq_2})に代入して \begin{eqnarray} \vec{q} &=&\vec{l} - \frac{<\vec{n},\vec{l}>}{<\vec{n},\vec{p}-\vec{l}>} (\vec{p}-\vec{l}) \\ &=&\frac{1}{<\vec{n},\vec{p}-\vec{l}>}\left[ -<\vec{n},\vec{l}>\vec{p}\;+<\vec{n},\vec{p}>\vec{l} \right]\\ &=&\frac{1}{<\vec{n},\vec{l}-\vec{p}>}\left[ <\vec{n},\vec{l}>\vec{p}\;-<\vec{n},\vec{p}>\vec{l} \right]\\ &=&\alpha\left[ <\vec{n},\vec{l}>\vec{p}\;-<\vec{n},\vec{p}>\vec{l} \right] \label{eq_4} \end{eqnarray} $\alpha\equiv 1/<\vec{n},\vec{l}-\vec{p}>$としたが、同次座標系を考えているので $\alpha$ は任意の実数として良い。 ここで、右辺第2項の第 $i$ 成分 $<\vec{n},\vec{p}>l_i$ を考える。 \begin{eqnarray} <\vec{n},\vec{p}>l_i &=&\sum_j\;n_j p_j l_i \\ &=&\sum_j\;l_i n_j p_j \\ &=&\left[\left(\vec{l}\cdot\vec{n}^{T}\right)\vec{p}\right]_i \end{eqnarray} 従って \begin{eqnarray} \vec{q} &=& \alpha\left[<\vec{n},\vec{l}>\vec{p}\;-\left(\vec{l},\vec{n}^{T}\right)\vec{p}\right] \\ &=& \alpha\left[ <\vec{n},\vec{l}>E-\left(\vec{l},\vec{n}^{T}\right)\right] \vec{p} \end{eqnarray} となる。$E$ は単位行列である。以上から \begin{eqnarray} \vec{q}&=&\alpha\;M\;\vec{p}\\ M&\equiv& <\vec{n},\vec{l}>E -\left(\vec{l},\vec{n}^{T}\right) \end{eqnarray} を得る。$M$が影行列である。

実装

影行列は207行目から222行目にかけて計算されている。$\alpha=1$ とした。また、床は $z=0$ ではなく $z=-0.005$ としてある。影をきれいに描画させるためである。

結果

開発環境

　Xcode Version 6.2

疑問

　$\alpha$ を負の数に取ると影が描画されない。何か勘違いしているのか？