GEMM矩阵相乘与深度学习

1. GEMM矩阵相乘

1.1 GEMM算法基础

  GEMM（General Matrix Multiplication），通用矩阵乘法。

则根据矩阵相乘的定义，可得C=AB+C的计算公式为：

从而可实现为：

for i=0:m-1
    for j=0:n-1
        for p=0:k-1
            C(i,j):=A(i,p)*B(p,j)+C(i,j)
        endfor
    endfor
endfor

  可以看出，根据定义实现的算法计算复杂度为O(mnk)，也就是O(n^3)。

1.2 GEMM算法优化

1.2.1 循环重排充分利用缓存

  应该首先注意到我们访问数据的模式。我们按照下图A的形式逐行遍历数据，按照下图B的形式逐列遍历数据。

  它们的存储也是行优先的，因此一旦我们找到 A[i, k]，则它在该行中的下一个元素A[i, k+1]已经被缓存了。接下来我们来看B中发生了什么：

  列的下一个元素并未出现在缓存中，即出现了缓存缺失（cache miss）。这时尽管获取到了数据，CPU也出现了一次停顿。获取数据后，缓存同时也被 B 中同一行的其他元素填满。我们实际上并不会使用到它们，因此它们很快就会被删除。多次迭代后，当我们需要那些元素时，我们将再次获取它们。我们在用实际上不需要的值污染缓存。

  我们需要重新修改loop，以充分利用缓存能力。如果数据被读取，则我们要使用它。这就是我们所做的第一项改变：循环重排序（loop reordering）。

  我们将i,j,k 循环重新排序为 i,k,j：

for i in0..M: 
	for k in0..K: 
		for j in0..N:

  答案仍然是正确的，因为乘/加的顺序对结果没有影响。而遍历顺序则变成了如下状态：

  循环重排序这一简单的变化，却带来了相当可观的加速。

1.2.2 平铺（Tiling）充分利用缓存

  要想进一步改进重排序，我们需要考虑另一个缓存问题。

  对于A中的每一行，我们针对B中所有列进行循环。而对于 B 中的每一步，我们会在缓存中加载一些新的列，去除一些旧的列。当到达A的下一行时，我们仍旧重新从第一列开始循环。我们不断在缓存中添加和删除同样的数据，即缓存颠簸（cache thrashing）。

  如果所有数据均适应缓存，则颠簸不会出现。如果我们处理的是小型矩阵，则它们会舒适地待在缓存里，不用经历重复的驱逐。庆幸的是，我们可以将矩阵相乘分解为子矩阵。要想计算 C 的r×c平铺，我们仅需要A的r行和B的c列。接下来，我们将 C 分解为6x16的平铺：

C(x, y) += A(k, y) * B(x, k);
C.update().tile(x, y, xo, yo, xi, yi, 6, 16)/*in pseudocode:
for xo in 0..N/16:
	 for yo in 0..M/6: 
		for yi in 6: for xi in 0..16:
			 for k in 0..K: C(...) = ...*/

  我们将x,y 维度分解为外侧的xo,yo和内侧的xi,yi。我们将为该6x16 block优化micro-kernel（即xi,yi），然后在所有block上运行micro-kernel（通过xo,yo进行迭代）。

1.2.3 展开（Unrolling）

  循环使我们避免重复写同样代码的痛苦，但同时它也引入了一些额外的工作，如检查循环终止、更新循环计数器、指针运算等。如果手动写出重复的循环语句并展开循环，我们就可以减少这一开销。例如，不对1个语句执行8次迭代，而是对4个语句执行2次迭代。

  这种看似微不足道的开销实际上是很重要的，最初意识到这一点时我很惊讶。尽管这些循环操作可能「成本低廉」，但它们肯定不是免费的。每次迭代2-3个额外指令的成本会很快累积起来，因为此处的迭代次数是数百万。随着循环开销越来越小，这种优势也在不断减小。

1.2.4 内存对齐

内存对齐的原则：任何K字节的基本对象的地址必须都是K的倍数。

  假设 cache line 为 32B。待访问数据大小为 64B，地址在 0x80000001，则需要占用 3 条 cache 映射表项；若地址在 0x80000000 则只需要 2 条。内存对齐变相地提高了 cache 命中率。

1.2.5 向量化

  大部分现代CPU支持SIMD（Single Instruction Multiple Data，单指令流多数据流）。在同一个CPU循环中，SIMD可在多个值上同时执行相同的运算/指令（如加、乘等）。如果我们在4个数据点上同时运行SIMD指令，就会直接实现4倍的加速。

  使用FMA（Fused Multiply-Add）。尽管乘和加是两种独立的浮点运算，但它们非常常见，有些专用硬件单元可以将二者融合为一，作为单个指令来执行。编译器通常会管理FMA的使用。

1.2.6 矩阵分块

  充分利用多核性能。

  GEMM的分块方式：

　　图中共有6中拆分方法，依次分析：

1. 拆分K，A矩阵拆成多列，B矩阵拆成多行。拆分M，A的一列拆分成三个block。拆分N，B的一行也可以拆分成多个小slice，每个slice放到寄存器中。遍历A的列，A的每个block乘以B的一行，得到矩阵C的一行的部分值。 拆分之后Aip为(mc, kc), Bpj为(kc, nr), Cij为(mc, nr)。
   
   
   如果kc*mc很小，那么Aip、Bpj、Cij都放入cache中，Aip（A矩阵的一个block）只需要被加载进cache一次，提高了cache命中率。对Ai和Bj进行pack使其内存连续。由于处理B矩阵是按照每个slice依次进行，所以这种划分更适合于列存储优先的矩阵乘。

2. 拆分K，A拆分成多列，B拆成多行。拆分N，B的一行被拆分为三个block。拆分M，A的一列被拆分为多个slice。A的一列乘以B的每个block，得到矩阵C的一列的部分值。
   
   

   类似于(1), 只是变成A的一行乘以B的一个block。所以这种划分更适合于行存储优先的矩阵乘。

3. 拆分N，B矩阵被拆分为多列。再拆分K，A拆分成多列，B拆成多个block。拆分M，A的一列被拆分为多个小block。A的一列乘以B的每个block，得到矩阵C的一列的部分值。
   
   

   这种划分与(2)的划分唯一的区别是，访问B矩阵是按列还是按行。很显然(2)的划分方式好于(3)。

4. 类似于普通矩阵乘，A的一行*B的一列。然后在K拆分，将A矩阵的每行划分为多个block, 将B矩阵的每列划分为多个block。
   
   　　　　　　　　　　　　　 　
   每次执行可以得到C矩阵一个block的值。当MNK非常大时，cache无法存下Cj，所以划分方法没有什么优势。

5. 类似于(4), 都是每次执行可以得到C矩阵的一个block的值。同样当MNK非常大时，cache无妨存下Cj，所以该方法没有什么优势。
   

6. 类似于（1）。不一样在于：对A矩阵的block遍历是按行访问。和（1）相比，最外层循环是i，而非p。遍历顺序为A0B0 + A1B1 + … + ApBp 才能得到C一行的最终值，因此需要不断的更新C这一行的值，可以使用一个顺序的cache数组存放 pack ，最后在unpack的时候累加。这与(1)中的按列访问A矩阵的block不同,（1)中Ci += Ai, 在最内层循环时只需要使用寄存器存储每次Ai*Bpj的结果。
   
   

   由于在循环外 cache 级别 unpack C 远复杂于在循环内 register 级别 unpack C，因此(1)比(6)好。

1.2.7 双缓冲

  计算的同时在不同级的存储介质中进行数据搬移，隐藏访存开销


将矩阵分块同时利用不同核进行计算，充分利用多核性能

2. GEMM与卷积计算

2.1 计算卷积的方法

  计算卷积的方法有很多种，常见的有以下几种方法：

• 滑窗：这种方法是最直观最简单的方法。 但是，该方法不容易实现大规模加速，因此，通常情况下不采用这种方法（但是也不是绝对不会用，在一些特定的条件下该方法反而是最高效的）。
• im2col：目前几乎所有的主流计算框架包括 Caffe, MXNet 等都实现了该方法. 该方法把整个卷积过程转化成了GEMM过程，而GEMM在各种 BLAS 库中都是被极致优化的，一般来说，速度较快。
• FFT：傅里叶变换和快速傅里叶变化是在经典图像处理里面经常使用的计算方法，但是，在 ConvNet中通常不采用，主要是因为在 ConvNet 中的卷积模板通常都比较小，例如 3×3 等，这种情况下，FFT 的时间开销反而更大，所以很少在CNN中利用FFT实现卷积。
• Winograd： Winograd 是存在已久最近被重新发现的方法，在很多场景中, Winograd方法都显示出较大的优势，目前cudnn中计算卷积就使用了该方法。

2.2 Img2col

2.2.1 CNN中张量的存储

  通常张量（如CNN中的张量）的存储顺序是NCHW、NHWC等。

  NCHW存储顺序，即如果HxW 图像的block为N，通道数为C，则具备相同N的所有图像是连续的，同一block内通道数C相同的所有像素是连续的

2.2.2 卷积运算转化为GEMM

  卷积是滤波器和输入图像块（patch）的点乘。如果我们将滤波器展开为2-D矩阵，将输入块展开为另一个2-D矩阵，则将两个矩阵相乘可以得到同样的数字。与CNN不同，近几十年来矩阵相乘已经得到广泛研究和优化，成为多个科学领域中的重要问题。将图像块展开为矩阵的过程叫做im2col（image to column）。我们将图像重新排列为矩阵的列，每个列对应一个输入块，卷积滤波器就应用于这些输入块上。

  下图展示了一个正常的3x3卷积：

  下图展示的是该卷积运算被实现为矩阵相乘的形式。右侧的矩阵是im2col的结果，它需要从原始图像中复制像素才能得以构建。左侧的矩阵是卷积的权重，它们已经以这种形式储存在内存中。

  出于视觉简洁考虑，此处将每个图像块作为独立的个体进行展示。而在现实中，不同图像块之间通常会有重叠，因而im2col可能导致内存重叠。生成im2col 缓冲（im2col buffer）和过多内存（inflated memory）所花费的时间必须通过GEMM实现的加速来抵消。

  使用im2col可以将卷积运算转换为矩阵相乘。现在我们可以使用更加通用和流行的线性代数库（如OpenBLAS、Eigen）对矩阵相乘执行高效计算，而这是基于几十年的优化和微调而实现的。