RBM受限玻尔兹曼机的公式推导及代码实现(matlab)

  考虑一组具有 m m m个样本的数据集 X = { x ( 1 ) , … , x ( m ) } \mathbb{X}=\{x^{(1)},\dots,x^{(m)}\} X={ x(1),…,x(m)}，独立地由真实数据生成分布 p d a t a ( x ) p_{data}(x) pdata​(x)生成。令 p m o d e l ( x ; θ ) p_{model}(x;\theta) pmodel​(x;θ)是一族由 θ \theta θ确定在相同空间上的概率分布，换言之， p m o d e l ( x ; θ ) p_{model}(x;\theta) pmodel​(x;θ)将任意输入 x x x映射到实数来估计真实概率 p d a t a ( x ) p_{data}(x) pdata​(x)。对 θ \theta θ的最大似然估计被定义为 θ M L = a r g m a x θ p m o d e l ( X ; θ ) = a r g m a x θ ∏ i = 1 m p m o d e l ( x ( i ) ; θ ) (1) \theta_{ML}=argmax_\theta p_{model}(\mathbb{X};\theta)=argmax_\theta\prod_{i=1}^m p_ {model}(x^{(i)};\theta) \tag{1} θML​=argmaxθ​pmodel​(X;θ)=argmaxθ​i=1∏m​pmodel​(x(i);θ)(1)   因为多个概率的乘积不便于计算，便转为求和形式 θ M L = a r g m a x θ ∑ i = 1 m l o g p m o d e l ( x ( i ) ; θ ) (2) \theta_{ML}=argmax_\theta\sum_{i=1}^mlogp_{model}(x^{(i)};\theta)\tag{2} θML​=argmaxθ​i=1∑m​logpmodel​(x(i);θ)(2)   重新缩放代价函数时，对优化结果没有影响，可以除以 m m m，从而以训练数据经验分布 p ^ d a t a \hat{p}_{data} p^​data​的期望作为准则 θ M L = a r g m a x θ E x ∼ p ^ d a t a l o g p m o d e l ( x ; θ ) (3) \theta_{ML}=argmax_\theta\mathbb{E}_{x\sim\hat{p}_{data}}logp_{model}(x;\theta) \tag{3} θML​=argmaxθ​Ex∼p^​data​​logpmodel​(x;θ)(3)   一种关于最大似然估计的观点是看作最小化训练集上的经验分布 p ^ d a t a \hat{p}_{data} p^​data​和模型分布 p m o d e l p_{model} pmodel​之间的差异，差异可以通过 K L KL KL散度定义为 D K L ( p ^ d a t a ∣ ∣ p m o d e l ) = E x ∼ p ^ d a t a [ l o g p ^ d a t a ( x ) − l o g p m o d e l ( x ) ] (4) D_{KL}(\hat{p}_{data}||p_{model})=\mathbb{E}_{x\sim\hat{p}_{data}}[log\hat{p}_{data}(x)-logp_{model}(x)] \tag{4} DKL​(p^​data​∣∣pmodel​)=Ex∼p^​data​​[logp^​data​(x)−logpmodel​(x)](4) 左边仅涉及数据生成过程，和模型无关，所以只需要最小化 θ M L = a r g m i n θ − E x ∼ p ^ d a t a l o g p m o d e l ( x ; θ ) (5) \theta_{ML}=argmin_\theta-\mathbb{E}_{x\sim\hat{p}_{data}}logp_{model}(x;\theta) \tag{5} θML​=argminθ​−Ex∼p^​data​​logpmodel​(x;θ)(5) 我们可以将最大似然看作使模型分布尽可能和经验分布 p ^ d a t a \hat{p}_{data} p^​data​相匹配，理想情况下能够匹配真实分布 p d a t a p_{data} pdata​，但我们无法直接知道这个真实分布。   考虑用受限玻尔兹曼机RBM来对 p m o d e l p_{model} pmodel​进行建模，RBM有两层，分别称为可见层(visible layer)和隐藏层(hidden layer)，可见层为可观察变量 v v v，隐藏层为潜变量 h h h。层内无连接，层间全连接，是一个二分网络结构，所以当给定可见层神经元状态时，隐藏层各神经元条件独立，反之亦然。可见层神经单元用来描述观察数据，隐藏层神经单元可以看作特征提取层。   就像普通的玻尔兹曼机，受限玻尔兹曼机也是基于能量的模型，其联合概率分布由能量函数指定（能量函数的概念最早来自于统计热力学家研究磁体的易辛模型，后来被Hinton借鉴发展为RBM模型） P θ ( v = v , h = h ) = 1 Z e x p ( − E ( v , h ) ) = 1 Z ( θ ) e x p ( ∑ i = 1 D ∑ j = 1 F W i j v i h j + ∑ i = 1 D v i b i + ∑ j = 1 F h j a j ) (6) P_\theta(\mathtt{v}=v,\mathtt{h}=h)=\frac{1}{Z}exp(-E(v,h)) \\ =\frac{1}{Z(\theta)}exp(\sum_{i=1}^D\sum_{j=1}^FW_{ij}v_ih_j+\sum_{i=1}^Dv_ib_i+\sum_{j=1}^Fh_ja_j)\tag{6} Pθ​(v=v,h=h)=Z1​exp(−E(v,h))=Z(θ)1​exp(i=1∑D​j=1∑F​Wij​vi​hj​+i=1∑D​vi​bi​+j=1∑F​hj​aj​)(6) P θ ( v = v ) = 1 Z ( θ ) ∑ h e x p ( v T W h + c T h + b T v ) (7) P_\theta(\mathtt{v}=v)=\frac{1}{Z(\theta)}\sum_h exp(v^TWh+c^Th+b^Tv) \tag{7} Pθ​(v=v)=Z(θ)1​h∑​exp(vTWh+cTh+bTv)(7) RBM的能量函数由下给出 E ( v , h ) = − b T v − c T h − v T W h (8) E(v,h)=-b^Tv-c^Th-v^TWh \tag{8} E(v,h)=−bTv−cTh−vTWh(8) 我们通过最大似然估计来确定RBM的参数 θ = a r g m a x θ L ( θ ) = a r g m a x θ 1 N ∑ i = 1 m l o g P θ ( v ( i ) ) (9) \theta=argmax_\theta L(\theta)=argmax_\theta\frac{1}{N}\sum_{i=1}^mlogP_\theta(v^{(i)})\tag{9} θ=argmaxθ​L(θ)=argmaxθ​N1​i=1∑m​logPθ​(v(i))(9) 可以通过随机梯度下降 确定参数，首先要求 L ( θ ) L(\theta) L(θ)对 W W W的导数 ∂ L ( θ ) ∂ W i j = 1 N ∑ n = 1 N ∂ ∂ W i j l o g ( ∑ h e x p [ v ( n ) T W h + c T h + b T v ( n ) ] ) − ∂ ∂ W i j l o g Z ( θ ) = E P d a t a [ v i h j ] − E P θ [ v i h j ] = 1 N ∑ n = 1 N [ P ( h j = 1 ∣ v ( n ) ) v i ( n ) − ∑ v P ( v ) P ( h j = 1 ∣ v ) v i ] (10) \frac{\partial L(\theta)}{\partial W_{ij}}=\frac{1}{N}\sum_{n=1}^N \frac{\partial}{\partial W_{ij}}log(\sum_hexp[v^{(n)T}Wh+c^Th+b^Tv^{(n)}])-\frac{\partial}{\partial W_{ij}}logZ(\theta)\\ =E_{P_{data}}[v_ih_j]-E_{P_\theta}[v_ih_j] \\ =\frac{1}{N}\sum_{n=1}^N[P(h_j=1|v^{(n)})v_i^{(n)}-\sum_vP(v)P(h_j=1|v)v_i] \tag{10} ∂Wij​∂L(θ)​=N1​n=1∑N​∂Wij​∂​log(h∑​exp[v(n)TWh+cTh+bTv(n)])−∂Wij​∂​logZ(θ)=EPdata​​[vi​hj​]−EPθ​​[vi​hj​]=N1​n=1∑N​[P(hj​=1∣v(n))vi(n)​−v∑​P(v)P(hj​=1∣v)vi​](10) ∂ L ( θ ) ∂ b i = E P d a t a [ v i ] − E P θ [ v i ] = ∑ n = 1 N [ v i n − ∑ v P ( v ) v i ] (11) \frac{\partial L(\theta)}{\partial b_{i}}=\mathbb{E}_{P_{data}}[v_i]-\mathbb{E}_{P_\theta}[v_i] \\ =\sum_{n=1}^N[v_i^{n}-\sum_vP(v)v_i] \tag{11} ∂bi​∂L(θ)​=EPdata​​[vi​]−EPθ​​[vi​]=n=1∑N​[vin​−v∑​P(v)vi​](11) ∂ L ( θ ) ∂ c j = E P d a t a [ h j ] − E P θ [ h j ] = ∑ n = 1 N [ P ( h j = 1 ∣ v ( n ) ) − ∑ v P ( v ) P ( h i = 1 ∣ v ) ] (12) \frac{\partial L(\theta)}{\partial c_{j}}=\mathbb{E}_{P_{data}}[h_j]-\mathbb{E}_{P_\theta}[h_j] \\ =\sum_{n=1}^N[P(h_j=1|v^{(n)})-\sum_vP(v)P(h_i=1|v)] \tag{12} ∂cj​∂L(θ)​=EPdata​​[hj​]−EPθ​​[hj​]=n=1∑N​[P(hj​=1∣v(n))−v∑​P(v)P(hi​=1∣v)](12)

  上式的前一项在全部数据集上求平均值即可，后一项等于 ∑ v , h v i h j P θ ( v , h ) \sum_{\mathtt{v},\mathtt{h}}v_ih_jP_\theta(\mathtt{v},\mathtt{h}) ∑v,h​vi​hj​Pθ​(v,h)，其中 Z Z Z是被称为配分函数的归一化常数 Z = ∑ v ∑ h e x p ( − E ( v , h ) ) (13) Z=\sum_v\sum_hexp(-E(v,h)) \tag{13} Z=v∑​h∑​exp(−E(v,h))(13) 计算配分函数 Z Z Z的朴素方法是对所有状态进行穷举求和，计算上是难以处理的，Long and Servedio(2010)正式证明配分函数 Z Z Z是难解的 。但是RBM的二分网络结构具有特定性质，因为可见层和隐藏层内部各神经元是条件独立的，所以条件分布 p ( h ∣ v ) p(\mathtt{h}|\mathtt{v}) p(h∣v)和 p ( v ∣ h ) p(\mathtt{v}|\mathtt{h}) p(v∣h)是因子的，并且计算和采样相对简单。 P ( h ∣ v ) = P ( h , v ) P ( v ) = 1 P ( v ) 1 Z e x p { b T v + c T h + v T W h } = 1 Z ′ e x p { c T h + v T W h } = 1 Z ′ e x p { ∑ j = 1 n h c j T h j + ∑ n h j = 1 v T W : , j h j } = 1 Z ′ ∏ j = 1 n h e x p { c j T h j + v T W : , j h j } (14) P(h|v)=\frac{P(h,v)}{P(v)} \\ = \frac{1}{P(v)}\frac{1}{Z}exp\{b^Tv+c^Th+v^TWh\} \\ = \frac{1}{Z'}exp\{c^Th+v^TWh\} \\ = \frac{1}{Z'}exp\{\sum_{j=1}^{n_h}c_j^Th_j+\sum_{n_h}^{j=1}v^TW_{:,j}h_j\} \\ = \frac{1}{Z'}\prod_{j=1}^{n_h}exp\{c_j^Th_j+v^TW_{:,j}h_j\} \tag{14} P(h∣v)=P(v)P(h,v)​=P(v)1​Z1​exp{ bTv+cTh+vTWh}=Z′1​exp{ cTh+vTWh}=Z′1​exp{ j=1∑nh​​cjT​hj​+nh​∑j=1​vTW:,j​hj​}=Z′1​j=1∏nh​​exp{ cjT​hj​+vTW:,j​hj​}(14)   在基于可见单元 v v v计算条件概率 P ( h ∣ v ) P(h|v) P(h∣v)时，可以将其视为常数，因为 P ( h ∣ v ) P(h|v) P(h∣v)的性质，我们可以将向量 h h h上的联合概率写成单独元素 h j h_j hj​（未归一化）上分布的乘积。 P ( h j = 1 ∣ v ) = P ~ ( h j = 1 ∣ v ) P ~ ( h j = 0 ∣ v ) + P ~ ( h j = 1 ∣ v ) = e x p { c j + v T W : , j } e x p { 0 } + e x p { c j + v T W : , j } = 1 1 + e x p { − c j − ∑ i v i W i j } = σ ( c j + v T W : , j ) (15) P(h_j=1|v)=\frac{\tilde{P}(h_j=1|v)}{\tilde{P}(h_j=0|v)+\tilde{P}(h_j=1|v)} \\ =\frac{exp\{c_j+v^TW_{:,j}\}}{exp\{0\}+exp\{c_j+v^TW_{:,j}\}} \\ =\frac{1}{1+exp\{-c_j-\sum_iv_iW_{ij}\}} \\ =\sigma(c_j+v^TW_{:,j}) \tag{15} P(hj​=1∣v)=P~(hj​=0∣v)+P~(hj​=1∣v)P~(hj​=1∣v)​=exp{ 0}+exp{ cj​+vTW:,j​}exp{ cj​+vTW:,j​}​=1+exp{ −cj​−∑i​vi​Wij​}1​=σ(cj​+vTW:,j​)(15) 其中 σ ( x ) = 1 1 + e x p ( − x ) \sigma(x)=\frac{1}{1+exp(-x)} σ(x)=