【C语言算法】前馈神经网络实现手写数字识别

一.简介

最近在学习神经网络相关的知识,为了巩固自己对相关知识的理解,尝试使用C语言来编写一个简单的神经网络.
前馈神经网络是一个最简单的神经网络模型,每层神经元采用全连接的方式与下一层神经元相连接,信号的传递方向是单向的.如图所示:

本文主要介绍使用C语言实现前馈神经网络的方法并通过训练实现一个识别手写数字的效果.
参考了以下文章:
神经网络之手写数字识别
机器学习笔记丨神经网络的反向传播原理及过程（图文并茂+浅显易懂）
神经网络与深度学习.邱锡鹏

二.矩阵运算

神经网络依赖于矩阵运算,首先实现前馈神经网络需要使用到的矩阵运算.
新建 mymatrix.h 文件:
定义数据类型和错误类型:

typedef double m_float;

// 定义错误类型
#define M_ERRTYPE_NULL			(0)
#define M_ERRTYPE_MEMERR		(1)
#define M_ERRTYPE_PARERR		(2)
#define M_ERRTYPE_CALCERR		(3)

声明如下函数:
其中参数m=要计算的矩阵,row=矩阵的行数,column=矩阵的列数

	// 初始化矩阵为0，成功返回0
	int m_init_matrix(m_float*m, int row, int column);

	// 初始化为单位矩阵
	int m_init_identity(m_float*m, int len);

	// 初始化矩阵为随机数范围是min到max，成功返回0
	int m_init_matrix_random(m_float*m, int row, int column);

	// 复制矩阵
	int m_matrix_copy(m_float*out, m_float*m,int row,int column);

	// 打印矩阵值
	int m_printf(m_float*m, int row, int column);

	// 两个矩阵相加,out可以和a或b是同一个变量
	int m_matrix_add(m_float* out, m_float* a, m_float* b, int row, int column);

	// 两个矩阵相减,out可以和a或b是同一个变量
	int m_matrix_sub(m_float* out, m_float* a, m_float* b, int row, int column);

	// 两个矩阵相乘,column_a==row_b
	int m_matrix_multiply(m_float* out, m_float* a, m_float* b, int row_a, int column_a, int column_b);

	// a的转置与b相乘,row_a==row_b
	int m_matrix_multiply_tr_a(m_float* out, m_float* a, m_float* b, int row_a, int column_a, int column_b);

	// a与b的转置相乘,column_a==column_b
	int m_matrix_multiply_tr_b(m_float* out, m_float* a, m_float* b, int row_a, int column_a, int row_b);

	// 求矩阵的hadamard积,out可以和a或b是同一个变量
	int m_matrix_hadamard(m_float* out, m_float* a, m_float* b, int row, int column);

	// 求矩阵的转置矩阵
	int m_matrix_transpose(m_float* out, m_float* a, int row, int column);


	// 取得矩阵第i行,第j列的数据
#define m_macro_matrix_value(m,column,i,j) (m)[(i)*(column)+(j)]

1.矩阵基本操作函数

// 初始化矩阵为0，成功返回0
int m_init_matrix(m_float* m, int row, int column)
{
	for (int i = 0; i < row * column; i++)
	{
		m[i] = 0;
	}
	return M_ERRTYPE_NULL;
}



// 初始化为单位矩阵
int m_init_identity(m_float* m, int len)
{
	for (int i = 0; i < len; i++)
	{
		for (int j = 0; j < len; j++)
		{
			if (i == j)
				m[i * len + j] = 1;
			else
				m[i * len + j] = 0;
		}
	}
	return M_ERRTYPE_NULL;
}




// 初始化矩阵为随机数
int m_init_matrix_random(m_float* m, int row, int column)
{
	static int seed = 0;
	if (seed == 0) seed = (int)time(NULL);
	for (int i = 0; i < row * column; i++)
	{
		srand(seed + i);
		seed = rand();
		m[i] = seed % 10000 / 10000.0;
	}
	return M_ERRTYPE_NULL;
}


// 复制矩阵
int m_matrix_copy(m_float* out, m_float* m, int row, int column)
{
	for (int i = 0; i < row * column; i++)
	{
		out[i] = m[i];
	}
	return M_ERRTYPE_NULL;
}



