ボルツマンマシンの定義

ボルツマンマシン（Boltzmann Machine; BM）では、状態ベクトルを \( \vec{x}=(x_1,\dots,x_N) \) とし、エネルギー関数を次のように定義します。

$$ E_\theta(\vec{x}) = -\sum_i b_i x_i -\sum_{i<j} w_{ij} x_i x_j \tag{1} $$

ここで、\(x_i \in \{0,1\}\)、\(b_i\) はバイアス、\(w_{ij}\) は結合重みを表します。これ以降は、パラメータをまとめて \(\theta=\{\vec{b}, \vec{w}\}\) で表記します。エネルギー関数に基づき、確率分布（ボルツマン分布）を導入します。定義は以下の通りです。

$$ P_\theta(\vec{x}) = \frac{1}{Z_\theta}\exp(-E_\theta(\vec{x})) \tag{2} $$

$$ Z_\theta = \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) \tag{3} $$

KL情報量の最小化

モデル分布をデータ分布に近づけるため、KL情報量を利用します。

$$ \begin{aligned} D_{KL}(P_{\text{data}}||P_\theta) &= \sum_{\vec{x}} P_{\text{data}}(\vec{x}) \log \frac{P_{\text{data}}(\vec{x})}{P_\theta(\vec{x})} \\ &= \sum_{\vec{x}} P_{\text{data}}(\vec{x}) \left( \log P_{\text{data}}(\vec{x}) – \log P_\theta(\vec{x}) \right) \\ &= \sum_{\vec{x}} P_{\text{data}}(\vec{x})\log P_{\text{data}}(\vec{x}) – \sum_{\vec{x}} P_{\text{data}}(\vec{x})\log P_\theta(\vec{x}) \quad (\text{∵ } \log \tfrac{a}{b}=\log a-\log b) \end{aligned} $$

第一項はモデルパラメータ \(\theta\) に依存しない定数です。したがって、KL情報量の最小化は以下の量の最大化と等価です。

$$ \max_{\theta} \quad f(\theta) := \sum_{\vec{x}} P_{\text{data}}(\vec{x})\log P_\theta(\vec{x}) \tag{6} $$

更新式の導出

\(f(\theta)\) を最大化するため、勾配上昇法を用います。したがって、以下の勾配を計算します。

$$ \nabla f(\theta) = \sum_{\vec{x}} P_{\text{data}}(\vec{x})\nabla \log P_\theta(\vec{x}) \tag{7} $$

まず、\(\nabla \log P_\theta(\vec{x})\) を計算します。

$$ \begin{aligned} \nabla \log P_\theta(\vec{x}) &= \nabla \log \frac{\exp(-E_\theta(\vec{x}))} {\sum_{\vec{x}}\exp(-E_\theta(\vec{x}))} \\ &= \nabla \left( \log \exp(-E_\theta(\vec{x})) – \log \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) \right) \quad (\text{∵ } \log \tfrac{a}{b}=\log a-\log b) \\ &= -\nabla E_\theta(\vec{x}) – \nabla \log \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) \quad (\text{∵ } \log e^a=a) \end{aligned} \tag{8} $$

ここで、第2項を式変形すると以下のようになります。

$$ \begin{aligned} \nabla \log \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) &= \frac{ \nabla \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) } { \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) } \quad (\text{∵ } \nabla\log g=\tfrac{\nabla g}{g}) \\ &= \frac{ \sum_{\vec{x}} \nabla \exp(-E_\theta(\vec{x})) } { \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) } \\ &= \frac{ \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) (-\nabla E_\theta(\vec{x})) } { \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) } \quad (\text{∵ } \nabla e^f=e^f\nabla f) \\ &= – \frac{ \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) \nabla E_\theta(\vec{x}) } { \sum_{\vec{x}} \exp(-E_\theta(\vec{x})) } \\ &= – \sum_{\vec{x}} P_\theta(\vec{x}) \nabla E_\theta(\vec{x}) \end{aligned} $$

これより、式(8)は以下のようになります。

$$ \nabla \log P_\theta(\vec{x}) = -\nabla E_\theta(\vec{x}) + \sum_{\vec{x}^{\,\prime}} P_\theta(\vec{x}^{\,\prime}) \nabla E_\theta(\vec{x}^{\,\prime}) $$

これを式(7)に代入すると、次の勾配が得られます。

$$ \begin{aligned} \nabla f(\theta) &= \sum_{\vec{x}} P_{\text{data}}(\vec{x}) \left( -\nabla E_\theta(\vec{x}) + \sum_{\vec{x}^{\,\prime}} P_\theta(\vec{x}^{\,\prime}) \nabla E_\theta(\vec{x}^{\,\prime}) \right) \\ &= – \sum_{\vec{x}} P_{\text{data}}(\vec{x}) \nabla E_\theta(\vec{x}) + \sum_{\vec{x}} P_{\text{data}}(\vec{x}) \sum_{\vec{x}^{\,\prime}} P_\theta(\vec{x}^{\,\prime}) \nabla E_\theta(\vec{x}^{\,\prime}) \\ &= – \sum_{\vec{x}} P_{\text{data}}(\vec{x}) \nabla E_\theta(\vec{x}) + \sum_{\vec{x}^{\,\prime}} P_\theta(\vec{x}^{\,\prime}) \nabla E_\theta(\vec{x}^{\,\prime}) \quad (\text{∵ } \sum_{\vec{x}}P_{\text{data}}(\vec{x})=1) \\ &= – \mathbb{E}_{P_{\text{data}}} [\nabla E_\theta(\vec{X})] + \mathbb{E}_{P_\theta} [\nabla E_\theta(\vec{X})] \end{aligned} $$

これで、勾配の形式が得られました。ここからは、\(\nabla E_\theta(\vec{X})\) について、具体的なパラメータで微分計算を行います。

パラメータ更新

バイアス項の微分

まず、バイアス \(b_i\) に関する微分を行います。

