目录
- 手写BP
- 前言
- 1. 数据加载
- 2. 前向传播
- 3. 反向传播
- 总结
手写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, float
labels.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-hot
return 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;
}
// 通过指定索引数组，并且指定图像、索引的起点和终点。获取对应索引的对应图像，返回matrix
Matrix 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{
// 对应的就是seed
static default_random_engine global_random_engine;
Matrix normal_distribution_matrix(int rows, int cols, float mean = 0.0f, float stddev = 0.0f){
// seed
normal_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 10
auto 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(); // 60000
int num_input = trainimages.cols(); // 784
int 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 dataloader
int 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 m
auto 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. 反向传播
我们先来回忆下反向传播的公式
根据上图来计算 L o s s Loss Loss 分别对 W 1 W_1 W1 B 1 B_1 B1 W 2 W_2 W2 B 2 B_2 B2 的梯度
-
1.隐藏层的输出为： H = r e l u ( X W 1 + B 1 ) H = relu(XW_1 + B_1) H=relu(XW1+B1)
-
2.输出层的预测概率为： P = s i g m o i d ( H W 2 + B 2 ) P = sigmoid(HW_2 + B_2) P=sigmoid(HW2+B2)
-
3.损失为： L = B i n a r y C r o s s E n t r o p y L o s s ( P , Y ) L = BinaryCrossEntropyLoss(P,Y) L=BinaryCrossEntropyLoss(P,Y)
-
4.计算 L L L 对 W 2 W_2 W2 和 B 2 B_2 B2 的梯度： ∂ L ∂ W 2 = H T ( P − Y ) \frac{\partial L}{\partial W_2}=H^{T}(P-Y) ∂W2∂L=HT(P−Y) ∂ L ∂ B 2 = r e d u c e _ s u m ( P − Y ) \frac{\partial L}{\partial B_{2}}= reduce\_sum(P-Y) ∂B2∂L=reduce_sum(P−Y)
-
5.计算 L L L 对 W 1 W_1 W1 和 B 1 B_1 B1 的梯度： ∂ L ∂ W 1 = X T ∂ L ∂ ( X W 1 + B 1 ) \frac{\partial L}{\partial W_{1}}=X^{T}\frac{\partial L}{\partial(X W_{1}+B_{1})} ∂W1∂L=XT∂(XW1+B1)∂L ∂ L ∂ B 1 = r e d u c e _ s u m ∂ L ∂ ( X W 1 + B 1 ) \frac{\partial L}{\partial B_{1}}=reduce\_sum\frac{\partial L}{\partial(X W_{1}+B_{1})} ∂B1∂L=reduce_sum∂(XW1+B1)∂L
-
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 10
auto 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){
// 为了拿到每一行预测的标签值，使用argmax
int 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(); // 60000
int num_input = trainimages.cols(); // 784
int 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 dataloader
int 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 m
auto 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 @ w2
auto 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 / dhidden
auto dhidden = nn::drelu(dhidden_act, hidden);
// x @ w1 + b1
auto 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 L L 分别对 W 1 W_1 W1 B 1 B_1 B1 W 2 W_2 W2 B 2 B_2 B2 的导数计算，后续就是利用 SGD 算法进行参数更新。最后在每个 epoch 结束后，计算测试集上的准确率并输出结果
总结
本次课程跟着杜老师手写了一遍 BP 反向传播算法，加深了对 BP 的认识，锻炼了自己的动手能力，提高了对 C++ 和神经网络的进一步理解。整个 BP 过程还是比较清晰的，无非是前传 => 计算Loss => 反传 => 更新，但还是要求对许多细节的把握，比如如何去优化代码提高性能，如何尽可能的减少矩阵相乘的次数等等，多想多写吧。期待下次手写 AutoGrad 课程😄
