配分函数|机器学习推导系列（二十六）

一、概述

对于有向概率图模型来说，由于图中存在天然的拓扑排序关系，所以有向概率图的因式分解的形式很容易写出来。而对于无向图来说就需要根据它图中的最大团来写成一个因式分解的形式，无向图模型在局部并没有表现出是一个概率模型，在整体上才表现地是一个概率模型，由此我们也就遇到了配分函数。在无向图模型的学习和评估问题中，我们会面对概率公式中的配分函数（Partition Function），往往这个配分函数是很难处理的。

对于连续或离散的高维随机变量 x ∈ R p    o r    { 0 , 1 , ⋯   , k } p x\in \mathbb{R}^{p}\; or\; \left \{0,1,\cdots ,k\right \}^{p} x∈Rpor{0,1,⋯,k}p，它可以表示成一个无向概率图，模型参数为$\theta$，它的概率公式也就可以写成以下形式：

$$P(x;\theta )=\frac{1}{Z(\theta )}\hat{P}(x;\theta )$$

其中$Z(\theta )$也就是配分函数，可以表示为：

$$Z(\theta )=\hat{P}(x;\theta )\mathrm{d}x$$

对于这个概率模型的参数估计，可以采用极大似然估计的方法，首先，我们有一些样本，表示为 X = { x i } i = 1 N X=\left \{x_{i}\right \}_{i=1}^{N} X={xi​}i=1N​，然后使用这些样本来做极大似然估计：

$$\begin{gathered} \hat{\theta }=\underset{\theta }{\mathrm{argmax}}P(X;\theta ) \\ =\underset{\theta }{\mathrm{argmax}}\prod _{i=1}^{N}P(x_i;\theta ) \\ =\underset{\theta }{\mathrm{argmax}}log\prod _{i=1}^{N}P(x_i;\theta ) \\ =\underset{\theta }{\mathrm{argmax}}\sum _{i=1}^{N}logP(x_i;\theta ) \\ =\underset{\theta }{\mathrm{argmax}}\sum _{i=1}^N\left(log\hat{P}(x_i;\theta )-logZ(\theta )\right) \\ =\underset{\theta }{\mathrm{argmax}}\sum _{i=1}^{N}log\hat{P}(x_i;\theta )-NlogZ(\theta ) \\ =\underset{\theta }{\mathrm{argmax}}\underset{记作l(\theta )}{\underbrace{\frac{1}{N}\sum _{i=1}^{N}log\hat{P}(x_i;\theta )-logZ(\theta )}} \end{gathered}$$

这里我们也就得到了目标函数$l(\theta )$：

$$l(\theta )=\frac{1}{N}\sum _{i=1}^{N}log\hat{P}(x_i;\theta )-logZ(\theta )$$

接下来使用梯度上升的方法来求解参数$\theta$，因此也就需要对$l(\theta )$求导：

$${\nabla }_{\theta }l(\theta )=\underset{①}{\underbrace{\frac{1}{N}\sum _{i=1}^N{\nabla }_{\theta }log\hat{P}(x_i;\theta )}}-\underset{②}{\underbrace{{\nabla }_{\theta }logZ(\theta )}}$$

这里我们首先看一下②这一项的求导：

$$\begin{gathered} ②={\nabla }_{\theta }logZ(\theta ) \\ =\frac{1}{Z(\theta )}{\nabla }_{\theta }Z(\theta ) \\ =\frac{P(x;\theta )}{\hat{P}(x;\theta )}{\nabla }_{\theta }\int \hat{P}(x;\theta )\mathrm{d}x \\ =\frac{P(x;\theta )}{\hat{P}(x;\theta )}\int {\nabla }_{\theta }\hat{P}(x;\theta )\mathrm{d}x \\ =\int \frac{P(x;\theta )}{\hat{P}(x;\theta )}{\nabla }_{\theta }\hat{P}(x;\theta )\mathrm{d}x \\ =\int P(x;\theta ){\nabla }_{\theta }log\hat{P}(x;\theta )\mathrm{d}x \\ =E_{P(x;\theta )}\left[{\nabla }_{\theta }log\hat{P}(x;\theta )\right] \end{gathered}$$

注意这里的$\frac{P(x;\theta )}{\hat{P}(x;\theta )}$之所以能够放到积分号里面，是因为对于任意$x$来说$\frac{P(x;\theta )}{\hat{P}(x;\theta )}$都是个常数。

