18 February 2016

1. 導入

これはボルツマンマシン,特に制約付きボルツマンマシンについてまとめたものです(Safariで見た方がいいかも).

近年,画像認識やその他もろもろの分野で深層学習が話題となっています. この研究の皮切りとなったものはボルツマンマシンを基礎とした確率的な深層学習モデルです.ボルツマンマシンは画像のノイズ除去(復元)や最適化問題へのアプローチとして用いられ,ニューラルネットワーク(以後 NN と表記)のモデルの1つであるホップフィールドネットワークを歴史的な母体としています.


2. ホップフィールド ネットワーク

ホップフィールドネットワーク( Hopfield Network: HN )は以下のページで説明します.

E-musu’s Tech Memorandum: Hopfield Network


3. ボルツマンマシン

HNでは極小値に陥ると抜け出すことができませんでした.そこで各ユニットを確率的に動作させることによって,極小値の問題を回避させるように工夫されたものがボルツマンマシン(Boltzmann Machine: BM)です. 各ユニット間の重みは対称\(w_{ij} = w_{ji}\)で,自己結合は無し\(w_{ii} = 0\)と仮定するのはHNと同じで,また,エネルギー関数も同じ形を持ちます.

HNでは\(I_{i}(\mathbf{x})\)の値に従って決定論的な状態変化を行っていました. BMでも各ユニットの出力は0または1を取りますが,ランダムに任意のユニット\(i\)を1つ選択し,以下の確率でユニットの状態を変化させる点で異なります. これにより,エネルギー関数の値が増加する方向へも状態変化が許されることになります.

\begin{equation} Prob\{x_{i}^{‘} = 1\} = \frac{1}{1+\exp \left( \frac{-I_{i}(\mathbf{x})}{T} \right)} \end{equation}

このときTは正のパラメータであり,温度を表したものです. \(T \rightarrow \infty\)ではユニットの出力は\(I_{i}(\mathbf{x})\)の値に関係なく,1になるか0になるかは1/2となります. 対して,\(T \rightarrow 0\)ではHNと同様の状態変化則をとるため,\(I_{i}(\mathbf{x})\)の値に従います. BMではTを大きな値にすることにより,HNの極小値問題から脱出できる可能性を持っています. ただし,温度を高くするとエネルギー最小の状態をとる確率が低くなってしまい,温度を低くするとエネルギー最小の状態になる確率が高くなるが,定常状態になるまで時間がかかるという問題があります。ちなみに後述するBMの学習方法ではT=1で学習を行っています.

ここで説明のため,しいき値\(h_{i}\)をバイアス\(b_{i}\)とします.バイアス\(b_{i}\)及び重み\(w_{ij}\)をまとめてパラメータ\(\Theta\)とすると,BMのエネルギー関数を以下のようになります.

\begin{equation} \Phi \left( \mathbf{x},\Theta \right) = -\sum_{i=1}{b_{i}x_{i}}-\sum_{\left( i,j \right)\in \epsilon }^{}{w_{ij}x_{i}x_{j}} \end{equation}

BMの状態変化則は確率に変化するものであり,ネットワークの状態は1時点前の状態に依存し変化します.つまり,BMの状態変化の過程はマルコフ連鎖をなしていると考えることができます. このとき,更新プロセスを十分繰り返すことによるユニットの状態パターンが以下のボルツマン分布(定常分布)に従うとされています.

\begin{equation} p( \mathbf{x} | \Theta )=\frac{1}{Z( \Theta )}\exp \{ -\Phi ( \mathbf{x},\Theta )\} (Z( \theta )は正規化定数) \label{eq:bol} \end{equation}

この分布はBMのネットワークの状態\(\mathbf{x}\)の確率分布を表し,エネルギー関数が小さい値を取る状態の生起確率が比較的高くなります. ここで,HNでは安定状態をエネルギー関数の極小値を取ることで用いて判断しましたが,これに関してエネルギー関数を評価(コスト)関数と見なすことができ,エネルギーをより小さくする\(\mathbf{x}\)ほどより良い\(\mathbf{x}\)と考えることができます. このことから, \((\ref{eq:bol})\)はある課題に対して望ましい\(\mathbf{x}\)がより高い確率で実現するような確率モデルと言えます.

観測データは個々に出現のし易さが異なるものであるため,いわばそれは確率的に生成されたものと考えることができます. このときBMを確率的な規則に従ってデータを生成するモデルとして捉えると,統計的機械学習により,HNのような連想記憶を含めた様々な応用が期待されます.そこで,BMでは学習を行うことでパラメータ\(\Theta\)(\(w_{ij}\)や\(b_{i}\))を調節し,観測データを生成した真のモデル(仮に存在すると考えて)への近似を目指します.


3.1 学習

BMの学習について,単純に最尤推定を考えてみます. ここでは\(p(\mathbf{x} | \Theta )\)を実際にBMが観測データを生成する確率と解釈するとき,以下の対数尤度関数を最大とするパラメータの値を求めます.

\begin{equation} L\left( \Theta \right)=\prod_{n=1}^{}{p\left( \mathbf{x}_{n}|\Theta \right)} \end{equation}

\begin{equation} \log L\left( \Theta \right)=\sum_{n=1}^{N}{\log p\left( \mathbf{x}_{n} | \Theta \right)=\sum_{n=1}^{N}{\left\{ -\Phi \left( \mathbf{x}_{n},\Theta \right)-\log Z\left( \Theta \right) \right\}}} \end{equation}

この対数尤度関数のパラメータである\(b_{i}\)の勾配は天下り的だが

\begin{eqnarray}\frac{\partial Z\left( \Theta \right)}{\partial b_{i}} &=& \frac{\partial }{\partial b_{i}}\sum_{\mathbf{x}}^{}{\exp \left\{ -\Phi \left( \mathbf{x},\Theta \right) \right\}} = \frac{\partial }{\partial b_{i}}\sum_{\mathbf{x}}^{}{\exp \left( \sum_{i=1}{b_{i}x_{i}} + \sum_{\left( i,j \right)\in \epsilon }^{}{w_{ij}x_{i}x_{j}} \right)}\\ &=& \sum_{\mathbf{x}}^{}{\exp \left( \sum_{i=1}{b_{i}x_{i}} + \sum_{\left( i,j \right)\in \epsilon }^{}{w_{ij}x_{i}x_{j}} \right)}x_{i}\\ &=& \sum_{\mathbf{x}}^{}{\exp \left\{ -\Phi \left( \mathbf{x},\Theta \right) \right\}}x_{i} \end{eqnarray}

より,

\begin{eqnarray}\frac{\partial \log L\left( \Theta \right)}{\partial b_{i}} & = & \sum_{n=1}^{N}{\left\{ x_{ni}-\frac{\partial }{\partial b_{i}}\log Z\left( \Theta \right) \right\}} = \sum_{n=1}^{N}{x_{ni}} - N \frac{\partial }{\partial b_{i}}\log Z\left( \Theta \right)\\ & = & \sum_{n=1}^{N}{x_{ni}}-N\frac{1}{Z\left( \Theta \right)}\frac{\partial Z\left( \Theta \right)}{\partial b_{i}}\\ & = & \sum_{n=1}^{N}{x_{ni}}-N\frac{1}{Z\left( \Theta \right)}\sum_{\mathbf{x}}^{}{\exp \left\{ -\Phi \left( \mathbf{x},\Theta \right) \right\}}x_{i}\\ & = & \sum_{n=1}^{N}{x_{ni}}-N\sum_{\mathbf{x}}^{}{x_{i}\frac{1}{Z\left( \Theta \right)}\exp \left\{ -\Phi \left( \mathbf{x},\Theta \right) \right\}}\\ & = & \sum_{n=1}^{N}{x_{ni}}-N\sum_{\mathbf{x}}^{}{x_{i}p\left( \mathbf{x} | \Theta \right)}\\ & = & \sum_{n=1}^{N}{x_{ni}-N \mbox{E}_{p\left( \mathbf{x} | \Theta \right)}\left[ x_{i} \right]} \end{eqnarray}

となります.

また,パラメータである\(w_{ij}\)の勾配も同様に

\begin{eqnarray}\frac{\partial Z\left( \Theta \right)}{\partial w_{ij}} & = & \frac{\partial }{\partial w_{ij}}\sum_{\mathbf{x}}^{}{\exp \left\{ -\Phi \left( \mathbf{x},\Theta \right) \right\}} = \frac{\partial }{\partial w_{ij}}\sum_{\mathbf{x}}^{}{\exp \left( \sum_{i=1}{b_{i}x_{i}} + \sum_{\left( i,j \right)\in \epsilon }^{}{w_{ij}x_{i}x_{j}} \right)}\\ & = & \sum_{\mathbf{x}}^{}{\exp \left( \sum_{i=1}{b_{i}x_{i}} + \sum_{\left( i,j \right)\in \epsilon }^{}{w_{ij}x_{i}x_{j}} \right)}x_{i}x_{j}\\ & = & \sum_{\mathbf{x}}^{}{\exp \left\{ -\Phi \left( \mathbf{x},\Theta \right) \right\}}x_{i}x_{j} \end{eqnarray}

より,

\begin{eqnarray}\frac{\partial \log L\left( \Theta \right)}{\partial w_{ij}} & = & \sum_{n=1}^{N}{\left\{ x_{ni}x_{nj}-\frac{\partial }{\partial w_{ij}}\log Z\left( \Theta \right) \right\}} = \sum_{n=1}^{N}{x_{ni}x_{nj}} - N \frac{\partial }{\partial w_{ij}}\log Z\left( \Theta \right)\\ & = & \sum_{n=1}^{N}{x_{ni}x_{nj}}-N\frac{1}{Z\left( \Theta \right)}\frac{\partial Z\left( \Theta \right)}{\partial w_{ij}}\\ & = & \sum_{n=1}^{N}{x_{ni}x_{nj}}-N\frac{1}{Z\left( \Theta \right)}\sum_{\mathbf{x}}^{}{\exp \left\{ -\Phi \left( \mathbf{x},\Theta \right) \right\}}x_{i}x_{j}\\ & = & \sum_{n=1}^{N}{x_{ni}x_{nj}}-N\sum_{\mathbf{x}}^{}{x_{i}x_{j}\frac{1}{Z\left( \Theta \right)}\exp \left\{ -\Phi \left( \mathbf{x},\Theta \right) \right\}}\\ & = & \sum_{n=1}^{N}{x_{ni}x_{nj}}-N\sum_{\mathbf{x}}^{}{x_{i}x_{j}p\left( \mathbf{x} | \Theta \right)}\\ & = & \sum_{n=1}^{N}{x_{ni}x_{nj}-N \mbox{E}_{p\left( \mathbf{x} | \Theta \right)}\left[ x_{i}x_{j} \right]} \end{eqnarray}

となります. このとき,それぞれの左辺2項目の\(\mbox{E}_{p ( \mathbf{x} | \Theta )}[ \cdot ]\)は期待値\(\sum_{\mathbf{x}}^{}{\left( \cdot \right)p\left( \mathbf{x} | \Theta \right)}\)を表しているのですが,特に分配関数を用いて,\(\sum_{\mathbf{x}}^{}{=\sum_{x_{1}=0,1}^{}{\sum_{x_{2}=0,1}^{} \cdots {\sum_{x_{N}=0,1}^{}{}}}}\)のように\(\mathbf{x}\)のすべての値の組み合わせ(\(2^{N}\)通り)についての和を計算する必要があるため,ユニットの数の増加と共に計算量が膨大なものとなってしまいます.

なお,対数尤度の最大点,つまり\(\frac{\partial \log L\left( \Theta \right)}{\partial b_{i}} = \frac{\partial \log L\left( \Theta \right)}{\partial w_{ij}} = 0\)のとき,次のような連立方程式が求まりますが,これは解析的に解くことは困難とされているため,\(w_{ij}\)と\(b_{i}\)について,求めた勾配方向に反復更新し,対数尤度関数において勾配を利用した頒布鵜処理を実行することになります.

\begin{equation} \frac{1}{N}\sum_{n=1}^{N}{x_{ni} = \mbox{E}_{p\left( \mathbf{x} | \Theta \right)}\left[ x_{i} \right]} \end{equation}

\begin{equation} \frac{1}{N}\sum_{n=1}^{N}{x_{ni}x_{nj} = \mbox{E}_{p\left( \mathbf{x} | \Theta \right)}\left[ x_{i}x_{j} \right]} \end{equation}

以上から,BMの学習に関しては最尤推定はなかなか骨が折れそうです. 特に分配関数を直接計算しなければならない点が非常にネックです. 最尤推定以外にもBMの計算ではギブスサンプリングを用いた近似法も存在しますが, ギブスサンプリングでもその精度を高めるためには十分な反復回数が必要となるため一般的には計算コストは高めです. しかしながら,データの生成確率\(p(\mathbf{x} | \Theta)\)に従うサンプルを生成することのできる貴重な手段と言えます.


3.2 隠れ変数の導入

BMに限らず,多くの機械学習モデルは(パラメータの値を変化させて)再現できる分布の豊富さが重要となります. なぜなら,パラメータの最適値を求めても生成モデル自体が真の分布とかけ離れたものであった場合,うまく真の分布を近似することは不可能だからです.

BMではエネルギー関数の形で分布が決定します.そのため,より複雑でパラメータ数の多いエネルギー関数を導入することで,様々な未知の分布に対しても,それらをうまく表現することができる柔軟な生成モデルにすることが期待されます. そのアプローチの1つとしてBMに隠れ変数を追加したモデルを考えます. このときユニットは可視変数\(\mathbf{v}\)と隠れ変数\(\mathbf{h}\)で分けることができ, 特に既存の変数を可視変数としているため\(\mathbf{v}\)はデータ\(\mathbf{x}\)と対応しています.

hidden

このときエネルギー関数\(\Phi\left( \mathbf{v}, \mathbf{h}, \Theta \right)\)及びネットワークの状態確率\(p\left( \mathbf{v}, \mathbf{h} | \Theta \right)\)は\(\mathbf{v}\), \(\mathbf{h}\)を順番に並べたベクトル\(\mathbf{z} = \{v_{1},v_{2} \cdots h_{j-1}, h_{j}\}\)を導入すると,以下のように表すことができます. \begin{equation} \Phi\left( \mathbf{v}, \mathbf{h}, \Theta \right) = \Phi \left( \mathbf{z},\Theta \right)\; =\; -\sum_{i=1}^{}{b_{i}z_{i}}-\sum_{\left( i,j \right)\in \epsilon }^{}{w_{ij}z_{i}z_{j}} \end{equation}

\begin{equation} p\left( \mathbf{v}, \mathbf{h} | \Theta \right) = p\left( \mathbf{z} | \Theta \right)=\frac{1}{Z\left( \Theta \right)}\exp \left( -\Phi \left( \mathbf{z},\Theta \right) \right) \end{equation}

ここで,先ほどと同じように最尤推定を考えてみます. しかしながら,データの生成確率\(p\left( \mathbf{v}, \mathbf{h} | \Theta \right)\)は観測データとは無関係である隠れ変数 が含まれているため,\(p\left( \mathbf{v} | \Theta \right)=\sum_{\mathbf{h}}^{}{p\left( \mathbf{v}, \mathbf{h} | \Theta \right)}\)のように隠れ変数\(\mathbf{h}\)を周辺化することにより可視変数\(\mathbf{v}\)のみの分布にします.

すると対数尤度関数は以下のようになります.

\begin{equation} \log L\left( \Theta \right) \propto \frac{1}{N}\sum_{n=1}^{}{\left[ \log \sum_{\mathbf{h}}^{}{\exp \left\{ -\Phi \left( \mathbf{v}_{n},\mathbf{h},\Theta \right) \right\}-\log Z\left( \Theta \right)} \right]} \end{equation}

また,パラメータである\(b_{i}\)及び\(w_{ij}\)の勾配は以下のように求まります.

\begin{equation} \frac{\partial \log L\left( \Theta \right)}{\partial b_{i}}=\sum_{\mathbf{v}, \mathbf{h}}^{}{x_{ni} p\left( \mathbf{h}| \mathbf{v}, \Theta \right) -N \mbox{E}_{p\left( \mathbf{v}, \mathbf{h} | \Theta \right)}\left[ x_{i} \right]} \end{equation}

\begin{equation} \frac{\partial \log L\left( \Theta \right)}{\partial w_{ij}}=\sum_{\mathbf{v}, \mathbf{h}}^{}{x_{ni}x_{nj} p\left( \mathbf{h}| \mathbf{v}, \Theta \right) -N \mbox{E}_{p\left( \mathbf{v}, \mathbf{h} | \Theta \right)}\left[ x_{i}x_{j} \right]} \end{equation}

このとき求めたそれぞれ勾配について1項目,2項目共に分配関数を含んでいるためさらに計算が困難なものとなっています.ここで任意の分布を表現しようと隠れ変数を導入し,エネルギー関数を複雑にすることは膨大な計算量とのトレードオフな関係と言えます.


4. 制約付きボルツマンマシン

