C++实现深度神经网络解决多分类逻辑回归问题

如何用C++实现一个深度神经网络(DNN, Deep Neural Network)处理多分类逻辑回归问题(MLR, Multiclass Logistic Regression)且不使用第三方库

唔…

今天就来尝试一下: 使用sigmoid\text{sigmoid}激活函数+Softmax\text{Softmax}逻辑回归+交叉熵损失函数+反向传播实现一个DNN。整个过程只使用C++标准库而不使用任何第三方库

建议食用本文前请掌握神经网络的基本知识,或者先忽略在"约定和初始化"部分所看到的陌生函数和名词,在后文中会有相应解释。

P.S.: 由于学习Go语言的缘故,我开始左大括号不换行了QwQ。

问题

给定kk维空间的KK个点集(在几何距离上相近)。每次查询一个新点,问这个点最有可能属于哪一个点集。

kk在一定范围内增大不会影响我们的解法,因此这里我们先从简处理k=2k = 2的情况,即22维平面的点集。


方案

如何表示一个划分方案?

用一个K×1K \times 1向量表示,第ii个数字表示这个点属于点集ii的概率。

对于初始点集作为给出的标准答案来说一个属于ii点集的点属于点集ii的概率当然是11


损失函数

如何评价一个划分方案的好坏?

在监督学习(提供标准答案)的情况下,常用的有两类评价函数,一种是平方差损失函数(MSEMSE),公式如下:

Loss=12Ki=1K(yiy^i)2Loss = \frac{1}{2K}\sum_{i=1}^K (y_i - \hat y_i)^2

这个1/21/2是为了求导方便。MSEMSE是一种很容易想到而且有效性显而易见的损失函数。不过另一种相对更好的评价函数是交叉熵损失函数(CECE),公式如下:

Loss=i=1Kyiln(y^i)Loss = -\sum_{i=1}^Ky_i\ln(\hat y_i)

注意这里i=1Kyi=1\sum_{i=1}^Ky_i = 1LossLoss越小的划分方案我们认为越优,直观理解就是对于yiy_i比较大的项,追求ln(1y^i)\ln(\frac{1}{\hat y_i})较小即1y^i\frac{1}{\hat y_i}较小即y^i\hat y_i较大,这与yiy_i是一致的。


代码实现:

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)(i,j)。这种方法因为内存连续性,常数比较优秀。但是不能很好的处理每层节点数不同的情况,无法自定义空间复杂度已经最大层数和最大节点数。

还有一种常见的方法使用矩阵表示所有的转移。这样的实现简洁不易出错,代码更易理解。但是灵活性较差,往往需要使用第三方库(例如OpenCV)。而手动实现矩阵则凭空增加了代码复杂度。

于是我选择用vectorstruct组合来完成这一任务。这种做法是最直观的(相对于矩阵方法更贴近于神经网络的物理意义而不是数学意义)。


定义一个神经网络有N+1N+1层,分别记为第0,1,2,...,N0,1,2,...,N层。

其中第00层为输入层,只有k=2k = 2个节点,输入形式为x/RANGE,y/RANGEx/\text{RANGE},y/\text{RANGE},其中xx为坐标范围,这样保证输入的每个数都与11在一个数量级上,防止sigmoid\text{sigmoid}函数将输入"抹平"。

所有00N1N-1层的激活函数都是sigmoid\text{sigmoid},但第NN层的激活函数为Softmax\text{Softmax},作为输出层,有恰好KK个节点。

ii层有L[i].n+1L[i].n+1个节点,其中00节点代表偏置,因此00节点的输出总是11,偏置由边权调整。以后用(i,j)(i,j)表示第ii层第jj个节点。

每个节点有hh表示未经激活函数处理的输入,xx表示经过激活函数处理的输出,dd表示反向传播的误差。以上内容按层记录,另外每层L[i].w[j][k]L[i].w[j][k]表示这样一条边: (i1,k)L[i].w[j][k](i,j)(i-1,k) \xrightarrow{L[i].w[j][k]} (i,j)

η\eta是学习率,代码中为etaInOut以及Ans是对外的接口,为了增强通用性,这些下标从00开始,注意这意味着接口和神经网络有11的下标偏差。


因此基本结构的代码实现如下:

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\mu = 0, \sigma = 0.01,即N(0,0.01)N(0,0.01)

所有的边权被随机赋值为一个接近于00而很小的数,而偏置常数项被设定为中值0.50.5

η\eta设置为0.010.10.01 \sim 0.1

我也不知道为啥这样玄学调参

注意ww的初始化要求比较高: 显然不能全部初始化为00,这样整个网络永远都是00;也不能全部初始化为一个数,这样同一层所有节点是完全对称的地位,结果也会完全一样,节点数目就形同虚设了;也不能全部随机初始化,因为wx\sum wx很大会导致sigmoid\text{sigmoid}函数值全都趋向于11