上面式子之所以这么变换的目的也就是要将这个导数写成关于$P(x;\theta )$的期望的形式，然而$P(x;\theta )$正是我们要求解的分布，是一个未知的分布，因此没办法精确求解，也就没办法计算这个梯度${\nabla }_{\theta }l(\theta )$，只能近似采样。如果没有②这一项，$l(\theta )$就可以采用梯度上升，但是由于存在配分函数就无法直接采用梯度上升了。以上就是问题所在。

二、随机最大似然（Stochastic Maximum Likelihood）

现在我们可以把${\nabla }_{\theta }l(\theta )$下面的形式：

$${\nabla }_{\theta }l(\theta )=\frac{1}{N}\sum _{i=1}^N{\nabla }_{\theta }log\hat{P}(x_i;\theta )-E_{P(x;\theta )}\left[{\nabla }_{\theta }log\hat{P}(x;\theta )\right]$$

对于我们的训练数据，我们可以把它们看做服从一个分布$P_{data}$，这个分布也就是训练数据的真实分布，是我们要去拟合的分布，然而也是一个我们永远不可能真正精确得到的分布，因为我们能够拿到的只有这个分布的若干样本（也就是训练数据）。有了$P_{data}$这个定义以后，我们也可以把上面梯度中减号左边的一项写成期望的形式：

$${\nabla }_{\theta }l(\theta )=E_{P_{data}}\left[{\nabla }_{\theta }log\hat{P}(x;\theta )\right]-E_{P(x;\theta )}\left[{\nabla }_{\theta }log\hat{P}(x;\theta )\right]$$

我们的目的也就是要让$P(x;\theta )$尽可能地逼近$P_{data}$，这个$P(x;\theta )$也就是我们的模型，所以我们把$P(x;\theta )$这个分布记作$P_{model}$，现在我们就有了两个定义的分布，即数据的真实分布$P_{data}$和用来拟合$P_{data}$的$P_{model}$：

$$\{\begin{array}{c}datadistribution:P_{data}\\ modeldistribution:P_{model}\triangleq P(x;\theta )\end{array}$$

现在${\nabla }_{\theta }l(\theta )$就可以写成${\nabla }_{\theta }log\hat{P}(x;\theta )$关于$P_{data}$和$P_{model}$的期望的形式：

$${\nabla }_{\theta }l(\theta )=\underset{positivephase}{\underbrace{E_{P_{data}}\left[{\nabla }_{\theta }log\hat{P}(x;\theta )\right]}}-\underset{negativephase}{\underbrace{E_{P_{model}}\left[{\nabla }_{\theta }log\hat{P}(x;\theta )\right]}}$$

这里分别定义等号左边和右边的部分为正相（positive phase）和负相（negative phase）。

我们期待使用梯度上升法来迭代地求解最优的$\theta$：

$${\theta }^{(t+1)}={\theta }^{(t)}+\eta {\nabla }_{\theta }l({\theta }^{(t)})$$

对于负相来说，我们无法对这个期望值积分，因此只能采用近似的方法，也就是在每次迭代时利用MCMC的方法（比如吉布斯采样）从$P_{model}=P(x;{\theta }^{(t)})$里采样得到样本，然后利用这些样本来近似计算负相这个积分值。这里采样得到的$M$个样本记作 { x ^ i } i = 1 m \left \{\hat{x}_{i}\right \}_{i=1}^{m} {x^i​}i=1m​，这些样本叫做幻想粒子（fantacy particles）:

$$\begin{array}{c}{\hat{x}}_1\sim P(x;{\theta }^{(t)})\\ ⋮\\ {\hat{x}}_m\sim P(x;{\theta }^{(t)})\end{array}\}fantacyparticles$$

有了这些样本就可以计算负相，也就是说现在就可以用梯度上升的方法来迭代求解参数$\theta$了，以下就是迭代求解的公式（这里也会从训练集中抽$m$个样本）：

$${\theta }^{(t+1)}={\theta }^{(t)}+\eta \left(\sum _{i=1}^m{\nabla }_{\theta }log\hat{P}(x_i;{\theta }^{(t)})-\sum _{i=1}^m{\nabla }_{\theta }log\hat{P}({\hat{x}}_i;{\theta }^{(t)})\right)$$

这个方法就叫做Gradient Ascent based on MCMC。

