conv2d之forward

1 直接卷积

在直接卷积中，卷积核（滤波器）通过在输入图像上滑动并与图像的每个位置进行逐元素相乘，然后将所有乘积相加得到输出的单个像素值。这个过程可以看作是在输入图像上进行的一种滑动窗口操作，其中卷积核的大小决定了窗口的大小。下图所示为一个输入进行卷积的过程，如果为批量输入，则结果也多个。维度变化为：

其中：

/*
不考虑输出的padding
input的维度：[batch, h, w, channel]
kernel的维度：[kerner_cnt, kh, kw, channel]
*/
int output_h = (h - kh) / sh + 1;
int output_w = (w - kw) / sw + 1;

for (int b = 0; b < batch; b ++)
    for (int i = 0; i < output_h; i ++)
        for (int j = 0; j < output_w; j ++)
            for (int d = 0; d < kernel_cnt; d ++) 
                for (int c = 0; c < channel; c ++)
                    for (int ph = 0; ph < kh; ph ++)
                        for (int pw = 0; pw < kw; pw ++)
                            res[b][i][j][d] += kernel[d][ph][pw][c] * input[b][i * sh + ph][j * sw +  pw][c];

2 Img2col+Gemm算法实现卷积

2.1 概述

最广泛的卷积算法之一就是基于 Img2col+GEMM 算法来实现卷积。它包括以下步骤：

将输入转换为一个大矩阵，其中每一行对应一个输出位置上的卷积 patch（展平后的输入块）。

另一个矩阵由所有卷积核展平后按列拼接得到。最后通过 GEMM 完成卷积，如图所示。

转换后的输入矩阵按行存放 patch：GEMM 中 的每一行与 的每一列做内积，等价于一个 patch 与一个卷积核做点乘。 的每一行对应一个输出位置， 的每一列对应一个卷积核，因此 的每个元素就是卷积结果。用于获取转换后矩阵的函数通常称为 im2col。

由于GPU体系结构的性质和线性代数库中GEMM方法的高效实现，卷积非常适合于放在GPU上计算。然而，这种方法需要大量的内存来存储转换后的矩阵，特别是转换后的输入矩阵，由于卷积中有滑动步长和Padding的存在，转换后的数据矩阵必须存储滑动中重复访问的元素，和一些Padding后额外数据。这在神经网络训练过程中是十分占用显存的，于是就有另一种GEMM实现。

2.2 计算过程

2.2.1 几何角度

• 把 个 的卷积核各自展平为列向量，按列拼成 ，记为 Weight。
• 把 的输入按卷积窗口展平：每个输出位置对应一行 patch，所有 patch 按行堆叠得到 ，记为 Input。设卷积步长为 ，则：

• 以上两步都是 Im2col；之后对每个 batch 执行 GEMM：

• 的每一列对应一个卷积核在整个特征图上的响应；遍历 batch 得到 ，再 reshape 为 。

优点：

• 通过向量化把原本不连续的内存访问变成连续的内存访问
• 将计算转换为高效的矩阵乘法，达到运算加速的目的。

缺点：

• 存储展平后的输入和卷积核需要开辟新的空间。
• 通常卷积步长会小于卷积核的大小，参与卷积运算的相邻块有重叠，向量化的代价是冗余存储。

2.2.2 公式角度

符号约定（与后文代码 [B,H,W,C]、[D,KH,KW,C] 一致）：

• 输入 ，卷积核 ，输出
• patch 展平顺序为 ，与 im2col_lhs_new 中 for h; for w; for c 一致

统一卷积公式（VALID、无 padding）：

以下分情形说明上式如何对应到 Im2col + GEMM 的 。

2.2.2.1 单输入单核单通道单步长

卷积计算的输入从二维的矩阵到四维张量，以及卷积核从二维矩阵到四维矩阵都对应不同大小的输出，最后会统一到一个算法里面，这里用一个函数来表示：

其中 是输入， 是卷积核， 是卷积的结果，即特征图（feature map），先看输入和卷积核都是矩阵的情况：

上图中蓝色的为输入，红色的是卷积核，灰色是结果

卷积的过程就是，将输入划分成若干个与卷积核相同大小的不同子集，再分别与卷积核点乘：

如上图所示，本栗子的输入可以划分为4个子子集（为什么是4个，这里涉及一个步长的变量，后面讨论），将它们各自与卷积核点乘，得到的结果就是红框上面的输出

现在，将上述的过程用公式向量表示，就拿第一个子集（将其表示为 ）来说吧：

然后将其展开，变成一个行向量：

对四个 patch 分别展平为行向量后，按行堆叠成 Im2col 矩阵（）：

