BP反向传播算法

前言

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

本次课程主要是讲解BP反向传播算法

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

在这里插入图片描述

1. MNIST

MNIST是一个手写数字识别数据集，它包含60,000个训练图像和10,000个测试图像。每个图像是28x28像素大小的灰度图像，表示0到9之间的一个数字。该数据集最初由美国国家标准与技术研究所（NIST）在20世纪80年代收集和创建，后来由MNIST数据库维护。该数据集已经成为机器学习和计算机视觉领域的标准基准数据集之一。(from chatGPT)

MNIST官网有关于手写数字数据集及其格式的详细介绍

MNIST数据集的格式为二进制格式，分为训练集和测试集两部分，每个部分都由一个图像文件和一个标签文件组成。训练集包含60,000个图像，测试集包含10,000个图像。图像文件和标签文件的名称分别为：

训练集图像文件：train-images-idx3-ubyte.gz
训练集标签文件：train-labels-idx1-ubyte.gz
测试集图像文件：t10k-images-idx3-ubyte.gz
测试集标签文件：t10k-labels-idx1-ubyte.gz

图像文件包含原始的灰度图像数据，标签文件包含每个图像对应的数字标签。

具体来说，图像文件的格式如下：

前4个字节：magic number(MSB first)
4个字节：数据集中图像的数量
4个字节：每个图像的行数
4个字节：每个图像的列数
后续的字节：表示每个图像的像素值，按照行优先的顺序排列

标签文件的格式如下：

前4个字节：magic number(MSB first)
4个字节：数据集中标签的数量
后续的字节：表示每个图像对应的标签，取值范围为0到9

因此，MNIST数据集中的每个图像都是一个28x28像素的灰度图像，由0到255之间的整数值表示每个像素的灰度级别。标签是一个0到9之间的整数值，表示每个图像对应的数字。

下面的示例代码用于加载MNIST数据集并返回对应的图像和标签文件：

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <vector>

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_columns;
};

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;
}

int main(){

    FILE* f = fopen("mnist/train-labels.idx1-ubyte", "rb");
    mnist_labels_header_t labels_header{0};
    fread(&labels_header, 1, sizeof(labels_header), f);
    printf("labels_header.magic_number = %X, number_of_items = %d\n", 
        inverse_byte(labels_header.magic_number), inverse_byte(labels_header.number_of_items));
    
    unsigned char label = 0;
    fread(&label, 1, sizeof(label), f);
    printf("First label is: %d\n", label);
    fclose(f);

    f = fopen("mnist/train-images.idx3-ubyte", "rb");
    mnist_images_header_t images_header{0};
    fread(&images_header, 1, sizeof(images_header), f);
    printf("images_header.magic_number = %X, number_of_images = %d, number_of_rows = %d, number_of_columns = %d\n", 
        inverse_byte(images_header.magic_number), 
        inverse_byte(images_header.number_of_images),
        inverse_byte(images_header.number_of_rows),
        inverse_byte(images_header.number_of_columns)
    );
    
    std::vector<unsigned char> image(inverse_byte(images_header.number_of_rows) * inverse_byte(images_header.number_of_columns));
    fread(image.data(), 1, image.size(), f);
    for(int i = 0;i < image.size(); ++i){
        if(image[i] == 0)
            printf("--- ");
        else
            printf("%03d ", image[i]);

        if((i + 1) % inverse_byte(images_header.number_of_columns) == 0)
            printf("\n");
    }
    fclose(f);
    return 0;
}

上述示例代码用于解析 MNIST 数据集中的标签和图像数据。代码首先定义了两个结构体：mnist_labels_header_t 和 mnist_images_header_t，分别表示标签和图像的头部信息，其中使用 packed 属性来确保结构体不会被编译器优化对齐。

接下来定义了一个 inverse_byte 函数，用于将4字节整数类型的大小端转换

主函数首先打开标签文件，读取标签头部信息并打印。然后读取第一个标签，并打印；接着打开图像文件，读取图像头部信息并打印；然后读取第一张图像，将其打印出来；最后关闭文件。

运行效果如下图所示：

在这里插入图片描述

2. 感知机

2.1 前言

感知机可以简单用下图表示：

在这里插入图片描述

上图表示了：

