如何用C++实现一个深度神经网络(DNN, Deep Neural Network)处理多分类逻辑回归问题(MLR, Multiclass Logistic Regression)且不使用第三方库?
唔…
今天就来尝试一下: 使用sigmoid激活函数+Softmax逻辑回归+交叉熵损失函数+反向传播实现一个DNN。整个过程只使用C++标准库而不使用任何第三方库。
建议食用本文前请掌握神经网络的基本知识,或者先忽略在"约定和初始化"部分所看到的陌生函数和名词,在后文中会有相应解释。
P.S.: 由于学习Go语言的缘故,我开始左大括号不换行了QwQ。
问题
给定k维空间的K个点集(在几何距离上相近)。每次查询一个新点,问这个点最有可能属于哪一个点集。
k在一定范围内增大不会影响我们的解法,因此这里我们先从简处理k=2的情况,即2维平面的点集。
方案
如何表示一个划分方案?
用一个K×1向量表示,第i个数字表示这个点属于点集i的概率。
对于初始点集作为给出的标准答案来说一个属于i点集的点属于点集i的概率当然是1。
损失函数
如何评价一个划分方案的好坏?
在监督学习(提供标准答案)的情况下,常用的有两类评价函数,一种是平方差损失函数(MSE),公式如下:
Loss=2K1i=1∑K(yi−y^i)2
这个1/2是为了求导方便。MSE是一种很容易想到而且有效性显而易见的损失函数。不过另一种相对更好的评价函数是交叉熵损失函数(CE),公式如下:
Loss=−i=1∑Kyiln(y^i)
注意这里∑i=1Kyi=1。Loss越小的划分方案我们认为越优,直观理解就是对于yi比较大的项,追求ln(y^i1)较小即y^i1较小即y^i较大,这与yi是一致的。
代码实现:
1 2
| 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;}
|
约定和初始化
为了处理好细节问题,现在这里对神经网络的结构做一下说明和实现。
一开始我试图用一个二维数组来表示神经网络,每个神经元对应一个坐标(i,j)。这种方法因为内存连续性,常数比较优秀。但是不能很好的处理每层节点数不同的情况,无法自定义空间复杂度已经最大层数和最大节点数。
还有一种常见的方法使用矩阵表示所有的转移。这样的实现简洁不易出错,代码更易理解。但是灵活性较差,往往需要使用第三方库(例如OpenCV
)。而手动实现矩阵则凭空增加了代码复杂度。
于是我选择用vector
和struct
组合来完成这一任务。这种做法是最直观的(相对于矩阵方法更贴近于神经网络的物理意义而不是数学意义)。
定义一个神经网络有N+1层,分别记为第0,1,2,...,N层。
其中第0层为输入层,只有k=2个节点,输入形式为x/RANGE,y/RANGE,其中x为坐标范围,这样保证输入的每个数都与1在一个数量级上,防止sigmoid函数将输入"抹平"。
所有0到N−1层的激活函数都是sigmoid,但第N层的激活函数为Softmax,作为输出层,有恰好K个节点。
第i层有L[i].n+1个节点,其中0节点代表偏置,因此0节点的输出总是1,偏置由边权调整。以后用(i,j)表示第i层第j个节点。
每个节点有h表示未经激活函数处理的输入,x表示经过激活函数处理的输出,d表示反向传播的误差。以上内容按层记录,另外每层L[i].w[j][k]表示这样一条边: (i−1,k)L[i].w[j][k](i,j)。
η是学习率,代码中为eta
。In
、Out
以及Ans
是对外的接口,为了增强通用性,这些下标从0开始,注意这意味着接口和神经网络有1的下标偏差。
因此基本结构的代码实现如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
| #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,即N(0,0.01)。
所有的边权被随机赋值为一个接近于0而很小的数,而偏置常数项被设定为中值0.5。
η设置为0.01∼0.1。
我也不知道为啥这样玄学调参
注意w的初始化要求比较高: 显然不能全部初始化为0,这样整个网络永远都是0;也不能全部初始化为一个数,这样同一层所有节点是完全对称的地位,结果也会完全一样,节点数目就形同虚设了;也不能全部随机初始化,因为∑wx很大会导致sigmoid函数值全都趋向于1。
神经元和前向传播
如何实现一个神经元?
之前提到过(i−1,k)L[i].w[j][k](i,j),至于为什么使用L[i].w[j][k]而不是L[i−1].w[k][j],是为了后面代码方便。接下来简记为wi,j,k,其他同理。
具体地,神经元(i,j)受到上一层神经元的影响是每个神经元的输出(xi−1,k)乘以边权:
hi,j=k=0∑ni−1wi,j,k⋅xi−1,k=Wi,j⋅Xi−1
可以看到,因为采用wi,j,k而不是wi−1,k,j,这里可以记为向量点积(inner_product
)。wi,j,0表示的是偏置常数项。实现的时候,偏置常数项的输出xi,0永远是1。
然后每个神经元将其受到的影响经过激活函数处理后输出:
xi,j=sigmoid(hi,j)=sigmoid(Wi,j⋅Xi−1)
其中sigmoid(x)=1+e−x1,也可以简记为σ(x)。注意到σ函数值域是(0,1),只能实现二分类(判断靠近0还是靠近1),一个它的拓展是Softmax函数,是一个序列到等长序列的映射:
hN,j→∑iehN,iehN,j=Outj
事实上就是经过ex放大以后求出每个数字占整体的比例,来得到一个概率分布。
代码实现:
运行神经网络:
1 2 3 4 5 6
| 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函数:
1
| 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都是一开始设定好的,这个神经网络根本没有学习能力。如何让神经网络具有学习能力?
关键就是通过大量数据调整所有wi,j,k的值,使得最后的Loss最小。
一个很常见的思路就是求导: 每次贪心地将每个wi,j,k尝试调大或者调小。至于具体调大还是调小,当然就是向着改变以后Loss会降低的方向调整。Loss降低的方向就是Loss关于wi,j,k的偏导数正向的反方向。即每次调整:
wi,j,k=wi,j,k−η∂wi,j,k∂Loss
这种方法就是梯度下降算法。η是学习率(就是下降的步长)。
先来尝试推导最后一层节点的偏导数。首先回忆之前的交叉熵损失函数以及Softmax函数:
Loss=−i=1∑KAnsiln(Out^i)Outj=∑iehN,iehN,j=SumehN,j
就可以得到:
∂hN,j∂Outi=−Sum2ehN,iehN,j=−OutiOutj∂hN,j∂Outj=Sum2Sumi=jehN,j=(1−Outj)Outj
∂hN,j∂Loss=−i=1∑KAnsiOuti1∂hN,j∂Outi=−Ansj(1−Outj)+i=j∑AnsiOutj=−Ansj+i=1∑KAnsiOutj=Outj−Ansj
最后一个等式是因为Ans是一个概率分布,∑i=1KAnsi=1。于是就有:
∂wN,j,k∂Loss=∂hN,j∂Loss∂wN,j,k∂hN,j=(Outj−Ansj)xN−1,k
利用链导法则推广到一般情况: 如果记di,j=−η∂hi,j∂Loss,那么就有Δwi,j,k=di,jxi−1,k,而:
di,k=−η∂hi,k∂Loss=−η∂hi+1,j∂Loss∂hi,k∂hi+1,j=di+1,j∂hi,k∂hi+1,j
只要能求出后半部分偏导数就做完了:
∂hi,k∂hi+1,j=∂hi,k∂∑txi,twi+1,j,t=∂hi,k∂σ(hi,k)wi+1,j,k=wi+1,j,kσ′(hi,k)
然后这里就需要用到σ(x)的美妙性质了:
σ′(x)=(1+e−x1)′=(1+e−x)2e−x=σ(x)(1−σ(x))
因此:
di,k=di+1,jwi+1,j,kσ(hi,k)(1−σ(hi,k))=di+1,jwi+1,j,kxi,k(1−xi,k)
大成功!!!
代码实现:
1 2 3 4 5 6 7
| 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(Nn2)的。而反向传播由于链式法则求导运算复杂得多,并且最重要的是反向传播不满足内存连续性。因此反向传播这个部分可以说是整个神经网络的复杂度瓶颈: 一个微小的改动可能造成巨大的常数差异。例如:
1
| L[i].d[k] += L[i+1].d[j]*L[i+1].w[j][k]*L[i].x[k]*(1-L[i].x[k])
|
如果按照之前的公式直接写为:
1
| 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代码:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
| #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
|
瞧!一个深度神经网络的核心代码居然只有30行左右!
训练
如何训练神经网络?
核心问题是如何造数据。这里自己YY了一下,采用这种造数据的方式: 随机一个点集数目K,对于每个点集,随机一个矩形区域,在里面生成大量点。然后反复选择其中一部分作为初始给出,另一部分作为检测数据。
下面是一个实例:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
| #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);
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();
} 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;
} } return 0; }
|
这里先离题说一个小Trick:
可以方便地实现多行注释,只要删除第一个/
就可以转变为注释状态,加上第一个/
就转变为不注释状态。
回到正题,当提供的数据在4000∗K左右时,
令人出乎意料的是NN({2,50,K},.1)作为最小的神经网络竟然表现最好,基本可以保持正确率在90%以上,且绝大多数情况下都是100%正确。
NN({2,200,K},.1)表现也很好,基本可以保持正确率在90%以上,且多数情况下是100%正确。
NN({2,30,30,K},.1)基本可以保持正确率在70%以上。
NN({2,100,100,K},.1)基本可以保持正确率在85%以上,经常可以达到95%甚至100%。
可能是由于训练数据集大小、神经网络对于这一类问题的适应性、随机的不确定性等原因,单隐层神经网络表现得比多隐层神经网络要好,而且隐层节点数适度减少并不一定会导致正确率下滑。
扫描二维码即可在手机上查看这篇文章,或者转发二维码来分享这篇文章:
