图解ncnn实现Winograd卷积

心已赠人 2023-01-17 07:42 155阅读 0赞

# 图解ncnn实现Winograd卷积 #

## 前言 ##

[上一篇][Link 1]讲了Winograd的基本实现思路, 提到ncnn的实现是**使用结论公式描述计算的流程, 用固定已知系数的加法代替矩阵乘法**, 除此之外, ncnn的实现还涉及到一系列的内存重排配合SIMD操作提高效率. 本文用多幅图描述ncnn实现Winograd的详细流程. 希望能够帮助大家更好的理解.

## 流程概述 ##

Y = A T \[ ( G g G T ) ⊙ ( B T d B ) \] A = A T \[ U ⊙ V \] A = A T M A (1) Y=A^T\[(GgG^T)⊙(B^TdB)\]A =A^T\[U⊙V\]A=A^TMA \\tag\{1\} Y=AT\[(GgGT)⊙(BTdB)\]A=AT\[U⊙V\]A=ATMA(1)

回顾一下公式1, 其中 g g g代表卷积核,  d d d代表卷积输入,  Y Y Y代表卷积输出, 其余都为已知矩阵. 二维Winograd卷积的实现分为4步, 分别是:

*  kernel\_transform对应公式中的 U = G g G T U=GgG^T U=GgGT
 *  input\_transform对应公式中的 V = B T d B V=B^TdB V=BTdB
 *  multi对应公式中的 M = U ⊙ V M=U⊙V M=U⊙V
 *  ouput\_transform对应公式中的 Y = A T M A Y=A^TMA Y=ATMA

本文分析的是ncnn的Winograd实现函数包括[conv3x3s1\_winograd64\_transform\_kernel\_neon5][conv3x3s1_winograd64_transform_kernel_neon5]和[conv3x3s1\_winograd64\_neon5][conv3x3s1_winograd64_transform_kernel_neon5], 对应的Winograd的 F ( 6 , 3 ) F(6,3) F(6,3), 即每个tile块的输出尺寸是\[6x6\], 卷积核尺寸是\[3x3\], 输入尺寸是\[8x8\].

## 步骤一: kernel\_transform ##

### step 1 ###

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center]

图1.kernel\_transform step1

对于一个nchw排布的kernel, 先根据公式 U = G g G T U=GgG^T U=GgGT对3x3的卷积核 g g g进行了转换得到8x8的尺寸. 其实这个过程是可以推理之前做好存储起来的, 所以即使用最慢的矩阵乘法也是可以的.  
 G = \[ 1 0 0 − 2 9 − 2 9 − 2 9 − 2 9 2 9 − 2 9 1 90 1 45 2 45 1 90 − 1 45 2 45 32 45 16 45 8 45 32 45 − 16 45 8 45 0 0 1 \] G T = \[ 1 − 2 9 − 2 9 1 90 1 90 32 45 32 45 0 0 − 2 9 2 9 1 45 − 1 45 16 45 16 45 0 0 − 2 9 − 2 9 2 45 2 45 8 45 8 45 1 \] (2) G = \\left\[ \\begin\{matrix\} 1 & 0 & 0\\\\ -\\frac\{2\}\{9\} & -\\frac\{2\}\{9\} & -\\frac\{2\}\{9\} \\\\ -\\frac\{2\}\{9\} & \\frac\{2\}\{9\} & -\\frac\{2\}\{9\} \\\\ \\frac\{1\}\{90\} & \\frac\{1\}\{45\} & \\frac\{2\}\{45\} \\\\ \\frac\{1\}\{90\} & -\\frac\{1\}\{45\} & \\frac\{2\}\{45\} \\\\ \\frac\{32\}\{45\} & \\frac\{16\}\{45\} & \\frac\{8\}\{45\} \\\\ \\frac\{32\}\{45\} & -\\frac\{16\}\{45\} & \\frac\{8\}\{45\} \\\\ 0 & 0 & 1\\\\ \\end\{matrix\} \\right\] \\qquad G^T = \\left\[ \\begin\{matrix\} 1 & -\\frac\{2\}\{9\} & -\\frac\{2\}\{9\} & \\frac\{1\}\{90\} & \\frac\{1\}\{90\} & \\frac\{32\}\{45\} & \\frac\{32\}\{45\} & 0\\\\ 0 & -\\frac\{2\}\{9\} & \\frac\{2\}\{9\} & \\frac\{1\}\{45\} & -\\frac\{1\}\{45\} & \\frac\{16\}\{45\} & \\frac\{16\}\{45\} & 0\\\\ 0 & -\\frac\{2\}\{9\} & -\\frac\{2\}\{9\} & \\frac\{2\}\{45\} & \\frac\{2\}\{45\} & \\frac\{8\}\{45\} & \\frac\{8\}\{45\} & 1\\\\ \\end\{matrix\} \\right\] \\tag\{2\} G=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡1−92−929019014532453200−9292451−4514516−451600−92−924524524584581⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤GT=⎣⎡100−92−92−92−9292−92901451452901−4514524532451645845324516458001⎦⎤(2)  
