Fast Convolutional Nets With fbfft: A GPU Performance Evaluation

12/24/2014 ∙ by Nicolas Vasilache, et al. ∙ Facebook 0

We examine the performance profile of Convolutional Neural Network training on the current generation of NVIDIA Graphics Processing Units. We introduce two new Fast Fourier Transform convolution implementations: one based on NVIDIA's cuFFT library, and another based on a Facebook authored FFT implementation, fbfft, that provides significant speedups over cuFFT (over 1.5x) for whole CNNs. Both of these convolution implementations are available in open source, and are faster than NVIDIA's cuDNN implementation for many common convolutional layers (up to 23.5x for some synthetic kernel configurations). We discuss different performance regimes of convolutions, comparing areas where straightforward time domain convolutions outperform Fourier frequency domain convolutions. Details on algorithmic applications of NVIDIA GPU hardware specifics in the implementation of fbfft are also provided.



There are no comments yet.


page 6

page 10

Code Repositories


Fast GPU kernels for convolutional networks.

view repo


Auto-tuner that builds optimized kernels for Convolution Layers on ConvNets

view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Deep convolutional neural networks (CNNs) have emerged as one of the most promising techniques to tackle large scale learning problems, whether in image and face recognition, audio and speech processing or natural language understanding. A convolutional layer within these networks provides useful properties such as translation equivariance of activations. A limiting factor for use of convolutional nets on large data sets was, until recently, their computational expense.

Krizhevsky et al. (2012)

demonstrated that training of large CNNs with millions of weights and massive data sets is tractable when graphics processing units (GPUs) are properly put to use. Since then, renewed interest in CNNs insufflated a fresh breath in various frameworks and implementations, including Torch (

Collobert et al. (2011a)

), Theano (

Bergstra et al. (2010)), cuda-convnet (Krizhevsky (2014)

) and Caffe (

Jia et al. (2014)). Many of these frameworks are based around codes for NVIDIA GPUs using CUDA (Garland et al. (2008)).

We discuss our contributions to convolution performance on these GPUs, namely using Fast Fourier Transform (FFT) implementations within the Torch framework. We summarize the theory behind training convolutional layers both in the time and frequency domain in Section 2. We then detail our implementations. The first is based on NVIDIA’s cuFFT and cuBLAS libraries (Section 3). We evaluate our relative performance to NVIDIA’s cuDNN library (Chetlur et al. (2014)) on over different configurations (Section 4). We significantly outperform cuDNN and other time domain convolution implementations for a wide range of problem sizes.

Our second implementation is motivated by limitations in using a black box library such as cuFFT in our application domain, which we describe. In reaction, we implemented a from-scratch open-source implementation of batched 1-D FFT and batched 2-D FFT, called Facebook FFT (fbfft), which achieves over speedup over cuFFT for the sizes of interest in our application domain. This implementation achieves GPU efficiency ratios of over in certain cases. We describe an ongoing effort to further improve the performance of our solution based on algorithmic tiling (Section 6) before we conclude. Our implementation is released as part of the fbcuda and fbcunn open-source libraries at

2 Convolution

Discrete convolution and cross-correlation are used in CNNs. We quickly summarize these and their implementation, with a formulation mirroring Mathieu et al. (2013). Forward propagation (fprop) inputs are a set of input feature planes . These are cross-correlated111Torch practice is that the forward pass is cross-correlation, hence the . with different filter kernel weights , producing output feature planes . Each input and output feature can be part of a minibatch , so we have and :

The feature planes are reduced (summed) pointwise. For back-propagation (bprop), the gradient of the loss with respect to outputs are convolved with the kernels:

Reduction is over here. Finally, the kernel weights are updated using the gradient of the loss with respect to the weights (accGrad):

Reduction is over here. For purposes of this paper, we use set symbols interchangeably to refer to their size: each input plane is a 2-D matrix of size , and each filter kernel is a 2-D matrix of size 2222-D can be extended to -D, .. The output planes are of size , and implement valid-only convolution, as per MATLAB terminology.

Input zero padding

and input mirror padding around the margins of the input () can be optionally added.333Input size , output size .

