ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

从零实现一个C++机器学习库:从Tensor到自动求导的完整实践

从零实现一个C++机器学习库:从Tensor到自动求导的完整实践 如果你和我一样在一门算法课或者一个实际项目中碰见了“C机器学习库开发”这种需求心里多半会冒出一连串问题Python 不香吗为什么非得用 C 去折腾这些事说实话我自己在动手之前也是这个想法但等到要把模型部署到没有 Python 环境的服务器上、要跑在资源受限的嵌入式设备里、或者要把训练好的模型嵌进现成的 C 业务系统中才发现一个轻量、可控、不依赖巨型框架的 C 机器学习库确实是绕不开的刚需。这篇文章就是基于我从零开发一个 C 机器学习库的完整经历聊一聊库的边界设计、数据结构怎么选、自动求导怎么实现、性能优化和工程化落地。内容适合已经有一定 C 基础、想从“调包”走向“实现”的读者也适合在学校做课程设计、毕业设计时参考。1. 先给库定位不做的东西和要做的东西同样重要1.1 为什么不用现成框架而要自己写一个 C 机器学习库如果你只是做实验、跑模型PyTorch 和 TensorFlow 无疑是更高效的选择。但如果你是在一个 C 为主语言的团队里做推理服务直接引第三方深度学习框架往往并不舒服依赖体积巨大光是 CUDA、MKL 这些运行时就能让安装包膨胀到几个 GB。框架的 C API 虽然存在但设计上始终以 Python 为第一公民很多高级特性封装得不够顺手。在嵌入式设备或特殊操作系统上很多框架无法编译运行。某些业务场景要求我们完全掌控模型结构和前向推理逻辑比如要做模型的权限管理、做特征对齐、做算子融合用框架反而很别扭。所以我当时的目标并不是“重新发明一个 TensorFlow”而是写一个满足特定场景的轻量级 C 机器学习库能够很方便地嵌入到现有工程里并且核心逻辑自己可控。明确这个边界后后面的一切设计才不会跑偏。1.2 功能清单第一版只做真正必要的模块在动笔之前我先列出了第一版的功能清单多维数组与基础矩阵运算加法、乘法、转置、逐元素运算。自动求导支持标量和张量级别的反向传播能完成梯度累积和拓扑排序。常见网络组件线性层全连接层、激活函数层ReLU、Sigmoid、Tanh。损失函数均方误差 MSE、交叉熵 CrossEntropy。优化器SGD、Momentum、Adam。同时明确不做卷积、循环神经网络、分布式训练、GPU 加速、模型序列化与转换工具。这些功能并不是不重要而是第一版的核心价值应该是先把“训练一个小型网络”的闭环跑通。等这个闭环稳定了后续再按实际需求扩展远比一开始就铺开要靠谱。1.3 模块划分一个让别人能看懂的结构库的代码组织直接决定后续扩展和维护的效率。我采用的目录结构是这样ml_lib/ ├── include/ml/ │ ├── tensor.h // 多维数组 │ ├── autograd.h // 自动求导核心 │ ├── layers.h // 网络层 │ ├── loss.h // 损失函数 │ ├── optimizer.h // 优化器 │ └── model.h // 模型容器 ├── src/ │ ├── tensor.cpp │ ├── autograd.cpp │ ├── layers.cpp │ ├── loss.cpp │ └── optimizer.cpp ├── tests/ │ ├── test_tensor.cpp │ ├── test_autograd.cpp │ └── test_layers.cpp ├── examples/ │ └── linear_regression.cpp └── CMakeLists.txt这个结构看起来很简单但每个模块之间的依赖关系必须清晰tensor 是基础autograd 依赖 tensorlayers 和 loss 依赖 autogradoptimizer 依赖 tensor 和 autogradmodel 负责把前向过程串起来。这个方向一旦依赖反过来后面写起来就会痛苦不堪。2. 数据结构的地基多维数组到底该怎么存2.1 别急着写类先把内存布局想清楚机器学习库最底层、最重要的数据结构就是张量。很多初学者会用std::vectorstd::vectorfloat直接表示二维矩阵因为写起来直观。但这种做法的性能问题很大每一行都是一个独立的堆分配多行数据在内存中并不是连续存放的CPU 缓存命中率会非常差。我当时在做一个大矩阵乘法的测试时对比了vectorvector和连续一维数组两种方式在 1024x1024 的矩阵规模下后者快了将近一倍。原因很简单连续内存可以一次加载进缓存而离散的二维结构会导致频繁 cache miss数据量越大差距越明显。这也就是热词里“多维数组 C 指针”背后真正要解决的问题如何在不牺牲访问效率的前提下让多维数组用起来仍像 a(i, j) 这样直观。答案是用一维数组存储配合指针和索引计算来模拟多维访问。2.2 Tensor 接口设计的核心实现下面这个简化版 Tensor 展示了最核心的思路template typename T class Tensor { public: Tensor() default; Tensor(int rows, int cols) : rows_(rows), cols_(cols), data_(rows * cols, T(0)) {} T operator()(int row, int col) { return data_[row * cols_ col]; // 行主序 } const T operator()(int row, int col) const { return data_[row * cols_ col]; } int rows() const { return rows_; } int cols() const { return cols_; } T* data() { return data_.data(); } const T* data() const { return data_.data(); } private: int rows_ 0; int cols_ 0; std::vectorT data_; // 一维连续存储 };这里最关键的是row * cols_ col这个映射。行主序的意思是存储时先存完第一行、再存第二行。MATLAB 和 Fortran 默认是列主序Python 的 NumPy 默认是行主序C 社区也大多采用行主序。选择行主序的另一个原因是它和 C/C 的二维数组在内存中的排列方式一致将来如果要对接 BLAS、OpenBLAS 这类底层库数据可以直接通过data()指针传进去几乎零拷贝。在这个基础之上如果要支持更高维度做法其实一样比如一个形状为 [batch, channel, height, width] 的四维张量访问索引(n, c, h, w)时偏移量就是size_t offset (((n * channels c) * height h) * width w);一层一层按顺序把维度乘进去。核心思想是多维数组在代码里看起来是“多维”的但在内存里始终是一根连续的线。2.3 生命周期管理移动语义和引用计数不能忽略Tensor 一旦承载了较大的数据量拷贝的代价就极高。如果每次运算都产生临时对象再深拷贝训练一个稍大一点的模型就会非常慢。所以我在实现 Tensor 时做了一个很关键的决定默认禁用拷贝构造只提供Clone()显式拷贝同时实现移动构造和移动赋值。Tensor(const Tensor) delete; Tensor operator(const Tensor) delete; Tensor(Tensor other) noexcept : rows_(other.rows_), cols_(other.cols_), data_(std::move(other.data_)) { other.rows_ 0; other.cols_ 0; }这样写之后每一次运算返回的临时 Tensor 都能被高效“偷走”资源而不是把整块内存再复制一遍。实际跑训练循环时这个决定的收益非常明显。如果你打算在机器学习库的底层用裸指针加手动new/delete我建议放弃这种想法std::vector的内存管理已经足够好而且避免内存泄漏的烦恼。另外C 里还有一个容易踩的坑当你在多线程环境下共享 Tensor 时std::vector并不保证线程安全。如果未来要做并行推理需要对数据并发访问加锁或者干脆在实现上让每个线程持有自己的 Tensor。第一版先别考虑太复杂保证单线程正确后面再优化并行。3. 自动求导引擎把反向传播从黑盒变成白盒3.1 计算图的最小设计节点、边和反向函数机器学习库的自动求导核心就是维护一张计算图。每一个参与运算的变量都是一个节点节点与节点之间的运算关系就是边。前向传播时数据沿着计算图从输入流向输出反向传播时梯度沿着相反方向流回每个节点。在 C 里实现这个机制最简单的做法是让每个节点保存三样东西当前数值、累积梯度、以及一个反向计算函数。class ValueNode : public std::enable_shared_from_thisValueNode { public: float value 0.0f; float grad 0.0f; std::vectorstd::shared_ptrValueNode children; std::functionvoid() backward_fn; void Backward() { // 这里简化处理先重置梯度再调用 backward_fn // 实际要按拓扑序进行图中依赖多个子节点时需注意顺序 grad 1.0f; backward_fn(); } };当然面向张量的真实实现比这个复杂得多。你要考虑 Tensor 的梯度、维度规约、分支判断但思想是完全一致构建一个有向无环图反向传播时从根节点出发调用每个节点的反向函数把梯度传给它依赖的节点。3.2 一个最容易出错的细节梯度累积与清零反向传播中特别容易掉进坑的一个地方是梯度累积。假设某个参数被多个路径引用比如在一个网络里同一个权重矩阵被不同分支使用那么它收到的梯度应该是来自这些所有分支的梯度之和。如果反向函数里直接写grad upstream_grad后覆盖先梯度就丢了。正确做法是每次反向传播前先调用ZeroGrad()把梯度清零然后在反向过程中执行grad upstream_grad。很多学习者在给自己的库写反向时都会在这里遇到问题具体的表现是 loss 不减反增或者训练一开始就 NaN。NaN 的原因往往是梯度累积导致数值溢出而不是学习率太大。3.3 计算图的内存回收训练完一定要释放中间节点自动求导引擎如果采用动态图设计每次前向都会创建大量中间节点。训练一个 batch 结束后如果不及时释放这些节点内存占用会持续增长跑十几个 epoch 就可能把服务器内存吃光。我的经验是两点训练模式下每个 batch 前向之前记录一下计算图根节点反向完成后把不再需要的中间节点都清掉推理模式下完全不构建计算图前向过程中的临时结果直接复用一块预分配的内存。这样既能省内存又能让推理速度明显提升。这里我建议在动手写自动求导之前先去熟悉你熟悉框架源码里autograd目录的运作方式或者用 Python 的 NumPy 写一个小规模的自动求导 Demo。把逻辑想通之后再落到 C会省很多返工时间。4. 把算法落到代码层、损失函数和优化器的设计4.1 网络层的抽象如何管理参数在 C 面向对象的思想下网络层最自然的抽象是基类加子类。基类Layer声明 forward 和 backward 的接口每一个具体层只关心自己的运算和参数。class Layer { public: virtual ~Layer() default; virtual Tensor Forward(const Tensor input) 0; virtual void Backward(const Tensor grad_output) 0; virtual void UpdateParameters(Optimizer optimizer) 0; virtual std::vectorTensor* Parameters() 0; };这里有个值得注意的点参数更新不应该由 Layer 自己写死。不同的优化器有不同的更新规则所以 Layer 只负责暴露自己的参数指针具体怎么更新交给 Optimizer 去处理。这样后续添加新的优化器非常方便不用改动已有层。权重的初始化也是一个重要细节。全连接层权重如果不做初始化第一轮前向就会出现梯度爆炸或消失。我写的第一版用的是最简单的随机均匀分布初始化范围是[-0.05, 0.05]。后来换成 Kaiming 初始化对 ReLU 激活函数的网络训练稳定性提升非常明显。如果你以后要支持更深的网络推荐实现 Xavier 和 Kaiming 两套初始化策略。4.2 损失函数的数值稳定性细节决定成败均方误差相对简单float MSELoss::Forward(const Tensor pred, const Tensor target) { float sum 0.0f; int n pred.rows() * pred.cols(); for (int i 0; i n; i) { float diff pred.data()[i] - target.data()[i]; sum diff * diff; } return sum / n; // 除以样本数量才是均值 }交叉熵损失就要小心了。如果直接按公式-log(softmax(x)[class])计算当softmax的输出非常接近 0 时log(0)会产生-inf。解决办法是把 softmax 和交叉熵合并成一个函数利用数值技巧先减去最大值再计算指数既保持了数学等价性又避免了溢出。float softmax_cross_entropy(const std::vectorfloat logits, int target) { float max_val *std::max_element(logits.begin(), logits.end()); float sum 0.0f; for (float v : logits) { sum std::exp(v - max_val); } float log_sum_exp max_val std::log(sum); return log_sum_exp - logits[target]; }这种合并写法在很多框架内部都是标配。第一次实现时如果分开写 softmax 再算 loss精度和稳定性都会打折扣。4.3 优化器的实现SGD、Momentum 与 Adam 的更新公式优化器就是把反向传播得到的梯度作用到参数上。SGD 最简单for (size_t i 0; i params.size(); i) { params[i]-data()[j] - learning_rate * params[i]-grad()[j]; }但 SGD 收敛速度慢容易在局部振荡。Momentum 引入了速度项更新公式为velocity[i] momentum * velocity[i] learning_rate * grad; param[i] - velocity[i];Adam 则同时维护一阶矩估计和二阶矩估计并且加入偏差修正。实现 Adam 时一个容易踩的坑是要为每个参数维护两个额外的状态向量而且要在初始化时把m和v都置为 0。如果状态没有正确保存或者更新顺序写错训练曲线很容易出现一开始正常、后面剧烈震荡的情况。下面是一个精简版的 Adam 更新核心for (int i 0; i params.size(); i) { float g grad[i]; m[i] beta1 * m[i] (1 - beta1) * g; v[i] beta2 * v[i] (1 - beta2) * g * g; float m_hat m[i] / (1 - std::pow(beta1, t)); float v_hat v[i] / (1 - std::pow(beta2, t)); param[i] - learning_rate * m_hat / (std::sqrt(v_hat) epsilon); }注意这里的epsilon建议取1e-8太小会导致分母为零太大则会让更新步长缩小得过多。5. 性能优化实录从“能跑”到“跑得快”5.1 先测量再谈优化很多人在写 C 机器学习库时最容易犯的错就是“凭感觉优化”觉得哪段代码慢就改哪段结果半天下来没有实际收益。正确做法是先用性能分析工具定位热点。Linux 上可以用perf不想用命令行工具的话可以在工程里用简单的计时器对关键函数打点。我在完成了第一版之后跑了一个三层感知机在 MNIST 上的训练发现矩阵乘法耗时占了整体训练时间的 73%。这个结果说明优化矩阵乘法的优先级最高其他函数再优化也影响有限。5.2 减少内存分配对象池与缓存复用训练循环里最容易被忽视的性能杀手是频繁的堆内存分配。每一次Tensor运算都会产生新的临时对象如果这些对象在循环中不断创建和释放malloc/free 的耗时甚至会超过计算本身。我当时的优化方案是给常见运算增加“复用缓冲区”机制。具体做法是在 Layer 内部为前向和反向各准备一块缓存 Tensor每次前向后结果写入缓存而不是重新分配。虽然代码写起来没有原来干净但训练效率提升非常明显。5.3 并行化OpenMP 让矩阵乘法立刻提速C 做并行最简单的方式是 OpenMP。只要在编译器里开了对应选项在循环前加一行指令就可以多线程执行#pragma omp parallel for for (int i 0; i rows; i) { for (int k 0; k inner; k) { float a_ik A(i, k); for (int j 0; j cols; j) { C(i, j) a_ik * B(k, j); } } }我用 4 核机器测试这样改造后矩阵乘法速度大约提升了 2.8 倍。要注意的问题是并行循环内的C(i, j)写入操作如果存在竞争会得到错误结果。上面的 i 循环并行时每个 i 对应不同的行不会互相冲突这样是安全的。如果对 j 循环做并行就可能同时写同一个C(i, j)需要小心。5.4 SIMD 和第三方库的取舍如果只是学习自己写矩阵乘法没问题。但如果要投入生产环境建议直接链接 OpenBLAS 或 Eigen。这些库不仅实现了 SIMD 指令级优化还有向量化、循环分块等高级技巧手写版本很难追上它们的性能。我在做性能对比时手写版矩阵乘法 1024x1024 大约需要 40ms而链接 OpenBLAS 后只需要 6ms 左右。最终我给自己的库加了一个可选的后端如果编译时检测到 OpenBLAS就调用它的接口如果没有就退回内置的实现。这样一个库既能用于学习研究也具备了一定的实际部署价值。5.5 内存对齐与缓存友好最后内存对齐是一个容易被忽略的细节。如果数据没有对齐SIMD 指令会退化成慢速路径。最简单的方法是用alignas(64)修饰 Tensor 内部的存储或者使用std::vector时做一些额外处理。64 字节对齐通常是为了配合 cache line 的大小这一步优化对矩阵乘法的性能提升有 10% 到 20% 的贡献。6. 工程化落地构建、测试、示例代码一个都不能少6.1 构建系统CMake 是不二选择C 工程没有统一的标准包管理器CMake 虽然不是最优雅的但它是事实标准兼容所有主流平台。最小可用的 CMakeLists.txt 可以这样写cmake_minimum_required(VERSION 3.14) project(ml_lib VERSION 0.1.0 LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) add_library(ml_lib STATIC src/tensor.cpp src/autograd.cpp src/layers.cpp src/loss.cpp src/optimizer.cpp ) target_include_directories(ml_lib PUBLIC include) # 可选开启 OpenMP find_package(OpenMP) if(OpenMP_CXX_FOUND) target_link_libraries(ml_lib PUBLIC OpenMP::OpenMP_CXX) endif() # 测试 enable_testing() add_subdirectory(tests)如果你是在 VS Code 里开发只需要在.vscode/tasks.json中配置好 cmake 构建任务在launch.json中设置调试程序路径就可以比较舒服地调 C 代码。热词里经常有人问“vscode 配置 C/C 环境”这块需要特别注意的就是编译器路径和 CMake 工具链的匹配问题Windows 上推荐装 MinGW-w64 或微软的 CMake 工具。6.2 单元测试不测试的库等于埋雷机器学习库的数值计算不允许出现一点错误否则训练出来的模型可能完全不可用。我用 GoogleTest 为每个模块写了基础测试。比如矩阵运算测试TEST(TensorTest, AddTwoMatrices) { Tensor a(2, 2, 1.0f); Tensor b(2, 2, 2.0f); Tensor c a b; EXPECT_NEAR(c(0, 0), 3.0f, 1e-6); EXPECT_NEAR(c(1, 1), 3.0f, 1e-6); }自动求导的测试是关键中的关键我采用的做法是把数值梯度和解析梯度做对比对每个参数给一个很小的扰动用差分近似计算梯度然后和反向传播的梯度比较。这个测试能抓出大量实现细节上的 bug。6.3 完整示例用自己写的库训练线性回归最后用一个可以直接运行的例子展示整个库的用法#include ml/ml_lib.h #include iostream int main() { // 构造模拟数据 y 2x 1 Tensor X(4, 1), y(4, 1); X(0,0)0.0f; y(0,0)1.0f; X(1,0)1.0f; y(1,0)3.0f; X(2,0)2.0f; y(2,0)5.0f; X(3,0)3.0f; y(3,0)7.0f; Linear model(1, 1); MSELoss loss; SGD optimizer(model.Parameters(), 0.01f); for (int epoch 0; epoch 500; epoch) { optimizer.ZeroGrad(); Tensor pred model.Forward(X); float loss_value loss.Forward(pred, y); model.Backward(loss.Backward()); optimizer.Step(); if (epoch % 50 0) { std::cout epoch epoch , loss loss_value std::endl; } } // 期望学到的权重接近 2偏置接近 1 std::cout weight: model.Parameters()[0]-data()[0] , bias: model.Parameters()[1]-data()[0] std::endl; return 0; }如果你能跑通这段代码并且最终输出的 weight 在 2 左右、bias 在 1 左右说明库的底层矩阵运算、自动求导、网络层和优化器已经形成闭环。这个最小示例是检验整个库是否可靠最好的入门测试。做完这个项目之后我特别深的一个体会是实现一个 C 机器学习库本质上是在 C 的工程化和机器学习理论的交叉点上做大量选择。每一层设计都不是孤立的底层 Tensor 的存储方式直接影响自动求导的写法自动求导的接口设计又决定了上层 Layer 能写得有多自然。还有一个比较实际的建议如果你是在校学生完全可以把这个库作为课程设计或毕业设计的切入点。先从最小的矩阵运算开始逐步加入自动求导和一个简单的 MLP跑通 MNIST 手写数字识别就会对机器学习底层机制形成远比调包深刻的理解。整个过程中我最想提醒你的还是那句话别贪多先把一个简单的闭环走通再考虑功能扩展。
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进