由于这部分ncnn的实现也是很简单粗暴的矩阵乘法, 因此就不做过多解释

static void conv3x3s1_winograd64_transform_kernel_neon5(const Mat& kernel, Mat& kernel_tm, int inch, int outch){ 
        //......
        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}
        };
        //......
        float tmp[8][3];
        for (int i = 0; i < 8; i++) { 
            tmp[i][0] = k0[0] * ktm[i][0] + k0[1] * ktm[i][1] + k0[2] * ktm[i][2];
            tmp[i][1] = k1[0] * ktm[i][0] + k1[1] * ktm[i][1] + k1[2] * ktm[i][2];
            tmp[i][2] = k2[0] * ktm[i][0] + k2[1] * ktm[i][1] + k2[2] * ktm[i][2];
        }
        //......
        for (int j = 0; j < 8; j++){ 
            float* tmpp = &tmp[j][0];
            for (int i = 0; i < 8; i++) { 
                kernel_tm0[j * 8 + i] = tmpp[0] * ktm[i][0] + tmpp[1] * ktm[i][1] + tmpp[2] * ktm[i][2];
            }
        }
        //......
    }

### step2 ###

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 1]

图2.kernel\_transform step2

这一步就是reshape的过程, 其实就是将原本的hw两个维度合并成一个w维度, 不改变元素在内存中的排列顺序. 看图应该就很直观

### step3 ###

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 2]

图3.kernel\_transform step3

这一步是一个reorder的过程, 是改变元素排列顺序的. 这个过程用文字描述就是: reorder前沿着n这个维度每4个为一组, 横着放到reorder后的一行上, 这一行由多个c的这样的组构成. 每一行代表\[8x8\]tile块的一个元素,

## 步骤二: input\_transform ##

### step1 ###

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 3]

图4.input\_transform step1

这一步包含两个操作, 第一个操作是将输入bottom\_blob划分成多个tile, 第二个操作是对每个tile根据公式 V = B T d B V=B^TdB V=BTdB完成transform. 由于transform前后每个tile都是8x8的尺寸, 所以图4中就将这两个操作合在一起了. 划分tile不多说, 重点讲讲这个transform的过程.

首先看计算 d B dB dB的矩阵乘

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 4]

我们观察矩阵 B B B, 根据其数值特点将其划分为4块

*  红色列有两个数值特点:
    
    1.  很多0和1
    2.  列内系数成对出现(不考虑正负号)
    
    分布在左右的红色两列需要分别计算. 在计算红色第一列时根据特点1可以利用[上一篇][Link 1]提到思路三 **“很多0和正负1, 其实我是可以不用乘的. 固定死写成加法即可”** 减少乘法操作数量, 同时根据特点2列内系数成对出现, 可以 **“乘法合并”, 先加减后乘除** , 进一步降低乘法操作数量.  
     ( d B ) 0 = 1 ∗ d 0 + 0 ∗ d 1 − 5.25 ∗ d 2 + 0 ∗ d 3 + 5.25 ∗ d 4 + 0 ∗ d 5 − d 6 ∗ 1 + 0 ∗ d 7 = d 0 − 5.25 ∗ d 2 + 5.25 ∗ d 4 − d 6 = d 0 − d 6 + 5.25 ∗ ( d 4 − d 2 ) (4) \\begin\{aligned\} (dB)\_0 &= 1\*d\_0+0\*d\_1-5.25\*d\_2+0\*d\_3+5.25\*d\_4+0\*d\_5-d\_6\*1+0\*d\_7\\\\ &=d\_0-5.25\*d\_2+5.25\*d\_4-d\_6\\\\ &=d\_0-d6+5.25\*(d4-d2) \\end\{aligned\} \\tag\{4\} (dB)0=1∗d0\+0∗d1−5.25∗d2\+0∗d3\+5.25∗d4\+0∗d5−d6∗1\+0∗d7=d0−5.25∗d2\+5.25∗d4−d6=d0−d6\+5.25∗(d4−d2)(4)
 *  绿色列的特点除了红色列具有的上述两个特点外还有第三个特点:
    
    1.  两列之间的差异只是部分符号相反
    
     ( d B ) 1 = d 1 + d 2 − 4.25 ∗ d 3 − 4.25 ∗ d 4 + d 5 + d 6 = ( d 2 + d 6 − 4.25 ∗ d 4 ) + ( d 1 + d 5 − 4.25 ∗ d 3 ) = a + b (5) \\begin\{aligned\} (dB)\_1 &= d\_1+d\_2-4.25\*d\_3-4.25\*d\_4+d\_5+d\_6\\\\ &=(d\_2+d\_6-4.25\*d\_4)+(d\_1+d\_5-4.25\*d\_3)\\\\ &=a+b\\\\ \\end\{aligned\} \\tag\{5\} (dB)1=d1\+d2−4.25∗d3−4.25∗d4\+d5\+d6=(d2\+d6−4.25∗d4)\+(d1\+d5−4.25∗d3)=a\+b(5)
    
     ( d B ) 2 = − d 1 + d 2 + 4.25 ∗ d 3 − 4.25 ∗ d 4 − d 5 + d 6 = ( d 2 + d 6 − 4.25 ∗ d 4 ) − ( d 1 + d 5 − 4.25 ∗ d 3 ) = a − b (6) \\begin\{aligned\} (dB)\_2 &= -d\_1+d\_2+4.25\*d\_3-4.25\*d\_4-d\_5+d\_6\\\\ &=(d\_2+d\_6-4.25\*d\_4)-(d\_1+d\_5-4.25\*d\_3)\\\\ &=a-b\\\\ \\end\{aligned\} \\tag\{6\} (dB)2=−d1\+d2\+4.25∗d3−4.25∗d4−d5\+d6=(d2\+d6−4.25∗d4)−(d1\+d5−4.25∗d3)=a−b(6)
    
    你看公式5和公式6所表示的计算绿色第一列和第二列, 是可以提炼出公共部分a和b, 因此又可以合并计算降低乘法操作数量. 蓝色列与黄色列同理.

再看计算 B T ( d B ) B^T(dB) BT(dB)的矩阵乘, 同理不赘述

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 5]

此时再看代码, 就一目了然了

static void conv3x3s1_winograd64_neon5(const Mat& bottom_blob, Mat& top_blob, const Mat& kernel_tm, const Mat& _bias, const Option& opt){
        //...
        //计算dB
        for (int m = 0; m < 8; m++) {
            // 红色两列
            tmp[0][m] = r0[0] - r0[6] + (r0[4] - r0[2]) * 5.25f;
            tmp[7][m] = r0[7] - r0[1] + (r0[3] - r0[5]) * 5.25f;
    
            // 绿色两列
            float tmp12a = (r0[2] + r0[6] - r0[4] * 4.25f);
            float tmp12b = (r0[1] + r0[5] - r0[3] * 4.25f);
            tmp[1][m] = tmp12a + tmp12b;
            tmp[2][m] = tmp12a - tmp12b;
    
            // 蓝色两列
            float tmp34a = (r0[6] + r0[2] * 0.25f - r0[4] * 1.25f);
            float tmp34b = (r0[1] * 0.5f - r0[3] * 2.5f + r0[5] * 2.f);
            tmp[3][m] = tmp34a + tmp34b;
            tmp[4][m] = tmp34a - tmp34b;
    
            // 黄色两列
            float tmp56a = (r0[6] + (r0[2] - r0[4] * 1.25f) * 4.f);
            float tmp56b = (r0[1] * 2.f - r0[3] * 2.5f + r0[5] * 0.5f);
            tmp[5][m] = tmp56a + tmp56b;
            tmp[6][m] = tmp56a - tmp56b;
    
            r0 += w;
        }
        //...
        //计算BT(dB)
        for (int m = 0; m < 8; m++) {
            const float* tmp0 = tmp[m];
    		// 红色两列
            r0_tm_0[0] = tmp0[0] - tmp0[6] + (tmp0[4] - tmp0[2]) * 5.25f;
            r0_tm_7[0] = tmp0[7] - tmp0[1] + (tmp0[3] - tmp0[5]) * 5.25f;
    		
            // 绿色两列
            float tmp12a = (tmp0[2] + tmp0[6] - tmp0[4] * 4.25f);
            float tmp12b = (tmp0[1] - tmp0[3] * 4.25f + tmp0[5]);
            r0_tm_1[0] = tmp12a + tmp12b;
            r0_tm_2[0] = tmp12a - tmp12b;
    
            // 蓝色两列
            float tmp34a = (tmp0[6] + tmp0[2] * 0.25f - tmp0[4] * 1.25f);
            float tmp34b = (tmp0[1] * 0.5f - tmp0[3] * 2.5f + tmp0[5] * 2.f);
            r0_tm_3[0] = tmp34a + tmp34b;
            r0_tm_4[0] = tmp34a - tmp34b;
    		
            // 黄色两列
            float tmp56a = (tmp0[6] + (tmp0[2] - tmp0[4] * 1.25f) * 4.f);
            float tmp56b = (tmp0[1] * 2.f - tmp0[3] * 2.5f + tmp0[5] * 0.5f);
            r0_tm_5[0] = tmp56a + tmp56b;
            r0_tm_6[0] = tmp56a - tmp56b;
            //...
        }
        //...
    }

