【CNN】——矩阵乘法优化

转载自：通用矩阵乘（GEMM）优化与卷积计算
转载自：AI移动端优化之Im2Col+Pack+Sgemm
代码参考：BBuf/how-to-optimize-gemm

简介：
本文简要介绍通用矩阵乘（GEMM，General Matrix Multiplication）优化的基本概念和方法、QNNPACK 对特定场景的矩阵乘的优化方法、以及用 GEMM 优化神经网络中卷积计算的一点方向。

1. native

第一种方式就是
通用矩阵乘（下文简称 GEMM）的一般形式是 = C=AB，其中 A 和 B 涵盖了各自转置的含义。图一是矩阵乘计算中为计算一个输出点所要使用的输入数据。三个矩阵的形状也如图所示。

2. 1*4的拆分

将输出的计算拆分为 1×4 的小块，即将维度拆分为两部分。计算该块输出时，需要使用矩阵的 1 行，和矩阵的 4 列。

下面是该计算的伪代码表示，这里已经将 1×4 中 N 维度的内部拆分进行了展开。这里的计算操作数仍然是 2 ，这一点在本文中不会有变化。这里的内存访问操作数尚未出现变化，仍然是 4 ，但接下来会逐步改进。

for (int m = 0; m < M; m++) {for (int n = 0; n < N; n += 4) {C[m][n + 0] = 0;C[m][n + 1] = 0;C[m][n + 2] = 0;C[m][n + 3] = 0;for (int k = 0; k < K; k++) {C[m][n + 0] += A[m][k] * B[k][n + 0];C[m][n + 1] += A[m][k] * B[k][n + 1];C[m][n + 2] += A[m][k] * B[k][n + 2];C[m][n + 3] += A[m][k] * B[k][n + 3];}}
}

3. 4*4拆分

类似地，我们可以继续拆分输出的维度，从而在内侧循环中计算 4×4 输出

for (int m = 0; m < M; m += 4) {for (int n = 0; n < N; n += 4) {C[m + 0][n + 0..3] = 0;C[m + 1][n + 0..3] = 0;C[m + 2][n + 0..3] = 0;C[m + 3][n + 0..3] = 0;for (int k = 0; k < K; k++) {C[m + 0][n + 0..3] += A[m + 0][k] * B[k][n + 0..3];C[m + 1][n + 0..3] += A[m + 1][k] * B[k][n + 0..3];C[m + 2][n + 0..3] += A[m + 2][k] * B[k][n + 0..3];C[m + 3][n + 0..3] += A[m + 3][k] * B[k][n + 0..3];}}
}

4. 4 * 4 * 4拆分

到目前为止，我们都是在输出的两个维度上展开，而整个计算还包含一个削减（Reduction）维度。图五展示了在计算 4×4 输出时，将维度拆分，从而每次最内侧循环计算出输出矩阵的 4x4 部分和。

5. im2col + pack + sgemm

转载自:AI移动端优化之Im2Col+Pack+Sgemm

在获取了输入特征图和卷积核的Im2Col变换矩阵之后其实就可以利用Sgemm计算出卷积的结果了。

但是如果直接使用矩阵乘法计算，在卷积核尺寸比较大并且输出特征图通道数也比较大的时候，我们会发现这个时候Im2Col获得矩阵是一个行非常多列非常少的矩阵，在做矩阵乘法的时候访存会变得比较差，从而计算效率边地。这是因为当代CPU架构的某些原因，导致程序在行方向上的处理始终比列方向慢一个档次，所以我们如果能想个办法使得行方向的元素变少，列方向的元素变多这样就可以有效改善缓存达到加速的目的，另外列方向上元素增多更容易让我们发挥Neon指令集的作用。所以，这就是本文将要提到的第一个优化技巧，数据打包(Pack)。

具体来说，对于卷积核我们进行4 × 4的Pack（所谓4× 4 的Pack就是在Im2Col获得的二维矩阵的高维度进行压缩，在宽维度进行膨胀，不知道我这种说法是否合适，参考一下下面的图应该很好理解），如下图所示：

这是一个3 × 3 的卷积核并且输出通道数为4，它经过Im2Col之后首先变成上图的上半部分，然后经过Pack 4 × 4 之后变成了上图的下半部分，即变成了每个卷积核的元素按列方向交织排列。

这部分的代码也比较容易实现，另外一个技巧是，每次在执行卷积计算时，对于Image的Im2col和Pack每次都会执行，但对于卷积核，Im2col和Pack在任意次只用做一次，所以我们可以在模型初始化的时候提前把卷积核给Pack好，这样就可以节省卷积核Im2Col和Pack耗费的时间，代码实现如下：

