Optimizing Memory Efficiency for Convolution Kernels on Kepler GPUs

TL;DR

This paper presents optimized convolution kernels for Kepler GPUs that significantly improve memory efficiency and performance, addressing the mismatch between GPU memory architecture and computation data width.

Contribution

It introduces a general model for memory-data width mismatch and develops two optimized convolution kernels, achieving substantial performance gains over cuDNN.

Findings

Achieved 5.16X performance improvement for the special case.

Achieved 35.5% average performance improvement for the general case.

Developed memory access and computation pattern optimizations.

Abstract

Convolution is a fundamental operation in many applications, such as computer vision, natural language processing, image processing, etc. Recent successes of convolutional neural networks in various deep learning applications put even higher demand on fast convolution. The high computation throughput and memory bandwidth of graphics processing units (GPUs) make GPUs a natural choice for accelerating convolution operations. However, maximally exploiting the available memory bandwidth of GPUs for convolution is a challenging task. This paper introduces a general model to address the mismatch between the memory bank width of GPUs and computation data width of threads. Based on this model, we develop two convolution kernels, one for the general case and the other for a special case with one input channel. By carefully optimizing memory access patterns and computation patterns, we design a communication-optimized kernel for the special case and a communication-reduced kernel for the general case. Experimental data based on implementations on Kepler GPUs show that our kernels achieve 5.16X and 35.5% average performance improvement over the latest cuDNN library, for the special case and the general case, respectively.

Tables

Table 1. Table 1: Best configurations of our general case convolution kernel for different filter sizes for Kepler K40m.

Filter size	3$\times$3	5$\times$5	7$\times$7
$W$	32	32	64
$H$	4	8	4
$F_{TB}$	64	32	32
$W_T$	16	8	8
$F_T$	4	8	8
$C_{SH}$	2	1	1

Equations

$$W_{SMB}=n\cdot W_{CD}.\text{vspace}-4pt$$

Full text

††conference: DAC ’17, June 18-22, 2017, Austin, TX, USA

{CCSXML}

<ccs2012> <concept> <concept_id>10010147.10010169.10010170.10010174</concept_id> <concept_desc>Computing methodologies Massively parallel algorithms</concept_desc> <concept_significance>500</concept_significance> </concept> </ccs2012>

\ccsdesc

[500]Computing methodologies Massively parallel algorithms

\printccsdesc

1 Introduction

Convolution is a fundamental operation in many image processing and computer vision applications. For example, image convolution is a key component in numerous basic image processing routines, such as edge detection [1], smoothing [1], template-based object detection [2], etc. Recently, convolutional neural networks (CNNs) [3] have become a powerful deep learning model which has been widely adopted in various computer vision applications, such as image recognition [4], image classification [5], object detection [6], etc. State-of-the-art CNNs typically have quite a few convolutional layers. Propagating through these convolutional layers is always a computation bottleneck in both the training and inference phases of CNNs.

With the rapid development of many-core parallel processors, new methods have been developed by leveraging high computation throughput and memory bandwidth of graphics processing units (GPUs) to accelerate convolution operations. These methods can be roughly classified into four categories: (1) general matrix multiplication (GEMM) based convolution [7, 8], (2) direct convolution [9, 10, 11], (3) fast Fourier transform (FFT) based convolution [12, 13, 14], and (4) the Winograd algorithm [15, 16].

Convolution can be easily converted into a multiplication of two matrices by unrolling all the involved convolution operations [7]. Highly optimized GEMM kernels (e.g., cuBLAS [17]) can be invoked to compute matrix multiplications. This is the default method in Caffe [18], a popular deep learning framework. Although good performance can be attained, it requires a huge amount of additional memory. Recently, cuDNN [8] adopted a GEMM-like method, in which sub-blocks of the input matrices are constructed in on-chip memory at run-time, and thus no additional memory is needed. A direct method was proposed in [9], but the reported performance is not good enough when there are more than 100 channels. In [10], optimization techniques were discussed for direct convolution on GPUs, but the proposed method was not compared with any public library. Cuda-convnet2 [11] also implemented direct convolution on GPUs, but there is no detailed document to describe its methodology or performance. FFT-based convolution [12, 13, 14] can reduce the arithmetic complexity compared with direct methods. However, the filters need to be padded to the same size as the input image, which incurs additional memory and computation time. In addition, in order to reuse the Fourier transform of the filters, the batch size should be big enough. Recent studies have shown that the Winograd algorithm can significantly reduce the arithmetic complexity for the 3$\times$3 filter [15, 16], at the cost of increased memory usage and filter size dependent specialized processing.