// 打印矩阵值
int m_printf(m_float* m, int row, int column)
{
	printf("matrix[%d*%d]:\r\n",row,column);
	for (int i = 0; i < row; i++)
	{
		printf("[");
		for (int j = 0; j < column; j++)
			printf("%f,", m[i * column + j]);
		printf("]\r\n");
	}
	return M_ERRTYPE_NULL;
}

2.矩阵加减法

// 两个矩阵相加,out可以和a或b是同一个变量
int m_matrix_add(m_float* out, m_float* a, m_float* b, int row, int column)
{
	for (int i = 0; i < row * column; i++)
	{
		out[i] = a[i] + b[i];
	}
	return M_ERRTYPE_NULL;
}

// 两个矩阵相减,out可以和a或b是同一个变量
int m_matrix_sub(m_float* out, m_float* a, m_float* b, int row, int column)
{
	for (int i = 0; i < row * column; i++)
	{
		out[i] = a[i] - b[i];
	}
	return M_ERRTYPE_NULL;
}

3.矩阵乘法

// 两个矩阵相乘
int m_matrix_multiply(m_float* out, m_float* a, m_float* b, int row_a, int column_a, int column_b)
{
	m_float value;
	for (int i = 0; i < row_a; i++)
	{
		for (int j = 0; j < column_b; j++)
		{
			value = 0;
			for (int k = 0; k < column_a; k++)
			{
				value += m_macro_matrix_value(a, column_a, i, k) *
					m_macro_matrix_value(b, column_b, k, j);
			}
			m_macro_matrix_value(out, column_b, i, j) = value;
		}
	}
	return M_ERRTYPE_NULL;
}



// a的转置与b相乘
int m_matrix_multiply_tr_a(m_float* out, m_float* a, m_float* b, int row_a, int column_a, int column_b)
{
	m_float value;
	for (int i = 0; i < column_a; i++)
	{
		for (int j = 0; j < column_b; j++)
		{
			value = 0;
			for (int k = 0; k < row_a; k++)
			{
				value += m_macro_matrix_value(a, column_a, k, i) *
					m_macro_matrix_value(b, column_b, k, j);
			}
			m_macro_matrix_value(out, column_b, i, j) = value;
		}
	}
	return M_ERRTYPE_NULL;
}



// a与b的转置相乘
int m_matrix_multiply_tr_b(m_float* out, m_float* a, m_float* b, int row_a, int column_a, int row_b)
{
	m_float value;
	for (int i = 0; i < row_a; i++)
	{
		for (int j = 0; j < row_b; j++)
		{
			value = 0;
			for (int k = 0; k < column_a; k++)
			{
				value += m_macro_matrix_value(a, column_a, i, k) *
					m_macro_matrix_value(b, column_a, j, k);
			}
			m_macro_matrix_value(out, row_b, i, j) = value;
		}
	}
	return M_ERRTYPE_NULL;
}




// 求矩阵的hadamard积,out可以和a或b是同一个变量
int m_matrix_hadamard(m_float* out, m_float* a, m_float* b, int row, int column)
{
	for (int i = 0; i < row * column; i++)
	{
		out[i] = a[i] * b[i];
	}
	return M_ERRTYPE_NULL;
}

三.神经网络实现

1.定义神经网络类

根据本人的理解,输入层只有自变量,所以定义神经网络类的时候没有把输入层计算在内,方便程序编写.
每一层神经元包含的变量有输入x;偏置b;参数w;输出a,每一层的输出a同时也是下一层的输入x.如果神经网络不具备学习功能,仅用这些参数就够了.
每层神经元用于学习的变量有净输出值z,这个值用于计算激活函数的梯度;损失值loss,这个变量用于反向传播时把误差向前传播.
变量loss是网络的整体误差,此误差越小越好.

// 定义神经网络最大层数(不包括输入层)
#define NERVE_MAX_LAYER_NUM		(10)