### step2 ###

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 6]

图5.input\_transform step2

第二步是一个内存重排的操作. 把每个tile的同一位置的元素连在一起形成一个块, 本文称其为"tiles块"(注意区分前面说的"tile块"). 同时tiles也代表tile的数量, 图5中等于4.

### step3 ###

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 7]

图6.input\_transform step3

第三步又是一个内存重排. 首先你可能会疑惑, 为什么第二步做了内存重排, 紧接着第三步又做了内存重排, 不能一步到位直接排成想要的样子么? 其实ncnn的代码中没有第二步, 而是对每个tile块矩阵乘之后直接排成了第二步图5里的bottom\_blob\_tm的样子. 我特地增加第二步是为了容易理解.

继续讲第三步的内存重排. 为了说明该步的操作细节, 将之前输入的4个tile块变成15个tile块. 其实重排的过程图6已经用不同的颜色块表示的比较清晰了. 一定要用语言来描述这个重排的过程话那就是: 对于bottom\_blob\_tm的每个tiles块, 沿着竖的h维度每8个元素划为一个"大块"(图6中红色的1-8表示), 横着放到bottom\_blob\_tm2的w维度上, 同时需要把其余c维度的"大块"也跟在w维度后面. tiles块剩下不足以构成8个元素的"大块", 就每4个元素构成一个"中块"(图6中绿色的17-20表示), 然后剩下不足以构成"中块"就以元素为单位同样操作(图6中蓝色的25-27表示).

## 步骤三: multi ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 8]

图7.multi step1

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 9]

图8.与图7等效的未transform的multi操作

然后进行两个矩阵的点对点相乘. 这两个矩阵分别是步骤一讲的卷积核kernel经过transform之后得到的矩阵kernel\_tm2和步骤二讲的输入bottom\_blob经过transform之后得到的bottom\_blob\_tm2, 也就是Winograd结论公式 Y = A T \[ ( G g G T ) ⊙ ( B T d B ) \] A = A T \[ U ⊙ V \] A = A T M A Y=A^T\[(GgG^T)⊙(B^TdB)\]A =A^T\[U⊙V\]A=A^TMA Y=AT\[(GgGT)⊙(BTdB)\]A=AT\[U⊙V\]A=ATMA里的 U U U和 V V V. 我们需要比较清楚的了解这两个transform之后的矩阵中每个元素对应于transform之前是什么位置的, 才能理解这个操作的正确性. 然后进一步理解经过transform之后实现点对点相乘操作的高效性. 如图7所示, kernel\_tm2里所涂色的4个颜色块表示原始矩阵kernel中"4个chw卷积核, 并且各个卷积核的8x8hw矩阵的第0个元素", bottom\_blob\_tm2里涂色的8个颜色块表示原始矩阵bottom\_blob中"8个不同8x8tile矩阵的第0个元素". 图7的操作可以按照以下描述: kernel\_tm2的4个颜色块分别与bottom\_blob\_tm2的8个颜色块进行乘法, 结果保存到一个长度为4\*8的临时变量中. 之所以没有把这个临时变量涂色涂完全, 是因为现在只计算了第0个c, 我们知道卷积的时候需要沿着c维度进行累加. 还需要图9所示将第1个c结果也累加起来, 这个临时变量才算是完全. 图7表示的经过transform后的multi操作等效于图8表示的未经transform的multi操作.  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 10]