Although FFT-based methods and the Winograd algorithm can be faster than direct methods in some cases, they are not universal. Direct convolution is still fundamental and considered the best in general. In this paper, we aim to improve the memory efficiency of direct convolution on GPUs, targeting at two cases: (1) a special case with one input channel, which appears in numerous image processing applications and the input layer of CNNs (for grayscale images), and (2) the general case for CNNs. Specifically, we introduce a general model to address the mismatch between the shared memory bank width of GPUs and computation data width of threads. Based on this model, by carefully optimizing the memory access patterns and computation patterns, we design a communication-optimized kernel for the special case and a communication-reduced kernel for the general case. Experimental data based on implementations on Kepler GPUs show that our convolution kernels achieve 5.16$\times$ and 35.5% average performance improvement compared with the latest cuDNN library, for the special case and the general case, respectively.

2 Problem Formulation

In this section, we present the problem formulation, which illustrates the main challenges and the general model that we propose to overcome such challenges.

2.1 GPU Memory Constraints and Modeling

Most GPU programs are memory bandwidth hungry. GPUs usually have a complex memory hierarchy subject to different constraints. Global memory (GM) accesses should be coalesced in order to reduce latency. Bank conflict should be avoided when accessing the shared memory (SM). For the constant memory (CM), all the access addresses within one warp should be identical to take full advantage of the broadcast mechanism. These are basic constraints that GPU programs should satisfy to achieve good performance.

The SM bank width, which also plays an important role in GPU performance, however, has received less attention from programmers and researchers. We elaborate this problem below. Let $W_{S\!M\!B}$ be the SM bank width. ($W_{S\!M\!B}$ is 8 (bytes) on Kepler and 4 on other GPU architectures.) Further, let $W_{C\!D}$ be the computation data width for each thread. For example, if each thread takes float as the minimum unit fort computation, then $W_{C\!D}=4$. The relation between $W_{S\!M\!B}$ and $W_{C\!D}$ can be described by

[TABLE]

If $n\!=\!1$, the SM bank width and computation data width are matched; otherwise, they are unmatched. Mismatch between $W_{S\!M\!B}$ and $W_{C\!D}$ frequently occurs in practice. Even when $W_{S\!M\!B}$ is 4, $W_{C\!D}$ can be 2 (for short or fp16) or 1 (for char). Fig. 1 illustrates the impact of a mismatch. Consider multiple threads reading from or writing to the SM. A conventional method shown in Fig. 1a is often used, where contiguous threads access contiguous elements, as it is easy to program. But, such a method may fail to fully utilize the available SM bandwidth. For example, if $n\!=\!2$, as shown in Fig. 1a, any two accesses that fall into the same bank have to be serialized. Yet, if we can double the computation data width through intelligent thread layout and computation pattern redesign so that $W_{S\!M\!B}\!=\!W_{C\!D}$, as in Fig. 1b, then each thread can obtain 2 elements together in a single access, resulting in a 2$\times$ improvement in the SM bandwidth.

To demonstrate the importance of this problem, we compare the performance of the single-precision GEMM on a Kepler K40m GPU, as shown in Fig. 2. MAGMA is highly optimized for Fermi and is faster than cuBLAS on the Fermi architecture [19]; however, it becomes 2.4$\times$ slower than cuBLAS on the Kepler architecture. The SM bank width of the Kepler architecture is twice of that of Fermi, causing a mismatch between $W_{C\!D}$ and $W_{S\!M\!B}$ for the MAGMA kernel that operates on float, which results in the loss of half of the SM bandwidth. Yet, a modification to the MAGMA kernel by matching $W_{C\!D}$ with $W_{S\!M\!B}$ saves 36% of the execution time on average.