卷积核同样展平为列向量 ：

然后将 与 相乘：

最后将变形，就得到了卷积的结果：

以上，就是输入和卷积核都是矩阵时的情况。当 时，可写成（单核单通道，记 ）：

其中 是卷积核高和宽，下标从 0 开始； 是输出位置。

实现代码

int output_h = (h - kh) / sh + 1;
int output_w = (w - kw) / sw + 1;
vec2d res(output_h, vec1d(output_w));
for (int i = 0; i < output_h; i ++)
    for (int j = 0; j < output_w; j ++)
        for (int ph = 0; ph < kh; ph ++)
            for (int pw = 0; pw < kw; pw ++)
                res[i][j] += kernel[ph][pw] * input[i + ph][j + pw];

2.2.2.2 单输入单核多通道单步长

现在将输入升级为多通道（例如图片一般都是三通道的），即变成一个三维张量，相对应的，之前的情况称之为单通道；同时，卷积核也必须升级成多通道的，为什么？因为卷积核的通道数必须与输入的通道数相等：

多通道

如上图所示，输入变成双通道后，卷积核也必须变成双通道的，但是，输出还是一个通道，下面看一下计算过程，先画一个详细的图示：

将两个通道的 patch 按 顺序展平（与代码 for h; for w; for c 一致），得到长度 18 的行向量，记第一个 patch 为 。图示中上下拼接的两个 子块分别对应 两个通道：

对四个 patch 重复上述操作，得到 ；卷积核展平为 ，展平顺序与输入相同。

然后，计算卷积结果（，）：

可以看到，多通道的情况与单通道的情况基本上是没有区别的；对应公式为：

代码实现

for (int i = 0; i < output_h; i ++)
    for (int j = 0; j < output_w; j ++)
        for (int c = 0; c < channel; c ++) 
            for (int ph = 0; ph < kh; ph ++)
                for (int pw = 0; pw < kw; pw ++)
                    res[i][j] += kernel[ph][pw][c] * input[i + ph][j + pw][c];

2.2.2.3 单输入多核多通道单步长

接下来，在多通道的基础上，进一步升级，使用多个卷积核，大致如下图：

多卷积核

可以看到，现在输出的结果变成多通道了，也就是说，卷积结果的通道数等于卷积核的个数。计算过程大致如下：

因为输入没有变化，所以划分出来的 Im2col 矩阵与上一节相同：

两个卷积核分别按 顺序展平为列向量。第一个卷积核（）记为 ：

第二个卷积核（）同样展平为 ，再按列拼成权重矩阵：

计算卷积结果：

从上式可知，每个卷积核独立与输入做卷积。将 reshape 为两个输出通道的特征图 （，）：

其中第一个下标 表示卷积核编号，后两个下标为空间位置。

再上一个从网上找来的动图，相信能更好的理解

图片来自网络

以上，就是多卷积核多通道的 Im2col+GEMM 过程，对应公式为：

其中 是卷积核序号。

代码实现

for (int d = 0; d < kernel_cnt; d ++) 
    for (int i = 0; i < output_h; i ++)
        for (int j = 0; j < output_w; j ++)
            for (int c = 0; c < channel; c ++)
                for (int ph = 0; ph < kh; ph ++)
                    for (int pw = 0; pw < kw; pw ++)
                        res[i][j][d] += kernel[d][ph][pw][c] * input[i + ph][j + pw][c];

2.2.2.4 多输入多核多通道单步长

有了前面多卷积核的计算过程做对比后，多输入的计算过程也是大同小异的，在前面的多卷积核的计算中，每个卷积核都独立与输入进行卷积，现在有多个输入也是一样，每一个输入独立与所有卷积核进行卷积计算得到一个输出，所以输入的数量与输出的数量一致

批量输入

还是来看栗子吧，设输入的 batch 为 2，对第一个输入采用与前面一样的操作，变换为一个矩阵：

第二个输入变换为：

然后将两个 batch 的矩阵按行拼接，在数学上等价于一次更大的 GEMM；文末 conv2d 实现中则是对每个 batch 独立调用 GEMM，结果相同：

如果非常熟悉矩阵乘法，肯定已经知道了输出的结果：

上面 至 是 与 相乘的结果， 至 是 与 相乘的结果。更一般地：

所以，最后的公式就是（）：

其中 是 batch 序号。

代码实现

for (int b = 0; b < batch; b ++)
    for (int i = 0; i < output_h; i ++)
        for (int j = 0; j < output_w; j ++)
            for (int d = 0; d < kernel_cnt; d ++) 
                for (int c = 0; c < channel; c ++)
                    for (int ph = 0; ph < kh; ph ++)
                        for (int pw = 0; pw < kw; pw ++)
                            res[b][i][j][d] += kernel[d][ph][pw][c] * input[b][i + ph][j +  pw][c];