$$ \begin{aligned} \frac{\partial E_\theta(\vec{x})}{\partial b_i} &= \frac{\partial}{\partial b_i} \left( -\sum_k b_k x_k – \sum_{k<\ell} w_{k\ell}x_kx_\ell \right) \\ &= -x_i \end{aligned} $$

重み項の微分

同様に、重み \(w_{ij}\) に関する微分を行います。

$$ \begin{aligned} \frac{\partial E_\theta(\vec{x})}{\partial w_{ij}} &= \frac{\partial}{\partial w_{ij}} \left( -\sum_k b_k x_k – \sum_{k<\ell} w_{k\ell}x_kx_\ell \right) \\ &= -x_ix_j \end{aligned} $$

勾配式への代入

これらを勾配式に代入します。

$$ \begin{aligned} \frac{\partial f}{\partial b_i} &= – \mathbb{E}_{\text{data}} \left[ \frac{\partial E_\theta(\vec{x})}{\partial b_i} \right] + \mathbb{E}_\theta \left[ \frac{\partial E_\theta(\vec{x})}{\partial b_i} \right] \\ &= \mathbb{E}_{\text{data}}[x_i] – \mathbb{E}_\theta[x_i] \end{aligned} $$

$$ \begin{aligned} \frac{\partial f}{\partial w_{ij}} &= – \mathbb{E}_{\text{data}} \left[ \frac{\partial E_\theta(\vec{x})}{\partial w_{ij}} \right] + \mathbb{E}_\theta \left[ \frac{\partial E_\theta(\vec{x})}{\partial w_{ij}} \right] \\ &= \mathbb{E}_{\text{data}}[x_ix_j] – \mathbb{E}_\theta[x_ix_j] \end{aligned} $$

繰り返しになりますが、対数尤度 \(f(\theta)\) を最大化するため勾配上昇法を用います。

$$ \theta \leftarrow \theta + \eta \nabla f(\theta) \tag{9} $$

したがって、バイアスと重みの更新式は以下のようになります。

$$ b_i \leftarrow b_i + \eta \left( \mathbb{E}_{\text{data}}[x_i] – \mathbb{E}_\theta[x_i] \right) \tag{10} $$

$$ w_{ij} \leftarrow w_{ij} + \eta \left( \mathbb{E}_{\text{data}}[x_ix_j] – \mathbb{E}_\theta[x_ix_j] \right) \tag{11} $$

これで、更新式まで導出できました。BMの学習では、パラメータはデータの期待値とモデルの期待値が一致するように更新されることが分かります。次の節からは、これらの更新式を実際にどのように計算するのかを説明します。

期待値の計算

更新式において、データ分布に関する期待値 \(\mathbb{E}_{\text{data}}\) は標本平均として容易に計算できます。一方、モデル分布に関する期待値 \(\mathbb{E}_\theta\) を厳密に計算するには、分配関数の計算が必要になります。すると、状態数 \(2^N\) の総和を含むため、現実的には計算が困難になります。

そこで、モデル分布 \(P_\theta(\vec{x})\) からサンプルを生成し、サンプル平均によって期待値を近似的に計算することを考えます。しかし、ここで一つ重要な問題があります。モデル分布

$$ P_\theta(\vec{x}) = \frac{1}{Z_\theta}\exp(-E_\theta(\vec{x})) $$

から直接サンプルしようとしても、分配関数 \(Z_\theta\) を計算できないため、現実的ではありません。つまり、「モデル分布からサンプルしたい」という要請は明確である一方で、「そのままではサンプルできない」という困難があるわけです。

この問題に対して、一般にはマルコフ連鎖モンテカルロ（Markov Chain Monte Carlo; MCMC）法を用います。MCMCでは、目的とする分布を定常分布として持つマルコフ連鎖を構成し、その連鎖を十分に長く回すことで、目的分布からのサンプルを得ます。

では、BMではどのようにマルコフ連鎖を構成すればよいのでしょうか。ここで重要になるのが、同時分布 \(P_\theta(\vec{x})\) は直接扱いにくい一方で、各変数の条件付き確率

$$ P(x_i \mid \vec{x}_{-i}) $$

は分配関数を含まない形で計算できるという点です。この性質を利用すると、他の変数を固定した条件で、ある一つの変数だけを順番に更新することができます。こうして各変数を条件付き確率に従って順番に更新していくサンプリング法が、ギブスサンプリング（Gibbs sampling） です。

具体的には、現在の状態 \(\vec{x}\) に対して、

$$ x_1 \sim P(x_1 \mid x_2,\dots,x_N) $$

$$ x_2 \sim P(x_2 \mid x_1,x_3,\dots,x_N) $$

$$ \vdots $$

$$ x_N \sim P(x_N \mid x_1,\dots,x_{N-1}) $$

のように、一変数ずつ更新していきます。この操作を繰り返すことで、状態列

$$ \vec{x}^{(1)},\vec{x}^{(2)},\vec{x}^{(3)},\dots $$

が得られます。理論的には、このマルコフ連鎖の定常分布は \(P_\theta(\vec{x})\) になります。したがって、十分に長く連鎖を回した後のサンプルを用いることで、

$$ \mathbb{E}_{P_\theta}[x_i], \qquad \mathbb{E}_{P_\theta}[x_ix_j] $$

のようなモデル期待値をサンプル平均によって近似できるようになります。以下では、そのために必要となる条件付き確率を具体的に導出します。

まず、条件付き確率の定義から始めます。

$$ \begin{aligned} P(X_i=1\mid \vec{x}_{-i}) &= \frac{P(x_i=1,\vec{x}_{-i})}{P(\vec{x}_{-i})} \\ &= \frac{P(x_i=1,\vec{x}_{-i})} {P(x_i=0,\vec{x}_{-i})+P(x_i=1,\vec{x}_{-i})} \quad (\text{∵ } x_i \in \{0,1\}) \end{aligned} $$

ここで、\(\vec{x}_{-i}\) は成分 \(x_i\) を除いた残りの成分全体を表す記法です。例えば、