上面介绍的内容中引入了几个不容易理解的名词：正相、负相和幻想粒子，现在可以直观地解释一下这几个名词的含义。${\nabla }_{\theta }l({\theta }^{(t)})$中存在正相和负相，对比目标函数$l(\theta )$的表达式可以直观地理解这样命名的用意。

正相的作用是让模型分布在训练样本处概率增大，也就是从$P_{data}$中采样，然后让这些样本在$P_{model}$中的概率增大。负相的作用是让$Z(\theta )$的值变小，也就是说要让从$P_{model}$中采样出来的幻想粒子在$P_{model}$中的概率减小，这些样本可以认为我们是不信任它的，称它们是fantasy的，因此要让它们的概率减小。正相和负相共同作用，最终结果就会让$P_{model}$逼近$P_{data}$，这个过程可以由下图表示：

可以想象如果$P_{model}$已经非常逼近$P_{data}$，那么采样得到的幻想粒子和从数据集中采样的样本就会非常一致，这时对这些样本既要增大它们的概率也要压低它们的概率，此时正相和负相的作用就会抵消，也就不会再产生梯度，训练也就必须停止。

三、对比散度

对于MCMC的方法，可以参考这两个链接：
①MCMC-1|机器学习推导系列（十五）
②MCMC-2|机器学习推导系列（十六）

上面的目标函数需要从$P_{data}$和$P_{model}$里面都采样$m$个样本，从$P_{data}$里面采样是很容易的，只需要从训练数据中抽取$m$个训练数据，而从$P_{model}$中采样就要采用MCMC的方法，具体的操作就是使用一个初始分布初始化马尔可夫链，然后等马尔可夫链随机游走到达平稳分布时进行采样，这里可以构建一条马尔可夫链然后从中采集$m$个样本，也可以构建$m$条马尔可夫链然后从每条马尔可夫链中采集一个样本（只不过这样比较消耗资源）。上述MCMC方法的问题是对于高维数据分布，很可能马尔可夫链的混合时间（或者叫燃烧期）会非常地长。

下图展示了采样多条马尔可夫链的吉布斯采样方法的过程，假设燃烧期是$k$，由于数据维度过高，很可能$k$就会非常大，我们可以认为是无穷大$\infty $。这里也要注意，吉布斯采样的过程是对高维随机变量的每一维依次采样，因此这里的图中的每个step其实都表示高维随机变量每一维采样的过程：

$$\begin{array}{c}\stackrel{mixingtime}{\overbrace{\underset{0-step}{◯}\to \underset{1-step}{◯}\to \cdots \underset{k-step}{◯}}}\to \cdots ◯\to {\hat{x}}_1\\ ⋮\\ \underset{0-step}{◯}\to \underset{1-step}{◯}\to \cdots \underset{k-step}{◯}\to \cdots ◯\to {\hat{x}}_m\end{array}\}GibbsSampling$$

对比散度（Contrastive Divergence）的方法可以避免燃烧期过长的问题，在上面的抽样过程中需要为0-step的${\hat{x}}_i$按照一个初始化分布进行初始化，而对比散度的方法就是使用$P_{data}$这个分布来作为初始化分布，具体的做法也就是从训练数据中抽取$m$个样本来初始化这$m$条马尔可夫链，也就是：

$${\hat{x}}_1=x_1,{\hat{x}}_2=x_2,\cdots ,{\hat{x}}_m=x_m$$

使用$P_{data}$初始化马尔可夫链以后只需要经过很短的几步就可以采集样本了，比如经过$k$步就开始采样，即使$k=1$也是可以的。总之，对比散度的方法与普通的直接MCMC采样的方法的区别就在于使用了$P_{data}$初始化马尔可夫链，然后不用等漫长的混合时间，只需要$k$步就可以采样了，无论$k$步时有没有达到平稳分布，而$k=1,2,3$等很小的数字也是可以的。这种对比散度的方法就叫做CD Learning，之前的方法叫做ML Learning。

接下来来看看对比散度这个名字的由来。首先，先看这个式子：

