手写BP

前言

手写AI推出的全新面向AI算法的C++课程 Algo C++，链接。记录下个人学习笔记，仅供自己参考。

本次课程主要是手写 BP 代码

课程大纲可看下面的思维导图

1. 数据加载

我们首先来实现下MNIST数据集的加载工作

#include <iostream>
#include <tuple>
#include <string>
#include <fstream>
#include "matrix.hpp"using namespace std;namespace Application{namespace io{// 不要内存对齐struct __attribute__((packed)) mnist_labels_header_t{unsigned int magic_number;unsigned int number_of_items;}struct __attribute__((packed)) mnist_images_header_t{unsigned int magic_number;unsigned int number_of_images;unsigned int number_of_rows;unsigned int number_of_cols;}unsigned int inverse_byte(unsigned int v){unsigned char* p = (unsigned char*)&v;std::swap(p[0], p[3]);std::swap(p[1], p[2]);return v;}/* 加载mnist数据集 */tuple<Matrix, Matrix> load_data(const string& image_file, const string& label_file){// file, mode = rb/wb/a+Matrix images, labels;fstream fimage(image_file, ios::binary | ios::in);fstream flabel(label_file, ios::binary | ios::in)mnist_images_header_t images_header;mnist_labels_header_t labels_header;fimage.read((char *)&images_header, sizeof(images_header));flabel.read((char *)&labels_header, sizeof(labels_header));// 大小端转换images_header.number_of_images = inverse_byte(images_header.number_of_images);labels_header.number_of_items  = inverse_byte(labels_header.number_of_items);images.resize(images_header.number_of_images, 28 * 28);// one hot, floatlabels.resize(labels_header.number_of_items, 10);// 中间存储作用, 它的大小等于文件中图像数据的大小std::vector<unsigned char> buffer(images.rows() * images.cols());fimage.read((char*)buffer.data(), buffer.size());// buffer 是unsigned char类型的图像数据，pixel// 现在需要转化到images矩阵中// 顺便把标准化给做了for(int i = 0; i < buffer.size() ++i)images.ptr()[i] = (buffer[i] / 255.0f - 0.1307f) / 0.3081f;// 开始处理label，ont-hot过程// labels是全零的矩阵buffer.resize(labels.rows());flabel.read((char*)buffer.data(), buffer.size());for(int i = 0; i < buffer.size(); ++i)labels.ptr(i)[buffer[i]] = 1;	// one-hotreturn make_tuple(images, labels);}}int run(){// 验证Matrix images, labels;tie(images, labels) = io::load_data("mnist/train-images.idx3-ubyte","mnist/train-labels.idx1-ubyte");cout << labels;return 0;}
}int main(){return Application::run();
}

上面示例代码演示了 mnist 数据集加载的示例，其中

• 定义了两个用于 mnist 数据集的结构体，分别是 mnist_labels_header_t 和 mnist_images_header_t。这两个结构体表示 mnist 标签和图像的头信息
• __attribute__((packed)) 是一个 GCC/Clang 的扩展，用于告诉编译器不要进行内存对齐，以便我们自己手动控制内存的布局和对齐。
• inverse_byte 用于大小端转换。在 mnist 数据集中，头文件中存储的是大端模式，而读取数据时需要将其转换为小端模式。inverse_byte 函数实现的方法是将 unsigned int 类型的变量按照字节顺序重新排列。
• load_data 函数用于加载 mnist 数据集，需要传入图像和标签的文件名。在函数内部，使用 fstream 打开文件并读取头信息，对头信息进行大小端转换，然后根据图像和标签的大小创建矩阵。之后，读取图像和标签数据并进行归一化，最后返回归一化后的矩阵。

标签处理的代码有点看不懂，来分析下

buffer.resize(labels.rows());
flabel.read((char*)buffer.data(), buffer.size());
for(int i = 0; i < buffer.size(); ++i)labels.ptr(i)[buffer[i]] = 1;