$$ \vec{x} = (x_1,\dots,x_{i-1},x_i,x_{i+1},\dots,x_N) $$

とすると、\(\vec{x}_{-i}\) は次で与えられます。

$$ \vec{x}_{-i} = (x_1,\dots,x_{i-1},x_{i+1},\dots,x_N) $$

したがって、\((\vec{x}_{-i},1)\) は「\(i\) 番目の成分だけを 1 に固定し、他の成分は \(\vec{x}_{-i}\) のままにした状態」を表します。同様に、\((\vec{x}_{-i},0)\) は「\(i\) 番目の成分だけを 0 に固定した状態」を表します。

ここで、ボルツマン分布を代入すると、以下のように分配関数を消去できます。

$$ \begin{aligned} P(X_i=1\mid \vec{x}_{-i}) &= \frac{ \frac{1}{Z_\theta}\exp(-E_\theta(\vec{x}_{-i},1)) }{ \frac{1}{Z_\theta}\exp(-E_\theta(\vec{x}_{-i},0)) + \frac{1}{Z_\theta}\exp(-E_\theta(\vec{x}_{-i},1)) } \\ &= \frac{\exp(-E_\theta(\vec{x}_{-i},1))} {\exp(-E_\theta(\vec{x}_{-i},0))+\exp(-E_\theta(\vec{x}_{-i},1))} \quad (\text{∵ } \tfrac{1}{Z_\theta} \text{ が分子分母で共通}) \\ &= \frac{1} {\exp\left(-E_\theta(\vec{x}_{-i},0)+E_\theta(\vec{x}_{-i},1)\right)+1} \\ &= \frac{1} {1+\exp\left(E_\theta(\vec{x}_{-i},1)-E_\theta(\vec{x}_{-i},0)\right)} \end{aligned} $$

ここで、分母のエネルギー差を計算します。エネルギー関数のうち \(x_i\) に依存する部分のみを取り出すと、

$$ \begin{aligned} E_\theta(\vec{x}) &= – b_i x_i – \sum_{j\neq i} w_{ij}x_ix_j + C(\vec{x}_{-i}) \\ &= – x_i \left( b_i+\sum_{j\neq i} w_{ij}x_j \right) + C(\vec{x}_{-i}) \end{aligned} $$

と書けます。ここで \(C(\vec{x}_{-i})\) は \(x_i\) に依存しない項です。

$$ \begin{aligned} E_\theta(\vec{x}_{-i},1) &= – \left( b_i+\sum_{j\neq i} w_{ij}x_j \right) + C(\vec{x}_{-i}) \\ E_\theta(\vec{x}_{-i},0) &= C(\vec{x}_{-i}) \end{aligned} $$

したがって、

$$ E_\theta(\vec{x}_{-i},1)-E_\theta(\vec{x}_{-i},0) = – \left( b_i+\sum_{j\neq i} w_{ij}x_j \right) $$

となります。これを条件付き確率に代入すると、次式を得ます。

$$ \begin{aligned} P(X_i=1\mid \vec{x}_{-i}) &= \frac{1} {1+\exp\left( -\left( b_i+\sum_{j\neq i} w_{ij}x_j \right) \right)} \\ &= \sigma \left( b_i+\sum_{j\neq i} w_{ij}x_j \right) \quad (\text{∵ } \sigma(z)=\tfrac{1}{1+e^{-z}}) \end{aligned} $$

これで、条件付き確率が得られました。以降は、この条件付き確率を用いて各変数を順番に更新し、モデル期待値を近似的に計算していきます。

具体的な計算方法（アルゴリズム）

ここまでで、BM の学習に必要な更新式

$$ b_i \leftarrow b_i + \eta \left( \mathbb{E}_{\text{data}}[x_i] – \mathbb{E}_\theta[x_i] \right), \qquad w_{ij} \leftarrow w_{ij} + \eta \left( \mathbb{E}_{\text{data}}[x_ix_j] – \mathbb{E}_\theta[x_ix_j] \right) $$

が得られました。次に問題となるのは、これらの期待値を実際にどのように計算するかです。この節では、その具体的な計算方法、すなわち学習アルゴリズムを説明します。

データ期待値の計算

まず、データ分布に関する期待値は標本平均として計算できます。学習データが

$$ \vec{x}^{(1)}, \vec{x}^{(2)}, \dots, \vec{x}^{(M)} $$

のように与えられているとします。このとき、データ期待値は次のように近似できます。

$$ \mathbb{E}_{\text{data}}[x_i] \approx \frac{1}{M} \sum_{m=1}^{M} x_i^{(m)} $$

$$ \mathbb{E}_{\text{data}}[x_ix_j] \approx \frac{1}{M} \sum_{m=1}^{M} x_i^{(m)}x_j^{(m)} $$

すなわち、データ期待値は「各データに対して値を計算し、その平均を取る」ことで求まります。これは単なる平均計算なので、比較的容易に実行できます。

モデル期待値の計算

一方で、モデル期待値の計算は難しくなります。例えば、

$$ \mathbb{E}_\theta[x_i] = \sum_{\vec{x}} P_\theta(\vec{x})x_i, \qquad \mathbb{E}_\theta[x_ix_j] = \sum_{\vec{x}} P_\theta(\vec{x})x_ix_j $$

ですが、ここで和を取るべき状態数は \(2^N\) 個あります。そのため、厳密な計算は事実上不可能です。

そこで、前節で導いたギブスサンプリングを用いてモデル分布からサンプル列を生成し、その平均によってモデル期待値を近似します。サンプル列を

$$ \vec{x}^{[1]}, \vec{x}^{[2]}, \dots, \vec{x}^{[K]} $$

とすると、モデル期待値は次のように近似されます。

$$ \mathbb{E}_\theta[x_i] \approx \frac{1}{K} \sum_{k=1}^{K} x_i^{[k]} $$

$$ \mathbb{E}_\theta[x_ix_j] \approx \frac{1}{K} \sum_{k=1}^{K} x_i^{[k]}x_j^{[k]} $$

