deep learning:DBN源码 ----- C++

其代码原型来自于yusugomori，DBN其主要就是RBM的堆叠，在代码实现里需要区分好第几个RBM，其每个RBM的隐层是下一个RBM的输入层。以及微调的原理就是根据已有的标签数据和最后一个隐含层生成的输出层进行调整。初学小娃。

DBN.h

class DBN {

public:
    int N;
    int n_ins;
    int *hidden_layer_sizes;
    int n_outs;
    int n_layers;
    HiddenLayer **sigmoid_layers;
    RBM **rbm_layers;
    LogisticRegression *log_layer;
    DBN(int, int, int*, int, int);
    ~DBN();
    void pretrain(int*, double, int, int);
    void finetune(int*, int*, double, int);
    void predict(int*, double*);
};

HiddenLayer.h

class HiddenLayer {

public:
    int N;
    int n_in;
    int n_out;
    double **W;
    double *b;
    HiddenLayer(int, int, int, double**, double*);
    ~HiddenLayer();
    double output(int*, double*, double);
    void sample_h_given_v(int*, int*);
};

LogisticRegression.h

class LogisticRegression {

public:
    int N;
    int n_in;
    int n_out;
    double **W;
    double *b;
    LogisticRegression(int, int, int);
    ~LogisticRegression();
    void train(int*, int*, double);
    void softmax(double*);
    void predict(int*, double*);
};

RBM.h

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*);
};

DBN.cpp

#include <iostream>
#include <cmath>
#include "HiddenLayer.h"
#include <stdlib.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];

  //初始化为多少层，在本代码是有2个RBM组成，所以是2个
    for(int i=0; i<n_layers; i++)
    {
        if(i == 0)
        {
        input_size = n_ins;//第一层就是原始输入层
        }
        else
        {
            input_size = hidden_layer_sizes[i-1];//第二层为第一个隐含层作为下一层的输入层
        }

    //构造隐含层
    sigmoid_layers[i] = new HiddenLayer(N, input_size, hidden_layer_sizes[i], NULL, NULL);

    //构造RBM层
    rbm_layers[i] = new RBM(N, input_size, hidden_layer_sizes[i],sigmoid_layers[i]->W, sigmoid_layers[i]->b, NULL);
    }
     //根据上面的构造可以发现构造出的是2层RBM结构
    //初始化logistic,很显然起输入层是最后一层隐含层，输出层是整个模型的输出层
    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;
}

//数据处理，注意是两层RBM，所以需要分开处理，另外要得到最后一层隐藏层。
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++)
    {
        for(int epoch=0; epoch<epochs; epoch++)
        {
            for(int n=0; n<N; n++)
            {
                for(int m=0; m<n_ins; m++) train_X[m] = input[n * n_ins + m];
                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;
}
//微调就是根据标签数据的标签来进行微调，过程和LR类似。
void DBN::finetune(int *input, int *label, double lr, int epochs)
{
    int *layer_input;
    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++)
        {
            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];
            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);
        }
    }
    delete[] layer_input;
    delete[] train_X;
    delete[] train_Y;
}

void DBN::predict(int *x, double *y)
{
    double *layer_input;
    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];
    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 += 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()
{
    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);//获得h0

    for(int step=0; step<k; step++)
    {
        if(step == 0)
        {
            gibbs_hvh(ph_sample, nv_means, nv_samples, nh_means, nh_samples);//获得V1,h1
        }
        else
        {
            gibbs_hvh(nh_samples, nv_means, nv_samples, nh_means, nh_samples);
        }
    }

    //更新权值，双向偏移量。由于hinton提出的CD-K，可以知道其v0代表的是原始数据x
    //h0即ph_sigm_out，h0近似等于对v0下h的概率
    //v1即代表的是经过一次转换后的x，近似等于对h0下v的概率。
    //h1同理。CD-K主要就是求出这个三个数据，便能够很好的近似计算梯度。至于为什么我也不知道。
    for(int i=0; i<n_hidden; i++)
    {
        for(int j=0; j<n_visible; j++)
        {
      //可以根据权重公式发现，其实P（hi=1|v）代表的就是h0，p(hi=1|Vyk)和Vyk代表的就是h1和V1.
            W[i][j] += lr * (ph_mean[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}
    };
    //构造RBN包括构造多层，隐藏层，RBM及LR
    DBN dbn(train_N, n_ins, hidden_layer_sizes, n_outs, n_layers);
    //预处理过程
    dbn.pretrain(*train_X, pretrain_lr, k, pretraining_epochs);
    //微调
    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;
}