$$\begin{gathered} \hat{\theta }=\underset{\theta }{\mathrm{argmax}}\frac{1}{N}\sum _{i=1}^{N}logP(x;\theta ) \\ =\underset{\theta }{\mathrm{argmax}}{E}_{P_{data}}\left[log{P}_{model}\right] \\ =\underset{\theta }{\mathrm{argmax}}\int P_{data}log{P}_{model}\mathrm{d}x \\ =\underset{\theta }{\mathrm{argmax}}\int P_{data}log{P}_{model}\mathrm{d}x-\underset{与\theta 无关}{\underbrace{\underset{\theta }{\mathrm{argmax}}\int P_{data}log{P}_{data}\mathrm{d}x}} \\ =\underset{\theta }{\mathrm{argmax}}\int P_{data}log\frac{P_{model}}{P_{data}}\mathrm{d}x \\ =\underset{\theta }{\mathrm{argmax}}-KL(P_{data}||P_{model}) \\ =\underset{\theta }{\mathrm{argmin}}KL(P_{data}||P_{model}) \end{gathered}$$

可以看到，原来的极大似然估计的方法其实是在最小化$P_{data}$和$P_{model}$的KL散度。而现在在对比散度的方法中，我们用$P_{data}$作为第$0$步的初始化分布，而在经过很长的步数后才能达到平稳分布，我们现在把第$0$步的分布记作$P^{(0)}$，最终的平稳分布记作$P^{(\infty )}$，而中间的第$k$步的马尔可夫链的分布记作$P^{(k)}$。由于现在$P_{model}$采样的样本来自$P^{(k)}$而非$P^{(\infty )}$，所以通过这些样本进行的参数估计就不是在最小化$KL(P^{(0)}||P^{(\infty )})$（也就是$KL(P_{data}||P_{model})$）了，而是在按照下面式子中的目标函数进行参数估计：

$$\hat{\theta }=argmin\underset{ContrastiveDivergence}{\underbrace{\left[KL(P^{(0)}||P^{(\infty )})-KL(P^{(k)}||P^{(\infty )})\right]}}$$

这个目标函数就是对比散度。使用CD-Learning的方法的算法如下：

$t+1$时刻:
①为正相从$P_{data}$中采样，$x_1,x_2,\cdots ,x_m$是采样的数据，也就是训练数据；
②为负相从$P_{model}=P(x;{\theta }^{(t)})$中采样，使用为正相采样的数据初始化马尔可夫链：
$${\hat{x}}_1=x_1,{\hat{x}}_2=x_2,\cdots ,{\hat{x}}_m=x_m$$
然后使用吉布斯采样从马尔可夫链中第$k$步的分布抽取$m$个样本，$k$可以是很小的数字，甚至是$1$：
$$\begin{array}{c}\underset{0-step}{◯}\to \underset{1-step}{◯}\to \cdots \underset{k-step}{◯}\to {\hat{x}}_1\\ ⋮\\ \underset{0-step}{◯}\to \underset{1-step}{◯}\to \cdots \underset{k-step}{◯}\to {\hat{x}}_m\end{array}\}GibbsSampling$$

四、受限玻尔兹曼机的学习

1. 表示

受限玻尔兹曼机在前一篇介绍了它的表示和推断问题，参考链接如下：受限玻尔兹曼机|机器学习推导系列（二十五）。

它的概率模型如下：

