RBM（受限玻尔兹曼机）原理及代码

本文转自http://blog.csdn.net/zdy0_2004/article/details/45798223，虽然源贴也是转载的，转自http://www.cnblogs.com/xiaokangzi/p/4492466.html

本文在转载时，加入了本人的一些勘误和补充了缺失的代码，另外，还有一点未搞明白的是，下面所附的源码中使用的函数

binomial(1, mean[i])，目前尚未搞清楚其具体实现和源码，虽然翻了一通Gibbs Sampling的理论，可仍没搞出来这个函数的实现。如果有人明白，还望不吝赐教，在此先行谢过。

基于能量模型 (EBM)

基于能量模型将关联到感兴趣的变量每个配置的标量能量。学习修改的能量函数使他它的形状具有最好的性能。例如，我们想的得到最好的参量拥有较低的能量。

EBM的概率模型定义通过能量函数的概率分布，如下所示：

　　　　　　　　

规则化系数 Z 称为分区函数和物理系统的能量模型相似。

　　　　　　　　

一种基于能量模型可以学习通过随机梯度下降的方法处理负对数似然训练数据的。至于 logistic 回归分析我们将首次作为负对数似然定义对数似然然后损失函数。

　　　　　　　　

使用随机梯度  更新参数权值,  是模型中的各种参数 .

 

EBMs 的隐藏神经元

在很多情况下, 我们看不到部分的隐藏单元  , 或者我们要引入一些不可见的参量来增强模型的能力.所以我们考虑一些可见的神经元(依然表示为 ) 和 隐藏的部分 . 我们可以这样写我们的表达式:

                                                                                                                                          (2)

在这种情况下,公式类似于 (1), 我们引入符号， *能量, 定义如下:

                                                                                                                                                         (3)

也可以这么写,

数据的负极大似然梯度表示为.

                                                                                                                                 (4)

注意上面的梯度表示为两个部分，涉及到正面的部分和负面的部分.正面和负面的不表示等式中每部分的符号，而是表示对模型中概率密度的影响. 第一部分增加训练数据的概率 (通过降低相应的*能量), 第二部分降低模型确定下降梯度通常是很困难的, 因为他涉及到计算. 这无非在所有配置下的期望 (符合由模型生成的概率分布 ) !

 

第一步是计算估计固定数量的模型样本的期望. 用来表示负面部分梯度的表示为负粒子, 表示为 . 梯度可以谢伟:

                                                                                                                  (5)

我们想 根据 取样元素 of   (例如. 我们可以做 蒙特卡罗方法). 通过上面的公式, 我们几乎可以使用随机粒子算法来学习EBM模型. 唯一缺少的就是如何提取这些负粒子T . 统计学上有许多抽样方法, 马尔可夫链蒙特卡罗方法特别适合用于模型如受限玻尔兹曼机 (RBM), 一种特殊的 EBM.

 

受限玻尔兹曼机 (RBM)

玻尔兹曼机（BMS）是一种特殊的对数线性马尔可夫随机场（MRF）的形式，即，其能量函数在其*参数的线性空间里。使他们强大到足以代表复杂的分布，我们考虑到一些变量是没有观察到（他们称为隐藏）。通过更多的隐藏变量（也称为隐藏的单位），我们可以增加的玻尔兹曼机的建模能力（BM）。受限玻尔兹曼机进一步限制BMS中那些可见-可见和隐藏-隐藏的连接。下面是一个RBM的图形描述。

RBM能量方程  定义为 :

                                  (6)

 表示隐藏单元和可见单元连接的权重， ,  是可见层和隐藏层的偏置.

转化成下面的*能量公式:

由于RBM的特殊结构, 可见和隐藏单元 是条件独立的. 利用这个特性, 我们可以得出:

 

二值化的RBMs 

在通常的情况下使用二值化单元 (和), 我们从公式. (6) and (2)得到, 概率版本的常用神经元激活函数:

                                        (7)

                                       (8)

二值RBM的*能量可以更简单的表示为:

                           (9)

 

二值RBM的参数更新方程

结合等式 (5) 和 (9), 我们得到下面的对数似然梯度方程:

                                                      (10)

 

RBM中的采样

样本  可以通过运行Markov chain收敛得到，使用gibbs采样作为转移操作.

N随机变量的Gibbs采样的联合概率可以通过一系列的采样得到 ，其中 包含  个 中其他的随机参数但不包括.

对于RBM，包含一组可见或者不可见单元，由于他们是条件独立的你可以执行块Gibbs采样。在这种背景下对可见单元同时采样来得到隐藏单元的值. 同样的可以对隐藏单元同时采样，得到可见单元的值。 一步Markov chainA 由下面公式得到:

其中 表示是指在Markov chain第n步 的所有隐藏单元. 意思是, 例如,  是根据概率随机选择为1（或者0）, 相似的,  是根据概率随机选择为1（或者0）.

这可以用下图说明：

 

当 , 样本 保证处于真实样本下 .

理论下，每个参数的获取都必须运行这样一个链知道收敛，不用说这样做代价是十分昂贵的。因此，大家为RBM设计了不少算法，为了在学习过程中得到在概率分布下的有效样本。

