diff --git a/docs/12_convolution/01_naive_conv/README.md b/docs/12_convolution/01_naive_conv/README.md index 0a02796..a935c2c 100644 --- a/docs/12_convolution/01_naive_conv/README.md +++ b/docs/12_convolution/01_naive_conv/README.md @@ -1 +1,182 @@ -# 卷积算子简易实现 \ No newline at end of file +# 卷积算子简易实现 + +## 1. 什么是卷积运算 + +卷积运算是一种特殊的线性运算,它是一种对两个实变函数的一种数学运算,表征函数之间的一种重要的运算方式。卷积运算的定义如下: + +$$ +(f*g)(t) = \int_{-\infty}^{\infty} f(\tau)g(t-\tau)d\tau +$$ + +其中, $f(t)$ 和 $g(t)$ 是两个实变函数,$*$表示卷积运算。卷积运算的定义是在实数域上的,但是在深度学习中,我们通常使用离散卷积运算。离散卷积运算的定义如下: + +$$ +(f*g)(t) = \sum_{\tau=-\infty}^{\infty} f(\tau)g(t-\tau) +$$ + +其中,$f(t)$ 和 $g(t)$ 是两个离散函数,$*$表示离散卷积运算。 + +在图象处理中,卷积运算是一种对图像进行滤波的操作,它可以提取图像的特征。卷积运算的过程如下: + +1. 将卷积核与图像进行逐元素相乘; +2. 将相乘的结果求和。 + +其中,卷积核是一个小的矩阵,它可以提取图像的特征。卷积核的大小通常是$3\times 3$或$5\times 5$。 +下图是一个卷积运算的示意图: + +![picture 0](images/9ee5c9b3e312cca4ecd53ef2759dae7d3cba5a083f65c353a6c85935bd7ff256.png) + +:::note + +卷积的运算过程就是相应区域的数字相乘后相加 + +::: + +有图像处理相关经验的同学应该知道,一个图像是有三个通道的。那么理论上卷积运算也应该是对三个通道都要进行卷积运算才行。所以我们就会有三个卷积核分别进行计算。 + +![picture 1](images/80eae29b069a39fd446c2551ab8f8cea396190707e2104d6c868a35eb3735154.png) + +卷积是一个特征提取的过程,不同的卷积核可以提取图像中不同的特征。 + +本文将介绍如何使用 CUDA 实现一个卷积算子。 + +## 2. 卷积算子的实现 + +CUDA 代码是从并发的角度来实现卷积算子的。具体来说,我们将使用一个线程来计算一个输出元素。这样,我们可以使用一个线程块来计算一个输出通道的所有输出元素。 + +首先我们先明确一下我们的任务,我们的任务是给定一个矩阵和一个卷积核,计算卷积运算的结果。矩阵的纬度是 $n \times c \times h \times w$,卷积核的纬度是 $k \times c \times r \times s$。其中, $n$ 是 batch size, $c$ 是通道数, $h$ 是数据高, $w$ 是数据宽, $k$ 是卷积核数量, $r$ 是卷积核高, $s$ 是卷积核宽。一些常用的符号定义如下: + +```cpp +int n; // batch szie +int c; // 通道数 +int h; // 数据高 +int w; // 数据宽 +int k; // 卷积核数量 +int r; // 卷积核高 +int s; // 卷积核宽 +int u; // 卷积在高方向上的步长 +int v; // 卷积在宽方向上的步长 +int p; // 卷积在高方向上的补边 +int q; // 卷积在宽方向上的补边 + +int out_h; // 卷积在高方向上的输出大小 +int out_w; // 卷积在宽方向上的输出大小 +``` + +卷积运算的结果的纬度是 $n \times k \times out_h \times out_w$。 $out_h$ 和 $out_w$ 是卷积运算的结果的高和宽它们的计算代码如下: + +```cpp +out_h = (h + 2 * p - r) / u + 1; +out_w = (w + 2 * q - s) / v + 1; +``` + +首先我们先想一下如何组织线程来进行计算,我们可以用一个线程来计算一个输出元素。这样,我们可以使用一个线程块来计算一个输出通道的所有输出元素。输入的矩阵是一个4维矩阵,我们可以用二维的 grid 和二维的 block: + +```cpp +// 定义线程块的大小 +const int blockDim_x = 16; +const int blockDim_y = 16; + +// 计算网格中线程块的数量 +const int gridDim_x = (out_h * out_w + blockDim_x - 1) / blockDim_x; +const int gridDim_y = (k + blockDim_y - 1) / blockDim_y; +``` + +gridDim_x 和 gridDim_y,分别表示在 x 和 y 方向上需要的线程块数量。 + +对于输出特征图的大小,我们需要考虑每个输出像素需要进行卷积操作,因此使用 out_h * out_w 计算总的输出像素数量。将总的输出像素数量除以每个线程块中的线程数量,即 blockDim_x 和 blockDim_y,可以得到在 x 和 y 方向上需要的线程块数量。由于 GPU 线程块的数量必须是整数,所以我们使用 `(total_pixels + blockDim - 1) / blockDim` 来确保线程块数量足够覆盖所有像素,并且没有多余的线程块。 + +接下来我们来看一下卷积运算的实现代码, 首先我们先需要获取线程的索引: + +```cpp +// 根据线程在二维网格中的位置计算出线程在一维线程空间中的全局唯一索引 +int x = blockIdx.x * blockDim.x + threadIdx.x; +int y = blockIdx.y * blockDim.y + threadIdx.y; +// 表示当前处理的 batch 中的索引 +int z = blockIdx.z; +``` + +检查当前线程是否超出了工作范围。如果超出了输出特征图的范围或者卷积核的数量,或者超出了 batch 的数量: + +```cpp +if (x >= out_h * out_w || y >= k || z >= n) +{ + return; +} +``` + +下面我们需要明确几个坐标,首先是这个线程计算的是输出特征图的哪一个像素,然后是这个像素对应的输入特征图的哪一个像素, + +```cpp +// 计算当前线程处理的输出特征图中的索引 +int pos_out_h = x / out_w; +int pos_out_w = x % out_w; + +// 计算输入数据的坐标 +int pos_ori_h = pos_out_h * u - p; +int pos_ori_w = pos_out_w * v - q; + +// 计算输入数据的偏移量,即输入数据中当前像素的起始位置 +int in_offset = z * c * h * w + pos_ori_h * w + pos_ori_w; +// 计算卷积核的偏移量,即卷积核中当前通道的起始位置 +int weight_offset = y * c * r * s; +// 计算输入数据中通道之间的偏移量 (输入数据的大小是 n * c * h * w) +int in_channel_offset = h * w; +// 计算卷积核中通道之间的偏移量 (卷积核的大小是 k * c * r * s) +int weight_channel_offset = r * s; +``` + +在理清楚了坐标之后,我们就可以开始计算卷积运算的结果了,计算过程很简单就是把卷积核和输入特征图对应位置的元素相乘然后求和: + +```cpp +float sum = 0.0f; + +for (int i = 0; i < r; i++) +{ + for (int j = 0; j < s; j++) + { + int pos_real_h = pos_ori_h + i; + int pos_real_w = pos_ori_w + j; + + // 只处理有效的数据点 + if (pos_real_h >= 0 && pos_real_w >= 0 && pos_real_w < w && pos_real_h < h) + { + int in_offset_tmp = in_offset; + int wei_offset_tmp = weight_offset; + for (int channel = 0; channel < c; channel++) + { + // 计算卷积和 + sum += in[in_offset_tmp + i * w + j] * weight[wei_offset_tmp + i * s + j]; + in_offset_tmp += in_channel_offset; + wei_offset_tmp += weight_channel_offset; + } + } + } +} +``` + +最后我们将计算结果写入到输出特征图中: + +```cpp +// 将计算结果写入到输出特征图中 +int out_offset = z * k * out_h * out_w + y * out_h * out_w + x; +out[out_offset] = sum; +``` + +编译运行命令如下: + +```bash +nvcc -o naive_conv naive_conv.cu +./naive_conv +``` + +## 3. 总结 + +本文介绗了卷积运算的定义和卷积运算的实现。卷积运算是一种特殊的线性运算,它是一种对两个实变函数的一种数学运算,表征函数之间的一种重要的运算方式。卷积运算的定义是在实数域上的,但是在深度学习中,我们通常使用离散卷积运算。离散卷积运算的定义是在离散函数上的。卷积运算是一种对图像进行滤波的操作,它可以提取图像的特征。卷积运算的过程是将卷积核与图像进行逐元素相乘,然后将相乘的结果求和。 + +CUDA 代码是从并发的角度来实现卷积算子的。具体来说,我们将使用一个线程来计算一个输出元素。这样,我们可以使用一个线程块来计算一个输出通道的所有输出元素。 + +## References + +1. https://space.keter.top/docs/deep_learning/%E5%9F%BA%E7%A1%80%E7%9F%A5%E8%AF%86/%E5%8D%B7%E7%A7%AF%E7%A5%9E%E7%BB%8F%E7%BD%91%E7%BB%9C +2. http://cs231n.stanford.edu/ diff --git a/docs/12_convolution/01_naive_conv/images/80eae29b069a39fd446c2551ab8f8cea396190707e2104d6c868a35eb3735154.png b/docs/12_convolution/01_naive_conv/images/80eae29b069a39fd446c2551ab8f8cea396190707e2104d6c868a35eb3735154.png new file mode 100644 index 0000000..559b569 Binary files /dev/null and b/docs/12_convolution/01_naive_conv/images/80eae29b069a39fd446c2551ab8f8cea396190707e2104d6c868a35eb3735154.png differ diff --git a/docs/12_convolution/01_naive_conv/images/9ee5c9b3e312cca4ecd53ef2759dae7d3cba5a083f65c353a6c85935bd7ff256.png b/docs/12_convolution/01_naive_conv/images/9ee5c9b3e312cca4ecd53ef2759dae7d3cba5a083f65c353a6c85935bd7ff256.png new file mode 100644 index 0000000..2901cf6 Binary files /dev/null and b/docs/12_convolution/01_naive_conv/images/9ee5c9b3e312cca4ecd53ef2759dae7d3cba5a083f65c353a6c85935bd7ff256.png differ diff --git a/docs/12_convolution/01_naive_conv/naive_conv.cu b/docs/12_convolution/01_naive_conv/naive_conv.cu new file mode 100644 index 0000000..c2297ac --- /dev/null +++ b/docs/12_convolution/01_naive_conv/naive_conv.cu @@ -0,0 +1,239 @@ +#include +#include +#include + +/* + * @param n: batch size + * @param c: 通道数 + * @param h: 输入数据高 + * @param w: 输入数据宽 + * @param k: 卷积核数量 + * @param r: 卷积核高 + * @param s: 卷积核宽 + * @param out_h: 输出数据高 + * @param out_w: 输出数据宽 + * @param u: 卷积在高方向上的步长 + * @param v: 卷积在宽方向上的步长 + * @param p: 卷积在高方向上的补边 + * @param q: 卷积在宽方向上的补边 + * @param in: 输入数据 + * @param weight: 卷积核 + * @param out: 输出数据 + */ +__global__ void +naive_conv2d_kernel(int n, int c, int h, int w, + int k, int r, int s, + int out_h, int out_w, + int u, int v, int p, int q, + float *in, float *weight, float *out) +{ + // 获取线程在三维网格中的位置 + int x = blockIdx.x * blockDim.x + threadIdx.x; + int y = blockIdx.y * blockDim.y + threadIdx.y; + int z = blockIdx.z; + + // 如果线程超出工作范围则退出 + if (x >= out_h * out_w || y >= k || z >= n) + { + return; + } + + // 当前线程处理的数据点在out_h、out_w上的坐标 + int pos_out_h = x / out_w; + int pos_out_w = x % out_w; + + // 计算输入数据的坐标 + int pos_ori_h = pos_out_h * u - p; + int pos_ori_w = pos_out_w * v - q; + + float sum = 0.0; + + int in_offset = z * c * h * w + pos_ori_h * w + pos_ori_w; + int weight_offset = y * c * r * s; + int in_channel_offset = h * w; + int weight_channel_offset = r * s; + + // 执行卷积操作 + for (int i = 0; i < r; i++) + { + for (int j = 0; j < s; j++) + { + int pos_real_h = pos_ori_h + i; + int pos_real_w = pos_ori_w + j; + + // 只处理有效的数据点 + if (pos_real_h >= 0 && pos_real_w >= 0 && pos_real_w < w && pos_real_h < h) + { + int in_offset_tmp = in_offset; + int wei_offset_tmp = weight_offset; + for (int channel = 0; channel < c; channel++) + { + // 计算卷积和 + sum += in[in_offset_tmp + i * w + j] * weight[wei_offset_tmp + i * s + j]; + in_offset_tmp += in_channel_offset; + wei_offset_tmp += weight_channel_offset; + } + } + } + } + + // 计算输出偏移 + int out_offset = z * k * out_h * out_w + y * out_h * out_w + x; + out[out_offset] = sum; +} + +// CPU 端的卷积计算 +void conv2d_cpu(float *in, float *pwei, float *out, int n, int c, int h, int w, int k, int r, int s, int u, int v, int p, int q) +{ + int out_h = (h + 2 * p - r) / u + 1; + int out_w = (w + 2 * q - s) / v + 1; + + for (int n_num = 0; n_num < n; n_num++) + { + for (int k_num = 0; k_num < k; k_num++) + { + for (int i = 0; i < out_h; i++) + { + for (int j = 0; j < out_w; j++) + { + double sum = 0.0; + int pos_h = i * u - p; + int pos_w = j * v - q; + + for (int c_num = 0; c_num < c; c_num++) + { + for (int kh_num = 0; kh_num < r; kh_num++) + { + for (int kwNum = 0; kwNum < s; kwNum++) + { + int pos_ori_h = pos_h + kh_num; + int pos_ori_w = pos_w + kwNum; + if (pos_ori_w >= 0 && pos_ori_h >= 0 && pos_ori_w < w && pos_ori_h < h) + { + sum += (double)(in[n_num * c * h * w + c_num * (w * h) + pos_ori_h * w + pos_ori_w] * pwei[k_num * r * s * c + c_num * r * s + kh_num * s + kwNum]); + } + } + } + } + + out[n_num * k * out_h * out_w + k_num * out_h * out_w + i * out_w + j] = (float)sum; + } + } + } + } +} + +int main() +{ + // 定义输入数据和卷积核的尺寸 + const int n = 2; // batch size + const int c = 2; // 通道数 + const int h = 10; // 数据高 + const int w = 10; // 数据宽 + const int k = 5; // 卷积核数量 + const int r = 3; // 卷积核高 + const int s = 3; // 卷积核宽 + const int out_h = (h - r + 2 * 0) / 1 + 1; // 输出高 + const int out_w = (w - s + 2 * 0) / 1 + 1; // 输出宽 + const int u = 1; // 卷积在高方向上的步长 + const int v = 1; // 卷积在宽方向上的步长 + const int p = 0; // 卷积在高方向上的补边 + const int q = 0; // 卷积在宽方向上的补边 + + // 分配内存并随机生成输入数据和卷积核 + float *in, *weight, *out; + float *in_device, *weight_device, *out_device; + + in = (float *)malloc(n * c * h * w * sizeof(float)); + weight = (float *)malloc(k * c * r * s * sizeof(float)); + out = (float *)malloc(n * k * out_h * out_w * sizeof(float)); + + cudaMalloc((void **)&in_device, n * c * h * w * sizeof(float)); + cudaMalloc((void **)&weight_device, k * c * r * s * sizeof(float)); + cudaMalloc((void **)&out_device, n * k * out_h * out_w * sizeof(float)); + + // 随机生成输入数据和卷积核 + for (int i = 0; i < n * c * h * w; ++i) + { + in[i] = (float)rand() / RAND_MAX; + } + for (int i = 0; i < k * c * r * s; ++i) + { + weight[i] = (float)rand() / RAND_MAX; + } + + // 将输入数据和卷积核拷贝到 GPU + cudaMemcpy(in_device, in, n * c * h * w * sizeof(float), cudaMemcpyHostToDevice); + cudaMemcpy(weight_device, weight, k * c * r * s * sizeof(float), cudaMemcpyHostToDevice); + cudaMemcpy(out_device, out, n * k * out_h * out_w * sizeof(float), cudaMemcpyHostToDevice); + + // 定义线程块的大小 + const int blockDim_x = 16; + const int blockDim_y = 16; + + // 计算线程块和网格的数量 + const int gridDim_x = (out_h * out_w + blockDim_x - 1) / blockDim_x; + const int gridDim_y = (k + blockDim_y - 1) / blockDim_y; + + // 定义线程块和网 + + dim3 blockDim(blockDim_x, blockDim_y); + dim3 gridDim(gridDim_x, gridDim_y, n); + + // 调用 kernel 函数 + naive_conv2d_kernel<<>>(n, c, h, w, k, r, s, out_h, out_w, u, v, p, q, in_device, weight_device, out_device); + // 同步 + cudaDeviceSynchronize(); + + // 将 GPU 计算的结果拷贝到 CPU + cudaMemcpy(out, out_device, n * k * out_h * out_w * sizeof(float), cudaMemcpyDeviceToHost); + + // CPU 端进行卷积计算 + float *out_cpu = (float *)malloc(n * k * out_h * out_w * sizeof(float)); + conv2d_cpu(in, weight, out_cpu, n, c, h, w, k, r, s, u, v, p, q); + + // 比较 GPU 和 CPU 计算结果是否一致 + bool pass = true; + for (int i = 0; i < n * k * out_h * out_w; ++i) + { + if (abs(out[i] - out_cpu[i]) > 1e-5) + { + pass = false; + std::cout << "Verification failed at " << i << "!" << std::endl; + std::cout << "GPU: " << out_cpu[i] << " CPU: " << out[i] << std::endl; + break; + } + } + + if (pass) + { + std::cout << "Verification passed!" << std::endl; + + int iter = 100; + cudaEvent_t start, stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + cudaEventRecord(start, 0); + for (int i = 0; i < iter; i++) + { + naive_conv2d_kernel<<>>(n, c, h, w, k, r, s, out_h, out_w, u, v, p, q, in_device, weight_device, out_device); + } + cudaEventRecord(stop, 0); + cudaEventSynchronize(stop); + float elapsedTime; + cudaEventElapsedTime(&elapsedTime, start, stop); + std::cout << "GPU time: " << 1000 * elapsedTime / iter << "us" << std::endl; + cudaEventDestroy(start); + cudaEventDestroy(stop); + } + + // 释放内存 + cudaFree(in_device); + cudaFree(weight_device); + cudaFree(out_device); + free(in); + free(weight); + free(out); + + return 0; +} \ No newline at end of file