~ read.
[.AI]C++下DNN的實現

[.AI]C++下DNN的實現

最近為家裡領導開發一套投資分析系統,需要用到DNN算法(深度學習神經網絡),於是乎,天賦樹上又有了新的分支...

如何用C++實現一個深度神經網絡(DNN, Deep Neural Network)處理多分類邏輯回歸問題(MLR, Multiclass Logistic Regression)且不使用第三方庫?

今天就來嘗試一下: 使用sigmoid+激活函數+Softmax邏輯回歸+交叉熵損失函數+反向傳播實現一個DNN。整個過程只使用C++標準庫而不使用任何第三方庫。

建議食用本文前請掌握神經網絡的基本知識,或者先忽略在「約定和初始化」部分所看到的陌生函數和名詞,在後文中會有相應解釋。

問題

給定𝑘維空間的𝐾個點集(在幾何距離上相近)。每次查詢一個新點,問這個點最有可能屬於哪一個點集。

𝑘在一定範圍內增大不會影響我們的解法,因此這裡我們先從簡處理𝑘=2的情況,即2維平面的點集。

方案

如何表示一個劃分方案?

用一個𝐾×1向量表示,第𝑖個數字表示這個點屬於點集𝑖的概率。

對於初始點集作為給出的標準答案來說一個屬於𝑖點集的點屬於點集𝑖的概率當然是1。

損失函數

如何評價一個劃分方案的好壞?

在監督學習(提供標準答案)的情況下,常用的有兩類評價函數,一種是平方差損失函數(𝑀𝑆𝐸),公式如下:

這個1/2是為了求導方便。𝑀𝑆𝐸是一種很容易想到而且有效性顯而易見的損失函數。不過另一種相對更好的評價函數是交叉熵損失函數(𝐶𝐸),公式如下:

注意這裡:越小的劃分方案我們認為越優,直觀理解就是對於𝑦𝑖比較大的項,追求$\ln(\frac{1}{\hat y_i})$較小即較小即較大,這與𝑦𝑖是一致的。

代碼實現:

inline db MSELoss(){db S = 0; for(int j = 0; j < K; S += (Out[j]-Ans[j])*(Out[j]-Ans[j]), j++); return S/2/K;}  
inline db CELoss(){db S = 0; for(int j = 0; j < K; S += CE(Ans[j],Out[j]), j++); return -S;}  

約定和初始化

為了處理好細節問題,現在這裡對神經網絡的結構做一下說明和實現。

一開始我試圖用一個二維數組來表示神經網絡,每個神經元對應一個坐標(𝑖,𝑗)。這種方法因為內存連續性,常數比較優秀。但是不能很好的處理每層節點數不同的情況,無法自定義空間複雜度已經最大層數和最大節點數。

還有一種常見的方法使用矩陣表示所有的轉移。這樣的實現簡潔不易出錯,代碼更易理解。但是靈活性較差,往往需要使用第三方庫(例如OpenCV)。而手動實現矩陣則憑空增加了代碼複雜度。

於是我選擇用vectorstruct組合來完成這一任務。這種做法是最直觀的(相對於矩陣方法更貼近於神經網絡的物理意義而不是數學意義)。

定義一個神經網絡有𝑁+1層,分別記為第0,1,2,…,𝑁層。

其中第0層為輸入層,只有𝑘=2個節點,輸入形式為𝑥/RANGE,𝑦/RANGE,,其中𝑥為坐標範圍,這樣保證輸入的每個數都與1在一個數量級上,防止sigmoid函數將輸入」抹平」。

所有0到𝑁−1層的激活函數都是sigmoid,但第𝑁層的激活函數為Softmax,作為輸出層,有恰好𝐾個節點。

第𝑖層有𝐿[𝑖].𝑛+1個節點,其中0節點代表偏置,因此0節點的輸出總是1,偏置由邊權調整。以後用(𝑖,𝑗)表示第𝑖層第𝑗個節點。

每個節點有ℎ表示未經激活函數處理的輸入,𝑥表示經過激活函數處理的輸出,𝑑表示反向傳播的誤差。以上內容按層記錄,另外每層𝐿[𝑖].𝑤[𝑗][𝑘]表示這樣一條邊: \[(i-1,k) \xrightarrow{L[i].w[j][k]} (i,j)\]

𝜂學習率,代碼中為eta。In、Out以及Ans是對外的接口,為了增強通用性,這些下標從0開始,注意這意味著接口和神經網絡有1的下標偏差。

因此基本結構的代碼實現如下:

#define db double
#define range(x) (x).begin(),(x).end()
#define sigmoid(x) (1./(1.+exp(-(x))))
#define CE(a,b) ((a)*log(max(b,1e-12))+(1-(a))*log(max(1-(b),1e-12)))
typedef vector<db> Arr; mt19937 _rd(1061109589); auto real_rd = std::bind(std::normal_distribution<db>(0.,.01),mt19937(998244353));  
struct NeuralNetwork{  
    struct Layer{
        Arr h,x,d; vector<Arr> w; int n;
        inline Layer(int _n, int c){
            n = _n; x = h = d = Arr(n+1); x[0] = 1; w.clear();
            for(int i = 0, j; i <= n; w[i][0] = 0.5, i++)
                for(w.push_back(Arr(c+1)), j = 1; j <= c; w[i][j] = real_rd(), j++);
        }
    };  vector<Layer> L; Arr Out, Ans, In; int N, K; db eta;
    inline NeuralNetwork(vector<int> Num, db _eta){L.clear(); eta = _eta; N = -1; K = 0; for(int&j:Num) AddLayer(j);}
}

這裡的real_rd是一個正態分布𝜇=0,𝜎=0.01,即𝑁(0,0.01)所有的邊權被隨機賦值為一個接近於0而很小的數,而偏置常數項被設定為中值0.5。

𝜂設置為0.01∼0.1。

注意𝑤的初始化要求比較高: 顯然不能全部初始化為0,這樣整個網絡永遠都是0也不能全部初始化為一個數,這樣同一層所有節點是完全對稱的地位,結果也會完全一樣,節點數目就形同虛設了;也不能全部隨機初始化,因為∑𝑤𝑥會導致sigmoid函數值全都趨向於1。

神經元和前向傳播

如何實現一個神經元?

之前提到過$(i-1,k) \xrightarrow{L[i].w[j][k]} (i,j)$至於為什麼使用𝐿[𝑖].𝑤[𝑗][𝑘]而不是𝐿[𝑖−1].𝑤[𝑘][𝑗],是為了後面代碼方便。接下來簡記為𝑤𝑖,𝑗,𝑘,其他同理。

具體地,神經元(𝑖,𝑗)受到上一層神經元的影響是每個神經元的輸出(𝑥𝑖−1,𝑘)乘以邊權:

可以看到,因為採用𝑤𝑖,𝑗,𝑘而不是𝑤𝑖−1,𝑘,𝑗,這裡可以記為向量點積(inner_product)。𝑤𝑖,𝑗,0表示的是偏置常數項。實現的時候,偏置常數項的輸出𝑥𝑖,0永遠是1。

然後每個神經元將其受到的影響經過激活函數處理後輸出:

其中$\text{sigmoid}(x) = \frac{1}{1+e^{-x}}$也可以簡記為𝜎(𝑥)。注意到𝜎函數值域是(0,1),只能實現二分類(判斷靠近0還是靠近1),一個它的拓展是Softmax函數,是一個序列到等長序列的映射:

事實上就是經過$e^x$放大以後求出每個數字佔整體的比例,來得到一個概率分布。

代碼實現:

inline void Run(){  
    for(int i = (copy(range(In),++L[0].x.begin()),1), j; i <= N; i++)
        for(j = 1; j <= L[i].n; L[i].x[j]=sigmoid(L[i].h[j]=inner_product(range(L[i-1].x),L[i].w[j].begin(),0.)), j++);
    Out.resize(K); for(int j = 0; j < K; Out[j] = exp(L[N].h[j+1]), j++); db S = accumulate(range(Out),0.); for(db&j:Out) j /= S;
}
inline int Judge(){return (int)(max_element(range(Out))-Out.begin());}  

Judge()就是輸出概率分布中最大值作為神經網絡最後的分類,只在檢測正確率的時候用到。這裡的前向傳播和最後的Softmax函數合併了,也可以單獨提取出Softmax函數:

inline void Softmax(){Out.resize(K); for(int j = 0; j < K; Out[j] = exp(L[N].h[j+1]), j++); db S = accumulate(range(Out),0.); for(db&j:Out) j /= S;}  

梯度下降反向傳播

現在每個𝑤都是一開始設定好的,這個神經網絡根本沒有學習能力。如何讓神經網絡具有學習能力?

關鍵就是通過大量數據調整所有𝑤𝑖,𝑗,𝑘的值,使得最後的𝐿𝑜𝑠𝑠最小。

一個很常見的思路就是求導: 每次貪心地將每個$w_{i,j,k}$嘗試調大或者調小。至於具體調大還是調小,當然就是向著改變以後𝐿𝑜𝑠𝑠會降低的方向調整。

𝐿𝑜𝑠𝑠降低的方向就是𝐿𝑜𝑠𝑠關於$w_{i,j,k}$的偏導數正向的反方向。即每次調整:

這種方法就是梯度下降算法。𝜂是學習率(就是下降的步長)。

先來嘗試推導最後一層節點的偏導數。首先回憶之前的交叉熵損失函數以及Softmax函數:

就可以得到:

最後一個等式是因為𝐴𝑛𝑠是一個概率分布:

於是就有:

利用鏈導法則推廣到一般情況:

如果記:

那麼就有:

而:

只要能求出後半部分偏導數就做完了:

然後這裡就需要用到𝜎(𝑥)的美妙性質了:'

因此:

DONE!!

代碼:

inline void Adjust(){  
    for(int j = 1; j <= K; L[N].d[j] = -eta*(Out[j-1]-Ans[j-1]), j++);
    for(int i = N-1, j, k; i; i--) for(fill(range(L[i].d),0.), j = 1; j <= L[i+1].n; j++)
        for(k = 0; k <= L[i].n; L[i].d[k] += L[i+1].d[j]*L[i+1].w[j][k]*L[i].x[k]*(1-L[i].x[k]), k++);
    for(int i = N, j, k; i; i--) for(j = 1; j <= L[i].n; j++)
        for(k = 0; k <= L[i-1].n; L[i].w[j][k] += L[i].d[j]*L[i-1].x[k], k++);
}

注意一下前向傳播和反向傳播複雜度都是$O(Nn^2)$的。而反向傳播由於鏈式法則求導運算複雜得多,並且最重要的是反向傳播不滿足內存連續性。因此反向傳播這個部分可以說是整個神經網絡的複雜度瓶頸: 一個微小的改動可能造成巨大的常數差異。例如:

L[i].d[k] += L[i+1].d[j]*L[i+1].w[j][k]*L[i].x[k]*(1-L[i].x[k])  

如果按照之前的公式直接寫為:

L[i].d[k] += L[i+1].d[j]*L[i+1].w[j][k]*sigmoid(L[i].h[k])*(1-sigmoid(L[i].h[k]))  

由於多調用了兩次sigmoid函數,而exp函數又十分緩慢,僅僅這一點差別可以讓運行時間相差數倍至數十倍(親測血的教訓)。

神經網絡

完整DNN代碼:

// NN.h
#include <bits/stdc++.h>
using namespace std;  
#define db double
#define range(x) (x).begin(),(x).end()
#define sigmoid(x) (1./(1.+exp(-(x))))
#define CE(a,b) ((a)*log(max(b,1e-12))+(1-(a))*log(max(1-(b),1e-12)))
typedef vector<db> Arr; mt19937 _rd(1061109589); auto real_rd = std::bind(std::normal_distribution<db>(0.,.01),mt19937(998244353));  
struct NeuralNetwork{  
    struct Layer{
        Arr h,x,d; vector<Arr> w; int n;
        inline Layer(int _n, int c){
            n = _n; x = h = d = Arr(n+1); x[0] = 1; w.clear();
            for(int i = 0, j; i <= n; w[i][0] = 0.5, i++)
                for(w.push_back(Arr(c+1)), j = 1; j <= c; w[i][j] = real_rd(), j++);
        }
    };  vector<Layer> L; Arr Out, Ans, In; int N, K; db eta; inline void AddLayer(int n){++N; L.emplace_back(n,K); K = n;}
    inline NeuralNetwork(vector<int> Num, db _eta){L.clear(); eta = _eta; N = -1; K = 0; for(int&j:Num) AddLayer(j);}
    inline void Run(){
        for(int i = (copy(range(In),++L[0].x.begin()),1), j; i <= N; i++)
            for(j = 1; j <= L[i].n; L[i].x[j]=sigmoid(L[i].h[j]=inner_product(range(L[i-1].x),L[i].w[j].begin(),0.)), j++);
        Out.resize(K); for(int j = 0; j < K; Out[j] = exp(L[N].h[j+1]), j++); db S = accumulate(range(Out),0.); for(db&j:Out) j /= S;
    }
    inline int Judge(){return (int)(max_element(range(Out))-Out.begin());}
    inline db CELoss(){db S = 0; for(int j = 0; j < K; S += CE(Ans[j],Out[j]), j++); return -S;}
    inline void Adjust(){
        for(int j = 1; j <= K; L[N].d[j] = -eta*(Out[j-1]-Ans[j-1]), j++);
        for(int i = N-1, j, k; i; i--) for(fill(range(L[i].d),0.), j = 1; j <= L[i+1].n; j++)
            for(k = 0; k <= L[i].n; L[i].d[k] += L[i+1].d[j]*L[i+1].w[j][k]*L[i].x[k]*(1-L[i].x[k]), k++);
        for(int i = N, j, k; i; i--) for(j = 1; j <= L[i].n; j++)
            for(k = 0; k <= L[i-1].n; L[i].w[j][k] += L[i].d[j]*L[i-1].x[k], k++);
    }
};
#undef db
#undef range

最後,關於訓練

