ゆるふわ Restricted Boltzmann Machine

Deep Learningで用いられるらしいということで, Restricted Boltzmann Machine(RBM)について調べたので概要とPythonによる実装例をまとめた. 主にAn Introduction to Restricted Boltzmann Machinesを参考にしているので, 数式の詳細はそちらをあたって欲しい.

概要

Restricted Boltzmann Machineは分布をモデル化するアルゴリズム.
目的は, 与えられた観測変数の集合 $\{\mathbf{x}^i\}_{i=1}^{n}$ からその確率分布 $P(\mathbf{x})$ を求めることだ.

RBMでは観測変数 $\mathbf{x}$ の他に隠れ変数 $\mathbf{h}$ を導入する. そして, $P(\mathbf{x}, \mathbf{h})$ を求めたのち, 周辺化によって $P(\mathbf{x})$ を求める.

RBMでは $P(\mathbf{x}, \mathbf{h})$ は次のように表される.

$P(\mathbf{x}, \mathbf{h}) = \frac{1}{Z}\exp(-E(\mathbf{x}, \mathbf{h}))$
ただし, $E(\mathbf{x}, \mathbf{h}) = -\mathbf{b}^T \mathbf{x} - \mathbf{c}^T\mathbf{h} - \mathbf{h}^TW\mathbf{x}$ .

Eはエネルギーと呼ばれるもので, エネルギーが小さい状態程起こる確率は高くなる. 水は低きに流れる的な精神. また, $\mathbf{b}, \mathbf{c}, W$ の3つはこのモデルのパラメーターで, これらをデータから最尤推定で求める(詳細は後述). そして $Z$ は正規化定数. 全部足しあわせたときに1になるように調整するためのもの.

最尤推定でモデルパラメーターを求めた後は, 得た分布 $P(\mathbf{x},\mathbf{h})$ を使って適当にやりたいことをやる. RBM完！

解釈

RBMは次のようなグラフィカルモデルで表される.

上側の円はそれぞれ $\mathbf{x}$ の１次元目の要素, ２次元目の要素,...に対応している. 下側の円はそれぞれ $\mathbf{h}$ の各次元の要素に対応している.
一応, 図があった方が見栄えがいいから載せてみたけど, この図はあんまり重要ではない. グラフィカルモデルを知らない人は, 次の２つの独立性

$P(\mathbf{h} | \mathbf{x}) = \prod_{i=1}^M P(h_i | \mathbf{x})$
$P(\mathbf{x} | \mathbf{h}) = \prod_{i=1}^D P(x_i | \mathbf{h})$

がRBMでは成り立つことだけ認識すればOK. なお, $M,D$ はそれぞれ $\mathbf{h},\mathbf{x}$ の次元数で, $h_i,x_i$ はそれぞれのi番目の要素を指す.

この文献では, RBMにおいて $\mathbf{h}$ の各要素 $h_i$ は, $\mathbf{x}$ に対してM回分類を行った際のi番目の分類結果を表していると解釈できると述べている. hの各要素が0,1の値を取る場合その結果の組み合わせは $2^M$ 通りとなり, 豊かな表現力を持つ( $\mathbf{h}$ を指定すると細かい $\mathbf{x}$ の集合が指定できる).

また, 先ほどの文献に示された他の解釈として,

$P(\mathbf{x}) = \sum_{\mathbf{h}} P(\mathbf{x} | \mathbf{h}) P(\mathbf{h})$

と書けることから, $P(\mathbf{x} | \mathbf{h})$ がガウス分布のとき, これは $2^M$ 個のガウス分布を用いた混合ガウス分布と解釈できる(ちょっとエネルギーの定義をいじったら $P(\mathbf{x} | \mathbf{h})$ がガウス分布になるっぽい. Greedy Layer-Wise Training of Deep Networksを読んでください).

たぶんこれが言いたいのは, かなりの多峰性を持つ分布も表現できる可能性があるということだろう. Kernel Density Estimationみたいなモデルだとデータ点の数しか峰は持てない.

応用

正直, あんまり調べてないからわからない. 一つの例としてはRBMはDeep Learningの部品として使われている. そこでは, $P(\mathbf{h} | \mathbf{x})$ を最大にする $\mathbf{h}$ が $\mathbf{x}$ のより高次の表現とみなされている(ex. $\mathbf{x}$ はピクセル値を表し $\mathbf{h}$ はエッジの存在を表す的な).

また, 分類問題にも使うことができて, 使い方の一つの例として $\mathbf{x}$ 側にラベルを表現する次元を一つ足してモデルを学習するらしい. そして, 新しいデータ点に対しては, $P(label | \mathbf{x})$ を最大にするようなラベルを割りあてるっぽい.

詳細

ここから先は実際にどうやってRBMのパラメーターを学習をするかという話. 概要を思い出すか見返して欲しい. なお, ここからは $\mathbf{x},\mathbf{h}$ の各要素が0,1の値を取る場合を考える. 連続値の場合はまだ勉強してない.

最尤推定では, 対数尤度関数

$\ln{L(\{\mathbf{x}^i\}_{i=1}^{n} | \theta)} = \sum_{i=1}^{n}\ln P(\mathbf{x}^i | \theta)$

を最大にするようなパラメーター $\theta$ を求める. ここでは, $\theta=\{W, \mathbf{b}, \mathbf{c}\}$ .

対数尤度を最大化するために確率的勾配降下法を用いる. そのためデータ点あたりの対数尤度の各パラメーターに関する偏微分を計算する必要がある. 例えばm番目のデータにおける $W$ のi行j列目に関する偏微分は

