GPU convolutions for neural networks

With all the popularity of deep learning, many researchers in the
field might wonder which framework is "right" to implement their
experiments. For plain neural networks, the main "work horse" is the
matrix multiplication, which can be accelerated a lot using graphics
processing units (GPU). For convolutional architectures, the matrix
multiplication is typically "replaced" by a convolution, and we would
also like to see them being fast(er) on GPU.


Neural net convolutions are somewhat special, since there filters are 3D and pool over input layers. Also, since they are usually applied to
many small "maps" at once, common FFT acceleration techniques do not
apply.


For my own implementations, I compared 3 convolution implementations:

• The convolutions that come with Theano (from git, 2011-1-14). This
  implementation is by far the most flexible, as we will see. It is
  based on the formely separate, now theano-integrated CudaNdarray
  library.
  
• Alex Krizhevsky, a PhD student in Toronto, wrote two publically available convolution routines. We already integrated the first
  version of his convolutions in CUV.
  
• Alex' new convolutions created for the cuda-convnet (svn, 2011-1-13)
  which are described as being "several times faster" than the first
  version.
  

Constraints

The (main) constraints of the three versions are quite different:

Implementation	Image Size	Memory-Ordering (row-major)	Other
Theano	any	(nImages, nChannels, imageH, imageW)	–
Alex old	square only	(nChannels, nImages, imageH*imageW)	nFilters%2==0
Alex new	square only	(nChannels, imageH*imageW, nImages)	nFilters%16==0

Regarding squared images, one can argue that in random image
collections the shapes vary, anyway, and for batch processing it is
necessary to square them.


The ordering is tricky. At first sight, Theano's ordering looks most
intuitive. However, all operations which are functions of all channels
of a single pixel are a bit tricky to optimize. Alex' old and new
orderings can both use efficient matrix-row operations for
cross-channel functions. The "Alex old" convolution has the
disadvantage that images in one batch are not in the columns or the
rows of a matrix, so that final "full" layers (for example in LeNet)
require reordering the matrix. The new convolutions have images in the
columns of a matrix, solving the reordering problem, even though this
ordering looks most un-intuitive.


I should also mention the "sparse" filter option in Alex' code, which
allows to convolve only certain maps with a filter. I'm not going into
detail since Theano does not have this feature and I want to compare
execution times.

Speed

In the following table, all operations were computed 10 times and the
(wall clock) times averaged. For Theano, I varied the 'version'
parameter, but found that the auto-selection (-1) selects the best
algorithm. I used a GTX480 and in an Intel Xeon X5650 (2.67 GHz).

Execution speed of convolution packages

Version	Image Size	Filter Size	Type	Time (ms)	Comment
Naive CPU	32,8,176,176	32,8,7,7	fwd	34200
			dimg	26800
			dflt	n/a
Alex new	32,8,176,176	32,8,7,7	fwd	75
			dimg	90
			dflt	55
			trn	0.3	transposing all input batch
			total	220.3
Alex old	32,8,176,176	32,8,7,7	fwd	101
			dimg	240	plus error padding (3 ms)
			dflt	115	plus summing over batch (.8 ms)
			total	459
Theano	32,8,176,176	32,8,7,7	fwd	268
			dimg	451
			dflt	281
			total	1000

Key:

Image Size

batch size, number of input maps, height, width

Filter Size

number of output maps, number of input maps, height, width

Type

fwd is the "forward pass" convolution, dimg is the
derivative w.r.t. the inputs and dflt is the derivative w.r.t. the
filters.

Discussion: I was quite surprised to see Theano is comparably slow.
It seems that Alex' new convolutions are indeed faster, albeit not
several times (for the tested case) (Update: With patches for small
batch sizes kindly provided by Alex, speed nearly doubled!). The
overhead of a transpose (to comply with the "weird" memory layout) is
negligible compared to the overall advantages. All GPU
implementations significantly outperform a naive CPU version (just
many nested for-loops). Note however that theano is able to generate
code for efficient CPU convolutions.


Combinations: Theano is quite flexible, but "Alex new" is
fast. How do we get the best of two worlds? It is interesting to
note that the memory layouts of both convolutions are transposed to
each other, and that for just 0.3 ms (in the above setting), we can
get from one to the other. So we can get speed or flexibility at
wish.

Maintenance concerns

Both implementations are not particularly very well documented, but
well tested. At least for CudaNdarray, there is a successor on the way. It seems to me that optimized code at this level is mostly
write-only anyway.

Easy Parallelization with C++0X lambda functions, Thread Building Blocks

Lambda functors in C++0X