ここで重要なのは、モデル期待値は「モデル分布から得られたサンプル列の平均」として計算する、という点です。

ギブスサンプリングの具体的な手順

ギブスサンプリングでは、現在の状態ベクトル

$$ \vec{x}=(x_1,\dots,x_N) $$

を用意し、各成分を 1 つずつ更新していきます。具体的には、ユニット \(i\) について

$$ P(X_i=1 \mid \vec{x}_{-i}) = \sigma \left( b_i+\sum_{j\neq i}w_{ij}x_j \right) $$

を計算し、この確率に従って \(x_i\) を 0 または 1 に更新します。

例えば、確率が

$$ P(X_i=1 \mid \vec{x}_{-i}) = 0.7 $$

であれば、70% の確率で \(x_i=1\)、30% の確率で \(x_i=0\) を選びます。

この操作を \(i=1,\dots,N\) について順番に行います。1 回すべてのユニットを更新する操作を 1 sweep と呼びます。これを繰り返すことで、

$$ \vec{x}^{[1]}, \vec{x}^{[2]}, \vec{x}^{[3]}, \dots $$

というサンプル列が得られます。

BM の学習アルゴリズム

以上をまとめると、BM の学習アルゴリズムは次のようになります。

1. 学習データからデータ期待値 \(\mathbb{E}_{\text{data}}[x_i]\), \(\mathbb{E}_{\text{data}}[x_ix_j]\) を標本平均として計算します。
2. 現在のパラメータ \(\theta\) のもとでギブスサンプリングを実行し、サンプル列 \(\vec{x}^{[1]},\dots,\vec{x}^{[K]}\) を生成します。
3. 生成されたサンプル列からモデル期待値 \(\mathbb{E}_\theta[x_i]\), \(\mathbb{E}_\theta[x_ix_j]\) をサンプル平均として計算します。
4. それらを更新式に代入し、 \(b_i\), \(w_{ij}\) を更新します。

計算量上の問題

このアルゴリズムは理論的には正しいのですが、実際には非常に時間がかかります。最大の理由は、ギブスサンプリングでは 1 変数ずつしか更新できない からです。

すなわち、状態ベクトルに \(N\) 個の変数がある場合、1 sweep だけでも \(N\) 回の条件付き確率計算とサンプリングが必要になります。さらに、モデル分布に十分近づくまでには多くの sweep を繰り返す必要があります。

したがって、BM では

• 1 回のサンプリングに時間がかかる
• 多くのサンプルを集める必要がある
• 結果として学習全体が非常に遅くなる

という問題が生じます。

この「モデル期待値を計算するためのサンプリングが遅い」という点が、BM を実用上扱いにくくしている本質的な理由です。次の節では、この問題を構造的な工夫によって改善したモデルである RBM を説明します。

制限ボルツマンマシン

前節では、ボルツマンマシン（BM）の学習則に現れるモデル期待値 \(\mathbb{E}_\theta[x_i x_j]\) が、状態数 \(2^N\) にわたる総和を含むため、直接計算が困難であることを見ました。

制限ボルツマンマシン（Restricted Boltzmann Machine; RBM）は、この計算を扱いやすくすることを目的として導入されるモデルです。RBM ではユニットを可視変数 \(\vec{v}=(v_1,\dots,v_D)\) と隠れ変数 \(\vec{h}=(h_1,\dots,h_F)\) の 2 層に分け、層内結合を持たない二部グラフ構造を仮定します。

具体的な構造制約は以下の通りです。

• 可視変数同士の結合を持ちません。
• 隠れ変数同士の結合を持ちません。
• 可視変数と隠れ変数の間にのみ結合を持ちます。

BM と同様に、RBM でもエネルギー関数を導入し、そのエネルギーに基づいて確率分布を定義し、対数尤度を最大化するように学習します。一方で、RBM ではこの構造制約によって条件付き独立性が成立するため、BM に比べてサンプリングや期待値計算が大幅に簡単になります。以後、RBM のパラメータは \(\theta=(\vec{a},\vec{b},W)\) とまとめて表します。

RBMのエネルギー関数と確率分布

RBM のエネルギー関数は以下で与えます。

$$ E_\theta(\vec{v},\vec{h}) = -\sum_{i=1}^{D} a_i v_i -\sum_{j=1}^{F} b_j h_j -\sum_{i=1}^{D}\sum_{j=1}^{F} v_i W_{ij} h_j \tag{13} $$

ここで、\(v_i \in \{0,1\}\) は可視ユニット、\(h_j \in \{0,1\}\) は隠れユニット、\(a_i\) は可視ユニットのバイアス、\(b_j\) は隠れユニットのバイアス、\(W_{ij}\) は可視ユニット \(i\) と隠れユニット \(j\) の結合重みを表します。

BM と同様に、エネルギー関数に基づいて同時分布を導入します。定義は以下の通りです。

$$ P_\theta(\vec{v},\vec{h}) = \frac{1}{Z_\theta}\exp(-E_\theta(\vec{v},\vec{h})) \tag{14} $$

ここで分配関数 \(Z_\theta\) は以下で定義します。

$$ Z_\theta = \sum_{\vec{v}} \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) \tag{15} $$

ここまでは BM と同様であり、「エネルギー関数を定め、その指数関数から確率分布を構成する」という考え方は共通です。一方で RBM では、学習データとして直接観測されるのは可視変数 \(\vec{v}\) のみであり、隠れ変数 \(\vec{h}\) は観測されません。そのため、可視変数の周辺分布 \(P_\theta(\vec{v})\) を直接扱う必要があります。

可視変数の周辺分布は以下で定義されます。

$$ P_\theta(\vec{v}) = \sum_{\vec{h}} P_\theta(\vec{v},\vec{h}) \tag{16} $$

式(14) を代入すると、次式を得ます。

$$ \begin{aligned} P_\theta(\vec{v}) &= \frac{1}{Z_\theta} \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) \quad  \end{aligned} \tag{17} $$

