前言

本文主要介绍最基础的神经网络，包括其结构，学习方法， $\text{C++}$ 的实现代码。$\text{Python}$ 的代码可以搜索互联网得到。

前排提示：本人涉及一丁点数学知识。

神经网络的结构

神经网络包括多个层次(Layer)。一般来说神经网络包括三个部分：输入层，隐藏层，输出层 。输入层和输出层都只有一个，而隐藏层可以有多个。顺序为 输入层，隐藏层，输出层。

每一个层次会有若干的神经元。我们称第 $L_l$ 层为输出层，也就是最后一层。对于第 $L$ 层（从 $0$ 开始计数），我们用 $n_L$ 表示其神经元的数量。

对于第 $L$ 层的第 $i$ 个神经元（从 $0$ 开始计数），我们用 $a_i^L$ 表示其输出到下一层神经元或者作为结果的值。对于生物体的神经元，$a_i^L$ 的取值只有 $0,1$ ，但是对于本文而言，计算机的神经元有 $a_i^L\in [0,1]$ 。同时，每个神经元还有一个偏移值 $b_i^L$ 。

除输入层外，每一个神经元都与上一层的所有神经元都有一个连接，对于第 $L(1<L)$ 的第 $i$ 个神经元，它与 $L-1$ 层的所有神经元都有连接，其中它与第 $j$ 个神经元的连接有权重 $w_{i,j}^L$ 。

下图展示了一个简单的神经网络的构成：

图中有 $n_0=3,n_1=n_2=2$ ，且第 $0$ 层为输入层，第 $2$ 层为输出层。

$\text{PS}$: 此图没有展示神经元的偏移值。

输出值的计算

对于输入层而言，其输出值 $a_i^0(i\in [0,n_0)\cap \mathbb{Z})$ ，其输出值与神经网络的输入相同。

对于隐藏层和输出层，其计算方式如下：
$$\begin{array}{rl}{a}_i^L& =f\left(z_i^L\right)\\ z_i^L& =b_i^L+\sum _{j=0}^{n_{L-1}}{w}_{i,j}^L{a}_j^{L-1}\end{array}$$
其中，$f(x)$ 函数一般为 $\mathrm{Sigmoid}$ 函数，常见于生物种群数量的自然增长。其定义和导数如下：
$$\begin{array}{rl}f(x)& =\frac{1}{1+e^{-x}}\\ f^{\prime }(x)& =f(x)(1-f(x))\end{array}$$
这个函数可以让我们的输出值取值范围固定在 $[0,1]$ 之间。

