BP神經網絡

起源：線性神經網絡與單層感知器

古老的線性神經網絡，使用的是單層Rosenblatt感知器。該感知器模型已經不再使用，但是你可以看到它的改良版：Logistic回歸。

 

可以看到這個網絡，輸入->加權->映射->計算分類誤差->迭代修改W、b，其實和數學上的回歸擬合別無二致。

 

Logistic回歸對該模型進行了改良：

線性神經網絡（回歸）使用的LMS（最小均方）的數學原理其實可由最大似然估計+假設誤差概率模型得到。（詳見Andrew Ng視頻）

在二類分類（誤差非0即1）情況下，適用於連續型數據的最小均方顯然不是很好的cost函數，會引起梯度過大。

仿照線性回歸假設誤差服從正態分布建立概率模型，Logistic回歸假設誤差服從二項分布建立概率模型。

Logistic函數的（0~1連續特性）在這里充當着，由輸入評估概率的角色，而不是像下面的BP網絡一樣，起的是高維空間非線性識別作用。

該手法同樣在RBM限制玻爾茲曼機中使用。

 

實際上，這兩種模型的起源都是最小二乘法的線性回歸。不同的是，早期的解決線性回歸使用的矩陣解方程組，求得參數。

而基於梯度下降使目標函數收斂的數學方法，在計算神經科學領域，就變成神經網絡了。

 

 

Part I ：BP網絡的結構與工作方式

 

BP網絡中使用隱層（HideLayer）設定，目的是通過全連接的網絡+非線性Sigmoid函數，瘋狂地通過特征空間映射來區分非線性數據。

注意，BP網絡的非線性處理能力主要來自於隱層結構，而不單單是非線性激活函數。

所以BP網絡的這種結構稱為多層感知器（MLP，Multi-layer Perceptron）

從數學上是很難解釋原理的。（與之相對的是SVM支持向量機，最大化間隔數學原理讓你竟無言以對）。

隱層層數和每層的神經元個數決定着網絡復雜度，已有數據表明(單、雙隱層), 每層神經元個數與樣本個數相近，網絡效率最高。

輸出層很特殊，如果需要多類分類的話，需要自己設計輸出方式。常用的是通過0/1二進制編碼來設計。

如00、01、10、11表示不同的類，這樣的方式需要log2(K)+1的神經元，盡可能減輕整個網絡的壓力，畢竟神經網絡用的是全連接。

BP網絡的工作方式分為兩部分：

①FP（Front Propagation）前向傳播：輸入->線性加權傳至隱層->在隱層Sigmoid並輸出->線性加權傳至輸入層

                                 ->輸入層使用Sigmoid(或線性函數)處理輸入，進行分類，計算誤差並累加到總誤差(目標函數)

②BP（Back Propagation） 反向傳播：從輸入層開始，由當前處理單條數據的誤差，通過梯度法更新W_ij、B_ij。

                                 由於網絡的特殊性，此時W_mi、B_mi的更新依賴於W_ij、B_ij，只能反向更新。

BP網絡訓練數據有兩種方法：

①單樣本串行<類似隨機梯度算法>：按順序/隨機輸入每個樣本，每次迭代只對一個樣本執行FP、BP。直至單個樣本誤差收斂，退出迭代。

②批樣本並行<類似批梯度算法>：按順序輸入每個樣本，每次迭代對每個按順序執行FP、BP，累計總誤差。

                                               執行完畢之后，算一次迭代，繼續從第一個樣本開始，進行第二次迭代。直至總誤差收斂，退出迭代。

由於串行每次迭代只用了一個樣本，因而總迭代次數應該是並行的M倍，否則誤差很大。

批樣本並行計算可以直觀看到總誤差，一般用來調參數，確認收斂情況。

而單樣本串行計算，則更多的是在調完參數后，觀察是否提升正確率。

 

 

Part II：BP過程的公式推導

定義ui,vi,u_j,v_j,分別是隱層、輸入層的I/O。

Logistic-Sigmoid函數的導數：$S^{'}(x)=S(x)(1-S(x))$

輸出層誤差：$e_{j}=d_{j}-v_{J}^{j}$,其中d是真值，v是預測值，由於輸出層自行設計，所以分類->真值之間需要加工處理。

輸出層總誤差（LMS目標函數）：$e=\frac{1}{2}\sum_{j=1}^{J}e_{j}^2$

輸出層第j個神經元的輸出：$v_{J}^{j}=S(u_{J}^{j})$

輸出層第j個神經元的輸入：$u_{J}^{j}=\sum_{i=1}^{I}(W_{ij}\cdot v_{I}^{i}+b_{ij})$

Wij的梯度：$\frac{\partial e}{\partial W_{ij}}=\frac{\partial e}{\partial e_{j}}\cdot\frac{\partial e_{j}}{\partial v_{J}^{j}}\cdot \frac{\partial v_{J}^{j}}{\partial u_{J}^{j}}\cdot \frac{\partial u_{I}^{i}}{\partial W_{ij}}==-e_j\cdot S'(u_{J}^{j})\cdot v_{I}^{i}$

求導使用的是鏈式法則（最好復習一遍高數）,它的鏈式：$e->e_{j}->v_{J}^{j}->u_{J}^{j}->W_{ij}$

其中$\frac{\partial e}{\partial e_{j}}=e_j ,\quad \frac{\partial e_{j}}{\partial v_{J}^{j}}=-1,\quad \frac{\partial v_{J}^{j}}{\partial u_{J}^{j}}=S'(u_{J}^{j}),\quad\frac{\partial u_{I}^{i}}{\partial W_{ij}}=v_{I}^{i}$

定義局部梯度:$\delta _{J}^{j}=-\frac{\partial e}{\partial e_{j}}\cdot\frac{\partial e_{j}}{\partial v_{J}^{j}}\cdot \frac{\partial v_{J}^{j}}{\partial u_{J}^{j}}=e_j\cdot S'(u_{J}^{j})$ （$\delta _{J}^{j}$將在W_{mi更新中使用）。}

於是：$W_{ij}^{new}=W_{ij}^{old}+\alpha\cdot \delta _{J}^{j}\cdot v_{I}^{i}$      $b_{ij}^{new}=b_{ij}^{old}+\alpha\cdot \delta _{J}^{j}$

按照同樣的方式，有$W_{mi}^{new}=W_{mi}^{old}+\alpha\cdot \delta _{I}^{i}\cdot M_{m}$

$\delta_{I}^{i}=-\frac{\partial e}{\partial u_{I}^{i}}$的推導比較有趣，其鏈式：$(e->e_{j}->v_{J}^{j}->u_{J}^{j}->v_{I}^{i})->u_{I}^{i}$

注意由於是全連接，所以單I神經元，與后一層全部J神經元有關，借用$\delta _{j}$的導數式，括號部分=$\sum_{j=1}^{J}W_{ij}\cdot \delta _{j}$

結合上面的推導，有：$\delta_{I}^{i}=\sum_{j=1}^{J}W_{ij}\cdot \delta _{j}\cdot \frac{\partial v_{I}^{i}}{\partial u_{I}^{i}}=\sum_{j=1}^{J}W_{ij}\cdot \delta _{j}\cdot S'(u_{I}^{i})$

這樣，即可通過先計算$\delta _{I}$ 、$\delta _{J}$這兩個局部梯度的向量，然后拼成W的完整梯度，進而使用梯度法。

 

Part III：參數調整與動量BP優化

   BP網絡非常吃參數，調整很復雜。因而推薦取100個樣本，先進行小規模調參。

   由於Sigmoid函數的有效定義域大概是[-3,3]，因而，首先對數據進行縮放，控制在[-1,1]內最佳。

   梯度法的參數調整一直是個麻煩。步長$\alpha$需要根據輸入數據的特征大小調整，如果特征數值過大或過小，都會導致Sigmoid函數爆掉而導致無法迭代收斂。

  在梯度方面，可以使用動量BP方法。動量BP法引入動量因子$\lambda$ $(0<\lambda<1)$，通常取值0.1~0.8

原更新量變成：$\Delta W_{ij}=\alpha\cdot \delta _{J}^{j}\cdot v_{I}^{i}\quad=>\;(1-\lambda)(\alpha\cdot \delta _{J}^{j}\cdot v_{I}^{i})+\lambda\Delta W_{ij}^{old}$

動量因子的使用，適當的考慮了前次更新：

①若前后兩次梯度方向相同，由於梯度值隨着目標函數的減小而減小，因而上一次的較大的梯度值會加大本次梯度混合值。

②若前后兩次梯度方向相反，表明可能在兩個位置有極值，上一次的梯度值可以抵消本次的部分梯度值，減小更新量。

動量因子的選取要看情況，通常先取0，然后，逐步增加。過大的動量因子，會導致抵消過大，無法收斂。

關於W參數初始值：已有論文表示，應當初始化隨機這個范圍的值$[-4*\frac{\sqrt{6}}{\sqrt{LayerInput+LayerOut}},4*\frac{\sqrt{6}}{\sqrt{LayerInput+LayerOut}}]$

這個和$\alpha$一樣，得根據數據情況，自行調節，過大或過小都會導致初始迭代失敗。

注意：BP網絡的W、b初始化非常重要，如果像Logistic回歸那樣直接設0,  最終迭代出來的網絡就和屎一樣。（比如下面的異或問題）

 

Part IV：測試與代碼

#include "cstdio"
#include "fstream"
#include "iostream"
#include "sstream"
#include "vector"
#include "math.h"
#include "stdlib.h"
using namespace std;
#define Dim dataSet[0].feature.size()
#define M dataSet.size()
#define alpha 0.6
#define delta 0.0000001
#define gamma 0.8
struct Data
{
    vector<double> feature;
    int y;
    Data(vector<double> feature,int y):feature(feature),y(y) {}
};
int HideLayerNum,OutputLayerNum,now_data=0;
vector<Data> dataSet,testSet;
vector<double> u_i,v_i,u_j,v_j;
vector<double> delta_j,delta_i,pdelta_j,pdelta_i;
vector< vector<double> > W_m_i,W_i_j,B_m_i,B_i_j,pW_i_j,pW_m_i;
double random(int f_in,int f_out)
{
    double ret1=rand()%((int)(sqrt(6)/sqrt(f_in+f_out)*100));
    double ret2=rand()%((int)(sqrt(6)/sqrt(f_in+f_out)*100));
    ret1/=100,ret2/=100;
    return ret1-ret2;
}
void read()
{
    ifstream fin("in.txt");
    string line;
    double fea;int cls;
    while(getline(fin,line))
    {
        stringstream sin(line);
        vector<double> feature;
        while(sin>>fea) feature.push_back(fea);
        cls=feature.back();feature.pop_back();
        dataSet.push_back(Data(feature,cls));
    }
    HideLayerNum=M-1; //隱層神經元個數=樣本數-1
    OutputLayerNum=1; //二類分類，一個輸出即可
   for(int i=0;i<Dim;i++) //初始化權W與偏置B
    {
        W_m_i.push_back(vector<double>(HideLayerNum,random(Dim,HideLayerNum)));
        pW_m_i.push_back(vector<double>(HideLayerNum,random(Dim,HideLayerNum)));
        B_m_i.push_back(vector<double>(HideLayerNum,random(Dim,HideLayerNum)));
    }
    for(int i=0;i<HideLayerNum;i++)
    {
       W_i_j.push_back(vector<double>(OutputLayerNum,random(HideLayerNum,OutputLayerNum)));
       pW_i_j.push_back(vector<double>(OutputLayerNum,random(HideLayerNum,OutputLayerNum)));
       B_i_j.push_back(vector<double>(OutputLayerNum,random(HideLayerNum,OutputLayerNum)));
    }
}
double sigmoid(double x) {return exp(x)/(1+exp(x));}
void buildInputLayer() //build inputLayer->hideLayer
{
    u_i.clear(); //re-calc
    v_i.clear();
    for(int i=0;i<HideLayerNum;i++)
    {
        double ret=0.0;
        for(int m=0;m<Dim;m++) ret+=(W_m_i[m][i]*dataSet[now_data].feature[m]+B_m_i[m][i]);
        u_i.push_back(ret);
        v_i.push_back(sigmoid(ret));
    }
}
double buildHideLayer() //build hideLayer->OutputLayer
{
    double error=0.0;
    u_j.clear();
    v_j.clear();
    for(int j=0;j<OutputLayerNum;j++)
    {
        double ret=0.0;
        for(int i=0;i<HideLayerNum;i++) ret+=(W_i_j[i][j]*v_i[i]+B_i_j[i][j]);
        u_j.push_back(ret);
        v_j.push_back(sigmoid(ret));
        error+=(dataSet[now_data].y-sigmoid(ret))*(dataSet[now_data].y-sigmoid(ret));
    }
    return error;
}
double FP()
{
    buildInputLayer();
    return buildHideLayer();
}
void BP()
{
    delta_i.clear();
    delta_j.clear();
    for(int j=0;j<OutputLayerNum;j++)  //calc delta_j=error*sigmoid'(u_j)=error*v_j(1-v_j)
    {
        double error=0.0;
        error+=(dataSet[now_data].y-v_j[j]); //all error;
        delta_j.push_back(v_j[j]*(1-v_j[j])*error);
    }
    for(int i=0;i<HideLayerNum;i++) //calc delta_i=Σ(delta_j*W_i_j)*sigmoid'(u_i)=Σ(delta_j*W_i_j)*v_i(1-v_i)
    {
        double ret=0.0;
        for(int j=0;j<OutputLayerNum;j++) ret+=W_i_j[i][j]*delta_j[j];
        delta_i.push_back(v_i[i]*(1-v_i[i])*ret);
    }
    for(int i=0; i<HideLayerNum; i++) //update W_i_j
        for(int j=0; j<OutputLayerNum; j++)
            {
                    double Delta=alpha*delta_j[j]*v_i[i];
                    W_i_j[i][j]+=(Delta*(1-gamma)+gamma*pW_i_j[i][j]);
                    pW_i_j[i][j]=(Delta*(1-gamma)+gamma*pW_i_j[i][j]);
                    B_i_j[i][j]+=alpha*delta_j[j];
            }
    for(int m=0; m<Dim; m++) //update W_m_i
        for(int i=0; i<HideLayerNum; i++)
            {
                double Delta=alpha*delta_i[i]*dataSet[now_data].feature[m];
                W_m_i[m][i]+=(Delta*(1-gamma)+gamma*pW_m_i[m][i]);
                pW_m_i[m][i]=(Delta*(1-gamma)+gamma*pW_m_i[m][i]);
                B_m_i[m][i]+=alpha*delta_i[i];
            }
}
void classify()
{
    /*
    ifstream fin("testdata.txt");
    dataSet.clear();
    double fea;int cls;
    string line;
    while(getline(fin,line))
    {
        stringstream sin(line);
        vector<double> feature;
        while(sin>>fea) feature.push_back(fea);
        cls=feature.back();feature.pop_back();
        dataSet.push_back(Data(feature,cls));
    }*/
    now_data=0;
    for(int s=0;s<dataSet.size();s++)
    {
        buildInputLayer();
        double ret=0.0;
        for(int j=0;j<OutputLayerNum;j++)
          for(int i=0;i<HideLayerNum;i++)
            ret+=(W_i_j[i][j]*v_i[i]+B_i_j[i][j]);
        printf("test%d: origin:%d test:%lf\n",now_data++,dataSet[s].y,sigmoid(ret));
    }
}
double iterProcess()
{
    double err=0.0;
    //random
    /*
    now_data=rand()%M;
    err=FP();
    BP();*/
    //batch
    now_data=0;
    for(int i=0;i<M;i++)
    {
        err+=FP();
        BP();
        now_data++;
    }
    return err;
}
int main()
{
    double oldLw,newLw;
    int iter=0;
    read();
    oldLw=iterProcess();
    newLw=iterProcess();
    while(fabs(oldLw-newLw)>delta)
    {
        oldLw=newLw;
        newLw=iterProcess();
        iter++;
        if(iter%1000==0) printf("iter:%d %lf\n",iter,newLw);
    }
    cout<<iter<<endl<<endl;
    classify();
}

自己測試，推薦使用著名的異或問題數據集，該數據只有四個點[0,0]、[1,0]、[0,1]、[1,1] ，分類[0,1,1,0]

線性分類器准確無法分類，原因參見：http://www.guokr.com/blog/793310/

使用$\alpha=0.6$、$\lambda=0.8$(動量因子)

如果你的BP網絡代碼沒錯的話，並行9000次,  單次迭代所有的樣本SSE（誤差平方和）大概是這樣。

 很好地解決了，這個線性不可分問題。