A popular convolution implementation is to unroll the data until the computation is in the form of a large matrix multiplication (Chellapilla et al. (2006)). This is the strategy followed by many implementors, since matrix multiplication is a well-tuned linear algebra primitive available on virtually any platform. While it is possible to provide instances of direct calculation that are faster than matrix unrolling (e.g., for large , Krizhevsky (2014)), it is challenging to provide an implementation that is faster for more than just a small subset of possible convolution problems.

Introducing strides in this form of convolution (

i.e., performing the convolution at every -th offset) is a popular way to reduce the computational cost at the expense of precision. The memory accesses required are very similar but with fewer reuse opportunities. On the other hand, by the convolution theorem, a convolution of two discrete signals can be performed with lower asymptotic complexity by performing the multiplication in the frequency domain. Applied to the forward pass, it becomes:

where denotes complex conjugation and is the pointwise product. The discrete Fourier basis used is the largest of the two components convolved and the output.444-dimensional or even bigger for performance (Section 3.2). Linearity of the DFT allows one to perform the sum above in the Fourier domain if desired. Applying the FFT then yields a procedure in lieu of the original . Similar transformations apply for the other two passes. We call this a frequency domain convolution, in contrast to time domain convolution via direct computation. Strided convolutions via FFT can be implemented efficiently to obtain good performance Brosch & Tam (2015). We do not consider those in this paper.

3 cuFFT Convolution Implementation

In this section we discuss implementation strategies using the NVIDIA cuFFT libraries and their efficiency.

3.1 FFT Convolution Details

We described the general formulation for the three types of convolutions in section 2. Here, we borrow the Torch naming convention: input for ; weight for ; output for ; gradOutput for ; gradInput for ; and gradWeight for

. All are stored as single-precision floating point 4-D tensors in row-major layout, and are stored in memory using the so-called

BDHW format. This is explicit in the expression , with input image row data as the innermost or most varying dimension.

Table 1 describes the in-order operations for FFT computation of the forward pass, using the and operators and Cgemm matrix multiplication. Similar implementations follow for the other two passes. The G prefix denotes gradients. The F suffix denotes -valued frequency domain tensors; the rest are over . The T suffix denotes transposed tensors.

Table 1: Implementation detail for forward propagation

Exact tensor dimensions are also given above. By taking advantage of the Hermitian symmetry property of the 2-D DFT for -valued inputs we only store about half the complex entries; the remaining can be obtained by complex conjugation. This results in array sizes such as

. We also perform interpolation by

zero-padding, which serves multiple purposes. First, it is necessary to handle boundary conditions.555In this case, we typically have and . Second, it is required to interpolate all operands over the same Fourier basis.666All tensors are zero-padded to before . Finally, padding has an impact on the FFT algorithm used in practice, as well as on the floating point operation count of non-FFT operations (Section 3.2).

Following the conversion into frequency domain, we perform transpositions to prepare the tensors for Cgemm matrix multiplication library calls. The transposition converts the BDHW layout into HWBD. The transposition is currently out-of-place and implemented using the Cgeam routine; we are also considering our own, in-place transposition routine. Cgemm library calls are performed on transposed tensors in the frequency domain. Casting the operation as a Cgemm call allows us to benefit from the heavily tuned cuBLAS routine. Eventually, we transpose the result back into the BDHW format and perform a 2-D inverse FFT. At this point, the resulting real tensor, always , is clipped to the appropriate final size: for fprop, for bprop, for accGrad.

3.2 cuFFT Design Space

We now discuss implementation aspects we explored. Multiple factors influence the computational efficiency of FFTs: transform size ,

’s prime factor decomposition, and whether batched or iterated single transforms are applied. In the deep learning domain, it is commonplace to deal with small sizes,

. If has undesirable properties, efficiency can drop by an order of magnitude.777