2.2.2.5 多输入多核多通道多步长

最后，讨论卷积步长。前面各节默认 ；当步长增大时，输出位置 对应的输入窗口起点变为 ：

上图是步长为1时的卷积结果，如果令步长为2，那么就变成了：

上图中的卷积结果为灰色的表示被表过，可以看到，步长为2时的情况相当于在步长为1的计算结果基础上每个一步取一个结果，同样，步长为3甚至更大也是依次类推

不过有时候输入的宽高不能恰好被卷积核大小和步长整除，这时需要根据实现策略做 padding 或截断。引入步长后的统一公式为：

Im2col 中 patch 的采样位置与上式一致：输出 对应输入窗口左上角 ，再按 顺序展平后与 做 GEMM。

代码实现

for (int b = 0; b < batch; b ++)
    for (int i = 0; i < output_h; i ++)
        for (int j = 0; j < output_w; j ++)
            for (int d = 0; d < kernel_cnt; d ++) 
                for (int c = 0; c < channel; c ++)
                    for (int ph = 0; ph < kh; ph ++)
                        for (int pw = 0; pw < kw; pw ++)
                            res[b][i][j][d] += kernel[d][ph][pw][c] * input[b][i * sh + ph][j * sw +  pw][c];

2.2.3 Img2col+Gemm的代码实现

// [b, h, w, c] -> [b, out_h * out_w, kh * kw * C]
vec3d im2col_lhs_new(vec4d v, int kh, int kw, int sh, int sw) {

    int B = v.size(), H = v[0].size(), W = v[0][0].size(), C = v[0][0][0].size();

    int out_h = (H - kh) / sh + 1;
    int out_w = (W - kw) / sw + 1;

    vec3d output(B, vec2d(out_h * out_w, vec1d(kh * kw * C)));

    for (int b = 0; b < B; b ++) {
        int row = 0;
        for (int h = 0; h < out_h; h ++)
            for (int w = 0; w < out_w; w ++) {
                int idx = 0;
                for (int i = 0; i < kh; i ++)
                    for (int j = 0; j < kw; j ++)
                        for (int c = 0; c < C; c ++)
                            output[b][row][idx ++] = v[b][sh * h + i][sw * w + j][c];
                row ++;
            }
    }
    return output;
}


// [D, kh, kw, C] -> [kh * kw * C, D]
vec2d im2col_rhs_new(vec4d v) {
    int Dout = v.size();
    int Hout = v[0].size();
    int Wout = v[0][0].size();
    int Cout = v[0][0][0].size();
        
    vec2d res(Hout * Wout * Cout, vec1d(Dout));

    for (int d = 0; d < Dout; d ++)
        for (int i = 0; i < Hout; i ++) 
            for (int j = 0; j < Wout; j ++)
                for (int k = 0; k < Cout; k ++)
                        res[i * Wout * Cout + j * Cout + k][d] = v[d][i][j][k];        
        
    return res;
}

// [out_h * out_w, kh * kw * C] dot [kh * kw * C, D]->[out_h * out_w, D]
vec2d gemm(vec2d a, vec2d b) {
    int n = a.size(), m = a[0].size(), l = b[0].size();
    vec2d c(n, vec1d(l));
    for (int i = 0; i < n; i ++)
        for (int k = 0; k < m; k ++)
            for (int j = 0; j < l; j ++)
                c[i][j] += a[i][k] * b[k][j];
    return c;
}

// [B, out_h, out_w, d]
vec4d conv2d(vec4d input, vec4d weight, int sh, int sw) {

    int H = input[0].size();
    int W = input[0][0].size();
    int D = weight.size();
    int kh = weight[0].size();
    int kw = weight[0][0].size();
    int out_h = (H - kh) / sh + 1;
    int out_w = (W - kw) / sw + 1;

    vec3d input_3d = im2col_lhs_new(input, kh, kw, sh, sw);
    vec2d weight_2d = im2col_rhs_new(weight);

    int B = input.size();

    vec4d ans = vec4d(B, vec3d(out_h, vec2d(out_w, vec1d(D))));
    cout << "维度： [" << B << ' ' << out_h << ' ' << out_w << ' ' << D << "]\n"; 
    for (int b = 0; b < B; b ++) {
        vec2d lhs = input_3d[b]; // [out_h * out_w, kh * kw * C]
        vec2d rhs = weight_2d; // [kh * kw * C, D]

        vec2d res = gemm(lhs, rhs); // [out_h * out_w, D]

        for (int d = 0; d < D; d ++)
            for (int i = 0; i < out_h; i ++)
                for (int j = 0; j < out_w; j ++)
                    ans[b][i][j][d] += res[i * out_w + j][d];
    }

    return ans;
}

