書接上文 從 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有什麼關系呢？請聽下回分解。