{ P ( h , v ) = 1 Z e x p { − E ( h , v ) } E ( h , v ) = − ( h T W v + α T v + β T h ) \left\{\begin{matrix} P(h,v)=\frac{1}{Z}exp\left \{-E(h,v)\right \}\\ E(h,v)=-\left (h^{T}Wv+\alpha ^{T}v+\beta ^{T}h\right ) \end{matrix}\right. {P(h,v)=Z1​exp{−E(h,v)}E(h,v)=−(hTWv+αTv+βTh)​

其中：

$$\begin{gathered} h={\left(\begin{array}{cccc}{h}_1& h_2& \cdots & h_m\end{array}\right)}^T \\ v={\left(\begin{array}{cccc}{v}_1& v_2& \cdots & v_n\end{array}\right)}^T \\ W={\left[w_{ij}\right]}_{m\times n} \\ \alpha ={\left(\begin{array}{cccc}{\alpha }_1& {\alpha }_2& \cdots & {\alpha }_n\end{array}\right)}^T \\ \beta ={\left(\begin{array}{cccc}{\beta }_1& {\beta }_2& \cdots & {\beta }_m\end{array}\right)}^T \end{gathered}$$

2. 具有隐变量的能量模型

对于具有隐变量的能量模型来说，我们有的数据是观测变量$v$的数据，假设数据集是$S$，$S$的规模是$N$，即$\left|S\right|=N$，另外用$\theta$表示参数$(W,\alpha ,\beta )$，那么数据的$log$似然就是：

$$l(\theta )=\frac{1}{N}\sum _{v\in S}logP(v)$$

对于概率$logP(v)$，来求它对$\theta$的梯度，首先对这个概率做一些变换：

l o g P ( v ) = l o g ∑ h P ( h , v ) = l o g ∑ h 1 Z e x p { − E ( h , v ) } = l o g ∑ h e x p { − E ( h , v ) } − l o g Z = l o g ∑ h e x p { − E ( h , v ) } ⏟ 记 作 ① − l o g ∑ h , v e x p { − E ( h , v ) } ⏟ 记 作 ② logP(v)=log\sum _{h}P(h,v)\\ =log\sum _{h}\frac{1}{Z}exp\left \{-E(h,v)\right \}\\ =log\sum _{h}exp\left \{-E(h,v)\right \}-logZ\\ =\underset{记作①}{\underbrace{log\sum _{h}exp\left \{-E(h,v)\right \}}}-\underset{记作②}{\underbrace{log\sum _{h,v}exp\left \{-E(h,v)\right \}}} logP(v)=logh∑​P(h,v)=logh∑​Z1​exp{−E(h,v)}=logh∑​exp{−E(h,v)}−logZ=记作① logh∑​exp{−E(h,v)}​​−记作② logh,v∑​exp{−E(h,v)}​​

接下来看①和②这两项对参数$\theta$的导数：

∂ ① ∂ θ = ∂ ∂ θ l o g ∑ h e x p { − E ( h , v ) } = − 1 ∑ h e x p { − E ( h , v ) } ∑ h e x p { − E ( h , v ) } ∂ E ( h , v ) ∂ θ = − ∑ h e x p { − E ( h , v ) } ∑ h e x p { − E ( h , v ) } ∂ E ( h , v ) ∂ θ = − ∑ h 1 Z e x p { − E ( h , v ) } 1 Z ∑ h e x p { − E ( h , v ) } ∂ E ( h , v ) ∂ θ = − ∑ h P ( h , v ) ∑ h P ( h , v ) ∂ E ( h , v ) ∂ θ = − ∑ h P ( h ∣ v ) ∂ E ( h , v ) ∂ θ ∂ ② ∂ θ = ∂ ∂ θ l o g ∑ h , v e x p { − E ( h , v ) } = − 1 ∑ h , v e x p { − E ( h , v ) } ∑ h , v e x p { − E ( h , v ) } ∂ E ( h , v ) ∂ θ = − ∑ h , v e x p { − E ( h , v ) } ∑ h , v e x p { − E ( h , v ) } ∂ E ( h , v ) ∂ θ = − ∑ h , v e x p { − E ( h , v ) } Z ∂ E ( h , v ) ∂ θ = − ∑ h , v P ( h , v ) ∂ E ( h , v ) ∂ θ \frac{\partial ①}{\partial \theta }=\frac{\partial}{\partial \theta }log\sum _{h}exp\left \{-E(h,v)\right \}\\ =-\frac{1}{\sum _{h}exp\left \{-E(h,v)\right \}}\sum _{h}exp\left \{-E(h,v)\right \}\frac{\partial E(h,v)}{\partial \theta }\\ =-\sum _{h}\frac{exp\left \{-E(h,v)\right \}}{\sum _{h}exp\left \{-E(h,v)\right \}}\frac{\partial E(h,v)}{\partial \theta }\\ =-\sum _{h}\frac{{\color{Red}{\frac{1}{Z}}}exp\left \{-E(h,v)\right \}}{{\color{Red}{\frac{1}{Z}}}\sum _{h}exp\left \{-E(h,v)\right \}}\frac{\partial E(h,v)}{\partial \theta }\\ =-\sum _{h}\frac{P(h,v)}{\sum _{h}P(h,v)}\frac{\partial E(h,v)}{\partial \theta }\\ =-\sum _{h}P(h|v)\frac{\partial E(h,v)}{\partial \theta }\\ \frac{\partial ②}{\partial \theta }=\frac{\partial}{\partial \theta }log\sum _{h,v}exp\left \{-E(h,v)\right \}\\ =-\frac{1}{\sum _{h,v}exp\left \{-E(h,v)\right \}}\sum _{h,v}exp\left \{-E(h,v)\right \}\frac{\partial E(h,v)}{\partial \theta }\\ =-\sum _{h,v}\frac{exp\left \{-E(h,v)\right \}}{\sum _{h,v}exp\left \{-E(h,v)\right \}}\frac{\partial E(h,v)}{\partial \theta }\\ =-\sum _{h,v}\frac{exp\left \{-E(h,v)\right \}}{Z}\frac{\partial E(h,v)}{\partial \theta }\\ =-\sum _{h,v}P(h,v)\frac{\partial E(h,v)}{\partial \theta } ∂θ∂①​=∂θ∂​logh∑​exp{−E(h,v)}=−∑h​exp{−E(h,v)}1​h∑​exp{−E(h,v)}∂θ∂E(h,v)​=−h∑​∑h​exp{−E(h,v)}exp{−E(h,v)}​∂θ∂E(h,v)​=−h∑​Z1​∑h​exp{−E(h,v)}Z1​exp{−E(h,v)}​∂θ∂E(h,v)​=−h∑​∑h​P(h,v)P(h,v)​∂θ∂E(h,v)​=−h∑​P(h∣v)∂θ∂E(h,v)​∂θ∂②​=∂θ∂​logh,v∑​exp{−E(h,v)}=−∑h,v​exp{−E(h,v)}1​h,v∑​exp{−E(h,v)}∂θ∂E(h,v)​=−h,v∑​∑h,v​exp{−E(h,v)}exp{−E(h,v)}​∂θ∂E(h,v)​=−h,v∑​Zexp{−E(h,v)}​∂θ∂E(h,v)​=−h,v∑​P(h,v)∂θ∂E(h,v)​

那么最终$logP(v)$对参数$\theta$的梯度为：

$$\begin{gathered} \frac{\partial }{\partial \theta }logP(v)=\frac{\partial ①}{\partial \theta }-\frac{\partial ②}{\partial \theta } \\ =\sum _{h,v}P(h,v)\frac{\partial E(h,v)}{\partial \theta }-\sum _{h}P(h|v)\frac{\partial E(h,v)}{\partial \theta } \end{gathered}$$

3. RBM的$log$似然梯度

接下来以求解$W$为例来看RBM的参数学习方法。上面有了$logP(v)$对参数$\theta$的梯度，类似地$logP(v)$对参数$w_{ij}$的梯度为：

$$\frac{\partial }{\partial w_{ij}}logP(v)=\sum _{h,v}P(h,v)\frac{\partial E(h,v)}{\partial w_{ij}}-\sum _{h}P(h|v)\frac{\partial E(h,v)}{\partial w_{ij}}$$

然后观察一下能量函数$E(h,v)$：

$$\begin{gathered} E(h,v)=-(h^{T}Wv+\underset{与W无关，记作\Delta }{\underbrace{{\alpha }^{T}v+{\beta }^{T}h}}) \\ =-(h^{T}Wv+\Delta ) \\ =-(\sum _{i=1}^m\sum _{i=1}^n{h}_i{w}_{ij}{v}_j+\Delta ) \end{gathered}$$

那么$\frac{\partial E(h,v)}{\partial w_{ij}}$也就很容易写出来：

$$\frac{\partial E(h,v)}{\partial w_{ij}}=-h_i{v}_j$$

代入$\frac{\partial }{\partial w_{ij}}logP(v)$就有：

$$\begin{gathered} \frac{\partial }{\partial w_{ij}}logP(v)=\sum _{h,v}P(h,v)(-h_i{v}_j)-\sum _{h}P(h|v)(-h_i{v}_j) \\ =\underset{记作①}{\underbrace{\sum _{h}P(h|v)h_i{v}_j}}-\underset{记作②}{\underbrace{\sum _{h,v}P(h,v)h_i{v}_j}} \end{gathered}$$

有一点要回忆一下，就是RBM无论隐变量还是观测变量都是二值的，取值只能是$0$或$1$。接下来对①和②继续进行变换，需要利用这一点：

$$\begin{gathered} ①=\sum _{h_1}\sum _{h_2}\cdots \sum _{h_i}\cdots \sum _{h_m}P(h_1,h_2,\cdots ,h_i,\cdots ,h_m|v)h_i{v}_j \\ =\sum _{h_i}P(h_i|v)h_i{v}_j \\ =P(h_i=1|v)1\cdot v_j+P(h_i=0|v)0\cdot v_j \\ =P(h_i=1|v)v_j \\ ②=\sum _h\sum _{v}P(h,v)h_i{v}_j \\ =\sum _h\sum _{v}P(v)P(h|v)h_i{v}_j \\ =\sum _{v}P(v)\underset{同①}{\underbrace{\sum _{h}P(h|v)h_i{v}_j}} \\ =\sum _{v}P(v)P(h_i=1|v)v_j \end{gathered}$$

因此也就有：

$$\frac{\partial }{\partial w_{ij}}logP(v)=P(h_i=1|v)v_j-\sum _{v}P(v)P(h_i=1|v)v_j$$

对于上面式子中的$P(h_i=1|v)$这个条件概率，我们是可以直接写出来的，这个概率在上一篇的推断问题中已经推导过了，可以参照本节开头的链接。

4. RBM的CD-k方法

所有样本的$log$似然$l(\theta )$对$w_{ij}$的梯度现在就可以表示为：

$$\begin{gathered} \frac{\partial l(\theta )}{\partial w_{ij}}=\frac{1}{N}\sum _{v\in S}\frac{\partial }{\partial w_{ij}}logP(v) \\ =\frac{1}{N}\sum _{v\in S}(P(h_i=1|v)v_j-\underset{E_{P(v)}\left[P(h_i=1|v)v_j\right]}{\underbrace{\sum _{v}P(v)P(h_i=1|v)v_j}}) \end{gathered}$$

这里红色的项需要对$v$进行积分，是untrackable的，这一项其实也就是关于$P(v)$的期望，因此需要借助MCMC（这里指CD-k的方法，因为RBM里的随机变量也是高维的，只用MCMC也会面临燃烧期过长的问题）的方法。

使用的采样方法还是吉布斯采样，这里具体的采样过程如下图所示，首先使用训练数据来初始化$v^{(0)}$，然后固定$v^{(0)}$来依次采样$h^{(0)}$的每一维，然后固定$h^{(0)}$再来依次采样$v^{(1)}$的每一维，按照如此流程进行$k$步，最终采样得到$v^{(k)}$这一个样本，这种固定一部分来采另一部分的方法叫做块吉布斯采样（block Gibbs Sampling）：

需要注意的是虽然按照吉布斯采样的方法需要依次采集每一个维度，但在RBM中由于模型的特殊性，在固定$v$或者$h$时，其余的随机变量都是相互独立的，因此实际操作中并行采每一维也是可以的。

另外我们一共有$N$个训练数据，因此使用每个训练数据初始化一次都可以采集到一个$v^{(k)}$，因此最终采集到的$v^{(k)}$一共有$N$个。

具体的CD-k for RBM算法为：

for each $v$：
  $v^{(0)}\leftarrow v$
  for $l=0,1,2,\cdots ,k-1$：
    for $i=1,2,\cdots ,m$：sample $h_i^{(l)}\sim P(h_i|v^{(l)})$
    for $j=1,2,\cdots ,n$：sample $v_j^{(l+1)}\sim P(v_j|h^{(l)})$
  for $i=1,2,\cdots ,m$,$j=1,2,\cdots ,n$：
    $\Delta w_{ij}\leftarrow \Delta w_{ij}+\frac{\partial }{\partial w_{ij}}logP(v)$

这里的$\Delta w_{ij}$初始化为$0$，这里的$\frac{\partial }{\partial w_{ij}}logP(v)$也就是：

$$\frac{\partial }{\partial w_{ij}}logP(v)\approx P(h_i=1|v^{(0)})v_j^{(0)}-P(h_i=1|v^{(k)})v_j^{(k)}$$

也就是说用每个训练数据采样得到的样本对应的$P(h_i=1|v^{(k)})v_j^{(k)}$来代替期望$E_{P(v)}\left[P(h_i=1|v)v_j\right]$，并且将所有训练数据计算得到的$\frac{\partial }{\partial w_{ij}}logP(v)$累计起来得到的$\Delta w_{ij}$作为$\frac{\partial l(\theta )}{\partial w_{ij}}$的近似值，即：

$$\frac{\partial l(\theta )}{\partial w_{ij}}=\frac{1}{N}\sum _{v\in S}\frac{\partial }{\partial w_{ij}}logP(v)\approx \frac{1}{N}\Delta w_{ij}$$

最后进行梯度上升迭代求解就可以了。