Gaussian-Bernouili RBMって
Restricted Boltzmann Maschine(RBM)と言えば可視層と隠れ層がBernouili分布に従いますよね。
本題のGaussian-Bernouili RBM（GBRBM）は可視層がGaussian分布で隠れ層がBernouili分布に従うRBMなわけです。
GBRBM学習って分散の学習が上手くいかない〜、とか、分散を1に固定して学習しない記事や実装コードばかり目にするので、
GBRBMを使おうとしてた時に相当困りました・・・（４年前）。
ということでコードを混ぜてGBRBMについて！！！
分散の学習について
悩んでる人は当時の僕だけじゃないはずですが、
一応前置きすると、可視特徴量の分散をきちんと再現できるわけでは無いです。
あくまで「途中でパラメータが発散してしまう！」などの、学習が破綻するようなレベルの問題をクリアするためのお話です。
結論から言うと、”学習時にCD法を行う際、きちんと隠れ層をサンプリングしてから可視層の分布を推定する”ってことです。
多くの記事やコードでは
$v_0$ ⇒ $p(h_0)$ ⇒ $p(v_1)$ ⇒ $p(h_1)$
のように、CD法を行ってます。
これを
$v_0$ ⇒ $p(h_0)$ ⇒ $h_0$ ⇒ $p(v_1)$ ⇒ $v_1$ ⇒ $p(h_1)$
のように、ちゃんと$h_0$をサンプリングしてから$p(v_1)$を算出すれば学習もちゃんと進みます。
「RBMの隠れ層の値は 1 or 0 程度の情報量しか持ってない」という前提で可視層が予測されるので
学習時はちゃんとサンプリングした方が良いそうです。
ただし、学習が終わってしまえば、デコード時は隠れ層の期待値からデコードして大丈夫（というかそうすべき？）です。
初期化
import numpy as np
class GBRBM(object):
def __init__(self, I, J):
self.W = np.random.randn(I, J)
self.b = np.zeros(I)
self.c = np.zeros(J)
self.z = np.zeros(I)
Wは接続重み、bとcはそれぞれ可視層のバイアスと隠れ層のバイアスです。
分散$\sigma^2_i$の学習ですが、非負制約を満たさなくてはいけないので
$$
z_i = \log\sigma^2_i
$$
として、$\sigma^2_i$の代わりに$z_i$を学習します。
ネットワーク定義
def propup(self, v):
v_ = v / np.exp(self.z)
tmp = np.dot(v_, self.W) + self.c
ph = self.sigmoid(tmp)
return ph
def propdown(self, h):
v_mean = np.dot(h, self.W.T) + self.b
return v_mean
def contrastive_divergence(self, v, k=1):
sigma = np.sqrt(np.exp(self.z))
ph0 = self.propup(v) # p(h|v)
phk = ph0
for step in range(k):
hk = np.random.binomial(size=phk.shape, n=1, p=phk) # sampling h
pvk = self.propdown(hk) # p(v|h)
vk = np.random.randn(pvk.shape[0], pvk.shape[1])*sigma + pvk
phk = self.propup(vk) # p(h|v)
return ph0, vk, phk
propupは$p(h|v)$を算出するブロックで、propdownは$p(v|h)$を算出するブロックです。
contrastive_divergenceはCD法によるサンプリングブロックで、
さっきも書きましたが$p(v|h)$を算出する際にphk（$p(h|v)$のこと）を使うのではなく
ちゃんとサンプリングしたhkを使います。
ただ、勾配計算にはどうせ期待値を使うのでサンプリングした値ではなく$p(h|v)$を返します。
勾配計算
def get_grad(self, v):
N = v.shape[0]
#--- CD_1 ---#
v0 = v
h0, v1, h1 = self.contrastive_divergence(v0)
var = np.exp(self.params['z'])
v0_, v1_ = v0/var, v1/var
W, b = self.W, self.b
#--- Calculate gradient params ---#
vw0 = np.dot(v0_.T, h0)
vw1 = np.dot(v1_.T, h1)
dW = (vw0 - vw1) / N
db = np.mean(v0_ - v1_, axis=0)
vz1 = 0.5*((v0 - b)**2) - np.dot(h0, W.T)*v0
vz2 = 0.5*((v1 - b)**2) - np.dot(h1, W.T)*v1
dz = np.mean((vz1 - vz2)/var, axis=0)
return dW, db, dz
ここはどの記事でもコードでもほぼ一緒だと思います。
全体コード
import numpy as np
class GBRBM(object):
def __init__(self, I, J):
self.W = np.random.randn(I, J)
self.b = np.zeros(I)
self.c = np.zeros(J)
self.z = np.zeros(I)
def propup(self, v):
v_ = v / np.exp(self.z)
tmp = np.dot(v_, self.W) + self.c
ph = self.sigmoid(tmp)
return ph
def propdown(self, h):
v_mean = np.dot(h, self.W.T) + self.b
return v_mean
def contrastive_divergence(self, v, k=1):
sigma = np.sqrt(np.exp(self.z))
ph0 = self.propup(v) # p(h|v)
phk = ph0
for step in range(k):
hk = np.random.binomial(size=phk.shape, n=1, p=phk) # sampling h
pvk = self.propdown(hk) # p(v|h)
vk = np.random.randn(pvk.shape[0], pvk.shape[1])*sigma + pvk
phk = self.propup(vk) # p(h|v)
return ph0, vk, phk
def get_grad(self, v):
N = v.shape[0]
#--- CD_1 ---#
v0 = v
h0, v1, h1 = self.contrastive_divergence(v0)
var = np.exp(self.params['z'])
v0_, v1_ = v0/var, v1/var
W, b = self.W, self.b
#--- Calculate gradient params ---#
vw0 = np.dot(v0_.T, h0)
vw1 = np.dot(v1_.T, h1)
dW = (vw0 - vw1) / N
db = np.mean(v0_ - v1_, axis=0)
vz1 = 0.5*((v0 - b)**2) - np.dot(h0, W.T)*v0
vz2 = 0.5*((v1 - b)**2) - np.dot(h1, W.T)*v1
dz = np.mean((vz1 - vz2)/var, axis=0)
return dW, db, dz
def train(self, v, nbatch=256, nepoch=100, lr=0.01):
N = v.shape[0]
for epoch in range(n_epoch):
perm = np.random.permutation(N)
for i in range(0, N-nbatch, nbatch):
v_batch = np.asarray(v[perm[i:i+n_batch]])
dW, db, dc, dz = self.get_grads(v_batch)
self.W -= lr * dW
self.b -= lr * db
self.c -= lr * dc
self.z -= lr * dz
def sigmoid(self, x):
return 1./(1. + np.exp(x))