实际上，神经网络使用的 $f(x)$ 不一定需要让输出值固定在 $[0,1]$ 之间，比如卷积神经网络(CNN)中有一个 $f(x)$ 是 $\mathrm{ReLU}$ 函数，其定义为：
$$f(x)=\{\begin{array}{ll}x& x>0\\ 0& \mathrm{otherwise}\end{array}$$

请留意 $z_i^L$ 这个值，我们暂且将其称为成为中间值。

神经网络的学习方法

整个神经网络的计算过程，可以认为是一个庞大的函数 $F$ ，其自变量的个数是输入层神经元的个数，而因变量的个数是输出层神经元的个数。而我们需要的让这个函数去拟合一个标准函数 $G$ 。

而神经网络的学习过程，实际上是一个用提供的若干个输入输出（我们称这些输入输出为训练集）不断调整权重 $w$ 以及偏移值 $b$ ，使得 $F$ 拟合 $G$ 的过程。

过拟合

实际上在神经网络的学习过程中，因为参数过多（或者其他原因），容易出现过拟合(overfitting)的现象。表现为在对于训练集中的数据，神经网络可以给出优秀的结果，但是对于训练集以外的数据，则出现很多错误。

一个形象的理解方法就是：当我们用给定的一些三次函数 $G$ 上的点，用一个 $6$ 次函数 $F$ 去拟合这个 $3$ 次函数 $G$ 的时候，对于给定的点，$F$ 都和 $G$ 的输出一致，但是本应该是 $0$ 的 $4,5,6$ 次项不全是 $0$ ，导致 $F$ 会出现一些弯曲，这时候就出现了过拟合。下图是一个过拟合的示意图：

其中绿色的点是训练集，$G(x)$ 是我们需要的函数，$F(x)$ 为神经网络拟合出的函数。

一般来说，我们会使用另一套输入输出（我们称之为测试集）来评判神经网络的学习效果。

在当前状态下，我们并不考虑过拟合问题。

反向传播方法(Backprop，简称BP)

反向传播方法是神经网络学习的核心算法。

对于训练集中的一个输入输出 $k$ ，我们用 $O_i$ 表示理想状态下，神经网络的输出层的第 $i$ 个神经元应该有的输出值。 这时候，我们可以计算出神经网络当前的函数 $F$ 和 $G$ 的误差，这里我们用类似方差的表达式表示其误差：
$$C_k=\sum _{i=0}^{N_{L_l}}(a_i^{L_l}-O_i{)}^2$$

实际上有其他表示误差的方式，本人不是很了解。

显然的，我们需要让 $C_k$ 越小越好，这时候，我们需要知道 $a_i^{L_l}$ 应该如何变化才能使得 $C_k$ 变小，可以想到运用导数计算出 $a$ 应该有的变化，也就是：
$$\frac{\partial C_k}{\partial a_i^{L_l}}=2\left(a_i^{L_l}-O_i\right)$$
想让 $C_k$ 变小，就需要让 $a_i^{L_l}$ 向着 $-\frac{\partial C_k}{a_i^{L_l}}$ 方向减小，至于如何做到这一点，我们可以调整 $b_i^{L_l},w_{i,\ast }^{L_l},a_{\ast }^{L_{l-1}}$。而调整这些参数，也可以使用求导的方式 （$0\le j<n_{l-1}$）：
$$\begin{array}{rl}\frac{\partial C_k}{\partial b_i^{L_l}}& =\frac{\partial z_i^{L_l}}{\partial b_i^{L_l}}\frac{\partial a_i^{L_l}}{\partial z_i^{L_l}}\frac{\partial C_k}{\partial a_i^{L_l}}=1\cdot f^{\prime }\left(z_i^{L_l}\right)\frac{\partial C_k}{\partial a_i^{L_l}}\\ \frac{\partial C_k}{\partial w_{i,j}^{L_l}}& =\frac{\partial z_i^{L_l}}{\partial w_{i,j}^{L_l}}\frac{\partial a_i^{L_l}}{\partial z_i^{L_l}}\frac{\partial C_k}{\partial a_i^{L_l}}=a_j^{L_l-1}{f}^{\prime }\left(z_i^{L_l}\right)\frac{\partial C_k}{\partial a_i^{L_l}}\\ \frac{\partial C_k}{\partial a_j^{L_l-1}}& =\sum _{i=0}^{n_{L_l}}\frac{\partial z_i^{L_l}}{\partial a_j^{L_l-1}}\frac{\partial a_i^{L_l}}{\partial z_i^{L_l}}\frac{\partial C_k}{\partial a_i^{L_l}}=\sum _{i=0}^{n_{L_l}-1}{w}_{i,j}^{L_l}{f}^{\prime }\left(z_i^{L_l}\right)\frac{\partial C_k}{\partial a_i^{L_l}}\end{array}$$
要调整这些参数，只需要向着导数的反方向增减这些变量即可，至于 $a_j^{L_l-1}$ 的调整，需要调整 $b_j^{L_l-1}$ 以及 $w_{j,\ast }^{L_l-1}$ ，其调整方式也是求导，且结构和上述的式子大致相同：
$$\begin{array}{rl}& \forall L\in [1,L_l-1]\cap \mathbb{Z}\\ & \frac{\partial C_k}{\partial b_i^L}=\frac{\partial z_i^L}{\partial b_i^L}\frac{\partial a_i^L}{\partial z_i^L}\frac{\partial C_k}{\partial a_i^L}=1\cdot f^{\prime }\left(z_i^L\right)\frac{\partial C_k}{\partial a_i^L}\\ & \frac{\partial C_k}{\partial w_{i,j}^L}=\frac{\partial z_i^L}{\partial w_{i,j}^L}\frac{\partial a_i^L}{\partial z_i^L}\frac{\partial C_k}{\partial a_i^L}=a_j^{L-1}{f}^{\prime }\left(z_i^L\right)\frac{\partial C_k}{\partial a_i^L}\\ & \frac{\partial C_k}{\partial a_j^L}=\sum _{i=0}^{n_{L+1}-1}\frac{\partial z_i^{L+1}}{\partial a_j^L}\frac{\partial a_i^{L+1}}{\partial z_i^{L+1}}\frac{\partial C_k}{\partial a_i^{L+1}}=\sum _{i=0}^{n_{L+1}-1}{w}_{i,j}^{L+1}{f}^{\prime }\left(z_i^{L+1}\right)\frac{\partial C_k}{\partial a_i^{L+1}}\end{array}$$
至此，我们得到了为了拟合数据 $k$ 而需要调整的数据的内容，以及调整之前的网络的信息，我们用两个数组（或者认为是矩阵，向量）将其整合一下：
$$\nabla C_k=\left[\begin{array}{c}\frac{\partial C_k}{\partial b_0^0}\\ \frac{\partial C_k}{\partial b_1^0}\\ ⋮\\ \frac{\partial C_k}{\partial b_{n_{L_l}-1}^{L_l}}\\ \frac{\partial C_k}{\partial w_{0,0}^1}\\ \frac{\partial C_k}{\partial w_{0,1}^1}\\ ⋮\\ \frac{\partial C_k}{\partial w_{n_{L_l}-1,n_{L_l-1}-1}^{L_l}}\end{array}\right],D=\left[\begin{array}{c}{b}_0^0\\ b_1^0\\ ⋮\\ b_{n_{L_l}-1}^{L_l}\\ w_{0,0}^1\\ w_{0,1}^1\\ ⋮\\ w_{n_{L_l}-1,n_{L_l-1}-1}^{L_l}\end{array}\right]$$
当我们调整网络的时候，可以使用一个参数 $\eta$ 调整这个数据对学习的影响，而调整方式，可以用 $D^{\prime }=D-\eta \cdot \nabla C_k$ 表示对这个数据的学习。可想而知 $\eta$ 越大，这个数据对网络的影响越强。而 $\eta$ 的调节，是一个需要讨论的问题，但是本人不太懂，所以本文中的 $\eta =1$ 。

Mini-Batch

如果对于每个数据，都进行一次调整 ，那么效率会变得很低，为此，我们可以使用一个小技巧，也就是 $\text{Mini-Batch}$ 。具体操作方式如下：

首先设定一个大小 $S$ ，本文中，我们设其为 $10$ 。接着，将数据集打乱，然后每 $S$ 个数据组成一个 $\mathrm{Batch}$ 。对于每一个 $\mathrm{Batch}$ ，设其中的数据编号为 $1,2,3,\dots ,S$ ，对于每一个数据 $i$ ，都用 $\text{BP}$ 计算出一个 $\nabla C_i$ ，接着将 $\mathrm{Batch}$ 中的调整信息整合，这里我们用平均值进行整合： $\nabla C=\frac{1}{S}{\sum }_{i=1}^S\nabla C_i$ ，接着用这个 $\nabla C$ 调整神经网络。

代码

$\text{NN}.h$：

#pragma once

#include <vector>
#include <cstdarg>
#include "gloconst.h"


class Network {
private:
	double **a, **b, **da, **db, **z;
	double ***w, ***dw;
	std::vector<int> network_size;
	/// <summary>
	/// 进行逆向传播
	/// </summary>
	/// <param name="output">学习的输出</param>
	void Backprop(double *output);
public:
	Network(std::vector<int> &network_size);
	~Network();
	
	/// <summary>
	/// 计算一个数据
	/// </summary>
	/// <param name="input">输入</param>
	void Calculate(double *input);
	/// <summary>
	/// 对模型进行训练
	/// </summary>
	/// <param name="data">数据 data[i][0] 表示第i个输入，data[i][1] 表示第i个输出</param>
	/// <param name="size">数据的数量</param>
	/// <param name="batch_size">mini-batch的大小</param>
	void Train(double ***data, int size, int batch_size = 100);
	/// <summary>
	/// 得到输出值a最大的节点的id
	/// </summary>
	/// <returns></returns>
	int GetMaxOutputNode();
};

$\text{NN.cpp}$：

#include "NN.h"
#include <algorithm>
#include <ctime>
#include <cassert>
#include <random>

unsigned seed = time(NULL);
std::default_random_engine gen(seed);
std::normal_distribution<double> dis(0, 1);

//生成正态分布的随机数 mu=0,sigma=1
inline double randn() {
	return dis(gen);
}
Network::Network(std::vector<int> &network_size) {
	
	this->network_size = network_size;
	a = new double *[network_size.size()];
	b = new double *[network_size.size()];
	da = new double *[network_size.size()];
	db = new double *[network_size.size()];
	z = new double *[network_size.size()];
	w = new double **[network_size.size() - 1];
	dw = new double **[network_size.size() - 1];

	for (int i = 0; i < network_size.size(); i++) { 
		int size_i = network_size[i];
		b[i] = new double[size_i], a[i] = new double[size_i], da[i] = new double[size_i], db[i] = new double[size_i], z[i] = new double[size_i];
		for (int j = 0; j < size_i; j++) b[i][j] = randn();
		if (i) {
			w[i] = new double *[size_i];
			dw[i] = new double *[size_i];
			for (int j = 0; j < size_i; j++) {
				int size2 = network_size[i - 1];
				w[i][j] = new double[size2];
				dw[i][j] = new double[size2];

				for (int k = 0; k < size2; k++)
					w[i][j][k] = randn();
			}
		}
	}
}

Network::~Network() {
	for (int i = 0; i < network_size.size(); i++) {
		int size_i = network_size[i];
		delete[] a[i], delete[] b[i], delete[] da[i], delete[] db[i], delete[] z[i];
		if (i) for (int j = 0; j < size_i; j++)
				delete[] w[i][j], delete[] dw[i][j];
	}
}

void Network::Calculate(double *input) {
	memcpy(a[0], input, sizeof(double) * network_size[0]);
	for (int i = 1; i < network_size.size(); i++) {
		int sizei = network_size[i], size_lst = network_size[i - 1];
		for (int j = 0; j < sizei; j++) {
			z[i][j] = b[i][j];
			for (int k = 0; k < size_lst; k++) z[i][j] += w[i][j][k] * a[i - 1][k];
			a[i][j] = Sigmoid(z[i][j]);
		}
	}
}

void Network::Backprop(double *output) {
	int lstid = network_size.size() - 1, sizei = network_size[network_size.size() - 1];
	for (int i = 0; i < sizei; i++) {
		da[lstid][i] = (a[lstid][i] - output[i]),
		db[lstid][i] = DSigmoid(z[lstid][i]) * da[lstid][i];
		for (int j = 0; j < network_size[lstid - 1]; j++)
			dw[lstid][i][j] = a[lstid - 1][j] * db[lstid][i];
	}
	for (int i = lstid - 1; i > 0; i--) {
		for (int j = 0; j < network_size[i]; j++) {
			da[i][j] = 0;
			for (int k = 0; k < network_size[i + 1]; k++)
				da[i][j] += w[i + 1][k][j] * db[i + 1][k];
			db[i][j] = DSigmoid(z[i][j]) * da[i][j];
			for (int k = 0; k < network_size[i - 1]; k++)
				dw[i][j][k] = a[i - 1][k] * db[i][j];
		}
	}
}

void array_add(double *dst, double *src, int sz) {
	for (int i = 0; i < sz; i++) dst[i] += src[i];
}
void array_mul(double *dst, double x, int sz) {
	for (int i = 0; i < sz; i++) dst[i] *= x;
}

void Network::Train(double ***data, int size, int batch_size) {
	batch_size = std::min(batch_size, size);
	std::random_shuffle(data, data + size);
	double **avr_db, ***avr_dw;
	avr_db = new double *[network_size.size()];
	avr_dw = new double **[network_size.size()];
	for (int i = 0; i < network_size.size(); i++) {
		avr_db[i] = new double[network_size[i]];
		avr_dw[i] = new double*[network_size[i]];
		if (i) for (int j = 0; j < network_size[i]; j++)
				avr_dw[i][j] = new double[network_size[i - 1]];
	}
	const double eta = 1;
	for (int st = 0; st < size; st += batch_size) {
		//clean the avr array
		for (int i = 0; i < network_size.size(); i++) {
			memset(avr_db[i], 0, sizeof(double) * network_size[i]);
			if (i) for (int j = 0; j < network_size[i]; j++)
					memset(avr_dw[i][j], 0, sizeof(double) * network_size[i - 1]);
		}
		for (int x = st; x < st + batch_size; x++) {
			Calculate(data[x][0]);
			Backprop(data[x][1]);
			for (int i = 0; i < network_size.size(); i++) {
				array_add(avr_db[i], db[i], network_size[i]);
				if (i) for (int j = 0; j < network_size[i]; j++)
						array_add(avr_dw[i][j], dw[i][j], network_size[i - 1]);
			}
			
		}
		for (int i = 1; i < network_size.size(); i++) {
			array_mul(avr_db[i], -eta / batch_size, network_size[i]);
			array_add(b[i], avr_db[i], network_size[i]);
			if (i) for (int j = 0; j < network_size[i]; j++) {
					array_mul(avr_dw[i][j], -eta / batch_size, network_size[i - 1]),
					array_add(w[i][j], avr_dw[i][j], network_size[i - 1]);
				}
		}
	}
	for (int i = 0; i < network_size.size(); i++) {
		delete[] avr_db[i];
		if (i) {
			for (int j = 0; j < network_size[i]; j++) delete[] avr_dw[i][j];
			delete[] avr_dw[i];
		}
	}
	delete[] avr_db, delete[] avr_dw;
}

int Network::GetMaxOutputNode() {
	int res = 0;
	double *lyr = a[network_size.size() - 1];
	int size = network_size[network_size.size() - 1];
	for (int i = 1; i < size; i++)
		res = (lyr[res] > lyr[i] ? res : i);
	return res;
}