书接上文 从 CNN 性能优化说起（一）。
上回说到贾杨清把 convolution 转化为矩阵乘。这个方法在图像领域曾经有人用过。扬清的 memo 里提到 MATLAB 里有一个操作 im2col 就是干这个的。不过在深度学习领域，扬清是第一个用这个方法加速 CNN 的。
那么矩阵乘应该怎么优化呢？这里用到的思路是通用的。在计算机专业本科的《计算机体系结构》里都有介绍的。优化数据库执行引擎、互联网后台服务、通用编译器以及深度学习编译器等系统用的也是这些路子。
CPU 提速思路
本科教材里介绍的 CPU 的设计里用来提升效能（用有限的晶体管）有如下几个：
- cache：利用程序对内存的访问”一般来说“是连续的，具有局部性的，于是预取（prefetch）被访问的内存区域之后的一段区域到 cache 里。这样接下来的内存访问访问 cache 就可以了。cache 位于 CPU 芯片内部，内存和CPU 是不同的芯片，所以访问 cache 更快。
- pipeline：把 CPU 里的晶体管分组。每条指令的执行顺序使用这些组。这样任何一组被用来执行一条指令之后，这条指令由下一组接手，而当前组不闲着，而是开始执行下一条指令。这个路子是利用大部分程序都在顺序执行一串指令（指令的局部性）。不过，如果程序里有选择分支或者循环，并不是简单地顺序执行一系列指令，这个路子就不灵了。
- vectorize：给 CPU 增加一些专门用来快速执行向量操作的指令。每条指令调动更多的晶体管，并发操作向量里的多个元素。和 pipeline 一样，vectorize 是并发思路在 CPU 里的最细粒度实现。
- multi-core：受限于制造工艺，能塞进一个 CPU 里的晶体管数量有限，所以一个 CPU 芯片里塞多个 CPU，甚至一台计算机里塞多个芯片。这里每个 CPU 叫一个 core。
接下来我们看看如何利用这些 CPU 特点加速矩阵乘。GPU 和 TPU 的体系结构虽然和 CPU 有很大区别，但是基本思路也都不超过上述范围。
利用 Cache：改变循环的嵌套顺序
我们写一个三重循环计算矩阵乘 A B = C 的时候，一般会写
for i in 0..Hc:
for j in 0..Wc:
for k in 0..Wa
C[i][j] = A[i][k] * B[k][j]
假设矩阵的内存布局是行主（row major），那么上述代码访问矩阵 A、B、C 的行和列顺序如下图。
图示来自 https://sahnimanas.github.io/post/anatomy-of-a-high-performance-convolution
请注意，这里对 A 和 C 的访问都可以利用 cache，因为内存布局是 row major，而对 A 和 C 都是一行一行地访问的。但是对 B 的访问是一列一列的。CPU prefetch B 的行，但是往往用不上。实际计算的时候还是经常得去访问内存。这样就慢了。
解决之道很简单。把上面三重循环的顺序调换一下：
for i in 0..Hc:
for k in 0..Wa
for j in 0..Wc:
C[i][j] = A[i][k] * B[k][j]
这样，对三个矩阵的访问顺序就都变成了一行一行的：
图示来自 https://sahnimanas.github.io/post/anatomy-of-a-high-performance-convolution
这样能快多少呢？根据矩阵大小不同，提速分别是几倍到几十十倍。
图示来自 https://sahnimanas.github.io/post/anatomy-of-a-high-performance-convolution
利用 Pipline：展开循环
在上面三重循环里的 step 只有一个相乘的指令，没有一串多条指令，也就没法利用 CPU 的 pipeline。一个直接的解决之道是展开循环 —— 把最内层循环展开，也就是手写成多条语句，以便编译器翻译成一串指令。
for i in 0..Hc:
for k in 0..Wa
C[i][0] = A[i][k] * B[k][0]
C[i][1] = A[i][k] * B[k][1]
C[i][2] = A[i][k] * B[k][2]
...
C[i][Wc-1] = A[i][k] * B[k][Wc-1]
这个想法不错，但是并不可行，因为展开的重复次数 Wc 是运行时才知道的，编译前写代码的时候并不能预知。
一个解决之道是 tiling —— 把 Wc 拆成某个常数 s（比如s=8）的倍数 —— Wc 是几个 s 的乘积虽然在编译时未知，但是 s 这个常数是已知的，就可以展开了。
矩阵乘问题的 tiling 思路是：结果 C 中的一个 r x s 的小方块需要用到 A 中的 r 行和 B 中的 s 列。这样就可以把 A 拆成若干个 r 行，把 B 拆成若干个 s 列。用 xt 和 yt 来索引 C 中的 tiles，用 x 和 y 来索引 tile 中的元素，则三重循环可以被重写如下：
for yt in 0..Hc/r:
for xt in 0..Wc/s:
for y in 0..r:
for k in 0..Wa:
for x in 0..s:
C[yt*r+y][xt*s+x] = A[yt*r+y][k] * B[k][xt*s+x]
这个最内层循环可以被展开。
for yt in 0..Hc/r:
for xt in 0..Wc/s:
for y in 0..r:
for k in 0..Wa:
C[yt*r+y][xt*s+0] = A[yt*r+y][k] * B[k][xt*s+0]
C[yt*r+y][xt*s+1] = A[yt*r+y][k] * B[k][xt*s+1]
...
C[yt*r+y][xt*s+s-1] = A[yt*r+y][k] * B[k][xt*s+s-1]
在原文 https://sahnimanas.github.io/post/anatomy-of-a-high-performance-convolution 里，作者把循环展开（unrolling）带来的性能提升解释为减少了循环 for x in 0..s 中对 x 和 s 的比较。其实如果只是减少这个比较，是不足以带来这么大的性能提升的。主要还是利用了 CPU 的 pipeline。
另外，需要强调一下，只做 tiling（不做循环展开）也是可以带来性能提升的，因为更好地利用了 cache。
利用 Vectorized 指令
Intel CPU 中有一类被称为 SIMD 的指令集，包括 SSE、SSE2、AVX，其中的每条指令可以做 8 个浮点数的算术运算。上文中的循环展开使得编译器更有机会发现可以被 SIMD 指令优化的代码。
利用多线程和Multicore
循环展开的另一种形式是多线程：不是依次执行 steps，而是并发执行。简单的说，把上面代码中某一层（相对外层）循环中的 for 改成 parallel for 即可。对于 CPU 优化，可以用 OpenMP Home 。
所有招数一起上
在原文 https://sahnimanas.github.io/post/anatomy-of-a-high-performance-convolution 最后，作者把上述招数都用上后发现原来的三重循环被优化了 100 多倍，达到了 OpenBLAS 的性能，接近了他使用的 CPU 的理论最高峰值计算能力。
Finally, we’re able to touch speeds of over 120 GFLOPs – respectably close to the peak performance of 160 GFLOPs, and matching production-level libraries like OpenBLAS.
小结
这些技法和TVM，XLA有什么关系呢？请听下回分解。