1. 前言

虽然推导过二值形式的RBM，但是对于可见层为实值的输入还是半知半解的。最近写个深度学习相关综述，看了一些关于RBM的文献，这里做一下对比总结。但是不倾向于对实值RBM的推导，而是相关代码的实现。

2. RBM回顾

RBM是具有特殊结构的玻尔兹曼机，可见单元到隐单元是全连接，但是层内是无连接的，也就是可见单元不会和可见单元连一起。RBM的这种二分结构保证了当给定可见单元的时候，隐藏层个单元是相互独立的，至于为什么，可以看我前面对概率图模型截图的两篇文章：概率无向图模型、概率有向图模型，这两篇文章不长，而且对理解RBM非常有帮助，毕竟RBM也属于概率无向能量图模型。
【笔者注】有句话要记住：“没有万能的模型，只有不断优化的理论和结构”。个人非常不建议这样的问题：“为什么我把XX结构用到XX数据上效果不好哇”，“XX这个结构用到XX上出现这个问题是什么情况”。多看论文，不都迎刃而解了么。

3. 二值RBM

3.1 参数更新理论

定义RBM的可见单元为$v_i$，对应偏置为$a_i$，隐层单元为$h_j$，对应偏置为$b_j$，那么RBM的可见层与隐藏层的联合概率分布就是
$$P(\mathbf{v},\mathbf{h})=\frac{exp(-E(\mathbf{v},\mathbf{h}))}{Z}$$
其中$E(\mathbf{v},\mathbf{h})$称为能量函数。当可见单元和隐单元都是二值的时候(值为0或者1)，能量函数定义为
$$E(\mathbf{v},\mathbf{h})=-\sum _{ij}{W}_{ij}{v}_i{h}_j-\sum _i{a}_i{v}_i-\sum _j{b}_j{h}_j$$
联合概率分布中的分母$Z$为配分函数(partition function)，用于归一化，可以确保$P(\mathbf{v},\mathbf{h})$是一个有效的概率分布${}^1$,其实也就是保证了和为1。这一项很难计算，因为存在对(指数级)连接方式的求和${}^2$
$$Z=\sum _{{\mathbf{v}}^{\prime },{\mathbf{h}}^{\prime }}E({\mathbf{v}}^{\prime },{\mathbf{h}}^{\prime })$$
随后对联合概率分布求解边缘分布，并且求解极大似然函数，便能够得到权重的更新方法
$$\Delta W_{ij}\propto <v_i{h}_j{>}_{data}-<v_i{h}_j{>}_{model}$$
其中$<\cdot {>}_{data}$代表对数据分布的期望，$<\cdot {>}_{model}$代表RBM能量模型的平衡分布 ，通过将训练集中的一个向量输入到可见单元，可以得到$<v_i{h}_j>$的一个无偏样本，可以并行计算隐单元
$$P(h_j=1|v)=\frac{1}{1+exp(-b_j-{\sum }_i{W}_{ij}{v}_i)}$$
对于从整个训练数据中的提取的’小批量’数据的每一个向量都要使用这个式子，主要是为了得到$<v_i{h}_j{>}_{data}$的经验估计(empirical estimate)。为了计算$<v_i{h}_j>$，我们需要从联合概率分布$p(\mathbf{v},\mathbf{h})$中获取无偏样本，方法就是吉布斯采样，通过迭代执行$P(\mathbf{h}|\mathbf{v})$采样以及$P(\mathbf{v}|\mathbf{h})$
$$P(v_i=1|\mathbf{h})=\frac{1}{1+exp(-a_i-{\sum }_j{W}_{ij}{h}_j)}$$
但是，在高维空间中执行吉布斯采样，需要很久的时间才能收敛。所以使用下式代替上面更新权重的式子
$$\Delta W_{ij}\propto <v_i{h}_j{>}_{data}-<v_i{h}_j{>}_{recon}$$
式中的第二项对应"重构"数据的分布。重构数据是通过从可见层上的一个数据向量开始，不断使用$P(h_j=1|\mathbf{v})$对隐层采样和使用$P(v_i=1|\mathbf{h})$对可见层采样，交错采样K次即可，偏置的学习方法为
$$\begin{array}{rl}\Delta a_i& \propto <v_i{>}_{data}-<v_j{>}_{recon}\\ \Delta b_j& \propto <h_j{>}_{data}-<h_j{>}_{recon}\end{array}$$
这个步骤不是极大似然学习法，因为极大似然学习很慢；如果我们近似地沿着另一个函数(对比散度)的梯度，学习效果依旧很好${}^2$${}^5$，其实就是交错进行K次完整步骤的吉布斯采样。

用向量的方法写RBM的联合概率分布，迭代采样中的条件概率分布如下${}^1$：
$$\begin{array}{rl}& P(\mathbf{v},\mathbf{h})=exp(v^T{b}_v+h^T{b}_h+v^{T}Wh)/Z\\ & P(H^{(j)}=1|v)=s(b_H+W^{T}v{)}^{(j)}\\ & P(V^{(i)}=1|h)=s(b_V+Wh{)}^{(i)}\end{array}$$
其中$s(x{)}^{(j)}=(1+exp(-x^{(j)}){)}^{-1}$是logistic函数，$x^{(j)}$是向量$x$的第$j$个成分。

3.2 参数更新代码

找代码的目的主要是为了看理论和代码的一致性，Hinton提供的代码结构分析如下：

1. 可见层$\to$ 隐藏层

   利用可见层对隐层采样

$$P(h_j=1|v)=\frac{1}{1+exp(-b_j-{\sum }_i{W}_{ij}{v}_i)}$$

  data = batchdata(:,:,batch);poshidprobs = 1./(1 + exp(-data*vishid - repmat(hidbiases,numcases,1)));  batchposhidprobs(:,:,batch)=poshidprobs;

​ 当然，由于隐层是二值形式，所以采用随机阈值激活方法

poshidstates = poshidprobs > rand(numcases,numhid);

​ 记录$<v_i{h}_j{>}_{data};<v_i{>}_{data};<h_j{>}_{data}$

 posprods    = data' * poshidprobs;poshidact   = sum(poshidprobs);posvisact = sum(data);

