BP神经网络公式推导及实现(MNIST)

BP神经网络的基础介绍见：http://blog.csdn.net/fengbingchun/article/details/50274471，这里主要以公式推导为主。

BP神经网络又称为误差反向传播网络，其结构如下图。这种网络实质是一种前向无反馈网络，具有结构清晰、易实现、计算功能强大等特点。



BP神经网络有一个输入层，一个输出层，一个或多个隐含层。每一层上包含了若干个节点，每个节点代表一个神经元，同一层上各节点之间无任何耦合连接关系，层间各神经元之间实现全连接，即后一层(如输入层)的每一个神经元与前一层(如隐含层)的每一个神经元实现全连接。网络按照监督学习的方式学习，当信息被输入网络后神经元受到刺激，激活值从输入层依次经过各隐含层节点，最后在输出层的各节点获得网络的输入响应。

BP神经网络的基本思想：BP神经网络的学习采用误差反向传播算法，BP算法是一种有监督的学习方法，其主要思想是把整个学习过程分为正向传播、反向(逆向)传播和记忆训练三个部分。正向传播时，输入样本从输入层输入，经各隐含层处理后传向输出层，每一层神经元的状态只影响下一层神经元的状态。如果在输出层得不到期望的输出，则转入误差的反向传播阶段，将输出误差以某种形式通过隐含层向输入层反传，并将误差分摊给各层的所有单元，从而获得各层单元的误差信号并将其作为修正各单元权值的依据。这种网络的信号正向传播与误差反向传播是反复交替进行的，权值的不断调整就是网络的记忆训练过程。网络的记忆训练过程一直进行到网络趋向收敛，即输出误差达到要求的标准。

三层BP神经网络的学习算法：为了使BP网络具有某种功能，完成某项任务，必须调整层间连接权值和节点阈值，使所有样品的实际输出和期望输出之间的误差稳定在一个较小的值之内。三层BP网络学习过程主要由四部分组成：(1)、输入模式顺传播(输入模式由输入层经隐含层向输出层传播计算)；(2)、输出误差逆传播(输出的误差由输出层经隐含层传向输入层)；(3)、循环记忆训练(模式顺传播与误差逆传播的计算过程反复交替循环进行)；(4)、学习结果判别(判定全局误差是否趋向极小值或是否已达到设定的最大迭代次数)。

(1)、输入模式顺传播：这一过程主要是利用输入模式求出它所对应的实际输出。

确定输入向量Xk：



式中，k=1,2,…,m；m是学习模式对数(训练模式对数)；n是输入层单元数。

确定期望输出向量Yk：



式中，k=1,2,…,m；m是学习模式对数(训练模式对数)；q为输出层单元数。

计算隐含层各神经元的激活值sj：



式中，n是输入层单元数；wij是输入层至隐含层的连接权值；θj是隐含层单元的阈值；j=1,2…p，p是隐含层单元数。

激活函数采用s型函数：



计算隐含层j单元的输出值：将上面的激活值即公式(3)代入激活函数即公式(4)中可得隐含层j单元的输出值：



阈值θj在学习过程中与权值wij一样也不断地被修正。

计算输出层第t个单元的激活值ot：



计算输出层第t个单元的实际输出值ct：



式中，wjt是隐含层至输出层的权值；θt是输出层单元阈值；j=1,2…p，p是隐含层单元数；xj为隐含层第j个节点的输出值；f是s型激活函数，t=1,2…,q，q为输出层单元数。

利用以上各公式就可以计算出一个输入模式的顺传播过程。

(2)、输出误差的逆传播：在第一步的模式顺传播计算中得到了网络的实际输出值，当这些实际的输出值与希望的输出值不一样或者误差大于所限定的数值时，就要对网络进行校正。

这里的校正是从前往后进行的，所以叫做误差逆传播，计算时是从输出层到隐含层，再从隐含层到输入层。

输出层的校正误差：



式中，t=1,2,…,q，q是输出层单元数；k=1,2,…,m，m是训练(学习)模式对数；ytk是希望输出；ctk是实际输出；f’(.)是对输出函数的导数。

隐含层各单元的校正误差：



式中，t=1,2,…,q，q是输出层单元数；j=1,2,…,p; p是隐含层单元数；k=1,2,…,m，m是训练(学习)模式对数。

对于输出层至隐含层连接权和输出层阈值的校正量：



式中，bjk是隐含层j单元的输出；dtk是输出层的校正误差；j=1,2…,p；t=1,2,…,q；k=1,2,…,m； α>0(输出层至隐含层学习率)。

隐含层至输入层的校正量：



式中，ejk是隐含层j单元的校正误差；xik是标准输入，i=1,2,…,n ，n是输入层单元数；0<β<1(隐含层至输入层学习率)。

(3)、循环记忆训练：为使网络的输出误差趋向于极小值，对于BP网输入的每一组训练模式，一般要经过数百次甚至上万次的循环记忆训练，才能使网络记住这一模式。这种循环记忆实际上就是反复重复上面介绍的输入模式顺传播和输出误差逆传播。

(4)、学习结果的判别：当每次循环记忆训练结束后，都要进行学习结果的判别。判别的目的主要是检查输出误差是否已经小到可以允许的程度。如果小到可以允许的程度，就可以结束整个学习过程，否则还要继续进行循环训练。

确定隐含层节点数：一般有3个经验公式：



式中，m为要设置的隐含层节点数；n为输入层节点数；l为输出层节点数；α为1至10之间的常数。

以下按照上面的公式实现的BP，通过MNIST库测试，识别率可以达到96.5%以上。

BP.hpp:

#ifndef _BP_HPP_
#define _BP_HPP_

namespace ANN {

#define num_node_input_BP	784 //输入层节点数
#define width_image_BP		28 //归一化图像宽
#define height_image_BP		28 //归一化图像高
#define num_node_hidden_BP	120 //隐含层节点数
#define num_node_output_BP	10 //输出层节点数
#define alpha_learning_BP	0.8 //输出层至隐含层学习率
#define beta_learning_BP	0.6 //隐含层至输入层学习率
#define patterns_train_BP	60000 //训练模式对数(总数)
#define patterns_test_BP	10000 //测试模式对数(总数)
#define iterations_BP		10000 //最大训练次数
#define accuracy_rate_BP	0.965 //要求达到的准确率

class BP {
public:
BP();
~BP();

void init(); //初始化，分配空间
bool train(); //训练
int predict(const int* data, int width, int height); //预测
bool readModelFile(const char* name); //读取已训练好的BP model

protected:
void release(); //释放申请的空间
bool saveModelFile(const char* name); //将训练好的model保存起来，包括各层的节点数，权值和阈值
bool initWeightThreshold(); //初始化，产生[-1, 1]之间的随机小数
bool getSrcData(); //读取MNIST数据
void calcHiddenLayer(const int* data); //计算隐含层输出
void calcOutputLayer(); //计算输出层输出
void calcAdjuctOutputLayer(const int* data); //计算输出层校正误差
void calcAdjuctHiddenLayer(); //计算隐含层校正误差
float calcActivationFunction(float x); //计算激活函数，对数S形函数
void updateWeightThresholdOutputLayer(); //更新输出层至隐含层权值和阈值
void updateWeightThresholdHiddenLayer(const int* data); //更新隐含层至输入层权值和阈值
float test(); //训练完一次计算一次准确率

private:
float weight1[num_node_input_BP][num_node_hidden_BP]; //输入层至隐含层连接权值
float weight2[num_node_hidden_BP][num_node_output_BP]; //隐含层至输出层连接权值
float threshold1[num_node_hidden_BP]; //隐含层阈值
float threshold2[num_node_output_BP]; //输出层阈值
float output_hiddenLayer[num_node_hidden_BP]; //顺传播，隐含层输出值
float output_outputLayer[num_node_output_BP]; //顺传播，输出层输出值
float adjust_error_outputLayer[num_node_output_BP]; //逆传播，输出层校正误差
float adjust_error_hiddenLayer[num_node_hidden_BP]; //逆传播，隐含层校正误差

int* data_input_train; //原始标准输入数据，训练
int* data_output_train; //原始标准期望结果，训练
int* data_input_test; //原始标准输入数据，测试
int* data_output_test; //原始标准期望结果，测试
};

}

#endif //_BP_HPP_

BP.cpp:

#include <assert.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <windows.h>

#include "BP.hpp"

namespace ANN {

BP::BP()
{
data_input_train = NULL;
data_output_train = NULL;
data_input_test = NULL;
data_output_test = NULL;
}

BP::~BP()
{
release();
}

void BP::release()
{
if (data_input_train) {
delete[] data_input_train;
}
if (data_output_train) {
delete[] data_output_train;
}
if (data_input_test) {
delete[] data_input_test;
}
if (data_output_test) {
delete[] data_output_test;
}
}

bool BP::initWeightThreshold()
{
srand(time(0) + rand());

for (int i = 0; i < num_node_input_BP; i++) {
for (int j = 0; j < num_node_hidden_BP; j++) {
weight1[i][j] = -1 + 2 * ((float)rand()) / RAND_MAX; //[-1, 1]
}
}

for (int i = 0; i < num_node_hidden_BP; i++) {
for (int j = 0; j < num_node_output_BP; j++) {
weight2[i][j] = -1 + 2 * ((float)rand()) / RAND_MAX;
}
}

for (int i = 0; i < num_node_hidden_BP; i++) {
threshold1[i] = -1 + 2 * ((float)rand()) / RAND_MAX;
}

for (int i = 0; i < num_node_output_BP; i++) {
threshold2[i] = -1 + 2 * ((float)rand()) / RAND_MAX;
}

return true;
}

static int reverseInt(int i)
{
unsigned char ch1, ch2, ch3, ch4;
ch1 = i & 255;
ch2 = (i >> 8) & 255;
ch3 = (i >> 16) & 255;
ch4 = (i >> 24) & 255;
return((int)ch1 << 24) + ((int)ch2 << 16) + ((int)ch3 << 8) + ch4;
}

static void readMnistImages(std::string filename, int* data_dst, int num_image)
{
std::ifstream file(filename, std::ios::binary);
assert(file.is_open());

int magic_number = 0;
int number_of_images = 0;
int n_rows = 0;
int n_cols = 0;
file.read((char*)&magic_number, sizeof(magic_number));
magic_number = reverseInt(magic_number);
file.read((char*)&number_of_images, sizeof(number_of_images));
number_of_images = reverseInt(number_of_images);
assert(number_of_images == num_image);
file.read((char*)&n_rows, sizeof(n_rows));
n_rows = reverseInt(n_rows);
file.read((char*)&n_cols, sizeof(n_cols));
n_cols = reverseInt(n_cols);
assert(n_rows == height_image_BP && n_cols == width_image_BP);

for (int i = 0; i < number_of_images; ++i) {
for (int r = 0; r < n_rows; ++r) {
for (int c = 0; c < n_cols; ++c) {
unsigned char temp = 0;
file.read((char*)&temp, sizeof(temp));
//data_dst[i * num_node_input_BP + r * n_cols + c] = (int)temp; //formula[1]
if (temp > 128) {
data_dst[i * num_node_input_BP + r * n_cols + c] = 1;
} else {
data_dst[i * num_node_input_BP + r * n_cols + c] = 0;
}
}
}
}
}

static void readMnistLabels(std::string filename, int* data_dst, int num_image)
{
std::ifstream file(filename, std::ios::binary);
assert(file.is_open());

int magic_number = 0;
int number_of_images = 0;
file.read((char*)&magic_number, sizeof(magic_number));
magic_number = reverseInt(magic_number);
file.read((char*)&number_of_images, sizeof(number_of_images));
number_of_images = reverseInt(number_of_images);
assert(number_of_images == num_image);

for (int i = 0; i < number_of_images; ++i) {
unsigned char temp = 0;
file.read((char*)&temp, sizeof(temp));
data_dst[i * num_node_output_BP + temp] = 1; //formula[2]
}
}

bool BP::getSrcData()
{
assert(data_input_train && data_output_train && data_input_test && data_output_test);

std::string filename_train_images = "D:/Download/MNIST/train-images.idx3-ubyte";
std::string filename_train_labels = "D:/Download/MNIST/train-labels.idx1-ubyte";
readMnistImages(filename_train_images, data_input_train, patterns_train_BP);
/*unsigned char* p = new unsigned char[784];
memset(p, 0, sizeof(unsigned char) * 784);
for (int j = 0, i = 59998 * 784; j< 784; j++, i++) {
p[j] = (unsigned char)data_input_train[i];
}
delete[] p;*/
readMnistLabels(filename_train_labels, data_output_train, patterns_train_BP);
/*int* q = new int[10];
memset(q, 0, sizeof(int) * 10);
for (int j = 0, i = 59998 * 10; j < 10; j++, i++) {
q[j] = data_output_train[i];
}
delete[] q;*/

std::string filename_test_images = "D:/Download/MNIST/t10k-images.idx3-ubyte";
std::string filename_test_labels = "D:/Download/MNIST/t10k-labels.idx1-ubyte";
readMnistImages(filename_test_images, data_input_test, patterns_test_BP);
readMnistLabels(filename_test_labels, data_output_test, patterns_test_BP);

return true;
}

void BP::init()
{
data_input_train = new int[patterns_train_BP * num_node_input_BP];
memset(data_input_train, 0, sizeof(int) * patterns_train_BP * num_node_input_BP);
data_output_train = new int[patterns_train_BP * num_node_output_BP];
memset(data_output_train, 0, sizeof(int) * patterns_train_BP * num_node_output_BP);
data_input_test = new int[patterns_test_BP * num_node_input_BP];
memset(data_input_test, 0, sizeof(int) * patterns_test_BP * num_node_input_BP);
data_output_test = new int[patterns_test_BP * num_node_output_BP];
memset(data_output_test, 0, sizeof(int) * patterns_test_BP * num_node_output_BP);

initWeightThreshold();
getSrcData();
}

float BP::calcActivationFunction(float x)
{
return 1.0 / (1.0 + exp(-x)); //formula[4] formula[5] formula[7]
}

void BP::calcHiddenLayer(const int* data)
{
for (int i = 0; i < num_node_hidden_BP; i++) {
float tmp = 0;
for (int j = 0; j < num_node_input_BP; j++) {
tmp += data[j] * weight1[j][i];
}

tmp -= threshold1[i]; //formula[3]
output_hiddenLayer[i] = calcActivationFunction(tmp);
}
}

void BP::calcOutputLayer()
{
for (int i = 0; i < num_node_output_BP; i++) {
float tmp = 0;
for (int j = 0; j < num_node_hidden_BP; j++) {
tmp += output_hiddenLayer[j] * weight2[j][i];
}

tmp -= threshold2[i]; //formula[6]
output_outputLayer[i] = calcActivationFunction(tmp);
}
}

void BP::calcAdjuctOutputLayer(const int* data)
{
for (int i = 0; i < num_node_output_BP; i++) {
adjust_error_outputLayer[i] = (data[i] - output_outputLayer[i]) *
output_outputLayer[i] * (1.0 - output_outputLayer[i]); //formula[8], f'(x)= f(x)*(1. - f(x))
}
}

void BP::calcAdjuctHiddenLayer()
{
for (int i = 0; i < num_node_hidden_BP; i++) {
float tmp = 0;
for (int j = 0; j < num_node_output_BP; j++) {
tmp += weight2[i][j] * adjust_error_outputLayer[j];
}

adjust_error_hiddenLayer[i] = tmp * (output_hiddenLayer[i] * (1.0 - output_hiddenLayer[i])); //formula[9]
}
}

void BP::updateWeightThresholdOutputLayer()
{
for (int i = 0; i < num_node_output_BP; i++) {
for (int j = 0; j < num_node_hidden_BP; j++) {
weight2[j][i] += (alpha_learning_BP * adjust_error_outputLayer[i] * output_hiddenLayer[j]); //formula[10]
}

threshold2[i] += (alpha_learning_BP * adjust_error_outputLayer[i]); //formula[11]
}
}

void BP::updateWeightThresholdHiddenLayer(const int* data)
{
for (int i = 0; i < num_node_hidden_BP; i++) {
for (int j = 0; j < num_node_input_BP; j++) {
weight1[j][i] += (beta_learning_BP * adjust_error_hiddenLayer[i] * data[j]); //formula[12]
}

threshold1[i] += (beta_learning_BP * adjust_error_hiddenLayer[i]); //formula[13]
}
}

float BP::test()
{
int count_accuracy = 0;

for (int num = 0; num < patterns_test_BP; num++) {
int* p1 = data_input_test + num * num_node_input_BP;
calcHiddenLayer(p1);
calcOutputLayer();

float max_value = -9999;
int pos = -1;

for (int i = 0; i < num_node_output_BP; i++) {
if (output_outputLayer[i] > max_value) {
max_value = output_outputLayer[i];
pos = i;
}
}

int* p2 = data_output_test + num * num_node_output_BP;
if (p2[pos] == 1) {
count_accuracy++;
}
Sleep(1);
}

return (count_accuracy * 1.0 / patterns_test_BP);
}

bool BP::saveModelFile(const char* name)
{
FILE* fp = fopen(name, "wb");
if (fp == NULL) {
return false;
}

int num_node_input = num_node_input_BP;
int num_node_hidden = num_node_hidden_BP;
int num_node_output = num_node_output_BP;
fwrite(&num_node_input, sizeof(int), 1, fp);
fwrite(&num_node_hidden, sizeof(int), 1, fp);
fwrite(&num_node_output, sizeof(int), 1, fp);
fwrite(weight1, sizeof(weight1), 1, fp);
fwrite(threshold1, sizeof(threshold1), 1, fp);
fwrite(weight2, sizeof(weight2), 1, fp);
fwrite(threshold2, sizeof(threshold2), 1, fp);

fflush(fp);
fclose(fp);

return true;
}

bool BP::readModelFile(const char* name)
{
FILE* fp = fopen(name, "rb");
if (fp == NULL) {
return false;
}

int num_node_input, num_node_hidden, num_node_output;

fread(&num_node_input, sizeof(int), 1, fp);
assert(num_node_input == num_node_input_BP);
fread(&num_node_hidden, sizeof(int), 1, fp);
assert(num_node_hidden == num_node_hidden_BP);
fread(&num_node_output, sizeof(int), 1, fp);
assert(num_node_output == num_node_output_BP);
fread(weight1, sizeof(weight1), 1, fp);
fread(threshold1, sizeof(threshold1), 1, fp);
fread(weight2, sizeof(weight2), 1, fp);
fread(threshold2, sizeof(threshold2), 1, fp);

fflush(fp);
fclose(fp);

return true;
}

int BP::predict(const int* data, int width, int height)
{
assert(data && width == width_image_BP && height == height_image_BP);

const int* p = data;
calcHiddenLayer(p);
calcOutputLayer();

float max_value = -9999;
int ret = -1;

for (int i = 0; i < num_node_output_BP; i++) {
if (output_outputLayer[i] > max_value) {
max_value = output_outputLayer[i];
ret = i;
}
}

return ret;
}

bool BP::train()
{
int i = 0;
for (i = 0; i < iterations_BP; i++) {
std::cout << "iterations : " << i;

float accuracyRate = test();
std::cout << ",    accuray rate: " << accuracyRate << std::endl;
if (accuracyRate > accuracy_rate_BP) {
saveModelFile("bp.model");
std::cout << "generate bp model" << std::endl;
break;
}

for (int j = 0; j < patterns_train_BP; j++) {
int* p1 = data_input_train + j * num_node_input_BP;
calcHiddenLayer(p1);
calcOutputLayer();

int* p2 = data_output_train + j * num_node_output_BP;
calcAdjuctOutputLayer(p2);
calcAdjuctHiddenLayer();

updateWeightThresholdOutputLayer();
int* p3 = data_input_train + j * num_node_input_BP;
updateWeightThresholdHiddenLayer(p3);
}
}

if (i == iterations_BP) {
saveModelFile("bp.model");
std::cout << "generate bp model" << std::endl;
}

return true;
}

}

test.cpp:

#include <iostream>
#include "BP.hpp"
#include <opencv2/opencv.hpp>

int test_BP();

int main()
{
test_BP();
std::cout << "ok!" << std::endl;
}

int test_BP()
{
//1. bp train
ANN::BP bp1;
bp1.init();
bp1.train();

//2. bp predict
ANN::BP bp2;
bool flag = bp2.readModelFile("bp.model");
if (!flag) {
std::cout << "read bp model error" << std::endl;
return -1;
}

int target[10] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
std::string path_images = "../../../../test-images/";

int* data_image = new int[width_image_BP * height_image_BP];

for (int i = 0; i < 10; i++) {
char ch[15];
sprintf(ch, "%d", i);
std::string str;
str = std::string(ch);
str += ".jpg";
str = path_images + str;

cv::Mat mat = cv::imread(str, 2 | 4);
if (!mat.data) {
std::cout << "read image error" << std::endl;
return -1;
}

if (mat.channels() == 3) {
cv::cvtColor(mat, mat, cv::COLOR_BGR2GRAY);
}

if (mat.cols != width_image_BP || mat.rows != height_image_BP) {
cv::resize(mat, mat, cv::Size(width_image_BP, height_image_BP));
}

memset(data_image, 0, sizeof(int) * (width_image_BP * height_image_BP));

for (int h = 0; h < mat.rows; h++) {
uchar* p = mat.ptr(h);
for (int w = 0; w < mat.cols; w++) {
if (p[w] > 128) {
data_image[h* mat.cols + w] = 1;
}
}
}

int ret = bp2.predict(data_image, mat.cols, mat.rows);
std::cout << "correct result: " << i << ",    actual result: " << ret << std::endl;
}

delete[] data_image;

return 0;
}

train结果如下图所示：



predict结果如下图所示，测试图像是从MNIST test集合中选取的：



GitHub：https://github.com/fengbingchun/NN