0. はじめに

RBMはDeep Learningで多層化する場合に、事前学習として用いられる手法です。
ソースコードはtheanoで公開されています。
http://deeplearning.net/tutorial/rbm.html
今回はこのコードをchainerに移植してみました。

RBMの場合、誤差関数はFree Energyの関数となり複雑(?)ですが
逆誤差伝搬のため重み行列、バイアスの変数で微分すると綺麗な計算式となり、
Contrastive Divergence(CD)処理の値から計算できます。
chainerで実装する場合、わざわざ誤差関数(loss)を作成して逆誤差伝搬(loss.backward)するよりも、
単純にcupy(numpy+GPU)で実装したほうが単純となりますが、あえて誤差関数を使って実装してみました。
cupy(numpy+GPU)だけのソースコードはまた機会がある場合に公開します。
といっても以下のソースコードもCDの処理はcupyを使って高速化してますが。。。

# -*- coding: utf-8 -*-

import os,sys
import argparse
import time

import numpy as np

import chainer
from chainer import computational_graph
from chainer import cuda, Variable

import chainer.functions as F
import chainer.links as L
import chainer.optimizers as O


import pickle
import gzip

parser = argparse.ArgumentParser()
parser.add_argument('--gpu', '-g', default=-1, type=int,
                    help='GPU ID (negative value indicates CPU)')
parser.add_argument('--pcd', '-p', default=1, type=int,
                    help='pcd_flag')
parser.add_argument('--kcd', '-k', default=1, type=int,
                    help='cd-k')


''' GPUで計算する時は、cupy = numpy + GPUで処理する。 '''
args = parser.parse_args()
if args.gpu >= 0:
    cuda.check_cuda_available()
xp = cuda.cupy if args.gpu >= 0 else np

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

