受限玻尔兹曼机代码实现（Python）

1、前言

上一篇受限玻尔兹曼机原理推导，我们讲了受限玻尔兹曼机的原理推导。本文，我们来看看受限玻尔兹曼机如何应用到推荐系统中。
数学基础：【概率论与数理统计知识复习-哔哩哔哩】

2、案例

由于我们仅仅讲了二值的受限玻尔兹曼机，所以本文也仅仅针对二值问题。受限玻尔兹曼机也常常用在推荐系统中。

如何推荐？一般情况下，用户给出它们所喜欢（厌恶）的电影（或其他东西），比如有5部电影，甲给了【0，1，1，0，1】，0表示不喜欢，1表示喜欢。我们通常将这些0，1值当作训练数据。也就是观测变量。

那么隐藏层是什么东西呢？我们没有明确的界定，因为在神经网络中，其是一种抽象出来的概念。如果你的五部电影种，三部是喜剧，两部是恐怖片。那么设置两个隐藏层，或许可以理解为对喜剧和恐怖片的喜好。

那么，当有一个新的用户来的时候，它对该五部电影给出了偏好评价，然后我们自然可以计算出它喜欢什么类型的电影；亦或者，我们再进行一次反向的，当用户给了它喜欢什么类型的电影，我们同样也可以给他推荐相关的电影了。

3、代码实现

import numpy as np

np.random.seed(2)
class RBM():
    def __init__(self,x_layer,h_layer):
        '''
        :param x_num: 可见层维度
        :param h_num: 隐藏层维度
        '''
        self.x_layer=x_layer #可见层的维度
        self.h_layer=h_layer #隐藏层的维度
        self.w=np.random.normal(0, 0.01, size=(self.x_layer, self.h_layer)) #从正态分布中随机采样w
        self.a=np.random.normal(0, 0.01, size=(self.h_layer,1)) #从正态分布中随机采样a
        self.b=np.random.normal(0, 0.01, size=(self.x_layer,1)) #从正态分布中随机采样b
        self.learning_rate=0.1 #学习率
    def train(self,x,K):
        '''
        :param x: 训练数据
        :param K: 使用k次吉布斯采样
        :return:
        '''
        x_num=x.shape[0] #样本的个数

        for _ in np.arange(10000): #梯度上升迭代10000次
            x0=x
            #################
            #CD-K吉布斯采样
            for _ in np.arange(K): #吉布斯采样K次

                P_h=self.sigmoid_Ph_x(x0) #从v0计算出P(h=1|v0)

                #从P(h=1|v0)采样出h0
                h0=np.random.binomial(1,p=P_h,size=(x_num,self.h_layer))
                #计算出P(v|h0)
                P_x=self.sigmoid_Px_h(h0)

                #采样出v
                x0=np.random.binomial(1,p=P_x,size=(x_num,self.x_layer))
            #################
            #真实数据的P(h=1|x)
            true_h =self.sigmoid_Ph_x(x)

            #采样数据的P(h=1|x)
            x_sample_h=self.sigmoid_Ph_x(x0)
            #w梯度
            w_GD=(x.T@true_h-x0.T@x_sample_h)/x_num

            #a梯度
            a_GD=np.mean(true_h-x_sample_h,axis=0).reshape(-1,1)

            #b梯度
            b_GD=np.mean(x-x0,axis=0).reshape(-1,1)

            #梯度下降
            self.w+=self.learning_rate*w_GD
            self.a+=self.learning_rate*a_GD
            self.b+=self.learning_rate*b_GD
    def sigmoid_Ph_x(self,x):
        '''
        计算P(h=1|x)
        :param x: 数据
        :return:
        '''
        H=x @ self.w+self.a.T
        result=1/(1+np.exp(-H))
        return result
    def sigmoid_Px_h(self,h):
        '''
        计算P(x=1|h)
        :param h:
        :return:
        '''
        H=(self.w @ h.T + self.b).T
        result=1/(1+np.exp(-H))
        return result

    def predict(self,x):
        '''
        预测
        :param x: 数据
        :return:
        '''
        P=self.sigmoid_Ph_x(x) #计算概率
        result=np.random.binomial(1,P) #返回概率最大所对应的值
if __name__ == '__main__':
    K=2 #吉布斯采样k次
    x=np.array([[1,1,1,0,1],
                [1,1,1,0,0],
                [1,0,0,1,1],
                [0,0,0,1,1],]) #训练数据
    rbm = RBM(5, 2) #初始化
    rbm.train(x,K) #训练
    new_x = np.array([[1,0,1,0,0],
                     [0,0,0,1,1]]) #测试数据
    rbm.predict(new_x) #预测