图9.multi step2

图9表示所有的c都累加到临时变量后才构成完全上色的完整的临时变量, 然后将这个临时变量保存到长的很类似bottom\_blob\_tm的top\_blob\_tm中. 图9应该描述的很清楚.

## 步骤四: ouput\_transform ##

Y = A T \[ ( G g G T ) ⊙ ( B T d B ) \] A = A T \[ U ⊙ V \] A = A T M A Y=A^T\[(GgG^T)⊙(B^TdB)\]A =A^T\[U⊙V\]A=A^TMA Y=AT\[(GgGT)⊙(BTdB)\]A=AT\[U⊙V\]A=ATMA

M A = \[ m 0 m 1 m 2 m 3 m 4 m 5 m 6 m 7 m 8 m 9 m 10 m 11 m 12 m 13 m 14 m 15 m 16 m 17 m 18 m 19 m 20 m 21 m 22 m 23 m 24 m 25 m 26 m 27 m 28 m 29 m 30 m 31 m 32 m 33 m 34 m 35 m 36 m 37 m 38 m 39 m 40 m 41 m 42 m 43 m 44 m 45 m 46 m 47 m 48 m 49 m 50 m 51 m 52 m 53 m 54 m 55 m 56 m 57 m 58 m 59 m 60 m 61 m 62 m 63 \] \[ 1 0 0 0 0 0 1 1 1 1 1 1 1 − 1 1 − 1 1 − 1 1 2 4 8 16 32 1 − 2 4 − 8 16 − 32 32 16 8 4 2 1 32 − 16 8 − 4 2 − 1 0 0 0 0 0 1 \] MA = \\left\[ \\begin\{matrix\} m\_0 & m\_1 & m\_2 & m\_3 & m\_4 & m\_5 & m\_6 & m\_7\\\\ m\_8 & m\_9 & m\_\{10\} & m\_\{11\} & m\_\{12\} & m\_\{13\} & m\_\{14\} & m\_\{15\}\\\\ m\_\{16\} & m\_\{17\} & m\_\{18\} & m\_\{19\} & m\_\{20\} & m\_\{21\} & m\_\{22\} & m\_\{23\}\\\\ m\_\{24\} & m\_\{25\} & m\_\{26\} & m\_\{27\} & m\_\{28\} & m\_\{29\} & m\_\{30\} & m\_\{31\}\\\\ m\_\{32\} & m\_\{33\} & m\_\{34\} & m\_\{35\} & m\_\{36\} & m\_\{37\} & m\_\{38\} & m\_\{39\}\\\\ m\_\{40\} & m\_\{41\} & m\_\{42\} & m\_\{43\} & m\_\{44\} & m\_\{45\} & m\_\{46\} & m\_\{47\}\\\\ m\_\{48\} & m\_\{49\} & m\_\{50\} & m\_\{51\} & m\_\{52\} & m\_\{53\} & m\_\{54\} & m\_\{55\}\\\\ m\_\{56\} & m\_\{57\} & m\_\{58\} & m\_\{59\} & m\_\{60\} & m\_\{61\} & m\_\{62\} & m\_\{63\}\\\\ \\end\{matrix\} \\right\] \\left\[ \\begin\{matrix\} 1 & 0 & 0 & 0 & 0 & 0\\\\ 1 & 1 & 1 & 1 & 1 & 1\\\\ 1 & -1 & 1 & -1 & 1 & -1\\\\ 1 & 2 & 4 & 8 & 16 & 32\\\\ 1 & -2 & 4 & -8 & 16 & -32\\\\ 32 & 16 & 8 & 4 & 2 & 1\\\\ 32 & -16 & 8 & -4 & 2 & -1\\\\ 0 & 0 & 0 & 0 & 0 & 1\\\\ \\end\{matrix\} \\right\] MA=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡m0m8m16m24m32m40m48m56m1m9m17m25m33m41m49m57m2m10m18m26m34m42m50m58m3m11m19m27m35m43m51m59m4m12m20m28m36m44m52m60m5m13m21m29m37m45m53m61m6m14m22m30m38m46m54m62m7m15m23m31m39m47m55m63⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡111113232001−12−216−1600114488001−18−84−40011161622001−132−321−11⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