Consequently, for applications that are sensitive to the SM bandwidth, memory access patterns and computation patterns should be reorganized to match $W_{C\!D}$ with $W_{S\!M\!B}$. That is, each thread should be designed such that it accesses and computes $n$ basic elements as a single unit. In this way, we can obtain an $n\times$ improvement in the SM bandwidth.

2.2 Data Sharing in Convolution

One key aspect in developing highly efficient convolution operations on GPUs is to maximize data sharing, which is also a key factor for communication reduction. Consider the general case of convolution in CNNs (see Fig. 3a). Fig. 3b illustrates a simple data reuse method, in which pixels can be reused in both the horizontal and vertical directions as indicated by the solid circles within the dashed boxes. A simple analysis shows that an input pixel can be used up to $K\!\times\!K\!\times\!F$ times, where $K$ is the filter size and $F$ is the number of filters. This feature should be fully exploited to reduce both GM and SM accesses. For this aim, elaborate memory access patterns and computation patterns need to be used, while still satisfying the basic constraints of the GPU memory hierarchy.

Another challenge arises from the SM bank width model presented in the previous subsection. As convolution is SM bandwidth bounded, when designing convolution kernels on GPUs, we must strive to match $W_{C\!D}$ and $W_{S\!M\!B}$, in order to fully utilize the SM bandwidth.

Taking all these requirements into consideration, the goal of this paper is to develop general solutions for convolution on GPUs, such that (i) the memory communication is reduced as much as possible, (ii) the basic constraints of the GPU memory hierarchy are satisfied, and (iii) the SM bank width and the computation data width are matched.

3 Convolution for Special Case

This section presents our convolution kernel for the special case, in which the input has only one channel ($C\!=\!1$ in Fig. 3a). This case arises at the first layer of CNNs (for grayscale images) and in many image processing applications. We first show how we design the thread layout and then discuss how we achieve optimal memory accesses.

3.1 Thread Layout

The goal of thread layout is to judiciously allocate computation to thread blocks (TBs) and individual threads to maximize both coarse-grained and fine-grained parallelism. Fig. 3b depicts the general concept for our parallelization methodology, in which each thread keeps $K\!\times\!K$ pixels of the input image used by a convolution operation in the register. However, directly applying this method would cause a problem. A thread cannot move right and down simultaneously. In other words, once a thread moves to the right to compute the next convolution, it loses some pixels needed by the convolution below. One could use additional registers to store such lost pixels, but a lot more registers would be required. To resolve this issue, we propose an alternative scheme to maximize both parallelism and 2D data sharing.

To achieve coarse-grained parallelism, we partition the input image into blocks of size $H\!\times\!W$ each (see Fig. 4). Such partitioning enables data sharing along the vertical direction, since one row of the input image can be used by the convolutions of $K$ rows. A TB with $W$ threads handles one image block. Different image blocks are assigned to different TBs so that they are computed in parallel. Each block needs some additional pixels outside its right and bottom boundaries to compute convolutions.

In terms of fine-grained parallelism, all the $W$ threads in a TB compute convolutions of one row in parallel, as illustrated in Fig. 5a. Once a row is finished, the $W$ threads move down to compute the next row. Thus we read a new row from the input image for each down movement. This process continues until reaching the bottom of the block.

3.2 Optimizing Memory Accesses

We now discuss how we schedule memory accesses in coordination with the thread layout design in Section 3.1. Since the filters in the special case (with only 1 input channel) are typically small, they can reside in the CM and no further scheduling is needed. We focus our discussion on accessing the input image in the GM to minimize GM communication.

We first consider the simple case of $W_{C\!D}\!=\!W_{S\!M\!B}$. For each row of the block, we first read it into the SM (including the needed pixels outside the block boundaries), and then the $W$ threads read their corresponding pixels into respective registers. This process allows horizontal data sharing and avoids redundant reads from the GM as adjacent threads share some common pixels. Hence, our 2D data sharing method works as follows: in the horizontal direction, the SM provides inter-thread data sharing; in the vertical direction, intra-thread data sharing is achieved through the private registers of the threads.