BM では各状態 \(\vec{x}\) の確率をそのまま扱えば十分でしたが、RBM では隠れ変数を周辺化しなければならない点が最初の大きな違いです。

学習目標

RBM でも、モデル分布をデータ分布に近づけることを目的とします。ただし、ここで直接扱うのは可視変数の分布です。この点は BM における「観測される変数の分布をモデル化する」という考え方と同じですが、RBM では観測変数が可視層だけに限定される点が異なります。

したがって、対数尤度の期待値を以下のように導入します。

$$ f(\theta) = \mathbb{E}_{P_{\text{data}}} [\log P_\theta(\vec{v})] $$

ここで、式(17)を代入すると、

$$ \begin{aligned} \log P_\theta(\vec{v}) &= \log \left( \frac{1}{Z_\theta} \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) \right) \\ &= \log \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) – \log Z_\theta \quad (\text{∵ } \log \tfrac{a}{b}=\log a-\log b) \end{aligned} \tag{18} $$

となります。BM では \(\log P_\theta(\vec{x})\) がエネルギー項と分配関数項だけで直接表されましたが、RBM では隠れ変数 \(\vec{h}\) を周辺化する必要があるため、式(18) のように \(\log\sum_{\vec{h}}\exp(-E_\theta(\vec{v},\vec{h}))\) が現れます。

更新式の導出

BM と同様に、学習目標で導入した目的関数 \(f(\theta)=\mathbb{E}_{P_{\text{data}}}[\log P_\theta(\vec{v})]\) を最大化するため、勾配上昇法を用います。したがって、以下の勾配が計算できれば十分です。

$$ \nabla f(\theta) = \sum_{\vec{v}} P_{\text{data}}(\vec{v}) \nabla\log P_\theta(\vec{v}) \tag{19} $$

まず、\(\nabla\log P_\theta(\vec{v})\) を計算します。

$$ \nabla \log P_\theta(\vec{v}) = \nabla \log \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) – \nabla \log Z_\theta $$

ここでも BM と同様に、「対数確率の勾配を二つの項に分けて計算する」という流れは同じです。ただし、RBM では第1項に隠れ変数に関する総和が入るため、その扱いが少し変わります。

第一項の計算

まず、

$$ \log \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) $$

を微分します。

$$ \begin{aligned} \nabla \log \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) &= \frac{ \nabla \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) }{ \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) } \quad (\text{∵ } \nabla \log g = \tfrac{\nabla g}{g}) \\ &= \frac{ \sum_{\vec{h}} \nabla \exp(-E_\theta(\vec{v},\vec{h})) }{ \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) } \\ &= \frac{ \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) (-\nabla E_\theta(\vec{v},\vec{h})) }{ \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) } \quad (\text{∵ } \nabla e^f=e^f\nabla f) \\ &= – \sum_{\vec{h}} \frac{ \exp(-E_\theta(\vec{v},\vec{h})) }{ \sum_{\vec{h}^{\,\prime}} \exp(-E_\theta(\vec{v},\vec{h}^{\,\prime})) } \nabla E_\theta(\vec{v},\vec{h}) \end{aligned} $$

ここで、条件付き分布

$$ P_\theta(\vec{h}\mid \vec{v}) = \frac{\exp(-E_\theta(\vec{v},\vec{h}))} {\sum_{\vec{h}^{\,\prime}} \exp(-E_\theta(\vec{v},\vec{h}^{\,\prime}))} $$

を導入します。

これを用いると、

$$ \begin{aligned} \nabla \log \sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) &= – \sum_{\vec{h}} P_\theta(\vec{h}\mid \vec{v}) \nabla E_\theta(\vec{v},\vec{h}) \\ &= – \mathbb{E}_{P_\theta(\vec{H}\mid \vec{V}=\vec{v})} [\nabla E_\theta(\vec{v},\vec{h})] \end{aligned} $$

となります。BM では第1項は単純に \(-\nabla E_\theta(\vec{x})\) でしたが、RBM では隠れ変数を周辺化しているため、条件付き分布 \(P_\theta(\vec{h}\mid\vec{v})\) による期待値に変わります。

第二項の計算

次に、\(\nabla \log Z_\theta\) を計算します。

分配関数は

$$ Z_\theta = \sum_{\vec{v}}\sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) $$

であるため、

$$ \begin{aligned} \nabla Z_\theta &= \sum_{\vec{v}}\sum_{\vec{h}} \nabla \exp(-E_\theta(\vec{v},\vec{h})) \\ &= \sum_{\vec{v}}\sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) (-\nabla E_\theta(\vec{v},\vec{h})) \quad (\text{∵ } \nabla e^f=e^f\nabla f) \end{aligned} $$

したがって、

$$ \begin{aligned} \nabla \log Z_\theta &= \frac{1}{Z_\theta}\nabla Z_\theta \\ &= \frac{1}{Z_\theta} \sum_{\vec{v}}\sum_{\vec{h}} \exp(-E_\theta(\vec{v},\vec{h})) (-\nabla E_\theta(\vec{v},\vec{h})) \\ &= – \sum_{\vec{v}}\sum_{\vec{h}} \frac{\exp(-E_\theta(\vec{v},\vec{h}))}{Z_\theta} \nabla E_\theta(\vec{v},\vec{h}) \\ &= – \sum_{\vec{v}}\sum_{\vec{h}} P_\theta(\vec{v},\vec{h}) \nabla E_\theta(\vec{v},\vec{h}) \quad (\text{∵ } P_\theta(\vec{v},\vec{h})=\tfrac{\exp(-E_\theta(\vec{v},\vec{h}))}{Z_\theta}) \\ &= – \mathbb{E}_{P_\theta(\vec{V},\vec{H})} [\nabla E_\theta(\vec{v},\vec{h})] \end{aligned} $$

この部分は、BM における分配関数の微分と全く同じ考え方です。違いは、状態空間が \(\vec{x}\) ではなく \((\vec{v},\vec{h})\) に広がっていることだけです。