Contrastive Divergence（中文有时称为对比散度）

(CD-k)

CD使用两个技巧来加速采样过程 :

• 由于我们最终想要得到  (真实的处于数据分布的样本), 我们把马尔科夫链初始化为一个训练样本（即，我们要使一个分布靠近p，那么我们的链应该已经收敛到最后的分布p）
• CD不需要等待链收敛，样本结果k步gibbs采样得到，在实践中，k=1工作的出奇的好。

Persistent CD

Persistent CD [Tieleman08] 使用另一种方法近似下的样本，他依赖于一个拥有持续性的马尔科夫链， (i.e.,不为每个可见样本重新计算链 ). 对于每一个参量更新，我们简单的运行k-步链生成新的样本。 链的状态需要保存用来计算以后步骤的参量更新.

一般的直觉是，如果参量更新对于链的混合速率足够小，马尔科夫链应该可以察觉到模型的改变。

 

 

下面我们需要用代码来实现RBM（使用c语言）

//return a random number (double) with a Uniform Distribution, the value is between min and max

double random_uniform_distribution(double min,double max)

{

    return min+(max-min)*rand()/(RAND_MAX+1.0);    //RAND_MAX is 0x7FFF

    //According to the time-based seed which was initialised by srand(), a random number is created.

    //If no seed was initialised by srand(), rand() would use the seed initialised by srand(1) automatically and a pseudo-random number would be created.

    //Do NOT initialise seed within this sub-function, because everytime when this function is called, the seed would be initialised newly. CPU runs repaidly, there may be no time-difference among several successive calling for the function, and the seeds initialised are more likely the same, so the phenomenon that the random numbers obtained are all the same may occur.

}

double sigmoid(double x)

{

    double exp_value;

    double return_value;

    /*** Exponential calculation ***/

    exp_value = exp((double) -x);

    /*** Final sigmoid value ***/

    return_value = 1 / (1 + exp_value);

    return return_value;

}

void RBM__construct(RBM* ptRBM,int N, int n_visible,int n_hidden, \

                    double **W,double *hbias, double *vbias) {

    int i, j;

    double a =1.0 / n_visible;

    ptRBM->N = N;

    ptRBM->n_visible = n_visible;

    ptRBM->n_hidden = n_hidden;

    if(W ==NULL) {

        ptRBM->W = (double **)malloc(sizeof(double*) * n_hidden);

        ptRBM->W[0] = (double *)malloc(sizeof(double) * n_visible * n_hidden);

        for(i=0; i<n_hidden; i++) ptRBM->W[i] = ptRBM->W[0] + i * n_visible;

        for(i=0; i<n_hidden; i++) {

            for(j=0; j<n_visible; j++) {

                ptRBM->W[i][j] =random_uniform_distribution(-a, a);

            }

        }

    } else {

        ptRBM->W = W;

    }

    if(hbias ==NULL) {

        ptRBM->hbias = (double *)malloc(sizeof(double) * n_hidden);

        for(i=0; i<n_hidden; i++) ptRBM->hbias[i] =0;

    } else {

        ptRBM->hbias = hbias;

    }

    if(vbias ==NULL) {

        ptRBM->vbias = (double *)malloc(sizeof(double) * n_visible);

        for(i=0; i<n_visible; i++) ptRBM->vbias[i] =0;

    } else {

        ptRBM->vbias = vbias;

    }

}

//由可见层得到隐藏样本

double RBM_propdown(RBM* ptRBM,int *h, int i,double b) {

    int j;

    double pre_sigmoid_activation =0.0;

    for(j=0; j<ptRBM->n_hidden; j++) {

        pre_sigmoid_activation += ptRBM->W[j][i] * h[j];

    }

    pre_sigmoid_activation += b;

    returnsigmoid(pre_sigmoid_activation);

}

void RBM_sample_v_given_h(RBM* ptRBM,int *h0_sample, double *mean,int *sample) {

    int i;

    for(i=0; i<ptRBM->n_visible; i++) {

        mean[i] = RBM_propdown(ptRBM, h0_sample, i, ptRBM->vbias[i]);

        sample[i] = binomial(1, mean[i]);

    }

}

//由隐藏样本得到可见样本

double RBM_propup(RBM* ptRBM,int *v, double *w,double b) {

    int j;

    double pre_sigmoid_activation =0.0;

    for(j=0; j<ptRBM->n_visible; j++) {

        pre_sigmoid_activation += w[j] * v[j];

    }

    pre_sigmoid_activation += b;

    returnsigmoid(pre_sigmoid_activation);

}

void RBM_sample_h_given_v(RBM* ptRBM,int *v0_sample, double *mean,int *sample) {

    int i;

    for(i=0; i<ptRBM->n_hidden; i++) {

        mean[i] = RBM_propup(ptRBM, v0_sample, ptRBM->W[i], ptRBM->hbias[i]);

        sample[i] = binomial(1, mean[i]);

    }

}

//Gibbs采样：

void RBM_gibbs_hvh(RBM* ptRBM,int *h0_sample, double *nv_means,int *nv_samples, \

                   double *nh_means,int *nh_samples) {

    RBM_sample_v_given_h(ptRBM, h0_sample, nv_means, nv_samples);

    RBM_sample_h_given_v(ptRBM, nv_samples, nh_means, nh_samples);

}

//运行CD-K并且更新权值：

void RBM_contrastive_divergence(RBM* ptRBM,int *input, double lr,int k) {

    int i, j, step;

    double *ph_mean = (double *)malloc(sizeof(double) * ptRBM->n_hidden);

    int *ph_sample = (int *)malloc(sizeof(int) * ptRBM->n_hidden);

    double *nv_means = (double *)malloc(sizeof(double) * ptRBM->n_visible);

    int *nv_samples = (int *)malloc(sizeof(int) * ptRBM->n_visible);

    double *nh_means = (double *)malloc(sizeof(double) * ptRBM->n_hidden);

    int *nh_samples = (int *)malloc(sizeof(int) * ptRBM->n_hidden);

    /* CD-k */

    RBM_sample_h_given_v(ptRBM, input, ph_mean, ph_sample);

    for(step=0; step<k; step++) {

        if(step ==0) {

            RBM_gibbs_hvh(ptRBM, ph_sample, nv_means, nv_samples, nh_means, nh_samples);

        } else {

            RBM_gibbs_hvh(ptRBM, nh_samples, nv_means, nv_samples, nh_means, nh_samples);

        }

    }

    for(i=0; i<ptRBM->n_hidden; i++) {

        for(j=0; j<ptRBM->n_visible; j++) {

            // this->W[i][j] += lr * (ph_sample[i] * input[j] - nh_means[i] * nv_samples[j]) / this->N;

            ptRBM->W[i][j] += lr * (ph_mean[i] * input[j] - nh_means[i] * nv_samples[j]) / ptRBM->N;

        }

        ptRBM->hbias[i] += lr * (ph_sample[i] - nh_means[i]) / ptRBM->N;

    }

    for(i=0; i<ptRBM->n_visible; i++) {

        ptRBM->vbias[i] += lr * (input[i] - nv_samples[i]) / ptRBM->N;

    }

    free(ph_mean);

    free(ph_sample);

    free(nv_means);

    free(nv_samples);

    free(nh_means);

    free(nh_samples);

}

//如何重建样本：就是  v->h->v

void RBM_reconstruct(RBM* ptRBM,int *v, double *reconstructed_v) {

    int i, j;

    double *h = (double *)malloc(sizeof(double) * ptRBM->n_hidden);

    double pre_sigmoid_activation;

    for(i=0; i<ptRBM->n_hidden; i++) {

        h[i] = RBM_propup(ptRBM, v, ptRBM->W[i], ptRBM->hbias[i]);

    }

    for(i=0; i<ptRBM->n_visible; i++) {

        pre_sigmoid_activation = 0.0;

        for(j=0; j<ptRBM->n_hidden; j++) {

            pre_sigmoid_activation += ptRBM->W[j][i] * h[j];

        }

        pre_sigmoid_activation += ptRBM->vbias[i];

        reconstructed_v[i] = sigmoid(pre_sigmoid_activation);

    }

    

    free(h);

}

void RBM__destruct(RBM* ptRBM) {

    if(ptRBM->W !=NULL) {

        if(ptRBM->W[0]) {

            free(ptRBM->W[0]);

        }

        free(ptRBM->W);

    }

    if(ptRBM->hbias !=NULL) {

        free(ptRBM->hbias);

    }

    if(ptRBM->vbias !=NULL) {

        free(ptRBM->vbias);

    }

}

void test_rbm(void) {

    srand(0);

    int i, j, epoch;

    double learning_rate =0.1;

    int training_epochs =1000;

    int k =1;

    int train_N =6;

    int test_N =2;

    int n_visible =6;

    int n_hidden =3;

    // training data

    int train_X[6][6] = {

        {1,1, 1,0, 0,0},

        {1,0, 1,0, 0,0},

        {1,1, 1,0, 0,0},

        {0,0, 1,1, 1,0},

        {0,0, 1,0, 1,0},

        {0,0, 1,1, 1,0}

    };

    // construct RBM

    RBM rbm;

    RBM__construct(&rbm, train_N, n_visible, n_hidden,NULL, NULL,NULL);

    // train

    for(epoch=0; epoch<training_epochs; epoch++) {

        for(i=0; i<train_N; i++) {

            RBM_contrastive_divergence(&rbm, train_X[i], learning_rate, k);

        }

    }

    // test data

    int test_X[2][6] = {

        {1,1, 0,0, 0,0},

        {0,0, 0,1, 1,0}

    };

    double reconstructed_X[2][6];

    // test

    for(i=0; i<test_N; i++) {

        RBM_reconstruct(&rbm, test_X[i], reconstructed_X[i]);

        for(j=0; j<n_visible; j++) {

            printf("%.5f ", reconstructed_X[i][j]);

        }

        printf("\n");

    }

    

    

    // destruct RBM

    RBM__destruct(&rbm);

}