pythonで制約ボルツマンマシン実装


# はじめに

こんにちは @ta-ka です．
この記事は [数学とコンピュータⅡ Advent Calendar 2017](https://qiita.com/advent-calendar/2017/math-and-computer2) の25日目の記事です．
pythonで制約ボルツマンマシンを実装しました．
教科書として[『深層学習』](http://amzn.asia/0hEJtVx)を使いました．

本記事の構成

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


# 制約ボルツマンマシン

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

<img width="300" alt="graph.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/08e34d3d-cfc0-4342-f592-385716a64d85.png">


## ボルツマン分布

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

```math
\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)$ は以下のように定義してします．

```math
\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)$ を求めます．

```math
\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}
```

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

```math
\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}
```

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

```math
\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)$より以下の式が得られます．

```math
\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)$ が得られます．

```math
\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)$ はパラメータ更新で利用します．

## 対数尤度関数

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

```math
\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}$ で偏微分します．

```math
\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$ 回繰り返します．
サンプリングした結果を使ってパラメータを更新します．

```math
\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で実装したプログラムを掲載します．

```py3
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$ 次元の元画像を復元しています．

<img width="600" alt="rbm.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/b35d1d86-d70a-cb98-090d-b990f6e6b0e6.png">


# おわりに

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