第一点： $a\times W_1 + b\times W_2$
第二点： $d = a c t i v a t i o n (c)$ 这里的 activation 是一个激活函数
激活函数通常是为了非线性映射，例如 $\frac{1}{1+e^{-x}}$ 或者 $re l u = ma x (0, x)$

重点：

任何两个节点的连接线是具有权重值的，例如 $W_1$ ， $W_2$
多个节点连接到一个节点，指这多个节点值加权求和后，经过激活函数的结果即： $\times W_1 + b \times W_2)$

2.2 感知机-矩阵表示

我们在理解感知机可以通过用矩阵的方式来理解，定义 $A=\left\{\begin{array}{cc}a&b\end{array}\right\}$ $W=\left\{\begin{array}{c}w1\\ w2\end{array}\right\}$ 则输出为
$c=\left\{a\quad b\right\}\times\left\{\begin{array}{c}w1\\ w2\end{array}\right\}=AW$

2.3 感知机-矩阵表示-多个样本

我们来看下多个样本的情况，当增加了一个样本后其实是在矩阵 $A$ 中新增了一行(注意权重是同一组只是样本换了)

在这里插入图片描述

定义 $A=\left\{\begin{array}{cc}a&b\end{array}\right\}$ $W=\left\{\begin{array}{c}w1\\ w2\end{array}\right\}$ 那么输出为
$\left\{\begin{array}{c}c\\ f\end{array}\right\}=\left\{\begin{array}{c}a&b\\ x&y\end{array}\right\}\times\left\{\begin{array}{c}w1\\ w2\end{array}\right\}=AW$

2.4 感知机-增加偏置

关于偏置的存在，考虑 $y = k x + b$ 直线公式，若 $b = 0$ ，则退化为 $y = k x$ ，此时表达的直接必定过零点，无法表达不过零点的直线，所以偏置在这里非常重要，感知机增加偏置后的图如下所示：

在这里插入图片描述

$c=\left\{\begin{array}{c}a&b\end{array}\right\}\times\left\{\begin{array}{c}w1\\ w2\end{array}\right\}+bias=AW+bias$

2.5 感知机-多个输出

当感知机存在多个输出时，如下图所示，

在这里插入图片描述

定义 $c=\left\{a\quad b\right\}\times\left\{\begin{array}{c}w11\\ w12\end{array}\right\}+bias1$ ， $q=\left\{\begin{array}{c}a&b\end{array}\right\}\times\left\{\begin{array}{c}w21\\ w22\end{array}\right\}+bias2$

将输出 $c$ 和 $q$ 合并为一个矩阵，输出如下：
$\begin{aligned} \left\{\begin{array}{cc}c&q\end{array}\right\}=\left\{\begin{array}{cc}a&b\end{array}\right\}\times\left\{\begin{array}{cc}w11&w21\\ w12&w22\end{array}\right\}+\left\{\begin{array}{cc}bias1&bias2\end{array}\right\}=AW+B \end{aligned}$

2.6 总结

新增一个样本， $A$ 增加一行
新增一个输出， $W$ 增加一列
$A$ 的行数是样本数， $A$ 的列数是特征数
$W$ 的行数是输入特征数， $W$ 的列数是输出特征数
可以认为 $A$ 经过 $W$ 映射为新的特征

2.7 关于广播

广播机制是非常重要的一种特性，它可以使得不同形状的矩阵在一些条件下能够进行数学运算。当运算中两个矩阵地形状不一致时，会自动扩展，以满足运算条件，这个过程就称为广播。

对于矩阵 $A$ 和 $B$ 的元素操作(如点乘、点加、点除等等)，广播约定了假设 $A$ 是 $1\times 5$ ， $B$ 是 $3\times 5$ ，则约定把 $A$ 在行方向复制 3 份后，再与 $B$ 进行元素操作。同理可以发生在列上，发生在 $B$ 上。

对于
$\left\{\begin{array}{cc}c&q\\f&m\end{array}\right\}=\left\{\begin{array}{cc}a&b\\x&y\end{array}\right\}\times\left\{\begin{array}{cc}w11&w21\\w12&w22\end{array}\right\}+\left\{\begin{array}{cc}bias1&bias2\end{array}\right\}=AW+B$
等价于
$\left\{\begin{array}{cc}c&q\\f&m\end{array}\right\}=\left\{\begin{array}{cc}a&b\\x&y\end{array}\right\}\times\left\{\begin{array}{cc}w11&w21\\w12&w22\end{array}\right\}+\left\{\begin{array}{cc}bias1&bias2\\bias1&bias2\end{array}\right\}=AW+B$