隠れ変数を導入したBMではエネルギー関数を複雑にすることでモデルの表現能力向上を目指しましたが,求めた勾配に含まれる分配関数により膨大な量の計算をする必要がありました. そこで提案されたのが制約付きボルツマンマシン(Restricted Beltzmann Machine :RBM )です. RBMは完全2部グラフ上に定義されたBMであり,可視変数\(\mathbf{v}\)のみの層(可視層)と隠れ変数\(\mathbf{h}\)のみの層(隠れ層)の2つの層に分かれています.

hidden

このときエネルギー関数\(\Phi \left( \mathbf{v}, \mathbf{h}, \Theta \right)\)及びネットワークの状態確率\(p\left( \mathbf{v}, \mathbf{h} | \Theta \right)\)は可視変数のバイアスを\(b_{i}\),隠れ変数のバイアスを\(c_{j}\)とすると以下のようになります. \begin{equation} \Phi \left( \mathbf{v}, \mathbf{h}, \Theta \right) = -\sum_{i}{b_{i}v_{i}}-\sum_{j}{c_{j}h_{j}}-\sum_{\left( i,j \right)\in \epsilon }^{}{w_{ij}v_{i}h_{j}} \end{equation}

\begin{equation} p\left( \mathbf{v}, \mathbf{h} | \Theta \right) =\frac{1}{Z\left( \Theta \right)}\exp \left( -\Phi \left( \mathbf{v}, \mathbf{h},\Theta \right) \right) \end{equation}

RBMでは同じ層のユニットを結ぶリンクが存在しないことが制約であり,計算量の削減へとつながる重要なポイントとなります. この制約により可視変数と隠れ変数は条件付き独立性の性質が持つことになります. これはグラフィカルモデルにおいて,片方の層の値を固定するともう片方の層の確率確率変数が互いに独立になることを表し,次の式が成り立ちます.

\begin{equation} p\left( \mathbf{h} | \mathbf{v},\Theta \right)=\prod_{j=1}^{}{p\left( h_{j} | \mathbf{v}, \Theta \right)} \label{eq:j1} \end{equation}

\begin{equation} p\left( \mathbf{v} | \mathbf{h},\Theta \right)=\prod_{i=1}^{}{p\left( v_{i} | \mathbf{h}, \Theta \right)} \label{eq:j2} \end{equation}

例えば\(p\left( \mathbf{h} | \mathbf{v},\Theta \right)\)は具体的には \begin{eqnarray} p ( \mathbf{h} | \mathbf{v}, \Theta ) & = & \frac{ p ( \mathbf{v}, \mathbf{h} | \Theta )}{\sum_{\mathbf{h}}{p(\mathbf{v}, \mathbf{h} | \Theta)} }\\ & = & \frac{\frac{1}{Z\left( \Theta \right)}\prod_{i}^{}{\exp \left( b_{i}v_{i} \right)\prod_{j}^{}{\exp \left( c_{j}h_{j}+\sum_{i}^{}{w_{ij}v_{i}h_{j}} \right)}}}{\frac{1}{Z\left( \Theta \right)}\prod_{i}^{}{\exp \left( b_{i}v_{i} \right)\prod_{j}^{}{\sum_{h_{j}}^{}{\exp \left( c_{j}h_{j}+\sum_{i}^{}{w_{ij}v_{i}h_{j}} \right)}}}}\\ & = & \prod_{j}^{}{\frac{\exp \left( c_{j}h_{j}+\sum_{i}^{}{w_{ij}v_{i}h_{j}} \right)}{\sum_{h_{j}}^{}{\exp \left( c_{j}h_{j}+\sum_{i}^{}{w_{ij}v_{i}h_{j}} \right)}}}\\ & = & \prod_{j}^{}{\frac{\exp \{ h_{j} ( c_{j}+\sum_{i}^{}{w_{ij}v_{i}} ) \}}{ 1 + \exp \left( c_{j}+\sum_{i}^{}{w_{ij}v_{i}} \right) }} = \prod_{j}^{}{p\left( h_{j} | \mathbf{v}, \Theta \right)} \end{eqnarray}

と求まります.このとき可視層\(\mathbf{v}\)からの入力を固定したとき\(h_{j} = 1\)となる条件付き確率\(p(h_{j} = 1 | \mathbf{v} , \Theta)\)は\(\sigma \left( x \right)=\frac{1}{1+\exp \left( -x \right)}\)とすると \begin{eqnarray} p(h_{j} = 1 | \mathbf{v} , \Theta) & = & \frac{\exp \left( c_{j}+\sum_{i}^{}{w_{ij}v_{i}} \right)}{1+\exp \left( c_{j}+\sum_{i}^{}{w_{ij}v_{i}} \right)}\\ & = & \frac{1}{1+\exp \{ -\left( c_{j}+\sum_{i}^{}{w_{ij}v_{i}} \right) \}}\\ & = & \sigma (c_{j} + \sum_{i}{w_{ij}v_{i}}) \end{eqnarray}

となります.また対称性から\(p\left( \mathbf{v} | \mathbf{h},\Theta \right)\)及び隠れ層からの入力を固定したとき\(v_{i} = 1\)となる条件付き確率\(p(v_{i}=1 | \mathbf{h}, \Theta) \)はそれぞれ次のようになります. \begin{eqnarray} p ( \mathbf{v} | \mathbf{h}, \Theta ) & = & \prod_{i}^{}{\frac{\exp \{ v_{i} ( b_{j}+\sum_{j}^{}{w_{ij}h_{j}} ) \}}{ 1 + \exp \left( b_{i}+\sum_{j}^{}{w_{ij}h_{j}} \right) }} = \prod_{i}^{}{p\left( v_{i} | \mathbf{h}, \Theta \right)} \end{eqnarray}

\begin{equation} p(v_{i}=1 | \mathbf{h}, \Theta) = \sigma (b_{i} + \sum_{j}{w_{ij}h_{j}}) \end{equation}