2. 隐藏层$\to$可见层

   利用隐层二值状态poshidstates对可见层重构
   $$P(v_i=1|\mathbf{h})=\frac{1}{1+exp(-a_i-{\sum }_j{W}_{ij}{h}_j)}$$

   negdata = 1./(1 + exp(-poshidstates*vishid' - repmat(visbiases,numcases,1)));
   

   记录$<v_i{h}_j{>}_{recon};<v_i{>}_{recon};<h_j{>}_{recon}$

   neghidprobs = 1./(1 + exp(-negdata*vishid - repmat(hidbiases,numcases,1)));  negprods  = negdata'*neghidprobs;neghidact = sum(neghidprobs);negvisact = sum(negdata); 
   

3. 更新模型参数

$$\begin{array}{rl}\Delta W_{ij}& \propto <v_i{h}_j{>}_{data}-<v_i{h}_j{>}_{recon}\\ \Delta a_i& \propto <v_i{>}_{data}-<v_j{>}_{recon}\\ \Delta b_j& \propto <h_j{>}_{data}-<h_j{>}_{recon}\end{array}$$

 vishidinc = momentum*vishidinc + ...epsilonw*( (posprods-negprods)/numcases - weightcost*vishid);visbiasinc = momentum*visbiasinc + (epsilonvb/numcases)*(posvisact-negvisact);hidbiasinc = momentum*hidbiasinc + (epsilonhb/numcases)*(poshidact-neghidact);vishid = vishid + vishidinc;visbiases = visbiases + visbiasinc;hidbiases = hidbiases + hidbiasinc;

代码中多出的项分别为动量项momentum，学习率epsilon，权重衰减weightcost

由此可见，上述理论和代码可以作为二值RBM的标准写法。

4. 实值RBM

很多情况RBM的输入向量是实值形式的，比如手写数字的处理中，先对数据进行归一化，然后使用二值可见单元的实值概率取代激活情况。注意下面讨论的是可见单元为实值，隐单元为二值的RBM的情况，关于实值-实值的情况以后再讨论，其实有人说，实值-实值的RBM已经退化成PCA了，后续继续研究。

4.1 参数更新理论

4.1.1 实值RBM的能量函数

在Hinton的课程和文献${}^2$中书写方式为
$$E(\mathbf{v},\mathbf{h})=\sum _i\frac{(v_i-a_i{)}^2}{2{\sigma }_i^2}-\sum _{ij}(W_{ij}\cdot \frac{v_i}{{\sigma }_i}\cdot h_j)-\sum _j{b}_j{h}_j$$
在文献${}^4$中带范围的书写方式为
$$E(\mathbf{v},\mathbf{h})=\sum _{i=1}^m\frac{(v_i-a_i{)}^2}{2{\sigma }^2}-\sum _{j=1}^n\sum _{i=1}^m(W_{ij}\cdot \frac{v_i}{\sigma }\cdot h_j)-\sum _{j=1}^n{b}_j{h}_j$$
在文献${}^6$中使用平方的写法为
ERBM(v,h∣W,bv,bh,{σi2})=∑i(biv−vi)22σi2−∑i,j1σi2Wijvihj−∑jbjhhjE_{RBM}(\mathbf{v,h|W,b^v,b^h,\{\sigma_i^2}\})=\sum_i \frac{(b_i^v-v_i)^2}{2\sigma_i^2}-\sum_{i,j}\frac{1}{\sigma_i^2}W_{ij}v_ih_j-\sum_jb_j^hh_j ERBM​(v,h∣W,bv,bh,{σi2​})=i∑​2σi2​(biv​−vi​)2​−i,j∑​σi2​1​Wij​vi​hj​−j∑​bjh​hj​
其中$a_i$是可见单元的偏置，$b_j$是隐单元的偏置，$\sigma$是可见单元的高斯噪声的标准差${}^2$，对称权重$W_{ij}$连接可见曾和隐藏层。在Hinton在Coursera上的课程中提到了这个方差使得RBM的采样出现了问题

所以在实际中，经常让$\sigma =1$，而且在文献${}^5$中也提到了: 实际中，我们经常把我们的数据缩放到零均值和单位方差，固定式中方差${\sigma }_i$为1能够学习地更好，即使我们期望模型能够达到更高的精度，因而这就出现了后来的

• 在文献${}^1$中写出了向量表示的形式

$$E(\mathbf{v},\mathbf{h})=\frac{||v|{|}^2}{2}-v^T{b}_v-h^T{b}_h-v^{T}Wh$$

• 在文献${}^3$中的提出的能量函数写法

$$E(\mathbf{v},\mathbf{h})=\frac{1}{2}\sum _i{v}_i^2-\sum _{i,j}{v}_i{h}_j{w}_{ij}-\sum v_i{b}_i-\sum _j{h}_j{c}_j$$

4.1.2 实值RBM的采样

采样方法依旧是使用k次吉布斯采样，依据是一层的一个单元关于另一层所有单元的条件概率分布函数。需要注意的是，由于我们采用的是实值RBM，可见层单元服从零均值、单位方差的分布，所以依据隐层采样可见层不再是$P(v_i=1|\mathbf{h})$，而是利用隐单元对可见单元的均值做一个线性贡献${}^2$${}^6$，在高斯分布中进行采样得到的，而高斯分布的参数也是由下式给出的
$$\begin{array}{rl}P(v_i|\mathbf{h})& =N(a_i+{\sigma }_i\sum _j{W}_{ij}{h}_j,{\sigma }_i^2)\\ P(h_j=1|\mathbf{v})& =\frac{1}{1+exp(-b_j-{\sum }_i{W}_{ij}\frac{v_i}{\sigma })}\end{array}$$
其中$N(\mu ,{\sigma }^2)$代表均值为$\mu$，方差为$\sigma$的正太分布(normal distribution)

因为方差$\sigma =1$，所以出现了下面两种写法

• 文献${}^3$中的，在假设可见单元具有$N(0,1)$分布下，第$i$个可见单元可以从下式采样

$$N(\sum _j{w}_{ij}{h}_j+b_i,1)$$

• 文献${}^5$中提到的，假设${\sigma }_i=1$，那么条件分布就是

$$\begin{array}{rl}p(h_j=1|\mathbf{v})& =f(b_j+\sum _i{v}_i{w}_{ij})\\ P(v_i|\mathbf{h})& =N(c_i+\sum _j{h}_j{w}_{ij},1)\end{array}$$

其中$f(\cdot )$是logistic函数(对数函数)，$N(\mu ,V)$是高斯分布。

权重和偏置的更新方法如下

4.2参数更新代码

保险起见，对于实值RBM，分析的代码应该不止一个，而且读者没必要具体分析代码的各个变量，只需要关注更新的大概方法即可。代码公布下载戳此处：代码1、代码2、代码3 .

4.2.1 positive阶段

关于实值RBM训练阶段的positive过程(可见层$\to$隐层)就直接贴了，和二值形式一样，先计算条件概率$P(h_j=1|\mathbf{v})$，然后使用随机阈值激活为01二值状态。

• 代码1来自${}^5$ ，虽然代码有时序关系存在，多了几项，但是依旧可以看出端倪

eta =  w*(data(:,:,1)./gsd)' + ...   %bottom-up connectionsrepmat(bj, 1, numcases) + ...      %static biases on unitbjstar;                            %dynamic biaseshposteriors = 1./(1 + exp(-eta));    %logistic%Activate the hidden units    hidstates = double(hposteriors' > rand(numcases,numhid)); %Calculate positive gradients (note w.r.t. neg energy)wgrad = hidstates'*(data(:,:,1)./gsd);bigrad = sum(data(:,:,1)' - ...repmat(bi,1,numcases) - bistar,2)./gsd^2;bjgrad = sum(hidstates,1)';