神经元和前向传播

如何实现一个神经元?

之前提到过(i1,k)L[i].w[j][k](i,j)(i-1,k) \xrightarrow{L[i].w[j][k]} (i,j),至于为什么使用L[i].w[j][k]L[i].w[j][k]而不是L[i1].w[k][j]L[i-1].w[k][j],是为了后面代码方便。接下来简记为wi,j,kw_{i,j,k},其他同理。

具体地,神经元(i,j)(i,j)受到上一层神经元的影响是每个神经元的输出(xi1,kx_{i-1,k})乘以边权:

hi,j=k=0ni1wi,j,kxi1,k=Wi,jXi1h_{i,j} = \sum_{k=0}^{n_{i-1}} w_{i,j,k} \cdot x_{i-1,k} = W_{i,j} \cdot X_{i-1}

可以看到,因为采用wi,j,kw_{i,j,k}而不是wi1,k,jw_{i-1,k,j},这里可以记为向量点积(inner_product)。wi,j,0w_{i,j,0}表示的是偏置常数项。实现的时候,偏置常数项的输出xi,0x_{i,0}永远是11

然后每个神经元将其受到的影响经过激活函数处理后输出:

xi,j=sigmoid(hi,j)=sigmoid(Wi,jXi1)x_{i,j} = \text{sigmoid}(h_{i,j}) = \text{sigmoid}(W_{i,j} \cdot X_{i-1})

其中sigmoid(x)=11+ex\text{sigmoid}(x) = \frac{1}{1+e^{-x}},也可以简记为σ(x)\sigma(x)。注意到σ\sigma函数值域是(0,1)(0,1),只能实现二分类(判断靠近00还是靠近11),一个它的拓展是Softmax\text{Softmax}函数,是一个序列到等长序列的映射:

hN,jehN,jiehN,i=Outjh_{N,j} \rightarrow \frac{e^{h_{N,j}}}{\sum_i e^{h_{N,i}}} = Out_j

事实上就是经过exe^x放大以后求出每个数字占整体的比例,来得到一个概率分布。


代码实现:

运行神经网络:

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\text{Softmax}函数合并了,也可以单独提取出Softmax\text{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;}

梯度下降反向传播

系好安全带,开始飙车了!

现在每个ww都是一开始设定好的,这个神经网络根本没有学习能力。如何让神经网络具有学习能力?

关键就是通过大量数据调整所有wi,j,kw_{i,j,k}的值,使得最后的LossLoss最小。

一个很常见的思路就是求导: 每次贪心地将每个wi,j,kw_{i,j,k}尝试调大或者调小。至于具体调大还是调小,当然就是向着改变以后LossLoss会降低的方向调整。LossLoss降低的方向就是LossLoss关于wi,j,kw_{i,j,k}的偏导数正向的反方向。即每次调整:

wi,j,k=wi,j,kηLosswi,j,kw_{i,j,k} = w_{i,j,k} - \eta \frac{\partial Loss}{\partial w_{i,j,k}}

这种方法就是梯度下降算法。η\eta是学习率(就是下降的步长)。

先来尝试推导最后一层节点的偏导数。首先回忆之前的交叉熵损失函数以及Softmax\text{Softmax}函数:

Loss=i=1KAnsiln(Out^i)Outj=ehN,jiehN,i=ehN,jSumLoss = - \sum_{i=1}^KAns_i\ln(\hat{Out}_i) \\ Out_j = \frac{e^{h_{N,j}}}{\sum_i e^{h_{N,i}}} = \frac{e^{h_{N,j}}}{Sum} \\

就可以得到:

OutihN,j=ehN,iehN,jSum2=OutiOutjOutjhN,j=SumijehN,jSum2=(1Outj)Outj\frac{\partial Out_i}{\partial h_{N,j}} = -\frac{e^{h_{N,i}}e^{h_{N,j}}}{Sum^2} = -Out_iOut_j \\ \frac{\partial Out_j}{\partial h_{N,j}} = \frac{Sum_{i \ne j} e^{h_{N,j}}}{Sum^2} = (1-Out_j)Out_j \\

LosshN,j=i=1KAnsi1OutiOutihN,j=Ansj(1Outj)+ijAnsiOutj=Ansj+i=1KAnsiOutj=OutjAnsj\frac{\partial Loss}{\partial h_{N,j}} = -\sum_{i=1}^K Ans_i\frac{1}{Out_i}\frac{\partial Out_i}{\partial h_{N,j}} = -Ans_j(1-Out_j) + \sum_{i \ne j} Ans_iOut_j \\ = -Ans_j + \sum_{i=1}^K Ans_iOut_j = Out_j - Ans_j

