Gaussian-Bernoulli Restricted Boltzmann Machine

先のCD learningの考察によって、ようやく意味がわかって、入力が連続値タイプのGaussian-Bernoulli RBMを実装できた。ちょこちょこ動いて、とってもかわいい!

GRBMDemo

Hinton先生のありがたいPracticeはほとんど反映されていない超簡易実装R版。混合分布推定に比べるとかなりシンプルで驚いた。

# Gaussian Bernoulli Restricted Boltzmann Machine

library(mvtnorm)

# 1. data generation
N = 200
D = 2

mu = list(NULL)
mu[[1]] = c(-10,-10)
mu[[2]] = c(-10,0)
mu[[3]] = c(-10,10)
mu[[4]] = c(0,-10)
mu[[5]] = c(0,0)
mu[[6]] = c(0,10)
mu[[7]] = c(10,-10)
mu[[8]] = c(10,0)
mu[[9]] = c(10,10)

S = diag(D)

K = length(mu)
X = matrix(0, N, D)
Z = numeric(N)

for(n in 1:N){
	Z[n] = sample(1:K, 1)
	X[n,] = rmvnorm(1, mu[[Z[n]]], S)
}

plot(X[,1], X[,2])

# 2. GRBM modeling

# 2-1. set up data, hyperparameters and function
v0 = scale(X)    # data set
L = 200          # number of hidden nodes
eta = 1e-2       # learning rate
itmax = 200      # max number of iteration

sigmoid = function(x){
	1.0/(1.0 + exp(-x))
}

# 2-2. set up parameters & variables
W = matrix(rnorm(D*L, 0, 1e-2), D, L)
b = rnorm(D,0,1e-2)
c = rnorm(D,0,1e-2)
R = diag(D) * 1e-2

# 2-3. execute learning
par(mfrow=c(2,2))
par(cex=0.6)
par(mar=rep(2,4))

rss = rep(NA, itmax)

for(it in 1:itmax){

	# a. calculate reconstructed samples and probability of hidden nodes; h0, v1, h1
	h0 = matrix(0, N, L)
	v1 = matrix(0, N, D)
	h1 = matrix(0, N, L)

	Rinv = solve(R)

	for(n in 1:N){
		h0[n,] = as.double(runif(L) < sigmoid(c + t(v0[n,]) %*% Rinv %*% W))
		v1[n,] = rmvnorm(1, b + W %*% h0[n,], R)
		h1[n,] = sigmoid(c + t(v1[n,]) %*% Rinv %*% W)	
	}
	
	# b. calculate statistics
	Evht0 = t(v0) %*% h0
	Evht1 = t(v1) %*% h1
	Ev0 = apply(v0, 2, mean)
	Ev1 = apply(v1, 2, mean)
	Eh0 = apply(h0, 2, mean)
	Eh1 = apply(h1, 2, mean)

	Evvt0 = Ev0 %*% t(Ev0) + cov(v0)
	Evvt1 = Ev1 %*% t(Ev1) + cov(v1)

	rss[it] = sum((v0 - v1)^2)

	# c. calculate derivatives of the contrastive divergence (CD)
	dW = - Rinv %*% (Evht0 - Evht1)
	db = - Rinv %*% (Ev0 - Ev1)
	dc = - (Eh0 - Eh1)

	tmp1 = -0.5 * Rinv %*% (Evvt0 - 2*Ev0%*%t(b) - 2*Evht0%*%t(W)) %*% Rinv
	tmp2 = -0.5 * Rinv %*% (Evvt1 - 2*Ev1%*%t(b) - 2*Evht1%*%t(W)) %*% Rinv
	dR = diag(diag(tmp1 - tmp2))

	# d. gradient-descent update
	W = W - (eta * max(abs(W)) / max(abs(dW))) * dW
	b = b - (eta * max(abs(b)) / max(abs(db))) * db
	c = c - (eta * max(abs(c)) / max(abs(dc))) * dc
	R = R - (eta * max(abs(R)) / max(abs(dR))) * dR

	# e. plot the progress
	plot(v0[,1], v0[,2], col="gray", main=sprintf("it=%d, v0=gray, v1=blue",it), xlim=c(-2,2), ylim=c(-2,2))
	points(v1[,1], v1[,2], col="blue")

	plot(1:itmax, rss, type="l", main=sprintf("it=%d, rss",it), ylim=c(0,200))
	abline(h=36, col="red")
}


rev.1 プログラムにミスがあったので修正。勾配法の更新量を調整して安定化。

MLE vs. Contrastive Divergence Learning

Restricted Boltzman Machineのプログラムを組んでみたのだけれど、上手く学習ができない。

真面目にCD learningの勉強をしなければと思い立ち、1次元ガウスの場合で計算をしてみた。以下、その結果。


下図は、ガウスの平均と精度パラメータの学習過程である。

KLvsCD_fig1

Gradient descent版の最尤推定と比較してみると、学習が遅いことがわかる。しかし、下図のように尤度でみると、毎イテレーションで、基本的に単調増加はしている。ただし、保証はされていない。

KLvsCD_fig2

おまけに、今回真面目に計算して分かったのだが、Contrastive Divergenceとは経験分布(データ)と真のモデル分布とのKullback Leibler Divergenceの超荒い近似であるようだ。RBMでは、潜在変数のGibbs Samplingが必要なため、さらにその近似精度が下がるであろうことは推測できる。

すなわち、学習には忍耐が必要・・・ということか。

上記の詳細はこちら


ということで、CD学習が正しく機能することは分かったので、もう少し頑張ってみようかな。