受限玻尔兹曼机的python参考实现

简介

众所周知，玻尔兹曼机好是好，但是太复杂了，所以在实际应用中不会使用。而受限玻尔兹曼机（Restricted Boltzmann machine, RBM）是它的一个带约束的版本，因此模型变得更加简单。受限玻尔兹曼机是一种生成式的机器学习模型，能够学习样本的概率分布。根据学到的概率分布，可生成符合分布的样本。例如可以用它学习人脸图片的概率分布，然后用它学到的概率分布来生成一张这个世界不存在的人脸。

总而言之，RBM是通过输入数据集学习概率分布的随机生成神经网络。它把网络中的节点分为两层：

• 可见层
• 隐藏层

每层有若干节点。
可见层$v=(v_1,v_2,...,v_i,...,v_n)$，隐藏层$h=(h_1,h_2,...,h_j,...h_m)$，层内的节点没有连接。层间结点两两相连。每一条连接都有一个权重$w_{ij}$。每个节点都是二值的随机变量 v i , h j ∈ { 0 ， 1 } v_i,h_j \in \{0，1\} vi​,hj​∈{0，1}。

RBM是一种随机的动力系统，因此用联合组态能量 表示系统的一种总体状态。定义如下：
$$E(v,h;\theta )=-\sum _{ij}{w}_{ij}{v}_i{h}_j-\sum _i{b}_i{v}_i-\sum _j{a}_j{h}_j$$
其中模型参数 θ = { w i j , a j , b i } \theta = \{w_{ij},a_j,b_i\} θ={wij​,aj​,bi​}。$a_j,b_i$分别是偏置单元。

那么根据联合组态能量，可定义$v,h$这两组随机向量的联合概率分布：
$$\begin{array}{rl}{p}_{\theta }(v,h)& =\frac{1}{z(\theta )}exp(-E(v,h;\theta ))\\ & =\frac{1}{z(\theta )}\prod _{ij}{e}^{w_{ij}{v}_i{h}_j}\prod _i{e}^{b_i{v}_i}\prod _j{e}^{a_j{h}_j}\end{array}$$
其中$Z(\theta )={\sum }_{v,h}exp(-E(v,h;\theta ))$。

RBM通过最大化似然函数来找到最优的参数$W,a,b$。

对于给定的训练样本 D = { v ^ ( 1 ) , v ^ ( 3 ) , . . . , v ^ ( N ) } D=\{\hat{v}^{(1)},\hat{v}^{(3)},...,\hat{v}^{(N)}\} D={v^(1),v^(3),...,v^(N)},其对数似然函数为：
$$L(D;\theta )=\frac{1}{N}\sum _{n=1}^N\log p({\hat{v}}^{(N)};\theta )$$
$p({\hat{v}}^{(N)})$可在联合分布的基础上求边缘分布得到。然后求偏导，使用梯度下降法求解参数。

$$\log p(v)=\log \sum _h\exp (-E(v,h))-\log \sum _{v^{\prime },h^{\prime }}\exp (-E(v^{\prime }{h}^{\prime }))$$
$$\frac{\partial \log p(v)}{\partial \theta }=...=E_{p(h|v)}\left[\frac{-\partial E(v,h)}{\partial \theta }\right]-E_{p(v^{\prime },h^{\prime })}\left[\frac{-\partial E(v^{\prime },h^{\prime })}{\partial \theta }\right]$$

具体地：
$$\begin{array}{rl}\frac{\partial \log p(v)}{\partial w_{ij}}& =E_{p(h|v)}\left(v_i{h}_j\right)-E_{p(v^{\prime },h^{\prime })}\left(v_i^{\prime }{h}_j^{\prime }\right)\\ \frac{\partial \log p(v)}{\partial a_i}& =E_{p(h|v)}\left(v_i\right)-E_{p(v^{\prime },h^{\prime })}\left(v_i^{\prime }\right)\\ \frac{\partial \log p(v)}{\partial b_i}& =E_{p(h|v)}\left(h_j\right)-E_{p(v^{\prime },h^{\prime })}\left(h_j^{\prime }\right)\end{array}$$

偏导中含有的期望很难计算。因此需要进行采样来估计。采样及优化过程使用CD-k算法（对比散度学习算法）。通常k取1，所以算法过程如下：

其中的$p(h=1|v)$和$p(v=1|h)$也可以很容易从联合概率密度函数中推导出来。详见《神经网络与深度学习》P299。

代码

本文对RBM代码进行了封装。对外暴露两个函数：

• train
• generate

前者用于从样本学习的一概率分布，后者用于从学到的分布上进行一次随机采样。完整代码 如下。

import numpy as np
from typing import *

class RBM:
    """
    受限玻尔兹曼机参考实现
    """
    def __init__(self, N:int, M:int):
        """
        初始化受限波尔兹曼机的参数
        :param N: 可见层神经元个数
        :param M: 隐藏层神经元个数
        """
        self.W = np.zeros((N, M))
        self.a = np.zeros((N,))
        self.b = np.zeros((M))

    def __sigmoid(self, x:np.ndarray)->np.ndarray:
        """
        非线性激活函数
        :param x: 输入的数组
        :return: 输出的数组
        """
        return 1 / (1 + np.exp(- (x)))

    def __sample(self, p:np.ndarray)->np.ndarray:
        """
        按一定概率采样一个0，1数组
        :param p: 概率向量
        :return: 0/1向量
        """
        values = []
        for i in range(len(p)):
            value = np.random.choice([0,1], p=[1 - p[i], p[i]])
            values.append(value)
        return np.array(values)


    def __sample_hidden(self, v:np.ndarray)->np.ndarray:
        """
        根据的可见状态采样一个隐藏状态，即按照p(h=1|v)的概率进行采样
        :param v: 可见状态向量
        :return: 隐藏状态向量
        """
        WT = np.transpose(self.W, [1,0])
        b = self.b
        return self.__sample(self.__sigmoid(WT.dot(v) + b))

    def __sample_visible(self, h:np.ndarray)->np.ndarray:
        """
        根据隐藏状态采样一个可见状态, 即按照p(v=1|h)的概率进行采样
        :param h: 隐藏状态向量
        :return: 可见状态向量
        """
        W = self.W
        a = self.a
        return self.__sample(self.__sigmoid(W.dot(h) + a))


    def train(self, samples:List[np.ndarray], leaning_rate:float=0.01, epoch:int = 3):
        """
        使用CD-1算法不断优化模型的参数
        :param samples: 训练样本,维度为 (batch-size, N)
        :param leaning_rate: 学习率
        :param epoch: 训练的回合数
        :return: 无
        """
        for i in range(epoch):
            for j in range(len(samples)):
                sample = samples[j]
                h = self.__sample_hidden(sample)
                v1 = self.__sample_visible(h)
                h1 = self.__sample_hidden(v1)
                self.W += leaning_rate * (
                        sample.reshape((len(sample), 1)).dot( h.reshape((1, len(h))) )
                        -
                        v1.reshape((len(v1),1)).dot( h1.reshape((1, len(h1))) )
                )
                self.a += leaning_rate * (sample - v1)
                self.b += leaning_rate * (h - h1)

    def generate(self)->np.ndarray:
        """
        随机生成一个符合模型分布的可观测样本
        :return: 样本
        """
        ramdom_h = [np.random.randint(0,2) for i in range(len(self.b))]
        return self.__sample_visible(np.array(ramdom_h))

实验结果

本文使用MNIST数据集优化模型参数，以使RBM模型能够学习到数据集的分布。
数据集下载地址 http://yann.lecun.com/exdb/mnist/
实验代码地址 https://gitee.com/clouddea/rbm
下面是实验代码。程序将手写数字的图片数据从文件中读取。原图片是值为0-255单通道的灰度图像，这里为了方便模型处理，将之转换成二值图像。二值图像数据作为RBM模型的可见层的状态。

每100张图像合成一批(即一个batch)送给模型训练一次。 训练结束后，令模型随机生成100张图片。

from RBM import  RBM
import numpy as np
import struct
import cv2
import matplotlib.pyplot as plt
from typing import *

def decode_idx3_ubyte(idx3_ubyte_file):
    # 这部分代码来自： https://blog.csdn.net/panrenlong/article/details/81736754
    # 因此略去
    return images

def image2vector(image:np.ndarray)->np.ndarray:
    image[image < 0] = 0
    image[image > 0] = 1
    # plt.imshow(image * 256)
    # plt.show()
    # plt.pause(0)
    return image.flatten()

def vector2image(vector:np.ndarray)->np.ndarray:
    vector = vector.reshape(28,28)
    vector[vector != 1] = 0
    vector[vector == 1] = 256
    return vector


if __name__ == '__main__':
    #实例化模型
    rbm = RBM(784, 100) # 784 = 28 * 28， 因为这里一张图像有784个像素
    #读取图片数据
    images = decode_idx3_ubyte('./train-images.idx3-ubyte')
    i = 0
    while i < len(images):
        # 构建一批数据
        images_list = images[i:i+100]
        images_batch:List[np.ndarray] = []
        for j in range(len(images_list)):
            image = images_list[j].astype(np.int)
            images_batch.append(image2vector(image))

        #学习一次
        rbm.train(images_batch)
        print("正在处理%s / %s"%(i, len(images)))
        i += 100
        if i >= 10000:
            break


    #模型自动生成一张图片
    for i in range(100):
        v = rbm.generate()
        plt.imshow(vector2image(v), 'gray')
        plt.savefig("./output/%s.png"%(i,))
        print("正在生成 %s.png"%(i,))

下面是上诉实验代码生成的图像。可以发现模型生成的图像与输入的图像具有相似的特征。即比较亮的像素集中在图像中部。部分图像能够看出是某个数字。这说明模型学习到了输入样本的分布特征。

下面几张是看起来比较明显的：
数字7：

数字6

总结

本文实现了RBM，并使用它学习了MNIST数据集的分布。从生成的图像来看模型的确能够学到样本的分布特征。但是生成图像依然不够明显，这可能是因为模型本身不够强大，也有可能训练得不够充分。