【theano-windows】学习笔记十五——受限玻尔兹曼机,p6 root

文件名：【theano-windows】学习笔记十五——受限玻尔兹曼机,p6 root 【theano-windows】学习笔记十五——受限玻尔兹曼机 前言

终于到了最喜欢的模型: 受限玻尔兹曼机(RBM)了, 发现关于RBM是如何从能量模型发展过来的介绍非常不错, 而关于详细理论证明, 可以去看我前面的受限玻尔兹曼机的一系列博客.

国际惯例, 参考博客,超级强推第二个博客, 证明过程很给力:

Restricted Boltzmann Machines (RBM)

Contrastive divergence multi-layer RBMs

理论 基于能量的模型EBM

这个强烈去看看我前面的两篇关于概率有向图和概率无向图模型的博客, 很简短. 基于能量的模型是将所有的感兴趣变量与一个标量能量关联起来. 学习对应着修改能量函数, 而使其具有应有的结构. 我们也经常希望这个模型具有很低的能量. 至于为什么? 我们经常举的一个例子, 将一个小球放到碗里, 这个球最终基本都会停留在碗底, 为什么？当然是势能最低啦，这就是我们说的最小化能量函数。

基于能量的概率模型通过能量函数定义了一个概率分布:

p(x)=e−E(x)Z 其中分母是归一化因子,与概率模型中的归一化因子意义相同 Z=∑xe−E(x) 优化过程就是针对训练集上的负对数似然进行随机梯度下降法训练, 和logistics回归一样,我们先定义一个对数似然, 然后损失函数就是负对数似然 L(θ,D)=1N∑x(i)∈Dlogp(x(i))loss=−L(θ,D) 其中θ是模型参数, D是训练集，最终更新方法依旧是梯度下降,针对模型参数求解−∂logp(x(i))∂θ

具有隐层的能量模型

你以为所有的模型都有隐层么？想想马尔科夫模型是怎么发展到隐马尔科夫模型的吧. 千万不要有这个错误的想法:所有的模型都有隐层。

上面说过能量模型是定义一个概率分布, 与所有的感兴趣变量相关, 当我们加入了除了输入单元x之外的隐单元后, 整个模型的关于输入单元的概率分布就变成了关于所有隐单元的边缘分布了 p(x)=∑hp(x,h)=∑he−E(x,h)Z 为了看起来与能量模型中的第一个式子看起来像是双胞胎, 引入了一个自由能函数

F(x)=−log∑h(e−E(x,h)) 这样p(x)就可以写成: p(x)=e−F(x)Z 其中Z=∑xe−F(x)

然后我们求一波导数:

−∂logp(x)∂θ=−∂(loge−F(x)−logZ)∂θ=∂(F(x)+logZ)∂θ=∂F(x)∂θ−e−F(x^)Z∗∂F(x^)∂θ=∂F(x)∂θ−∑x^p(x^)∂F(x^)∂θ 注意这里为什么第一项是x, 而第二项是x^ ？可以看上面第二个博客的证明, 也可以去看看我前面关于概率无向图模型的简介, 它是输入样本所有可能的情况. 这两项经常被称为positive phase 和negative phase, 我经常称作正向推断和反向重构, 它们反应的是对模型概率密度的影响. 第一项通过降低对应自由能量去增加训练集的概率, 而第二项降低了模型生成样本的概率.

由于第二项的存在, 我们很难直接去计算梯度, 它是关于输入的所有可能组合情况的期望, 第二项的简写表达式为Ep(∂F(x)∂θ), 这个其实就是期望的另外种表示. 既然无法知道所有的情况, 就可以利用采样的方法, 选取固定数量的模型样本估计这个期望值, 用于估计负相位梯度的样本称为负颗粒(negative particles),记为ℵ, 梯度的新写法就是:

−∂logp(x)∂θ≈∂F(x)∂θ−1ℵ∑x^∈ℵ∂F(x^)∂θ 而这个固定数量的样本是如何选取的？一般来说用的是蒙特卡洛方法. 我前面有博客介绍过这个采样方法. 后面我们说的吉布斯采样也是属于马尔科夫链蒙特卡洛方法的一种.

玻尔兹曼机

随机变量X(观测变量)和隐变量H, 它们的联合分布是与能量函数有关的玻尔兹曼分布:

P(X=x,H=h)=e−E(x,h)Z 其中归一化常量是: Z=∑x,he−E(x,h) 在一般的玻尔兹曼机中,能量函数是一个二次多项式(quadratic polynomial),因为模型中隐层和可见层之间, 以及它们各自的单元都互相连接, 说明它们被同等对待了, 那么假设z=(x,h),那么,能量函数为 E(z)=E(x,h)=−∑ibizi−∑ijziwijzjs.t.zi∈{0,1} 如果X=x是观测的, H是隐藏的, 那么对数似然就是所有可能的H组合的和 P(X=x)=∑hP(X=x,H=h)=∑he−E(x,h)Z 在非限制玻尔兹曼机中分子中的h和分母中的x的求和是无法解决的问题, 因而它们的梯度也无法求得: ∂−logP(x)∂θ=∂−log∑he−E(x,h)Z∂θ=∂(−log∑he−E(x,h)+logZ)∂θ=∂−log∑he−E(x,h)∂θ+∂logZ∂θ 然后我们分别求第一项和第二项 ∂−log∑he−E(x,h)∂θ=−(1∑he−E(x,h)⋅∂∑he−E(x,h)∂θ)=−(1∑he−E(x,h)⋅∑h(−e−E(x,h)⋅∂E(x,h)∂θ))=−∑h(−e−E(x,h)∑he−E(x,h)⋅∂E(x,h)∂θ)=∑h(P(x,h)P(x)⋅∂E(x,h)∂θ)=∑h(P(h|x)⋅∂E(x,h)∂θ) 同理求第二项 ∂logZ∂θ=1Z⋅∂Z∂θ=1Z⋅∂∑x,he−E(x,h)∂θ=1Z⋅∑x,h(e−E(x,h)⋅∂E(x,h)∂θ)=∑x,h(e−E(x,h)Z⋅∂E(x,h)∂θ)=∑x,h(p(x,h)⋅∂E(x,h)∂θ) 所以最终的 ∂−logP(x)∂θ=∑h(P(h|x)⋅∂E(x,h)∂θ)−∑x,h(p(x,h)⋅∂E(x,h)∂θ) 玻尔兹曼机是对数线性马尔科夫随机场的特定形式, 参数在能量中是线性的, 为了足以强大以表示复杂分布, 为其添加了隐变量, 具有更多隐单元使得玻尔兹曼机的建模能力更强。

限制玻尔兹曼机 能量相关内容

前面那个玻尔兹曼机是隐层和可见层之间和各自内部的单元之间都有连接, 也就是说所有的神经元之间都有连接, 那么如果让隐单元自己内部和可见层自己内部的单元连接权重为0, 其实也就是没连接的时候, 就成了限制玻尔兹曼机了. 它的好处在于在可见层单元X给定的情况下,各个隐单元Hi是相互独立的,反之也成立, 给定隐单元的时候, 可见层单元也相互独立。

给定输入vi和二项隐单元(binomial unit)hj,参数是(bj,Wj⋅),分别代表隐单元的偏置和连接到隐单元j的权重向量, 那么就可以得到:

能量函数:−bj⋅hj−∑iviWijhj P(hj=1|v)=ebj+∑iWijvi1+ebj+∑iWijvi=sigmoid(bj+∑iWijvi)

如果给定的是固定方差的高斯单元hj, 参数为(aj,bj,Wj⋅), 那么也可以得到

