【总结】关于玻尔兹曼机(BM)、受限玻尔兹曼机(RBM)、深度玻尔兹曼机(DBM)、深度置信网络(DBN)理论总结和代码实践

前言

前段时间，团队被老师安排了一个涉及玻尔兹曼机及相关变体的任务。具体内容就是学习相关理论知识及代码实现。
所以就想着写一篇博客来总结一下繁杂的知识点（尤其是背后的数学推导公式），如果你是玻尔兹曼机及相关变体的初学者，并对这些知识非常感兴趣，那么恭喜你，发现了这篇文章。我会用非常白话、通俗易懂的方式向读者解释这些知识逻辑。
在此特别感谢王志强、董倩等研究生，对此篇博客做出的贡献。

玻尔兹曼机(BM)

玻尔兹曼机(Boltzman Machine)是一个随机动力系统，每个变量的状态都以一定的概率受到其他变量的影响。
玻尔兹曼机可以用概率无向图模型来描述，一个具有K个节点的玻尔兹曼机满足以下三个性质：

1. 二值化。每个节点的状态值只有0和1。
2. 一个玻尔兹曼机包括两类节点，一类是可观察的节点有N个，一类是不可观察的节点，即隐藏节点，有(K-N)个。
3. 节点之间是全连接的。每个节点都和其他节点连接。
4. 每两个变量之间的互相影响是对称的。这里的对称和上面无向其实是一个概念，说白了就是已知A点的状态值，那么求B的状态值和已知B的状态值，求A的状态值的影响是相等的。如果你还是没有理解这句话，碰巧你又了解过一点概率论的知识，那么你可以将上述理解为 P(B|A) = P(A|B)。

上图就是一个有六个节点的玻尔兹曼机。其中有三个可观察的节点，我已经标了黄色，还有三个不可观测的节点，即隐藏节点，我已经标了灰色。

波尔兹曼分布推导过程

玻尔兹曼机中，随机向量X的联合概率，也就是节点的状态值，是满足玻尔兹曼分布的。

玻尔兹曼分布是描述粒子处于特定状态下的概率，是关于状态能量E(x)与系统温度T的函数。一个粒子处于状态α的概率P(α)是关于状态能量E(x)与系统温度T的函数。

别急，我相信你读完这个定义一定是懵逼，心里已经开始呐喊，TMD啥叫特定状态啊？啥叫特定状态下的概率啊？啥是状态能量啊？…
不急，我这就用 “人话” 翻译一遍。
特定状态就是说，节点的状态值为1，还是0。x=α
特定状态就是说，节点状态值为1或者0时的概率。P(x=α)
状态能量就是说，粒子本身具有的能量。
玻尔兹曼分布就是计算P(x=α)时具体的概率函数与系统的状态能量E(x)和系统温度T的函数有关。具体表达式为：

通常玻尔兹曼分布还有另一个表达式，实则是上式的等价处理，如下：

E(x)为能量函数，T为系统温度，Z为配分函数，实则就是一个归一化因子。
从玻尔兹曼分布的定义，我们可以发现系统的两个不同状态的概率之比仅与系统能量有关：
$$P(a)/P(b)=exp((E1-E2)/KT)$$
具体的推导公式如下图：

吉布斯采样

我采用的是Python下的tensooflow语言来实现的玻尔兹曼机。请注意，我们使用的版本是tensorflow2。
在讲代码之前，有必要讲一下玻尔兹曼机的训练过程。
在玻尔兹曼机中，配分函数Z通常很难计算，因此，联合概率分布P(x)一般通过马尔科夫链蒙特卡洛方法(MCMC方法)来做近似计算。
玻尔兹曼机采用了基于吉布斯采样的样本生成方法来训练的。

吉布斯采样
玻尔兹曼机的吉布斯采样过程为：随机选择一个变量Xi，然后根据其全条件概率P(Xi|X-i)来设置其状态值，即以P(Xi=1|X-i)的概率将变量Xi设为1，否则为0。在固定的温度T下，运行足够时间后，玻尔兹曼机会达到热平衡状态。此时，任何全局状态的概率都服从玻尔兹曼分布P(x)，只和系统的能量有关，和初始状态无关。

如果你听不懂上述定义，(针对没有概率论和热力学基础的读者)。那你只需要了解，因为玻尔兹曼机的联合概率函数中配分函数Z很难处理，就采用另一种方法来将节点的状态值概率趋近于玻尔兹曼分布。

玻尔兹曼机可以解决两类问题，一类是搜索问题：当给定变量之间的连接权重时，需要找到一组二值向量，使得整个网络的能量最低。另一类是学习问题，当给定变量的多组观测值时，学习网络的最优权重。

受限玻尔兹曼机(RBM)

全连接的玻尔兹曼机在理论上固然有趣，但是由于其复杂性，目前为止并没有广泛运用。实际应用中，用得更多的是基于玻尔兹曼机改造的一个版本——受限玻尔兹曼机(RBM)，其网络架构如下：

玻尔兹曼机没有层的概念，它所有的节点都是全连接的。
但受限玻尔兹曼机有层的概念。它有两层，一层称为显层，用于观测和输入，一层为隐藏层，用于提取特征。
受限玻尔兹曼机相比玻尔兹曼机，层间的节点还是采用了对称的全连接的方式连接，但是层内的节点相互独立，相互不受影响。
因为层内节点相互独立，那么由Bayes条件独立定理，受限玻尔兹曼机可以并行地对所有的显层变量或隐藏层变量同时进行采样，从而更快达到热平衡。

能量函数

由于受限玻尔兹曼机变成了层的结构，所以受限玻尔兹曼机的能量函数也变成了由三部分组成，
一个是显层节点偏置乘以显层随机可观测变量部分，
一个是连接权重与显层随机可观测变量和隐层随机可观测变量相乘部分，
一个是隐层节点偏置乘以隐层随机可观测变量偏置部分。

CD学习算法

由于受限玻尔兹曼机的特殊结构，G·Hinton提出了一种比吉布斯采样更加有效的学习算法，即对比散度学习算法，又称为CD学习算法。
通过对CD学习算法的学习，我发现这个CD算法就是在吉布斯采样的基础上作出的一点改进，即在处理玻尔兹曼机时，运行无穷次的吉布斯采样改进为运行K次即可。
以前处理玻尔兹曼机时，吉布斯采样是一直对这个玻尔兹曼机处理，直到这个波尔兹曼机收敛。G·Hinton提出，在受限玻尔兹曼机中，不需要等到受限玻尔兹曼机完全收敛，只需要K步吉布斯采样，这时模型就非常好了。所以CD算法又称K步吉布斯采样法。

代码实现受限玻尔兹曼机

受限玻尔兹曼机有两个偏置项，隐藏层的偏置项有助于RBM在前向传递中获得非零激活值，而可见层的偏置项有助于受限玻尔兹曼机学习后向传递中的重建。
在正向传递中，每个输入数据乘以一个独立的权重，然后相加后再加上一个偏置项，最后将结果传递到激活函数来产生输出。
用对比散度计算正反向梯度，然后更新偏置和权重
因为最开始受限玻尔兹曼机权重是随机初始化的，所以重建结果和原始输入差距通常会比较大，这个差距可看作是重建误差，训练受限玻尔兹曼机是通过在可见层和隐藏层之间迭代学习不断正向反向传播，直至达到某个误差的最小值
我们采用了Python环境下Numpy工具库来撰写代码，具体注释已在代码中标注，不做过多讲解。

import numpy
class RBM:
    def __init__(self, n_visible, n_hidden):
        self.n_visible = n_visible
        self.n_hidden = n_hidden
        self.bias_a = np.zeros(self.n_visible)  # 可视层偏移量
        self.bias_b = np.zeros(self.n_hidden)  # 隐藏层偏移量
        self.weights = np.random.normal(0, 0.01, size=(self.n_visible, self.n_hidden))
        self.n_sample = None

    def encode(self, v):
        # 编码，即基于v计算h的条件概率：p(h=1|v)
        return sigmoid(self.bias_b + v @ self.weights)

    def decode(self, h):
        # 解码(重构)：即基于h计算v的条件概率：p(v=1|h)
        return sigmoid(self.bias_a + h @ self.weights.T)

    def gibbs_sample(self, v0, max_cd):
        # gibbs采样, 返回max_cd采样后的v以及h值
        v = v0
        for _ in range(max_cd):
            # 首先根据输入样本对每个隐藏层神经元采样。二项分布采样，决定神经元是否激活
            ph = self.encode(v)
            h = np.random.binomial(1, ph, (self.n_sample, self.n_hidden))
            # 根据采样后隐藏层神经元取值对每个可视层神经元采样
            pv = self.decode(h)
            v = np.random.binomial(1, pv, (self.n_sample, self.n_visible))
        return v
        
    def update(self, v0, v_cd, eta):
        # 根据Gibbs采样得到的可视层取值(解码或重构)，更新参数
        ph = self.encode(v0)
        ph_cd = self.encode(v_cd)
        self.weights += eta * (v0.T @ ph - v_cd.T @ ph)  # 更新连接权重参数
        self.bias_b += eta * np.mean(ph - ph_cd, axis=0)  # 更新隐藏层偏移量b
        self.bias_a += eta * np.mean(v0 - v_cd, axis=0)  # 更新可视层偏移量a
        return

    def train(self, data, max_step, max_cd=2, eta=0.1):
        # 训练主函数,采用对比散度算法(CD算法)更新参数
       
        assert data.shape[1] == self.n_visible, "输入数据维度与可视层神经元数目不相等"
        self.n_sample = data.shape[0]
        for i in range(max_step):
            v_cd = self.gibbs_sample(data, max_cd)
            self.update(data, v_cd, eta)
            error = np.sum((data - v_cd) ** 2) / self.n_sample / self.n_visible * 100
            if i == (max_step-1):  # 将重构后的样本与原始样本对比计算误差
                print("可视层(隐藏层)状态误差比例:{0}%".format(round(error, 2)))
    def predict(self, v):
        # 输入训练数据，预测隐藏层输出
        ph = self