前言
终于到了最喜欢的模型: 受限玻尔兹曼机(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 是训练集，最终更新方法依旧是梯度下降,针对模型参数求解
具有隐层的能量模型
你以为所有的模型都有隐层么？想想马尔科夫模型是怎么发展到隐马尔科夫模型的吧. 千万不要有这个错误的想法:所有的模型都有隐层。
上面说过能量模型是定义一个概率分布, 与所有的感兴趣变量相关, 当我们加入了除了输入单元
x
之外的隐单元后, 整个模型的关于输入单元的概率分布就变成了关于所有隐单元的边缘分布了
为了看起来与能量模型中的第一个式子看起来像是双胞胎, 引入了一个自由能函数
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 , 而第二项是
由于第二项的存在, 我们很难直接去计算梯度, 它是关于输入的所有可能组合情况的期望, 第二项的简写表达式为
Ep(∂F(x)∂θ)
, 这个其实就是期望的另外种表示. 既然无法知道所有的情况, 就可以利用采样的方法, 选取固定数量的模型样本估计这个期望值, 用于估计负相位梯度的样本称为负颗粒(negative particles),记为
ℵ
, 梯度的新写法就是:
−∂logp(x)∂θ≈∂F(x)∂θ−1ℵ∑x^∈ℵ∂F(x^)∂θ
而这个固定数量的样本是如何选取的？一般来说用的是蒙特卡洛方法. 我前面有博客介绍过这个采样方法. 后面我们说的吉布斯采样也是属于马尔科夫链蒙特卡洛方法的一种.
玻尔兹曼机
随机变量
X
(观测变量)和隐变量
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 是隐藏的, 那么对数似然就是所有可能的
P(X=x)=∑hP(X=x,H=h)=∑he−E(x,h)Z
在非限制玻尔兹曼机中分子中的 h 和分母中的
∂−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
, 其实也就是没连接的时候, 就成了限制玻尔兹曼机了. 它的好处在于在可见层单元
给定输入 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−∑iviWijhj
-
P(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−∑iviWijhj
- P(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 np import theano import theano.tensor as T from theano.tensor.shared_randomstreams import RandomStreams import cPickle,gzip from PIL import Image import pylab import 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_xy shared_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的随机函数用于将神经元按照概率激活为二值单元
#定义RBM class 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_visible self.n_hidden=n_hidden if 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=input if not input: self.input=T.matrix('input') self.W=W self.hbias=hbias self.vbias=vbias self.trng=trng self.params=[self.W,self.hbias,self.vbias]为了进行吉布斯采样, 我们需要定义positive phase和negative phase
前向计算过程: 可见层->隐层##########前向计算,从可见层到隐层################ #激活概率 def propup(self,vis): pre_sigmoid_activation=T.dot(vis,self.W)+self.hbias return [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.vbias return [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.hbias vbias_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_size index=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) #训练RBM train_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.params l=[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 Image import pylab dataset='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
扫一扫
1268
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