cuFFT implements FFTs with the ubiquitous Cooley-Tukey algorithm (Cooley & Tukey (1965)) which takes advantage of trigonometric equalities to recursively decompose and reuse computations. This is further discussed in the Supplement. Decomposition is built on specialized kernels of fixed sizes which correspond to the prime factor decomposition of . cuFFT implements specialized building blocks for radix sizes , and for sizes where , it can use more efficient kernels exploiting the conjugate symmetry property. When does not admit a prime factor decomposition using those radices only, the expensive Bluestein algorithm is used (Bluestein (1970)). Because our results are used in the time domain, we can in fact zero-pad the image and kernel to perform the FFT at any larger size that may be handled more efficiently. Exploiting more efficient, larger sizes should be balanced against the extra cost introduced in the subsequent transposition and matrix multiplication steps. Table 4’s last case is one in which the best tradeoff is not easily guessed. cuFFT also has batched mode optimizations when multiple FFTs of the same size are being performed.

3.3 cuBLAS Design Space

The cuBLAS library also comes with different implementations for batched and single operation modes. We had the choice between implementation options:

  • for larger batches over small matrices, the cublasCgemmBatched library call;

  • for smaller batches over larger matrices, multiple cublasCgemm calls from the host;

  • for intermediate batch and matrix sizes, devices of compute capability 3.5 and higher support dynamic parallelism which allows CUDA kernels to launch other kernels. This can be beneficial for many launches over small matrices.

Note that the discussion above applies to multiplications after transposition. So the matrix size is either , or and the number of such matrices is . Vendor libraries are usually optimized for throughput and not latency, so we expect it to be more efficient for larger sizes along critical dimensions (i.e., image size for the batch case and , or for the multiple kernel case). Due to build system limitations we were not able to experiment with the dynamic parallelism strategy; we leave this for future work.

At the system level, we use CUDA streams and buffering of all CUDA resources and intermediate buffers to remove synchronization points across convolutions. We are mindful of memory consumption; to address this we keep one single buffered copy of each type of tensor involved. This behavior is tailored for a bulk synchronous execution of layers on a GPU and is not adapted for multiple asynchronous convolutions on the same GPU. The buffers are automatically expanded as required and reused as much as possible.

3.4 Autotuning

We combine the above implementation with a simple autotuning strategy. We devise a strategy selection mechanism that runs once for each problem size and caches the fastest strategy out of a few dozen for later reuse. The autotuning strategy explores different possible Fourier basis sizes that can be decomposed in powers for which cuFFT has an efficient implementation. In other words, for an FFT dimension of size , we explore the sizes where . When the input size is a power of , the search space is reduced to a single point. In addition to Fourier basis sizes, we weigh in various cuBLAS calls and asynchronous modes.

4 cuFFT Convolution Performance

4.1 Performance versus cuDNN: 8,232 configurations

We compare our cuFFT convolution results against NVIDIA’s cuDNN 1.0 library (Chetlur et al. (2014)), which contains one of the fastest, general purpose convolution methods for the GPU, using matrix unrolling. It has decent performance for many problem sizes thanks to heavy autotuning of cuBLAS codes for different problems. It is a strong baseline for this reason.

Image CNNs to date have for the most part used square input images and filters, though rectangular filters are valid for other problems (notably text CNNs, Collobert et al. (2011b)). Thus, we restrict ourselves to a 5-D problem domain . Much of this space is not used in practice. Some areas are perhaps over-emphasized (large , small ) due to current engineering concerns. We evaluate cuDNN vs cuFFT-based convolution for Table 2’s configurations.888Parameterized on output rather than input size because the implied will be valid for any choice of .

Minibatch size () , , ,
Input filters () , , , , , ,
Output filters () , , , , , ,
Kernel h/w () , , , , ,
Output h/w () , , , , , ,
Table 2: Configuration elements evaluated

Figures 6-6 are performance summaries of cuFFT convolution versus cuDNN on a NVIDIA Tesla K40m, averaged across all three passes. The -axis problem size corresponds to the minibatch size multiplied by number of input and output planes (); each one of these is a pass reduction dimension. Many possible combinations of may map to the same problem size. cuDNN performance varies to a greater degree than cuFFT across passes. This is due to the asymmetry of convolution sizes in each pass, and the fact that a larger convolution kernel (as seen with gradient accumulation) is essentially free in the Fourier domain. Averaging the three passes together provides a proxy for overall performance. The -axis corresponds to output height/width. For deeper layers in image CNNs, output size will decrease while will increase, so depth corresponds to moving from the upper right to the lower left of the graph. Black areas in the chart are due to failed cuFFT runs, due to memory pressure or undetermined potential cuFFT 6.5 issues.