计算的过程跟前面计算 d B dB dB的思路是一样的, 不赘述了. 也是 **“0和1不用乘”** 以及 **“提炼公共部分”** , 至于这个公共部分, 就看代码吧

// 计算MA
    for (int m = 0; m < 8; m++) {
        // 3对6个公共部分
        float tmp024a = output0_tm_1[0] + output0_tm_2[0];
        float tmp135a = output0_tm_1[0] - output0_tm_2[0];
    
        float tmp024b = output0_tm_3[0] + output0_tm_4[0];
        float tmp135b = output0_tm_3[0] - output0_tm_4[0];
    
        float tmp024c = output0_tm_5[0] + output0_tm_6[0];
        float tmp135c = output0_tm_5[0] - output0_tm_6[0];
    
        // 结果都是用带上系数的公共部分组合而成
        tmp[0][m] = output0_tm_0[0] + tmp024a + tmp024b + tmp024c * 32;
        tmp[2][m] = tmp024a + tmp024b * 4 + tmp024c * 8;
        tmp[4][m] = tmp024a + tmp024b * 16 + tmp024c + tmp024c;
    
        tmp[1][m] = tmp135a + tmp135b + tmp135b + tmp135c * 16;
        tmp[3][m] = tmp135a + tmp135b * 8 + tmp135c * 4;
        tmp[5][m] = output0_tm_7[0] + tmp135a + tmp135b * 32 + tmp135c;
        // ...
    }
    
    // 计算AT(MA)
    for (int m = 0; m < 6; m++){
        const float* tmp0 = tmp[m];
    
        float tmp024a = tmp0[1] + tmp0[2];
        float tmp135a = tmp0[1] - tmp0[2];
    
        float tmp024b = tmp0[3] + tmp0[4];
        float tmp135b = tmp0[3] - tmp0[4];
    
        float tmp024c = tmp0[5] + tmp0[6];
        float tmp135c = tmp0[5] - tmp0[6];
    
        output0[0] = bias0 + tmp0[0] + tmp024a + tmp024b + tmp024c * 32;
        output0[2] = bias0 + tmp024a + tmp024b * 4 + tmp024c * 8;
        output0[4] = bias0 + tmp024a + tmp024b * 16 + tmp024c + tmp024c;
    
        output0[1] = bias0 + tmp135a + tmp135b + tmp135b + tmp135c * 16;
        output0[3] = bias0 + tmp135a + tmp135b * 8 + tmp135c * 4;
        output0[5] = bias0 + tmp0[7] + tmp135a + tmp135b * 32 + tmp135c;
    
        output0 += outw;
    }

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 11]

图10.output transform

## 总结 ##

本文主要讲了ncnn实现Winograd的算法, 该实现主要包含两个部分, 一个是Winograd公式的矩阵之间的优化计算方式, 另外一个是为了SIMD操作优化而做的内存重排. 前者本文是通过展示矩阵乘公式以及附上ncnn实现代码完成了对应说明, 后者本文通过多幅图完成了形象说明. 希望能帮助大家更好的理解Winograd卷积

[Link 1]: https://blog.csdn.net/hsqyc/article/details/116136385?spm=1001.2014.3001.5501
[conv3x3s1_winograd64_transform_kernel_neon5]: https://github.com/Tencent/ncnn/blob/master/src/layer/arm/convolution_3x3.h
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center]: /images/20221021/5139c1d3776e444587632a784579dffe.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221021/bb32ddfd81e84d4e95f1cd221e7198b0.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 2]: /images/20221021/f9bf2d4e17004e689663036750c9f28c.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 3]: /images/20221021/c04dc92a900d404587b3bcdae65d2c8b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 4]: /images/20221021/24a4cf8225e544ea88f8bf5e34a570e0.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 5]: /images/20221021/fdac4a794b11486aabe475e3f0306281.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 6]: /images/20221021/ea42eab612984c56bd052a4b7d1f0a43.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 7]: /images/20221021/cf7943aecc8640f9ac2f137d156d3f72.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 8]: /images/20221021/b1657f168e734a43ae7dbba913abb8d5.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 9]: /images/20221021/620edc62888840e3b77a5a9dc377f6c0.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 10]: /images/20221021/265536b5b1e9430ea60c23a00e033a4b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2hzcXlj_size_16_color_FFFFFF_t_70_pic_center 11]: /images/20221021/528f3cd5eca24dce9432dff1fb6c927a.png