A simple analysis shows that each pixel in a block is read from the GM only once, which is, of course, the theoretical lower bound. For the entire image, only those pixels which are needed by a block and outside the block boundaries are read more than once. But, the proportion of such halo pixels is small. As a result, this method is (almost) communication-optimal for GM accesses.

When considering the SM bank width model presented in Section 2.1, the above thread layout and memory access schedule are suitable only when $n\!=\!1$, i.e., $W_{C\!D}\!=\!W_{S\!M\!B}$. On the Kepler architecture, $n\!=\!2$ if we use float as the basic computation unit. Following the general idea depicted in Fig. 1b, we propose to have each thread read, write, and compute $n$ pixels together (using built-in data types such as float2 or float4). Each thread is responsible for $n$ contiguous output pixels in each row and $n\!\times\!H$ output pixels in the block (for one output feature map). The number of threads in a TB is reduced to $\frac{W}{n}$. With this approach, each thread needs a few more registers ($O(K\!\times\!(n\!-\!1))$) to store $K\!\times\!(K\!+\!n\!-\!1)$ pixels that are used by the convolutions for $n$ contiguous output pixels. Fig. 5b shows our modified convolution method specifically for the Kepler architecture.

3.3 Implementation

Algorithm 1 outlines the flow of our special case convolution method at the TB level. The algorithm starts by reading $K$ image rows of the block into the SM (line 1). After that, the first $K\!-\!1$ rows are read into the threads’ registers (line 3). Then convolution is performed on all the rows within the block iteratively in a loop (lines 4-11). For each row, the data are first read from the SM into the threads’ registers (line 6), and then each thread computes convolutions for all the filters (lines 7 and 8). We use a prefetching mechanism to overlap computations and GM accesses. Before the threads’ computation tasks, the next image row of the block is prefetched into the threads’ registers (line 5). Although this operation may take a long time, it can be overlapped with convolution computations, since they have no data dependency. After prefetching is finished, the prefetched data are written into the SM (line 10).

The above algorithm is quite memory efficient. When computing convolutions, the involved pixels are in registers so the latency can be ignored. The filters are fetched from the CM. In our method, all the threads within a warp always compute convolutions using the same filter at the same time, so they always access the identical address, which is the best case for the CM. As the filters are quite small in the special case, a high hit rate of the constant cache can be expected. When accessing the SM and GM, contiguous threads always read or write contiguous addresses (at the granularity of $n$ pixels as a single unit), so both coalesced GM access and conflict-free SM access are achieved. Our experimental results in Section 5 support our analysis here.

4 Convolution for General Case

This section presents our convolution kernel for the general case, where the input has multiple channels (see Fig. 3a). Note that the method for the special case cannot be applied here for the following reasons. For the special case, we can keep the needed pixels in the threads’ registers since the filters are small. Thus, we can finish one convolution at once with $K\!\times\!K$ fused multiply-add (FMA) operations. With multiple channels, the involved pixels of one convolution cannot entirely reside in the registers. Hence, the computation of one convolution must be divided into multiple steps, and the intermediate results should be accumulated in the registers. In addition, the filters (proportional to the number of channels) become larger and may no longer fit in the CM. Instead, the GM needs to store both the filters and the input image. Our basic idea for the general case is inspired by the blocked GEMM method for GPUs [19], but we optimize memory communication by maximally sharing data.

4.1 Thread Layout

Similar to the special case, each input channel of the input image is partitioned into blocks of size $H\!\times\!W$ each (see Fig. 4). We use a 2D TB layout which is similar to that adopted by the blocked GEMM method [19]. Since in the general case, a TB cannot be responsible for all the filters, we divide the computation into a 2D TB layout of size $TB_{X}\!\times\!TB_{Y}$. In the X dimension, a TB is responsible for $F_{TB}$ contiguous filters, where $TB_{X}=\left\lceil{F}/{F_{TB}}\right\rceil$. In the Y dimension, a TB is responsible for $C$ image blocks at the same location of all the $C$ channels. Within a TB, we use a 2D thread layout of size $T_{X}\!\times\!T_{Y}$. Each thread is responsible for $W_{T}$ output pixels and $F_{T}$ filters, where $T_{X}=\frac{F_{TB}}{F_{T}}$ and $T_{Y}=\frac{W\!\times\!H}{W_{T}}$. Each thread keeps $F_{T}\!\times\!W_{T}$ pixels in the register to store the intermediate convolution results for the $W_{T}$ pixels and $F_{T}$ filters. In the following subsections, our discussion will focus on optimizing memory access patterns to reduce memory communication.