void ConvolutionLayerSgemm::convolutionTransformKernel(float *const &kernel, const int &kernelW, const int &kernelH, float* &dest, const int &inChannel,const int &outChannel){int kernelSize = kernelH * kernelW;int ccOutChannel = 0;int ccRemainOutChannel = 0;int Stride = 0;ccOutChannel = outChannel >> 2;ccRemainOutChannel = ccOutChannel << 2;for(int cc = 0;  cc < ccOutChannel; cc ++){int c = cc << 2;const float* k0 = kernel + c * inChannel * kernelSize;const float* k1 = kernel + (c + 1) * inChannel * kernelSize;const float* k2 = kernel + (c + 2) * inChannel * kernelSize;const float* k3 = kernel + (c + 3) * inChannel * kernelSize;Stride = 4 * kernelSize * inChannel;float* destptr = dest + (c / 4) * Stride;for(int i = 0; i < inChannel * kernelSize; i++){destptr[0] = k0[0];destptr[1] = k1[0];destptr[2] = k2[0];destptr[3] = k3[0];destptr += 4;k0 += 1;k1 += 1;k2 += 1;k3 += 1;}}for(int cc = ccRemainOutChannel; cc < outChannel; cc++){int c = cc;const float* k0 = kernel + c * inChannel * kernelSize;Stride = 4 * kernelSize * inChannel;float* destptr = dest + (c / 4 + c % 4) * Stride;for(int i = 0; i < inChannel * kernelSize; i++){destptr[0] = k0[0];destptr += 1;k0 += 1;}}}

注意这个地方如果Pack的时候有拖尾部分，也就是说outChannel不能被4整除时，那么直接按照原始顺序排列即可，后面在Sgemm方法计算的时候需要注意一下。下图展示了一个5 × 3 × 3 的卷积核经过Im2Col再经过Pack4x4之后获得的结果，可以看到拖尾那一行是直接复制的，不够的长度直接置0即可。
5.2 数据的pack
接下来我们继续讲一下对于输入数据的Pack，我这里以对输入数据进行8x8的Pack为例子来讲解输入数据Pack之后应该是什么样子，以及如何和刚才已经4x4Pack好的卷积核执行矩阵乘法以完成整个卷积过程。

下面展示了一个输入维度为[ 1 , 1 , 7 , 4 ]的特征图使用Im2Col进行展开并Pack8x8之后获得结果（卷积核维度为[ 1 , 1 , 3 , 3 ] ），其中左边代表的是Im2Col后的结果，右边是将其进一步Pack8x8后获得的结果。这个地方由于outHeight×outWidth不能被8整除，所以存在拖尾部分无法进行Pack，即结果图中的第二和第三列，直接顺序放就可以了。

5.3 带Pack的矩阵乘法
所以就目前的情况来说，矩阵A就是我们的卷积核经过Im2Col+Pack4x4获得的输出矩阵，矩阵B就是我们的输入特征图经过Im2Col+Pack8x8之后获得的输出矩阵，然后矩阵C就是经过矩阵乘法之后获得的输出特征图了。在实现矩阵乘法的时候，以矩阵A为基准一次处理4行，并且在列方向分别处理4或者8个元素，这个可以结合上面的输入特征图的Im2Col+Pack8x8的示意图来想。最后，对于矩阵A不够4的拖尾部分，进行暴力处理即可。这个地方考虑的是当代典型的CNN架构一般通道数都是4 44的倍数，暂时没有实现拖尾部分的Neon优化，下面先看一下这个算法实现的整体部分，复刻NCNN。

这部分的代码实现如下：

