高速卷积：im2col

参考 from 机器之心：如何实现高速卷积

Introduction

对于一个卷积运算，例如我们所熟悉的高斯模糊，通过高斯滤波器和目标图像进行卷积操作来得到我们所需要的一个模糊的结果。很容易想到就是朴素卷积，如果简单考虑单通道的图像的话：

for ow in 0..output_width:
	for oh in 0..output_height:
		for kr in 0..kernel_row:
			for kc in 0..kernel_col:
				output_img[ow, oh] += kernel[kr, kc] * input_img[ow + kr, oh + kc]

显然这是一个4个嵌套的庞大的循环，花费的时间可想而知，是非常可怕的。这其中涉及到的运算$2fc\times fr\times ow\times oh-1$ 次乘法和加法的运算。

不过，在卷积神经网络，以及在此基础上实现的各种神经网络来说，这种程度的开销显然是不能接受的，而事实上，对于一个类似的卷积，在神经网络的一层中其实只需要毫秒级别的开销（例如在caffe框架上面）。因此，山人自有妙计，要用魔法打败魔法。

当我们谈到卷积的时候，往往会联想到的另一个运算就是矩阵相乘，即 matmul，或者是叫做通用矩阵相乘（Generalized Matrix Multiplication, GEMM），对于熟悉深度学习的人来说必不陌生。相比于卷积来说，矩阵相乘非常快速，因为涉及到的仅仅只是对应位置的点乘，对两个相乘的矩阵来说都只需要遍历一次即可。而卷积从矩阵相乘的角度来看，事实上就是随着模板移动的循环矩阵相乘的过程。外部循环不可避免，那只能从内部循环找到突破口。首先肯定能想到的一点就是：虽然每次相乘的图像局部不同，但是模板是不变的；另外一点就是：模板和局部的点乘相加，其实相当于矩阵乘法中的行列相乘。

所以，如果可以将每一次局部相乘累加的过程转换成两个矩阵中的一行和一列的相乘的话，那么，我们就可以把整个卷积过程转换成一个矩阵相乘的过程，大幅度提高卷积速度。而另一方面，对于卷积核来说，每次运算都是不变的，因此，最终可以将卷积过程转换为：
$$\left[\begin{array}{ccccccc}{p}_{0,0}& p_{0,1}& \dots & p_{0,kc}& p_{1,0}& \dots & p_{kr,kc}\\ p_{0,s}& p_{0,s+1}& \dots & p_{0,kc+s}& p_{1,s}& \dots & p_{kr+s,kc+s}\\ ⋮& ⋮& \ddots & & & & ⋮\\ ⋮& ⋮& & \ddots & & & ⋮\\ ⋮& ⋮& & & \ddots & & ⋮\\ ⋮& ⋮& & & & \ddots & ⋮\\ p_{oh-s+1,ow-s+1}& p_{oh-s+1,ow-s+2}& \dots & p_{oh-s,ow}& p_{oh-s+1,ow-s+1}& \dots & p_{oh,ow}\end{array}\right]\left[\begin{array}{c}{k}_{0,0}\\ k_{0,1}\\ ⋮\\ k_{0,kc}\\ k_{1,0}\\ ⋮\\ k_{kr,kc}\end{array}\right]$$
其中：

• $s$：步长
• $ow,oh$：输出下标
• $kr,kc$：模板下标

例如一个 $3\times 3$ 的卷积核，步长为 $1$，$5\times 5$ 的图像的 im2col 的结果为：
$$\left[\begin{array}{ccccccccc}{p}_{0,0}& p_{0,1}& p_{0,2}& p_{1,0}& p_{1,1}& p_{1,2}& p_{2,0}& p_{2,1}& p_{2,2}\\ p_{0,1}& p_{0,2}& p_{0,3}& p_{1,1}& p_{1,2}& p_{1,3}& p_{2,1}& p_{2,2}& p_{2,3}\\ p_{0,2}& p_{0,3}& p_{0,4}& p_{1,2}& p_{1,3}& p_{1,4}& p_{2,2}& p_{2,3}& p_{2,4}\\ p_{1,0}& p_{1,1}& p_{1,2}& p_{2,0}& p_{2,1}& p_{2,2}& p_{3,0}& p_{3,1}& p_{3,2}\\ p_{1,1}& p_{1,2}& p_{1,3}& p_{2,1}& p_{2,2}& p_{2,3}& p_{3,1}& p_{3,2}& p_{3,3}\\ p_{1,2}& p_{1,3}& p_{1,4}& p_{2,2}& p_{2,3}& p_{2,4}& p_{3,2}& p_{3,3}& p_{3,4}\\ p_{2,0}& p_{2,1}& p_{2,2}& p_{3,0}& p_{3,1}& p_{3,2}& p_{4,0}& p_{4,1}& p_{4,2}\\ p_{2,1}& p_{2,2}& p_{2,3}& p_{3,1}& p_{3,2}& p_{3,3}& p_{4,1}& p_{4,2}& p_{4,3}\\ p_{2,2}& p_{2,3}& p_{2,4}& p_{3,2}& p_{3,3}& p_{3,4}& p_{4,2}& p_{4,3}& p_{4,4}\end{array}\right]\left[\begin{array}{c}{k}_{0,0}\\ k_{0,1}\\ k_{0,2}\\ k_{1,0}\\ k_{1,1}\\ k_{1,2}\\ k_{2,0}\\ k_{2,1}\\ k_{2,2}\end{array}\right]$$
然后，我们就可以通过矩阵相乘的方法（例如 eigen 等一些线性代数库里的实现）快速实现卷积