FFT convolutions make large kernel sizes inexpensive, which make the performance of all three passes roughly equal (Table 4). On the other hand, zero-padding to penalizes smaller kernels compared to cuDNN. For kernels (Figure 6), cuFFT performance is poor compared to cuDNN. The overhead of multiple kernel launches, streaming memory in and out multiple times, and zero-padding to the input size often outweigh the algorithmic advantage of FFT. However, for the largest problem sizes, convolution via FFT can still be advantageous, with top speed faster than cuDNN. kernels (Figure 6) show an increasing dominance of the FFT strategy, with top speed faster. The tendency is confirmed for larger kernel sizes: at , maximum speedup is over cuDNN.

Figure 1: kernel (K40m)
Figure 2: kernel (K40m)
Figure 3: kernel (K40m)
Figure 4: kernel (K40m)
Figure 5: kernel (K40m)
Figure 6: kernel (K40m)

4.2 CNN Performance

In table 3, we show performance for real CNNs, AlexNet (Krizhevsky et al. (2012)) and OverFeat fast (Sermanet et al. (2014)), comparing against cuDNN and cuda-convnet2 (ccn2) kernels in Torch. The first layer uses cuDNN for the cuFFT runs because it is strided, but all other layers use cuFFT. The timings include all convolutional layers of the network.

AlexNet cuFFT
OverFeat fast cuFFT
Table 3: AlexNet and OverFeat fast performance (K40, )

Table 4 shows the performance of the cuDNN and our cuFFT convolution implementation for some representative layer sizes, assuming all the data is present on the GPU. Our speedups range from to over cuDNN. Unsurprisingly, larger , smaller all contribute to reduced efficiency with the FFT. More surprisingly, we experience noticeable speedups on small kernels as long as the input tensor remains of small size. The optimal FFT sizes that autotuning finds are reported in columns and ; note L5 padding being found by the autotuner. Column has the trillion equivalent time-domain reductions per second (single-precision floating point multiply-adds) achieved by our implementation on a NVIDIA Tesla K40m on CUDA 6.5. This number represents the throughput a time-domain kernel needs to achieve in order to match our implementation; it is computed as . This is a metric to compare relative efficiency across problem and padding sizes. In the cases L2, L3 and L4, a time domain convolution would need to exceed the K40m peak of Tflop/sec in order to match our throughput.

L1 Params:
L2 Params:
L3 Params:
L4 Params:
L5 Params:
Table 4: Representative layer performance (, K40m)

5 Implementation

This section presumes familiarity with GPU architecture. Refer to the Supplement for details.

When designing high-performance libraries, multiple objectives must be balanced against each other: memory latency/bandwidth tradeoffs, maximizing locality without sacrificing too much parallelism, good instruction mix, register usage and mapping strategy of computation and data to memories and compute elements. A key principle is to design a set of leaf kernels with well-tuned in-register performance and reduce the larger problem to a combination of these kernels by data and loop tiling (Irigoin & Triolet (1988)) and recursive decompositions (Gunnels et al. (2001)). Since vendors have to sustain high performance for a large class of application domains, there exist parameter configurations for which a carefully tuned approach significantly outperforms vendor-tuned libraries (Shin et al. (2010)). For common deep learning use, convolutional layers consist of many batched small 2-D convolutions. These are tiny relative to DSP and HPC standards and put us in a regime where (a) we fall outside of the highly tuned regime, (b) feature dimensions are often smaller than GPU warp sizes and can often fit exclusively in registers rather than in shared memory (SMEM), and (c) we are very sensitive to latencies. We determined that it is possible to obtain better efficiency than the existing batched cuFFT mode for CNNs.

5.1 Limitations of cuFFT 

