前言
卷积神经网络是很多任务尤其是计算机视觉任务的基础，但很大程度上，模型需要的大量卷积计算限制了模型的可用性。因此，如何快速的完成卷积操作就至关重要。
此处的卷积是指图像处理领域的卷积操作，且数据通常为多通道的二维数组，卷积核的长宽相等。
常用的优化方法包括三个方面：
- 硬件：堆数据堆模型再堆硬件是提升性能最常见的手段，本文假设读者都是买不起煤气灶的穷人
- 模型复杂度：降低模型复杂度包括降低参数冗余的花式的卷积设计，如ShuffleNet中的分组卷积+通道混合，MobileNet中的Depthwise和Pointwise分解卷积等，以及模型的裁剪、量化和稀疏化。
- 框架计算速度：目前主流的深度学习框架在加快计算上，一方面是利用数据SIMD的特性，进行硬件上并行化（SSE、Neon、线程并行），又或者通过一些矢量化手段（如Caffe、MXNet中的im2col）来充分利用软硬件的特点，实现更高的计算速度，另一方面则是像FFT、Strassen算法以及Winograd算法等在卷积计算原理上不同的方法，从而减少了计算量。
本文介绍的 Winograd 是存在已久最近被重新发现的方法（The Coppersmith-Winograd Matrix Multiplication Algorithm），在大部分场景中，Winograd 方法都显示和较大的优势，目前TF Lite、Tencent NCNN、Ali MNN 中计算卷积就使用了该方法。详述该方法并进行测试的是CVPR 2016中的一篇文章《Fast Algorithms for Convolutional Neural Networks》，本文主要以该文章来进行方法的介绍和讲解。
贴两张ARM报告的图来说明一下问题，换言之，我们的问题是需要在模型（耗时操作主要是小型卷积核如3x3的卷积）和设备确定的情况下，在Inference阶段尽可能快地计算卷积结果。
卷积问题定义
符号
$Image_{H\cdot W\cdot C}$:表示一张多通道（通常为三通道的）图片，使用$I$表示，下标视情况省略，当$C$为1时，表示灰度图，当$W$为1时，表示一维数据
$Kernel_{R\cdot R}$：表示一个大小为$R$的二维卷积核，使用$K$表示
$Height/Row,Width/Col,Channel/Depth$：分别表示一张图片的行，列，通道数
$I_{OH\cdot OW \cdot OC}=Conv(I_{H\cdot W \cdot C},K_R)$：表示使用卷积核$K$对图像$I$进行卷积操作，得到一张新的图像，为了和论文中一致，我们也使用$F(I_{OH\cdot OW \cdot OC},K_R)$来表示上述操作。
一维卷积
- 形式:$I_{OH}=Conv(I_{H},K_R)$或$F(I_{OH},K_R)$
输出的图像大小取决于卷积核大小、Stride以及Padding策略，在这里我们假设Stride都是1，没有Padding，且只统计乘法运算，则
- 计算量：$FLO=OH\cdot R=(H-R+1)\cdot R$
二维卷积
- 形式:$I_{OH\cdot OW \cdot OC}=Conv(I_{H\cdot W\cdot C},K_R)$或$F(I_{OH\cdot OW \cdot OC},K_R)$
多个通道的二维卷积遵循层内连乘，层间累加的方法，需要的卷积核数量为$C\cdot OC$:
- 计算量：$FLO=(OH\cdot OW \cdot OC)\cdot (R\cdot R)\cdot C$
具体而言，输出图片的每个通道都是通过一组卷积核在所有通道上相乘累加得到的。
有兴趣的童鞋可以自己算一下需要的加法操作次数，智障作者算了一下应该是$OP=(OH\cdot OW \cdot OC)\cdot ((R\cdot R-1)\cdot C+C-1)=(OH\cdot OW \cdot OC)\cdot (R\cdot R\cdot C-1)$次，和乘法基本相同。
具体讨论算法的时候我们假设输入输出的通道数都是1，即$FLO=(OH\cdot OW)\cdot (R\cdot R)$。
- 形式:$I_{OH\cdot OW}=Conv(I_{H\cdot W},K_R)$或$F(I_{OH\cdot OW},K_R)$
注意点
在各种开源框架中，CNN中的conv2d层执行的并不是数学上的卷积计算，而是数学上的互相关计算，具体定义的区别请参考《深度学习》Chap9.1或者这个。
在传统的图像处理领域，卷积核的参数是已知的（且往往是对称的），而在卷积神经网络中则成为了待定参数。实际上，虽然说将卷积核上下翻转并左右翻转（即旋转180度）才是真正的卷积操作，但是也可以认为卷积层做的就是卷积，只是特征矩阵是倒序存储的。又或者说由于卷积核的参数是可变的，这样做不但提高了效率也不影响结果。
跑题：im2col实现
im2col是一种非常容易理解的矢量化（Vectorization）手段，基于im2col和GEMM（ General Matrix Multiplication ）的方法可以获得较正常卷积计算较高的加速比，具体而言就是把我们每次进行卷积操作时涉及到的元素展开一个列向量，最终得到一个$I_{(R\cdot R)\cdot (OH\cdot OW)}$的矩阵，卷积核则拆成$K_{1\cdot (R\cdot R)}$的向量。
按照同样的方法把卷积核展开，最后卷积操作就可以用一个矩阵乘法来表示，计算完成后，再使用col2im将结果转换为图片。
这里输入是3个通道，输出是1个，相比写4个for循环，im2col大大加快了计算速度（矩阵乘法可并行，数据在内存中的存储连续，Cache命中率提高），但是没有减少计算量，且内存占用几乎倍增，另外，生成多通道矩阵依然需要写5个循环（NHWCRR），算是一种用空间来换取时间的做法。关于Caffe中im2col和col2im的实现可以看这里，具体的实现细节其实和我们理解的过程是有较大不同的。
需要强调的一点是，由于我们的硬件设计就是良好支持并行化的，所以算得少和算得快并不能简单的认为是一回事。
Caffe原作者贾扬清对其框架中卷积实现的吐槽可以看这里，摘录一段如下：
In the last few months chatting with people about Caffe, a common comment I got was: “Caffe’s convolution has some memory issues.”
While this is true in some sense, I am not sure whether it is truly an issue - rather, it is a graduate-student level design choice when I was writing the Caffe framework in just 2 months’ budget with a looming thesis deadline. It turns out to have its own pros （faster than any trivial implementation unless you optimize really seriously）and cons （large memory consumption） A more detailed explanation follows, if you are interested.
……
Thus, I took a simpler approach: reduce the problem to a simpler one, where others have already optimized it really well.
关于傅里叶变换计算卷积的方法在卷积神经网络中并不是很常用（因为一般还没直接算快），可以参看这个，Strassen则是在矩阵较大时具有较好的加速效果，具体的方法细节之后再详细说
（因为作者也没完全看懂╮(￣▽￣)╭）。
Winograd
方法
在我们的硬件受限的情况下（例如只有单核CPU），如何才能变得更高更快更强呢，由于计算一次乘法所需要的时钟周期要比计算一次加法大很多，对于乘法操作我们能够想象的下界是：输入数据的每个元素至少参与一次乘法。这也是我们力求达到的目标，即乘法次数尽可能的靠近下界。
形式化的表述，针对一个输出长度为$m$，卷积核大小为$r$的卷积运算，其所需要的最小乘法次数与输入数据的长度相同：
$$\mu(F(I_m, K_r))=m+r-1=H$$
每个维度相对独立，因此拓展到二维的情况也是一样的：
$$\begin{aligned} \mu(F(I_{m \cdot n}, K_{r \cdot s})) =\mu(F(I_{m}, K_{r})) \mu(F(I_{n}, K_{s})) =(m+r-1)(n+s-1) \end{aligned}=H \cdot W$$
该理论最早由Shmuel Winograd于1980年提出，是本文最重要的理论，没有之一。具体而言，给定一个确定的卷积问题，我们可以找到一种变换，来接近（甚至到达）乘法次数的理论下界。根据该理论，在卷积核较小的情况下，我们可以获得相对普通卷积计算而言较大的加速比。
$F(I_2,K_3)$与$F(I_{2\cdot 2},K_{3 \cdot 3})$
举个例子，$F(I_2,K_3)$表示输入信号$d=\left[ \begin{array}{llll}{d_{0}} & {d_{1}} & {d_{2}} & {d_{3}}\end{array}\right]^{T}$，卷积核$g=\left[ \begin{array}{lll}{g_{0}} & {g_{1}} & {g_{2}}\end{array}\right]^{T}$的一维卷积操作。那么问题就可以表示为如下的形式：
$$F(I_2,K_3) = \left[ \begin{array}{lll}{d_{0}} & {d_{1}} & {d_{2}} \\ {d_{1}} & {d_{2}} & {d_{3}}\end{array}\right] \left[ \begin{array}{l}{g_{0}} \\ {g_{1}} \\ {g_{2}}\end{array}\right]=\left[ \begin{array}{c}{r_0} \\ {r_1}\end{array}\right]$$
如果是一般的矩阵乘法，则需要6次乘法和4次加法，如下：
$$\begin{array}{l}{r_{0}=\left(d_{0} \cdot g_{0}\right)+\left(d_{1} \cdot g_{1}\right)+\left(d_{2} \cdot g_{2}\right)} \\ {r_{1}=\left(d_{1} \cdot g_{0}\right)+\left(d_{2} \cdot g_{1}\right)+\left(d_{3} \cdot g_{2}\right)}\end{array}$$
但是，卷积运算中输入信号转换成的矩阵不是任意矩阵，其中有规律地分布着大量的重复元素，因此卷积转换成的矩阵乘法比一般矩阵乘法的问题域更小，这就让优化存在了可能。
Winograd的具体做法是，
$$F(I_2,K_3) = \left[ \begin{array}{lll}{d_{0}} & {d_{1}} & {d_{2}} \\ {d_{1}} & {d_{2}} & {d_{3}}\end{array}\right] \left[ \begin{array}{l}{g_{0}} \\ {g_{1}} \\ {g_{2}}\end{array}\right]=\left[ \begin{array}{c}{m_{1}+m_{2}+m_{3}} \\ {m_{2}-m_{3}-m_{4}}\end{array}\right]$$
其中，
$$\begin{array}{ll}{m_{1}=\left(d_{0}-d_{2}\right) g_{0}} & {m_{2}=\left(d_{1}+d_{2}\right) \frac{g_{0}+g_{1}+g_{2}}{2}} \\ {m_{4}=\left(d_{1}-d_{3}\right) g_{2}} & {m_{3}=\left(d_{2}-d_{1}\right) \frac{g_{0}-g_{1}+g_{2}}{2}}\end{array}$$
在神经网络的推理阶段，卷积核$g$上的元素是固定的，共3次加法与2次乘法，因此$g$相关的运算可以提前算好，预测阶段只需计算一次，可以忽略，所以一共所需的运算次数4次乘法和8次加法/减法，和下界一致。 速度提升比例为1.5倍。当然所需要的存储空间也相应的增加，由原来的$r$增加为$H$，即从卷积核大小变成了输入数据的宽度。
可是这怎么就能减少计算量呢，会不会得到的结果不一样？有兴趣的童鞋可以自行把所有符号代入化简，针对这个例子，得到的结果是完全相同的！当然在具体计算的时候可能会存在浮点精度误差。因此在这个例子里面完全就是仅仅使用不同的数学表示就达到了减少计算量的目的，真是太巧妙了。
一维卷积形式化表述
上述例子可以用通用矩阵乘法（GEMM）和元素级乘法（EWMM）的混合矩阵变换来形式化表述。
一维卷积的形式化表述如下：
$$Y=A^{T}\left[(G g) \odot\left(B^{T} d\right)\right]$$
详细解释如下：
- $g$：卷积核 ，维度$r\cdot 1$
- $d$：输入信号，维度$H\cdot 1$
- $G$：Filter transform矩阵，用于将卷积核变换到另一个空间，维度$(m+r-1)\cdot r=H \cdot r$
- $B^T$：Input transform矩阵，用于将输入数据变换到另一个空间，维度$(m+r-1)\cdot (m+r-1)=H \cdot H$
- $A^T$：Output transform矩阵，用于将数据转换回输出空间，维度$m\cdot (m+r-1)=m*H$
- $Y$：卷积结果，维度$m\cdot 1$
其中，$\odot$表示对应位置相乘，暗示此处两个矩阵的维度是一致的，上面的例子里面具体的参数如下，
$$
B^{T}=\left[ \begin{array}{cccc}{1} & {0} & {-1} & {0} \\ {0} & {1} & {1} & {0} \\ {0} & {-1} & {1} & {0} \\ {0} & {1} & {0} & {-1}\end{array}\right]
$$
$$
G=\left[ \begin{array}{ccc}{1} & {0} & {0} \\ {\frac{1}{2}} & {\frac{1}{2}} & {\frac{1}{2}} \\ {\frac{1}{2}} & {-\frac{1}{2}} & {\frac{1}{2}} \\ {0} & {0} & {1}\end{array}\right]
$$
$$
A^{T}=\left[ \begin{array}{llll}{1} & {1} & {1} & {0} \\ {0} & {1} & {-1} & {-1}\end{array}\right]
$$
$$
g=\left[ \begin{array}{lll}{g_{0}} & {g_{1}} & {g_{2}}\end{array}\right]^{T},d=\left[ \begin{array}{llll}{d_{0}} & {d_{1}} & {d_{2}} & {d_{3}}\end{array}\right]^{T}
$$
整个计算过程在逻辑上可以分为4步：
- Input transform
- Filter transform
- Hadamard product（ 哈达玛积 ）
- Output transform
此处，$A^{T}$即$m$前面的系数，$B^{T}$即数据$d$前面的系数，$G$即卷积核$g$前面的系数，由于卷积核相关参数是提前计算好的，虽然我们将卷积的过程表述成了矩阵乘法的形式，但是此处只有$\odot$包含了乘法，其他的系数全是$\pm 1$，因此只有加减法。之后会提到具体如何实现，有兴趣的同学可以先自己思考一下下。
二维卷积形式化表述
文中关于一维卷积向二维卷积的扩展只有寥寥数语：
A minimal 1D algorithm $F(m, r)$ is nested with itself to obtain a minimal 2D algorithm.
The nesting technique can be generalized for non-square filters and outputs,$F(m × n, r × s)$,
by nesting an algorithm for $F(m, r)$ with an algorithm for $F(n, s)$.
文中给出的二维卷积形式化表述如下：
$$
Y=A^{T}\left[\left[G g G^{T}\right] \odot\left[B^{T} d B\right]\right] A
$$
详细解释如下：
- $g$：卷积核 ，维度$r\cdot r$
- $d$：输入信号，维度$H\cdot H$
- $G$：Filter transform矩阵，用于将卷积核变换到另一个空间，维度$(m+r-1)\cdot r=H \cdot r$
- $B^T$：Input transform矩阵，用于将输入数据变换到另一个空间，维度$(m+r-1)\cdot (m+r-1)=H \cdot H$
- $A^T$：Output transform矩阵，用于将数据转换回输出空间，维度$m\cdot (m+r-1)=m*H$
- $Y$：卷积结果，维度$m\cdot m$
对于$F(I_{2\cdot 2},K_{3 \cdot 3})$，其中所有矩阵中的参数，包括$G$，$B^T$，$A^T$都是和$F(I_2,K_3)$一样的。
依然只有$\odot$包含了乘法，乘法的次数为16次，相比于标准卷积的次数$36=2\times 2 \times 3 \times 3$，速度提升比例为2.25倍。其中Input transform包括$32=4\times 4 \times 2$次加法，Filter transform包括28次浮点数操作（预先计算），Output transform（Inverse transform）包括$24=2\times 2 \times 3 \times 2$次加法。
等等，怎么就nested with itself了，参数还是一样的？
那么问题来了：
- 二维卷积的形式化表述是否正确，如何证明
- 通常输入的图像尺寸都较大，如何使用Winograd对其进行卷积计算，难道要实现$F(I_{448\cdot 448}, K_{3 \cdot 3})$吗，三维的卷积（卷积神经网络中的实际情况）该如何实现
- 上面形式化表示中所使用的矩阵参数该如何获得
- 形式化表示中的矩阵乘法看起来似乎比常规卷积做了更多的乘法操作，真正的算法实现是如何转换成加法的
二维卷积形式化推导
我们先来解决第一个问题，要完成形式化的推导，我们需要先理解文章中的nested with itself到底是什么意思。下列图片来自ARM在Embedded Vision Summit 2018上的Slides，里面的符号表示会略有不同（ 用$k$来表示输入，$w$表示权重，$r$表示输出 ）。
我们以$F(I_{2\cdot 2},K_{3 \cdot 3})$为例，输入图片是$4\times 4$的，输出是$2\times 2$的，卷积核如下：
$$
W = \left[\begin{array}{lll}{w_{0}} & {w_{1}} & {w_{2}} \\ {w_{3}} & {w_{4}} & {w_{5}} \\ {w_{6}} & {w_{7}} & {w_{8}}\end{array}\right]
$$
根据我们之前对标准卷积的运算次数推导，$FLO=(OH\cdot OW \cdot OC)\cdot (R\cdot R)\cdot C=36$，$OP=(OH\cdot OW \cdot OC)\cdot (R\cdot R\cdot C-1)=32$，即需要36次乘法和32次加法。
现在我们按照im2row的形式展开卷积流程，可以得到如下的矩阵运算：
仔细观察可以发现，左侧矩阵中的部分元素是重复出现的，我们按照相同的颜色对其进行标记，并以此来进行矩阵和向量的分块操作：
使用更加简洁的表述，我们得到了如下所示的分块运算，现在问题的表述和$F(I_2,K_3)$完全一致了，不同的是我们的每一对元素操作都是$F(I_2,K_3)$，这就是$F(I_{2\cdot 2},K_{3 \cdot 3})$的堆叠实现：
形式化的表述如下：
$$
\begin{aligned} \left[\begin{array}{lll}{K_0} & {K_1} & {K_2} \\ {K_1} & {K_2} & {K_3} \end{array}\right] \left[\begin{array}{l}{W_0} \\ {W_1} \\ {W_2} \end{array}\right] &= \left[\begin{array}{l}{R_0} \\ {R_1} \end{array}\right] = \left[\begin{array}{l}{K_0W_0+K_1W_1+K_2W_2} \\ {K_1W_0+K_2W_1+K_3W_2} \end{array}\right] \\ \\ &= \left[\begin{array}{l}{F_{(2,3)}(D_0,W_0)+F_{(2,3)}(D_1,W_1)+F_{(2,3)}(D_2,W_2)} \\ {F_{(2,3)}(D_1,W_0)+F_{(2,3)}(D_2,W_1)+F_{(2,3)}(D_3,W_2)} \end{array}\right] \end{aligned}
$$
其中，$D_i$是$K_i$对应的输入序列，也即卷积输入的第$i$行:
$$
D = d^T=\left[\begin{array}{llll} {k_0} & {k_4} & {k_8} & {k_{12}} \\ {k_1} & {k_5} & {k_9} & {k_{13}} \\ {k_2} & {k_6} & {k_{10}} & {k_{14}} \\ {k_3} & {k_7} & {k_{11}} & {k_{15}} \end{array}\right] = \left[\begin{array}{l} D_0 & D_1 & D_2 & D_3 \end{array}\right]
$$
我们一共使用6个$F(I_2,K_3)$来计算$F(I_{2\cdot 2},K_{3 \cdot 3})$，并额外增加了8次加法，之前提到$F(I_2,K_3)$一共所需的运算次数为4次乘法和8次加法，所以计算时期共计24次乘法与48次加法， 速度提升比例为1.5倍。在卷积核预处理阶段，需要进行8次乘法和12次加法。
有兴趣的童鞋可以思考一下为什么是48次加法，因为数据有部分重复，数据变换（Input transform）只需要做4次，此处共计16次加法，而Output transform则是6次，共计24次加法，另外还有8次额外的加法（虽然只有4个加号却是8次）。
之前曾经提到，最小的乘法次数应该和输入的数据规模相等，因此这还不是最高的加速比（最高应该是$2.25=36\div 16$倍），还记得我们之前的一维卷积形式化表述吗，下面的图片中的所有表述和$F(I_2,K_3)$是完全一致的，但是每个元素都用矩阵和向量替换了：
按照这种方式，我们可以得到$F(I_{2\cdot 2},K_{3 \cdot 3})$的嵌套实现，形式化的表述如下：
$$
\begin{aligned} \left[ \begin{array}{c}{R_0} \\ {R_1}\end{array}\right] &= \left[ \begin{array}{c}{K_0 W_0 + K_1 W_1 + K_2 W_2} \\ {K_1 W_0 + K_2 W_1 + K_3 W_2} \end{array} \right] \\ \\ &= \left[\begin{array}{l}{F_{(2,3)}(D_0,W_0)+F_{(2,3)}(D_1,W_1)+F_{(2,3)}(D_2,W_2)} \\ {F_{(2,3)}(D_1,W_0)+F_{(2,3)}(D_2,W_1)+F_{(2,3)}(D_3,W_2)} \end{array}\right] \\ &= \left[ \begin{array}{c} {A^{T}\left[(G W_0) \odot\left(B^{T} D_0 \right)\right] + A^{T}\left[(G W_1) \odot\left(B^{T} D_1 \right)\right] + A^{T}\left[(G W_2) \odot\left(B^{T} D_2 \right)\right]} \\ {A^{T}\left[(G W_0) \odot\left(B^{T} D_1 \right)\right] + A^{T}\left[(G W_1) \odot\left(B^{T} D_2 \right)\right] + A^{T}\left[(G W_2) \odot\left(B^{T} D_3 \right)\right]} \end{array} \right] \\ \\ &=A^{T}\left[\left[G [W_0 \ W_1 \ W_2 ] G^{T}\right] \odot\left[B^{T} [D_0 \ D_1 \ D_2 \ D_3] B\right]\right]A \\ \\ &=A^{T}\left[\left[G w G^{T}\right] \odot\left[B^{T} d B\right]\right] A \\ \\ &\textit{(…w => g…)} \\ \\ &=A^{T}\left[\left[G g G^{T}\right] \odot\left[B^{T} d B\right]\right] A \end{aligned}
$$
中间的一步变化很关键，$\left[(G W_i) \odot\left(B^{T} D_j \right)\right]$是一个长度为4的列向量,$A^{T}\left[(G W_i) \odot\left(B^{T} D_j \right)\right]$则是一个长度为2的列向量，$A^{T}\left[(G W_0) \odot\left(B^{T} D_0 \right)+ (G W_1) \odot\left(B^{T} D_1 \right) + (G W_2) \odot\left(B^{T} D_2 \right)\right]$ 方括号内对应位置相乘再相加，相当于在每组相点乘结果构成的行向量上做卷积。最后的结果是一个长度为2的列向量。
实际上两种表述的维度并不相同，前者的维度是$I_{4\cdot 1}$，后者是$I_{2\cdot 2}$。
此处的推导过于复杂，作者并不会，有兴趣的童鞋可以去试一下这段代码，其作用是代入$F(I_{2\cdot 2},K_{3 \cdot 3})$的参数和符号，进行两种表达的展开，最后证明相等:
from sympy import Symbol, Matrix,pprint,simplify
import numpy as np
BT = Matrix([
[1, 0, -1, 0],
[0, 1, 1, 0],
[0, -1, 1, 0],
[0, 1, 0, -1]
])
G = Matrix([
[2, 0, 0],
[1, 1, 1],
[1, -1, 1],
[0, 0, 2]
])
G=G/2
AT = Matrix([
[1, 1, 1, 0],
[0, 1, -1, -1]
])
g = Matrix(3, 3, [Symbol(f'g{i}') for i in range(3*3)])
d = Matrix(4, 4, [Symbol(f'd{i}') for i in range(4*4)])
m = Matrix(4, 4, [Symbol(f'm{i}') for i in range(4*4)])
print('GgGT:')
GgGT=G*g*(G.T)
pprint(GgGT)
print('BTdB')
BTdB=BT*d*(BT.T)
pprint(BTdB)
RET1=(AT*np.multiply(GgGT,BTdB)*(AT.T))
print('AT * [GgGT em BTdB] * A ,shape={shape},first:'.format(shape=RET1.shape))
pprint(simplify(RET1[0]))
R0=(AT*np.multiply(G*g[0,:].T,BT*d[0,:].T)+AT*np.multiply(G*g[1,:].T,BT*d[1,:].T)+AT*np.multiply(G*g[2,:].T,BT*d[2,:].T))
R1=(AT*np.multiply(G*g[0,:].T,BT*d[1,:].T)+AT*np.multiply(G*g[1,:].T,BT*d[2,:].T)+AT*np.multiply(G*g[2,:].T,BT*d[3,:].T))
print('R0 ,shape={shape},first:'.format(shape=R0.shape))
pprint(simplify(R0[0]))
print('Diff of first:')
pprint(simplify(RET1[0]-R0[0]))
print('Diff of all:')
print(simplify(RET1-np.vstack((R0.T,R1.T))))
经过这样一番推导（虽然最关键的一步跳过了，嘻嘻嘻）我们就可以得到$F(I_{2\cdot 2},K_{3 \cdot 3})$的嵌套实现，
$$
F(I_{2\cdot 2},K_{3 \cdot 3})= A^{T} \left[ U \odot V \right] A
$$
其中，$U = G g G^{T},V = B^{T} d B$，一共需要 16次乘法和56次加法（$V = B^{T} d B$过程32次加法、$M=U \odot V$过程16次乘法、$Y=A^TMA$过程24次加法）。和一维卷积类似的，所需要的存储空间也相应的增加，由原来的$r\cdot r$增加为$H\cdot W$，即从卷积核大小变成了输入数据的宽度。
扩展到$F(I_{4\cdot 4},K_{3 \cdot 3})$
算法本身可以被扩展到更大的输入输出尺寸，从而可能得到更高的加速比，但是参数数量和加法操作次数也相应的大大增加。
要计算$F(I_{4\cdot 4},K_{3 \cdot 3})$，我们的输入是$6\times 6$的图片，常规卷积的操作需要计算$144=4\times 4\times 3\times3$次乘法，但是使用Winograd只需要$36=6\times 6$次乘法，速度提升比例为4倍。具体参数如下：
$$
B^{T}=\left[\begin{array}{rrrrrr}{4} & {0} & {-5} & {0} & {1} & {0} \\ {0} & {-4} & {-4} & {1} & {1} & {0} \\ {0} & {4} & {-4} & {-1} & {1} & {0} \\ {0} & {-2} & {-1} & {2} & {1} & {0} \\ {0} & {2} & {-1} & {-2} & {1} & {0} \\ {0} & {4} & {0} & {-5} & {0} & {1}\end{array}\right]
$$
$$
G=\left[\begin{array}{rrr}{\frac{1}{4}} & {0} & {0} \\ {-\frac{1}{6}} & {-\frac{1}{6}} & {-\frac{1}{6}} \\ {-\frac{1}{6}} & {\frac{1}{6}} & {-\frac{1}{6}} \\ {\frac{1}{24}} & {\frac{1}{12}} & {\frac{1}{6}} \\ {\frac{1}{24}} & {-\frac{1}{12}} & {\frac{1}{6}} \\ {0} & {0} & {1}\end{array}\right]
$$
$$
A^{T}=\left[\begin{array}{rrrrrr}{1} & {1} & {1} & {1} & {1} & {0} \\ {0} & {1} & {-1} & {2} & {-2} & {0} \\ {0} & {1} & {1} & {4} & {4} & {0} \\ {0} & {1} & {-1} & {8} & {-8} & {1}\end{array}\right]
$$
The number of additions and constant multiplications required by the minimal Winograd transforms increases quadratically with the tile size. Thus for large tiles, the complexity of the transforms will overwhelm any savings in the number of multiplications.
The magnitude of the transform matrix elements also increases with increasing tile size. This effectively reduces the numeric accuracy of the computation, so that for large tiles, the transforms cannot be computed accurately.
但是浮点操作次数也大大增加了，其中Input transform包括$144=12\times (6+6)$次浮点数操作（注意此时的系数已经不全是1了，因此也会包括乘法），Filter transform包括$72=8\times (3+6)$次浮点数操作（预先计算），Output transform（Inverse transform）包括$100=10\times (6+4)$次浮点数操作。
此外，转换矩阵的规模增大导致了计算精度误差增加，不过作者认为卷积神经网络对精度的要求其实比较低，因此在附录中讨论了$F(I_{6\cdot 6},K_{3 \cdot 3})$的可能性，其参数如下：
$$
B^{T}=\left[\begin{array}{rrrrrrrr} {1} & {0} & {-21/4} & {0} & {21/4} & {0} & {-1} & {0} \\ {0} & {1} & {1} & {-17/4} & {-17/4} & {1} & {1} & {0} \\ {0} & {-1} & {1} & {17/4} & {-17/4} & {-1} & {1} & {0} \\ {0} & {1/2} & {1/4} & {-5/2} & {-5/4} & {2} & {1} & {0} \\ {0} & {-1/2} & {1/4} & {5/2} & {-5/4} & {-2} & {1} & {0} \\ {0} & {2} & {4} & {-5/2} & {-5} & {1/2} & {1} & {0} \\ {0} & {-2} & {4} & {5/2} & {-5} & {-1/2} & {1} & {0} \\ {0} & {-1} & {0} & {21/4} & {0} & {-21/4} & {0} & {1} \end{array}\right], \\ G=\left[\begin{array}{rrr} {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} \end{array}\right], \\ A^{T}=\left[\begin{array}{rrrrrrrr} {1} & {1} & {1} & {1} & {1} & {1} & {1} & {0} \\ {0} & {1} & {-1} & {2} & {-2} & {1/2} & {-1/2} & {0} \\ {0} & {1} & {1} & {4} & {4} & {1/4} & {1/4} & {0} \\ {0} & {1} & {-1} & {8} & {-8} & {1/8} & {-1/8} & {0} \\ {0} & {1} & {1} & {16} & {16} & {1/16} & {1/16} & {0} \\ {0} & {1} & {-1} & {32} & {-32} & {1/32} & {-1/32} & {1} \\ \end{array}\right]
$$
相比直接卷积324次的乘法操作，$F(I_{6\cdot 6},K_{3 \cdot 3})$只需要64次，加速比达到5.06倍（然而并没有，实际有3就很不错了）。我们目前讨论的都是$3\times 3$的卷积，只是tile大小不同。
基于Winograd的卷积计算算法流程
第二个问题，其实我们并不会真的去实现$F(I_{448\cdot 448}, K_{3 \cdot 3})$，而是将图片划分成多个相同大小部分重叠的tile，在此基础上使用如$F(I_{2\cdot 2},K_{3 \cdot 3})$的方式来计算，最后合并统计结果。对于三维卷积，实际上是和标准卷积一样，逐层做二维卷积，再每层对应位置结果相加。但除此之外，针对多个卷积核还有更加巧妙的做法。
算法流程
tile在此处代表一个图像块，具体划分方式如下：
- 每个
tile的大小为$H\times H=(m+r-1)(m+r-1)$ - 每个
tile互相重叠的长度为$r-1$ - 每个通道可以划分成$P=\lceil H / m\rceil\lceil W / m\rceil $个
tile，注意此处的$H$表示的是通道的大小，和之前的不同，之后的算法表述也会出现类似的情况
之所以有部分重叠，是因为我们的卷积结果其实小于原图的大小，而每一边减小的大小就是$r-1$。
下图展示了标准卷积和Winograd $F(I_{2\cdot 2},K_{3 \cdot 3})$的区别,标准的卷积过程：
Winograd $F(I_{2\cdot 2},K_{3 \cdot 3})$:
之前提到，计算公式如下所示：
$$
F(I_{m\cdot m},K_{r \cdot r})= A^{T} \left[ U \odot V \right] A \\U = G g G^{T}\\V = B^{T} d B
$$
- Winograd的输入为固定大小，因此需要先对图像做Padding，每个
tile的步长即为非重叠的长度$m$ - $U$已经提前计算好，$V$每个
tile计算一次 - 最后输出的图像部分是不需要的，需要做裁剪，裁剪宽度和Padding宽度相同
前述只讨论了一些比较简单的情况，事实上在CNN中，由于输入的特征图只需要变换一次，而却会被多个滤波器复用，所以输入变换过程的额外开销会被平摊——卷积的滤波器（也即输出通道）越多，那么输入变换产生的额外开销的影响就越小。
多通道多卷积核
在图片$i$，通道$c$，卷积核$k$上做标准卷积可以用如下的形式表示：
$$
{Y_{i, k, x, y}=\sum_{c=1}^{C} \sum_{v=1}^{R} \sum_{u=1}^{S} D_{i, c, x+u, y+v} G_{k, c, u, v}}
$$
我们可以使用星号来表示2D互相关的过程，
$$
Y_{i, k}=\sum_{c=1}^{C} D_{i, c} * G_{k, c}
$$
我们使用$(\tilde{x}, \tilde{y})$来表示tile的坐标，则Winograd卷积可以被表示成如下形式：
$$
{\qquad \begin{aligned} Y_{i, k, \widetilde{x}, \widetilde{y}} &=\sum_{c=1}^{C} D_{i, c, \widetilde{x}, \widetilde{y}} * G_{k, c} \\ &=\sum_{c=1}^{C} A^{T}\left[U_{k, c} \odot V_{c, i, \widetilde{x}, \bar{y}}\right] A \\ &=A^{T}\left[\sum_{c=1}^{C} U_{k, c} \odot V_{c, i, \widetilde{x}, \widetilde{y}}\right] A \end{aligned}}
$$
因此，我们可以在变换的空间就对通道的那个维度求和，不需要等到输出变换（Inverse Transform）。
$$
M_{k, i, \widetilde{x}, \widetilde{y}}=\sum_{c=1}^{C} U_{k, c} \odot V_{c, i, \widetilde{x}, \widetilde{y}}
$$
按照图像/tile坐标$(i,\tilde{x}, \tilde{y})$的方式将数据全部展开至一个维度，我们可以将之前的遍历操作转换得到如下的矩阵乘法表示，其中$(\xi, \nu)$遍历的是每个tile大小的区域，表示一个逐点相乘的操作：
$$
M_{k, b}^{(\xi, \nu)}=\sum_{c=1}^{C} U_{k, c}^{(\xi, \nu)} V_{c, b}^{(\xi, \nu)}
$$
该表述实际上就是矩阵乘法，在下文的图中可以跟清晰的看出来，现在计算过程的表述如下，整个过程一共要重复tile大小的次数：
$$
M^{(\xi, \nu)}=U^{(\xi, \nu)} V^{(\xi, \nu)}
$$
综上，文章中描述的算法流程如下所示，其中的$K$表示卷积核的数量，也是输出的特征图的深度（变换以后的，真实的卷积核有$K\times C$个），$C$表示图片通道数，$P$代表tile的个数：
在文章《 Sparse Winograd Convolutional neural networks on small-scale systolic arrays 》中描述了多通道多卷积核的Winograd计算流程，和单通道不同的是，此处我们依然是使用GEMM来计算整个卷积过程，因此针对im2col中GEMM部分的优化在此处依然适用。最后我们是从所有的卷积结果里面拿出了一个tile应有的输出。
算法整体依然分为四步：
- Filter Transform
- Input transform
- Hadamard product $->$ Batched-GEMM（批量矩阵乘法）
- Output transform
Winograd中和矩阵乘法不同的地方主要是有一个元素级乘法（EWMM）的操作，而在这种使用GEMM表示的计算流程中，巧妙的将多通道和多卷积核使用矩阵表示，参与GEMM操作的矩阵尺寸分别为$\lceil H / m\rceil\lceil W / m\rceil \times C$和$C \times K$，一共做了$(m+r-1)^2=H\times W$个矩阵乘法操作（此处为Winograd输入大小，和图像大小不同），通道的那个维度在矩阵乘法中就被消掉了。而元素级乘法（EWMM）的操作实际上就在$(m+r-1)^2=H\times W$那么多个矩阵乘法操作中进行。无论是单个GEMM操作内部，还是多个GEMM操作都可以使用GEMM的并行技巧。
参数推导（劝退节）-TODO
当你在凝视深渊的时，Winograd也在凝视着你。
Shmuel Winograd, 1936年生，MIT本硕（1959年）+纽约大学博士（1968年）。此后一直在IBM当研究员，直到退休。IEEE Fellow，ACM Fellow，美国科学院院士。
欧几里得定理（辗转相除法）
欧几里得定理（ Euclidean algorithm ），又名辗转相除法，用来快速求两个数的最大公约数（ Greatest Common Divisor ）。
定义：$GCD(a,b)=GCD(b,a\% b)$
证明：
- 若$a<b$ ，则计算一次后两者交换
- $a=k\times b+r\quad r=a\%b $
- 若$d$为$a,b$的公约数，则$a \% d=0,b\%d=0$，则$(a-k\times b)\%d=r\%d=0$，故$d$也是$r$的公约数，故$GCD(a,b)$和$GCD(b,r)$相等（为了便于迭代我们都把小的数写在右边）。
- 当$b=0$时，$GCD(a,b)=a$
- 故只需要迭代上述过程至右侧数为0即可
int gcd(int a，int b) {
return b == 0 ? a : gcd(b, a % b);
}顺便附上求最小公倍数的方法：$LCM(a,b)=a\div GCD(a,b) \times b$
裴蜀定理（贝祖定理）
裴蜀定理（贝祖定理，Bézout’s identity ）得名于法国数学家艾蒂安·裴蜀。又称为贝祖定理和贝祖等式。
在数论中，裴蜀定理是一个关于最大公约数的定理，具体而言是一个关于未知数 x 和 y 的线性丢番图方程。
丢番图方程（Diophantine Equation）：有一个或者几个变量的整系数方程，它们的求解仅仅在整数范围内进行。
定义：若$a$，$b$是整数，则存在整数$x$，$y$使得$ax+by=GCD(a,b)$,且$GCD(a,b)$是$ax+by$集合中的最小正整数（若$a\ne 0$ 或 $b\ne 0$）。
表述2：若$a$，$b$是整数，$ax+by=c$有整数解的充要条件是$c\%GCD(a,b)=0$。
推论：$GCD(a,b)=1$的充要条件是存在整数解使$ax+by=1$。
该定理也可以推广到多个数的情况。这个定理仔细想想会觉得非常有道理，因为$abs(a-b)=k\times GCD(a,b)$必然成立，我们不可能组合出更小的数了，实在是找不到反例，但是证明比较复杂，就先跳过吧。
扩展欧几里得定理
扩展欧几里得定理（Extended Euclidean algorithm）可以用来求上述定理表述2中的整数$x$，$y$。
我们已经知道有解的充要条件，在该条件满足的情况下，我们只需要求$ax+by=GCD(a,b)$的解即可。
首先，我们考虑边界情况 ，当$b=0$时，原方程变为$ax=GCD(a,0)$，则
$$
\begin{cases} x=1\\ y=0 \end{cases}\tag{1}
$$
这个边界条件就是算法的终止条件，现在，根据欧几里得定理，$GCD(a,b)=GCD(b,a\% b)$，
我们设$a’=b,b’=a \% b$，则$a′x′+b′y′=GCD(a′,b′)=GCD(b,a \% b)=GCD(a,b)$，
交换一下符号，
$$ax+by=a′x′+b′y′=bx′+(a\%b)y′=bx′+(a−\lfloor\frac{a}{b}\rfloor b)y′=GCD(a,b)$$
整理可得：
$$
ax+by=ay′+b(x′−\lfloor\frac{a}{b}\rfloor y′)
$$
对比系数，可得：
$$
\begin{cases} x=y’ \\ y=x’−\lfloor\frac{a}{b}\rfloor y′ \end{cases}\tag{2}
$$
这里是说$x′,y′$分别是$a’,b’$对应的解，更新以后我们得到了更上一层的解，对比更新和终点这两个式子，我们就可以得到方程$ax+by=GCD(a,b)$的迭代解法，先递归到最后一层，再逐层往回迭代：
int exgcd(int a, int b, int &x, int &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
//r is GCD(a,b)
int r = exgcd(b, a%b, x, y);
int t = x; x = y; y = t - a/b*y;
return r;
}中国剩余定理
中国剩余定理/中国同余定理（Chinese Remainder Theorem） 来源于《孙子算经》 ，是数论的重要内容 。
今有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二，问物几何？
答曰：二十三。
术曰：三三数之剩二，置一百四十；五五数之剩三，置六十三，七七数之剩二，置三十，并之。得二百三十三，以二百一十减之，即得。凡三三数之剩一，则置七十；五五数之剩一，则置二十一；七七数之剩一，则置十五；一百六以上以一百五减之即得。
使用数学语言描述如下：$x \equiv a(mod\ m)$成为两个数模$m$同余，若$a<m$，则可以认为$x\%m=a $
问题：已知同余方程组
$$
\begin{cases} x \equiv 2(mod\ 3) \\ x \equiv 3(mod\ 5) \\ x \equiv 2(mod\ 7) \end{cases}
$$
求满足条件的$x$
求解：假设$m_1, m_2, …, m_n$两两互质，则对于任意整数$a_1, a_2, …, a_n$，同余方程组
$$
\begin{cases} x \equiv a_1(mod\ m_1) \\ x \equiv a_2(mod\ m_2) \\ … \\ x \equiv a_n(mod\ m_n) \end{cases}
$$
存在整数解，且解$X,Y$必满足$X \equiv Y (mod\ N)$，其中$N=\prod^n_{i=1}m_i$。
令$N_i=\frac{N}{m_i}$，$t_i$是方程$N_it_i\equiv 1 (\mod m_i)$的解，则$x$满足$ x \equiv \sum^n_{i=1} a_i N_i t_i (mod\ N) $。
证明：
对$\forall k \neq i$有$a_iN_it_i\equiv 0(\mod m_k)$，$a_iN_it_i\equiv a_i (\mod m_i)$，可以认为$a_iN_it_i$只对方程$i$有贡献，对其他方程没有影响。
具体解法：$t_i$可以通过扩展欧几里得来求解。我们构造如下贝组等式：
$N_it_i+m_iy=1$，显然$N_i$和$m_i$必然互质，则我们构造了形如$ax+by=GCD(a,b)$的方程，可以用扩欧来求解了。
顺带一提，此处的$t_i$其实就是$N_i$的模$m_i$乘法逆元。模逆元也叫模反元素，通常在算法中可以用来计算模数下的大数除法。详情参见维基百科模反元素。
扩展中国剩余定理
实际上并不存在这个定理，这个问题是CRT针对扩展欧几里得的构造扩展问题。
若$m_1, m_2, …, m_n$不满足两两互质，则针对不互质的两个数，可以将
$$
\begin{cases} x \equiv a_1(mod\ m_1) \\ x \equiv a_2(mod\ m_2) \\\end{cases}
$$
变形为
$$
\begin{cases} x =a_1+m_1y_1 \\ x =a_2+m_2y_2 \\\end{cases}
$$
即$m_1y_1-m_2y_2=a_2-a_1$，这是个形如扩展欧几里得的式子，如果无解（根据裴蜀定理判断），则原方程组无解，否则可构造一个可行解$x\equiv x_1 (\mod m)$，其中$m=LCM(m_1,m_2)$。
原方程组则等价于$x\equiv x_1 (\mod m)$，这样就消去了一个方程，和剩下的方程联立直到满足$m_1, m_2, …, m_n$满足两两互质即可。
多项式概念扩展
可以将整数的欧几里得除法定义扩展到多项式，用长除法来表示：
$$
\begin{align} \frac{x^4 + x^3 - x^2 + 3x + 2}{x^3 + 3x^2 + 3x + 2} &= \frac{-2x^3 - 4x^2 + x + 2}{x^3 + 3x^2 + 3x + 2} + (x) \\ &= \frac{2x^2 + 7x + 6}{x^3 + 3x^2 + 3x + 2} + (x-2) \end{align}
$$
$q(x)$为除数，$r(x)$为余数， 长除法除到最后余数的最高次幂小于除数的最高次幂即可。
$$
{x^4 + x^3 - x^2 + 3x + 2} = (x-2)\,\underbrace{(x^3 + 3x^2 + 3x + 2)}_{q(x)} +\underbrace{(2x^2 + 7x + 6)}_{r(x)}
$$
可以按照同样的方式来定义GCD：
$$
x^2 + 7x + 6 = (x + 1)(x + 6)\
x^2 − 5x − 6 = (x + 1)(x − 6)
$$
则两多项式的GCD为$(x+1)$，此时的互质是说两个多项式除零次多项式（常数）外不再有其他的公因式 。
类似的，贝祖定理、扩展欧几里得定理，中国剩余定理的解$ x \equiv \sum^n_{i=1} a_i N_i t_i (mod\ N) $也可以被扩展到多项式：
$$
CRT(x)\equiv \left(\sum_{i=0}^na^{(i)}(x)N^{(i)}(x)t^{(i)}(x)\right)(mod\ {N(x)})
$$
在多项式中有一个度（degree）的概念，即多项式中包含的最高次项的次数。
卷积与多项式乘法
拉格朗日多项式插值法
Cook-Toom Algorithm
Winograd Algorithm for $F(I_{2\cdot 2},K_{3 \cdot 3})$
Winograd Algorithm for $F(I_{4\cdot 4},K_{3 \cdot 3})$
实验
总流程
我们计算的总流程可以用下图（来源：2019 GMTC 全球大前端技术大会，阿里MNN组）表示，其中所有矩阵中的参数都是提前计算好的，在我们拿到卷积核的参数时可以先对其进行预处理，针对数据的变换和输出的变换则实时计算：
计算优化
在扩展到$F(I_{4\cdot 4},K_{3 \cdot 3})$时，我们并没有给出具体的计算方式，只给了矩阵的参数。和$F(I_{2\cdot 2},K_{3 \cdot 3})$类似的，其中很多计算都可以被复用。我们观察矩阵，给相关联的元素标上色：
括号中代表了所需的最少加法次数和乘法次数。原文中的表述是Floating Point Instruction，猜测是代表乘加指令，即$a \leftarrow a+b \times c$表示一个指令。
Input transform（12I，6M12A）:
$$
\begin{align} d’_0 &= 4d_0 - 5d_2 + d_4 \\ d’_1 &= -(4d_1 - d_3) - (4d_2 - d_4) \\ d’_2 &= (4d_1 - d_3) - (4d_2 - d_4) \\ d’_3 &= -2(d_1 - d_3) - (d_2 - d_4) \\ d’_4 &= 2(d_1 - d_3) - (d_2 - d_4) \\ d’_5 &= 4d_1 - 5d_3 + d_5 \end{align}
$$
Filter transform（8I，6M6A2S）:
$$
\begin{align} g’_0 &= \frac{1}{4}g_0 \\ g’_1 &= -\frac{1}{6}(g_0 + g_2) - \frac{1}{6}g_1 \\ g’_2 &= -\frac{1}{6}(g_0 + g_2) + \frac{1}{6}g_1 \\ g’_3 &= (\frac{1}{24}g_0 + \frac{1}{6}g_2) + \frac{1}{12}g_1 \\ g’_4 &= (\frac{1}{24}g_0 + \frac{1}{6}g_2) - \frac{1}{12}g_1 \\ g’_5 &= g_2 \end{align}
$$
EWMM（6M0A）:
$$
m_i=d’_ig’_i
$$
Output transform（10I，3M10A）:
$$
\begin{align} y_0 &= (m_1 + m_2) + (m_3 + m_4) + m_0 \\ y_1 &= (m_1 - m_2) + 2(m_3 - m_4) \\ y_2 &= (m_1 + m_2) + 4(m_3 + m_4) \\ y_3 &= (m_1 - m_2) + 8(m_3 - m_4) + m_5 \end{align}
$$
扩展到二维，其中Input transform（$B^TdB$）包括$144=12\times (6+6)$次浮点数操作（包括乘法），Filter transform（$GgG^T$）包括$72=8\times (3+6)$次浮点数操作（预先计算），Output transform（$A^TMA$）（Inverse transform）包括$100=10\times (6+4)$次浮点数操作。其中括号里面的数表示了变换矩阵的两个维度。
算法复杂度分析
Winograd算法的点乘复杂度如下，其中假设了输入输出的大小一致，均为$ H \times W $，输入通道数为$C$、卷积核组数/输出通道数为$K$，批大小为$N$：
$$
X=N\lceil H / m\rceil\lceil W / n\rceil C K(m+R-1)(n+S-1)
$$
当tile大小为$1\times 1$时，可以看出Winograd与普通卷积无异，即与标准的im2col做法没有区别。
简明起见，我们假设图像尺寸是tile大小的整数倍，并且卷积核和tile都是正方形。
我们使用$\alpha,\beta,\gamma,\delta$来分别表示点乘的乘法复杂度、输入变换（Data）、权重变换（Filter）、输出变换（Inverse）的运算复杂度，
则点乘的乘法复杂度为：
$$
\begin{align} X &= \frac{\alpha}{m^2} NHWCK \\ &= \alpha’ NHWCK \end{align}
$$
变换过程（Data、Filter、Inverse）的运算复杂度为：
$$
\begin{align} T(D) &= \frac{\beta}{m^2} NHWC \\ T(F) &= \gamma CK \\ T(I) &= \frac{\delta}{m^2} NHWK \end{align}
$$
将以上数值都进行归一化，保持和$X$一致的量纲，可以得到整个卷积过程的算法复杂度如下：
$$
\begin{align} L &=X+T(D)+T(F)+T(I)\\&= (1 + \frac{T(D)}{X} + \frac{T(F)}{X} + \frac{T(I)}{X}) X \\ &= (1 + \frac{\beta’}{K} + \frac{\gamma’}{P} + \frac{\delta’}{C}) \alpha’ NHWCK \end{align}
$$
其中$P$表示每个Channel的tile数量，
$$
\begin{align} \alpha’ &= \frac{\alpha}{m^2} = \frac{(m+r-1)^2}{m^2} \\ \beta’ &= \frac{\beta}{(m+r-1)^2} = \frac{\beta}{\alpha} \\ \gamma’ &= \frac{\gamma}{(m+r-1)^2} = \frac{\gamma}{\alpha} \\ \delta’ &= \frac{\delta}{(m+r-1)^2} = \frac{\delta}{\alpha} \\P&=\frac{NHW}{m^2} \end{align}
$$
- $T(F)$只需要最初的时候做一次，通常也不考虑进卷积计算的复杂度
- $T(D)$只对卷积输入做一次，其计算复杂度会被输出通道数量$K$平摊
- $T(I)$只对卷积输出做一次，其计算复杂度会被输入通道数量$C$平摊
- 显然，随着$C、K$增大，$T(D)、T(I)$对整体运算复杂度的影响就会减小，Winograd带来的加速效果越明显；
- 同理，当$C、K$太小时，$T(D)、T(I)$对整体运算复杂度的影响会大于Winograd带来的加速效果，其计算效率反而不如直接卷积
tile大小与各个参数的关系如下，其中，直接卷积可以看成tile大小为3：
内存占用与加速比
$F(I_{2\cdot 2},K_{3 \cdot 3})$的加速比上限是2.25，$F(I_{4\cdot 4},K_{3 \cdot 3})$的加速比上限是4，$F(I_{6\cdot 6},K_{3 \cdot 3})$的加速比上限是5.06。
常见的Winograd实现与普通卷积的效果比较如下：
ARM在Embedded Vision Summit 2018上的Slides中指出了一些常见模型上的加速情况，实现的流程和我们之前说的一致：
不过ARM Compute Library的优化做的还不够好，所以真实情况下的加速比还有可能做的更高：
实现
矩阵参数获取
CVPR论文作者提供了一份生成矩阵参数的代码，可以参考这里。
端侧推理
以$F(I_{2\cdot 2},K_{3 \cdot 3})$为例，我们来看看代码是如何实现的，代码来自Tencent NCNN。通常针对特定的平台会使用硬件加速指令，此处的代码只是单纯的C++实现。可以看到一个很重要的优化手段就是for-loop展开，并且在真正计算的部分，中间结果我们尽量不使用二维数组，提高访问效率。另外，就像我们之前提到的，部分计算可以被重用。
- Filter Transform
//3x3卷积，stride=1，F(2,3)
static void conv3x3s1_winograd23_transform_kernel(const Mat& kernel, Mat& kernel_tm, int inch, int outch)
{
kernel_tm.create(4*4, inch, outch);
// G
const float ktm[4][3] = {
{ 1.0f, 0.0f, 0.0f},
{ 1.0f/2, 1.0f/2, 1.0f/2},
{ 1.0f/2, -1.0f/2, 1.0f/2},
{ 0.0f, 0.0f, 1.0f}
};
#pragma omp parallel for
for (int p = 0; p<outch; p++)
{
for (int q = 0; q<inch; q++)
{
const float* kernel0 = (const float*)kernel + p*inch * 9 + q * 9;
//每一行都有一个U
float* kernel_tm0 = kernel_tm.channel(p).row(q);
// transform kernel
const float* k0 = kernel0;
const float* k1 = kernel0 + 3;
const float* k2 = kernel0 + 6;
// h
float tmp[4][3]; // tmp = (G*g)T
for (int i=0; i<4; 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];
}
// U
// kernel_tm0 = G*g*GT = U
for (int j=0; j<4; j++)
{
float* tmpp = &tmp[j][0];
for (int i=0; i<4; i++)
{
kernel_tm0[j*4 + i] = tmpp[0] * ktm[i][0] + tmpp[1] * ktm[i][1] + tmpp[2] * ktm[i][2];
}
}
}
}
}
- Input transform
static void conv3x3s1_winograd23(const Mat& bottom_blob, Mat& top_blob, const Mat& kernel_tm, const Mat& _bias, const Option& opt)
{
int w = bottom_blob.w;
int h = bottom_blob.h;
int inch = bottom_blob.c;
int outw = top_blob.w;
int outh = top_blob.h;
int outch = top_blob.c;
// pad to 2n+2, winograd F(2,3)
Mat bottom_blob_bordered = bottom_blob;
outw = (outw + 1) / 2 * 2;
outh = (outh + 1) / 2 * 2;
w = outw + 2;
h = outh + 2;
Option opt_b = opt;
opt_b.blob_allocator = opt.workspace_allocator;
copy_make_border(bottom_blob, bottom_blob_bordered, 0, h - bottom_blob.h, 0, w - bottom_blob.w, 0, 0.f, opt_b);
const float* bias = _bias;
// BEGIN transform input
Mat bottom_blob_tm;
{
//total size 4x4->2x2
int w_tm = outw / 2 * 4;
int h_tm = outh / 2 * 4;
//tile number
int nColBlocks = h_tm/4; // may be the block num in Feathercnn
int nRowBlocks = w_tm/4;
const int tiles = nColBlocks * nRowBlocks;
bottom_blob_tm.create(4*4, tiles, inch, 4u, opt.workspace_allocator);
// BT
// const float itm[4][4] = {
// {1.0f, 0.0f, -1.0f, 0.0f},
// {0.0f, 1.0f, 1.00f, 0.0f},
// {0.0f, -1.0f, 1.00f, 0.0f},
// {0.0f, -1.0f, 0.00f, 1.0f}
// };
#pragma omp parallel for num_threads(opt.num_threads)
for (int q=0; q<inch; q++)
{
const float* img = bottom_blob_bordered.channel(q);
float* out_tm0 = bottom_blob_tm.channel(q);
for (int j = 0; j < nColBlocks; j++)
{
const float* r0 = img + w * j * 2;
const float* r1 = r0 + w;
const float* r2 = r1 + w;
const float* r3 = r2 + w;
for (int i = 0; i < nRowBlocks; i++)
{
float d0[4],d1[4],d2[4],d3[4];
float w0[4],w1[4],w2[4],w3[4];
float t0[4],t1[4],t2[4],t3[4];
// load
for (int n = 0; n < 4; n++)
{
d0[n] = r0[n];
d1[n] = r1[n];
d2[n] = r2[n];
d3[n] = r3[n];
}
// w = B_t * d
for (int n = 0; n < 4; n++)
{
w0[n] = d0[n] - d2[n];
w1[n] = d1[n] + d2[n];
w2[n] = d2[n] - d1[n];
w3[n] = d3[n] - d1[n];
}
// transpose d to d_t
{
t0[0]=w0[0]; t1[0]=w0[1]; t2[0]=w0[2]; t3[0]=w0[3];
t0[1]=w1[0]; t1[1]=w1[1]; t2[1]=w1[2]; t3[1]=w1[3];
t0[2]=w2[0]; t1[2]=w2[1]; t2[2]=w2[2]; t3[2]=w2[3];
t0[3]=w3[0]; t1[3]=w3[1]; t2[3]=w3[2]; t3[3]=w3[3];
}
// d = B_t * d_t
for (int n = 0; n < 4; n++)
{
d0[n] = t0[n] - t2[n];
d1[n] = t1[n] + t2[n];
d2[n] = t2[n] - t1[n];
d3[n] = t3[n] - t1[n];
}
// save to out_tm
for (int n = 0; n < 4; n++)
{
out_tm0[n ] = d0[n];
out_tm0[n+ 4] = d1[n];
out_tm0[n+ 8] = d2[n];
out_tm0[n+12] = d3[n];
}
r0 += 2;
r1 += 2;
r2 += 2;
r3 += 2;
out_tm0 += 16;
}
}
}
}- Hadamard product
bottom_blob_bordered = Mat();
// BEGIN dot
Mat top_blob_tm;
{
int w_tm = outw / 2 * 4;
int h_tm = outh / 2 * 4;
int nColBlocks = h_tm/4; // may be the block num in Feathercnn
int nRowBlocks = w_tm/4;
const int tiles = nColBlocks * nRowBlocks;
top_blob_tm.create(16, tiles, outch, 4u, opt.workspace_allocator);
int nn_outch = outch >> 2;
int remain_outch_start = nn_outch << 2;
#pragma omp parallel for num_threads(opt.num_threads)
for (int pp=0; pp<nn_outch; pp++)
{
int p = pp * 4;
Mat out0_tm = top_blob_tm.channel(p);
Mat out1_tm = top_blob_tm.channel(p+1);
Mat out2_tm = top_blob_tm.channel(p+2);
Mat out3_tm = top_blob_tm.channel(p+3);
const Mat kernel0_tm = kernel_tm.channel(p);
const Mat kernel1_tm = kernel_tm.channel(p+1);
const Mat kernel2_tm = kernel_tm.channel(p+2);
const Mat kernel3_tm = kernel_tm.channel(p+3);
for (int i=0; i<tiles; i++)
{
float* output0_tm = out0_tm.row(i);
float* output1_tm = out1_tm.row(i);
float* output2_tm = out2_tm.row(i);
float* output3_tm = out3_tm.row(i);
float sum0[16] = {0.0f};
float sum1[16] = {0.0f};
float sum2[16] = {0.0f};
float sum3[16] = {0.0f};
int q = 0;
for (; q+3<inch; q+=4)
{
const float* r0 = bottom_blob_tm.channel(q).row(i);
const float* r1 = bottom_blob_tm.channel(q+1).row(i);
const float* r2 = bottom_blob_tm.channel(q+2).row(i);
const float* r3 = bottom_blob_tm.channel(q+3).row(i);
const float* k0 = kernel0_tm.row(q);
const float* k1 = kernel1_tm.row(q);
const float* k2 = kernel2_tm.row(q);
const float* k3 = kernel3_tm.row(q);
for (int n=0; n<16; n++)
{
sum0[n] += r0[n] * k0[n];
k0 += 16;
sum0[n] += r1[n] * k0[n];
k0 += 16;
sum0[n] += r2[n] * k0[n];
k0 += 16;
sum0[n] += r3[n] * k0[n];
k0 -= 16 * 3;
sum1[n] += r0[n] * k1[n];
k1 += 16;
sum1[n] += r1[n] * k1[n];
k1 += 16;
sum1[n] += r2[n] * k1[n];
k1 += 16;
sum1[n] += r3[n] * k1[n];
k1 -= 16 * 3;
sum2[n] += r0[n] * k2[n];
k2 += 16;
sum2[n] += r1[n] * k2[n];
k2 += 16;
sum2[n] += r2[n] * k2[n];
k2 += 16;
sum2[n] += r3[n] * k2[n];
k2 -= 16 * 3;
sum3[n] += r0[n] * k3[n];
k3 += 16;
sum3[n] += r1[n] * k3[n];
k3 += 16;
sum3[n] += r2[n] * k3[n];
k3 += 16;
sum3[n] += r3[n] * k3[n];
k3 -= 16 * 3;
}
}
for (; q<inch; q++)
{
const float* r0 = bottom_blob_tm.channel(q).row(i);
const float* k0 = kernel0_tm.row(q);
const float* k1 = kernel1_tm.row(q);
const float* k2 = kernel2_tm.row(q);
const float* k3 = kernel3_tm.row(q);
for (int n=0; n<16; n++)
{
sum0[n] += r0[n] * k0[n];
sum1[n] += r0[n] * k1[n];
sum2[n] += r0[n] * k2[n];
sum3[n] += r0[n] * k3[n];
}
}
for (int n=0; n<16; n++)
{
output0_tm[n] = sum0[n];
output1_tm[n] = sum1[n];
output2_tm[n] = sum2[n];
output3_tm[n] = sum3[n];
}
}
}
#pragma omp parallel for num_threads(opt.num_threads)
for (int p=remain_outch_start; p<outch; p++)
{
Mat out0_tm = top_blob_tm.channel(p);
const Mat kernel0_tm = kernel_tm.channel(p);
for (int i=0; i<tiles; i++)
{
float* output0_tm = out0_tm.row(i);
float sum0[16] = {0.0f};
int q = 0;
for (; q+3<inch; q+=4)
{
const float* r0 = bottom_blob_tm.channel(q).row(i);
const float* r1 = bottom_blob_tm.channel(q+1).row(i);
const float* r2 = bottom_blob_tm.channel(q+2).row(i);
const float* r3 = bottom_blob_tm.channel(q+3).row(i);
const float* k0 = kernel0_tm.row(q);
const float* k1 = kernel0_tm.row(q+1);
const float* k2 = kernel0_tm.row(q+2);
const float* k3 = kernel0_tm.row(q+3);
for (int n=0; n<16; n++)
{
sum0[n] += r0[n] * k0[n];
sum0[n] += r1[n] * k1[n];
sum0[n] += r2[n] * k2[n];
sum0[n] += r3[n] * k3[n];
}
}
for (; q<inch; q++)
{
const float* r0 = bottom_blob_tm.channel(q).row(i);
const float* k0 = kernel0_tm.row(q);
for (int n=0; n<16; n++)
{
sum0[n] += r0[n] * k0[n];
}
}
for (int n=0; n<16; n++)
{
output0_tm[n] = sum0[n];
}
}
}
}
bottom_blob_tm = Mat();
// END dot- Output transform
// BEGIN transform output
Mat top_blob_bordered;
top_blob_bordered.create(outw, outh, outch, 4u, opt.workspace_allocator);
{
// AT
// const float itm[2][4] = {
// {1.0f, 1.0f, 1.0f, 0.0f},
// {0.0f, 1.0f, -1.0f, 1.0f}
// };
int w_tm = outw / 2 * 4;
int h_tm = outh / 2 * 4;
int nColBlocks = h_tm/4; // may be the block num in Feathercnn
int nRowBlocks = w_tm/4;
#pragma omp parallel for num_threads(opt.num_threads)
for (int p=0; p<outch; p++)
{
Mat out_tm = top_blob_tm.channel(p);
Mat out = top_blob_bordered.channel(p);
const float bias0 = bias ? bias[p] : 0.f;
for (int j=0; j<nColBlocks; j++)
{
float* outRow0 = out.row(j*2);
float* outRow1 = out.row(j*2+1);
for(int i=0; i<nRowBlocks; i++)
{
float* out_tile = out_tm.row(j*nRowBlocks + i);
float s0[4],s1[4],s2[4],s3[4];
float w0[4],w1[4];
float d0[2],d1[2],d2[2],d3[2];
float o0[2],o1[2];
// load
for (int n = 0; n < 4; n++)
{
s0[n] = out_tile[n];
s1[n] = out_tile[n+ 4];
s2[n] = out_tile[n+ 8];
s3[n] = out_tile[n+12];
}
// w = A_T * W
for (int n = 0; n < 4; n++)
{
w0[n] = s0[n] + s1[n] + s2[n];
w1[n] = s1[n] - s2[n] + s3[n];
}
// transpose w to w_t
{
d0[0] = w0[0]; d0[1] = w1[0];
d1[0] = w0[1]; d1[1] = w1[1];
d2[0] = w0[2]; d2[1] = w1[2];
d3[0] = w0[3]; d3[1] = w1[3];
}
// Y = A_T * w_t
for (int n = 0; n < 2; n++)
{
o0[n] = d0[n] + d1[n] + d2[n] + bias0;
o1[n] = d1[n] - d2[n] + d3[n] + bias0;
}
// save to top blob tm
outRow0[0] = o0[0];
outRow0[1] = o0[1];
outRow1[0] = o1[0];
outRow1[1] = o1[1];
outRow0 += 2;
outRow1 += 2;
}
}
}
}
// END transform output
// cut result pad
copy_cut_border(top_blob_bordered, top_blob, 0, top_blob_bordered.h - top_blob.h, 0, top_blob_bordered.w - top_blob.w, opt);
}$F(I_{4\cdot 4},K_{3 \cdot 3})$的实现就更加复杂，可以看参考资料中Tencent NCNN的实现。
跑题： NHWC、HCHW和NC/4HW4 -TODO
我们在说NHWC时，实际上是说数据在内存上的排布策略，以及探讨该策略带来的访存性能变化。
NHWC与HCHW的区别用一张图就可以说清楚，其中我们的N=1，C=3，H=1，W=6：
NHWC 的访存局部性更好（每三个输入像素即可得到一个输出像素），NCHW 则必须等所有通道输入准备好才能得到最终输出结果，需要占用较大的临时空间。
在 CNN 中常常见到 1x1 卷积（例如：用于移动和嵌入式视觉应用的 MobileNets），也是每个输入 channel 乘一个权值，然后将所有 channel 结果累加得到一个输出 channel。如果使用 NHWC 数据格式，可以将卷积计算简化为矩阵乘计算，即 1x1 卷积核实现了每个输入像素组到每个输出像素组的线性变换。
TensorFlow 为什么选择 NHWC 格式作为默认格式？因为早期开发都是基于 CPU，使用 NHWC 比 NCHW 稍快一些（不难理解，NHWC 局部性更好，cache 利用率高）。
NCHW 则是 Nvidia cuDNN 默认格式，使用 GPU 加速时用 NCHW 格式速度会更快（也有个别情况例外）。
最早接触到NC/4HW4是因为阿里的端侧推理框架MNN，不得不说这个框架真是太复杂了，随便拿一个技术点就能扯上半天。
NC/4HW4 的布局是为了更好地利用 SIMD 加速。NC/4HW4第一个4表示把原Feature map的通道按4分组不够补0，然后每组内的4个Feature map按照RGBA交织排列。
总结
- Winograd算法存在一定的计算精度损失。不过CNN模型需要的计算精度实际上很低，例如有用fp16、int8实现CNN的方法，也有用更低bit数甚至binary计算实现的方法，它们都有不错的ImageNet分类精度。
tile越大，加速比越大，卷积核大小也类似，但所需的加法、transform和存储的代价，计算精度的损失都会迅速增大。 因此通常只会采用较小的tile和卷积核大小，最常见的实现为$F(I_{2\cdot 2},K_{3 \cdot 3})$，$F(I_{4\cdot 4},K_{3 \cdot 3})$，$F(I_{6\cdot 6},K_{3 \cdot 3})$。- Winograd算法可以用矩阵形式来表示，但是具体实现时，并不意味着要调用矩阵运算的接口，为了更快的计算速度，通常会直接将计算展开，故代码量较大，且对于不同的
tile大小需要专门定制的代码（好在也就那么几种），通常卷积核的大小为$2\times 2$到$7\times 7$。在划分方式上，$F(I_{1\cdot 1},K_{r \cdot r})$和im2col没有本质区别。
参考资料
- 卷积神经网络中的Winograd快速卷积算法
- Winograd卷积原理
- 源于《孙子算经》的Cudnn
- Paper-CVPR 2016-Fast Algorithms for Convolutional Neural Networks
- Book-Shmuel Winograd.1980 Arithmetic complexity of computations
- 开源背后 | 面对端侧推理引擎的挑战，阿里工程师如何应对
- Tencent FeatherCNN实现
- Tencent NCNN 3x3实现
- 深度学习轻量级推理及加速
- 【移动端DL框架】当前主流的移动端深度学习框架一览
- Video-Fast Algorithms for Convolutional Neural Networks by Andrew Lavin and Scott Gray
- Video-Even Faster CNNs Exploring the New Class of Winograd Algorithms
- PPT-算法解析
- Paper-Sparse Winograd Convolutional neural networks on small-scale systolic arrays
- Code-参数生成 andravin/wincnn
- PPT-Fast Convolution
- NCNN winograd详解（一）
- 深度学习（十一）——Winograd（1）
![[ML]-统计学习-1（卡尔曼滤波与粒子滤波）](/medias/featureimages/7.jpg)
