Python
機械学習
ニューラルネットワーク
ボルツマンマシン

はじめに

こんにちは @ta-ka です.
この記事は 数学とコンピュータⅡ Advent Calendar 2017 の25日目の記事です.
pythonで制約ボルツマンマシンを実装しました.
教科書として『深層学習』を使いました.

本記事の構成

  • はじめに
  • 制約ボルツマンマシン
    • ボルツマン分布
    • 条件付き分布
    • 対数尤度関数
    • パラメータ更新
  • CD法
  • 実装
  • 結果
  • おわりに

制約ボルツマンマシン

下図のような構造を持つ無向グラフを考えます.
$v_{i}$ を可視変数,$h_{j}$ を隠れ変数と呼び,$v_{i}, h_{j}$ は $0$ か $1$ の値をとります.

graph.png

ボルツマン分布

ボルツマン分布を下式で定義します.
指数部分の値が大きいほど,$\boldsymbol v, \boldsymbol h$ の生起確率が高くなるような分布です.

\begin{align}
p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta)
&= \cfrac{1}{Z(\boldsymbol \theta)} \exp \bigl(- \Phi(\boldsymbol v, \boldsymbol h, \boldsymbol \theta) \bigr) \tag{1}
\end{align}

エネルギー関数 $\Phi(\boldsymbol v, \boldsymbol h, \boldsymbol \theta)$,正規化項 $Z(\boldsymbol \theta)$ は以下のように定義してします.

\begin{align}
\Phi(\boldsymbol v, \boldsymbol h, \boldsymbol \theta)
&= - \sum_{i} a_{i} v_{i} - \sum_{j} b_{j} h_{j} - \sum_{i} \sum_{j} w_{ij} v_{i} h_{j} \tag{2} \\
Z(\boldsymbol \theta)
&= \sum_{\boldsymbol v} \sum_{\boldsymbol h} \exp \bigl(- \Phi(\boldsymbol v, \boldsymbol h, \boldsymbol \theta) \bigr) \tag{3}
\end{align}

条件付き分布

可視変数 $\boldsymbol v$ が与えられた時の隠れ変数 $\boldsymbol h$ の条件付き分布 $p(\boldsymbol h \mid \boldsymbol v, \boldsymbol \theta)$ を求めます.

\begin{align}
p(\boldsymbol h \mid \boldsymbol v, \boldsymbol \theta)
&= \cfrac{p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta)}{p(\boldsymbol v \mid \boldsymbol \theta)} \\
&= \cfrac{p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta)}{\sum_{\boldsymbol h} p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta)} \tag{4}
\end{align}

分子を以下のように計算します.

\begin{align}
p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta)
&= \cfrac{1}{Z(\boldsymbol \theta)} \exp \bigl(- \Phi(\boldsymbol v, \boldsymbol h, \boldsymbol \theta) \bigr) \\
&= \cfrac{1}{Z(\boldsymbol \theta)} \exp \left( \sum_{i} a_{i} v_{i} + \sum_{j} b_{j} h_{j} + \sum_{i} \sum_{j} w_{ij} v_{i} h_{j} \right) \\
&= \cfrac{1}{Z(\boldsymbol \theta)} \exp \left( \sum_{i} a_{i} v_{i} \right) \exp \left( \sum_{j} h_{j} \Bigl( b_{j} + \sum_{i} w_{ij} v_{i} \Bigr) \right) \\
&= \cfrac{1}{Z(\boldsymbol \theta)} \exp \left( \sum_{i} a_{i} v_{i} \right) \prod_{j} \exp \left(h_{j} \Bigl( b_{j} + \sum_{i} w_{ij} v_{i} \Bigr) \right) \tag{5}
\end{align}

分母を以下のように計算します.

\begin{align}
\sum_{\boldsymbol h} p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta)
&= \sum_{\boldsymbol h} \cfrac{1}{Z(\boldsymbol \theta)} \exp \bigl(- \Phi(\boldsymbol v, \boldsymbol h, \boldsymbol \theta) \bigr) \\
&= \cfrac{1}{Z(\boldsymbol \theta)} \exp \left( \sum_{i} a_{i} v_{i} \right) \sum_{\boldsymbol h} \prod_{j} \exp \left(h_{j} \Bigl( b_{j} + \sum_{i} w_{ij} v_{i} \Bigr) \right) \\
&= \cfrac{1}{Z(\boldsymbol \theta)} \exp \left( \sum_{i} a_{i} v_{i} \right) \prod_{j} \left(1 + \exp \Bigl( b_{j} + \sum_{i} w_{ij} v_{i} \Bigr) \right) \tag{6}
\end{align}

式$(4)$, $(5)$, $(6)$より以下の式が得られます.

