现代深度学习库对大多数操作都具有生产级的、高度优化的实现，这并不奇怪。但这些库究竟是什么魔法？他们如何能够将性能提高100倍？究竟怎样才能“优化”或加速神经网络的运行呢？在讨论高性能/高效DNNs时，我经常会问(也经常被问到)这些问题。
在这篇文章中，我将尝试带你了解在DNN库中卷积层是如何实现的。它不仅是在模型中最常见的和最重的操作，我还发现卷积高性能实现的技巧特别具有代表性——一点点算法的小聪明，非常多的仔细的调优和低层架构的开发。
我在这里介绍的很多内容都来自Goto等人的开创性论文：Anatomy of a high-performance matrix multiplication，该论文为OpenBLAS等线性代数库中使用的算法奠定了基础。
最原始的卷积实现
“过早的优化是万恶之源”——Donald Knuth
在进行优化之前，我们先了解一下基线和瓶颈。这是一个朴素的numpy/for循环卷积：
'''
Convolve `input` with `kernel` to generate `output`
input.shape = [input_channels, input_height, input_width]
kernel.shape = [num_filters, input_channels, kernel_height, kernel_width]
output.shape = [num_filters, output_height, output_width]
'''
for filter in 0..num_filters
for channel in 0..input_channels
for out_h in 0..output_height
for out_w in 0..output_width
for k_h in 0..kernel_height
for k_w in 0..kernel_width
output[filter, channel, out_h, out_h] +=
kernel[filter, channel, k_h, k_w] *
input[channel, out_h + k_h, out_w + k_w]这是一个6层嵌套的for循环(如果你迭代多个输入批次，则为7个)。我们还没有研究步幅，膨胀，或其他参数。如果我在这里输入MobileNet的第一层的大小，然后用普通的C语言运行它，它会花费惊人的22秒！使用最积极的编译器优化，如' -O3 '或' -Ofast '，它减少到2.2秒。但这对于第一层来说仍然非常慢。
如果我使用Caffe运行相同的层呢？这台电脑只用了18毫秒。这比100倍的加速还要快！整个网络在我的CPU上运行大约100毫秒。
瓶颈是什么，我们应该从哪里开始优化?
最内部的循环执行两个浮点运算(乘法和加法)，对于我使用的大小，它执行了大约8516万次，也就是说，这个卷积需要1.7亿个浮点运算(MFLOPs)。根据英特尔的数据，我的CPU的最高性能是每秒800亿次，也就是说，理论上它可以在0.002秒内完成这项工作。显然，我们离这个目标还很远，而且很明显，这里的原始处理能力已经足够了。
理论峰值没有达到(从来没有)的原因是内存访问也需要时间—如果不能快速获得数据，那么仅仅快速处理数据是不够的。事实证明，上面嵌套的for循环使得数据访问模式非常困难，这使得缓存利用率很低。
正如你将看到的，在整个讨论过程中反复出现的一个问题是，我们如何访问正在操作的数据，以及这些数据如何与存储方式相关联。
卷积底层实现
FLOP/s
我们对“性能”或速度的度量是吞吐量，以每秒浮点计算次数度量。具有更多浮点操作的更大操作自然会运行得更慢，因此FLOP/s速率可以使用更一致的方式来比较性能。
我们也可以用它来了解我们离CPU的最高性能有多近。在我的笔记本电脑CPU上：
- 有2个phsyical core
- 每个核的频率为2.5 GHz，即每秒2.5×10^9个CPU周期
- 在每个周期，它可以处理32FLOPs(使用AVX和FMA还会更多)
其峰值性能为 GFLOP/s。这是我的CPU的理论峰值。同样，对于单个内核，这个数字是80GFLOP/s。
存储顺序和行主序
虽然我们从逻辑上把矩阵/图像/张量看作多维的，但它们实际上存储在线性的一维计算机内存中。我们必须定义一个约定，该约定规定如何将这些多维数据展开到线性存储中，反之亦然。
大多数现代DL库使用行主序存储。这意味着同一行的连续元素彼此相邻存储。更一般地说，对于多维，行主序意味着当线性扫描内存时，第一个维度的变化最慢。
那么维度本身的顺序呢？通常对于四维张量，比如CNNs，你会听到NCHW, NHWC等存储顺序。我将在这篇文章中假设NCHW——如果我有N块HxW图像的C通道，那么所有具有相同N个通道的图像都是重叠的，在该块中，同一通道C的所有像素都是重叠的，以此类推。
Halide
这里讨论的许多优化都需要在底层使用神秘的C语法，甚至是程序集进行干预。这不仅使代码难以阅读，还使尝试不同的优化变得困难，因为我们必须重新编写整个代码。Halide是c++中的一种嵌入式语言，它帮助抽象这些概念，并被设计用来帮助编写快速图像处理代码。通过分解算法(要计算什么)和计划(如何/何时计算)，可以更容易地试验不同的优化。我们可以保持算法不变，并使用不同的策略。
我将使用Halide来表示这些较低级别的概念，但是你应该能够理解足够直观的函数名，以便理解。
从卷积到GEMM
我们上面讨论的简单卷积已经很慢了，一个更实际的实现只会因为步长、膨胀、填充等参数而变得更加复杂。我们可以继续使用基本的卷积作为一个工作示例，但是，正如你看到的，从计算机中提取最大性能需要许多技巧—在多个层次上进行仔细的微调并充分利用现有计算机体系结构的非常具体的知识。换句话说，如果我们希望解决所有的复杂性，这将是一项艰巨的任务。
我们能不能把它转化成一个更容易解决的问题？也许矩阵乘法？
矩阵乘法，或matmul，或Generalized Matrix Multiplication (GEMM)是深度学习的核心。它用于全连接层、RNNs等，也可用于实现卷积。
毕竟，卷积是带有输入padding的滤波器的点积。如果我们把滤波器放到一个二维矩阵中，把输入的小patch放到另一个矩阵中，然后把这两个矩阵相乘，就会得到相同的点积。与CNNs不同，矩阵乘法在过去几十年里得到了大量的研究和优化，在许多科学领域都是一个关键问题。
上面将图像块放到一个矩阵中的操作称为im2col ，用于图像到列。我们将图像重新排列成矩阵的列，使每一列对应一个应用卷积滤波器的patch。
考虑这个普通的，直接的3x3卷积：
下面是与矩阵乘法相同的操作。正确的矩阵是im2col的结果——它必须通过复制原始图像中的像素来构造。左边的矩阵有conv权值，它们已经以这种方式存储在内存中。
注意，矩阵乘积直接给出了conv输出——不需要额外的“转换”到原始形式。
为了清晰起见，我将每个patch都单独显示在这里。然而，在现实中，不同的图像块之间往往存在一定的重叠，因此im2col会产生一定的内存重复。生成这个im2col缓冲区和膨胀的内存所花费的时间，必须通过GEMM实现的加速来抵消。
利用im2col，我们已经将卷积运算转化为矩阵乘法。我们现在可以插入更多通用的和流行的线性代数库，如OpenBLAS、Eigen等，利用几十年的优化和仔细的调优，有效地计算这个matmul。
如果我们要证明im2col转换所带来的额外工作和内存是合理的，那么我们需要一些非常严重的加速，所以让我们看看这些库是如何实现这一点的。这也很好地介绍了在系统级进行优化时的一些通用方法。
虽然以前有过不同形式的convolution-as-GEMM思想，但Caffe是第一个在CPU和GPU上对通用convs使用这种方法的深度学习库之一，并显示了较大的加速。这里：https://github.com/yangqing//wiki/convoluin-Caffe:a-memo一个非常有趣的阅读来自于Yanqing Jia本人(Caffe的创始人)关于这个决定的起源，以及关于“临时”解决方案的想法。
加速GEMM
原始的方法
在这篇文章的其余部分，我将假设GEMM被执行为
和之前一样，首先让我们对基本的，课本上的矩阵乘法进行计时：
for i in 0..M:
for j in 0..N:
for k in 0..K:
C[i, j] += A[i, k] * B[k, j]使用Halide:
Halide::Buffer C, A, B;
Halide::Var x, y;
C(x,y) += A(k, x) *= B(y, k); // loop bounds, dims, etc. are taken care of automatically最里面的一行执行2个浮点运算(乘法和加法)，并执行\(M∗N∗K\)次，因此这个GEMM的FLOPs是\(2MNK\)。
我们来测量一下它在不同矩阵大小下的性能：
我们的性能才刚刚达到顶峰的10% ！虽然我们将研究使计算更快的方法，但一个反复出现的主题是，如果我们不能快速获得数据，仅仅快速计算数据是不够的。由于内存对于较大的矩阵来说是一个越来越大的问题，因此性能会逐渐下降。你最后看到的急剧下降，表示当矩阵变得太大而无法放入缓存时，吞吐量突然下降—你可以看到系统阻塞。
缓存
RAM是一个大而慢的存储器。CPU缓存的速度要快几个数量级，但要小得多，因此正确使用它们至关重要。但是没有明确的指令说“加载数据以缓存”。它是一个由CPU自动管理的进程。
每次从主存中获取数据时，CPU都会自动将数据及其相邻的内存加载到缓存中，希望利用引用的局部性。
你应该注意的第一件事是我们访问数据的模式。我们在A上按行遍历，在B上按列遍历。
它们的存储也是行主序的，所以一旦找到A[i, k]，行中的下一个元素A[i, k+1]就已经缓存了。酷。但看看B会发生什么：
- 该列的下一个元素没有出现在缓存中—在获取数据的时候，我们得到一个cache miss和CPU stalls。
- 一旦数据被获取，缓存也被填充在同一行B的其他元素。我们实际上不会使用它们，所以它们很快就会被驱逐。经过几次迭代之后，当实际需要它们时，我们将再次获取它们。我们正在用不需要的值污染缓存。
我们需要重新设计循环来利用这种缓存能力。如果正在读取数据，我们不妨利用它。这是我们要做的第一个更改：循环重新排序。
让我们重新排序循环，从i,j,k到i,k,j：
for i in 0..M:
for k in 0..K:
for j in 0..N:我们的答案仍然是正确的，因为乘法/加法的顺序无关紧要。遍历顺序现在看起来是这样的：
这个简单的改变，只是重新排序了一下循环，给了一个相当快的加速：
Tiling
为了进一步改进重新排序，我们还需要考虑一个缓存问题。
对于A的每一行，我们循环遍历整个B。在B中每进行一步，我们将加载它的一些新列并从缓存中删除一些旧列。当我们到达A的下一行时，我们从第一列开始重新开始。我们重复地从缓存中添加和删除相同的数据，这叫做抖动。
如果我们所有的数据都能放入缓存，就不会发生抖动。如果我们使用更小的矩阵，他们就可以幸福地生活在一起，而不会被反复驱逐。谢天谢地，我们可以分解子矩阵上的矩阵乘法。计算一个C中的小的r×c块，只需要A中的r行和B中的C列。让我们把C分成6x16的小块。
C(x, y) += A(k, y) * B(x, k);
C.update()
.tile(x, y, xo, yo, xi, yi, 6, 16)
/*
in pseudocode:
for xo in 0..N/16:
for yo in 0..M/6:
for yi in 6:
for xi in 0..16:
for k in 0..K:
C(...) = ...
*/我们已经把x,y的维度分解成外部的xo,yo和内部的xi,yi。我们将努力为较小的6x16块(标记为xi,yi)优化一个微内核，并在所有块上运行该微内核(由xo,yo迭代)。
Vectorization & FMA
大多数现代cpu都支持SIMD，或者Single Instruction Multiple Data。顾名思义，SIMD可以在相同的CPU周期内对多个值同时执行相同的操作/指令(如add、multiply等)。如果我们可以一次运行4个数据点上的SIMD指令，那么就可以实现4倍的加速。
因此，当我们计算处理器的峰值速度时，我们“有点”作弊，而是参考了这种向量化的性能。这对于像向量这样的数据非常有用，我们必须对每个向量元素应用相同的指令。但是我们仍然需要设计内核来正确地利用这一点。
我们在计算峰值故障时使用的第二个“hack”是FMA-Fused Multiply-Add。虽然乘法和加法被算作两个独立的浮点运算，但它们是如此常见，以至于可以使用专用的硬件单元来“融合”它们，并将它们作为一条指令执行。使用它通常由编译器处理。
在Intel cpu上，我们可以使用SIMD(称为AVX & SSE)在一条指令中处理多达8个浮点数。编译器优化通常能够自己识别向量化的机会，但为了确保这一点，我们将亲自动手。
C.update()
.tile(x, y, xo, yo, xi, yi, 6, 16)
.reorder(xi, yi, k, xo, yo)
.vectorize(xi, 8)
*/*
in pseudocode:
for xo in 0..N/16:
for yo in 0..M/6:
for k in 0..K:
for yi in 6:
for vectorized xi in 0..16:
C(...) = ...
*/*
Threading
到目前为止，我们只使用了一个CPU内核。我们有多个可用的内核，每个内核可以同时物理地执行多个指令。一个程序可以把自己分成多个线程，每个线程可以运行在一个单独的内核上。
C.update()
.tile(x, y, xo, yo, xi, yi, 6, 16)
.reorder(xi, yi, k, xo, yo)
.vectorize(xi, 8)
.parallel(yo)
*/*
in pseudocode:
for xo in 0..N/16 in steps of 16:
for parallel yo in steps of 6:
for k in 0..K:
for yi in 6:
for vectorized xi in 0..16 in steps of 8:
C(...) = ...
*/*你可能会注意到，对于非常小的大小，性能实际上会下降，因为在较小的工作负载下，线程的工作时间更少，而彼此同步的时间更多。在线程方面还有很多其他类似的问题，它们本身可能需要进一步深入研究。
Unrolling
循环让我们避免了一遍又一遍地编写同一行的痛苦，同时引入了一些额外的工作，比如检查循环终止、更新循环计数器、指针算法等。相反，如果我们手工编写重复循环语句并展开循环，我们可以减少这种开销。例如，我们可以运行包含4个语句的2个迭代，而不是1个语句的8个迭代。
起初，我感到惊讶的是，如此反复无常、看似微不足道的管理费用，实际上却很重要。然而，尽管这些循环操作可能“便宜”，但它们肯定不是免费的。如果你还记得这里迭代的数量是以百万为单位的，那么每次迭代额外的2-3条指令的成本将很快增加。当循环开销变得相对较小时，好处会逐渐减少。
展开是另一个优化，现在几乎完全由编译器来处理，除了在像我们这样的微内核中，我们更喜欢控制。
C.update()
.tile(x, y, xo, yo, xi, yi, 6, 16)
.reorder(xi, yi, k, xo, yo)
.vectorize(xi, 8)
.unroll(xi)
.unroll(yi)
*/*
in pseudocode:
for xo in 0..N/16:
for parallel yo:
for k in 0..K:
C(xi:xi+8, yi+0)
C(xi:xi+8, yi+1)
...
C(xi:xi+8, yi+5)
C(xi+8:xi+16, yi+0)
C(xi+8:xi+16, yi+1)
...
C(xi+8:xi+16, yi+5)
*/*我们现在可以达到60GFLOP/s。
全部放到一起
上面的步骤涵盖了一些最常用的转换来提高性能。它们通常以不同的方式组合在一起，从而产生越来越复杂的策略来计算相同的任务。
这是来自Halide 的一个更仔细优化的策略。
matrix_mul(x, y) **+=** A(k, y) ***** B(x, k);
out(x, y) **=** matrix_mul(x, y);
out.tile(x, y, xi, yi, 24, 32)
.fuse(x, y, xy).parallel(xy)
.split(yi, yi, yii, 4)
.vectorize(xi, 8)
.unroll(xi)
.unroll(yii);
matrix_mul.compute_at(out, yi)
.vectorize(x, 8).unroll(y);
matrix_mul.update(0)
.reorder(x, y, k)
.vectorize(x, 8)
.unroll(x)
.unroll(y)
.unroll(k, 2);
总结一下，它是这样做的：
- 分裂成32x24的小块。进一步将每个tile分成8x24个子tile
- 在临时缓冲区(
matrix_mul)中计算8x24的matmul，使用类似的重新排序、向量化和展开 - 使用矢量化、展开等方法将临时的
matrix_mul复制回out。 - 在所有32x24块上并行化这个过程
最后，我们能够达到超过120GFLOPs的速度—相当接近160 GFLOPs的峰值性能，并且能够匹配OpenBLAS等生产级库。使用类似的im2col微调代码，然后是gemm，相同的卷积现在运行时间为~20ms。如果你对这方面的研究感兴趣，可以尝试使用你自己的不同策略—作为一名工程师，我总是听到关于缓存、性能等方面的说法，看到它们的真实效果可能是有益的和有趣的。
注意，这种im2col+gemm方法是大多数深度学习库中流行的通用方法。然而，定制化是关键——对于特定的常用大小、不同的体系结构(GPU)和不同的操作参数(如膨胀、分组等)，这些库可能会再次使用针对这些情况的类似技巧或假设进行更定制化的实现。这些微内核是经过反复试验和错误的高度迭代过程构建的。程序员通常只对什么应该/不应该工作得很好有一种直觉，并且/或者必须基于结果考虑解释。听起来很适合深度学习研究，对吧？
卷积的优势
传统神经网络都是采用矩阵乘法来建立输入和输出之间的关系，假如我们有 M 个输入和 N个输出，那么在训练过程中，我们需要 M×N 个参数去刻画输入和输出的关系 。当 M 和 N都很大，并且再加几层的卷积网络，这个参数量将会大的离谱。
那么对于这种问题，卷积网络怎么应对呢？？
卷积运算主要通过三个重要思想来改进上述面来的问题：稀疏连接、参数共享、平移不变性。 一一来看
稀疏连接
可能有些人不懂这种稀疏连接是怎么实现的？先来说说卷积操作，以一个二维矩阵为输入（可以看作是一个单通道图片的像素值），卷积产生的稀疏连接根本原因就是这块的核函数，一般的核函数的大小远小于输入的大小。
以下图例：卷积操作可以看做是一种滑窗法，首先，输入维度是4×4，输入中红色部分，先和核函数中的元素对应相乘，就是输出中左上角的元素值s1,即 \(s1 = a×k1+b×k2+e×k3+f×k4\)。
假设步长为2，那么输入中的红色部分将移动到下图的位置，继续卷积运算：\(s2 = c×k1+d×k2+g×k3+h×k4\)。
这么一直计算下去就完成了卷积操作，输入4×4的矩阵通过卷积变成为2×2的矩阵，其中，输出值s1只和输入中的a、b、e、f有关，和其他元素无关，刚好对应到前面所说的，卷积神经的特点就是每一层神经元只响应前一层的局部范围内的神经元。
继续用这个图来看稀疏连接，s2只和x1、x2、x3有关，感觉这个图看着更明显点！！！！
参数共享
参数共享是指在一个模型的多个函数中使用相同的参数，它是卷积运算带来的固有属性。
在全连接中，计算每层的输出时，权重矩阵中的元素只作用于某一个输入元素一次；
而在卷积神经网络中，卷积核中的每一个元素将作用于每一个局部输入的特定位置上。根据参数共享的思想，我们只需要学习一组参数集合，而不需要针对每一个位置的每一个参数来进行优化学习，从而大大降低了模型的存储需求。
平移不变性
平移不变性：如果一个函数的输入做了一些改变，那么输出也跟着做出同样的改变，这就时平移不变性。
平移不变性是由参数共享的物理意义所得。在计算机视觉中，假如要识别一个图片中是否有一只猫，那么无论这只猫在图片的什么位置，我们都应该识别出来，即就是神经网络的输出对于平移不变性来说是等变的。
一下图为例：s2只和x1、x2、x3有关，
输入层数据向左移动一位：其中响应层的s2也只和x1、x2、x3有关，只是位置发生了稍微的变化，
卷积与全连接的关系
在经典分类网络，比如LeNet、AlexNet中，在前面的卷积层提取特征之后都串联全连接层来做分类。但是近些年来，越来越多的网络，比如SSD，FasterRCNN的RPN，MTCNN中的PNet，都使用卷积层来代替全连接，也一样可以做到目标分类的效果，而且
- 更灵活，不需要限定输入图像的分辨率；
- 更高效，只需要做一次前向计算。
本文首先对全连接层和卷积层关系做分析，然后比较全连接层和卷积层的优缺点，自然也就搞清楚了为什么用卷积层替代全连接层这个问题了。
全连接和卷积的关系
全连接层和卷积层只要设置好了对应的参数，可以在达到相同输入输出的效果，在这个意义上，在数学上可以认为它们是可以相互替换的。
下面我们以输入10 x 10 x 3的特征图，输出10 x 10 x 1的特征图来证明。
全连接层是一种核很大的卷积层
全连接层怎么做？全连接层输入10 x 10 x 3的特征图，首先将其reshape成一维的300个输入神经元，然后每个输出神经元的值都是这300个输入神经元的线性组合，最后将100个输出reshape成10 x 10 x 1的形状。
卷积层怎么做？卷积层可以直接用100个10 x 10 x 3的滤波器，分别直接贴着输入的10 x 10 x 3的输入特征图做滤波，得到100个一维的输出，然后把这个100个一维输出reshape成10 x 10 x 1的形状。这100个输出，每个也都是由输入特征图上300个像素点线性组合而来，在数学上和上面的全连接层理论上可以达到一样的效果。这里用到的卷积和通常我们看到的卷积唯一有点特殊的是，为了达到全连接（采集到所有输入特征图像素的信息）的效果，它的分辨率和输入特征图分辨率一样。
所以，由上面的讨论可知，全连接成可以用（分辨率和输入特征图相同的）卷积代替。
卷积层是一种稀疏的全连接层
卷积层怎么做？假设卷积核大小为3 x 3的方形核，stride = 1， 为了输出分辨率不变padding用SAME方式。那么卷积核的形状是3 x 3 x 3，输出通道数为1，只用一个这样的卷积核按照常规卷积来做就行，输出10 x 10 x 1的特征图很容易。
全连接层怎么做？全连接层输入10 x 10 x 3 = 300个神经元，卷积的时候每次卷积只是连接了其中的一个3 x 3 x 3 = 27的子集，那么可以在做全连接的时候，除了这27个神经元设置连接关系，其余的 300 - 27 = 273个连接系数直接设置为0就可以了。做 10 x 10 x 1 = 100次这样的全连接，就可以得到100个输出神经元，再reshape成10 x 10 x 1的形状就可以了。
所以，由上面的讨论可以得到，卷积层只是全连接层的一个子集，把全连接的某些连接系数设置为0，就可以达到和卷积相同的效果。
卷积替代全连接的优点
由上一节的讨论可以知道，其实卷积层和全连接层在理论上是可以相互替代的，那么，为什么我们看到的趋势是全连接被卷积层替代，而不是卷积层被全连接层替代呢？这要从卷积层相比于全连接层两个工程上的优点来讲：
对输入分辨率的限制
如果网络后面有全连接层，而全连接层的输入神经元个数就是固定的，那么反推上层卷积层的输出是固定的，继续反推可知输入网络的图片的分辨率是固定的。例如，LetNet由于由全连接层，输入就只能是28 x 28的。如果网络中的全连接层都用卷积层替代，网络中只有卷积层，那么网络的输出分辨率是随着输入图片的分辨率而来的，输出图中每一个像素点都对应着输入图片的一个区域（可以用stride,pooling来反算）。
计算效率比较
同样以LeNet来做例子，如果一个图片是280 x 280的分辨率，为了识别图片中所有的数字（为了简单，假设每个数字都是在这个大图划分为10 x 10的网格中），那么为了识别这100个位置数字，那么至少需要做100次前向；而全卷积网络的特点就在于输入和输出都是二维的图像，并且输入和输出具有相对应的空间结构，我们可以将网络的输出看作是一张heat-map，用热度来代表待检测的原图位置出现目标的概率，只做一次前向就可以得到所有位置的分类概率。
总结
首先在理论上论证了卷积层和全连接层的可互换性质，然后详细分析了在实践中用卷积层代替全连接层的两个好处，第一个是去掉全连接层对网络输入图像分辨率的限制；第二个好处是全卷积网络只需要做一次前向运算就可以获得一张目标所在位置的heat-map，节约了计算。