能量函数: a2jh2j−bjhj−∑iviWijhjP(hj|v)=1Ze−a2jh2j+bjhj+∑iviWijhj=1Ze−(bj−μ)22σ2 P(hj|v)=N(hj;μ,σ2) with σ2=12a2j,μ=bj+∑iviWij2a2j

当a特别小的时候, 会容易崩掉(blow up), 因而经常使用a+ϵ代替a

对于softmax输出单元hj,参数为(Wij,bj):

能量函数: −bjhj−∑iviWijhjP(hj=1|v)=ebj+∑iWijvi∑j′ebj′+∑iWij′vi=softmax(bj+∑iWijvi)

参数更新方法如下:

假设归一化常量为Z=∑x,ye−E(x,y), 那么

P(x,y)P(y|x)=e−E(x,y)Z=e−E(x,y)∑ye−E(x,y) 对于任意的能量函数,我们都可以推导出: ∂−logP(x)∂θ=∂−log∑yP(x,y)∂θ=∂−log∑ye−E(x,y)Z∂θ=−Z∑ye−E(x,y)∂∑ye−E(x,y)Z∂θ=−Z∑ye−E(x,y)∑y(−e−E(x,y)∂E(x,y)∂θ⋅Z)−∑ye−E(x,y)⋅∂Z∂θZ2=∑y(e−E(x,y)∑ye−E(x,y)∂E(x,y)∂θ)+1Z∂Z∂θ=∑yP(y|x)∂E(x,y)∂θ−∑xyP(x,y)∂E(x,y)∂θ=E[∂E(x,y)∂θ∣x]−E[∂E(x,y)∂θ] 在使用对比散度算法做梯度更新的时候:

二值输入与二值输出的时候:

权重更新 正向阶段贡献: P(y0j=1|x0)×1×x0i+(1−P(y0j=1|x0))×0×x0i=P(y0j=1|x0)x0i反向阶段共享:P(y1j=1|x1)×1×x1i+(1−P(y1j=1|x1))×0×x1i=P(y1j=1|x1i)x1i偏置更新 正向阶段贡献:P(y0i=1|x0)反向阶段贡献:P(y1i=1|x1)

高斯输入与二值输出的时候:

权重和偏置更新同上参数aj: 正向阶段贡献:2ai(x0i)2反向阶段贡献:2ai(x1i)2

二值输入和softmax输出的时候:

与二值情况相同, 只不过P(yi=1|x)使用的是softmax而非sigmoid 理论梳理

前面可能说的有点乱, 下面直接按照教程梳理一下玻尔兹曼机的参数更新过程. 假设:

可见层V={v0,v1,⋯,vj,⋯},偏置b={b1,b2,⋯,bj,⋯}隐藏层H={h1,h2,⋯,hi,⋯}, 偏置c=c1,c2,⋯,ci,⋯连接权重为W={w11,w12,⋯,wji,⋯}

那么能量函数为:

E(v,h)=−b′v−c′h−h′Wv 自由能函数为 F(v)=−log∑hie−E(v,h)=−b′v−∑ilog∑hiehi(ci+wiv) 特别地, 在二值RBM中, 经常通过以下激活函数得到神经元的激活概率: P(hi=1|v)=sigmoid(ci+Wiv)P(vj=1|h)=sigmoid(bj+W′jh) 然后自由能函数简化为: F(v)=−b′v−∑ilog∑hiehi(ci+wiv)=−b′v−∑ilog{e0∗(ci+wiv)+e1∗(ci+wiv)}=−b′v−∑ilog(1+e(ci+wiv)) 梯度更新方法为: −∂logP(v)∂Wij−∂logP(v)∂bj−∂logP(v)∂ci=Ev[p(hi|v)⋅vj]−v(i)j⋅sigmoid(Wi⋅v(i)+ci)=Ev[p(vj|h)−v(i)j]=Ev[p(hi|v)]−sigmoid(Wi⋅v(i))

判别式玻尔兹曼机