勾配の整理

以上をまとめると、

$$ \nabla \log P_\theta(\vec{v}) = – \mathbb{E}_{P_\theta(\vec{H}\mid \vec{V}=\vec{v})} [\nabla E_\theta(\vec{v},\vec{h})] + \mathbb{E}_{P_\theta(\vec{V},\vec{H})} [\nabla E_\theta(\vec{v},\vec{h})] $$

となります。

期待値による表現

以上の結果を式(19) に代入すると、目的関数の勾配は次のように期待値で表せます。

$$ \nabla f(\theta) = – \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})} [\nabla E_\theta(\vec{v},\vec{h})] + \mathbb{E}_{P_\theta(\vec{V},\vec{H})} [\nabla E_\theta(\vec{v},\vec{h})] \tag{20} $$

この式は、データ分布におけるエネルギー勾配と、モデル分布におけるエネルギー勾配の差として勾配が表されることを意味します。

式(20) は、BM における勾配式

$$ – \mathbb{E}_{P_{\text{data}}}[\nabla E_\theta(\vec{X})] + \mathbb{E}_{P_\theta}[\nabla E_\theta(\vec{X})] $$

に対応しています。

パラメータ更新

次に、エネルギー関数を各パラメータで微分し、式(20) に代入して更新式を得ます。

$$ E_\theta(\vec{v},\vec{h}) = -\sum_i a_i v_i -\sum_j b_j h_j -\sum_{i,j} v_iW_{ij}h_j $$

可視バイアスの微分

$$ \begin{aligned} \frac{\partial E_\theta(\vec{v},\vec{h})}{\partial a_i} &= \frac{\partial}{\partial a_i} \left( -\sum_k a_k v_k -\sum_j b_j h_j -\sum_{k,j} v_k W_{kj} h_j \right) \\ &= -v_i \end{aligned} $$

隠れバイアスの微分

$$ \begin{aligned} \frac{\partial E_\theta(\vec{v},\vec{h})}{\partial b_j} &= \frac{\partial}{\partial b_j} \left( -\sum_i a_i v_i -\sum_\ell b_\ell h_\ell -\sum_{i,\ell} v_i W_{i\ell} h_\ell \right) \\  &= -h_j \end{aligned} $$

重みの微分

$$ \begin{aligned} \frac{\partial E_\theta(\vec{v},\vec{h})}{\partial W_{ij}} &= \frac{\partial}{\partial W_{ij}} \left( -\sum_k a_k v_k -\sum_\ell b_\ell h_\ell -\sum_{k,\ell} v_k W_{k\ell} h_\ell \right) \\  &= -v_i h_j \end{aligned} $$

勾配式への代入

重みに関する勾配は

$$ \begin{aligned} \frac{\partial f}{\partial W_{ij}} &= – \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})} \left[ \frac{\partial E_\theta(\vec{v},\vec{h})}{\partial W_{ij}} \right] + \mathbb{E}_{P_\theta(\vec{V},\vec{H})} \left[ \frac{\partial E_\theta(\vec{v},\vec{h})}{\partial W_{ij}} \right] \\ &= – \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})} [-v_i h_j] + \mathbb{E}_{P_\theta(\vec{V},\vec{H})} [-v_i h_j] \quad (\text{∵ } \tfrac{\partial E_\theta(\vec{v},\vec{h})}{\partial W_{ij}}=-v_i h_j) \\ &= \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})} [v_i h_j] – \mathbb{E}_{P_\theta(\vec{V},\vec{H})} [v_i h_j] \end{aligned} $$

同様に、可視バイアスについては

$$ \frac{\partial f}{\partial a_i} = \mathbb{E}_{P_{\text{data}}}[v_i] – \mathbb{E}_{P_\theta(\vec{V},\vec{H})}[v_i] $$

隠れバイアスについては

$$ \frac{\partial f}{\partial b_j} = \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})}[h_j] – \mathbb{E}_{P_\theta(\vec{V},\vec{H})}[h_j] $$

となります。

更新式

ここでも対数尤度を最大化するため、BM と同様に勾配上昇法を用います。

$$ \theta \leftarrow \theta + \eta \nabla f(\theta) $$

したがって、重みの更新式は以下で与えられます。

$$ W_{ij} \leftarrow W_{ij} + \eta \left( \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})}[v_i h_j] – \mathbb{E}_{P_\theta(\vec{V},\vec{H})}[v_i h_j] \right) $$

同様に、可視バイアスと隠れバイアスの更新は以下です。

$$ a_i \leftarrow a_i + \eta \left( \mathbb{E}_{P_{\text{data}}}[v_i] – \mathbb{E}_{P_\theta(\vec{V},\vec{H})}[v_i] \right) $$

$$ b_j \leftarrow b_j + \eta \left( \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})}[h_j] – \mathbb{E}_{P_\theta(\vec{V},\vec{H})}[h_j] \right) $$

BM において更新式が

$$ \mathbb{E}_{\text{data}}[x_i x_j] – \mathbb{E}_\theta[x_i x_j] $$

の形をとったのに対し、RBM では

$$ \mathbb{E}_{\text{data}}[v_i h_j] – \mathbb{E}_\theta[v_i h_j] $$

の形に置き換わります。すなわち、BM におけるユニット間相関が、RBM では可視・隠れ間相関として表れます。式の形は非常によく似ていますが、RBM では観測変数と潜在変数の相関を学習している点が特徴です。

条件付き独立性

RBM が BM よりも扱いやすい最大の理由は、構造制約により条件付き独立性が成立することです。BM では一般に全てのユニットが相互作用するため、ギブスサンプリングでも 1 ユニットずつ順番に更新する必要がありました。一方で RBM では、可視層と隠れ層の二部構造によって、同一層内のユニットが条件付き独立になります。

隠れ変数の条件付き分布

可視変数 \(\vec v\) を固定したときの隠れ変数の条件付き分布を考えます。条件付き分布の定義より、

$$
P_\theta(\vec h \mid \vec v)
= \frac{P_\theta(\vec v,\vec h)}{P_\theta(\vec v)}
\propto P_\theta(\vec v,\vec h)
\propto \exp\bigl(-E_\theta(\vec v,\vec h)\bigr)
$$

です。ここで、RBM のエネルギー関数
$$
E_\theta(\vec v,\vec h)
= -\sum_i a_i v_i – \sum_j b_j h_j – \sum_{i,j} W_{ij} v_i h_j
$$
を、\(\vec v\) にのみ依存する項と \(\vec h\) を含む項に分けると、

$$
E_\theta(\vec v,\vec h)
=
\underbrace{-\sum_i a_i v_i}_{=:C(\vec v)}
-\sum_j h_j \left( b_j + \sum_i W_{ij} v_i \right)
$$

と書けるので、

$$
\begin{aligned}
P_\theta(\vec h \mid \vec v)
&\propto \exp\bigl(-E_\theta(\vec v,\vec h)\bigr) \\
&= \exp\left(-C(\vec v) + \sum_j h_j\left(b_j+\sum_i W_{ij}v_i\right)\right) \\
&= \exp\bigl(-C(\vec v)\bigr)
\prod_j \exp\left(h_j\left(b_j+\sum_i W_{ij}v_i\right)\right) \\
&\propto \prod_j \exp\left(h_j\left(b_j+\sum_i W_{ij}v_i\right)\right).
\end{aligned}
$$

したがって、正規化をすると以下のように表せます。

$$
\begin{aligned}
P_\theta(\vec h \mid \vec v)
&=
\frac{
\prod_j \exp\left(h_j\left(b_j+\sum_i W_{ij}v_i\right)\right)
}{
\sum_{\vec h}
\prod_j \exp\left(h_j\left(b_j+\sum_i W_{ij}v_i\right)\right)
} \\
&=
\frac{
\prod_j \exp\left(h_j\left(b_j+\sum_i W_{ij}v_i\right)\right)
}{
\prod_j \sum_{h_j\in\{0,1\}}
\exp\left(h_j\left(b_j+\sum_i W_{ij}v_i\right)\right)
} \\
&=
\prod_j
\frac{
\exp\left(h_j\left(b_j+\sum_i W_{ij}v_i\right)\right)
}{
\sum_{h_j\in\{0,1\}}
\exp\left(h_j\left(b_j+\sum_i W_{ij}v_i\right)\right)
}.
\end{aligned}
$$

よって、

$$
P_\theta(\vec h \mid \vec v)=\prod_j P_\theta(h_j\mid \vec v) \tag{21}
$$

が成り立ちます。すなわち、\(\vec v\) を固定すると、各隠れユニットは条件付き独立になります。これは、RBM では隠れユニット同士の相互作用項が存在しないためです。

さらに、各 \(j\) について
$$
\lambda_j := b_j+\sum_i W_{ij}v_i
$$
とおくと、上式より

$$
P_\theta(h_j\mid \vec v)
=
\frac{\exp(h_j\lambda_j)}{\sum_{h_j\in\{0,1\}}\exp(h_j\lambda_j)}
=
\frac{\exp(h_j\lambda_j)}{1+\exp(\lambda_j)}
$$

となります。したがって、\(h_j=1\) の場合には

$$
\begin{aligned}
P_\theta(H_j=1\mid \vec v)
&=
\frac{\exp(\lambda_j)}{1+\exp(\lambda_j)} \\
&=
\frac{1}{1+\exp(-\lambda_j)} \\
&=
\sigma(\lambda_j).
\end{aligned}
$$

ゆえに、

$$
P_\theta(H_j=1\mid \vec v)
=
\sigma\left(b_j+\sum_i W_{ij}v_i\right) \tag{22}
$$

を得ます。

可視変数の条件付き分布

同様に、隠れ変数 \(\vec{h}\) を固定したとき、可視変数の条件付き分布も分解されます。したがって、次が成り立ちます。

$$ P_\theta(\vec{v}\mid \vec{h}) = \prod_i P_\theta(v_i\mid \vec{h}) \tag{23} $$

各可視ユニットの条件付き確率は次式で与えられます。

$$ P_\theta(V_i=1\mid \vec{h}) = \sigma \left( a_i+\sum_j W_{ij}h_j \right) \tag{24} $$

BMとの比較

BM では、一般に \(P(X_i\mid \vec{X}_{-i})\) を 1 ユニットずつ更新する必要がありました。一方 RBM では、式(21) および式(23) が成立するため、同一層内のユニットを並列に更新できます。

この違いは計算量に直結します。BM では 1 ユニットずつ順番にサンプリングするため時間がかかりましたが、RBM では可視層または隠れ層をまとめて更新できるため、ギブスサンプリングが大幅に高速化されます。ここが、RBM が BM よりも実用的である最大の理由です。

RBMにおける具体的な計算方法（アルゴリズム）

ここまでで、RBM の学習に必要な更新式が得られました。次に、BM の場合と同様に、これらの期待値を実際にどのように計算するのかを説明します。

データ期待値の計算

データ側の期待値は、BM の場合と同様に標本平均として計算します。学習データが

$$ \vec{v}^{(1)}, \vec{v}^{(2)}, \dots, \vec{v}^{(M)} $$

と与えられているとすると、例えば可視バイアスのデータ期待値は

$$ \mathbb{E}_{\text{data}}[v_i] \approx \frac{1}{M} \sum_{m=1}^{M} v_i^{(m)} $$

のように計算できます。

一方、重みの更新に現れる

$$ \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})}[v_i h_j] $$

は、データ \(\vec{v}^{(m)}\) を与えたときの隠れユニットの条件付き期待値を使って計算できます。すなわち、

$$ \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})}[v_i h_j] \approx \frac{1}{M} \sum_{m=1}^{M} v_i^{(m)} P_\theta(H_j=1\mid \vec{v}^{(m)}) $$

