Hatena::ブログ(Diary)

Risky Dune

2012-12-12 ゆるふわ Restricted Boltzmann Machine

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

概要

Restricted Boltzmann Machineは分布をモデル化するアルゴリズム.
目的は, 与えられた観測変数の集合からその確率分布を求めることだ.

RBMでは観測変数の他に隠れ変数を導入する. そして, を求めたのち, 周辺化によってを求める.

RBMではは次のように表される.


ただし, .

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

最尤推定でモデルパラメーターを求めた後は, 得た分布を使って適当にやりたいことをやる. RBM完!

解釈

RBMは次のようなグラフィカルモデルで表される.
f:id:saket:20121212184655p:image:w360
上側の円はそれぞれの1次元目の要素, 2次元目の要素,...に対応している. 下側の円はそれぞれの各次元の要素に対応している.
一応, 図があった方が見栄えがいいから載せてみたけど, この図はあんまり重要ではない. グラフィカルモデルを知らない人は, 次の2つの独立性


がRBMでは成り立つことだけ認識すればOK. なお, はそれぞれの次元数で, はそれぞれのi番目の要素を指す.

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

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

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

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

応用

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

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

詳細

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

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

を最大にするようなパラメータを求める. ここでは, .

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

となる. これは頑張って計算すれば普通に導出できる. 二つの変数の違いに気をつけて欲しい. も同時分布を周辺化することで求めることが出来る.

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

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

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

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

Contrastive Divergence learning

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

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

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

実装

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



また,


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

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

また, 学習の経過を見るために逐一対数尤度関数を計測している(データの次元が低いため可能).
また初期値の決定法などは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:
    """
    今回は入力, 隠れ層ともに2値の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 : 2次元の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):
        """
        全てのデータに関して1度ずつ勾配を求めて, パラメーターを更新する
        
        mini_batch_size : パラメーター更新1回あたりにいくつデータ点を用いるか
        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の1点ずつバージョン.
        可読性は高いが遅い. デバッグ用
        """
        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):
        """
        ある1点からサンプリングするデバッグ用
        """
        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を使った場合. だんだん, パラメーターがデータにフィッティングし, 対数尤度が上昇していることがわかる. ちなみにそれぞれの手法の上下関係は試行ごとに変化するため, この図で優劣を判断するべきではない.

f:id:saket:20121212184654p:image:w720
また, 次の表示結果は, 上はそれぞれ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 ]

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

この例では, 最適化の尺度として対数尤度を用いているが現実には計算不可能であることが多い. そのため代わりの尺度として, それぞれのデータ点からに従いを生成し, そこから更にに従いを生成した上で, それととの距離を測る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は低く見積もられてしまうようだ(要確認).

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

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

まとめ

とりあえずRestricted Boltzmann Machineの概要はわかったので, 次はDeep Belief Netを実装しようかどうか悩んでいる. 先に面白そうなデータセットを探すべきかもしれない.

RBMに対して調べないといけないなぁと思ってるのは, 入力を連続値にしたケースは扱えるようになっておきたい. 他にはもうちょっとましなデータで実験すべきかなぁとも思う.

参考文献