\begin{align}
p(\boldsymbol h \mid \boldsymbol v, \boldsymbol \theta)
&= \prod_{j} \cfrac{\exp \Bigl(h_{j} (b_{j} + \sum_{i} w_{ij} v_{i}) \Bigr)}{1 + \exp \Bigl( b_{j} + \sum_{i} w_{ij} v_{i} \Bigr)} \tag{7} \\
p(h_{j} = 1 \mid \boldsymbol v, \boldsymbol \theta)
&= \cfrac{\exp \Bigl(b_{j} + \sum_{i} w_{ij} v_{i} \Bigr)}{1 + \exp \Bigl( b_{j} + \sum_{i} w_{ij} v_{i} \Bigr)} \\
&= {\rm sigmoid} \left(b_{j} + \sum_{i} w_{ij} v_{i} \right) \tag{8}
\end{align}

同様に, 隠れ変数 $\boldsymbol h$ が与えられた時の可視変数 $\boldsymbol v$ の条件付き分布 $p(\boldsymbol v \mid \boldsymbol h, \boldsymbol \theta)$ が得られます.

\begin{align}
p(\boldsymbol v \mid \boldsymbol h, \boldsymbol \theta)
&= \prod_{i} \cfrac{\exp \Bigl(v_{i} (a_{i} + \sum_{j} w_{ij} h_{j}) \Bigr)}{1 + \exp \Bigl( a_{i} + \sum_{j} w_{ij} h_{j} \Bigr)} \tag{9} \\
p(v_{i} = 1 \mid \boldsymbol h, \boldsymbol \theta)
&= \cfrac{\exp \Bigl(a_{i} + \sum_{j} w_{ij} h_{j} \Bigr)}{1 + \exp \Bigl( a_{i} + \sum_{j} w_{ij} h_{j} \Bigr)} \\
&= {\rm sigmoid} \left(a_{i} + \sum_{j} w_{ij} h_{j} \right) \tag{10}
\end{align}

式$(8)$, $(10)$ はパラメータ更新で利用します.

対数尤度関数

対数尤度関数は以下のように計算されます.

\begin{align}
\ln L &= \ln \prod_{n} p(\boldsymbol v_{n} \mid \boldsymbol \theta) \\
&= \sum_{n} \ln p(\boldsymbol v_{n} \mid \boldsymbol \theta) \\
&= \sum_{n} \left[ \sum_{i} a_{i} v_{ni} + \sum_{j} \ln \left(1 + \exp \Bigl( b_{j} + \sum_{i} w_{ij} v_{ni} \Bigr) \right) - \ln Z(\boldsymbol \theta) \right] \tag{11}
\end{align}

パラメータ更新

パラメータの更新式を求めるために
式$(11)$の対数尤度関数を $w_{ij}$, $a_{i}$, $b_{j}$ で偏微分します.

\begin{align}
\cfrac{\partial \ln L}{\partial w_{ij}}
&= \sum_{n} \cfrac{v_{ni} \exp \Bigl( b_{j} + \sum_{i} w_{ij} v_{ni} \Bigr)}{1 + \exp \Bigl( b_{j} + \sum_{i} w_{ij} v_{ni} \Bigr)} - \cfrac{N}{Z(\boldsymbol \theta)} \sum_{\boldsymbol v} \sum_{\boldsymbol h} v_{i} h_{j} \exp \left( \sum_{i} a_{i} v_{i} + \sum_{j} b_{j} h_{j} + \sum_{i} \sum_{j} w_{ij} v_{i} h_{j} \right) \\
&= \sum_{n} v_{ni} {\rm sigmoid} \left(b_{j} + \sum_{i} w_{ij} v_{ni} \right) - N \sum_{\boldsymbol v} \sum_{\boldsymbol h} v_{i} h_{j} p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta) \\
&= \sum_{n} v_{ni} p(h_{nj} = 1 \mid \boldsymbol v_{n}, \boldsymbol \theta) - N \sum_{\boldsymbol v} \sum_{\boldsymbol h} v_{i} h_{j} p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta) \tag{12} \\ \\

\cfrac{\partial \ln L}{\partial a_{i}}
&= \sum_{n} v_{ni} - \cfrac{N}{Z(\boldsymbol \theta)} \sum_{\boldsymbol v} \sum_{\boldsymbol h} v_{i} \exp \left( \sum_{i} a_{i} v_{i} + \sum_{j} b_{j} h_{j} + \sum_{i} \sum_{j} w_{ij} v_{i} h_{j} \right) \\
&= \sum_{n} v_{ni} - N \sum_{\boldsymbol v} \sum_{\boldsymbol h} v_{i} p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta) \tag{13} \\ \\