// 定义一个神经网络
typedef struct {
	// 网络层数,不包括输入层
	int layer_num;

	// 数据源个数
	// 第0层神经元接收的是外部输入数据,
	// 第i(i>0)层神经元接收的是前一层神经元的输出数据
	int src_num[NERVE_MAX_LAYER_NUM+1];

	// 外部输入数据,初始化时分配内存
	m_float* src_data;

	// 矫正值,用于计算误差,初始化时分配内存
	m_float* y;

	// 各层神经元数目
	int layer_cell_num[NERVE_MAX_LAYER_NUM];

	// 神经元输入数据,此变量不申请内存,
	// 第0层指向src_table,第i(i>0)层指向a[i-1]
	m_float* src_table[NERVE_MAX_LAYER_NUM+1];

	// 各层神经元在激活函数之前的输出值,在初始化时分配内存
	m_float* z_table[NERVE_MAX_LAYER_NUM];

	// 各层神经元的输出数据,在初始化时分配内存
	m_float* a_table[NERVE_MAX_LAYER_NUM];

	// 各层神经元的偏置,在初始化时分配内存
	m_float* b_table[NERVE_MAX_LAYER_NUM];

	// 各层神经元的权重,在初始化时分配内存
	// 对于第i层来说w是一个[layer_cell_num[i]*src_num[i]]的二维数组
	// 其中第j行存储的是第j个神经元对数据源的权重
	m_float* w_table[NERVE_MAX_LAYER_NUM];

	// 各层神经元的误差,在初始化时分配内存
	m_float* loss_table[NERVE_MAX_LAYER_NUM];

	m_float loss;// 误差
}nerve_obj;

2.激活函数ReLU

ReLU(x)=max(0,x),此函数在x小于0时恒等于0,在x>0时为x本身.
经我测试,激活函数ReLU比sigmoid学习效果要好,但斜率不能太大,猜测可能是太大的斜率产生了过大的数据,使之后的误差评价函数产生了指数爆炸.
所以在此给了0.1的斜率,实测没有产生指数爆炸.

// 激活函数
static int nerve_relu(m_float *a,m_float *z,int len)
{
	for (int i = 0; i < len; i++)
	{
		a[i] = z[i] > 0.0 ? z[i]*0.1 : 0;
	}
	return 0;
}


// 激活函数的导函数乘以一个数l,z是自变量
static int nerve_deriv_relu(m_float* out, m_float* l,m_float* z, int len)
{
	for (int i = 0; i < len; i++)
	{
		out[i] = z[i] > 0 ? l[i]*0.1 : 0;
	}
	return 0;
}

3.激活函数Softmax

Softmax函数用于评价输出y=c的概率.
Softmax(y=c|x[1…n])=exp(x[c])/(exp(x[1])+…+exp(x[n])).
在编码中使用了对数公式,在一定程度上避免了指数爆炸问题.

// softmax 函数,也可以看作是最后一层的激活函数
static m_float nerve_softmax(m_float *out,m_float *a,int len)
{
	m_float sum = 0,max;
	max = a[0];
	for (int i = 0; i < len; i++)
	{
		if (max < a[i])
			max = a[i];
	}
	for (int i = 0; i < len; i++)
	{
		sum += exp(a[i]-max);
	}
	for (int i = 0; i < len; i++)
	{
		out[i] = exp(a[i]-max) / sum;
	}
	return sum;
}

4.交叉熵损失函数

手写数字识别是分类问题,使用交叉熵损失函数来评价整体误差.
其中变量y为校正值,a为输出值,
对于手写数字识别任务来说,校正值只可能有一个为1,其余都为0,假设这个为1的项是第c项,则恒有误差loss=-log(a[c]).

// 误差计算函数
static m_float nerve_loss(m_float *y,m_float *a,int len)
{
	m_float t = 0;
	for (int i = 0; i < len; i++)
	{
		if(y[i]>0)
			t+= y[i] * log(a[i]);
	}
	return -t;
}

5.神经网络基本操作

实现神经网络初始化,输入输出接口等操作.其中nerve_get_max_index输出的是输出矩阵中最大值的序号,对于手写数字识别任务来说,此序号就是识别的数字.

// 初始化神经网络
int nerve_init(nerve_obj* n, int layer, int src_num, int* cell_nums)
{
	memset(n, 0, sizeof(nerve_obj));

	n->layer_num = layer;
	n->src_num[0] = src_num;
	n->src_data = calloc(n->src_num[0], sizeof(m_float));
	n->y = calloc(cell_nums[n->layer_num - 1], sizeof(m_float));
	n->src_table[0] = n->src_data;
	for (int i = 0; i < n->layer_num; i++)
	{
		n->layer_cell_num[i] = cell_nums[i];
		n->src_num[i+1]= cell_nums[i];
		n->z_table[i] = calloc(n->layer_cell_num[i], sizeof(m_float));
		n->a_table[i] = calloc(n->layer_cell_num[i], sizeof(m_float));
		n->b_table[i] = calloc(n->layer_cell_num[i], sizeof(m_float));
		n->loss_table[i] = calloc(n->layer_cell_num[i], sizeof(m_float));
		n->w_table[i] = calloc(n->layer_cell_num[i] * n->src_num[i], sizeof(m_float));
		m_init_matrix_random(n->w_table[i], n->layer_cell_num[i], n->src_num[i]);
		m_init_matrix_random(n->b_table[i], 1, n->layer_cell_num[i]);
		n->src_table[i + 1] = n->a_table[i];
	}
	return 0;
}



// 获取输入缓冲区
m_float* nerve_get_input_ptr(nerve_obj* n)
{
	return n->src_data;
}

// 获取输出缓冲区
m_float* nerve_get_output_ptr(nerve_obj* n)
{
	return n->a_table[n->layer_num-1];
}

// 获取校正值缓冲区
m_float* nerve_get_y_ptr(nerve_obj* n)
{
	return n->y;
}

// 获取最大输出值的序号
int nerve_get_max_index(nerve_obj* n)
{
	int index = 0;
	m_float* a = nerve_get_output_ptr(n);
	m_float max=a[0];
	for (int i = 0; i < n->layer_cell_num[n->layer_num - 1]; i++)
	{
		if (max < a[i])
		{
			max = a[i];
			index = i;
		}
	}
	return index;
}

6.正向传播计算

正向传播就是从第0层开始,依次计算z=w*x+b;a=f(z).其中f(.)为激活函数.
由于最后一层的激活函数为softmax,所以分开计算.

// 计算神经网络
int nerve_calc(nerve_obj* n)
{
	for (int i = 0; i < n->layer_num; i++)
	{
		m_matrix_multiply(n->z_table[i], n->w_table[i], n->src_table[i],
			n->layer_cell_num[i], n->src_num[i], 1);
		m_matrix_add(n->z_table[i], n->z_table[i], n->b_table[i], n->layer_cell_num[i], 1);
		if(i<n->layer_num-1)
			nerve_relu(n->a_table[i], n->z_table[i], n->layer_cell_num[i]);
		else
			nerve_softmax(n->a_table[i], n->z_table[i], n->layer_cell_num[i]);
	}
}

7.反向传播学习

对于交叉熵评价函数,其误差loss=(a-y),若是其他评价函数,需要对此函数进行修改.

// 计算最后一层的误差2L(z)/2(z)
static int nerve_last_loss(m_float *out,m_float *a,m_float *y,int len)
{
	m_matrix_sub(out, a, y, 1, len);
	return 0;
}

在本程序中学习率a只能取负值.
最后一层的误差使用公式loss=(a-y)计算得出;
其余层的误差使用公式loss[n]=f’(z)(tr_wloss[n+1])得出,其中tr_w为w的转置矩阵,f’(.)为激活函数的导函数.
更新权重使用公式w=w+alosstr_x;b=b+a*loss.其中tr_x为x的转置矩阵.
在编码中可以使用函数m_matrix_multiply_tr_b来计算新权重,这里使用的多重循环.

// 反向传播学习,a学习率
int nerve_calc_bp(nerve_obj* n,m_float a)
{
	m_float* w;

	// 计算预测误差
	n->loss=nerve_loss(n->y, nerve_get_output_ptr(n),n->layer_cell_num[n->layer_num-1]);

	// 计算最后一层的误差
	nerve_last_loss(n->loss_table[n->layer_num - 1], nerve_get_output_ptr(n),n->y, n->layer_cell_num[n->layer_num - 1]);
	
	// 计算其余层的误差
	for (int i = n->layer_num-1; i > 0; i--)
	{
		// 第i层的误差权重和
		m_matrix_multiply_tr_a(n->loss_table[i - 1], n->w_table[i], n->loss_table[i],
			n->layer_cell_num[i], n->src_num[i], 1);
		// 使用激活梯度函数可能是因为梯度消失的原因不能收敛
		nerve_deriv_relu(n->loss_table[i - 1], n->loss_table[i - 1], n->z_table[i - 1], n->layer_cell_num[i - 1]);
	}

	// 更新权重
	for (int i = 0; i < n->layer_num; i++)
	{
		w = n->w_table[i];
		for (int j = 0; j < n->layer_cell_num[i]; j++)
		{
			for (int k = 0; k < n->src_num[i]; k++)
			{
				m_macro_matrix_value(w, n->src_num[i],j,k)+= 
					a * n->src_table[i][k] * n->loss_table[i][j];
			}
			n->b_table[i][j] += a * n->loss_table[i][j];
		}
	}
}