// sgemm (int M, int N, int K, float *A, float *B, float *C)
// A (M x K)
// B (K x N)
// C (M x N)//int M = outChannel;int N = outHeight * outWidth;int K = kernelSize * inChannel;int ccOutChannel = outChannel >> 2;int ccRemainOutChannel = ccOutChannel << 2;#if USE_OMP#pragma omp parallel for num_threads(OMP_THREAD)
#endiffor(int cc = 0; cc < ccOutChannel; cc++){int c = cc << 2;float *destptr0 = dest + c * outSize;float *destptr1 = dest + (c + 1) * outSize;float *destptr2 = dest + (c + 2) * outSize;float *destptr3 = dest + (c + 3) * outSize;int i = 0;// N = outHeight*outWidthfor(; i + 7 < N; i = i+8){const float *ptrB = src_im2col_pack + (i / 8) *  packHeight * packWidth;
#if __aarch64__throw Exception(1, "Error: armv8 temporarily not supported!", __FILE__, __LINE__, __FUNCTION__);
#elseconst float *ptrA = kernel_im2col_pack + (c / 4) * kernelPackHeight * kernelPackWidth;
#endif#if USE_NEON#if __aarch64__throw Exception(1, "Error: armv8 temporarily not supported!", __FILE__, __LINE__, __FUNCTION__);
#elseasm volatile();
#endif#elsefloat sum0[8] = {0};float sum1[8] = {0};float sum2[8] = {0};float sum3[8] = {0};int j = 0;// K = kernelSize * inChannel// 同时计算4行，同时在每一列计算8个输出for(; j + 7 < K; j = j + 8){for(int n = 0; n < 8; n++){sum0[n] += ptrA[0] * ptrB[n];sum1[n] += ptrA[1] * ptrB[n];sum2[n] += ptrA[2] * ptrB[n];sum3[n] += ptrA[3] * ptrB[n];ptrA += 4;sum0[n] += ptrA[0] * ptrB[n + 8];sum1[n] += ptrA[1] * ptrB[n + 8];sum2[n] += ptrA[2] * ptrB[n + 8];sum3[n] += ptrA[3] * ptrB[n + 8];ptrA += 4;sum0[n] += ptrA[0] * ptrB[n + 16];sum1[n] += ptrA[1] * ptrB[n + 16];sum2[n] += ptrA[2] * ptrB[n + 16];sum3[n] += ptrA[3] * ptrB[n + 16];ptrA += 4;sum0[n] += ptrA[0] * ptrB[n + 24];sum1[n] += ptrA[1] * ptrB[n + 24];sum2[n] += ptrA[2] * ptrB[n + 24];sum3[n] += ptrA[3] * ptrB[n + 24];ptrA += 4;sum0[n] += ptrA[0] * ptrB[n + 32];sum1[n] += ptrA[1] * ptrB[n + 32];sum2[n] += ptrA[2] * ptrB[n + 32];sum3[n] += ptrA[3] * ptrB[n + 32];ptrA += 4;sum0[n] += ptrA[0] * ptrB[n + 40];sum1[n] += ptrA[1] * ptrB[n + 40];sum2[n] += ptrA[2] * ptrB[n + 40];sum3[n] += ptrA[3] * ptrB[n + 40];ptrA += 4;sum0[n] += ptrA[0] * ptrB[n + 48];sum1[n] += ptrA[1] * ptrB[n + 48];sum2[n] += ptrA[2] * ptrB[n + 48];sum3[n] += ptrA[3] * ptrB[n + 48];ptrA += 4;sum0[n] += ptrA[0] * ptrB[n + 56];sum1[n] += ptrA[1] * ptrB[n + 56];sum2[n] += ptrA[2] * ptrB[n + 56];sum3[n] += ptrA[3] * ptrB[n + 56];ptrA -= 28;}ptrA += 32;ptrB += 64;}// K = kernelSize * inChannel * 4// 如果是pack4x4那么末尾一定是4的倍数for(; j < K; j++){for(int n = 0; n < 8; n++){sum0[n] += ptrA[0] * ptrB[n];sum1[n] += ptrA[1] * ptrB[n];sum2[n] += ptrA[2] * ptrB[n];sum3[n] += ptrA[3] * ptrB[n];}ptrA += 4;ptrB += 8;}for(int n = 0; n < 8; n++){destptr0[n] = sum0[n];destptr1[n] = sum1[n];destptr2[n] = sum2[n];destptr3[n] = sum3[n];}#endifdestptr0 += 8;destptr1 += 8;destptr2 += 8;destptr3 += 8;}

5.4 pack矩阵乘法效果

在更大特征图上的测试结果：

有一个直观的猜测就是当输入通道数和卷积核尺寸的乘积也就是卷积核Im2Col获得的矩阵的宽度比较大时，本文介绍的Im2Col+Pack+Sgemm能带来加速效果。在输入特征图通道数较小并且特征图较大时甚至会比手工优化更慢，因为这个时候对访存的优化效果很小，因为这个时候Im2Col的矩阵行数和列数差距并不大，这个时候手工优化带来的Cache Miss并不会对卷积的计算带啦很大的影响。因此，此算法并非所有情况适用，需要灵活判断。

6. 视频

coming …

参考

https://jackwish.net/2019/gemm-optimization.html