\cfrac{\partial \ln L}{\partial b_{j}}
&= \sum_{n} \cfrac{v_{ni} \exp \Bigl( b_{j} + \sum_{i} w_{ij} v_{ni} \Bigr)}{1 + \exp \Bigl( b_{j} + \sum_{i} w_{ij} v_{ni} \Bigr)} - \cfrac{N}{Z(\boldsymbol \theta)} \sum_{\boldsymbol v} \sum_{\boldsymbol h} h_{j} \exp \left( \sum_{i} a_{i} v_{i} + \sum_{j} b_{j} h_{j} + \sum_{i} \sum_{j} w_{ij} v_{i} h_{j} \right) \\
&= \sum_{n} v_{ni} {\rm sigmoid} \left(b_{j} + \sum_{i} w_{ij} v_{ni} \right) - N \sum_{\boldsymbol v} \sum_{\boldsymbol h} h_{j} p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta) \\
&= \sum_{n} v_{ni} p(h_{nj} = 1 \mid \boldsymbol v_{n}, \boldsymbol \theta) - N \sum_{\boldsymbol v} \sum_{\boldsymbol h} h_{j} p(\boldsymbol v, \boldsymbol h \mid \boldsymbol \theta) \tag{14}
\end{align}

式$(12)$, $(13)$, $(14)$の第 $2$ 項は $v_{i}$, $h_{j}$ のすべての状態における重み付き和であり,計算コストが大きくなります.
例えば可視変数,隠れ変数が合わせて $100$ 個ある場合,状態の通り数は $2^{100}$ 通りになります.

CD法

上述のように式$(12)$, $(13)$, $(14)$の第 $2$ 項は計算コストが大きいため,
ギブスサンプリングを使って近似的に計算することを考えます.
CD法では,初期値 $\boldsymbol v^{(0)} = \boldsymbol v_{n}$ とし,サンプリングを $T$ 回繰り返します.
サンプリングした結果を使ってパラメータを更新します.

\begin{align}
\Delta w_{ij} &= \epsilon \left(v_{i}^{(0)} p_{j}^{(0)} - v_{i}^{(T)} p_{j}^{(T)} \right) \tag{15} \\
\Delta a_{i} &= \epsilon \left(v_{i}^{(0)} - v_{i}^{(T)} \right) \tag{16} \\
\Delta b_{j} &= \epsilon \left(p_{j}^{(0)} - p_{j}^{(T)} \right) \tag{17}
\end{align}

ただし,$p_{j}^{(T)} = p(h_{j} = 1 \mid \boldsymbol v^{(T)}, \boldsymbol \theta)$ です.
(CD法のサンプリングについては理解することができなかったので,後日勉強して追記します.)

実装

RBMによる二値画像の自己符号化器を実装しました.
python3で実装したプログラムを掲載します.

import numpy as np

class RBM:

    def __init__(self, n_v, n_h):
        self.w = np.random.randn(n_h, n_v)
        self.b = np.random.randn(n_h, 1)
        self.a = np.random.randn(n_v, 1)

    def train(self, V, epsilon, epoch, T):
        for epo in range(epoch):
            for (n, v_0) in enumerate(V.T):
                v_0 = np.copy(v_0).reshape(-1, 1)
                p_h_0 = np.copy(self.sigmoid(self.w.dot(v_0) + self.b))

                v, p_h = self.encode_decode(np.copy(v_0), T)
                self.update(v_0, v, p_h_0, p_h, epsilon)

    def encode_decode(self, v, T):
        for t in range(T):
            # visible
            p_h = self.sigmoid(self.w.dot(v) + self.b)
            h = (np.random.rand(n_h, 1) < p_h).astype('float64')

            # hidden
            p_v = self.sigmoid(self.w.T.dot(h) + self.a)
            v = (np.random.rand(n_v, 1) < p_v).astype('float64')

        return (v, p_h)

    def update(self, v_0, v, p_h_0, p_h, epsilon):
        self.w += epsilon * (v_0.T * p_h_0 - v.T * p_h)
        self.a += epsilon * (v_0 - v)
        self.b += epsilon * (p_h_0 - p_h)

    def sigmoid(self, x):
        return 1.0 / (1.0 + np.exp(-x))


if __name__ == "__main__":
    # init parameters
    N = 20
    side = 8
    n_v = side ** 2
    n_h = 32
    V = np.round(np.random.rand(n_v, N))
    T = 5

    # train
    rbm = RBM(n_v, n_h)
    rbm.train(V, 0.1, 200, T)

結果

$20$ サンプルで学習し,復元した画像のうち $10$ サンプルを下図に示します.
左がランダムに生成した $8 \times 8$ の元画像です.
右がRBMにより復元した画像です.隠れ変数の数を $32$ に設定し,$64$ 次元の元画像を復元しています.

rbm.png

おわりに

RBMによる二値画像の自己符号化器を実装できました.
CD法のサンプリングを勉強して理解できたら追記します.