class RBM(chainer.Chain):
    '''
    n_visbile:可視層の次元
    n_hidden:隠れ層の次元
    l.W:重み行列
    l.b:hbias
    l.a:vbaias
    real: 可視層が0,1か実数なのかで場合分け。 初期層は実数を想定しているが、中間層は0,1しか考えていない。
    chainerで h = l(v)とかけるように一般的なRBMに出てくる重み行列を転置している。
    '''
    def __init__(self,n_visible,n_hidden,real=0,k=1,pcd_flag=0):
        super(RBM,self).__init__(
            l=L.Linear(n_visible,n_hidden),
            # chainer 1.5のmanualにこの使い方は非推奨とあったのでadd_param使用。
#            a=L.Parameter(xp.zeros(n_visible,dtype=xp.float32)),
        )
        self.l.add_param("a",(n_visible),dtype=xp.float32)

        ''' l.aを初期化、chainerの使用で初期化しないと学習しない？ '''
        self.l.a.data.fill(0)

        self.n_visible = n_visible
        self.n_hidden = n_hidden
        self.real = real
        self.k = k
        self.pcd_flag = pcd_flag

    def __call__(self,v_data,v_prev_data=None):

        batch_size = v_data.shape[0]

        if self.pcd_flag == 0:
            v_prev_data = v_data
        elif self.pcd_flag ==1 and v_prev_data is None:
            v_prev_data = v_data
        vh_data = self.constrastive_divergence(v_prev_data,self.k)
        v = Variable(v_data)
        vh = Variable(vh_data.astype(np.float32))

        '''
        http://deeplearning.net/tutorial/rbm.htmlの式(5)から
        誤差関数(loss)は(可視層のデータのfree energy)と(可視層のデータのCD-kのfree energy)の差となる。
        loss = (self.free_energy(v) - self.free_energy(vh)) / batch_size
        v->vhへの写像は学習に入れないため、CD-kの処理完了後にVariable化
        '''

        loss = (self.free_energy(v) - self.free_energy(vh)) / batch_size
        return loss

    def free_energy(self,v):
        ''' Function to compute the free energy '''
        '''
        入力値vは必ずVariable化
        本来行単位でSUMを取ってから列単位(バッチ数)でSUMを取るべき処理だけど
        結局SUMを取るので一括でSUM
        '''
        batch_size = v.data.shape[0]
        n_visible = self.n_visible
        real = self.real
        if real == 0:
            '''
            可視層が[0,1]のとき
            vbias_term = -1 * SUM((a(i) * v(i))
            '''
            vbias_term = F.sum(F.matmul(v,self.l.a))
        else:
            '''
            可視層が実数のとき
            vbias_term = -0.5 * SUM((v(i)-a(i)) * (v(i)-a(i)))
            chainerではバッチ数*入力層で扱うため、各行からbiasを引くため
            無理やり m＊nを使用
            '''
            m = Variable(xp.ones((batch_size,1),dtype=xp.float32))
            n = F.reshape(self.l.a,(1,n_visible))
            v_ = v - F.matmul(m,n)
            vbias_term = -F.sum(0.5 * v_ * v_)

        wx_b = self.l(v)
        hidden_term = F.sum(F.log(1+F.exp(wx_b)))
        return -vbias_term-hidden_term

    def propup(self,vis):
        '''
        可視層データから隠れ層の確率分布を計算
        隠れ層は[0,1]なので1となる確率を返す。
        入力値visは必ずxp(np)
        '''
        pre_sigmoid_activation = xp.dot(vis,self.l.W.data.T) + self.l.b.data
        return sigmoid(pre_sigmoid_activation)

    def propdown(self,hid):
        '''
        隠れ層データから可視層の確率分布を計算
        可視層が[0,1]の時は1となる確率を返す。
        可視層が実数の時は正規分布の平均を返す。分散は1固定。
        入力値hidは必ずxp(np)
        '''
        real = self.real
        if real == 0:
            pre_sigmoid_activation = xp.dot(hid,self.l.W.data) + self.l.a.data
            v_mean = sigmoid(pre_sigmoid_activation)
        else:
            v_mean = xp.dot(hid,self.l.W.data) + self.l.a.data
        return v_mean

    def sample_h_given_v(self,v0_sample):
        '''
        v0_sample → h1_sampleをgibbs sampling
        [0,1]の値をrandom作成して、popupで計算した1となる確率と比較
        h1_mean[i] >  p[i] なら h1_sample[i] = 1
        h1_mean[i] <= p[i] なら h1_sample[i] = 0
        '''
        h1_mean = self.propup(v0_sample)
        h1_sample = xp.random.binomial(size=h1_mean.shape,n=1,p=h1_mean)
        return h1_mean,h1_sample

    def sample_v_given_h(self,h0_sample):
        '''
        h0_sample → v1_sampleをgibbs sampling
        可視層が[0,1]の場合
        [0,1]の値をrandom値 p[i]を作成して、popdownで計算した1となる確率と比較
        v1_mean[i] >  p[i] なら v1_sample[i] = 1
        v1_mean[i] <= p[i] なら v1_sample[i] = 0
        可視層が実数の場合
        v1_sample[i]は平均v1_mean[i]、分散1の正規分布となるので
        平均0、分散1の正規分布からrandom値 u[i]を作成して、popdownで計算した値を追加
        v1_sample[i] = v1_mean[i] + u[i]
        '''
        v1_mean = self.propdown(h0_sample)
        if self.real == 0:
            v1_sample = xp.random.binomial(size=v1_mean.shape,n=1,p=v1_mean)
        else:
            batch_number = h0_sample.shape[0]
            v1_sample = v1_mean + xp.random.randn(batch_number,self.n_visible)
        return v1_mean,v1_sample

    def gibbs_hvh(self,h0_sample):
        ''' h->v->hをgibbs sampling '''
        v1_mean,v1_sample = self.sample_v_given_h(h0_sample)
        h1_mean,h1_sample = self.sample_h_given_v(v1_sample)
        return v1_mean,v1_sample,h1_mean,h1_sample

    def gibbs_vhv(self,v0_sample):
        ''' v->h->vをgibbs sampling '''
        h1_mean,h1_sample = self.sample_h_given_v(v0_sample)
        v1_mean,v1_sample = self.sample_v_given_h(h1_sample)
        return h1_mean,h1_sample,v1_mean,v1_sample

    def constrastive_divergence(self,v0_sample,k=1):
        ''' CD-kの処理、GPU化させるにはcupy必要 '''
        vh_sample = v0_sample
        for step in range(k):
            ph_mean,ph_sample,vh_mean,vh_sample = self.gibbs_vhv(vh_sample)
        return vh_sample

    def reconstruct(self, v):
        h = sigmoid(xp.dot(v,self.l.W.data.T) + self.l.b.data)
        reconstructed_v = sigmoid(xp.dot(h,self.l.W.data) + self.l.a.data)
        return reconstructed_v

chainerの場合、chainer.links.Linearを使って各層の写像の重み行列(W),バイアス(b)を定義できますが
RBMの場合追加でバイアスが必要となるため
l.add_param("a",(n_visible),dtype=xp.float32)
でパラメータを追加しました。
l.W:重み行列
l.b:hbias
l.a:vbaias
となります。
可視層から隠れ層への写像をchainer.links.Linearで管理させるため
一般的なRBMの本に記載されている重み行列を転置しています。(Wij→Wji)

1. MNISTで実験

以下のソースコードをMNISTのデータで実験してみました。
GPUのサーバがなかったのでCPUで学習させたため、学習回数は50回です。
今回は隠れ層を500次元とて、学習係数0.01,モメンタム0.5,PCD-10での結果を使用しています。

以下がoriginalのデータ。実際は60000件あるが、150件のみ表示

以下が上記の画像をreconstructした画像。ほぼ再現ができています。

件のみ表示。

そして次が重み行列を可視化した画像である。隠れ層は500次元で各ベクトルの基底に関して
28 * 28の画像に変換したものです。