3 Implicit GEMM算法实现卷积

Implicit GEMM是一种用于实现卷积操作的方法，与Img2col+GEMM相似，同样利用了矩阵乘法的性质来加速卷积计算。不同点在于，在内存使用量方面，Implicit GEMM不需要任何额外的存储空间。

在Img2col+GEMM实现中，实现的步骤为：

• 输入输出矩阵转换，用另一块内存空间保存转换后的输入输出矩阵。
• 将输入输出矩阵进行GEMM运算，然后进行输出。

而在 Implicit GEMM中，矩阵转换发生在计算中，也就不需要内存来保存转换后的矩阵，而是在GEMM中进行输入输出的坐标映射，从而实现卷积运算。

Direct Convolution

Im2col转换后的GEMM实现卷积

在 Implicit GEMM中，主要思想是坐标的映射，怎样使GEMM在load时正确的取到对应坐标元素是关键。参考以下矩阵乘：

for (int i = 0; i < M; i++) {
    for (int j = 0; j < N; j++) {
        result[i][j] =0;
        for (int k = 0; k < K; k++) {
            result[i][j] += matrix1[i][k] * matrix2[k][j];
        }
    }
}

如果把上图中的Filter作为matrix1，Input作为matrix2，那么GEMM可以变成这样来实现卷积

for (int i = 0; i < K; i++) {
    for (int j = 0; j < N*Oh*Ow; j++) {
        output[i][j] =0;
        for (int k = 0; k < C*R*S; k++) {
            result[i][j] += filter[i][k] * input[k][j];
        }
    }
}

这个矩阵乘已经可以实现卷积运算，但是注意到Input、Filter和Output并不是二维数据（如果没有经过Im2col转换），这时就要进行坐标变换。

现在假设Input的数据排列格式为NCHW，Filter为KCRS，则输出为NKOhOw，卷积步长为STRIDE。

• Output坐标映射。

按照GEMM来说，Output由i和j两个循环变量决定，则应由这两个变量计算出NKOhOw四个方向的坐标。

n = j  /（Oh*Ow）                      //N维度坐标
k = i                                  //K维度坐标
oh =( j %（Oh*Ow）) / Ow               //Oh维度坐标
ow =( j %（Oh*Ow）) % Ow               //Ow维度坐标

• Filter坐标映射。

由于为KCRS排列，若把CRS当作一个维度，则Filter可以是一个K行，CRS列的矩阵。此时的循环变量i控制在（0，K-1）内，代表当前处理的是第k个卷积核，k控制在（0，CRS-1）内，代表在当前卷积核内积过程。

k = i                                  //K维度坐标
c = k / (R*S)                          //C维度坐标
r = k % (R*S) / S                      //R维度坐标
s = k % (R*S) % S                      //S维度坐标

• Input坐标映射。

Input同样有四个维度，由循环变量j和k，还有Filter坐标中的r和s决定。这里的H和W的与当前计算结果的Oh和Ow相关，还有内积过程中的R和S相关。具体公式可表达为：Input.H/W = Output.Oh/Ow * STRIDE - Padding + R/S。则四个维度方向的表示为

n = j  /（Oh*Ow）                      //N维度坐标
c = k / (R*S)                          //C维度坐标
h = oh * STRIDE + r                    //H维度坐标
w = ow * STRIDE + s                    //W维度坐标

接下来我们重新改写GEMM实现的卷积算法

//Implicit GEMM Convolution

for (int i = 0; i < K; i++) {
    for (int j = 0; j < N*Oh*Ow; j++) {
        int on = j/(Oh*Ow);                    //N维度坐标
        int oh = (j%(Oh*Ow))/Ow;               //Oh维度坐标
        int ow = (j%(Oh*Ow))%Ow;               //Ow维度坐标
        output[on][i][oh][ow] =0;
        for (int k = 0; k < C*R*S; k++) {
            int ic = k/(R*S);                  //C维度坐标
            int ir = k%(R*S)/S;                //R维度坐标
            int is = k%(R*S)%S;                //S维度坐标
            int ih = oh*STRIDE + ir;           //H维度坐标
            int iw = ow*STRIDE + is;           //W维度坐标
            output[on][i][oh][ow] += filter[i][ic][ir][is] * input[on][ic][ih][iw];
        }
    }
}