RBMについて勉強したので自分用にまとめる．Qiitaの記事を書く練習も兼ねて．
下記のように構成しようと思います．RBMを多層に積み重ねたDBMについても解説しようと思います．
- 導入（この投稿）
- RBM（Restricted Boltzmann Machine）
- DBM（Deep Boltzmann Machine）
RBMはいろんなライブラリに実装されてたり，いたるところでサンプルコードを見かけますが，DBMはHintonさんがGitで公開してるやつくらいしか見かけないので，自分が書いたものをここで紹介できればと思います．
導入
RBM（Restricted Boltzmann Machine）は深層学習の基礎モデルの1つです．観測データを生成する生成モデルを獲得するために使われます．
RBMはBM（Boltzmann Machine）の，BMはMRF（Markov Random Field）の特殊な形ですので，MRFとBMについて簡単に解説します．
マルコフ確率場（Markov Random Field; MRF）
MRFは無向グラフ$G = (\Omega, E)$上に定義される確率モデルです．$i$番目のノードに確率変数$X_{i} \in \{0, 1 \}$を対応させ，その実現値$x_{i}$を並べたベクトル$\bf x$$=\{x_{1}, x_{2}, ..., x_{n}\}$とします．そして，グラフ上に次のようなエネルギー関数を定義します．
\Phi({\bf x}) = - \sum_{\Omega} \phi_{i}(x_{i}) - \sum_{E} \psi_{ij}(x_{i}, x_{j}) \tag{1}
式(1)において関数$\phi_{i}$や$\psi_{ij}$は適当に設計されます．
MRFは式(1)のエネルギー関数を用いた以下の**ボルツマン分布（Boltzmann distribution）**として定義されます．
p({\bf X}={\bf x}) = \frac{1} {Z}\exp(-\Phi({\bf x})) \tag{2}
式(2)において$Z$はすべての実現値$\bf x$に関する確率の総和を1にするための分配関数（partiton function）であり，
Z = \sum_{{\bf x}}\exp(-\Phi({\bf x})) \tag{3}
です．式(2)からわかるように，MRFは$n$個の確率変数の同時分布であり，エネルギー関数の値をより低くする${\bf x}$がより高確率になります．
ボルツマンマシン（Boltzmann Machine; BM）
BMはMRFの特殊な形です．具体的には，式(1)における関数$\phi_{i}$，$\psi_{ij}$をそれぞれ$\phi_{i} = b_{i}x_{i}$，$\psi_{ij} = w_{ij}x_{i}x_{j}$としたものです．これによりエネルギー関数は
\Phi({\bf x}; \theta) = - \sum_{\Omega} b_{i}x_{i} - \sum_{E} w_{ij}x_{i}x_{j} \tag{4}
となります．ここで，パラメータ$\theta = \{{\bf b}, {\bf W}\}$はそれぞれ各ノードに割り当てられたバイアスと各結合に与えられた重みです．式(4)のエネルギー関数を用いてBMは以下のように定式化されます．
p({\bf X}={\bf x}; \theta) = \frac{1} {Z(\theta)}\exp(-\Phi({\bf x}; \theta)) \tag{5}
BMは生成モデルの一種なので，観測されたデータ点を生成する確率を最大にするパラメータ$\theta$を学習により求めます．
確率変数$X_{i}$は観測されたデータ点に直接対応している変数（観測された値を代入する変数）であるため，以降では**可視変数（visible variable）と呼ぶことにします．逆に，直接対応していない変数を隠れ変数（hidden variable）**といいます．また，可視変数であることを明確にしたいときは確率変数$X_{i}$のかわりに$V_{i}$（visibleのv）を用いることにします．
可視変数のみのBMとその学習
まずは隠れ変数をもたない，可視変数のみからなるBMの学習の方法について見てみましょう．すべての変数が可視変数なので，${\bf X} = \{X_{i} | i \in \Omega \}$のかわりに${\bf V} = \{V_{i} | i \in \Omega \}$を用います．
$n$個の要素をもつ観測データ点がそれぞれ独立な分布から$N$個生成されたとします．たとえば，手書き数字の画像とそのラベルからなるデータセットのMNISTを例にとると，$n = 28$pixel$* 28$pixel $= 784$個の要素をもつ観測データ点（手書き文字画像）が$N$個あたえられた状態です．
（MNISTは画像に対応するラベル値もセットになっていますが，議論を簡単にするために，いまはラベル値については考えないものとします．）
BMの学習では一般に最尤推定（maximum likelihood estimation; MLE）が用いられます．最尤推定では，あたえられた観測データ点の集合$\mathcal{D} = \{{\bf v}^{(\mu)} | \mu = 1, 2, ..., N\}$について尤度関数（likelihood function）
l_{\mathcal{D}}(\theta) = \prod_{\mu = 1}^{N} p({\bf V} = {\bf v}^{(\mu)} | \theta) \tag{6}
を定義します．式(6)において$p({\bf V} = {\bf v}^{(\mu)} | \theta)$はパラメータを$\theta$とするBMが観測データ点${\bf v}^{(\mu)}$を実際に生成する確率を表しています．観測データ点はそれぞれ独立な分布から生成されているので，それらの積（式（6））は，今手元にある観測データ点の集合$\mathcal{D}$をBMが生成する確率を表しているといえます．BMの学習はパラメータ$\theta$をうまく調節してやることで，手元にある観測データ点の集合$\mathcal{D}$を最も高い確率で生成する分布を獲得することであると解釈できます．
ここで，式(6)の最大値を求めるにあたり式(6)をパラメータ$\theta$で微分してやる必要がありますが，掛け算の微分は難しいので，その対数をとった対数尤度関数（log-likelihood function）
L_{\mathcal{D}}(\theta) = \ln{L_{\mathcal{D}}(\theta)} = \sum_{\mu = 1}^{N} \ln{p({\bf V} = {\bf v}^{(\mu)} | \theta)} \tag{7}
を考えます．対数尤度関数が最大となる点ではパラメータ$\theta$に関する勾配が$0$になるので，勾配が$0$となる$\theta$の値を求めます．式(7)の対数尤度関数のパラメータ$\theta = \{{\bf b}, {\bf W}\}$に関する勾配はそれぞれ
\frac{\partial L_{\mathcal{D}}(\theta)} {\partial b_{i}} = \sum_{\mu = 1}^{N}v_{i}^{(\mu)} - NE_{p({\bf V} | \theta)}[V_{i}] \tag{8}
\frac{\partial L_{\mathcal{D}}(\theta)} {\partial w_{ij}} = \sum_{\mu = 1}^{N}v_{i}^{(\mu)}v_{j}^{(\mu)} - NE_{p({\bf V} | \theta)}[V_{i}V_{j}] \tag{9}
となります．ここで，$E_{p({\bf V} | \theta)}[\cdot]$はBMに関する期待値であり，
E_{p({\bf V} | \theta)}[f({\bf V})] = \sum_{{\bf v}}f({\bf v})p({\bf v} | \theta) \tag{10}
であたえられます．これは確率変数$\bf V$の可能な実現値$\bf v$すべての組み合わせに関する総和を計算することにより得られます．
式(8)，(9)から，勾配が$0$となる対数尤度関数の最大点では
\frac{1} {N}\sum_{\mu = 1}^{N}v_{i}^{(\mu)} = E_{p({\bf V} | \theta)}[V_{i}] \tag{11}
\frac{1} {N}\sum_{\mu = 1}^{N}v_{i}^{(\mu)}v_{j}^{(\mu)} = E_{p({\bf V} | \theta)}[V_{i}V_{j}] \tag{12}
が成り立ちます．この方程式は**ボルツマンマシンの学習方程式（learning equation of Boltzmann machine）**と呼ばれ，求めたいパラメータ$\theta$はこの連立方程式の解です．式(11)，(12)の左辺は観測データ点集合の平均をとればよいので容易に計算できます．しかし，右辺はBMの実現可能な値すべての組み合わせについて計算する必要があるため，解析的に解くことは容易ではありません．
式(7)の対数尤度関数はパラメータ$\theta$に関して凹関数になっているので，原理的には勾配上昇法を用いて対数尤度関数が最大点となるパラメータを求めることができます．$n$個の確率変数をもつBMの期待値計算は$2^{n}$個の項の和を計算することになるため，組み合わせ爆発を起こしてしまいます．
隠れ変数をもつBMとその学習
次に，観測データ点に直接対応しない隠れ変数をもつBMについて考えます．$(n + m)$個の確率変数のうち，最初の$n$個を可視変数に，残りの$m$個を隠れ変数に割り当てます．また，可視変数を${\bf V} = \{V_{i} | i \in \mathcal{V} \}$，隠れ変数を${\bf H} = \{H_{i} | i \in \mathcal{H} \}$で表し，式(4)，式(5)をそれぞれ以下のように書き換えます．
\Phi({\bf V}, {\bf H}; \theta) = - \sum_{\mathcal{V}} b_{i}v_{i} - \sum_{\mathcal{H}} c_{i}h_{i} - \sum_{E} w_{ij}x_{i}x_{j} \tag{13}
p({\bf V}={\bf v}, {\bf H}={\bf h}; \theta) = \frac{1} {Z(\theta)}\exp(-\Phi({\bf V}, {\bf H}; \theta)) \tag{14}
式(13)において$\bf c$は隠れ変数に対応するノードにあたえられたバイアスです．また，$\sum_{E} w_{ij}x_{i}x_{j}$における$x_{i}$，$x_{j}$は，ノード$i$，$j$が可視変数のノードか隠れ変数のノードかに応じて適当に変換されます．変数の表記が少し変わっただけで，定義そのものは式(4)，式(5)と変わっていないことに注意してください．隠れ変数をもつBMでも最尤推定法を用いて学習しますが．以降の議論は，KLダイバージェンスと呼ばれる概念を導入する必要があるため，ここでは避けたいと思いますが，隠れ変数を導入した場合でも勾配上昇法の計算過程で組み合わせ爆発が起こることに変わりはありません．
次節以降で，勾配上昇法の計算過程で組み合わせ爆発が起こることを防ぐために考えられたRBMや，期待値の近似手法であるコントラスティブ・ダイバージェンス（contarstive divergence; CD）法について解説します．