The relatively new gcc-4.5 release supports lambda expressions,
which – in contrast to boost.lambda, boost.bind and the like – provide
easy capturing of variables in context in the lambda expression.
This finally makes the STL algorithms usable, such as

struct add_val{
  add_val(double d):val(d){}
  void operator()(const double& d){
    d+=val;
  }
  double val;
};


int main(){
  std::vector<double> v;
  // ... fill vector

  double inc = 3.0;

  // the old way of doing things
  add_val av(inc);
  std::for_each(v.begin(),v.end(),av);

  // using a boost.lambda function
  std::for_each(v.begin(),v.end(), _1+=inc);

  // using a c++0x lambda function
  std::for_each(v.begin(),v.end(), [=](double& d){d+=inc;});
}

The "old way" primarily has inconveniences for the programmer:

• We need to define a struct outside the scope, even for tiny
  functionality.
  
• We need to explicitly capture variables from the scope of the
  surrounding code (here: inc) in the struct.
  

Boost.Lambda tried to resolve this problem, in an elegant way I
believe. However, there are still shortcomings of this approach:

• Variables have unintuitive names (_1, _2)
  
• The code in the lambda expression is not really a block of code. It is an expression, where parts may be evaluated
  surprisingly.
  
• Also, since this is an expression, conditionals and loops must be
  expressed in an awkward (that is, non-C++) way using if_, for_
  and so on.
  

The new C++0X lambda syntax gets rid of all these problems. While the
syntax looks a bit strange at first, it is much more readable than
boost.lambda constructs. The block of code is not out of scope, all
names of the surrounding code can be used as copy ([=]) or
reference ([&]).


You might wonder, why we do not just use boost.foreach and get away
with writing

#include <boost/foreach.hpp>
#define foreach BOOST_FOREACH
// ...as before...
foreach(double& d, v){
  d+=inc;
}

// in c++0x, not yet implemented, this will be
for(double& d : v){
  d+=inc;
}

… which is an idiom quite well-known in other languages. The fun
part is, that we cannot change what for does, but we can change the
implementation of for_each. This is one of the things that the Intel (R) Threading Building Blocks library does.

Parallelizing your for-loop by changing two lines

A nice trick now is to change your (side-effect-free) for-loops like
this:

#include <tbb/parallel_for_each.h>

// before
foreach(double& d, v){
  d+=inc;
}

// after
tbb::parallel_for_each(v.begin(),v.end(),[=](double& d){
  d+=inc;
});

… and everything in this loop is automatically run in parallel.

Similar things can be done with OpenMP:

int end = v.size();
#pragma omp parallel for
for(int i=0;i<end;i++){
   v[i]+=inc;
}

but this is already much more intrusive. Furthermore, the loop index
must be an int and the last index must be known. While
variables may be passed as copy (using private (inc)), these
variables are private to the thread, not to the "wrapped" function. Of
course, OpenMP gives you much more fine-grained control over
parallelization as well.

Benchmarking GPUs using the CUV library

Since I started programming GPUs, a few generations of these cards
have been released, most recently the GTX580. I always wondered how
good these actually are and how programs written for one GPU scale to
the next generation.


We have a test suite available here, which is at least of importance
to us: The CUV library. Apart from unit-tests checking correctness of
the implementation, the library also has a few "tests" which measure
execution speed on GPU and CPU for comparison.


Instead of inventing new tests for the benchmark, we simply reuse the
speed tests which come with CUV, as they are probably relevant
use-cases that the programmers optimized for, anyway. A small perl
script that now resides in the scripts/ directory of CUV now
identifies speed tests by their name (*_speed), runs them and parses
their output. We simply collect all values in a defined order and save
them to a file. A larger number of these files can then be analyzed
using a python script included in the same directory, which uses the
superb matplotlib library to draw a bar chart. We use a reference GPU
to compare relative timings.


To cut a long story short, here are the results comparing GTX285,
GTX295, GX2-9800, GTX480 and GTX580:

As expected, Fermi generation cards perform a lot better than the
older generation, the GTX580 also improves on GTX480. Some operations,
which are apparently not implemented well, perform worse with newer
generation cards. The real work horses of our library have definitely
improved a lot over generations, even though we did not spend time on
optimizing them for the later cards.


We're not the first ones to compare these cards of course. Legit Reviews compares the frame rate rendered by the GTX480 and GTX580,
finding only marginal improvements. Brightsideofnews runs many
rendering related benchmarks as well (3DMark Vantage, Unigine Heaven,
Pripyat), with mixed results. However, we're more interested in
general purpose computing (GPGPU) and in writing algorithms for GPU
hardware which then scale with the new hardware, as it becomes
available.