となります。ここが BM と異なる点であり、RBM では隠れ変数が観測されないため、データから直接相関を数えるのではなく、条件付き確率を介して補います。

モデル期待値の計算

一方、モデル期待値

$$ \mathbb{E}_{P_\theta(\vec{V},\vec{H})}[v_i h_j] $$

は、BM と同様に厳密計算が困難です。そのため、モデル分布からサンプルを生成し、サンプル平均によって近似します。

ここで RBM の利点が現れます。BM では 1 ユニットずつしか更新できませんでしたが、RBM では式(21) と式(23) により、隠れ層は可視層を固定すると並列更新でき、可視層は隠れ層を固定すると並列更新できます。したがって、ギブスサンプリングは

1. 可視ベクトル \(\vec{v}\) を固定し、各隠れユニットを並列にサンプリングする
2. 得られた隠れベクトル \(\vec{h}\) を固定し、各可視ユニットを並列にサンプリングする

という 2 段階で実行できます。

このようにして得られたサンプル列

$$ (\vec{v}^{[1]},\vec{h}^{[1]}), (\vec{v}^{[2]},\vec{h}^{[2]}), \dots, (\vec{v}^{[K]},\vec{h}^{[K]}) $$

を用いると、モデル期待値は

$$ \mathbb{E}_{P_\theta(\vec{V},\vec{H})}[v_i h_j] \approx \frac{1}{K} \sum_{k=1}^{K} v_i^{[k]} h_j^{[k]} $$

として近似できます。

RBM の学習アルゴリズム

以上をまとめると、RBM の学習アルゴリズムは次のようになります。

1. 学習データから、可視ユニットの平均 \(\mathbb{E}_{\text{data}}[v_i]\) を計算します。
2. 各データ \(\vec{v}^{(m)}\) に対して \(P_\theta(H_j=1\mid \vec{v}^{(m)})\) を計算し、データ側の相関 \(\mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})}[v_i h_j]\) を求めます。
3. 現在のパラメータ \(\theta\) のもとでギブスサンプリングを実行し、モデルサンプル \((\vec{v}^{[k]},\vec{h}^{[k]})\) を生成します。
4. そのサンプル列から、モデル期待値 \(\mathbb{E}_{P_\theta(\vec{V},\vec{H})}[v_i h_j]\) をサンプル平均として計算します。
5. 更新式に代入し、 \(W_{ij}\), \(a_i\), \(b_j\) を更新します。

まとめ

RBM では、可視変数と隠れ変数の間だけに結合を許すことで、BM よりも単純な構造を得ます。

この構造制約により、学習の勾配は次の形に整理されます。

$$ \mathbb{E}_{\text{data}}[v_i h_j] – \mathbb{E}_{\text{model}}[v_i h_j] $$

さらに、条件付き独立性により式(22) および式(24) が成り立つため、BM に比べてサンプリングと期待値計算がはるかに扱いやすくなります。特に、BM では 1 ユニットずつしか更新できなかったのに対し、RBM では可視層・隠れ層をまとめて並列更新できる点が大きな利点です。

もっとも、RBM においてもモデル期待値 \(\mathbb{E}_{P_\theta(\vec{V},\vec{H})}[v_i h_j]\) を厳密に求めることは難しく、依然としてサンプリングが必要です。次の節では、このモデル期待値をさらに効率よく近似するために提案されたコントラスティブ・ダイバージェンス（CD法）について説明します。

コントラスティブ・ダイバージェンス（CD法）

前節では、RBM の学習勾配が

$$ \frac{\partial f}{\partial W_{ij}} = \mathbb{E}_{P_{\text{data}}(\vec{v})P_\theta(\vec{h}\mid \vec{v})}[v_i h_j] – \mathbb{E}_{P_\theta(\vec{V},\vec{H})}[v_i h_j] $$

となることを導きました。ここで問題となるのが第二項であり、BM同様、直接計算は不可能です。

そこで、モデル分布からサンプルを生成し、サンプル平均によって期待値を近似することを考えます。このためにギブスサンプリングを用います。

しかし、ここでも新たな問題が生じます。ギブスサンプリングは理論的にはモデル分布に収束しますが、そのためには非常に多くのステップが必要になります。

つまり、

• 各更新で長いマルコフ連鎖を回す必要があります。
• 学習全体が非常に遅くなります。

という問題があります。

この問題を緩和するために提案されたのが コントラスティブ・ダイバージェンス（Contrastive Divergence; CD） です。

コントラスティブ・ダイバージェンスの基本的な考え方

CD 法の基本的なアイデアは次の通りです。

通常はモデル分布に十分近づくまで長いギブスサンプリングを行う必要があります。

これに対して CD 法では、

• サンプリングを データ点から開始します。
• 数ステップのギブスサンプリングだけ行います。
• その結果をモデルサンプルとして用います。

という近似を行います。

このとき、ギブス更新を \(k\) 回行う CD 法を CD-\(k\) と呼びます。実際には \(k=1\) の CD-1 がよく用いられます。

CD-1 のアルゴリズム

CD-1 の学習手順は次の通りです。

1. データ \(\vec{v}^{(0)}\) を入力します。
2. 隠れユニットをサンプリングします。

$$ P(H_j=1\mid \vec{v}^{(0)}) = \sigma \left( b_j+\sum_i W_{ij}v_i^{(0)} \right) $$

3. 得られた \(\vec{h}^{(0)}\) から可視ユニットを再構成します。

$$ P(V_i=1\mid \vec{h}^{(0)}) = \sigma \left( a_i+\sum_j W_{ij}h_j^{(0)} \right) $$

4. 再構成された可視ベクトルを \(\vec{v}^{(1)}\) とします。
5. もう一度隠れユニットをサンプリングします。

$$ P(H_j=1\mid \vec{v}^{(1)}) = \sigma \left( b_j+\sum_i W_{ij}v_i^{(1)} \right) $$

このように、CD-1 では「データから 1 ステップだけ再構成したサンプル」を用いて、モデル期待値を近似します。

本記事の作成者

鹿内怜央