核心問題是如何造數據。這裡自己YY了一下,採用這種造數據的方式: 隨機一個點集數目𝐾,對於每個點集,隨機一個矩形區域,在裡面生成大量點。然後反復選擇其中一部分作為初始給出,另一部分作為檢測數據。

下面是一個實例:

#include <bits/stdc++.h>
#include "NN.h" 
using namespace std;  
#define TESTS 5
#define RANGE 100000
#define SET_RANGE 100
#define NUM_PER_SET_MIN 4000
#define NUM_PER_SET_MAX 5000
#define db double
#define rint register int
inline int rf(){int r;int s=0,c;for(;!isdigit(c=getchar());s=c);for(r=c^48;isdigit(c=getchar());(r*=10)+=c^48);return s^45?r:-r;}  
struct P{int x,y,i;}; vector<P> Set[32], Homework, Exam; mt19937 rd(time(0)); int K = rd()%30+1, T = TESTS, _; db Loss, Score;  
int main(){  
    for(rint i = K, j, LX, LY; i; i--)
        for(LX = rd()%RANGE+1, LY = rd()%RANGE+1, j = rd()%(NUM_PER_SET_MAX-NUM_PER_SET_MIN)+NUM_PER_SET_MIN; j--; )
            Set[i].push_back(P{rd()%SET_RANGE+LX,rd()%SET_RANGE+LY,i-1});
    for(rint i, j; T--; printf("%.6lf %.6lf\n",Loss/Exam.size(),Score/Exam.size())){
        for(vector<P>().swap(Homework), vector<P>().swap(Exam), i = 1; i <= K; i++)
            for(shuffle(Set[i].begin(),Set[i].end(),rd), j = 0; j < Set[i].size(); j++)
                (j<NUM_PER_SET_MIN-10?Homework:Exam).push_back(Set[i][j]);
        NeuralNetwork NN({2,50,K},.1);
//      NeuralNetwork NN({2,100,K},.1);
//      NeuralNetwork NN({2,200,K},.1);
//      NeuralNetwork NN({2,30,30,K},.1);
//      NeuralNetwork NN({2,50,50,K},.1);
//      NeuralNetwork NN({2,100,100,K},.1);
        Loss = Score = 0;
        shuffle(Homework.begin(),Homework.end(),rd); _ = 0;
        for(auto p:Homework){
            NN.In = {(db)p.x/RANGE,(db)p.y/RANGE};
            NN.Ans = Arr(K); NN.Ans[p.i] = 1;
            NN.Run(); NN.Adjust();
            /*
            printf("Training #%d:\n",++_);
            for(j = 0; j < K; printf("%4.4lf ",NN.Out[j]), j++); puts("");
            for(j = 0; j < K; printf("%4.4lf ",NN.Ans[j]), j++); puts("");
            printf("Output = %d, Answer = %d\n\n",NN.Judge(),p.i);
            //*/
        }
        shuffle(Exam.begin(),Exam.end(),rd); _ = 0;
        for(auto p:Exam){
            NN.In = {(db)p.x/RANGE,(db)p.y/RANGE};
            NN.Ans = Arr(K); NN.Ans[p.i] = 1;
            NN.Run(); Loss += NN.CELoss(); Score += NN.Judge()==p.i; 
            /*
            printf("Testing #%d:\n",++_);
            for(j = 0; j < K; printf("%4.4lf ",NN.Out[j]), j++); puts("");
            for(j = 0; j < K; printf("%4.4lf ",NN.Ans[j]), j++); puts("");
            printf("Output = %d, Answer = %d\n\n",NN.Judge(),p.i);
            //*/
        }
    }   return 0;
}

當提供的數據在4000∗𝐾左右時,

令人出乎意料的是𝑁𝑁({2,50,𝐾},.1)最小的神經網絡竟然表現最好,基本可以保持正確率在90%以上,且絕大多數情況下都是100%正確。

𝑁𝑁({2,200,𝐾},.1)表現也很好,基本可以保持正確率在90% 以上,且多數情況下是100%正確。

𝑁𝑁({2,30,30,𝐾},.1)基本可以保持正確率在70%以上。

𝑁𝑁({2,100,100,𝐾},.1)基本可以保持正確率在85%以上,經常可以達到95%甚至100%。

可能是由於訓練數據集大小、神經網絡對於這一類問題的適應性、隨機的不確定性等原因,單隱層神經網絡表現得比多隱層神經網絡要好,而且隱層節點數適度減少並不一定會導致正確率下滑。

一个小Trick:

//*
//*/

可以方便地實現多行注釋,只要刪除第一個/就可以轉變為注釋狀態,加上第一個/就轉變為不注釋狀態。

———————————
© 2020. All rights reserved. Setup on Dedicate Server in Europe . Built by Mille Yin.