3. BP

BP(Back Propagation)误差反向传播算法，使用反向传播算法的多层感知器又称为 BP 神经网络。BP 是当前人工智能主要采用的算法，例如所知道的 CNN、GAN、NLP的Bert、Transformer，都是 BP 体系下的算法框架。

理解 BP 对于理解网络如何训练很重要

在这里我们采用最简单的思路理解BP。确保能够理解并且复现

在这里插入图片描述

使用BP的训练流程：

1.计算隐藏层的输出： $H = relu(XW_1+B_1)$
2.计算输出层的预测概率： $P = sigmoid(HW_2+B_2)$
3.计算损失： $L = B ina ry C ross E n t ro p y L oss (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.拿到梯度后，对每一个参数应用优化器进行更新迭代

部分核心代码如下：

// 开始循环所有的batch
for(int ibatch = 0; ibatch < num_batch_per_epoch; ++ibatch, ++total_batch){
    
    // 前向传播
    auto x           = data::choice_rows(trainimages, image_indexs, ibatch * batch_size, batch_size);
    auto y           = data::choice_rows(trainlabels, image_indexs, ibatch * batch_size, batch_size);
    auto hidden      = x.gemm(input_to_hidden) + hidden_bias;
    auto hidden_act  = nn::relu(hidden);
    auto output      = hidden_act.gemm(hidden_to_output) + output_bias;
    auto probability = nn::sigmoid(output);
    float loss 		 = nn::compute_loss(probability, y);
    
    // 反向传播
    // C = AB
    // dA =  G * BT
    // dB = AT * G
    // loss部分求导，loss对output求导
    auto doutput           = (probability - y) * (1 / (float)batch_size);
    
    // 第二个Linear求导
    auto doutput_bias      = data::reduce_sum_by_row(output);
    auto dhidden_to_output = hidden_act.gemm(doutput, true);
    auto hidden_act        = doutput.gemm(hidden_to_output, false, true);
    
    // 第一个Linear输出求导
    auto dhidden           = nn::drelu(dhidden_act, hidden);
    
    // 第一个Linear求导
    auto dinput_to_hidden  = x.gemm(dhidden, true);
    auto dhidden_bias      = data::reduce_sum_by_row(dhidden);
    
    // 调用优化器来调整更新参数
    optim.update_params(
    	(&input_to_hidden,  &hidden_bias,  &hidden_to_output,  &output_bias),
        (&dinput_to_hidden, &dhidden_bias, &dhidden_to_output, $doutput_bias),
    	lr, momentum
    );
}

4. 动量SGD

对于参数更新方向等于 $- l r * g r a d$ ，我们定义 $D = - l r * g r a d$

而梯度下降时，我们有： $\theta^+=\theta+D$

假设梯度方向固定沿着右边取值相同，则每个时刻的推进都是均匀的，如下图所示：

在这里插入图片描述

对于动量 Momentum，则是基于物理上的惯性设计，定义动量系数 $m$

定义 $t_1$ 时刻的累计梯度量： $D_{t0}\cdot m+ D_{t1}$ 其中 $D_{t0} = 0$

$A$ 就是动量 SGD 的参数更新方向 $\theta^+=\theta+A$

假设梯度方向固定沿着右边取值相同，则每个时刻的推进都有惯性作用，也可以连续下降的区域，会具有更快的下降速度。若在梯度方向不同时，也会存在正负抵消，从而更小心翼翼的前进，如下图所示：

在这里插入图片描述

5. BP示例代码

main.cpp

#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <cmath>
#include <tuple>
#include <iomanip>
#include <stdarg.h>
#include <memory.h>
#include <random>
#include <algorithm>
#include <chrono>
#include "matrix.hpp"

using namespace std;

namespace Application{

    static default_random_engine global_random_engine;

    namespace logger{

        #define INFO(...)  Application::logger::__printf(__FILE__, __LINE__, __VA_ARGS__)

        void __printf(const char* file, int line, const char* fmt, ...){

            va_list vl;
            va_start(vl, fmt);

            // None   = 0,     // 无颜色配置
            // Black  = 30,    // 黑色
            // Red    = 31,    // 红色
            // Green  = 32,    // 绿色
            // Yellow = 33,    // 黄色
            // Blue   = 34,    // 蓝色
            // Rosein = 35,    // 品红
            // Cyan   = 36,    // 青色
            // White  = 37     // 白色
            /* 格式是： \e[颜色号m文字\e[0m   */
            printf("\e[32m[%s:%d]:\e[0m ", file, line);
            vprintf(fmt, vl);
            printf("\n");
        }
    };

    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_columns;
        };

        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){

            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);
            labels.resize(labels_header.number_of_items, 10);

            std::vector<unsigned char> buffer(images.rows() * images.cols());
            fimage.read((char*)buffer.data(), buffer.size());

            for(int i = 0; i < buffer.size(); ++i)
                images.ptr()[i] = (buffer[i] / 255.0f - 0.1307f) / 0.3081f;
                //images.ptr()[i] = (buffer[i] - 127.5f) / 127.5f;

            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;   // onehot
            return make_tuple(images, labels);
        }

        void print_image(const float* ptr, int rows, int cols){

            for(int i = 0;i < rows * cols; ++i){

                //int pixel = ptr[i] * 127.5 + 127.5;
                int pixel = (ptr[i] * 0.3081f + 0.1307f) * 255.0f;
                if(pixel < 20)
                    printf("--- ");
                else
                    printf("%03d ", pixel);

                if((i + 1) % cols == 0)
                    printf("\n");
            }
        }

        bool save_model(const string& file, const vector<Matrix>& model){

            ofstream out(file, ios::binary | ios::out);
            if(!out.is_open()){
                INFO("Open %s failed.", file.c_str());
                return false;
            }

            unsigned int header_file[] = {0x3355FF11, model.size()};
            out.write((char*)header_file, sizeof(header_file));

            for(auto& m : model){
                int header[] = {m.rows(), m.cols()};
                out.write((char*)header, sizeof(header));
                out.write((char*)m.ptr(), m.numel() * sizeof(float));
            }
            return out.good();
        }

        bool load_model(const string& file, vector<Matrix>& model){

            ifstream in(file, ios::binary | ios::in);
            if(!in.is_open()){
                INFO("Open %s failed.", file.c_str());
                return false;
            }

            unsigned int header_file[2];
            in.read((char*)header_file, sizeof(header_file));

            if(header_file[0] != 0x3355FF11){
                INFO("Invalid model file: %s", file.c_str());
                return false;
            }

            model.resize(header_file[1]);
            for(int i = 0; i < model.size(); ++i){
                auto& m = model[i];
                int header[2];
                in.read((char*)header, sizeof(header));
                m.resize(header[0], header[1]);
                in.read((char*)m.ptr(), m.numel() * sizeof(float));
            }
            return in.good();
        }
    };

    namespace data{

        int argmax(float* ptr, int size){
            return std::max_element(ptr, ptr + size) - ptr;
        }

        Matrix choice_rows(const Matrix& m, const vector<int>& indexs, int begin=0, int size=-1){

            if(size == -1) size = indexs.size();
            Matrix out(size, m.cols());
            for(int i = 0; i < size; ++i){
                int mrow = indexs[i + begin];
                int orow = i;
                memcpy(out.ptr(orow), m.ptr(mrow), sizeof(float) * m.cols());
            }
            return out;
        }

        Matrix reduce_sum_by_row(const Matrix& value){
            Matrix out(1, value.cols());
            auto optr = out.ptr();
            auto vptr = value.ptr();
            for(int i = 0; i < value.numel(); ++i, ++vptr)
                optr[i % value.cols()] += *vptr;
            return out;
        }
    };

    namespace tools{

        vector<int> range(int end){
            vector<int> out(end);
            for(int i = 0; i < end; ++i)
                out[i] = i;
            return out;
        }

        double timenow(){
            return chrono::duration_cast<chrono::microseconds>(chrono::system_clock::now().time_since_epoch()).count() / 1000.0;
        }
    };

    namespace nn{

        Matrix relu(const Matrix& input){
            Matrix out(input);
            for(int i = 0; i < out.numel(); ++i)
                out.ptr()[i] = std::max(0.0f, out.ptr()[i]);
            return out;
        }

        Matrix drelu(const Matrix& grad, const Matrix& x){
            Matrix out = grad;
            auto optr = out.ptr();
            auto xptr = x.ptr();
            for(int i = 0; i < out.numel(); ++i, ++optr, ++xptr){
                if(*xptr <= 0)
                    *optr = 0;
            }
            return out;
        }

        Matrix sigmoid(const Matrix& input){
            Matrix out(input);
            float eps = 1e-5;
            for(int i = 0; i < out.numel(); ++i){
                float& x = out.ptr()[i];

                /* 处理sigmoid数值稳定性问题 */
                if(x < 0){
                    x = exp(x) / (1 + exp(x));
                }else{
                    x = 1 / (1 + exp(-x));
                }

                /* 保证x不会等于0或者等于1 */
                x = std::max(0.0f + eps, std::min(1.0f - eps, x));
            }
            return out;
        }

        float compute_loss(const Matrix& probability, const Matrix& onehot_labels){

            float eps = 1e-5;
            float sum_loss  = 0;
            auto pred_ptr   = probability.ptr();
            auto onehot_ptr = onehot_labels.ptr();
            int numel       = probability.numel();
            for(int i = 0; i < numel; ++i, ++pred_ptr, ++onehot_ptr){
                auto y = *onehot_ptr;
                auto p = *pred_ptr;
                p = max(min(p, 1 - eps), eps);
                sum_loss += -(y * log(p) + (1 - y) * log(1 - p));
            }
            return sum_loss / probability.rows();
        }

        float eval_test_accuracy(const Matrix& probability, const Matrix& labels){

            int success = 0;
            for(int i = 0; i < probability.rows(); ++i){
                auto row_ptr = probability.ptr(i);
                int predict_label = std::max_element(row_ptr, row_ptr + probability.cols()) - row_ptr;
                if(labels.ptr(i)[predict_label] == 1)
                    success++;
            }
            return success / (float)probability.rows();
        }
    };

    namespace random{

        Matrix create_normal_distribution_matrix(int rows, int cols, float mean=0.0f, float stddev=1.0f){

            normal_distribution<float> norm(mean, stddev);
            Matrix out(rows, cols);
            for(int i = 0; i < rows; ++i){
                for(int j = 0; j < cols; ++j)
                    out.ptr(i)[j] = norm(global_random_engine);
            }
            return out;
        }
    };

    namespace optimizer{

        struct SGDMomentum{
            vector<Matrix> delta_momentums;

            // 提供对应的参数params，和对应的梯度grads，进行参数的更新
            void update_params(const vector<Matrix*>& params, const vector<Matrix*>& grads, float lr, float momentum=0.9){

                if(delta_momentums.size() != params.size())
                    delta_momentums.resize(params.size());

                for(int i =0 ; i < params.size(); ++i){
                    auto& delta_momentum = delta_momentums[i];
                    auto& param          = *params[i];
                    auto& grad           = *grads[i];

                    if(delta_momentum.numel() == 0)
                        delta_momentum.resize(param.rows(), param.cols());
                    
                    delta_momentum = momentum * delta_momentum - lr * grad;
                    param          = param + delta_momentum;
                }
            }
        };
    };
    
    int run(){

        Matrix trainimages, trainlabels;
        Matrix valimage, 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");
        
        int num_images  = trainimages.rows();
        int num_input   = trainimages.cols();
        int num_hidden  = 1024;
        int num_output  = 10;
        int num_epoch   = 10;
        float lr        = 1e-1;
        int batch_size  = 256;
        float momentum  = 0.9f;
        int num_batch_per_epoch = num_images / batch_size;
        auto image_indexs       = tools::range(num_images);

        // 凯明初始化，fan_in + fan_out
        // W1 B1
        Matrix input_to_hidden  = random::create_normal_distribution_matrix(num_input,  num_hidden, 0, 2.0f / sqrt((float)(num_input + num_hidden)));
        Matrix hidden_bias(1, num_hidden);

        // W2 B2
        Matrix hidden_to_output = random::create_normal_distribution_matrix(num_hidden, num_output, 0, 1.0f / sqrt((float)(num_hidden + num_output)));
        Matrix output_bias(1, num_output);

        optimizer::SGDMomentum optim;
        auto t0 = tools::timenow();
        int total_batch = 0;
        for(int epoch = 0; epoch < num_epoch; ++epoch){

            if(epoch == 8){
                lr *= 0.1;
            }

            // 打乱索引
            // 0, 1, 2, 3, 4, 5 ... 59999
            // 199, 20, 1, 9, 10, 6, ..., 111
            std::shuffle(image_indexs.begin(), image_indexs.end(), global_random_engine);
            
            // 开始循环所有的batch
            for(int ibatch = 0; ibatch < num_batch_per_epoch; ++ibatch, ++total_batch){

                // 前向过程
                // trainimages -> X(60000, 784)
                // idx = image_indexs[0:256] -> 乱的
                // X = trainimages[idx]
                auto x           = data::choice_rows(trainimages,   image_indexs, ibatch * batch_size, batch_size);
                auto y           = data::choice_rows(trainlabels,   image_indexs, ibatch * batch_size, batch_size);
                auto hidden      = x.gemm(input_to_hidden) + hidden_bias;
                auto hidden_act  = nn::relu(hidden);
                auto output      = hidden_act.gemm(hidden_to_output) + output_bias;
                auto probability = nn::sigmoid(output);
                float loss       = nn::compute_loss(probability, y);

                // 反向过程
                // C = AB
                // dA = G * BT
                // dB = AT * G
                // loss部分求导，loss对output求导
                auto doutput           = (probability - y) * (1 / (float)batch_size);

                // 第二个Linear求导
                auto doutput_bias      = data::reduce_sum_by_row(doutput);
                auto dhidden_to_output = hidden_act.gemm(doutput, true);
                auto dhidden_act       = doutput.gemm(hidden_to_output, false, true);

                // 第一个Linear输出求导
                auto dhidden           = nn::drelu(dhidden_act, hidden);

                // 第一个Linear求导
                auto dinput_to_hidden  = x.gemm(dhidden, true);
                auto dhidden_bias      = data::reduce_sum_by_row(dhidden);

                // 调用优化器来调整更新参数
                optim.update_params(
                    {&input_to_hidden,  &hidden_bias,  &hidden_to_output,  &output_bias},
                    {&dinput_to_hidden, &dhidden_bias, &dhidden_to_output, &doutput_bias},
                    lr, momentum
                );

                if((total_batch + 1) % 50 == 0){
                    auto t1 = tools::timenow();
                    auto batchs_time = t1 - t0;
                    t0 = t1;
                    INFO("Epoch %.2f / %d, Loss: %f, LR: %f [ %.2f ms / 50 batch ]", epoch + ibatch / (float)num_batch_per_epoch, num_epoch, loss, lr, batchs_time);
                }
            }

            // 模型对测试集进行测试，并打印精度
            auto test_hidden      = nn::relu(valimage.gemm(input_to_hidden) + hidden_bias);
            auto test_probability = nn::sigmoid(test_hidden.gemm(hidden_to_output) + output_bias);
            float accuracy        = nn::eval_test_accuracy(test_probability, vallabels);
            float test_loss       = nn::compute_loss(test_probability, vallabels);
            INFO("Test Accuracy: %.2f %%, Loss: %f", accuracy * 100, test_loss);
        }

        INFO("Save to model.bin .");
        io::save_model("model.bin", {input_to_hidden, hidden_bias, hidden_to_output, output_bias});

        for(int i = 0; i < valimage.rows(); ++i){

            auto input = data::choice_rows(valimage, {i});
            auto test_hidden      = nn::relu(input.gemm(input_to_hidden) + hidden_bias);
            auto test_probability = nn::sigmoid(test_hidden.gemm(hidden_to_output) + output_bias);

            int ilabel = data::argmax(test_probability.ptr(), test_probability.cols());
            float prob = test_probability.ptr()[ilabel];

            io::print_image(input.ptr(), 28, 28);
            INFO("Predict %d, Confidence = %f", ilabel, prob);

            printf("Pass [Enter] to next.");
            getchar();
        }
        return 0;
    }
};

int main(){
    return Application::run();
}