在 MNIST 数据集中，标签是一个整数值，代表着手写数字的真实值。但是，在神经网络训练中，我们通常需要将标签转化为 one-hot 编码，以便于计算误差。例如，对于数字 5，它的 ont-hot 编码为[0,0,0,0,0,1,0,0,0,0]，也就是在第 6 个位置上为 1，其余位置都为0

这段代码中，首先通过 resize 函数将 buffer 的大小设置为标签行数，也就是样本数。然后从 flabel 中读取 buffer.size() 个字节到 buffer 中。由于标签是一个字节大小的整数，所以 buffer 中存储的是一连串的整数值。接下来的 for 循环中，遍历 buffer 中的所有元素，将标签对应的位置设置为 1，其余位置设置为 0

labels.ptr(i) 返回 labels 中第i行的地址，然后用 [buffer[i]] 访问该行中第 buffer[i] 列的元素，将该元素的值设置为1。这样就完成了标签数据的转化，使其变为了适合神经网络训练的形式

2. 前向传播

来看下整个前向传播过程，包括初始化超参数和权重、偏置矩阵，矩阵相乘进行前传并求取 Loss

#include <iostream>
#include <tuple>
#include <string>
#include <fstream>
#include <random>
#include <cmath>
#include <string.h>
#include <algorithm>    // shuffle
#include "matrix.hpp"using namespace std;namespace Application{namespace tools{vector<int> range(int end){vector<int> output(end);for(int i = 0; i < end; ++i)output[i] = i;return  output;}// 通过指定索引数组，并且指定图像、索引的起点和终点。获取对应索引的对应图像，返回matrixMatrix slice(const Matrix& images, const vector<int>& indexs, int start, int size){Matrix output(size, images.cols());for(int i = 0; i < size; ++i){int image_row_index = indexs[start + i];// 把images[image_row_index]对应的一行，复制到output[i]行memcpy(output.ptr(i), images.ptr(image_row_index), sizeof(float) * images.cols());}return output;}};namespace random{// 对应的就是seedstatic default_random_engine global_random_engine;Matrix normal_distribution_matrix(int rows, int cols, float mean = 0.0f, float stddev = 0.0f){// seednormal_distribution<float> norm(mean, stddev);Matrix output(rows, cols);for(int i = 0; i < rows * cols; ++i)output.ptr()[i] = norm(global_random_engine);return output;}// inplace -> 现场修改void shuffle_array(vector<int>& indexs){std::shuffle(indexs.begin(), indexs.end(), global_random_engine);}};namespace nn{Matrix relu(cosnt Matrix& x){Matrix output = x;auto p = output.ptr();for(int i = 0; i < x.numel(); ++i, ++p){// max(float, int)*p = std::max(*p, 0.0f);}return output;            }float sigmoid_cross_entropy_loss(const Matrix& output, const Matrix& y){static auto sigmoid = [](float x){return 1.0f / (1.0f + exp(-x));}// output => n x 10// y      => n x 10auto po = output.ptr();auto py = y.ptr();float loss = 0;float eps = 1e-5;for(int i = 0; i < output.numnel(); ++i, ++po, ++py){float prob = sigmoid(*po);prob = max(min(prob, 1-eps), eps);// loss += -(y * log(p) + (1-y) * log(1-p));loss += (*py) * log(prob) + (1 - *py) * log(1 - prob);}return -loss / output.rows();}};int run(){Matrix trainimages, trainlabels;Matrix valimages, vallabels;tie(trainimages, trainlabels) = io::load_data("mnist/train-images.idx3-ubyte", "mnist/train-labels.idx1-ubyte");tie(valimage, vallabels)      = io::load_data("mnist/t10k-images.idx3-ubyte",  "mnist/t10k-labels.idx1-ubyte");// hidden = images @ w1 + b1// output = hidden @ w2 + b2// loss   = lossfn(output, onehot_label)// loss.backward()int num_images = trainimages.rows();  // 60000int num_input  = trainimages.cols();  // 784int num_hidden = 1024;int num_output = 10;int num_epoch  = 10;float lr       = 1e-1;int batch_size = 256;float momentum = 0.9f;// drop last -> pytorch dataloaderint num_batch_per_epoch = num_images / batch_size;// 训练流程// 1. 随机抓取一个batch，shuffle=True//      怎么实现随机呢？//      重点是，要求，每一个epoch，抓取的随机性不同//     按照索引抓取，index[0,1,2,3,4,5,6,7]//     shuffle(indexs) ->[6,5,4,0,1,3,2,7] auto image_indexs = tools::range(num_images);// 确定矩阵的大小，并且对矩阵做初始化(凯明初始化fan_in,fan_out)float w1_gain =  2.0f / std::sqrt((float)num_input + num_hidden); // 对应激活函数        Matrix w1 = random::normal_distribution_matrix(num_input, num_hidden, 0, w1_gain);Matrix b1 = random::normal_distribution_matrix(1, num_hidden);        float w2_gain =  1.0f / std::sqrt((float)num_hidden + num_output); // 对应激活函数        Matrix w2 = random::normal_distribution_matrix(num_hidden, num_output, 0, w1_gain);Matrix b2 = random::normal_distribution_matrix(1, num_output);for(int epoch = 0; epoch < num_epoch; ++epoch){random::shuffle_array(image_indexs);for(int ibatch = 0; ibatch < num_batch_per_epoch; ++ibatch){// image_indexs[ibatch * batch_size : (ibatch + 1) * batch_size]auto x = tools::slice(trainimages, image_indexs, ibatch * batch_size, batch_size);cout << x.rows() << ", " << x.cols() << endl;auto y = tools::slice(trainlabels, image_indexs, ibatch * batch_size, batch_size);// n x m + 1 x mauto hidden = x.gemm(w1) + b1;auto hidden_act = nn::relu(hidden);auto output     = hidden_act.gemm(w2) + b2;float loss      = nn::sigmoid_cross_entropy_loss(output, y);cout << "Loss: " << loss << endl;}}// 验证// Matrix t = random::normal_distribution_matrix(3, 3);// Matrix s = tools::slice(t, {0, 1, 2}, 1, 2);// cout << t;// cout << s;return 0;}
}int main(){return Application::run();
}

