如何用C++实现一个深度神经网络(DNN, Deep Neural Network)处理多分类逻辑回归问题(MLR, Multiclass Logistic Regression)且不使用第三方库？
唔……
今天就来尝试一下: 使用$\text{sigmoid}$激活函数+$\text{Softmax}$逻辑回归+交叉熵损失函数+反向传播实现一个DNN。整个过程只使用C++标准库而不使用任何第三方库。
建议食用本文前请掌握神经网络的基本知识，或者先忽略在”约定和初始化”部分所看到的陌生函数和名词，在后文中会有相应解释。
P.S.: 由于学习Go语言的缘故，我开始左大括号不换行了QwQ。
问题
给定$k$维空间的$K$个点集(在几何距离上相近)。每次查询一个新点，问这个点最有可能属于哪一个点集。
$k$在一定范围内增大不会影响我们的解法，因此这里我们先从简处理$k = 2$的情况，即$2$维平面的点集。
方案
如何表示一个划分方案？
用一个$K \times 1$向量表示，第$i$个数字表示这个点属于点集$i$的概率。
对于初始点集作为给出的标准答案来说一个属于$i$点集的点属于点集$i$的概率当然是$1$。
损失函数
如何评价一个划分方案的好坏？
在监督学习(提供标准答案)的情况下，常用的有两类评价函数，一种是平方差损失函数($MSE$)，公式如下:
这个$1/2$是为了求导方便。$MSE$是一种很容易想到而且有效性显而易见的损失函数。不过另一种相对更好的评价函数是交叉熵损失函数($CE$)，公式如下:
注意这里$\sum_{i=1}^Ky_i = 1$。$Loss$越小的划分方案我们认为越优，直观理解就是对于$y_i$比较大的项，追求$\ln(\frac{1}{\hat y_i})$较小即$\frac{1}{\hat y_i}$较小即$\hat y_i$较大，这与$y_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)$。这种方法因为内存连续性，常数比较优秀。但是不能很好的处理每层节点数不同的情况，无法自定义空间复杂度已经最大层数和最大节点数。
还有一种常见的方法使用矩阵表示所有的转移。这样的实现简洁不易出错，代码更易理解。但是灵活性较差，往往需要使用第三方库(例如OpenCV)。而手动实现矩阵则凭空增加了代码复杂度。
于是我选择用vector和struct组合来完成这一任务。这种做法是最直观的(相对于矩阵方法更贴近于神经网络的物理意义而不是数学意义)。
定义一个神经网络有$N+1$层，分别记为第$0,1,2,…,N$层。
其中第$0$层为输入层，只有$k = 2$个节点，输入形式为$x/\text{RANGE},y/\text{RANGE}$，其中$x$为坐标范围，这样保证输入的每个数都与$1$在一个数量级上，防止$\text{sigmoid}$函数将输入”抹平”。
所有$0$到$N-1$层的激活函数都是$\text{sigmoid}$，但第$N$层的激活函数为$\text{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) \xrightarrow{L[i].w[j][k]} (i,j)$。
$\eta$是学习率，代码中为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是一个正态分布$\mu = 0, \sigma = 0.01$，即$N(0,0.01)$。
所有的边权被随机赋值为一个接近于$0$而很小的数，而偏置常数项被设定为中值$0.5$。
$\eta$设置为$0.01 \sim 0.1$。
我也不知道为啥这样玄学调参
注意$w$的初始化要求比较高: 显然不能全部初始化为$0$，这样整个网络永远都是$0$；也不能全部初始化为一个数，这样同一层所有节点是完全对称的地位，结果也会完全一样，节点数目就形同虚设了；也不能全部随机初始化，因为$\sum wx$很大会导致$\text{sigmoid}$函数值全都趋向于$1$。
神经元和前向传播
如何实现一个神经元？
之前提到过$(i-1,k) \xrightarrow{L[i].w[j][k]} (i,j)$，至于为什么使用$L[i].w[j][k]$而不是$L[i-1].w[k][j]$，是为了后面代码方便。接下来简记为$w_{i,j,k}$，其他同理。
具体地，神经元$(i,j)$受到上一层神经元的影响是每个神经元的输出($x_{i-1,k}$)乘以边权:
可以看到，因为采用$w{i,j,k}$而不是$w{i-1,k,j}$，这里可以记为向量点积(inner_product)。$w{i,j,0}$表示的是偏置常数项。实现的时候，偏置常数项的输出$x{i,0}$永远是$1$。
然后每个神经元将其受到的影响经过激活函数处理后输出:
其中$\text{sigmoid}(x) = \frac{1}{1+e^{-x}}$，也可以简记为$\sigma(x)$。注意到$\sigma$函数值域是$(0,1)$，只能实现二分类(判断靠近$0$还是靠近$1$)，一个它的拓展是$\text{Softmax}$函数，是一个序列到等长序列的映射:
事实上就是经过$e^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()就是输出概率分布中最大值作为神经网络最后的分类，只在检测正确率的时候用到。这里的前向传播和最后的$\text{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;}
|
梯度下降反向传播
系好安全带，开始飙车了!
现在每个$w$都是一开始设定好的，这个神经网络根本没有学习能力。如何让神经网络具有学习能力？
关键就是通过大量数据调整所有$w_{i,j,k}$的值，使得最后的$Loss$最小。
一个很常见的思路就是求导: 每次贪心地将每个$w{i,j,k}$尝试调大或者调小。至于具体调大还是调小，当然就是向着改变以后$Loss$会降低的方向调整。$Loss$降低的方向就是$Loss$关于$w{i,j,k}$的偏导数正向的反方向。即每次调整:
这种方法就是梯度下降算法。$\eta$是学习率(就是下降的步长)。
先来尝试推导最后一层节点的偏导数。首先回忆之前的交叉熵损失函数以及$\text{Softmax}$函数:
就可以得到:
最后一个等式是因为$Ans$是一个概率分布，$\sum_{i=1}^K Ans_i = 1$。于是就有:
利用链导法则推广到一般情况: 如果记$d{i,j} = - \eta \frac{\partial Loss}{\partial h{i,j}}$，那么就有$\Delta w{i,j,k} = d{i,j}x_{i-1,k}$，而:
只要能求出后半部分偏导数就做完了:
然后这里就需要用到$\sigma(x)$的美妙性质了:
因此:
大成功！！！
代码实现:
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(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
|
#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(\lbrace2,50,K\rbrace,.1)$作为最小的神经网络竟然表现最好，基本可以保持正确率在$90 \%$以上，且绝大多数情况下都是$100 \%$正确。
$NN(\lbrace2,200,K\rbrace,.1)$表现也很好，基本可以保持正确率在$90 \%$以上，且多数情况下是$100 \%$正确。
$NN(\lbrace2,30,30,K\rbrace,.1)$基本可以保持正确率在$70 \%$以上。
$NN(\lbrace2,100,100,K\rbrace,.1)$基本可以保持正确率在$85 \%$以上，经常可以达到$95 \%$甚至$100 \%$。
可能是由于训练数据集大小、神经网络对于这一类问题的适应性、随机的不确定性等原因，单隐层神经网络表现得比多隐层神经网络要好，而且隐层节点数适度减少并不一定会导致正确率下滑。
扫描二维码即可在手机上查看这篇文章，或者转发二维码来分享这篇文章:
