八. 实战：CUDA-BEVFusion部署分析-学习spconv的优化方案(Explicit GEMM conv)
目录
前言
自动驾驶之心推出的 《CUDA与TensorRT部署实战课程》，链接。记录下个人学习笔记，仅供自己参考
本次课程我们来学习下课程第八章——实战：CUDA-BEVFusion部署分析，一起来学习 spconv 的优化方案（Explicit GEMM conv）
Note：之前在学习杜老师的课程中有简单记录过 Sparse Convolution 的一些基础知识，感兴趣的可以看下：复杂onnx解决方案（以sparseconv为例）
课程大纲可以看下面的思维导图
0. 简述
本小节目标：理解 im2col，Explicit GEMM conv 是什么，以及理解 spconv 中使用 Gather 和 ScatterAdd 做优化的用意是什么
这节给大家讲解第八章节第 4 小节，学习 spconv 的优化方案，那这一小节我们从 Explicit GEMM Conv 显式 GEMM Conv 方式去看 spconv 是怎么优化的
1. 什么是Explicit GEMM Conv
在学 spconv 优化方案之前我们得先去理解 Explicit GEMM Conv 是什么，那要去理解 Explicit GEMM Conv 之前我们得先去理解 im2col 是什么，一层一层递进，
首先我们看一下 GEMM 是什么，GEMM 全称是 GEneral Matrix Multiplication 通用矩阵乘法，我们在第二章 CUDA 中也有涉及到，那 Explicit GEMM Conv 就是显式矩阵乘法卷积，通常指把 conv 计算的输入输出通过 im2col 的方式转换成 matrix 来进行矩阵乘法，可以高效的实现卷积（CUDA 中常用的方法，因为 CUDA 非常擅长 MxN=P 的矩阵乘法）
Conv 卷积操作过程中有权重，有输入输出维度，且维度可以是多维。我们在 CUDA、cuDNN、cuTLASS 去做卷积的时候，内部实际的计算其实是一个叫 im2col 的方式
我们可以把卷积的 input 输入铺平成一个二维的 Matrix 矩阵，这是因为任何维度的东西可以按照某一个方向，某一个顺序展开成一个二维的。比如三维的 CxHxW 可以变成二维的，四维的 NxCxHxW 也可以变成二维的
那同理 conv 的权重也是可以铺平成二维矩阵的，卷积的 weight 维度一般是 C_outxC_inxK_hxK_w，其中 Kh 和 Kw 为 Kernel 的高和宽，四个维度如果我们按照一定的方式去展开的话也是可以变成一个二维的
那所以说我们通过 im2col 可以把 conv 的输入和权重都变成二维的 Matrix，那变成二维之后，这个其实就是开始变成 CUDA 比较擅长的领域，也就是矩阵乘法
那矩阵乘法我们在第二章中讲过可以使用 shared memory 等多种优化方式，那所以就意味着我们最终需要做的是什么呢，就是把 tensor 级别上的卷积转换成 MxN=P 这么一个矩阵乘法去做计算，那这个其实就是 Explicit GEMM Conv
下面是来自 ChatGPT 关于 Explicit GEMM Conv 的解释：
卷积操作在本质上是一种多维数组的加权求和操作，然而，卷积操作本身在计算上可能并不高效，特别是在计算机硬件上执行时。为了克服这一问题，卷积操作可以被重新构造为矩阵乘法操作，这一方法在效率上可以得到很大的提升
因此，Explicit GEMM Conv 就是指通过将输入数据和卷积核重排成适合矩阵乘法的形式，然后利用 GEMM 算法来执行卷积操作
2. im2col
im2col 字面的意思是 image to column 即将图像转换为列，这个过程涉及将输入图像数据转换为列矩阵的形式，以便可以通过通用矩阵乘法（GEMM）来实现卷积操作。
使用 im2col 方法，我们可以将图像中每个卷积窗口的像素展开成一个列向量。将这些列向量堆叠起来，我们就可以得到一个大的矩阵，其中每列代表一个卷积窗口。然后卷积核也被展开成一个行向量，并与上述矩阵进行乘法运算，实现了卷积操作。
这样，原本复杂的卷积操作就被转换成了两个矩阵之间的乘法，这是现代计算硬件（特别是使用 cuBLAS 库）擅长的操作
OK，下面我们来看 im2col 的案例，看一下卷积是怎么一步步变成矩阵乘法运算的
我们先来看一个最简单的案例，如下图所示：
在上图中我们的输入 input 维度是 1x5x5，filter 维度是 1x1x3x3，输出 output 维度是 1x3x3，stride 等于 1，padding 等于 0，我们先不看 batch 维度
从下图中可以看出卷积从头到尾 filter 一共需要滑动 9 次，每滑动一次 input 里面的 9 个元素都需要跟 filter 中的 9 个元素做一个乘加，得到一个输出数据
那这也就意味着整个计算其实等价于一个 1x9 乘以一个 9x1 的矩阵乘法，那么所以说我们可以把这个计算换成下图这个样子：
图中上面是 convolution，下面我们把它展平，我们把 filter 中的 9 个数据给它展开成一行，同理把 input 中每次滑动窗口的 9 个数据按列展开，那么它就变成了一行乘以一列，这样就得到了输出 9 个数据中的第 1 个数据，对于输出我们也将它按行展开
以上就是 filter 在 input 中滑动第一次的一个计算，那么同理第二次的计算中，filter 本身是不变的，变的是 input 中跟 filter 相乘的那一块，那么依此类推 filter 滑动 9 次之后就得到了完整的输出。那通过这样一个过程我们就可以把 conv 的计算等价成一个 1x9 乘以 9x9 的矩阵乘法
OK，我们理解了这个之后我们再稍微复杂一点，我们把 input 的 channel 设置成 2，那么 filter 的输入 channel 也就是 2，那么整个过程就变成了下面的样子：
filter 它之前不是 9 个数据吗，那么现在我们其实有两个 channel，因此它按行展开就有 18 个数据，同理 input 中滑动窗口的 18 个数据需要在列上做展开，那么最终二者计算得到输出中的一个数值，如下图所示：
那 filter 它最终滑动了多少次呢，那还是 9 次，这也就意味着 input 的宽度还是 9，不同的是高度变成了 18
那么最后我们把所有数据全展平后我们就得到一个 18x9 的矩阵，那整个 conv 过程就是 1x18 乘以 18x9 得到 1x9 的矩阵运算，那 output 的大小不变，只不过参与计算的数据变多了
OK，我们再复杂化一点，我们把 filter 的个数增加，增加成 3 个，那整个过程就变成了下面的样子：
如果说只有一个 filter 我们按行展开就行了，那现在多个 fiter 意味着每个 filter 就是一行，最终展开变成了如下的样子：
从图中我们能看出不再是 vector 和 matrix 相乘了，而是 matrix 和 matrix 相乘，这是因为我们的 filter 也变成了一个矩阵，维度是 3x18，与 input 中的 18x9 维的数据相乘得到最终的 output 输出
以上就是 im2col 的简单分析了，下面我们从公式的角度来分析下
我们规定下：
- I C IC IC， I H IH IH， I W IW IW 是 input 的 c c c， h h h， w w w 的大小
- K H KH KH， K W KW KW 是 kernel 的 h h h， w w w 的大小
- O C OC OC， O H OH OH， O W OW OW 是 output 的 c c c， h h h， w w w 的大小
那么 filter 每滑动一次时，input 中与 filter 参与计算的数量是：
- I C × K H × K W = 18 IC\times KH\times KW =18 IC×KH×KW=18
filter 滑动的次数是：
- ( I H ? K H + s t r i d e ) × ( I W ? K W + s t r i d e ) = 9 (IH-KH+stride)\times(IW-KW+stride)=9 (IH?KH+stride)×(IW?KW+stride)=9
filter 的个数是：
- O C OC OC
由此可见，我们是可以把参与 conv 计算的激活值、权重以及输出给放在一个矩阵里，用矩阵乘法的方式去计算 conv 即 M × N = P M\times N=P M×N=P
M M M 激活值的大小：
- O C × ( I C × K H × K W ) OC\times (IC\times KH\times KW) OC×(IC×KH×KW)
N N N 权重的大小：
- ( I C × K H × K W ) × ( O H × O W ) (IC\times KH\times KW) \times (OH\times OW) (IC×KH×KW)×(OH×OW)
P P P 输出的大小：
- O C × ( O H × O W ) OC\times (OH\times OW) OC×(OH×OW)
那其实我们在做计算的时候就非常适合 CUDA 加速了，这个模式大家可能都已经看过很多遍了，就是每一个线程负责一个点，那每一个点其实就是利用 shared memory 做一个乘加
Explicit GEMM Conv 就是显性的分配额外的空间用来做 im2col 的处理，将 N 维的数据转为 2 维，并用优化过的矩阵乘法算法来加速
3. spconv是如何使用Explicit GEMM Conv的
谈到 spconv 的加速我们不得不聊一个 repo，就是 https://github.com/traveller59/spconv
这个 repo 实现的 spconv 是相当不错的，现在都有很多人在用它去加速，它这里面采用的方式有 Explicit GEMM Conv，对稀疏矩阵乘法做了很大程度上的加速，最新的 spconv2.3 已经支持 int8
我们可以看到这里面作者其实就是充分利用了 CUDA 的一些特性去做一个加速，那现在整个更新已经到 CUDA-12.0 了，还是比较新的，我们可以直接用 pip 方式去安装
我们可以看他现在的 spconv 已经更新到 2.3，可以支持 int8 的量化，那里面有一些使用案例，大家感兴趣的可以看下
同时我们之前不是做了 CenterPoint 的环境搭建吗，那 CenterPoint 中 SCN 网络的 spconv 加速其实就是用它这个库来进行加速的
所以由此可见这个 spconv 它影响力还是比较大，那大家安装之后可以根据它的 README example 先跑一跑一些示例，看看里面的 spconv 加速是怎么实现的
OK，我们下面再来看 spconv 是如何使 Explicit GEMM Conv 做加速的，官方提供了一个文档用来说明 spconv 算法的加速，具体在 spconv/docs/spconv2_algo.pdf 中，我们一起来看看他是怎么做的
这里面的 input 输入数据是 5x5 的，filter 是 3x3 的，其中 input 的 25 个数据只有 5 个数据是有效的，其它全为 0，那它来跟这个 filter 来做计算得到一个 output，其中 output 中有 6 个有效数据
那么 filter 在这里通过滑动窗口来进行计算，那滑动窗口中只要有一个数据在 filter 里面，也就是 input 滑动窗口中只要有一个数据被 filter 盖住了，那它就参与计算
那这里我们可以看到 filter 有 6 次 的滑动会捕捉到信息，那这里面有个点大家不能忽视，就是我们要对 input 和 filter 通过 im2col 展开，那展开的话其实我们能发现 input 中存在大量的无效数据也就是 0
也就是说我们转换成矩阵运算之后，虽然可以把很多计算给省掉，但是依旧存在很多无用的计算，主要原因在于我们的 input 是稀疏的，转换成的矩阵也是稀疏的
那所以作者就说 Too much zeros!，这里面的 0 太多了，我们需要想办法把这些 0 给去除掉
那作者是怎么做的呢，这里面主要是两个步骤即 Gather 和 ScatterAdd，就是用 Gather 和 Scale 的方式来做一个压缩和扩充，如下图所示：
我们每一列可能 0 特别多，那么我们就可以通过 Gather 把里面的数据给压缩，把一列中所有的 0 全给去除掉就只保留没有 0 的部分，比如图中 5 个数据有 2 个是 0，通过 Gather 我们就可以把这 2 个 0 给去除掉
接着我们让它再跟权重去做计算，那无用的计算就会少很多，通过这种方式得到输出值。此外我们还需要保证输入和输出的维度大小一致，因此我们还要把之前删除的 0 再加上去，通过 ScatterAdd 的方式来完成
所以总结下来 traveller59 作者
traveller 提供的解决方案是使用 Gather 和 ScatterAdd。**通过 Gather 将 input 中参与计算但没有意义的 0 点数据去除掉，**只留下有用的数据，这样 GEMM 的计算会变得更加 Dence，GEMM 计算完了以后，再通过 ScatterAdd 将 0 点添加进去。那这个就是作者提供的用 Explicit GEMM Conv 做 spconv 加速的一个方案
这里博主有些困惑，那在 filter 的 9 次滑动中其实 6 次是有效的，但不知道为什么图中只计算了 5 次，博主按照自己的理解，绘制了一个草图，描述了整个过程，如下图所示：
4. 使用Explicit GEMM Conv处理spconv的优缺点
OK，我们讲完 Explicit GEMM Conv 之后，我们思考一下它有什么缺点，我们说 Explicit 代表显性的意思，那显性去做 im2col 就意味着我们需要分配额外的空间，那额外的开销它其实是不能忽略的，在某种意义上是会产生很大的延迟的
那这部分其实我们也需要去考虑是否能够去除掉，因为它导致会有很多没有必要的 memory R/W，那这个地方就很容易成为瓶颈，那这个就是显式做 GEMM Conv 的一个问题
那有了这个问题，那么肯定就会有很多人去想怎么去解决它，那么所以有一个相应的概念叫做 Implicit GEMM Conv，隐式的 GEMM Conv
这里先做个预告，如上图所示，我们可以看下显式 GEMM Conv 把 N 维的 tensor 给转换成矩阵，然后通过 GEMM 计算得到一个值，那这个是显式的
那隐式是怎么做的呢，隐式我们可以发现它没有通过 im2col 这个操作去生成矩阵，那它直接就是在 tensor 这个维度上去做矩阵乘法
矩阵乘法中的每一个数据我们不再给给它分配空间了，我们直接通过索引的方式去寻找这个矩阵方法所需要的计算点到底对应 tensor 中的哪一个点，那这样通过索引的方式去寻找计算点可以直接跳过 im2col 这个步骤
那这个方式其实也是 spconv 可以选择的一种加速方式，那这部分具体怎么做，我们下节课再给大家讲解
OK，本次课程到这里就结束了
5. 拓展-conv加速
5.1 Introduction
这里我们将展示一种将卷积运算转换为矩阵乘法的方法。这种方法的优点是计算速度更快，但内存使用量更大。我们采用 im2col 操作，将输入图像转换为矩阵，然后将该矩阵与重塑后的 Kernel 相乘。最后我们再通过 col2im 运算将乘积矩阵重塑为图像
5.2 im2col
按照正常方式我们需要使用大量 for 循环来实现卷积，虽然这有助于理解，但速度不够快，这里我们将学习如何以矢量化方式实现卷积。
首先，如果我们仔细观察，卷积运算基本上是 Kernel 与移动窗口选择的局部区域之间的点积，而移动窗口选择的局部区域大小与我们的 Kernel 相同，如果我们在内存中扩展所有可能的窗口，并将点乘作为矩阵乘法来执行，会发生什么情况？答案是速度提高 200 倍或更多，但会消耗更多内存
比如，如果输入值为 [227x227x3]，需要用步长为 4、填充为 0 的 11x11x3 滤波器进行卷积，那么我们就需要在输入值中提取 [11x11x3] 个像素块，然后将每个像素块拉伸为大小为 11*11*3=363 的列向量
以输入 227 计算，步长为 4，填充为 0，则宽度和高度上都有 ((227-11)/4)+1=55 个位置，从而得到大小为 [363x3025] 的输出矩阵 X_col。在这里，每一列都是一个拉伸的感受野，总共有 55*55=3025 个
总结一下我们如何计算 im2col 的输出大小：
[img_height, img_width, img_channels] = size(img);
newImgHeight = floor(((img_height + 2*P - ksize) / S)+1);
newImgWidth = floor(((img_width + 2*P - ksize) / S)+1);
cols = single(zeros((img_channels*ksize*ksize),(newImgHeight * newImgWidth)));
CONV 层的权重也同样被拉伸成行。例如，如果有 96 个大小为 [11x11x3] 的滤波器，那么矩阵 W_row 的大小为 [96x363]，其中 11x11x3=363
图像和 Kernel 转换成矩阵后，卷积可以通过简单的矩阵乘法实现，在我们的例子中就是 W_col[96x363] 乘以 X_col[363x3025]，得到一个矩阵 [96x3025]，需要将其重塑为 [55x55x96]
最后的重塑也可也通过一个名为 col2im 的函数来实现
值得注意的是，im2col 的某些实现会将结果转置，如果是这种情况，则必须改变矩阵乘法的顺序
5.3 Forward graph
为了帮助我们使用 im2col 进行卷积，并推导出反向传播，我们以图形的形式展示了使用 im2col 进行的卷积，如下图所示。这里的输入张量是单一的 3 通道 4x4 图像，它将进入一个卷积层，卷积层的 S:1 P:0 K:2 F:1（输出量）
5.4 Backward graph
使用 im2col 技术，计算图与 FC 层详细，有着相同的公式 f ( x , θ , β ) = ( x . θ T ) + β f(x,\theta,\beta)=(x.\theta^T)+\beta f(x,θ,β)=(x.θT)+β，不同的是现在我们有了 reshape、transpose、im2col 块
关于反向传播过程中的 reshape 和 transpose，你只需要使用另一种 reshape 或 transpose 来反转它们的操作即可，需要注意的是，如果在前向传播过程中使用的是行为主的 reshape，那么在反向传播过程中也需要使用行为主的 reshape
唯一需要注意的是 im2col 的反向传播操作，它不能使用简单的 reshape 来实现，因为 patches 实际上可能会重叠（取决于步长），因此需要对 patches 相交处的梯度求和
5.5 Python example for forward propagation
def conv_forward_naive(x, w, b, conv_param):
"""
A naive implementation of the forward pass for a convolutional layer.
The input consists of N data points, each with C channels, height H and width
W. We convolve each input with F different filters, where each filter spans
all C channels and has height HH and width HH.
Input:
- x: Input data of shape (N, C, H, W)
- w: Filter weights of shape (F, C, HH, WW)
- b: Biases, of shape (F,)
- conv_param: A dictionary with the following keys:
- 'stride': The number of pixels between adjacent receptive fields in the
horizontal and vertical directions.
- 'pad': The number of pixels that will be used to zero-pad the input.
Returns a tuple of:
- out: Output data, of shape (N, F, H', W') where H' and W' are given by
H' = 1 + (H + 2 * pad - HH) / stride
W' = 1 + (W + 2 * pad - WW) / stride
- cache: (x, w, b, conv_param)
"""
out = None
pad_num = conv_param['pad']
stride = conv_param['stride']
N,C,H,W = x.shape
F,C,HH,WW = w.shape
H_prime = (H+2*pad_num-HH) // stride + 1
W_prime = (W+2*pad_num-WW) // stride + 1
out = np.zeros([N,F,H_prime,W_prime])
#im2col
for im_num in range(N):
im = x[im_num,:,:,:]
im_pad = np.pad(im,((0,0),(pad_num,pad_num),(pad_num,pad_num)),'constant')
im_col = im2col(im_pad,HH,WW,stride)
filter_col = np.reshape(w,(F,-1))
mul = im_col.dot(filter_col.T) + b
out[im_num,:,:,:] = col2im(mul,H_prime,W_prime,1)
cache = (x, w, b, conv_param)
return out, cache
5.6 Python example for backward propagation
def conv_backward_naive(dout, cache):
"""
A naive implementation of the backward pass for a convolutional layer.
Inputs:
- dout: Upstream derivatives.
- cache: A tuple of (x, w, b, conv_param) as in conv_forward_naive
Returns a tuple of:
- dx: Gradient with respect to x
- dw: Gradient with respect to w
- db: Gradient with respect to b
"""
dx, dw, db = None, None, None
x, w, b, conv_param = cache
pad_num = conv_param['pad']
stride = conv_param['stride']
N,C,H,W = x.shape
F,C,HH,WW = w.shape
H_prime = (H+2*pad_num-HH) // stride + 1
W_prime = (W+2*pad_num-WW) // stride + 1
dw = np.zeros(w.shape)
dx = np.zeros(x.shape)
db = np.zeros(b.shape)
# We could calculate the bias by just summing over the right dimensions
# Bias gradient (Sum on dout dimensions (batch, rows, cols)
#db = np.sum(dout, axis=(0, 2, 3))
for i in range(N):
im = x[i,:,:,:]
im_pad = np.pad(im,((0,0),(pad_num,pad_num),(pad_num,pad_num)),'constant')
im_col = im2col(im_pad,HH,WW,stride)
filter_col = np.reshape(w,(F,-1)).T
dout_i = dout[i,:,:,:]
dbias_sum = np.reshape(dout_i,(F,-1))
dbias_sum = dbias_sum.T
#bias_sum = mul + b
db += np.sum(dbias_sum,axis=0)
dmul = dbias_sum
#mul = im_col * filter_col
dfilter_col = (im_col.T).dot(dmul)
dim_col = dmul.dot(filter_col.T)
dx_padded = col2im_back(dim_col,H_prime,W_prime,stride,HH,WW,C)
dx[i,:,:,:] = dx_padded[:,pad_num:H+pad_num,pad_num:W+pad_num]
dw += np.reshape(dfilter_col.T,(F,C,HH,WW))
return dx, dw, db
5.7 Im2col and Col2im sources in python
该实现将接收三维张量 [channels, rows, cols]，并创建二维矩阵 [rows=(new_h*new_w), cols=(kw*kw*C)]，注意该算法将输出上图的转置版本
def im2col(x,hh,ww,stride):
"""
Args:
x: image matrix to be translated into columns, (C,H,W)
hh: filter height
ww: filter width
stride: stride
Returns:
col: (new_h*new_w,hh*ww*C) matrix, each column is a cube that will convolve with a filter
new_h = (H-hh) // stride + 1, new_w = (W-ww) // stride + 1
"""
c,h,w = x.shape
new_h = (h-hh) // stride + 1
new_w = (w-ww) // stride + 1
col = np.zeros([new_h*new_w,c*hh*ww])
for i in range(new_h):
for j in range(new_w):
patch = x[...,i*stride:i*stride+hh,j*stride:j*stride+ww]
col[i*new_w+j,:] = np.reshape(patch,-1)
return col
def col2im(mul,h_prime,w_prime,C):
"""
Args:
mul: (h_prime*w_prime*w,F) matrix, each col should be reshaped to C*h_prime*w_prime when C>0, or h_prime*w_prime when C = 0
h_prime: reshaped filter height
w_prime: reshaped filter width
C: reshaped filter channel, if 0, reshape the filter to 2D, Otherwise reshape it to 3D
Returns:
if C == 0: (F,h_prime,w_prime) matrix
Otherwise: (F,C,h_prime,w_prime) matrix
"""
F = mul.shape[1]
if(C == 1):
out = np.zeros([F,h_prime,w_prime])
for i in range(F):
col = mul[:,i]
out[i,:,:] = np.reshape(col,(h_prime,w_prime))
else:
out = np.zeros([F,C,h_prime,w_prime])
for i in range(F):
col = mul[:,i]
out[i,:,:] = np.reshape(col,(C,h_prime,w_prime))
return out
def col2im_back(dim_col,h_prime,w_prime,stride,hh,ww,c):
"""
Args:
dim_col: gradients for im_col,(h_prime*w_prime,hh*ww*c)
h_prime,w_prime: height and width for the feature map
strid: stride
hh,ww,c: size of the filters
Returns:
dx: Gradients for x, (C,H,W)
"""
H = (h_prime - 1) * stride + hh
W = (w_prime - 1) * stride + ww
dx = np.zeros([c,H,W])
for i in range(h_prime*w_prime):
row = dim_col[i,:]
h_start = (i / w_prime) * stride
w_start = (i % w_prime) * stride
dx[:,h_start:h_start+hh,w_start:w_start+ww] += np.reshape(row,(c,hh,ww))
return dx
5.8 Smaller example
为了让问题变简单，我们以 X[3x3] 与 W[2x2] 的卷积为例进行说明
总结
这节课程我们学习了 Explicit GEMM Conv 去加速 spconv 的方案，我们首先介绍了 im2col，卷积的计算实际上是可以转换为矩阵乘法运算的，接着我们分析了使用 Explicit GEMM Conv 方法加速 spconv，主要是通过 Gather 和 ScatterAdd 两个步骤来完成的，最后我们分析了 Explicit GEMM Conv 的缺点主要是需要额外分配内存空间，内存的读写非常耗时，因此引出了 Implicit GEMM Conv 的加速方案
OK，以上就是第 4 小节有关 Explicit GEMM conv 优化方案的全部内容了，下节我们将去学习 spconv 另一种优化方案即 Implicit GEMM conv，敬请期待😄
下载链接
参考
本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符，请联系我的编程经验分享网邮箱：chenni525@qq.com进行投诉反馈，一经查实，立即删除！