上面示例代码演示了前向传播的过程，包括初始化神经网络参数、进行前向传播计算、计算损失函数

• 首先，对神经网络的参数进行了初始化，包括输入层到隐藏层的权重矩阵 w1 和偏置向量 b1，隐藏层到输出层的权重矩阵 w2 和偏置向量 b2。这些参数都是通过调用 random::normal_distribution_matrix() 函数生成的，该函数使用了高斯分布来初始化，而 w1 和 w2 的初始化还使用了凯明初始化
• 然后，代码进入了训练流程，首先进行了数据的随机抽取。随机抽取的方法是使用了 tools::range() 函数生成一个索引数组，然后使用 random::shuffle_array() 函数对该索引数组进行打乱。打乱的目的是保证每个 epoch 中数据的随机性不同。
• 对于每个 batch，代码首先使用 tools::slice() 函数从训练集抽取出 batch_size 个样本，然后进行前向传播计算。
  • hidden = images @ w1 + b1
  • output = hidden @ w2 + b2
• 在得到输出层的数据后，代码计算了该 batch 的损失函数。其具体实现在 nn::sigmoid_cross_entropy_loss() 函数中，采用交叉熵损失，它首先将 output 中的每个元素通过 sigmoid 函数，然后根据公式计算 loss 并返回。

3. 反向传播

我们先来回忆下反向传播的公式

根据上图来计算 $Loss$ 分别对 $W_1$ $B_1$ $W_2$ $B_2$ 的梯度

• 1.隐藏层的输出为：$H=relu(X{W}_1+B_1)$

• 2.输出层的预测概率为：$P=sigmoid(H{W}_2+B_2)$