Because the cuFFT library is a black box, zero-padding999This is different from the FFTW compatibility padding mode for in-place transforms. has to be explicitly embedded in the input and output arrays. The consequence is that one may need to allocate a duplicate, larger memory region (only once) and copy data from non-padded tensors to padded tensors. This memory consumption and spurious copies affect latency significantly. Instead, we devised an implementation for batched 1-D FFT and 2-D FFT of sizes 2-256 and reaches up to 78% efficiency at 97.5% occupancy. We also implemented an IFFT kernel based on our FFT kernel.

In our implementation we use clipping to conditionally load a value if reading within bounds or a constant () otherwise. This is an approach used in automatic code generation tools such as Halide (Ragan-Kelley et al. (2013)) and relies on aggressive if-conversion properties of the CUDA compiler. It allows for more efficient control flow rather than using explicit loop prologues and epilogues. This mechanism does not require any additional memory allocation and is zero-copy; this is particularly desirable in the latency sensitive mode.

Additionally, since cuFFT and cuBLAS are closed source, it is impossible to take advantage of algorithmic simplifications that may be available. For instance, in the forward pass of our computation as shown in Table 1, the result of the first cuFFT call is of the form . With we return it in the form where the two innermost data dimensions are transposed. This allows us to remove a full data transposition from each of the FFT kernels. Another domain-specific optimization we have yet to explore is eliminating bit reversal portions of the FFT and IFFT. This can be done by performing the FFT with decimation in frequency (DIF) and the IFFT with decimation in time (DIT), discussed in the Supplement.

5.2 Warp-level 1-D FFT and 2-D FFT for size

For batched FFT of power of two sizes we view a single warp as a small distributed system with lockstep collective communication capabilities and we program it in a bulk-synchronous fashion (Valiant (1990)). We implement DIF and enforce the following invariants for the steps:

  • each warp thread originally loads one real element of the input vector and locally computes one complex twiddle factor (i.e. a root of unity);

  • at each step, all warp threads exchange data with another thread in the warp in parallel and produce a new value;

  • then, all warp threads exchange twiddle factors with another thread in the warp in parallel, and produce a new value.

The two bulk-synchronous exchanges can be written each with one warp-wide instruction. After the steps, the FFT is computed and stored in a distributed and bit reversed manner within register across a warp. For sizes , bit reversal can be implemented with a single warp shuffle.

We either load twiddle factors from device memory or compute them with the sincosf function only once, and subsequently swap them within registers. This greatly reduces the reliance on either memory bandwidth or on the special functional unit at the expense of a few additional registers. The decision between explicitly loading twiddle factors from device memory or computing them is a tradeoff between arithmetic intensity and memory bandwidth. For sizes and the arithmetic pipeline is the bottleneck. Loading twiddle factors from memory for these two special sizes results in a performance increase of and respectively.

The discussion above applies to 1-D FFT and to each independent FFT within a larger 2-D FFT. A -D Fourier transform is separable and can be implemented with sets of multiple 1-D FFT with transpositions between each of these sets. In 2-D FFT -to-, the first set comprises FFTs and the second set comprises FFTs by Hermitian symmetry. Following standard techniques Lyons (1996) we further pack real FFTs into a single complex FFT . The extra term in the quantity makes the computation ill-balanced and can bring down performance by lowering occupancy. We chose to dimension our kernels to have size and introduce additional control flow to handle the border case. This results in additional performance. We implement the transposition in SMEM across warps following Ruetsch & Micikevicius (2009). Data is already resident in registers so our main concerns are limiting SMEM usage to keep occupancy high, and limiting load/stores by using vector instructions to avoid saturating the load-store unit (LSU).

5.3 1-D FFT and 2-D FFT for size

