詳解卷積中的Winograd加速演算法
6. 何時開啟WinoGrad卷積
和Sgemm用於卷積一樣，我們也需要思考WinoGrad在何種情況下是適用的，或者說是有明顯加速的。這篇文章介紹的WinoGrad卷積是針對NCHW這種記憶體排布的，然後我們來看一下NCNN在基於NCHW這種記憶體排佈下，是在何種情況下啟用WinoGrad()？
通過檢視NCNN的原始碼(https://github.com/Tencent/ncnn/blob/master/src/layer/arm/convolution_arm.cpp)可以發現，只有在輸入輸出通道均>=16，並且特徵圖長寬均小於等於120的條件下才會啟用WinoGrad卷積。
那麼這個條件是如何得出的，除了和手工優化的conv3x3s1
https://github.com/msnh2012/Msnhnet/blob/master/src/layers/arm/MsnhConvolution3x3s1.cpp)在不同條件下做速度對比測試之外，我們也可以感性的分析一下。
第一，WinoGrad演算法設計到幾個矩陣變換，如果計算量不大，這幾個矩陣變換的成本佔計算總成本的比例就越大，所以WinoGrad應當是在計算量比較大時才能有效，如VGG16。
第二，當計算量比較大的時候，又要考慮到Cache命中率的問題，這個時候WinoGrad訪存可能會比直接手動優化更差，導致速度上不去。
詳解卷積中的Winograd加速演算法
1天前閱讀400「GiantPandaCV導語」：這篇文章為大家介紹一下用來加速卷積運算的WinoGrad演算法的原理，工程實現以及相關優化思路，如果你對卷積加速演算法感興趣可以看看這篇文章。演算法的完整實現請到MsnhNet的github倉庫檢視，地址為：https://github.com/msnh2012/Msnhnet
1. 為什麼會引入WinoGrad？
做過ACM/OI的朋友大家應該對FFT並不陌生，我們知道對於兩個序列的乘法通過FFT可以從原始O(n^2)複雜度變成O(nlogn)，所以我們就會想著FFT這個演算法是否可以應用到我們計算卷積中來呢？當然是可以的，但是FFT的計算有個問題哦，會引入複數。而移動端是不好處理複數的，對於小卷積核可能減少的計算量和複數運算帶來的降速效果是不好說誰會主導的。所以在這種情況下，針對卷積的WinoGrad演算法出現了，它不僅可以類似FFT一樣降低計算量，它還不會引入複數，使得卷積的運算加速成為了可能。因此，本文嘗試從工程實現的角度來看一下WinoGrad，希望對從事演算法加速的小夥伴有一些幫助。
2. 為什麼會有這篇文章？
最近嘗試給MsnhNet做卷積的WinoGrad實現，然後開始瞭解這個演算法，並嘗試參考著NCNN來理解和動手寫一下。參考了多篇優秀的講解文章和NCNN原始碼，感覺算是對這個演算法有了較為清楚的認識，這篇文章就記錄一下我在實現並且步長為的WinoGrad卷積時的一些理解。這篇文章的重點是WinoGrad卷積的實現，關於WinoGrad卷積裡面的變化矩陣如何推導可以看樑德澎作者的文章：詳解Winograd變換矩陣生成原理 （聽說後續他會做個視訊來仔細講講QAQ)，現在就假設我們知道了WinoGrad的幾個變換矩陣。如果你不知道也沒關係，因為有一個Python工具包可以直接幫我們計算，地址為：https://github.com/andravin/wincnn 。然後現在我們就要用拿到的這幾個矩陣來實現WinoGrad演算法，聽起來比較簡單，但我們還是得一步步理清楚是不。
3. WinoGrad演算法原理
WinoGrad演算法起源於1980年，是Shmuel Winograd提出用來減少FIR濾波器計算量的一個演算法。它指出，對於輸出個數為，引數個數為的FIR濾波器，不需要次乘法計算，而只需要次乘法計算即可。
下面是一個經典例子，以1維卷積為例，輸入訊號，卷積核，則卷積可以寫成如下矩陣乘法形式：
式子1如果這個計算過程使用普通的矩陣乘法，則一共需要「6次乘法和4次加法」 。
但是，我們仔細觀察一下，卷積運算中輸入訊號轉換得到的矩陣不是任意矩陣，其有規律的分佈著大量的重複元素，例如第一行的和，卷積轉換成的矩陣乘法比一般乘法的問題域更小，所以這就讓優化存為了可能。
然後WinoGrad的做法就是：
式子2其中，
等式3我們知道，在CNN的推理階段，卷積核上的元素是固定的，所以上式中和相關的式子可以提前算好，在預測階段只用計算一次，可以忽略。所以這裡一共需要「4次乘法加4次加法」。
相比於普通的矩陣乘法，使用WinoGrad演算法之後乘法次數減少了，這樣就可以達到加速的目的了。
這個例子實際上是「1D的WinoGrad演算法」，我們將上面的計算過程寫成矩陣的形式如下：
式子4其中，表示element-wise multiplication（Hadamard product）對應位置相乘。其中，
相關矩陣解釋- ：表示卷積核
- ：表示輸入訊號
- ：卷積核變換矩陣，尺寸為
- ：輸入變換矩陣，尺寸
- ：輸出變換矩陣，尺寸
所以整個計算過程可以分為4步：
- 輸入變換
- 卷積核變換
- 外積
- 輸出變換
然後我們將1D的WinoGrad擴充套件到2D，就可以實現卷積的加速了，那麼如何從1維擴充套件到2維呢？公式如下：
式子5其中，為的卷積核，為的影象塊，我們把上面的擴充套件到，先寫成矩陣乘法的方式：
F(2x2,3x3) 圖片來自https://www.cnblogs.com/shine-lee/p/10906535.html上圖表示我們將卷積核的元素拉成了一列，將輸入訊號每個滑動視窗中的元素拉成了一行。注意圖中紅線分成的矩陣塊，每個矩陣塊中重複元素的位置與一維相同，即：
二維和一維的WinoGrad矩陣關係然後，令，即影象視窗中的第0行元素，然後表示第行，，然後可以推導：
2D WinoGrad矩陣形式計算推導在上面的推導中，表示長度為4的和長度為的卷積結果，結果為長度為2的列向量，其中和均為長度為4的列向量。
進一步，可以看成3對長度為4的列向量兩兩對應位置相乘再相加，結果為長度為4的列向量，也可以看成是4組長度為3的行向量的點積運算。
同樣，也是3對長度為4的列向量的內積運算。
然後類似1D WinoGrad演算法，我們考慮兩者的重疊部分和，剛好對應1D WinoGrad中的每一行在的對應行上進行1維卷積，基於上面推導的1D WinoGrad公式，行向量的卷積只需要將所有左乘的變換矩陣轉置後變成右乘即可。
然後上面的推導就做完了。
下圖表示2D WinoGrad的示意圖：
2D WinoGrad示意圖這個時候，WinoGrad演算法的乘法次數為，而如果直接卷積乘法次數為，「降低了2.25倍的乘法計算複雜度」。
4. 從工程角度來看WinoGrad
下面我們就從一個實際例子來說，如何利用WinoGrad來實現並且步長為1的卷積運算。基於上面介紹的2D WinoGrad的原理，我們現在只需要分4步即可實現WnoGrad演算法：
- 第一步就是對輸入卷積核的變換：
- 第二步就是對輸入資料的變換：
- 第三步就是對M矩陣的計算：
- 最後一步就是結果的計算：
接下來我們就以WinoGrad實現並且步長為1的卷積計算為例子，來理解一下WinoGrad的工程實現。
4.1 對輸入卷積核進行變換
這一步就是對卷積核進行變化，公式為：，其中表示輸出通道標號，表示輸入通道標號，一個對應卷積核的一個。由於我們要實現的是，因此是一個的矩陣，我們不難寫出這部分程式碼（其中，矩陣可以通過https://github.com/andravin/wincnn 這個工具進行計算）：
// 矩陣G
const float ktm[8][3] = {
{1.0f, 0.0f, 0.0f},
{-2.0f / 9, -2.0f / 9, -2.0f / 9},
{-2.0f / 9, 2.0f / 9, -2.0f / 9},
{1.0f / 90, 1.0f / 45, 2.0f / 45},
{1.0f / 90, -1.0f / 45, 2.0f / 45},
{1.0f / 45, 1.0f / 90, 1.0f / 180},
{1.0f / 45, -1.0f / 90, 1.0f / 180},
{0.0f, 0.0f, 1.0f}
};
const int kernelTmSize = inChannel * 8 * 8;
#if USE_OMP
#pragma omp parallel for num_threads(OMP_THREAD)
#endif
for(int outc = 0; outc < outChannel; outc++){
for(int inc = 0; inc < inChannel; inc++){
const float* kernel0 = (const float*)kernel + outc * inChannel * 9 + inc * 9;
float *kernel_tm0 = kernel_tm + outc * kernelTmSize + inc * 64;
//需要變換的卷積核
const float* k0 = kernel0;
const float* k1 = kernel0 + 3;
const float* k2 = kernel0 + 6;
float tmpG[8][3]; // tmp = G*g
for(int i = 0; i < 8; i++){
tmpG[i][0] = k0[0] * ktm[i]