• 3.损失为：$L=BinaryCrossEntropyLoss(P,Y)$

• 4.计算 $L$ 对 $W_2$ 和 $B_2$ 的梯度：$\frac{\partial L}{\partial W_2}=H^T(P-Y)$ $\frac{\partial L}{\partial B_2}=reduce\_sum(P-Y)$

• 5.计算 $L$ 对 $W_1$ 和 $B_1$ 的梯度：$\frac{\partial L}{\partial W_1}=X^T\frac{\partial L}{\partial (X{W}_1+B_1)}$ $\frac{\partial L}{\partial B_1}=reduce\_sum\frac{\partial L}{\partial (X{W}_1+B_1)}$

• 6.拿到梯度后，对每一个参数应用优化器进行更新迭代

来看下整个反向传播的过程，根据 Loss 分别对 w1、b1、w2、b2 求取梯度，并更新参数

#include <iostream>
#include <tuple>
#include <string>
#include <fstream>
#include <random>
#include <cmath>
#include <string.h>
#include <algorithm>    // shuffle
#include "matrix.hpp"using namespace std;namespace Application{namespace nn{Matrix drelu(const Matrix& g, cosnt Matrix& x){Matrix output = g;auto px = x.ptr();auto pg = output.ptr();for(int i = 0; i < x.numel(); ++i, ++px, ++pg){if(*px < 0)*pg = 0;}return output;}Matrix relu(cosnt Matrix& x){Matrix output = x;auto p = output.ptr();for(int i = 0; i < x.numel(); ++i, ++p){// max(float, int)*p = std::max(*p, 0.0f);}return output;            }Matrix sigmoid(const Matrix& x){Matrix output = x;auto p = output.ptr();for(int i = 0; i < x.numel(); ++i, ++p){float t = *p;*p = 1.0f / (1.0f + exp(-t));}return output;}float sigmoid_cross_entropy_loss(const Matrix& prob, const Matrix& y){// static auto sigmoid = [](float x){return 1.0f / (1.0f + exp(-x));}auto prob = sigmoid(output);// output => n x 10// y      => n x 10auto po = prob.ptr();auto py = y.ptr();float loss = 0;float eps = 1e-5;for(int i = 0; i < prob.numel(); ++i, ++po, ++py){float prob = *po;// prob = max(min(prob, 1-eps), eps);// loss += -(y * log(p) + (1-y) * log(1-p));loss += (*py) * log(p) + (1 - *py) * log(1 - p);}return -loss / prob.rows();}float eval_accuracy(const Matrix& testset, const Matrix& y,const Matrix& w1, const Matrix& b1, const Matrix& w2, const Matrix& b2){auto hidden = nn::relu(testset.gemm(w1) + b1);auto prob   = nn::sigmoid(hidden.gemm(w2) + b2);int correct = 0;for(int i = 0; i < prob.rows(); ++i){// 为了拿到每一行预测的标签值，使用argmaxint predict_label = prob.argmax(i);// 因为y是one-hot，所以对应预测标签如果是1，则正确if(y(i, predict_label) == 1)correct++;}return correct / (float)prob.rows();}};int run(){Matrix trainimages, trainlabels;Matrix valimages, vallabels;tie(trainimages, trainlabels) = io::load_data("mnist/train-images.idx3-ubyte", "mnist/train-labels.idx1-ubyte");tie(valimages, vallabels)      = io::load_data("mnist/t10k-images.idx3-ubyte",  "mnist/t10k-labels.idx1-ubyte");// hidden = images @ w1 + b1// output = hidden @ w2 + b2// loss   = lossfn(output, onehot_label)// loss.backward()int num_images = trainimages.rows();  // 60000int num_input  = trainimages.cols();  // 784int num_hidden = 1024;int num_output = 10;int num_epoch  = 10;float lr       = 1e-1;int batch_size = 256;float momentum = 0.9f;// drop last -> pytorch dataloaderint num_batch_per_epoch = num_images / batch_size;  // 60000 / 256 = 234(iter)// 训练流程// 1. 随机抓取一个batch，shuffle=True//      怎么实现随机呢？//      重点是，要求，每一个epoch，抓取的随机性不同//     按照索引抓取，index[0,1,2,3,4,5,6,7]//     shuffle(indexs) ->[6,5,4,0,1,3,2,7] auto image_indexs = tools::range(num_images);// 确定矩阵的大小，并且对矩阵做初始化(凯明初始化fan_in,fan_out)float w1_gain =  2.0f / std::sqrt((float)num_input + num_hidden); // 对应激活函数        Matrix w1 = random::normal_distribution_matrix(num_input, num_hidden, 0, w1_gain);Matrix b1 = random::normal_distribution_matrix(1, num_hidden);        float w2_gain =  1.0f / std::sqrt((float)num_hidden + num_output); // 对应激活函数        Matrix w2 = random::normal_distribution_matrix(num_hidden, num_output, 0, w1_gain);Matrix b2 = random::normal_distribution_matrix(1, num_output);for(int epoch = 0; epoch < num_epoch; ++epoch){random::shuffle_array(image_indexs);for(int ibatch = 0; ibatch < num_batch_per_epoch; ++ibatch){    // 一个epoch// image_indexs[ibatch * batch_size : (ibatch + 1) * batch_size]auto x = tools::slice(trainimages, image_indexs, ibatch * batch_size, batch_size);cout << x.rows() << ", " << x.cols() << endl;auto y = tools::slice(trainlabels, image_indexs, ibatch * batch_size, batch_size);// n x m + 1 x mauto hidden      = x.gemm(w1) + b1;auto hidden_act  = nn::relu(hidden);auto output      = hidden_act.gemm(w2) + b2;auto probability = nn::sigmoid(output)float loss       = nn::sigmoid_cross_entropy_loss(probability, y);cout << "Loss: " << loss << endl;// backward过程// dloss / doutput = (sigmoid(output) - y) / output.rows// C = AB// G = dloss / dC// dloss / dA = G @ B^T// dloss / dB = A^T @ G// output = hidden_act @ w2auto doutput     = (probability - y) / (float)output.rows();auto db2         = doutput.reduce_sum_by_row();auto dhidden_act = doutput.gemm(w2, false, true);auto dw2         = hidden_act.gemm(doutput, true);// dloss / dhiddenauto dhidden     = nn::drelu(dhidden_act, hidden);// x @ w1 + b1auto db1         = dhidden.reduce_sum_by_row();auto dw1         = x.gemm(dhidden, true);// 更新w1 = w1 - lr * dw1;b1 = b1 - lr * dw1;}float accuracy = nn::eval_accuracy(valimages, vallabels, w1, b1, w2, b2);printf("%d. accuracy: %.2f %%\n", epoch, accuracy * 100);}// 验证// Matrix t = random::normal_distribution_matrix(3, 3);// Matrix s = tools::slice(t, {0, 1, 2}, 1, 2);// cout << t;// cout << s;return 0;}
}int main(){return Application::run();
}

在上面的示例代码中，主要是求取梯度的过程比较繁琐，需要明确 $L$ 分别对 $W_1$ $B_1$ $W_2$ $B_2$ 的导数计算，后续就是利用 SGD 算法进行参数更新。最后在每个 epoch 结束后，计算测试集上的准确率并输出结果

总结

本次课程跟着杜老师手写了一遍 BP 反向传播算法，加深了对 BP 的认识，锻炼了自己的动手能力，提高了对 C++ 和神经网络的进一步理解。整个 BP 过程还是比较清晰的，无非是前传 => 计算Loss => 反传 => 更新，但还是要求对许多细节的把握，比如如何去优化代码提高性能，如何尽可能的减少矩阵相乘的次数等等，多想多写吧。期待下次手写 AutoGrad 课程😄