当玻尔兹曼机用于分类的时候,假设:

输入单元为P,索引j隐单元为L,索引i,偏置为C标签单元为Y,索引k, 偏置为B输入单元和隐单元的连接权重为V输出单元和隐单元的连接权重为W

那么输入单元与隐单元的能量函数为−VijPjLi−CiLi,

输出单元和隐单元的能量函数为−WikYkLi−BkYk

标签y的计算与输入P相关, 而非隐单元L active(Yk)=Bk+∑i⎛⎝softplus(Wki+Ci+∑j(Vij⋅p(Pj))⎞⎠ 其中softplus(x)=log(1+ex)

那么输出

P(Y∣inputs)=softmax(active(Y)) 训练过程:

第一种方法是直接梯度下降:

计算cost=−logP(Y=onehot(k)∣inputs)

第二种方法是对比散度(contrastive divergence):

将X=[Y,P]当做一个大的层让x=[onehot(k),p(P)]作为隐层的输入将(X,L)当做RBM训练

关于这一部分可以看前面的一个博客判别模型的玻尔兹曼机论文源码解读

代码 训练

【注】此代码是为了更加方便理解而对官方教程代码进行了理解后的修改, 建议看完以后立即去看标准的官方教程

首先引入必要的库文件

import numpy as npimport theanoimport theano.tensor as Tfrom theano.tensor.shared_randomstreams import RandomStreamsimport cPickle,gzipfrom PIL import Imageimport pylabimport os

随后是读数据集的函数

#定义读数据的函数,把数据丢入到共享区域def load_data(dataset):data_dir,data_file=os.path.split(dataset)if os.path.isfile(dataset):with gzip.open(dataset,'rb') as f:train_set,valid_set,test_set=cPickle.load(f)#共享数据集def shared_dataset(data_xy,borrow=True):data_x,data_y=data_xyshared_x=theano.shared(np.asarray(data_x,dtype=theano.config.floatX),borrow=borrow)shared_y=theano.shared(np.asarray(data_y,dtype=theano.config.floatX),borrow=borrow)return shared_x,T.cast(shared_y,'int32')#定义三个元组分别存储训练集,验证集,测试集train_set_x,train_set_y=shared_dataset(train_set)valid_set_x,valid_set_y=shared_dataset(valid_set)test_set_x,test_set_y=shared_dataset(test_set)rval=[(train_set_x,train_set_y),(valid_set_x,valid_set_y),(test_set_x,test_set_y)]return rval

接下来就是定义整个RBM模型:

首先是初始化该有的权重和偏置,一个权重两个偏置, numpy的随机函数用于初始化权重和偏置, theano的随机函数用于将神经元按照概率激活为二值单元

#定义RBMclass RBM(object):def __init__(self,rng=None,trng=None,input=None,n_visible=784,n_hidden=500,W=None,hbias=None,vbias=None):self.n_visible=n_visibleself.n_hidden=n_hiddenif rng is None:rng=np.random.RandomState(1234)if trng is None:trng=RandomStreams(rng.randint(2**30))#初始化权重和偏置 if W is None:initW=np.asarray(rng.uniform(low=-4*np.sqrt(6./(n_hidden+n_visible)),high=4*np.sqrt(6./(n_hidden+n_visible)),size=(n_visible,n_hidden)),dtype=theano.config.floatX)W=theano.shared(value=initW,name='W',borrow=True)if hbias is None:inithbias=np.zeros(n_hidden,dtype=theano.config.floatX)hbias=theano.shared(value=inithbias,name='hbias',borrow=True)if vbias is None:initvbias=np.zeros(n_visible,dtype=theano.config.floatX)vbias=theano.shared(value=initvbias,name='vbias',borrow=True)self.input=inputif not input:self.input=T.matrix('input')self.W=Wself.hbias=hbiasself.vbias=vbiasself.trng=trngself.params=[self.W,self.hbias,self.vbias]

为了进行吉布斯采样, 我们需要定义positive phase和negative phase 前向计算过程: 可见层->隐层

##########前向计算,从可见层到隐层#################激活概率def propup(self,vis):pre_sigmoid_activation=T.dot(vis,self.W)+self.hbiasreturn [pre_sigmoid_activation,T.nnet.sigmoid(pre_sigmoid_activation)]#二值激活def sample_h_given_v(self,v0_samples):pre_sigmoid_h1,h1_mean=self.propup(v0_samples)h1_sample=self.trng.binomial(size=h1_mean.shape,n=1,p=h1_mean,dtype=theano.config.floatX)return [pre_sigmoid_h1,h1_mean,h1_sample]

反向计算:隐层->可见层

##########反向计算,从隐层到可见层#################激活概率def propdown(self,hid):pre_sigmoid_activation=T.dot(hid,self.W.T)+self.vbiasreturn [pre_sigmoid_activation,T.nnet.sigmoid(pre_sigmoid_activation)]#二值激活def sample_v_given_h(self,h0_samples):pre_sigmoid_v1,v1_mean=self.propdown(h0_samples)v1_sample=self.trng.binomial(size=v1_mean.shape,n=1,p=v1_mean,dtype=theano.config.floatX)return [pre_sigmoid_v1,v1_mean,v1_sample]

一次吉布斯采样: 可见层->隐层->可见层

##########吉布斯采样#################可见层->隐层->可见层def gibbs_vhv(self,v0_samples):pre_sigmoid_h1,h1_mean,h1_sample=self.sample_h_given_v(v0_samples)pre_sigmoid_v1,v1_mean,v1_sample=self.sample_v_given_h(h1_sample)return [pre_sigmoid_v1,v1_mean,v1_sample,pre_sigmoid_h1,h1_mean,h1_sample]

从第二章节: 具有隐层的能量模型来看, 梯度更新就是对自由能量函数球梯度, 当然如果我们自己写梯度的话,可以按照RBM理论梳理中来做, 事实上在matlab就是这样玩的, 手动写梯度更新方法, 而theano提供了自动求导方法, 所以我们就直接定义自由能函数, 然后对其求导即可 自由能量函数的定义:

############自由能量函数###############def free_energy(self,v_samples):wx_b=T.dot(v_samples,self.W)+self.hbiasvbias_term=T.dot(v_samples,self.vbias)#第一项hidden_term=T.sum(T.log(1+T.exp(wx_b)),axis=1)#第二项return -hidden_term-vbias_term

然后是梯度更新

############梯度更新#################def get_cost_updates(self,lr=0.1,k=1):([pre_sigmoid_nvs,nv_means,nv_samples,pre_sigmoid_nhs,nh_means,nh_samples],updates)=\theano.scan(self.gibbs_vhv,outputs_info=[None,None,self.input,None,None,None],n_steps=k,name='gibbs_vhv')chain_end=nv_samples[-1]cost=T.mean(self.free_energy(self.input))-T.mean(self.free_energy(chain_end))gparams=T.grad(cost,self.params,consider_constant=[chain_end])for gparam,param in zip(gparams,self.params):updates[param]=param-gparam*T.cast(lr,dtype=theano.config.floatX)##################期望看到交叉熵损失##############monitor_cost=self.get_reconstruction_cost(pre_sigmoid_nvs[-1])return monitor_cost,updates

当然我们看到了一个函数叫get_reconstruction_cost, 想想RBM是干什么? 是模型的估计出来的样本分布和真实样本分布更加接近, 而度量这个的最好方法就是交叉熵了, 这就是我们在更新的时候希望看到的类似于loss的东东:

########非持续性对比散度,重构误差#########def get_reconstruction_cost(self,pre_sigmoid_nv):cross_entropy=T.mean(T.sum(self.input*T.log(T.nnet.sigmoid(pre_sigmoid_nv))+\(1-self.input)*T.log(1-T.nnet.sigmoid(pre_sigmoid_nv)),axis=1))return cross_entropy

接下来我们就可以训练了, 一下采用15次CD算法训练(虽然一次就已经足够了)

#初始化并训练玻尔兹曼机def test_rbm(learning_rate=0.1,training_epochs=15,dataset='mnist.pkl.gz',batch_size=20,n_chains=20,n_samples=10,n_hidden=500):#读取数据集datasets=load_data(dataset)train_set_x,train_set_y=datasets[0]test_set_x,test_set_y=datasets[2]n_train_batches=train_set_x.get_value(borrow=True).shape[0]//batch_sizeindex=T.lscalar()x=T.matrix('x')rng=np.random.RandomState(1234)trng=RandomStreams(rng.randint(2**30))#初始化一个RBM实例rbm=RBM(input=x,n_visible=28*28,n_hidden=n_hidden,rng=rng,trng=trng)#更新参数cost,updates=rbm.get_cost_updates(lr=learning_rate,k=15)#训练RBMtrain_rbm=theano.function([index],cost,updates=updates,givens={x:train_set_x[index*batch_size:(index+1)*batch_size]},name='train_rbm')for epoch in range(training_epochs):mean_cost=[]for batch_index in range(n_train_batches):mean_cost+=[train_rbm(batch_index)]print ('Training epoch %d,cost is ' % epoch,np.mean(mean_cost))save_file=open('best_model_rbm.pkl','wb')model=rbm.paramsl=[model[0].get_value(),model[1].get_value(),model[2].get_value()]cPickle.dump(l,save_file)

接下来看看训练及结果

test_rbm()'''('Training epoch 0,cost is ', -221.54436)('Training epoch 1,cost is ', -193.33549)('Training epoch 2,cost is ', -188.36243)('Training epoch 3,cost is ', -185.52127)('Training epoch 4,cost is ', -183.83319)('Training epoch 5,cost is ', -182.85869)('Training epoch 6,cost is ', -181.76605)('Training epoch 7,cost is ', -180.4059)('Training epoch 8,cost is ', -179.306)('Training epoch 9,cost is ', -178.59416)('Training epoch 10,cost is ', -178.47371)('Training epoch 11,cost is ', -177.25623)('Training epoch 12,cost is ', -177.09161)('Training epoch 13,cost is ', -176.70366)('Training epoch 14,cost is ', -175.92307)''' 测试

先读一张图片, 还是使用之前教程用过的代码

#取一张图片预测from PIL import Imageimport pylabdataset='mnist.pkl.gz'datasets=load_data(dataset)test_set_x,test_set_y=datasets[2]test_set_x=test_set_x.get_value()test_data=test_set_x[12:13]img=test_data.reshape(28,28)pylab.imshow(img)pylab.show()

初始化一个RBM, 并对参数赋值

a=cPickle.load(open('best_model_rbm.pkl'))x=T.matrix('x')rbm_test=RBM(input=x,n_visible=28*28,n_hidden=500)rbm_test.params[0].set_value(a[0])rbm_test.params[1].set_value(a[1])rbm_test.params[2].set_value(a[2])

执行1000次的吉布斯采样

([pre_sigmoid_nvs,nv_means,nv_samples,pre_sigmoid_nhs,nh_means,nh_samples],updates)=\theano.scan(rbm_test.gibbs_vhv,outputs_info=[None,None,x,None,None,None],n_steps=1000,name='gibbs_vhv')recon=nv_samples[-1]recon_fn=theano.function([x],recon,updates=updates)

可视化看看结果

b=recon_fn(test_data)b=b.reshape(28,28)pylab.imshow(b)pylab.show()

【注】若感觉程序有问题, 请反馈在评论区, 毕竟我是python渣 code链接:链接: https://pan.baidu.com/s/1hsiYqNi 密码: wu2d