可以发现此代码保留了$\sigma$，由可见层得到隐层状态的公式为
$$P(h_j=1|\mathbf{v})=sigmoid(w\ast \frac{data}{\sigma }+b_j)$$
还有一点需要注意的是，从能量函数可以看出，方差项在权重$w$和可见层偏置$b_i$都出现过，说明这两项的更新方法与$\sigma$有关
<v,h>data=h{0,1}∗dataσ<v>data=∑i(vi−bi)σ2<h>data=∑h{0,1}<v,h>_{data}=h_{\{0,1\}}*\frac{data}{\sigma}\\ <v>_{data}=\frac{\sum_i (v_i-b_i)}{\sigma^2}\\ <h>_{data}=\sum h_{\{0,1\}} <v,h>data​=h{0,1}​∗σdata​<v>data​=σ2∑i​(vi​−bi​)​<h>data​=∑h{0,1}​

• 代码2来自${}^7$

 %pass 3-way term + gated biases + hidbiases through sigmoid poshidprobs = 1./(1 + exp(-yvisfeat*hidfac'  ...-ypastfeatB*hidfacB' - repmat(hidbiases,numcases,1)));%-data*vishid - repmat(hidbiases,numcases,1)));  %Activate the hidden units    hidstates = single(poshidprobs > rand(numcases,numhid));

此代码未保留$\sigma$即$\sigma =1$，那么更新形式就较为简单
$$P(h_j=1|\mathbf{v})=sigmoid(w\ast data+b_j)$$

• 代码3摘自${}^8$，非完整摘取，只有与隐层状态计算相关的关键部分

%Term calculated by summing over hiddens
hidinp = data*vishid + effhidbias;
exphidinp = exp(hidinp); %cache this computation
poshidprobsall(:,:,cc) = exphidinp./(1+exphidinp);
poshidprobs(idx,:) = poshidprobsall(idx,:,cc);   
posprods(:,:,cc) = (data(idx,:)./gsd)'*poshidprobs(idx,:); %smoothed: probs, not binary
poshidact(cc,:) = sum(poshidprobs(idx,:),1); %col vector; again smoothed
posvisact(cc,:) = sum(data(idx,:),1)./gsd^2; %row vector

这里的激活函数依旧是sigmoid函数，只不过变了一个形式
$$sigmoid(x)=\frac{1}{1+e^{(-x)}}=\frac{e^x}{1+e^x}$$
但是有一个小疑问就是，偏置的更新变了，在上述代码5中利用的是h{0,1}h_{\{0,1\}}h{0,1}​更新，但是此代码使用 KaTeX parse error: Expected 'EOF', got '}' at position 10: P(h_i=1|v}̲ 更新
$$\begin{gathered} <v,h{>}_{data}=P(h=1|v)\ast \frac{data}{\sigma } \\ <v{>}_{data}=\frac{{\sum }_i(v_i-b_i)}{{\sigma }^2} \\ <h{>}_{data}=\sum _{j}P(h_j=1|v) \end{gathered}$$
就我个人来说，写代码的时候，经常采用的是使用smooth version 即概率 $P(h=1|v) $进行更新，这样就可以发现可见层$\to$隐层单元的更新，与二值形式是一样的

4.2.2 negative阶段

反向阶段比较模糊的从$N(\mu ,{\sigma }^2)$中采样的具体写法,其实比较建议先看我关于吉布斯采样翻译和matlab实现的那一部分代码，基本就知道如何根据均值和方差采样了，而且也有效果展示，可以明白为什么使用这个代码能够得到符合次分布的采样点。

• 看看代码1的保留$\sigma$的写法

 topdown = gsd.*(hidstates*w);%This is the mean of the Gaussian %Instead of properly sampling, negdata is just the mean%If we want to sample from the Gaussian, we would add in%gsd.*randn(numcases,numdims);negdata =  topdown + ...            %top down connectionsrepmat(bi',numcases,1) + ...    %static biasesbistar';                        %dynamic biases

可以发现$P(v_i|\mathbf{h})$的计算方法很简单
P(vi∣h)=σ∗h{0,1}∗w+biP(v_i|\mathbf{h})=\sigma*h_{\{0,1\}}*w+b_i P(vi​∣h)=σ∗h{0,1}​∗w+bi​
对于偏置的更新代码：

eta =  w*(negdata./gsd)' + ...     %bottom-up connectionsrepmat(bj, 1, numcases) + ...  %static biases on unit (no change)bjstar;                        %dynamic biases (no change)hposteriors = 1./(1 + exp(-eta));   %logistic%Calculate negative gradientsnegwgrad = hposteriors*(negdata./gsd); %not using activationsnegbigrad = sum( negdata' - ...repmat(bi,1,numcases) - bistar,2)./gsd^2;negbjgrad = sum(hposteriors,2);