最后一个等式是因为AnsAns是一个概率分布,i=1KAnsi=1\sum_{i=1}^K Ans_i = 1。于是就有:

LosswN,j,k=LosshN,jhN,jwN,j,k=(OutjAnsj)xN1,k\frac{\partial Loss}{\partial w_{N,j,k}} = \frac{\partial Loss}{\partial h_{N,j}} \frac{\partial h_{N,j}}{\partial w_{N,j,k}} = (Out_j - Ans_j) x_{N-1,k}

利用链导法则推广到一般情况: 如果记di,j=ηLosshi,jd_{i,j} = - \eta \frac{\partial Loss}{\partial h_{i,j}},那么就有Δwi,j,k=di,jxi1,k\Delta w_{i,j,k} = d_{i,j}x_{i-1,k},而:

di,k=ηLosshi,k=ηLosshi+1,jhi+1,jhi,k=di+1,jhi+1,jhi,kd_{i,k} = - \eta \frac{\partial Loss}{\partial h_{i,k}} = - \eta \frac{\partial Loss}{\partial h_{i+1,j}}\frac{\partial h_{i+1,j}}{\partial h_{i,k}} = d_{i+1,j} \frac{\partial h_{i+1,j}}{\partial h_{i,k}}

只要能求出后半部分偏导数就做完了:

hi+1,jhi,k=txi,twi+1,j,thi,k=σ(hi,k)wi+1,j,khi,k=wi+1,j,kσ(hi,k)\frac{\partial h_{i+1,j}}{\partial h_{i,k}} = \frac{\partial \sum_{t} x_{i,t}w_{i+1,j,t}}{\partial h_{i,k}} = \frac{\partial \sigma(h_{i,k})w_{i+1,j,k}}{\partial h_{i,k}} = w_{i+1,j,k}\sigma'(h_{i,k})

然后这里就需要用到σ(x)\sigma(x)的美妙性质了:

σ(x)=(11+ex)=ex(1+ex)2=σ(x)(1σ(x))\sigma'(x) = (\frac{1}{1+e^{-x}})' = \frac{e^{-x}}{(1+e^{-x})^2} = \sigma(x)(1 - \sigma(x))

因此:

di,k=di+1,jwi+1,j,kσ(hi,k)(1σ(hi,k))=di+1,jwi+1,j,kxi,k(1xi,k)d_{i,k} = d_{i+1,j} w_{i+1,j,k} \sigma(h_{i,k})(1-\sigma(h_{i,k})) = d_{i+1,j}w_{i+1,j,k}x_{i,k}(1-x_{i,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)O(Nn^2)的。而反向传播由于链式法则求导运算复杂得多,并且最重要的是反向传播不满足内存连续性。因此反向传播这个部分可以说是整个神经网络的复杂度瓶颈: 一个微小的改动可能造成巨大的常数差异。例如:

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
// 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

瞧!一个深度神经网络的核心代码居然只有30行左右!


训练

如何训练神经网络?

核心问题是如何造数据。这里自己YY了一下,采用这种造数据的方式: 随机一个点集数目KK,对于每个点集,随机一个矩形区域,在里面生成大量点。然后反复选择其中一部分作为初始给出,另一部分作为检测数据。

下面是一个实例:

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

这里先离题说一个小Trick:

1
2
//*
//*/

可以方便地实现多行注释,只要删除第一个/就可以转变为注释状态,加上第一个/就转变为不注释状态。


回到正题,当提供的数据在4000K4000 * K左右时,

令人出乎意料的是NN({2,50,K},.1)NN(\lbrace2,50,K\rbrace,.1)作为最小的神经网络竟然表现最好,基本可以保持正确率在90%90 \%以上,且绝大多数情况下都是100%100 \%正确。

NN({2,200,K},.1)NN(\lbrace2,200,K\rbrace,.1)表现也很好,基本可以保持正确率在90%90 \%以上,且多数情况下是100%100 \%正确。

NN({2,30,30,K},.1)NN(\lbrace2,30,30,K\rbrace,.1)基本可以保持正确率在70%70 \%以上。

NN({2,100,100,K},.1)NN(\lbrace2,100,100,K\rbrace,.1)基本可以保持正确率在85%85 \%以上,经常可以达到95%95 \%甚至100%100 \%

可能是由于训练数据集大小、神经网络对于这一类问题的适应性、随机的不确定性等原因,单隐层神经网络表现得比多隐层神经网络要好,而且隐层节点数适度减少并不一定会导致正确率下滑。


扫描二维码即可在手机上查看这篇文章,或者转发二维码来分享这篇文章:


文章作者: Magolor
文章链接: https://magolor.cn/2019/08/02/2019-08-02-blog-01/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Magolor
扫描二维码在手机上查看或转发二维码以分享Magolor的博客