RBMではこれらの条件付き確率により例えば可視層の(ユニットの)値を固定し,対応する隠れユニットにサンプル値\( \in \{0, 1\}\)を\(p(h_{j} | \mathbf{v} , \Theta)\)に従って得ることができます. また,\((\ref{eq:j1})\), \((\ref{eq:j2})\)より各層の確率は同時に計算できるため, それぞれの条件付き分布\(p\left( \mathbf{v} | \mathbf{h},\Theta \right), p\left( \mathbf{h} | \mathbf{v},\Theta \right)\)を用いることで層ごとに交互のサンプリングを行うことできます. これを繰り返すことはネットワークの状態確率\(p(\mathbf{v}, \mathbf{h} | \Theta )\)に従うサンプルを生成することに繋がります.

hidden


4.1 学習

RBMにおける対数尤度関数は以下のようになります.

\begin{eqnarray}\log L\left( \Theta \right) & = & \sum_{n=1}^{}{\log \sum_{\mathbf{h}}^{}{p\left(\mathbf{v}_n,\mathbf{h} | \Theta \right)}}\\ & = & \sum_{n=1}^{}{\log \sum_{\mathbf{h}}^{}{\exp \{ -\Phi \left( \mathbf{v}_{n},\mathbf{h} \right) \}-\log Z\left( \Theta \right)}}\\ & = & \sum_{n=1}^{}{\log \sum_{\mathbf{h}}^{}{\exp \{ -\Phi \left( \mathbf{v}_{n},\mathbf{h} \right) \}}} - N \log \sum_{\mathbf{v},\mathbf{h}}^{}{\exp \{ -\Phi \left( \mathbf{v},\mathbf{h} \right) \}} \end{eqnarray}

このとき勾配は

\begin{eqnarray}\frac{\partial \log L\left( \Theta \right)}{\partial \Theta} & = & \sum_{n=1}^{}{\frac{\partial}{\partial \Theta}\left[\log \sum_{\mathbf{h}}^{}{\exp \{ -\Phi \left( \mathbf{v}_n,\mathbf{h} \right) \}}\right]} - N\frac{\partial}{\partial \Theta}\left[\log \sum_{\mathbf{v},\mathbf{h}}^{}{\exp \{ -\Phi \left( \mathbf{v},\mathbf{h} \right) \}}\right]\\ & = & - \sum_{n=1}^{}{\frac{\sum_{\mathbf{h}}^{}{\left[ \exp \left\{ -\Phi \left( \mathbf{v}_{n},\mathbf{h},\Theta \right) \right\}\frac{\partial \Phi \left( \mathbf{v}_{n},\mathbf{h},\Theta \right)}{\partial \Theta } \right]}}{\sum_{\mathbf{h}}^{}{\exp \left\{ -\Phi \left( \mathbf{v}_{n},\mathbf{h},\Theta \right) \right\}}}}+N\frac{\sum_{\mathbf{v},\mathbf{h}}^{}{\left[ \exp \left\{ -\Phi \left( \mathbf{v},\mathbf{h},\Theta \right) \right\}\frac{\partial \Phi \left( \mathbf{v},\mathbf{h},\Theta \right)}{\partial \Theta } \right]}}{\sum_{\mathbf{v},\mathbf{h}}^{}{\exp \left\{ -\Phi \left( \mathbf{v},\mathbf{h},\Theta \right) \right\}}}\\ & = & - \sum_{n=1}^{}{\sum_{\mathbf{h}}^{}{\frac{\partial \Phi \left( \mathbf{v}_n, \mathbf{h} | \Theta \right)}{\partial \Theta }p\left( \mathbf{h} | \mathbf{v}_n, \Theta \right)}} + N \sum_{\mathbf{v},\mathbf{h}}^{}{\frac{\partial \Phi \left( \mathbf{v},\mathbf{h},\Theta \right)}{\partial \Theta }p\left( \mathbf{v},\mathbf{h}|\Theta \right)} \end{eqnarray}

となります. そして\(b_{i}, w_{ij}, c_{j}\)に注目したとき,次のように求まります.

\begin{eqnarray}\frac{\partial \log L\left( \Theta \right)}{\partial b_{i}} & = & - \sum_{n=1}^{}{\sum_{\mathbf{h}}^{}{\frac{\partial \Phi \left( \mathbf{v}_n, \mathbf{h} | \Theta \right)}{\partial b_{i} }p\left( \mathbf{h} | \mathbf{v}_n, \Theta \right)}} + N \sum_{\mathbf{v},\mathbf{h}}^{}{\frac{\partial \Phi \left( \mathbf{v},\mathbf{h},\Theta \right)}{\partial b_{i} }p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\\ & = & \sum_{n=1}^{}{\sum_{\mathbf{h}}^{}{v_{ni}p\left( \mathbf{h} | \mathbf{v}_n, \Theta \right)}} - N \sum_{\mathbf{v},\mathbf{h}}^{}{v_{i}p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\\ & = & \sum_{n=1}^{}{v_{ni}} - N \mbox{E}_{p\left( \mathbf{v},\mathbf{h} | \Theta \right)}\left[ v_{i} \right] \label{eq:gb} \end{eqnarray}

\begin{eqnarray}\frac{\partial \log L\left( \Theta \right)}{\partial w_{ij}} & = & - \sum_{n=1}^{}{\sum_{\mathbf{h}}^{}{\frac{\partial \Phi \left( \mathbf{v}_n, \mathbf{h} | \Theta \right)}{\partial w_{ij} }p\left( \mathbf{h} | \mathbf{v}_n, \Theta \right)}} + N \sum_{\mathbf{v},\mathbf{h}}^{}{\frac{\partial \Phi \left( \mathbf{v},\mathbf{h},\Theta \right)}{\partial w_{ij} }p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\\ & = & \sum_{n=1}^{}{\sum_{\mathbf{h}}^{}{v_{ni}h_{j}p\left( \mathbf{h} | \mathbf{v}_n, \Theta \right)}} - N \sum_{\mathbf{v},\mathbf{h}}^{}{v_{i}h_{j}p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\\ & = & \sum_{n=1}^{}{v_{ni}p\left( h_{j} = 1 | \mathbf{v}_n, \Theta \right)} - N \mbox{E}_{p\left( \mathbf{v},\mathbf{h} | \Theta \right)}\left[ v_{i}h_{j} \right] \label{eq:gw} \end{eqnarray}

\begin{eqnarray}\frac{\partial \log L\left( \Theta \right)}{\partial c_{j}} & = & - \sum_{n=1}^{}{\sum_{\mathbf{h}}^{}{\frac{\partial \Phi \left( \mathbf{v}_n, \mathbf{h} | \Theta \right)}{\partial c_{j} }p\left( \mathbf{h} | \mathbf{v}_n, \Theta \right)}} + N \sum_{\mathbf{v},\mathbf{h}}^{}{\frac{\partial \Phi \left( \mathbf{v},\mathbf{h},\Theta \right)}{\partial c_{j} }p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\\ & = & \sum_{n=1}^{}{\sum_{\mathbf{h}}^{}{h_{j}p\left( \mathbf{h} | \mathbf{v}_n, \Theta \right)}} - N \sum_{\mathbf{v},\mathbf{h}}^{}{h_{j}p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\\ & = & \sum_{n=1}^{}{p\left( h_{j}=1 | \mathbf{v}_{n}, \Theta \right) -N \mbox{E}_{p\left( \mathbf{v},\mathbf{h} |\Theta \right)}\left[ h_{j} \right]} \label{eq:gc} \end{eqnarray}

このとき完全2部グラフ上に定義されていることにより,隠れ変数を導入しただけの場合と比較して1項目は容易に計算することができ,計算量は抑えられていることが分かります. RBMの2項目は例によって\(\mbox{E}_{p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\left[ \cdot \right] = \sum_{\mathbf{v},\mathbf{h}}^{}{\left( \cdot \right)p\left( \mathbf{v},\mathbf{h} | \Theta \right)}\)と分配関数を2重に含むものとなり計算は困難です.

ここで,RBMでは(条件付き確率が容易に計算可能であるから)層ごとの交互のサンプリングを行うことは簡単に行うことができるため,\(p(\mathbf{v}, \mathbf{h} | \Theta )\)に従うサンプルを生成することができます. ギブスサンプリングを用いた計算方法が考えてみます. RBMのNNではギブスサンプリングを十分な回数の繰り返すことによって,任意の初期分布によらずある定常分布に収束します.しかしこの場合マルコフ連鎖が定常分布に収束するまでサンプリングを繰り返さなければ,定常分布からのサンプルを得ることができません.また,定常分布への収束はパラメータを更新するたびに必要となるため,パラメータを収束させるまでに必要な計算時間が長くなりすぎる問題が存在します. そこでより計算が早く,精度も(それなりに)良い方法が求められます.


4.2 コントラスティブダイバージェンス

RBMではパラメータ推定手法としてコントラスティブダイバージェンス法(Contrastive Divergence Method: CD法)と呼ばれる手法が用いられます. これは収束するまで十分な回数の繰り返す必要のあるギブスサンプリングをわずかk回のみ行うという単純な方法です.ギブスサンプリングの回数がk回であることを明示するためにCD-k法と表記される場合があります.

CD法はまず\(\mathbf{v}^{(0)}\)に観測データをセットすることで初期化を行います. これにより可視ユニットにはそれぞれ0または1が出力値として設定されていることになります. このとき,\(\mathbf{v}^{(0)}\)を0ステップ目の可視層の状態を表すものと考えます. 次に条件付き確率\(p(\mathbf{h}^{(1)}=1|\mathbf{v}^{(0)}, \Theta)\)を計算し,その確率に従って\(\mathbf{h}^{(1)}\)についてサンプルを生成します. これにより隠れユニットに0または1が出力値として設定されることになります. 今度は反対方向に\(p(\mathbf{v}^{(1)}=1|\mathbf{h}^{(1)}, \Theta)\)を計算し,その確率に従って \(\mathbf{v}^{(1)}\)についてサンプルを生成し,可視ユニットの値を書き換えます. ここで\(\mathbf{v}^{(0)}\)から\(\mathbf{v}^{(1)}\)へ可視層の状態の更新がされたと言えます.

hidden

このときの\(\mathbf{v}^{(0)} \rightarrow \mathbf{h}^{(1)} \rightarrow \mathbf{v}^{(1)} \)のような一連の流れをギブスサンプリングの1ステップとしてGibbs Stepと呼びます. このGibbs Stepをk回続けた結果に得られるサンプル\(\mathbf{v}^{(k)}, \mathbf{h}^{(k)}\)は\(k \rightarrow \infty\)ならば\(p(\mathbf{v}, \mathbf{h}|\Theta)\)からサンプリングされたもの同等と考えられています.

RBMにおける対数尤度関数について,勾配 \((\ref{eq:gb})\), \((\ref{eq:gw})\), \((\ref{eq:gc})\) はそれぞれ2項目に\(\mbox{E}_{p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\left[ \cdot \right] = \sum_{\mathbf{v},\mathbf{h}}^{}{\left( \cdot \right)p\left( \mathbf{v},\mathbf{h} | \Theta \right)}\)となる可視変数と隠れ変数の同時確率による期待値が存在していましたが, CD法はGibbs Stepにより得られたサンプル\(\mathbf{v}^{(k)} = \{v_{1}^{(k)}, \cdots , v_{N}^{(k)}\}\)を用いて同時確率を次の式で近似します.

\begin{equation} p\left( \mathbf{v},\mathbf{h} | \Theta \right) \approx q(\mathbf{v}^{(k)})p\left( \mathbf{h}|\mathbf{v}^{(k)}, \Theta \right) \end{equation}

このとき,\(q(\mathbf{v})\)は経験分布と呼ばれ,観測データから出現確率を決める確率分布を表すものです. これを用いると観測データの平均は経験分布に関する期待値として以下のように表すことができます.

\begin{equation} \frac{1}{N}\sum_{n=1}^{}{x_{ni}} = \sum_{\mathbf{v}}^{}{x_{i}q\left( \mathbf{v} \right)} \end{equation}

CD法は通常k=1が用いられており,とても荒い近似をすることになりますがうまく学習できることが分かっています.経験的に. なぜこれでもうまく学習できるのかということについては複数の視点からの解析をなされています.

以上を踏まえた上で例えば \((\ref{eq:gb})\)の\(b_{i}\)におけるの勾配の2項目は次のように求めることができます.

\begin{eqnarray}N \mbox{E}_{p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\left[ v_{i}^{(k)} \right] & = & N\sum_{\mathbf{v},\mathbf{h}}^{}{v_{i}^{(k)}p\left( \mathbf{v},\mathbf{h} | \Theta \right)} = N\sum_{\mathbf{v},\mathbf{h}}^{}{v_{i}^{(k)}q(\mathbf{v}^{(k)})p\left( \mathbf{h} | \mathbf{v}^{(k)}, \Theta \right)}\\ & = & N\sum_{\mathbf{v}}^{}{v_{i}^{(k)}q(\mathbf{v}^{(k)})}\\ & = & \sum_{n=1}{v_{ni}^{(k)}} \end{eqnarray}

また,\((\ref{eq:gw})\)の\(w_{ij}\)におけるの勾配の2項目は次のように求めることができます.

\begin{eqnarray}N \mbox{E}_{p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\left[ v_{i}^{(k)}h_{j} \right] & = & N\sum_{\mathbf{v},\mathbf{h}}^{}{v_{i}^{(k)}h_{j}p\left( \mathbf{v},\mathbf{h} | \Theta \right)} = N\sum_{\mathbf{v},\mathbf{h}}^{}{v_{i}^{(k)}h_{j}q(\mathbf{v}^{(k)})p\left( \mathbf{h} | \mathbf{v}^{(k)}, \Theta \right)}\\ & = & N \sum_{\mathbf{h}}^{}{\left\{ \frac{1}{N}\sum_{n=1}^{}{v_{ni}^{(k)}h_{j}p\left( \mathbf{h} | \mathbf{v}_{n}^{(k)}, \Theta \right)} \right\}}\\ & = & \sum_{n=1}^{}{v_{ni}^{(k)}\sum_{\mathbf{h}}^{}{h_{j}p\left( \mathbf{h} | \mathbf{v}_{n}^{(k)}, \Theta \right)}}\\ & = & \sum_{n=1}^{}{v_{ni}^{(k)}p\left( h_{j} = 1 | \mathbf{v}_{n}^{(k)}, \Theta \right)} \end{eqnarray}

最後に \((\ref{eq:gc})\)の\(c_{j}\)におけるの勾配の2項目は

\begin{eqnarray}N \mbox{E}_{p\left( \mathbf{v},\mathbf{h}|\Theta \right)}\left[ h_{j} \right] & = & N\sum_{\mathbf{v},\mathbf{h}}^{}{h_{j}^{(k)}p\left( \mathbf{v},\mathbf{h} | \Theta \right)} = N\sum_{\mathbf{v},\mathbf{h}}^{}{h_{j}q(\mathbf{v}^{(k)})p\left( \mathbf{h} | \mathbf{v}^{(k)}, \Theta \right)}\\ & = & N \sum_{\mathbf{h}}^{}{\left\{ \frac{1}{N}\sum_{n=1}^{}{h_{j}p\left( \mathbf{h} | \mathbf{v}_{n}^{(k)}, \Theta \right)} \right\}}\\ & = & \sum_{n=1}^{}{\sum_{\mathbf{h}}^{}{h_{j}p\left( \mathbf{h} | \mathbf{v}_{n}^{(k)}, \Theta \right)}}\\ & = & \sum_{n=1}^{}{p\left( h_{j} = 1 | \mathbf{v}_{n}^{(k)}, \Theta \right)} \end{eqnarray}

と求めることができます. 以上からRBMにおける対数尤度関数のパラメータ\(b_{i}, c_{j}, w_{ij}\)の勾配は次のよう求まります.

\begin{eqnarray} \frac{\partial \log L\left( \Theta \right)}{\partial b_{i}} & = & \sum_{n=1}^{}{v_{ni} - \sum_{n=1}^{}{v_{ni}^{(k)}}} \\ & = & v_{i} - v_{i}^{(k)} \end{eqnarray}

\begin{eqnarray} \frac{\partial \log L\left( \Theta \right)}{\partial w_{ij}} & = & \sum_{n=1}^{}{v_{ni}p\left( h_{j}=1 | \mathbf{v}_{n}, \Theta \right) - \sum_{n=1}^{}{v_{ni}^{(k)}p\left( h_{j} = 1 | \mathbf{v}_{n}^{(k)}, \Theta \right)}} \\ & = & v_{i}^{(0)}p(h_{j}=1|\mathbf{v}, \Theta) - v_{i}^{(k)}p(h_{j}=1|\mathbf{v}^{(k)}, \Theta) \end{eqnarray}

\begin{eqnarray} \frac{\partial \log L\left( \Theta \right)}{\partial c_{j}} & = & \sum_{n=1}^{}{p\left( h_{j}=1 | \mathbf{v}_{n}, \Theta \right) - \sum_{n=1}^{}{p\left( h_{j} = 1 | \mathbf{v}_{n}^{(k)}, \Theta \right)}} \\ & = & p(h_{j}=1|\mathbf{v}, \Theta) - p(h_{j}=1|\mathbf{v}^{(k)}, \Theta) \end{eqnarray}

また,学習係数を\(\eta\)とすると更新式は以下のようになります.

\begin{equation} b_{i} \leftarrow b_{i} + \eta \{ v_{i}^{(0)} - v_{i}^{(k)}\} \end{equation}

\begin{equation} w_{ij} \leftarrow w_{ij} + \eta \{ v_{i}^{(0)}p(h_{j}=1|\mathbf{v}, \Theta) - v_{i}^{(k)}p(h_{j}=1|\mathbf{v}^{(k)}, \Theta) \} \end{equation}

\begin{equation} c_{j} \leftarrow c_{j} + \eta \{ p(h_{j}=1|\mathbf{v}^{(0)}, \Theta) - p(h_{j}=1|\mathbf{v}^{(k)}, \Theta) \} \end{equation}

CD法の処理手順を以下に示します.

hidden

ちなみにRBMの学習方法としてCD法をさらに発展させた持続的CD法(Persistent CD: PCD)も存在します(説明は省きます).

最後に,ウソを言ってたらスイマセン・・・


4.3 Pythonでの実装

RBM.py

# -*- coding: utf-8 -*-
import numpy as np

def sigmoid(x):
    return 1. / (1 + np.exp(-x))

class RBM(object):

    def __init__(self, input, n_visible, n_hidden, np_rng):
        self.n_visible = n_visible
        self.n_hidden = n_hidden

        a = 1. / n_visible
        self.W = np.array(np_rng.uniform(low=-a, high=a, size=(n_visible, n_hidden)))

        self.np_rng = np_rng
        self.input = input
        self.hbias = np.zeros(n_hidden)
        self.vbias = np.zeros(n_visible)

    def cd(self, eta=0.1, k=1):
        ph_mean, ph_sample = self.sample_h(self.input)
        chain_start = ph_sample

        for step in range(k):
            if step == 0:
                nv_means, nv_samples, nh_means, nh_samples = self.gibbs_samp(chain_start) 
            else:
                nv_means, nv_samples, nh_means, nh_samples = self.gibbs_samp(nh_samples)

        self.W     += eta * (np.dot(self.input.T, ph_mean) - np.dot(nv_samples.T, nh_means))
        self.vbias += eta * np.mean(self.input - nv_samples, axis=0)
        self.hbias += eta * np.mean(ph_mean - nh_means, axis=0)

    def sample_h(self, v0_sample):
        pre_sigmoid_activation = np.dot(v0_sample, self.W) + self.hbias
        h1_mean = sigmoid(pre_sigmoid_activation)
        h1_sample = self.np_rng.binomial(size=h1_mean.shape, n=1, p=h1_mean)
        return [h1_mean, h1_sample]


    def sample_v(self, h0_sample):
        pre_sigmoid_activation = np.dot(h0_sample, self.W.T) + self.vbias
        v1_mean = sigmoid(pre_sigmoid_activation)
        v1_sample = self.np_rng.binomial(size=v1_mean.shape, n=1, p=v1_mean)
        return [v1_mean, v1_sample]

    def gibbs_samp(self, h0_sample):
        v1_mean, v1_sample = self.sample_v(h0_sample)
        h1_mean, h1_sample = self.sample_h(v1_sample)
        return [v1_mean, v1_sample, h1_mean, h1_sample]

    def get_cross_entropy(self):
        pre_sigmoid_activation_h = np.dot(self.input, self.W) + self.hbias
        sigmoid_activation_h = sigmoid(pre_sigmoid_activation_h)
        pre_sigmoid_activation_v = np.dot(sigmoid_activation_h, self.W.T) + self.vbias
        sigmoid_activation_v = sigmoid(pre_sigmoid_activation_v)
        cross_entropy =  - np.mean(
            np.sum(self.input * np.log(sigmoid_activation_v) + (1 - self.input) * np.log(1 - sigmoid_activation_v), axis=1))
        return cross_entropy

    def reconstruct(self, v):
        h = sigmoid(np.dot(v, self.W) + self.hbias)
        mRecon = sigmoid(np.dot(h, self.W.T) + self.vbias)
        sRecon = self.np_rng.binomial(size=mRecon.shape, n=1, p=mRecon)
        return mRecon, sRecon


data = np.array([[1,1,1,0,0,0],
                 [0,0,1,1,0,0],
                 [0,0,0,1,1,1]])

rbm = RBM(input=data, n_visible=len(data[0]), n_hidden=4, np_rng=np.random.RandomState(123))

# train
training_epochs=1000
learning_rate=0.1
k=1
for epoch in range(training_epochs):
    rbm.cd(eta=learning_rate, k=k)
    cost = rbm.get_cross_entropy()
    print 'Training epoch: %d, cost: ' % epoch, cost

test = np.array([[0,0,1,1,1,1],
                 [1,1,0,0,0,1]])
mRecon, sRecon = rbm.reconstruct(test)
print sRecon

出力


.
.
Training epoch: 996, cost:  0.0794359546865
Training epoch: 997, cost:  0.0784908779092
Training epoch: 998, cost:  0.0784908779092
Training epoch: 999, cost:  0.0784908779092
[[0 0 0 1 1 1]
 [1 1 1 0 0 0]]

プログラムは想起したパターンを出力しています. 0 : □,1 : ■としたとき,次の3つのパターンを学習させました.

■■■□□□ | □□■■□□ | □□□■■■

学習後,新たにパターンを与え動かすと学習したパターンを想起します. 以下は入力パターン → 想起したパターンと対応しています.

□□■■■■ → □□□■■■ | ■■□□□■ → ■■■□□□



Categories

Tags



blog comments powered by Disqus