$\frac{\partial \ln L(\mathbf{x^m}|\theta)}{\partial w_{ij}} = P(h_i=1 | \mathbf{x}^m)x^m_j - \sum_{\mathbf{x' \in \{0,1\}^D}} P(\mathbf{x'})P(h_i=1|\mathbf{x'})x_j'$

となる. これは頑張って計算すれば普通に導出できる. 二つの変数 $\mathbf{x}, \mathbf{x}'$ の違いに気をつけて欲しい. $P(h_i=1 | \mathbf{x})$ も同時分布を周辺化することで求めることが出来る.

問題は最後の項が全ての $\mathbf{x}$ に関する総和となっているところ. 明らかに $\mathbf{x}$ の次元数が大きくなると計算不可能になる. これを近似するのがRBMの学習における要点である.

真っ先に考えられるのは, Gibbsサンプリングを用いて分布 $P(\mathbf{x})$ から適当に何点かサンプリングし, それらに関して平均を取ることである. Gibbsサンプリングに関しては適当に手順だけ.

あるデータ点 $\mathbf{x}^{m}$ を元に, 現在のパラメーターで定められる分布 $P(\mathbf{h} | \mathbf{x}^{m})$ から $\mathbf{h}^{(0)}$ をサンプリングする

次に $P(\mathbf{x} | \mathbf{h}^{(0)})$ から, $\mathbf{x}^{(1)}$ をサンプリングする

$P(\mathbf{h} | \mathbf{x}^{(1)})$ から $\mathbf{h}^{(1)}$ をサンプリングする

$P(\mathbf{x} | \mathbf{h}^{(1)})$ から, $\mathbf{x}^{(2)}$

この手順をk回繰り返して $\mathbf{x}^{(k)}$ を得る

このようにして得た $\mathbf{x}^{(k)}$ は, kを大きく取ると始点に関係なく分布 $P(\mathbf{x})$ に従う(筈). 詳しくはPRMLのサンプリングの章を読んで欲しい. 結構わかりやすい. なお, 分布 $P(\mathbf{x},\mathbf{h})$ から直接サンプリングするのは, 正規化定数Zを求めるのが困難なため難しい.

サンプリングで分布 $P(\mathbf{x})$ に従う観測点の集合が得られれば上記の総和を近似できる. 他のパラメーターも同様にして勾配を求めてやればよい. その勾配を用いてパラメーターを最適化することができる.

Contrastive Divergence learning

だが, 残念ながら上の方法は計算量が非常に大きい. 上記のGibbsサンプリングは出発点 $\mathbf{x}^{m}$ への依存性を断ち切り, $\mathbf{x}^{(k)}$ が $P(\mathbf{x})$ に従うようにするため, 比較的大きなkを取る必要があると考えられる.

そこで, Contrastive Divergence learningは比較的小さなkを用いる. ふざけんなって感じだがk=1の場合でも, 経験上多くの場合うまくいくことが知られているっぽい. その上, 複数の点をサンプリングして用いるのではなく, １点のみをサンプリングして総和を近似する. そのため更新式は次のようになる.

$\frac{\partial \ln L(\mathbf{x}^m|\theta)}{\partial w_{ij}} = P(h_i=1 | \mathbf{x}^m )x^m_j - P(h_i=1|\mathbf{x^{(k)}})x^{(k)}_j'$

このようにして得た更新式は厳密には対数尤度を最大化せず, Contrastive Divergenceという対数尤度に対してバイアスが乗った尺度を最大化することになる. そのため, この手法はContrastive Divergence learningと呼ばれる.

実装

細かい議論は後に回し, 実装例を示す. とりあえず必要な数式を以下に列挙しておく.

$\frac{\partial \ln L(\mathbf{x}^m|\theta)}{\partial w_{ij}} = P(h_i=1 | \mathbf{x}^m )x^m_j - P(h_i=1|\mathbf{x^{(k)}})x^{(k)}_j'$
$\frac{\partial \ln L(\mathbf{x}^m|\theta)}{\partial b_{j}} = x^{m}_j - x^{(k)}_j$
$\frac{\partial \ln L(\mathbf{x}^m|\theta)}{\partial c_{i}} = P(h_i=1 | \mathbf{x}^m) - P(h_i=1 | \mathbf{x}^{(k)})$

また,

$P(h_i=1 | \mathbf{x}) = \sigma \Bigl(\sum_{j=1}^{M} W_{ij}x_j + c_i \Bigr)$
$P(x_j=1 | \mathbf{h}) = \sigma \Bigl(\sum_{i=1}^{D} W_{ij}h_i + b_j \Bigr)$

これらの式を使えば, サンプリングや勾配の計算を行える.

実装例では, 入力次元D=4, 隠れ層の次元M=3としている. 今回のプログラムでは, まず適当に分布 $P(\mathbf{x})$ を生成し, そこからデータ点を発生させ, そのデータ点から $P(\mathbf{x})$ を学習させている.

また, 学習の経過を見るために逐一対数尤度関数を計測している(データの次元が低いため可能).
また初期値の決定法などはA Practical Guide to Training Restricted Boltzmann Machinesを参考にしている.
以下はソースコードだが, 先にその下の結果を読んだ方がいいと思う.

import scipy as sp

sp.set_printoptions(precision=4, suppress=True)

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

class RBM:
    """
    今回は入力, 隠れ層ともに２値のRBMを扱う

    v : 観測ベクトル(xじゃなくてすみません)
    h : 隠れ層を表すベクトル
    W[i,j] : i番目の隠れ層とj番目の入力との組に対応する重み

    E(v,h) = -dot(b,v) - dot(c,h) - sp.dot( sp.dot(h,W), x )
    というエネルギー関数を持つモデル
    """
    def __init__(self, dim_in, dim_hidden, V, k=1, is_PCD=False):
        """
        @type V : ２次元のsp.array
        @param V : 観測ベクトルの集合. V[j] はi番目のデータ
        """
        self.dim_hidden = dim_hidden
        self.dim_in = dim_in
        
        #パラメーターの初期化はA Practical Guide to Training Restricted Boltzmann Machinesを参照せよ
        self.W = 0.01*sp.random.randn(dim_hidden, dim_in)
        self.b = sp.log(sp.mean(V, axis=0) / (1-sp.mean(V, axis=0)))

        #self.b = sp.zeros(dim_in)
        self.c = sp.zeros(dim_hidden)
        self.is_PCD = is_PCD

        self.V = V
        self.num_in = len(self.V)
        self.k = k

    def train(self, mini_batch_size=10, learning_rate=0.01):
        """
        全てのデータに関して１度ずつ勾配を求めて, パラメーターを更新する
        
        mini_batch_size : パラメーター更新１回あたりにいくつデータ点を用いるか
        learning_rate : どれぐらいの割合で, Wの10^(-3)倍ぐらいのオーダーがいいらしい

        mini_batch_sizeがどれぐらいがいいかは
        A Practical Guide to Training Restricted Boltzmann Machines
        で議論されている.

        小さいとデータ点あたりの計算効率があがるが（特にPythonでは）,
        大きすぎると, 学習効率が悪くなるらしい.
        これはRBMに特有の問題というよりは,
        確率的勾配法一般の問題なのかな？
        """
        #データをmini_batch_sizeごとに切り分け
        V_batches = [self.V[i*mini_batch_size:(i+1)*mini_batch_size] for i in range(self.num_in / mini_batch_size)]
        V_sample = V_batches[0]

        for V_batch in V_batches:
            if self.is_PCD:
                V_sample = self.GibbsSampling(V_in=V_sample, k=self.k)
            else:
                V_sample = self.GibbsSampling(V_in=V_batch, k=self.k)

            grad_W, grad_b, grad_c = self.grad(V_batch, V_sample)

            self.W += learning_rate * grad_W
            self.b += learning_rate * grad_b
            self.c += learning_rate * grad_c
            #print sp.sum((grad_W-fake_W)**2), sp.sum((grad_b-fake_b)**2), sp.sum((grad_c-fake_c)**2)

    def grad(self, V_in, V_sample):
        """
        1点のサンプリングデータを元にパラメーター更新を行う.
        複数のデータ点に関して同時に計算してるから読みづらい

        V_in : このデータ点の集合に関して勾配を計算する
        V_sample : サンプリングした点の集合 (V_inと同じサイズ)
        """
        #p_H_given_v[i,m] =  (H_i = 1 | V_in[m])
        p_H_given_V_in = sigmoid(sp.dot(self.W, V_in.T) + self.c.reshape(self.dim_hidden, 1))
        p_H_given_V_sample = sigmoid(sp.dot(self.W, V_sample.T) + self.c.reshape(self.dim_hidden, 1))
        grad_W = sp.dot(p_H_given_V_in, V_in) - sp.dot(p_H_given_V_sample, V_sample)
        grad_b = sp.sum(V_in, axis=0) - sp.sum(V_sample, axis=0)
        grad_c = sp.sum(p_H_given_V_in, axis=1) - sp.sum(p_H_given_V_sample, axis=1)

        num_in = len(V_in)
        grad_W = grad_W / num_in
        grad_b = grad_b / num_in
        grad_c = grad_c / num_in

        return grad_W, grad_b, grad_c

    def grad_one(self, v_in, v_sample):
        """
        gradの１点ずつバージョン.
        可読性は高いが遅い. デバッグ用
        """
        p_h_given_v_in = sigmoid(sp.dot(self.W, v_in) + self.c)
        p_h_given_v_sample = sigmoid(sp.dot(self.W, v_sample) + self.c)
        grad_W = sp.outer(p_h_given_v_in, v_in) - sp.outer(p_h_given_v_sample, v_sample)
        grad_b = v_in - v_sample
        grad_c = p_h_given_v_in - p_h_given_v_sample

        return grad_W, grad_b, grad_c

    def GibbsSampling(self, V_in, k):
        """
        V_inのそれぞれの点から, GibbsSamplingにおけるstep kで得られる
        点を計算する

        V_sampleのそれぞれの点は, 確率を使ったbinary値ではなく,
        直接, 確率の値を使った方が収束性能がよいという議論がある
        """
        num_in = len(V_in)
        V_now = V_in
        
        for idx_k in range(k):
            p_H_given_V = sigmoid(sp.dot(self.W, V_now.T) + self.c.reshape(self.dim_hidden, 1))
            #H[i] : i 番目のサンプル
            H = p_H_given_V.T > sp.random.uniform(size=(num_in, self.dim_hidden))
            H = H.astype(float)

            #Vの生成
            #p_V_given_H[m,i] =  (V_i = 1 | H[m]) (dim_in * m)
            p_V_given_H = sigmoid(sp.dot(H, self.W) + self.b)

            V_now = p_V_given_H > sp.random.uniform(size=(num_in, self.dim_in))
            V_now = V_now.astype(float)
            
        return V_now

    def sampling_one(self, v, k):
        """
        ある１点からサンプリングするデバッグ用
        """
        p_h_v = sigmoid(sp.dot(self.W, v) + self.c)
        h = p_h_v > sp.random.uniform(size=dim_hidden)
        p_v_h = sigmoid(sp.dot(h, self.W) + self.b)
        v_sample = p_v_h > sp.random.uniform(size=dim_in)

        return v_sample, p_h_v, p_v_h
            
    def log_likelihood(self, V):
        """
        尤度関数
        P(self.V | self.W, self.b, self.c)を計算する
        これは, 分配関数の計算にvの全ての場合に関しての総和を必要とするため,
        一般には計算できない. vの状態数が小さいときのみ有効
        """
        all_v = []
        
        for i in range(2**self.dim_in):
            all_v.append([(i>>con) % 2 for con in range(self.dim_in)])

        l = sp.exp(sp.dot(all_v, self.b))
        r = sp.multiply.reduce(1+sp.exp(
            self.c + sp.dot(all_v, self.W.T)),axis=1)
        all_prob = l*r
        Z = sum(all_prob)
        #print all_prob / Z
        self.all_prob = all_prob / Z

        #対数尤度を求める
        result = -self.num_in*sp.log(Z) + sp.sum(sp.dot(V, self.b)) + \
                 sp.sum(sp.log(1 + sp.exp(self.c + sp.dot(V, self.W.T)) ))

        return result

if __name__ == "__main__":
    #観測次元4, 隠れ変数3でやってみる
    dim_in = 4
    dim_hidden = 3
    #学習データの点数
    num_data = 1000
    #学習率
    lerning_rate=0.01

    #入力データの全パターンを列挙
    all_v = []
    for i in range(2**dim_in):
        all_v.append([(i>>con) % 2 for con in range(dim_in)])
    all_v = sp.array(all_v)

    #分布を生成
    dist = sp.random.uniform(size=len(all_v))
    dist = dist / sum(dist)
    #累積分布関数を求める
    dist_accum = sp.add.accumulate(dist)

    #データの生成
    sample_idx = sp.searchsorted(dist_accum, sp.random.uniform(size=num_data))
    V_in = all_v[sample_idx]

    rbm_k1 = RBM(dim_in, dim_hidden, V_in)
    rbm_k5 = RBM(dim_in, dim_hidden, V_in, k=5)
    rbm_pcd = RBM(dim_in, dim_hidden, V_in, is_PCD=True)

    rbms = [rbm_k1, rbm_k5, rbm_pcd]
    labels = ["k1", "k5", "pcd"]
    
    import pylab
    scores = []

    #3000回同じデータセットを使う
    for i in range(5000):
        tmp = []
        
        for idx_rbm, rbm in enumerate(rbms):
            rbm.train()
            score = rbm.log_likelihood(V_in)
            tmp.append(score)

        scores.append(tmp)

    #分布の表示
    print "true:"
    print dist
    print "estimated:"
    print rbm_k1.all_prob

    scores = sp.array(scores)

    for idx_rbm in range(len(rbms)):
        pylab.plot(scores[:,idx_rbm], label=labels[idx_rbm])
    pylab.legend()
    pylab.show()

次のグラフは横軸が学習の繰り返し回数. 縦軸が対数尤度である. 青線がサンプリングにおいてk=1とした場合, 緑はk=5, 赤は後で説明するPersisten Contrastive Divergenceを使った場合. だんだん, パラメーターがデータにフィッティングし, 対数尤度が上昇していることがわかる. ちなみにそれぞれの手法の上下関係は試行ごとに変化するため, この図で優劣を判断するべきではない.

また, 次の表示結果は, 上はそれぞれx=[0,0,0,0], [0,0,0,1]...になる確率で, 下はそれに対応する推定した確率分布である.

true:
[ 0.056 0.0119 0.0591 0.019 0.0831 0.0473 0.0851 0.0855 0.0333
0.1102 0.1178 0.0156 0.0422 0.0763 0.1157 0.0421]
estimated:
[ 0.0556 0.0129 0.0553 0.0157 0.0819 0.0458 0.0869 0.0716 0.0332
0.1205 0.1208 0.0161 0.0343 0.098 0.1115 0.04 ]

それなりに大体正しく推定できていることがわかる.

この例では, 最適化の尺度として対数尤度を用いているが現実には計算不可能であることが多い. そのため代わりの尺度として, それぞれのデータ点 $\mathbf{x}^m$ から $P(\mathbf{h} | \mathbf{x}^m)$ に従い $\mathbf{h}^m$ を生成し, そこから更に $P(\mathbf{x} | \mathbf{h}^m)$ に従い $\hat{\mathbf{x}}^m$ を生成した上で, それと $\mathbf{x}^m$ との距離を測るreconstruction errorなるものがよく使われるらしい.

また, 正規化定数Zも計算不能であるため, もしそれが必要な場合は(条件付き分布のみを使うなら不要)適当にサンプリングして求めるんだろう. たぶん.

バリエーションと議論

この節の議論の多くはA Practical Guide to Training Restricted Boltzmann Machinesのものである.

一番重要な議論は, Contrastive Divergenceの適切性だろう. これは, kが大きくなれば対数尤度と一致するが, そうでない場合は, バイアスが存在する. そのため, 試行回数を増やすと対数尤度が減少する場合がある. 著しい不具合があるときはkを大きくしたりサンプリング方法を工夫したりといった措置が求められる.

工夫の一つは上のプログラムにも実装されている, サンプリングをする上で, データ点を初期値とするのではなく, 前回のサンプリング結果を初期値とする手法があり, それはPersistent Contrastive Divergence learningと呼ばれている. かなり性能がいいらしいが, 上のプログラムではそれ程差が出ていない. 問題が簡単すぎるのか, パラメーターが不適切なのか, 実装が間違っているのか・・・. より性能のいいParallel Temperingという手法もあるらしい.

上と少し関係するのがReconstruction Errorを学習進度の尺度として用いる場合, それは対数尤度と一致しないため, これをあんまり信用しすぎてはいけないらしい. パラメーターの値が大きくなり, サンプリングのmixing rate(定常状態に収束していく速度)が遅くなると, Reconstruction Errorは低く見積もられてしまうようだ(要確認).

また, 確率勾配法を使うのに１回のパラメーター更新あたりどれぐらいのデータ点を使うか(mini-batch size)に関しても議論がある. 複数のデータを同時に使う場合は行列演算により, GPUやPythonでは効率的に計算できる. しかし, 学習効率はあまり良くないらしいので, 大きくしすぎてはいけないらしい.

他には, サンプリングした値で総和を近似するのではなく, 確率値を代わりに使う, つまり $P(\mathbf{x} | \mathbf{h}^{(1)})$ から0,1のベクトルをサンプリングする代わりに, 確率の値をベクトルの要素とする技があるらしく, 学習効率をあげるらしいが, 全然動かなかった. たぶん手法の解釈が間違ってるんだと思う.