With size as our building block, we extend our strategy to larger sizes. We use the same single warp approach to compute a full 1-D FFT. The main difference is that the computation is now distributed across multiple registers across threads in a warp ( Fourier coefficients and twiddle factors in registers per thread). Because we perform a full FFT per warp, a performance cross-over where cuFFT wins happens after register usage limits occupancy too much. We outperform 1-D cuFFT for , with a hard register limit at ( and similarly for 2-D FFT). This is still well within our application domain. The following modifications handle multiple registers per thread:

  • Hermitian symmetry allows us to perform half the computation. There is a tradeoff between adding control-flow divergence and performing less work. At , benefits from reduced computations dominate divergence losses;

  • we take advantage of trigonometric symmetries and twiddle factor distribution to compute only a fraction of the roots of unity needed for each FFT, distributed with register to register copies;

  • twiddle factor re-balancing across a warp and across registers requires a different implementation. We managed to implement it fully within registers;

  • bit reversal occurs across registers and across warps. The high-order bits represent the register while the low-order bits represent the warp. Without a sophisticated implementation, this results in indirect addressing of registers which is costly. We implement a simple bit reversal in SMEM, which is an occupancy bottleneck at for 1-D FFT.

In the 2-D FFT case, the intermediate transpose becomes significantly more expensive. We experimented with various strategies to keep occupancy high, including partial transpositions within a warp to use minimal amounts of SMEM.

5.4 Discussion

We report the relative performance of our implementation compared to cuFFT for various batch and input sizes of interest. The number of batches to consider depends on the dimension of CNN layers as well as any multi-GPU parallelization strategy that may be involved. At typical sizes of interest, is between and faster. We tried up to million batches and at larger sizes gains stabilize around but efficiency goes down as more and more memory is used.

Figure 7: fbfft-1D FFT and IFFT (K40m, cuFFT 6.5 @ 1x)
Figure 8: fbfft-2D FFT and IFFT (K40m, cuFFT 6.5 @ 1x)

Figure 8 shows the performance in the 1-D case. These numbers do not exercise our implicit zero-copy padding, so we expect additional gains when we incorporate our FFT in the convolution. Our implementation outperforms cuFFT for all cases of interest, more dramatically so for smaller batch sizes. Small batch sizes also correspond to the latency sensitive regime in Figures 6-6 for which the cuFFT based implementation performs quite worse than cuDNN. We achieve efficiency at occupancy for size at batch size , as reported by nvvp.

Figure 8 shows the performance in the 2-D case. Relative performance gains for sizes are more modest than in the 1-D case, even losing to cuFFT at size and small batch sizes. The magnitude of the relative gains at various batch sizes drops faster than in the 1-D case. Looking at the performance of the FFT, we obtain speedup over cuFFT at batches. The same ratio is not obtained until batches in 1-D FFT.101010This is not unexpected because these two computations perform the same number of flops when accounting for Hermitian symmetry, plus the fact that the efficiency of cuFFT increases while remains high but almost constant. When coupled with the tiling strategy in Section 6, we emphasize that the sizes of interest are actually -, and depend on but not input . Batch sizes can vary on the whole spectrum.

We interfaced into our convolution module and ran experiments with kernels for the different convolution passes over inputs of sizes . For problem size, we used . By swapping our FFT implementation we observed an overall mean speedup of

with standard deviation

and geometric mean

The minimum speedup was , despite sometimes performing more computations with which can only interpolate to a power of . These experiments exercise the zero-copy padding and lower memory footprints of compared to cuFFT but do not yet reflect additional optimizations such as tiling and bit twiddling elision.

6 Current Limitations and Future Work

In our current implementation, heavily relies on shuffle instructions. In spite of a good efficiency, we only utilize of the available memory bandwidth. This is due to the load and store instructions in our kernel competing with the shuffle instructions for the Load-Store Unit (LSU). As a consequence, our first bottleneck is the number of instructions issued on the LSU. For instance, on Kepler (capability ), the throughput for -bit floating point multiply-add operations is per cycle but the throughput for shuffles is only . In the future we will investigate and release faster implementations as they become available.

Temporary memory overhead requirements are a common issue when performing convolutions in the Fourier domain. In this first implementation, we introduced the following memory buffers to support our implementation:

  • for each of input, output and weight tensors we store buffer for the frequency array and buffer for its complex transpose. These buffers store the Fourier representation and are generally limited by the weight tensor which is independent of the mini-batch size. Because of the global memory pressure we introduce, we reuse buffers at each layer and pass on the opportunity to (1) reuse 2 FFT results in each hidden layer, reducing the cost of forward FFTs by ; and (2) asynchronously precompute FFTs of the weight tensors and their gradients to better fill the gpu utilization pipeline,

  • when using cuFFT we additionally pad the input, weight and output tensors explicitly to the best performing common fft size

  • when using cuFFT additional temporary memory is reserved by each cufftPlan

  • with padding is implicit but and no temporary memory buffer is needed until we reach size . On the other hand, only supports square convolutions whose size is a power of . As a consequence, too much padding could occur and adversely affect both performance and memory consumption. The tiling strategy we describe next is a good way to circumvent the problem.

