动手用CUDA实现CNN

项目地址： github.com/hijkzzz/neu…

作为深度学习领域的初学者，很多人会好奇TensorFlow和PyTorch等深度学习框架是怎样实现高效反向传播算法的。毫无例外，这些系统都使用了CUDA进行并行计算加速，在此我用CUDA实现了一个简单的CNN网络，并在无扭曲的MNIST数据集上实现了99.23%的准确率。

预备知识

CUDA

这里推荐 《CUDA编程极简入门教程》

矩阵求导

推荐知乎大神分享的《矩阵求导术》

设计

存储

thrust.github.io/

在实现神经网络之前，我们需要设计一个存储类，用于保存GPU上的参数和数据。这里称为Storage类，为了方便实现，我们直接使用CUDA提供的thrust::device_vector（类似于std::vector）管理显存上的动态数组。并增加一个std::vector保存Storage的形状，可以理解为TensorFlow中Tensor的形状。

矩阵乘法

神经网络的实现大量用到矩阵乘法，所以CUDA并行加速的一个关键在于实现高效的并行矩阵乘法。这里我直接使用了《CUDA编程极简入门教程》中Shared Memory加速的矩阵乘法。实际上还可以继续优化，使效率大大提升。

全连接层

设X为输入数据矩阵，其中每一行为一个样本。W为参数矩阵，b为偏置向量，L为样本平均损失。* 表示矩阵乘法，而非逐元素相乘，^T表示转置：

全连接

前向传播
Y = X * W
反向传播
dL/dX = dL/dY * W^T
dL/dW = X^T * dL/dY

偏置

前向传播
Y = X + b
反向传播
dL/db = sum(dL/dY, 0) 逐样本梯度求和
复制代码

卷积层

hal.inria.fr/file/index/…

为了方便用矩阵乘法实现卷积，我参考了Caffe的卷积原理，即im2col：

基本的思想是把卷积运算展成矩阵乘法，所以可以用并行加速的矩阵乘法高效实现卷积。设F为卷积核参数，且形状为： channel_out*channel_in*kernel_width*kernel_height ，X为一个输入样本形状为 channel_in*width*height ，b为偏置向量。

卷积
前向传播
col = im2col(im) 根据im2col展开输入图
Y = F * col
反向传播
dL/dF = dL/dY * col^T
dL/d_col = F^T * dL/dY
dL/d_im = col2im(dL/d_col)

偏置
前向传播
Y = X + b 逐通道相加
反向传播
dL/db = sum(sum(X, 2), 1) 对整个通道进行规约
复制代码

Maxpool

Maxpool的反向传播需要记录池化前元素的位置，然后把反向梯度直接传回

激活函数

激活函数的前向反向传播都是一样的

ReLU
前向传播
Y = relu(X)
反向传播
dL/dX = relu'(X) element_mul dL/dY 逐元素相乘
其中relu'(x) = 1 if x > 0 else 0 

Sigmoid
前向传播
Y = sigmoid(X)
反向传播
dL/dX = sigmoid'(X) element_mul dL/dY 逐元素相乘
其中 sigmoid'(x) = sigmoid(x) * (1 - sigmoid(x))
复制代码

Softmax

在工程实现上：为了防止Softmax的分母溢出，一般使用LogSoftmax代替。设定1_n为全为1的列向量

Logsoftmax
正向传播
Y = log_softmax(X) = x - log(exp(X) * 1_n) * 1_n^T

由前言中矩阵求导的方法可得
反向传播
dL/dX = dL/dY - (dL/dY * 1_n * exp(x)) / (exp(x) * 1_n)

复制代码

NLLLoss

NLLLoss是平均负的对数似然损失，为了配合LogSoftmax使用而实现。设Y为样本标签矩阵，每一行为一个样本。N为样本数量

前向传播
L = mean(sum(-log_P element_mul Y, 1), 0)

反向传播
用矩阵乘法，L可表示为 L = 1_n^T * ((-log_P element_mul Y) * 1_k) / N
由矩阵求导术可得
dL/d(log_P) = -Y / N

NLLLoss+LogSoftmax为我们常见的Softmax损失
将dL/d(log_P)带入LogSoftmax梯度中可得softmax损失的梯度： softmax(X) - Y
复制代码

RMSProp

为了实现单独的优化器，我们需要在反向传播的时候把梯度保存下来，然后用RMSProp算法进行统一的滑动平均计算新梯度。同理可以很方便的实现Adam等优化器。

实现

源码结构

src
    cuda        CUDA源码
    minist      MNIST DEMO
test
    cuda        CUDA源码单元测试
CMakeLists.txt  CMake编译脚本
复制代码

由于篇幅有限，所以这里只能去看GitHub上的实际代码。每个层都封装为了一个类，并且可调用connect函数连接层与层。

Debug/调优

可以通过CUDA提供的Visual Profiler可以很方便的看出程序的性能瓶颈。

在我的实验中发现80%的执行时间都在等待显卡I/O，所以通过Pinned Memory以及合并传输/内存分配等方式使运行效率提升了数十倍。 其次是矩阵乘法还有较大的优化空间，不过总的来说在GTX1070上数十秒便可以跑完MNIST的6W个样本，基本实现了我的目标。

编程一天，调试两天，Debug是开发的一个困难而且重要的环节，掌握适当的 工具 的方法将事半功倍。CUDA提供的Nsight、cuda-memcheck都是很好的工具。当然printf+注释大法也是屡试不爽。

测试

网络结构
conv 1 32 5 relu
maxpool 2
conv 32 64 5 relu
maxpool 2
conv 64 128 3 relu
fc 4 * 128 128 relu
fc 128 10 relu
softmax
nllloss

调参
shuffle = true
batch_size = 128
learning_rate = 0.003
L2 = 0.0001
beta = 0.99

准确率
1 epoch 93%
10 epochs 99.12%
30 epochs 99.23%
10s / epoch(GTX1070)
复制代码