C++ Implementation

这个是 caffe 框架下的实现，具体过程如上面所示，本人对核心函数的一些涉及到卷积的相关参数作了适当的注释，可以参考一下：

#include<iostream>

using namespace std;

bool is_a_ge_zero_and_a_lt_b(int a,int b)
{
  if(a>=0 && a <b) return true;
    return false;
}
/*
data:
    data_im - Data of the input image
    channels - Number of channels of the input image
    data_col - im2col's result
image and kernel's properties:
    height - Image's height
    width - Image's width
    kernel_h - Height of kernel
    kernel_w - Width of kernel
parameters:
    pad_h - padding along height
    pad_w - padding along width
    stride_h - stride along height
    stride_w - Stride along width
    dilation_h - Dilation along height
    dilation_w - Dilation along width
*/

void im2col_cpu(const float* data_im, const int channels,
    // basic properties of input image and kernel
    const int height, const int width, const int kernel_h, const int kernel_w,
    // which means length of kernel's padding from boundry of the input image
    const int pad_h, const int pad_w,
    // which means stride of kernel's moving
    const int stride_h, const int stride_w,
    // which means stride between pixels inside the very region in the input image
    const int dilation_h, const int dilation_w,
	// 输出结果
    float* data_col) {


  // calculate output height and output width according to parameters
  // details
  //    First, we calculate height with padding
  //    And then calculate kernel's exact height during calculation
  //    Finally we determine the output height by dividing the difference of the two results above with stride along height
  //    Alternatively calculating width of output
  // 具体实现视dilation的定义范围而不同
  // 如果 dilation 从1开始，即1表示没有膨胀
  const int output_h = (height + 2 * pad_h - (dilation_h * (kernel_h - 1/*here is 1 when dilation not begin with 0, see following explanation*/) + 1)) / stride_h + 1;
  const int output_w = (width + 2 * pad_w - (dilation_w * (kernel_w - 1) + 1)) / stride_w + 1;
  // 如果 dilation 从0开始，即0表示没有膨胀
  // const int output_h = (height + 2 * pad_h - (dilation_h * (kernel_h - 1) + kernel_h)) / stride_h + 1;
  // const int output_w = (width + 2 * pad_w - (dilation_w * (kernel_w - 1) + kernel_w)) / stride_w + 1;

  const int channel_size = height * width;
  // channels
  for (int channel = channels; channel--; data_im += channel_size) {
    // kernel's row
    for (int kernel_row = 0; kernel_row < kernel_h; kernel_row++) {
      // kernel's column
      for (int kernel_col = 0; kernel_col < kernel_w; kernel_col++) {
        int input_row = -pad_h + kernel_row * dilation_h;
        // convolution's reuslt's row
        for (int output_rows = output_h; output_rows; output_rows--) {
          // 这里利用矩阵存储的连续性，通过指针的移动，减少变量的使用
          if (!is_a_ge_zero_and_a_lt_b(input_row, height)) {
            // convolution's result's column
            for (int output_cols = output_w; output_cols; output_cols--) {
              *(data_col++) = 0;
            }

          } else {

            int input_col = -pad_w + kernel_col * dilation_w;

            for (int output_col = output_w; output_col; output_col--) {

              if (is_a_ge_zero_and_a_lt_b(input_col, width)) {

                *(data_col++) = data_im[input_row * width + input_col];

              } else {

                *(data_col++) = 0;

              }

              input_col += stride_w;

            }
          }
          input_row += stride_h;
        }
      }
    }
  }
}