四.测试

1.统计最近n次的平均识别误差

// 统计最近n次的误差
#define TRAIN_LOSS_NUM	(200)

// 计算最近TRAIN_LOSS_NUM次的误差
m_float calc_loss(m_float loss)
{
	static int num = 0, ptr=0;
	static m_float loss_table[TRAIN_LOSS_NUM];
	m_float loss_ave, loss_sum;

	if (num < TRAIN_LOSS_NUM)
	{
		loss_table[num] = loss;
		num++;
	}
	else
	{
		loss_table[ptr] = loss;
		ptr++;
		if (ptr >= num) ptr = 0;
	}
	loss_sum = 0;
	for (int i = 0; i < num; i++)
		loss_sum += loss_table[i];
	loss_ave = loss_sum / num;
	return loss_ave;
}

2.定义监督学习的标签

m_float drain_dst_table[10][10] = {
	{ 1,0,0,0,0,0,0,0,0,0 },
	{ 0,1,0,0,0,0,0,0,0,0 },
	{ 0,0,1,0,0,0,0,0,0,0 },
	{ 0,0,0,1,0,0,0,0,0,0 },
	{ 0,0,0,0,1,0,0,0,0,0 },
	{ 0,0,0,0,0,1,0,0,0,0 },
	{ 0,0,0,0,0,0,1,0,0,0 },
	{ 0,0,0,0,0,0,0,1,0,0 },
	{ 0,0,0,0,0,0,0,0,1,0 },
	{ 0,0,0,0,0,0,0,0,0,1 },
};

3.编写主函数

int main()
{
	// 神经网络类
	nerve_obj nerve = { 0 };
	// 输入,输出,校正值
	m_float* bmp_x = 0, * out = 0, * z_true;
	// 训练集
	const char* train_data, * train_p;
	// 测试集
	const char* test_data, * test_p;
	// 两层神经网络,第一层88个神经元,第二层10个神经元
	int cells[] = { 2,88,10 };
	// 序号1,序号2,测试成功的次数
	int index, index2, test_true;
	// 加载数据集
	train_data = data_read_file(".\\hand\\mnist_train.csv");
	test_data = data_read_file(".\\hand\\mnist_test.csv");
	// 初始化神经网络
	nerve_init(&nerve, cells[0], 28 * 28, &cells[1]);
	bmp_x = nerve_get_input_ptr(&nerve);
	out = nerve_get_output_ptr(&nerve);
	z_true = nerve_get_y_ptr(&nerve);
	printf("初始化完成\r\n");

	for (int m = 0; m < 1; m++)
	{
		train_p = train_data;
		for (int i = 0; i < 59999; i++)
		{
			// 读取一行数据,也就是一个训练图片
			index = data_tr_buff(bmp_x, train_p);
			train_p = data_get_next_line(train_p);
			memcpy(z_true, drain_dst_table[index], sizeof(m_float) * 10);
			// 计算神经网络
			nerve_calc(&nerve);
			// 使用校正值来学习
			nerve_calc_bp(&nerve, -0.05);
			// 显示训练次数和误差,正常情况下误差loss会越来越小
			printf("time=%d,loss=%f           \r", i, calc_loss(nerve.loss));

		}
		printf("\n训练完成\r\n");



		test_p = test_data;
		test_true = 0;
		for (int i = 0; i < 9999; i++)
		{
			// 读取一行数据
			index = data_tr_buff(bmp_x, test_p);
			test_p = data_get_next_line(test_p);
			// 计算神经网络
			nerve_calc(&nerve);
			// 获取识别的数字
			index2 = nerve_get_max_index(&nerve);
			// 如果识别的数字与标签数字相等,则识别成功
			//if ((index == index2)&&(out[index2]>0.85))
			if (index == index2)
				test_true++;
			printf("time=%d           \r", i);
		}
		printf("测试识别正确 %d 次\r\n", test_true);
	}
	return 0;
}

4.测试结果

对于60000个数据的训练集,1000个数据的测试集mnist来说,训练之后的测试成功次数基本都在9000次以上,也就是识别成功率在90%以上.