源自:http://blog.youkuaiyun.com/marvin521/article/details/8886971
经过三年的狂刷理论,觉得是时候停下来做些有用的东西了,因此决定开博把他们写下来,一是为了整理学过的理论,二是监督自己并和大家分享。先从DeepLearning谈起吧,因为这个有一定的实用性(大家口头传的“和钱靠的很近”),国内各个大牛也都谈了不少,我尽量从其他方面解释一下。
DeepLearning算是多伦多大学Geoffery hinton教授第二春吧,第一春就是传统神经网络,由于传统的多层感知机很容易陷入局部最小,直接用反向传播算法(Back Propagation)求取的分类效果并不如意,原因一是特征是手工的,二就是局部最小问题。而DeepLearning引入了概率图模型里的生成模型,他可以直接自动的从训练集里提取所需要的特征,典型的模型为有限制玻尔兹曼机(Restricted Boltzmann Machines,简称RBM),自动提取的特征解决了手工特征考虑不周的因素,而且很好的初始化了神经网络权重,接着可以采用反向传播算法进行分类,实验得出了很好的效果。 因此DeepLearning被喻为下一代神经网络。今天的话题就来讨论下RBM:
再说RBM之前,我们先提一下基于能量的模型(Engery based model),能量方法来源于热动力学,分子在高温中运动剧烈,能够克服局部约束(分子之间的一些物理约束,比如键值吸引力等),在逐步降到低温时,分子最终会排列出有规律的结构,此时也是低能量状态。受此启发,早期的模拟退火算法就是在高温中试图跳出局部最小。随机场作为物理模型之一,也引入了此方法。在马尔科夫随机场(MRF)中能量模型主要扮演着两个作用:一、全局解的度量(目标函数);二、能量最小时的解(各种变量对应的配置)为目标解,能否把最优解嵌入到能量函数中至关重要,决定着我们具体问题求解的好坏。统计模式识别主要工作之一就是捕获变量之间的相关性,同样能量模型也要捕获变量之间的相关性,变量之间的相关程度决定了能量的高低。把变量的相关关系用图表示出来,并引入概率测度方式就够成了概率图模型的能量模型,其实实际中也可以不用概率表示,比如立体匹配中直接用两个像素点的像素差作为能量,所有像素对之间的能量和最小时的配置即为目标解。RBM作为一种概率图模型,引入概率就是为了方便采样,因为在CD(contrastive divergence)算法中采样部分扮演着模拟求解梯度的角色。
先来看一下RBM模型定义(图一)和能量的定义(公式一):
(公式一)
其中v表示实际样本节点,以下叫可视节点,h表示隐藏节点数据,theta={W,b,a},W表示网络权重,b表示可视节点偏移,a表示隐藏节点偏移。这些参数是需要我们求解的,一般开始时我们随机初始化这些参数,对于16*16的二值图像,可视节点有256个,每个节点可能的状态有两个(0,1);隐藏节点的数目由自己制定,每个节点的状态也有两个(0,1)。当每个节点都有一个指定的状态时,此时的RBM有一个对应的能量。此时所有节点对应某个状态的组合成为一个配置。穷举所有配置对应的能量之和就构成了归一化常量Z(配分函数-partition function)。每个配置能量除以常量Z即为当前配置的概率。在一般计算时,为了避免数值计算问题,我们常常引入极家族函数(对数函数和指数函数),那么模型概率定义为(公式二):
由于RBM是一种生成模型(generative model),我们的目标使得概率p(x)最大时对应的参数theta,常规思路应该使用最大似然法,在求最大似然函数的极值时,我们需要用到梯度法,梯度推导如(图二)所示(大家不要被公式吓倒,CD算法不用它,=、=,只是证明求梯度有点难度):
公式虽然这么写,但是实际中计算归一化分割函数Z时不可取的,因为组合状态是NP-complete的。但是这个公式给出了Contrastive divergence算法的原型(公式三):
第一项(数据相关期望项):利用可视节点v1乘上权重,然后用sigmoid函数计算出隐藏节点h每个节点的概率,通过此概率对隐藏节点采样,得出隐藏节点状态h1,然后用能量公式计算出所有节点的期望能量(<v1*h1>)作为第一项。
第二项(模型期望项):根据第一项中隐藏节点h1乘上权重,然后用sigmoid函数计算出可视节点v2的概率,在此通过采样得出v2 状态,同样用v2计算出h2的状态,用能量公式计算出所有节点的期望能量(<v2*h2>)。套用上面公式,就可以求出deltaW,
有了类梯度,那么就可以做权重更新咯,迭代就行了,b,a类似求解。
这是用CD算法模拟梯度,关键在于使用什么采样方法,基于MCMC的各种采样方法(如Gibbs,MH,拒绝采样等)请酌情使用,深度理解采样方法有助于我们理解为什么CD可以模拟梯度,用hinton的话来说:这是个惊人的事实。我概括了这么一句话来给我们一个直观的理解:数据期望分布是我们整个数据的分布中心,可以看成最小点,模型期望是采样出来的,具有一定的随机性,分布在中心点周围,二者相减正好可以模拟梯度。其实求解参数theta的方法不止这一种,比如在Russlan的论文中也给出另外一种方法:用随机逼近过程来估计模型期望,用变分方法来估计数据相关期望等。碍于作者水平有限,就不讲了。
另外对DBN这个生成模型(generative model)对参数求导时,得出的结果也和CD公式相似,DBN的训练方法也可以用类似方法训练,具体可参考russlan的代码。
DBM训练方法其实就是层层通过对RBM分别训练。对于russlan论文中提到的其他衍生应用和技巧大家可以根据自己所需详细阅读。比如估计分割函数(Z)是图模型里一个重要的基石部分;autoencoders的神奇作用其实是模拟了神经科学里神经分布式描述,相当于找到了低维流形;深度网络的微调部分(fine-tune),其实通过上述方法做了无监督预训练后,结合标签部分用BP等方法做了权重微更新,值得一提的识别可以用生成模型(generative model),也可以是判别模型(discriminative model),针对具体应用大家可以试着采用不同方法。
我想我要说的基本说完了,整个深度学习除了需要概率图模型和神经网络基础支撑外,也用到了大量的trick。
哦,对了,hinton也说过,神经网络几乎可以模拟现在绝大多数的机器学习方法,包括SVM,Adaboost,当然也包括概率图模型的那一套(其中有结构化学习),这有点把图模型收编的味道,作者视野有限,不敢不认同,也不敢认同。就讲到这里吧,作者写的轻浮一些,公式都是从别人paper上扣下来的,
而且作者水平有限,可能有些会说错,希望大家多多指正。
最后附加一点:概率图模型的能量和概率之间的关系可以参阅gibbs measure和Hammersley–Clifford theorem,他们俩在能量模型到概率之间的转换上扮演着很重要的角色。还有为什么最大化P(x)或者p(v)可以使得能量最小?这个问题一开始我没搞明白,后来一个做事严谨网友(北流浪子)的花了好长时间找到答案,并且无私的与俺分享:对(图二)中的公式(5.7)两边同时取对数得到:lnp(v)=-FreeEnergy(v)-lnZ,这样当p(v)最大时,自由能FreeEnergy(v)最小,而自由能最小决定了系统达到统计平衡,详细可以参考:Gibbs free energy。
参考文献:
[1] Learning Deep Generative models.Ruslan Salakhutdinov
[2] Learning Deep Architectures for AI.Yoshua Bengio
[3] A Tutorial on Energy-Based Models. yann lecun
下面贴出RBM C++版本的代码,一些大牛写的,结合上篇博文来加深大家对RBM理论的理解。。。
RBM类定义声明:
- class RBM {
- public:
- int N;
- int n_visible;
- int n_hidden;
- double **W;
- double *hbias;
- double *vbias;
- RBM(int, int, int, double**, double*, double*);
- ~RBM();
- void contrastive_divergence(int*, double, int);
- void sample_h_given_v(int*, double*, int*);
- void sample_v_given_h(int*, double*, int*);
- double propup(int*, double*, double);
- double propdown(int*, int, double);
- void gibbs_hvh(int*, double*, int*, double*, int*);
- void reconstruct(int*, double*);
- };
从上面声明中可以很直观的看出和上篇文章公式符号正好完美对应。下面是代码实现部分:
- #include <iostream>
- #include <math.h>
- #include "RBM.h"
- using namespace std;
- double uniform(double min, double max) {
- return rand() / (RAND_MAX + 1.0) * (max - min) + min;
- }
- int binomial(int n, double p) {
- if(p < 0 || p > 1) return 0;
- int c = 0;
- double r;
- for(int i=0; i<n; i++) {
- r = rand() / (RAND_MAX + 1.0);
- if (r < p) c++;
- }
- return c;
- }
- double sigmoid(double x) {
- return 1.0 / (1.0 + exp(-x));
- }
- RBM::RBM(int size, int n_v, int n_h, double **w, double *hb, double *vb) {
- N = size;
- n_visible = n_v;
- n_hidden = n_h;
- if(w == NULL) {
- W = new double*[n_hidden];
- for(int i=0; i<n_hidden; i++) W[i] = new double[n_visible];
- double a = 1.0 / n_visible;
- for(int i=0; i<n_hidden; i++) {
- for(int j=0; j<n_visible; j++) {
- W[i][j] = uniform(-a, a);
- }
- }
- } else {
- W = w;
- }
- if(hb == NULL) {
- hbias = new double[n_hidden];
- for(int i=0; i<n_hidden; i++) hbias[i] = 0;
- } else {
- hbias = hb;
- }
- if(vb == NULL) {
- vbias = new double[n_visible];
- for(int i=0; i<n_visible; i++) vbias[i] = 0;
- } else {
- vbias = vb;
- }
- }
- RBM::~RBM() {
- for(int i=0; i<n_hidden; i++) delete[] W[i];
- delete[] W;
- delete[] hbias;
- delete[] vbias;
- }
- void RBM::contrastive_divergence(int *input, double lr, int k) {
- double *ph_mean = new double[n_hidden];
- int *ph_sample = new int[n_hidden];
- double *nv_means = new double[n_visible];
- int *nv_samples = new int[n_visible];
- double *nh_means = new double[n_hidden];
- int *nh_samples = new int[n_hidden];
- /* CD-k */
- sample_h_given_v(input, ph_mean, ph_sample);
- for(int step=0; step<k; step++) {
- if(step == 0) {
- gibbs_hvh(ph_sample, nv_means, nv_samples, nh_means, nh_samples);
- } else {
- gibbs_hvh(nh_samples, nv_means, nv_samples, nh_means, nh_samples);
- }
- }
- for(int i=0; i<n_hidden; i++) {
- for(int j=0; j<n_visible; j++) {
- W[i][j] += lr * (ph_sample[i] * input[j] - nh_means[i] * nv_samples[j]) / N;
- }
- hbias[i] += lr * (ph_sample[i] - nh_means[i]) / N;
- }
- for(int i=0; i<n_visible; i++) {
- vbias[i] += lr * (input[i] - nv_samples[i]) / N;
- }
- delete[] ph_mean;
- delete[] ph_sample;
- delete[] nv_means;
- delete[] nv_samples;
- delete[] nh_means;
- delete[] nh_samples;
- }
- void RBM::sample_h_given_v(int *v0_sample, double *mean, int *sample) {
- for(int i=0; i<n_hidden; i++) {
- mean[i] = propup(v0_sample, W[i], hbias[i]);
- sample[i] = binomial(1, mean[i]);
- }
- }
- void RBM::sample_v_given_h(int *h0_sample, double *mean, int *sample) {
- for(int i=0; i<n_visible; i++) {
- mean[i] = propdown(h0_sample, i, vbias[i]);
- sample[i] = binomial(1, mean[i]);
- }
- }
- double RBM::propup(int *v, double *w, double b) {
- double pre_sigmoid_activation = 0.0;
- for(int j=0; j<n_visible; j++) {
- pre_sigmoid_activation += w[j] * v[j];
- }
- pre_sigmoid_activation += b;
- return sigmoid(pre_sigmoid_activation);
- }
- double RBM::propdown(int *h, int i, double b) {
- double pre_sigmoid_activation = 0.0;
- for(int j=0; j<n_hidden; j++) {
- pre_sigmoid_activation += W[j][i] * h[j];
- }
- pre_sigmoid_activation += b;
- return sigmoid(pre_sigmoid_activation);
- }
- void RBM::gibbs_hvh(int *h0_sample, double *nv_means, int *nv_samples, \
- double *nh_means, int *nh_samples) {
- sample_v_given_h(h0_sample, nv_means, nv_samples);
- sample_h_given_v(nv_samples, nh_means, nh_samples);
- }
- void RBM::reconstruct(int *v, double *reconstructed_v) {
- double *h = new double[n_hidden];
- double pre_sigmoid_activation;
- for(int i=0; i<n_hidden; i++) {
- h[i] = propup(v, W[i], hbias[i]);
- }
- for(int i=0; i<n_visible; i++) {
- pre_sigmoid_activation = 0.0;
- for(int j=0; j<n_hidden; j++) {
- pre_sigmoid_activation += W[j][i] * h[j];
- }
- pre_sigmoid_activation += vbias[i];
- reconstructed_v[i] = sigmoid(pre_sigmoid_activation);
- }
- delete[] h;
- }
- void test_rbm() {
- srand(0);
- 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(train_N, n_visible, n_hidden, NULL, NULL, NULL);
- // train
- for(int epoch=0; epoch<training_epochs; epoch++) {
- for(int i=0; i<train_N; i++) {
- rbm.contrastive_divergence(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(int i=0; i<test_N; i++) {
- rbm.reconstruct(test_X[i], reconstructed_X[i]);
- for(int j=0; j<n_visible; j++) {
- printf("%.5f ", reconstructed_X[i][j]);
- }
- cout << endl;
- }
- }
- int main() {
- test_rbm();
- return 0;
- }
干脆把运行结果也贴出来,给那些终极极品思考者提供一些方便
0.98472 0.67248 0.99120 0.01000 0.01311 0.01020
0.01021 0.00720 0.99525 0.65553 0.98403 0.00497
考虑到大家有可能对深度学习的识别有点模糊,因此决定写一个短博客,简单介绍下如何识别,结合本系列的第一篇博文提到的深度学习之所以叫深度,其中之一的原因是多层RBM模仿了人脑多层神经元对输入数据进行层层预处理(值得一提的是并不是每层都是RBM,DBN就是个例外),即深层次的数据拟合,多个RBM连接起来构成DBM(deep boltzmann machines),每层RBM的节点数自己指定,这需要一些经验,DBN和DBM在训练上没有区别,都是逐层使用CD算法训练,也叫贪心预训练算法,DBN和DBM的区别如(图一)所示。
图一
对于二者都使用同一个算法来训练,看起来毫无区别,但是DBM有一个优势,由于RBM是无向的,这就决定了无论给定可视节点还是隐藏节点,各个节点都是独立的,可由图模型的马尔科夫性看出。作为无向图的DBM天生具有一些优秀的基因,比如当人看到一个外观性质,知道它是什么物体,同样你告诉他物体名字,他可以知道物体的外观应该是什么样子。这种互相推理的关系正好可以用无向图来表示。这种优势也顺理成章的延伸出了autoencoder(大家所谓的自编码神经网络)和栈式神经网络,最终输出的少量节点是可以推理(重建)出原来样本,也起到了降维的作用,无形中也找到了特征(编码),autoencoder的效果如图二所示。但是DBN中有些层是有向的,就不具有这种优势。
图二
二者逐层预训练后,结合样本标签,使用BP算法进行权重微调,说白了就是在预训练后的权重基础上使用BP算法进行训练,这样得出的权重更好些。。。
下面贴出部分DBN代码,大家可以看出总体思路是按照构建DBN网络(刚构建后的每层的权重是随机生成的,从代码也能看出),贪心层层预训练,权重微调,预测(识别)这个步骤来的。另外代码中softmax其实是多变量的逻辑回归函数,注意我发的下面的代码中权重微调使用的是逻辑回归,不是BP:
- #include <iostream>
- #include <math.h>
- #include "HiddenLayer.h"
- #include "RBM.h"
- #include "LogisticRegression.h"
- #include "DBN.h"
- using namespace std;
- double uniform(double min, double max) {
- return rand() / (RAND_MAX + 1.0) * (max - min) + min;
- }
- int binomial(int n, double p) {
- if(p < 0 || p > 1) return 0;
- int c = 0;
- double r;
- for(int i=0; i<n; i++) {
- r = rand() / (RAND_MAX + 1.0);
- if (r < p) c++;
- }
- return c;
- }
- double sigmoid(double x) {
- return 1.0 / (1.0 + exp(-x));
- }
- // DBN
- DBN::DBN(int size, int n_i, int *hls, int n_o, int n_l) {
- int input_size;
- N = size;
- n_ins = n_i;
- hidden_layer_sizes = hls;
- n_outs = n_o;
- n_layers = n_l;
- sigmoid_layers = new HiddenLayer*[n_layers];
- rbm_layers = new RBM*[n_layers];
- // construct multi-layer
- for(int i=0; i<n_layers; i++) {
- if(i == 0) {
- input_size = n_ins;
- } else {
- input_size = hidden_layer_sizes[i-1];
- }
- // construct sigmoid_layer
- sigmoid_layers[i] = new HiddenLayer(N, input_size, hidden_layer_sizes[i], NULL, NULL);
- // construct rbm_layer
- rbm_layers[i] = new RBM(N, input_size, hidden_layer_sizes[i],\
- sigmoid_layers[i]->W, sigmoid_layers[i]->b, NULL);
- }
- // layer for output using LogisticRegression
- log_layer = new LogisticRegression(N, hidden_layer_sizes[n_layers-1], n_outs);
- }
- DBN::~DBN() {
- delete log_layer;
- for(int i=0; i<n_layers; i++) {
- delete sigmoid_layers[i];
- delete rbm_layers[i];
- }
- delete[] sigmoid_layers;
- delete[] rbm_layers;
- }
- void DBN::pretrain(int *input, double lr, int k, int epochs) {
- int *layer_input;
- int prev_layer_input_size;
- int *prev_layer_input;
- int *train_X = new int[n_ins];
- for(int i=0; i<n_layers; i++) { // layer-wise
- for(int epoch=0; epoch<epochs; epoch++) { // training epochs
- for(int n=0; n<N; n++) { // input x1...xN
- // initial input
- for(int m=0; m<n_ins; m++) train_X[m] = input[n * n_ins + m];
- // layer input
- for(int l=0; l<=i; l++) {
- if(l == 0) {
- layer_input = new int[n_ins];
- for(int j=0; j<n_ins; j++) layer_input[j] = train_X[j];
- } else {
- if(l == 1) prev_layer_input_size = n_ins;
- else prev_layer_input_size = hidden_layer_sizes[l-2];
- prev_layer_input = new int[prev_layer_input_size];
- for(int j=0; j<prev_layer_input_size; j++) prev_layer_input[j] = layer_input[j];
- delete[] layer_input;
- layer_input = new int[hidden_layer_sizes[l-1]];
- sigmoid_layers[l-1]->sample_h_given_v(prev_layer_input, layer_input);
- delete[] prev_layer_input;
- }
- }
- rbm_layers[i]->contrastive_divergence(layer_input, lr, k);
- }
- }
- }
- delete[] train_X;
- delete[] layer_input;
- }
- void DBN::finetune(int *input, int *label, double lr, int epochs) {
- int *layer_input;
- // int prev_layer_input_size;
- int *prev_layer_input;
- int *train_X = new int[n_ins];
- int *train_Y = new int[n_outs];
- for(int epoch=0; epoch<epochs; epoch++) {
- for(int n=0; n<N; n++) { // input x1...xN
- // initial input
- for(int m=0; m<n_ins; m++) train_X[m] = input[n * n_ins + m];
- for(int m=0; m<n_outs; m++) train_Y[m] = label[n * n_outs + m];
- // layer input
- for(int i=0; i<n_layers; i++) {
- if(i == 0) {
- prev_layer_input = new int[n_ins];
- for(int j=0; j<n_ins; j++) prev_layer_input[j] = train_X[j];
- } else {
- prev_layer_input = new int[hidden_layer_sizes[i-1]];
- for(int j=0; j<hidden_layer_sizes[i-1]; j++) prev_layer_input[j] = layer_input[j];
- delete[] layer_input;
- }
- layer_input = new int[hidden_layer_sizes[i]];
- sigmoid_layers[i]->sample_h_given_v(prev_layer_input, layer_input);
- delete[] prev_layer_input;
- }
- log_layer->train(layer_input, train_Y, lr);
- }
- // lr *= 0.95;
- }
- delete[] layer_input;
- delete[] train_X;
- delete[] train_Y;
- }
- void DBN::predict(int *x, double *y) {
- double *layer_input;
- // int prev_layer_input_size;
- double *prev_layer_input;
- double linear_output;
- prev_layer_input = new double[n_ins];
- for(int j=0; j<n_ins; j++) prev_layer_input[j] = x[j];
- // layer activation
- for(int i=0; i<n_layers; i++) {
- layer_input = new double[sigmoid_layers[i]->n_out];
- for(int k=0; k<sigmoid_layers[i]->n_out; k++) {
- // linear_output = 0.0; //原代码中删除此句
- for(int j=0; j<sigmoid_layers[i]->n_in; j++) {
- linear_output = 0.0; //原代码中添加此句
- linear_output += sigmoid_layers[i]->W[k][j] * prev_layer_input[j];
- }
- linear_output += sigmoid_layers[i]->b[k];
- layer_input[k] = sigmoid(linear_output);
- }
- delete[] prev_layer_input;
- if(i < n_layers-1) {
- prev_layer_input = new double[sigmoid_layers[i]->n_out];
- for(int j=0; j<sigmoid_layers[i]->n_out; j++) prev_layer_input[j] = layer_input[j];
- delete[] layer_input;
- }
- }
- for(int i=0; i<log_layer->n_out; i++) {
- y[i] = 0;
- for(int j=0; j<log_layer->n_in; j++) {
- y[i] += log_layer->W[i][j] * layer_input[j];
- }
- y[i] += log_layer->b[i];
- }
- log_layer->softmax(y);
- delete[] layer_input;
- }
- // HiddenLayer
- HiddenLayer::HiddenLayer(int size, int in, int out, double **w, double *bp) {
- N = size;
- n_in = in;
- n_out = out;
- if(w == NULL) {
- W = new double*[n_out];
- for(int i=0; i<n_out; i++) W[i] = new double[n_in];
- double a = 1.0 / n_in;
- for(int i=0; i<n_out; i++) {
- for(int j=0; j<n_in; j++) {
- W[i][j] = uniform(-a, a);
- }
- }
- } else {
- W = w;
- }
- if(bp == NULL) {
- b = new double[n_out];
- } else {
- b = bp;
- }
- }
- HiddenLayer::~HiddenLayer() {
- for(int i=0; i<n_out; i++) delete W[i];
- delete[] W;
- delete[] b;
- }
- double HiddenLayer::output(int *input, double *w, double b) {
- double linear_output = 0.0;
- for(int j=0; j<n_in; j++) {
- linear_output += w[j] * input[j];
- }
- linear_output += b;
- return sigmoid(linear_output);
- }
- void HiddenLayer::sample_h_given_v(int *input, int *sample) {
- for(int i=0; i<n_out; i++) {
- sample[i] = binomial(1, output(input, W[i], b[i]));
- }
- }
- // RBM
- RBM::RBM(int size, int n_v, int n_h, double **w, double *hb, double *vb) {
- N = size;
- n_visible = n_v;
- n_hidden = n_h;
- if(w == NULL) {
- W = new double*[n_hidden];
- for(int i=0; i<n_hidden; i++) W[i] = new double[n_visible];
- double a = 1.0 / n_visible;
- for(int i=0; i<n_hidden; i++) {
- for(int j=0; j<n_visible; j++) {
- W[i][j] = uniform(-a, a);
- }
- }
- } else {
- W = w;
- }
- if(hb == NULL) {
- hbias = new double[n_hidden];
- for(int i=0; i<n_hidden; i++) hbias[i] = 0;
- } else {
- hbias = hb;
- }
- if(vb == NULL) {
- vbias = new double[n_visible];
- for(int i=0; i<n_visible; i++) vbias[i] = 0;
- } else {
- vbias = vb;
- }
- }
- RBM::~RBM() {
- // for(int i=0; i<n_hidden; i++) delete[] W[i];
- // delete[] W;
- // delete[] hbias;
- delete[] vbias;
- }
- void RBM::contrastive_divergence(int *input, double lr, int k) {
- double *ph_mean = new double[n_hidden];
- int *ph_sample = new int[n_hidden];
- double *nv_means = new double[n_visible];
- int *nv_samples = new int[n_visible];
- double *nh_means = new double[n_hidden];
- int *nh_samples = new int[n_hidden];
- /* CD-k */
- sample_h_given_v(input, ph_mean, ph_sample);
- for(int step=0; step<k; step++) {
- if(step == 0) {
- gibbs_hvh(ph_sample, nv_means, nv_samples, nh_means, nh_samples);
- } else {
- gibbs_hvh(nh_samples, nv_means, nv_samples, nh_means, nh_samples);
- }
- }
- for(int i=0; i<n_hidden; i++) {
- for(int j=0; j<n_visible; j++) {
- W[i][j] += lr * (ph_sample[i] * input[j] - nh_means[i] * nv_samples[j]) / N;
- }
- hbias[i] += lr * (ph_sample[i] - nh_means[i]) / N;
- }
- for(int i=0; i<n_visible; i++) {
- vbias[i] += lr * (input[i] - nv_samples[i]) / N;
- }
- delete[] ph_mean;
- delete[] ph_sample;
- delete[] nv_means;
- delete[] nv_samples;
- delete[] nh_means;
- delete[] nh_samples;
- }
- void RBM::sample_h_given_v(int *v0_sample, double *mean, int *sample) {
- for(int i=0; i<n_hidden; i++) {
- mean[i] = propup(v0_sample, W[i], hbias[i]);
- sample[i] = binomial(1, mean[i]);
- }
- }
- void RBM::sample_v_given_h(int *h0_sample, double *mean, int *sample) {
- for(int i=0; i<n_visible; i++) {
- mean[i] = propdown(h0_sample, i, vbias[i]);
- sample[i] = binomial(1, mean[i]);
- }
- }
- double RBM::propup(int *v, double *w, double b) {
- double pre_sigmoid_activation = 0.0;
- for(int j=0; j<n_visible; j++) {
- pre_sigmoid_activation += w[j] * v[j];
- }
- pre_sigmoid_activation += b;
- return sigmoid(pre_sigmoid_activation);
- }
- double RBM::propdown(int *h, int i, double b) {
- double pre_sigmoid_activation = 0.0;
- for(int j=0; j<n_hidden; j++) {
- pre_sigmoid_activation += W[j][i] * h[j];
- }
- pre_sigmoid_activation += b;
- return sigmoid(pre_sigmoid_activation);
- }
- void RBM::gibbs_hvh(int *h0_sample, double *nv_means, int *nv_samples, \
- double *nh_means, int *nh_samples) {
- sample_v_given_h(h0_sample, nv_means, nv_samples);
- sample_h_given_v(nv_samples, nh_means, nh_samples);
- }
- void RBM::reconstruct(int *v, double *reconstructed_v) {
- double *h = new double[n_hidden];
- double pre_sigmoid_activation;
- for(int i=0; i<n_hidden; i++) {
- h[i] = propup(v, W[i], hbias[i]);
- }
- for(int i=0; i<n_visible; i++) {
- pre_sigmoid_activation = 0.0;
- for(int j=0; j<n_hidden; j++) {
- pre_sigmoid_activation += W[j][i] * h[j];
- }
- pre_sigmoid_activation += vbias[i];
- reconstructed_v[i] = sigmoid(pre_sigmoid_activation);
- }
- delete[] h;
- }
- // LogisticRegression
- LogisticRegression::LogisticRegression(int size, int in, int out) {
- N = size;
- n_in = in;
- n_out = out;
- W = new double*[n_out];
- for(int i=0; i<n_out; i++) W[i] = new double[n_in];
- b = new double[n_out];
- for(int i=0; i<n_out; i++) {
- for(int j=0; j<n_in; j++) {
- W[i][j] = 0;
- }
- b[i] = 0;
- }
- }
- LogisticRegression::~LogisticRegression() {
- for(int i=0; i<n_out; i++) delete[] W[i];
- delete[] W;
- delete[] b;
- }
- void LogisticRegression::train(int *x, int *y, double lr) {
- double *p_y_given_x = new double[n_out];
- double *dy = new double[n_out];
- for(int i=0; i<n_out; i++) {
- p_y_given_x[i] = 0;
- for(int j=0; j<n_in; j++) {
- p_y_given_x[i] += W[i][j] * x[j];
- }
- p_y_given_x[i] += b[i];
- }
- softmax(p_y_given_x);
- for(int i=0; i<n_out; i++) {
- dy[i] = y[i] - p_y_given_x[i];
- for(int j=0; j<n_in; j++) {
- W[i][j] += lr * dy[i] * x[j] / N;
- }
- b[i] += lr * dy[i] / N;
- }
- delete[] p_y_given_x;
- delete[] dy;
- }
- void LogisticRegression::softmax(double *x) {
- double max = 0.0;
- double sum = 0.0;
- for(int i=0; i<n_out; i++) if(max < x[i]) max = x[i];
- for(int i=0; i<n_out; i++) {
- x[i] = exp(x[i] - max);
- sum += x[i];
- }
- for(int i=0; i<n_out; i++) x[i] /= sum;
- }
- void LogisticRegression::predict(int *x, double *y) {
- for(int i=0; i<n_out; i++) {
- y[i] = 0;
- for(int j=0; j<n_in; j++) {
- y[i] += W[i][j] * x[j];
- }
- y[i] += b[i];
- }
- softmax(y);
- }
- void test_dbn() {
- srand(0);
- double pretrain_lr = 0.1;
- int pretraining_epochs = 1000;
- int k = 1;
- double finetune_lr = 0.1;
- int finetune_epochs = 500;
- int train_N = 6;
- int test_N = 3;
- int n_ins = 6;
- int n_outs = 2;
- int hidden_layer_sizes[] = {3, 3};
- int n_layers = sizeof(hidden_layer_sizes) / sizeof(hidden_layer_sizes[0]);
- // 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, 1, 0, 0},
- {0, 0, 1, 1, 1, 0}
- };
- int train_Y[6][2] = {
- {1, 0},
- {1, 0},
- {1, 0},
- {0, 1},
- {0, 1},
- {0, 1}
- };
- // construct DBN
- DBN dbn(train_N, n_ins, hidden_layer_sizes, n_outs, n_layers);
- // pretrain
- dbn.pretrain(*train_X, pretrain_lr, k, pretraining_epochs);
- // finetune
- dbn.finetune(*train_X, *train_Y, finetune_lr, finetune_epochs);
- // test data
- int test_X[3][6] = {
- {1, 1, 0, 0, 0, 0},
- {0, 0, 0, 1, 1, 0},
- {1, 1, 1, 1, 1, 0}
- };
- double test_Y[3][2];
- // test
- for(int i=0; i<test_N; i++) {
- dbn.predict(test_X[i], test_Y[i]);
- for(int j=0; j<n_outs; j++) {
- cout << test_Y[i][j] << " ";
- }
- cout << endl;
- }
- }
- int main() {
- test_dbn();
- return 0;
- }
- <pre></pre>
- <p>程序运行结果,是个二维
今天就来讨论下deep learning怎么来处理real valued data。对于图像来说,二值图像毕竟是少数,更多的还是实值图像。对于这样的情况,RBM已经无法很好的处理它们,因此需要改进它,对于了解计算机视觉的人而言,想必高斯混合背景模型大家已不陌生,高斯混合模型可以很好的对实值图像建模,OpenCV中早就用高斯混合背景模型来分割物体。接下来要引出的高斯有限制玻尔兹曼机(Gaussian Restricted Boltzmann machine-GRBM)和高斯混合模型有一定的等价性。
在上一节中,作者提到过对于RBM这个无向图模型而言,无论是给定隐藏节点还是可视节点,剩下的那层内节点之间互相独立。那么对于RBM这个图模型描述的联合分布可以分解成(公式一)的形式:
(公式一)
其中P(h)可看成高斯混合系数,如果P(v|h)是高斯分布,那么这个条件下的RBM就和高斯混合分布具有等价性。因此GRBM合理登场 ,下面开始进入GRBM正题,首先当然要假定P(v|h)服从高斯分布,不然没法进展下去咯,先来看下GRBM的能量公式定义(公式二):
从能量公式中可以很明显的看出有高斯分布的影子,当然还有其他类型的能量公式,关键看你做什么假设,假设不一样,能量公式也不一样,假设的好就会产生一篇好文章哦。其中参数要求的为:。求取算法和此系列文章一方法类似,也是用CD算法,但是对于GRBM的CD算法,需要多啰嗦几句,首先我们把博文一中的求取数据期望项和模型期望项的过程改成官方称呼:分别对应为正阶段(positive phase)和负阶段(negative phase)。
正阶段的目的就是找到给定数据的情况下,隐藏节点的配置(暂称它为"好配置"),并且使得能量降低。负阶段的目的当然就是抬升"好配置"周围的配置的能量,可能有点拗口,推荐看下博文一中的参考文献:Energy based model tutorial,可以把它总结成一句话:把和实际样本有关的配置能量降低,把fancy的样本的能量抬高,这样得到的权重是对样本的最佳描述,当然样本越多越好,所以大数据这么火。有些人会有疑问,样本多会不会导致过拟合?其实生成模型的魅力之一就在这,样本越多越好。正阶段中隐藏节点h的状态求取仍然用sigmoid函数得出概率,然后采样;但是负阶段中可视节点V的状态的概率就不是用sigmoid函数了,而是用高斯函数,二者的公式分别如(公式三)和(公式四)表示:
(公式三,f为sigmoid函数)
(公式四,高斯分布函数)
对于参数的求取,过程只是多了sigma而已,如果你熟悉高斯混合模型,自然对sigma的实际作用也有些直观认识,庆幸的是sigma仍然也可以在CD算法中与其他参数一样同时求出来。下面给出GRBM中CD算法的部分代码(matlab),其中invfstdInc为sigma,不要用用这段代码,而是用它解释自己对GRBM实现细节的疑问,这段代码中,CD算法引入了momentum(冲量)技巧,可以让算法收敛快些:
- %%%%%%%%% START POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- data = batchdata(:,:,batch); %nxd
- pos_hidprobs = 1./(1 + exp(-(data./Fstd)*vhW - Ones*hb)); %p(h_j =1|data)
- pos_hidstates = pos_hidprobs > rand( size(pos_hidprobs) );
- pos_prods = (data./Fstd)'* pos_hidprobs;
- pos_hid_act = sum(pos_hidprobs);
- pos_vis_act = sum(data)./(fstd.^2); %see notes on this derivation
- %%%%%%%%% END OF POSITIVE PHASE %%%%%%%%%
- for iterCD = 1:params.nCD
- %%%%%%%%% START NEGATIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- negdataprobs = pos_hidstates*vhW'.*Fstd+Ones*vb;
- negdata = negdataprobs + randn(n, d).*Fstd;
- neg_hidprobs = 1./(1 + exp(-(negdata./Fstd)*vhW - Ones*hb )); %updating hidden nodes again
- pos_hidstates = neg_hidprobs > rand( size(neg_hidprobs) );
- end %end CD iterations
- neg_prods = (negdata./Fstd)'*neg_hidprobs;
- neg_hid_act = sum(neg_hidprobs);
- neg_vis_act = sum(negdata)./(fstd.^2); %see notes for details
- %%%%%%%%% END OF NEGATIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- errsum = errsum + sum(sum( (data-negdata).^2 ));
- if epoch > params.init_final_momen_iter,
- momentum=params.final_momen;
- else
- momentum=params.init_momen;
- end
- %%%%%%%%% UPDATE WEIGHTS AND BIASES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- vhWInc = momentum*vhWInc + r/n*(pos_prods-neg_prods) - r*params.wtcost*vhW;
- vbInc = momentum*vbInc + (r/n)*(pos_vis_act-neg_vis_act);
- hbInc = momentum*hbInc + (r/n)*(pos_hid_act-neg_hid_act);
- invfstd_grad = sum(2*data.*(Ones*vb-data/2)./Fstd,1) + sum(data' .* (vhW*pos_hidprobs') ,2)';
- invfstd_grad = invfstd_grad - ( sum(2*negdata.*(Ones*vb-negdata/2)./Fstd,1) + ...
- sum( negdata'.*(vhW*neg_hidprobs') ,2 )' );
- invfstdInc = momentum*invfstdInc + std_rate(epoch)/n*invfstd_grad;
- if params.SPARSE == 1 %nair's paper on 3D object recognition
- %update q
- if batch==1 && epoch == 1
- q = mean(pos_hidprobs);
- else
- q_prev = q;
- q = 0.9*q_prev+0.1*mean(pos_hidprobs);
- end
- p = params.sparse_p;
- grad = 0.1*params.sparse_lambda/n*sum(pos_hidprobs.*(1-pos_hidprobs)).*(p-q)./(q.*(1-q));
- gradW =0.1*params.sparse_lambda/n*(data'./Fstd'*(pos_hidprobs.*(1-pos_hidprobs))).*repmat((p-q)./(q.*(1-q)), d,1);
- hbInc = hbInc + r*grad;
- vhWInc = vhWInc + r*gradW;
- end
- ptot = ptot+mean(pos_hidprobs(:));
- vhW = vhW + vhWInc;
- vb = vb + vbInc;
- hb = hb + hbInc;
- invfstd = 1./fstd;
- invfstd = invfstd + invfstdInc;
- fstd = 1./invfstd;
- fstd = max(fstd, 0.005); %have a lower bound!
参考文献:
Guassian-Bernoulli Deep Boltzmann Machine. KyungHyun Cho
结构决定功能和协同处理:鲁棒有限制玻尔兹曼机(RoBM)
上一篇博文中提到,GRBM训练样本越多越好,样本多蕴含的分布更具有一般性,这对其他模型也适用。但是实际样本数据中往往有大量的噪声,这或多或少的影响了GRBM的性能,工业界一般都设有清洗数据的岗位,用人脑去除噪声数据。试想:人脑为什么具有如此强大的抗噪和容错能力?其实生命科学中有一句经典的总结:结构决定功能。不同的结构的网络具有不同的功能,下面要引入的鲁棒玻尔兹曼机(Robust boltzmann machine-RoBM)就是一个很好的例证。RoBM由GRBM和RBM组成,确切的说还有一个成分:选择门。对于这个模型我不想说它怎么训练的咯,因为GRBM和RBM的训练方法前面都说过了,我只想强调下“结构决定功能”怎么体现在概率图模型上。RoBM的模型如(图一)所示:
(图一)
在给定V-tilde(tilde就是V头上的波浪号哈,表示输入节点像素)时,S和V是独立的,同样G和H是独立的,对于GRBM和RBM的节点独立前面也已经讨论过。由于这些独立性再结合对于我们问题有利的假设(如高斯分布),能量函数定义就出来了,(对于能量函数的定义,也可以用概率图模型的clque(团)来解释,其实也是变量间的独立关系决定)。RoBM的能量函数如(公式一)所示:
(公式一)
上面公式的第一行就是牵扯到变量S,V和V-tilde的选择门项,这项其实就决定选择门的功能,控制信号的输入,比如当S等于0时,随你V和V-tilde怎么变,这块能量不会太高;当S等于1时,gamma的平方控制着V和V-tilde的差异程度,如果二者差异太大,这块能量就会偏高。
第二行就是关于噪声S的RBM能量项(和RBM的能量公式是一样的),它控制着模型和噪声的关系,S是一个指示变量,表示该样本是否“应该”有噪声(“是否应该”这是个配置,有系统来完成)。
第三行当然就是干净数据V的GRBM能量项了(和GRBM的能量公式是一样的),他控制着模型和干净数据的关系。
第四行就是噪声的分布,假设噪声成高斯分布。

参考文献:
Robust Boltzmann Machines for Recognition and Denoising.Yi chuanTang,Ruslan,Geoffery Hinton