4.2 Optimizing Memory Accesses

Fig. 6 depicts our convolution method for the general case at the TB level. To improve the GM efficiency, we store $C_{SH}$ channels of image blocks (including the needed halo pixels outside the block) and the filters in the SM (shown as the two blue boxes in Fig. 6). When reading filter values from the GM to the SM, since the block is transposed, padding (the gray box in Fig. 6) is required for the SM to avoid bank conflict. For the image blocks, pixels are directly read into the SM without transposition, so padding is not needed.

To increase data sharing within a thread, the $W_{T}$ output pixels computed by one thread are contiguous along the horizontal direction. This is a major difference from the blocked GEMM method [19] where contiguous output pixels are computed by contiguous threads. Computing $W_{T}$ contiguous output pixels by one thread involves reading $(W_{T}\!+\!K\!-\!1)\times\!K\!\times\!C$ pixels from the SM, instead of $W_{T}\!\times\!K\!\times\!K\!\times\!C$ if they are computed by different threads. As $C$ may be large, it is impossible to put all of these involved pixels in the register. So we only keep a row of $W_{T}\!+\!K\!-\!1$ input pixels in each thread’s register, and the convolution results are accumulated iteratively. The $W_{T}\!+\!K\!-\!1$ input pixels are used in $K$ rounds of computation of $W_{T}$ output pixels. A round of computation refers to an FMA operation in a convolution.

When reading filter values from the SM, we use the same thread layout as that used in the blocked GEMM method [19]. However, in order to meet the requirement of the SM bank width model, each thread should read $n$ contiguous values as a single unit along the horizontal direction, as illustrated in Fig. 6 (the upper blue box). This method is conflict-free as contiguous threads in the X dimension read contiguous units from the SM.

Our thread layout and memory access patterns also avoid bank conflict when accessing the SM for image blocks. As the X dimension of the 2D thread layout is assigned along the feature direction, $T_{X}$ contiguous threads in the X dimension access the identical address of the SM for image blocks, which benefits from the broadcast mechanism of the SM.

The only issue of this approach is at the writing back phase. Writing the results back to the GM is not coalesced, as contiguous threads in the X dimension compute different output feature maps. However, we have found that in the general case convolution, the writing back phase consumes very little time, so we do not optimize the uncoalesced writing back operations. If one wants to make the writing back coalesced, the SM can be used as a buffer to reorganize the data layout. However, this would lead to additional cost including the SM latency and TB barriers.

4.3 Implementation

Algorithm 2 outlines our general case convolution method at the TB level. The algorithm starts by clearing the results (line 3) and loading the first $C_{SH}$ channels of image blocks and filters into the SM (lines 4 and 5). After that, a loop iteratively accumulates the results for all the channels (lines 7–19). For each channel, a thread needs to conduct $K$ rows of computations (lines 11–15), and for each row, $K$ rounds of computation are conducted (lines 13–15). The image data are loaded into each thread’s register only for each row (line 12), and these data are used by $K$ rounds of computation. The filter data are loaded into the register in each round (line 14). We also use a prefetching method to overlap GM accesses and computation (lines 8, 9, 17 and 18). After the intermediate results of all the channels are accumulated, the final results are written back to the GM (line 20).

Fig. 6 illustrates the first two rounds of computation for thread $(0,0)$. The thread first loads a row of $W_{T}\!+\!K\!-\!1$ pixels from the SM into the register. In the first round, it loads $F_{T}$ filter values from the first row of the SM storing the filters, and then updates the intermediate results by multiplying the $F_{T}$ filter values and the first $W_{T}$ pixels (the green dashed lines). In the second round, the $F_{T}$ filter values are loaded from the second row of the SM, but the pixels are not loaded again, as they are already in the row of $W_{T}\!+\!K\!-\!1$ pixels with an offset (the purple dashed lines).