可以发现是先利用一次$P(v_i|\mathbf{h})$得到新的隐状态激活概率，然后
$$\begin{gathered} <v,h{>}_{recon}=\sum _{ij}\frac{h_j\ast v_i}{\sigma } \\ <v{>}_{recon}=\sum _i\frac{v_i-b_i}{{\sigma }^2} \end{gathered}$$

• 代码2中，设置$\sigma =1$

 negdata = (yfeat.*yhid)*visfac' + ...(yfeatpastA)*visfacA' + ...repmat(visbiases,numcases,1);  

由于采用三通道，所以会多出几项，但是也很容易看出一般情况的更新方法
P(vi∣h)=h{0,1}∗w+biP(v_i|\mathbf{h})=h_{\{0,1\}}*w+b_i P(vi​∣h)=h{0,1}​∗w+bi​

• 代码3中，依旧保留了$\sigma$

 negdata(idx,:) =  gsd.*(hidstates(idx,:)*vishid(:,:,cc)') + ...effvisbiases(idx,:,cc);        

同样符合上式。结合三个代码可以发现，从高斯分布中采样的方法为
P(vi∣h)=N(ai+σi∑jWijhj,σi2)=σ∗h{0,1}∗w+ai\begin{aligned} P(v_i|\mathbf{h})&=N(a_i+\sigma_i\sum_jW_{ij}h_j,\sigma_i^2)\\ &=\sigma*h_{\{0,1\}}*w+a_i \end{aligned} P(vi​∣h)​=N(ai​+σi​j∑​Wij​hj​,σi2​)=σ∗h{0,1}​∗w+ai​​

4.3 实值RBM总结

综合上述理论及代码，我们可以得到RBM的条件分布采样方法以及参数更新公式
P(vi∣h)=σ∗h{0,1}∗w+aiP(hj=1∣v)=11+exp(−bj−∑iWijviσ)\begin{aligned} P(v_i|\mathbf{h})&=\sigma*h_{\{0,1\}}*w+a_i\\ P(h_j=1|\mathbf{v})&=\frac{1}{1+exp(-b_j-\sum_iW_{ij}\frac{v_i}{\sigma})}\\ \end{aligned} P(vi​∣h)P(hj​=1∣v)​=σ∗h{0,1}​∗w+ai​=1+exp(−bj​−∑i​Wij​σvi​​)1​​
可见层$\to$隐层时候
<v,h>data=h{0,1}∗dataσ<v>data=∑i(vi−bi)σ2<h>data=∑jhj\begin{aligned} <v,h>_{data}&=h_{\{0,1\}}*\frac{data}{\sigma}\\ <v>_{data}&=\frac{\sum_i (v_i-b_i)}{\sigma^2}\\ <h>_{data}&=\sum_j{h_j} \end{aligned} <v,h>data​<v>data​<h>data​​=h{0,1}​∗σdata​=σ2∑i​(vi​−bi​)​=j∑​hj​​
隐层$\to$可见层时候
$$\begin{array}{rl}<v,h{>}_{recon}& =\sum _{ij}\frac{h_j\ast v_i}{\sigma }\\ <v{>}_{recon}& =\sum _i\frac{v_i-b_i}{{\sigma }^2}\\ <h{>}_{recon}& =\sum _j{h}_j\end{array}$$
参数更新方法
$$\begin{array}{rl}\Delta a_i& \propto <v{>}_{data}-<v{>}_{recon}\\ \Delta b_j& \propto <h{>}_{data}-<h{>}_{recon}\\ \Delta W_{ij}& \propto <vh{>}_{data}-<vh{>}_{model}\end{array}$$

Reference
[1] The Recurrent Temporal Restricted Boltzmann Machine
[2] Two Distributed-State Models For Generating High-Dimensional Time Series
[3] Robust Generation of Dynamical Patterns in Human Motion by a Deep Belief Nets
[4] Modeling Human-Skeleton Motion Patterns Using Conditional Deep Boltzmann Machine
[5] Modeling Human Motion Using Binary Latent Variables
[6] Temporal Autoencoding Improves Generative Models of Time Series
[7] Factored Conditional Restricted Boltzmann Machines for Modeling Motion Style
[8] Implicit mixtures of Conditional Restricted Boltzmann Machines