Additionally, we recently developed an in-place transposed batched CGEMM which permits the removal of the complex transposed buffer. For this problem, a tool like MaxAS Lavin (2015) could be valuable.

provides the most gains over cuFFT at sizes -. A tiling strategy for the input can be used to exploit this advantage. When the kernel is significantly smaller than the input, we can decompose a large convolution into several smaller ones. For simplicity, we consider 1D convolution on a single input plane, as it can trivially be extended. Let be an input of size , a kernel of size and . We write for the vector formed by contiguous elements of : . Let . From the definition of the convolution, we have:

So the convolution of the input of size can be computed with convolutions with inputs of size . The cost of the convolution goes down from to . From this formula, we see that the optimal is of the order of , to get the complexity . This strategy allows us to speed up forward and backward propagation. Tiling can also be used to reduce memory cost for temporary storage by not running all the tiles in parallel (just the tiles which do run in parallel need their scratch space), at the potential expense of parallelism or efficiency.

For the gradient accumulation, we cannot reuse this strategy, since it involves a larger convolution between an input of size and a kernel of size . However, we have a similar formula:

And so

We have a few other optimizations that are planned as well. Since much of the data we have is already available in registers or in shared memory, we are implementing our own in-place, in-register transpose via recursive decomposition. The pointwise multiplications in the Fourier domain, especially with tiling, are rather small, so our own matrix multiplication routines integrated with the rest of the convolution kernel code might win over cuBLAS, and prevent the need for multiple CUDA kernel launches and their associated overhead. Finally, as mentioned earlier, bit reversal portions can be eliminated with the FFT using DIF and the IFFT using DIT.

7 Conclusion

To summarize, we achieve significant gains in CNNs using FFTs, with a cuFFT convolution implementation achieving speedups over cuDNN for common sizes. In reaction to cuFFT and cuBLAS limitations in the context of our specific application domain, we developed our own FFT implementation, fbfft, which is more suited to deep learning problem sizes (large batches, small feature planes). itself is faster than cuFFT transforms for these problems of interest. For convolution, it is faster than the cuFFT as well, with a mean of for sizes that we wish to exploit.

Given our new efficient primitive for size - convolution, we are continuing work on bit twiddling, transposition and pointwise multiplication optimizations, and continuing work on tiling to make the computational advantage at that size apply to larger convolution problems. These will all allow for reduced training time and use of ever larger and deeper CNNs.


We would like to thank Julien Demouth from NVIDIA who suggested further improvements are still possible by virtue of the current implementation being LSU throughput-bound rather than memory-bound.


8 Supplement

8.1 cuFFT Convolution Performance Breakdown

We show a breakdown of cuFFT convolution performance for the steps indicated in Table 1. The timings do not add up to of the reported performance in the previous table because we do not report additional copies needed for zero-padding here. We also enforce force extra synchronizations to isolate the contribution of each operation. Abstracting from these details, the FFT and IFFT take up a significant amount of compute resources, which we address in Section 5.

Table 5: cuFFT convolution performance breakdown (K40m, )

In the particular case of L1, the FFTs take more than of the runtime. This is due to the wasteful interpolation of the kernel tensor from a up to , which is the minimal size to compute the FFT of the input array without interpolation loss. In such cases, the tiling strategy we are developing (see section 6) will result in large additional performance gains.

8.2 FFT : Decimation in Time vs Frequency

A Fourier transform projects and -valued functions onto a harmonic orthogonal basis. The discrete Fourier transform of a vector is the vector:

where is the

-root of unity. The traditional radix-2 Cooley-Tukey algorithm recursively decomposes the computation between an odd and even part:

This decomposition is called decimation in time (DIT). An alternate decomposition performs decimation in frequency (DIF):

When is a power of , both decimations recursively decompose into a perfectly balanced tree and take advantage of the symmetry properties of the roots of unity. The dataflow graph for the radix-2 FFT has a butterfly shape and is a good way of visualizing the computations. There is a symmetry between DIT and DIF in both the order of operations applied and in whether the input or the output order is shuffled (Figure 9).

Figure 9: DIT output ordered (left); DIF input ordered (right) (Burrus (2008))

8.3 GPU Programming

There are a variety of references available that describe CUDA and NVIDIA’s various GPU architectures (Garland et al. (2008)) which we won’t discuss in detail, but the implementation of very much depends upon specifics of the Kepler GPU architecture.

NVIDIA GPUs execute code at the granularity of a warp which is defined as a set of threads in all existing architectures; each thread is assigned a lane within the warp. These threads execute in a SIMT (single instruction, multiple thread) fashion, meaning that a warp is an atomic unit of execution. It holds a single program counter (PC) and can thus only execute a single instruction at a time across all of its threads. Collections of warps are brought together in blocks or CTAs, which together share a region of fast shared memory resident on chip. Blocks themselves can only exchange data via much slower global memory, resident on the GPU or in the host CPU’s address space.

Individual threads within a warp are free to take divergent paths, but since a single PC is present, each branch in the execution will be serialized. Threads that aren’t participating in the branch in question are disabled. In other words, if all threads were to take divergent code paths, we would obtain only of the computational efficiency.

Divergent code paths are hard to avoid, but the NVIDIA instruction set has means to reduce their cost (Giles (2014)). One is with predicated instructions, which are used for small branches, in which all warp threads execute both parts of the branch, with non-participating threads having no side effects.

Block threads have access to a register file, with up to registers per thread for Kepler. Registers are allocated statically by the CUDA compiler. An important performance factor when writing CUDA kernels is that data should be kept in registers as much as possible to avoid communications. Registers in CUDA are “addressable”: it is possible to declare a static array within registers and operate on its elements. The limitation is that all addressing should be performed using statically determined constants so the compiler can translate these accesses to a register number known at compile time. Indirect addressing is also supported but results in copies to a local region within global memory, which essentially constitutes register spilling. Even with the presence of caches, using local memory usually comes with a performance hit.111111There are bleeding edge cases where a little local memory consumption helps performance; for instance, when restricting the number of registers per thread to increase occupancy. As a consequence, we design our kernels using aggressive inlining, template parameters and unrolling directives to make all register accesses statically addressable.

The Kepler architecture introduced specialized shuffle instructions to exchange data between registers within a warp synchronously, which avoids round-trips to shared or global memory. Interestingly, these shuffle instructions allow the dynamic indexing of an array held in registers, as long as the array is distributed in a cyclic fashion across registers in each thread within a warp.

  float arr[3];
  // This simulates a linear array float realArr[96]:
  // arr[0] holds elements 0-31 (lane i holds element i)
  // arr[1] holds elements 32-63 (lane i holds element 32 + i)
  // arr[2] holds elements 64-95 (lane i holds element 64 + i)
  // Example: all warp threads read value held at realArr[34]
  float val = __shfl(arr[1], 2); // ‘1‘ must be statically known
                                 // ‘2‘ can be dynamic

Many warps run in parallel and can be switched by the GPU hardware at each cycle. When enough parallelism is available (measured in occupancy of the GPU as a first approximation), long latency operations are hidden thanks to fast context switching. Registers and shared memory come in finite quantities on each GPU compute multiprocessor. These limited resources are partitioned by the compiler and the hardware amongst computations at the level of a CUDA kernel. Increased usage of registers or of shared memory can reduce GPU occupancy, which limits the ability to hide long latency operations. Reduced occupancy does not necessarily result in performance loss (Volkov (2010)). There are often non-obvious performance tradeoffs in increasing or decreasing threads per block, shared memory per block or registers per thread that are hard to discover. This problem is one of the many reasons why designing a one-size-fits-all implementation that aims to be efficient for any problem is difficult.