Compared with direct GEMM-based convolution methods, our method reduces GM communication by approximately $\frac{1}{K}$, since one image row is used by the convolution of $K$ rows. The SM communication for fetching image pixels is reduced by $\frac{W_{T}+K-1}{W_{T}\cdot K}$.

5 Experimental Results

We have implemented our proposed methods and conducted experiments on a Kepler K40m GPU with peak performance of 4290 giga floating-point operations per second (GFlop/s) for single-precision. Our code is compiled with compute capability 3.5. As we aim at direct convolution, we compare our kernels with the GEMM-based convolution provided by cuDNN [8] (version 5.1).

5.1 Results of Special Case

Through design space exploration, we determined that the best block size for the special case convolution kernel is $W\!=\!256$ and $H\!=\!8$. The performance comparison between our kernel and cuDNN for different convolution parameters (image size $N$, filter size $K$, and number of filters $F$) is shown in Fig. 7. For the 1$\times$1 filter, actually there is no data sharing; however, our kernel still obtains an average 6.16$\times$ performance gain, due to the well-designed communication-optimal kernel. For the 3$\times$3 and 5$\times$5 filters, our kernel obtains 6.43$\times$ and 2.90$\times$ average performance gains over cuDNN, respectively. The average performance gain of the three filters we have tested is 5.16$\times$. The performance is lower when $F\!=\!1$, due to the low overlap between communication and computation, as the computation workload is quite low for $F\!=\!1$. However, our kernel can be more than 10$\times$ faster than cuDNN when $F\!=\!1$.

For the 3$\times$3 filter, we have also implemented another kernel in which $W_{C\!D}$ and $W_{S\!M\!B}$ are unmatched, i.e., the basic unit for computation is float. As seen from Fig. 7b, the performance is reduced by 19% if $W_{C\!D}$ and $W_{S\!M\!B}$ are unmatched. It can be expected that the performance degradation will be higher for the general case if $W_{C\!D}$ and $W_{S\!M\!B}$ are unmatched, as the SM is used to store both the input image and the filters in the general case. If we compare the unmatched kernel with cuDNN, even if $W_{C\!D}$ and $W_{S\!M\!B}$ are unmatched, our parallelization strategy is still much better than cuDNN for the special case.

5.2 Results of General Case

Table 1 lists the best configurations of our general case convolution kernel for different filter sizes for the Kepler K40m GPU. The performance comparison between our kernel and cuDNN for different convolution parameters (image size $N$, filter size $K$, number of channels $C$, and number of filters $F$) is shown in Fig. 8. For the three filter sizes we tested, we get 30.5%, 45.3%, and 30.8% average improvements over cuDNN. Only when the image is very small (32$\times$32), our kernel may be a little slower than cuDNN. In all the other cases, our kernel is always faster than cuDNN. The average performance improvement of the three filter sizes is 35.5%. The highest performance we have achieved is 2020 GFlop/s, which is 47% of the hardware peak performance.

6 Conclusions

In this paper, we introduced a general model to address the mismatch between the SM bank width and computation data width of threads. Based on this model, we designed and optimized two convolution kernels on GPUs. By carefully optimizing the thread layout and memory access patterns, we attained 5.16$\times$ and 35.5% average performance improvements over the latest cuDNN library, for the special case and the general case, respectively.

Although we have only implemented our convolution kernels on the Kepler architecture, our proposed ideas can be applied to other applications and architectures. For example, one of the recent development trends of CNNs is to use shorter data types, such as half-precision floating-points and 16- or 8-bit fixed-points, to reduce both the storage requirement and execution time. For these data types, mismatch between the SM bank width and the computation data width exists even for architectures with 4-byte SM bank width. As a result, our proposed model and method will benefit applications using these data types.

7 Acknowledgments

This project was supported by the National Science Foundation under grants 1640081, 1217906, 1629914 and 1617735, and the Nanoelectronics Research Corporation (NERC), a wholly owned subsidiary of the Semiconductor Research Corporation (SRC), through Extremely Energy Efficient Collective Electronics (EXCEL), an SRC-NRI Nanoelectronics Research Initiative under Research Task IDs 2698.004 and 2698.005.