A short note on automatic differentiation

Do you remember your undergraduate calculus course? Frankly speaking, I don’t. But I do remember a big blackboard completely covered in chalk scribbles. That was a nightmare! Hopefully, those times are gone. I mean, why would we care about derivatives now?

Here at ELEKS, we do care about derivatives. We create and deliver precise models of various dynamic processes. The application field doesn’t matter: it may be physics, finances, or biology. What really matters is once you need to simulate something – most likely you’ll need to compute derivatives. And you need to compute them fast and with great precision. Thus, we created ADEL – a C++ template based-library that fits our needs. But let’s start with some background information.

The requirements listed above do not allow for solutions such as symbolic derivation or finite difference schemes. The first approach generates the exact formulae for derivatives of all levels and directions. This is a very memory/time consuming method. It produces enormous formulae. While the second one computes the numerical approximations of the value of a derivative. These suffer from various round-offs, cancellation and discretization errors.

Fortunately, there is a middle ground – automatic differentiation (AD). The way it computes derivatives is essentially through a chain rule. It generates evaluations (but not formulii) of the derivatives. As the result, there are none of the drawbacks(accuracy loss) of numerical differentiation where the step size needs to be carefully chosen to avoid errors. All intermediate expressions are evaluated as soon as possible; this saves memory, and removes the need for later simplication. Moreover AD is not as complicated as symbolic representation.

There are two modes of automatic differentiation – the forward mode and the reverse mode – both use different directions of the chain rules to propagate the derivative information. There are two main implementations of automatic differentiations:

Operator overloading is straightforward and easy to understand – one overloads the operators not only to compute function values but also to compute derivative information or build up the computational graph. 

Source transformation works like a compiler, which reads in the code that computes the function, analyzes the code and exports a source code that computes the derivatives. This implementation normally generates more optimized programs but involves a lot of implementation efforts.

ADEL contains the implementation of operator overloaded versions of both forward and reverse modes. This is an open source project; one may find all the required materials at: https://github.com/eleks/ADEL.

The main feature of ADEL’s forward mode is careful handmade optimization. Template based implementation with loop unrolling is compiled into very efficient sequential code. This implementation suits GPU computing architecture very well. The support of CUDA makes the solution truly unique.

The reverse mode of the ADEL library is special in the way it stores the computational graph. It uses the stack of the application for this purpose. Each overloaded operation creates the data structure that is placed on the top of the stack; the life-time of these structures is just enough to compute the derivative data. This makes the algorithm extremely efficient since the most of the operations could be inlined.

We compared our implementation with some of already existing solutions in order to see how it stands up against the others. We implemented Newton’s method for testing the accuracy of our tools. For testing the performance, a random expressions generator was used. It is capable of generating some monstrous things (your professor would never dream of such things). Here is an example of a test function:

template<typename ADType>

ADType TestFunction(const ADType(&x)[6]) {

  y[0] = x[1]; y[1] = -2.730; y[2] = 4.555;

  for (int i = 0; i < 1000; i++) {

    ADType y[3];

    y[0] = x[1]+x[5]+4.624/x[0]+(x[3]-y[0])/x[4]+(atan(-y[2])-x[2])*y[1];

    y[1] = (sin(y[1])-y[2])/x[2]/x[4]*(x[3]+log(x[0])/(y[0]-x[5]))+x[1];

    y[2] = atan(x[5])*(-y[2])/x[3]/(x[4]/(x[2]/x[1]+acos(y[1])))*y[0];

  }

  return 2.658/log(x[0])*(log(y[2])/x[4])*(x[1]+x[5])*(y[0]/y[1])*x[3];

};

The main characteristic of the test case is the number of independent (x) and dependent (y) variables as well as total number of functions/operators. The given task was to compute the gradient and hessian of the return value with respect to all independent variables. In the example above, we have 6 independent and 3 temporary (dependent) variables that produce a single output. 

We tested the following third-party AD implementations:

FADBAD (http://www.fadbad.com/fadbad.html) – a lightweight, user-friendly library. We used a combined mode: the gradient was computed via the reverse mode, while the hessian – via the forward mode.

ADOL-C (http://www.coin-or.org/projects/ADOL-C.xml) – a huge math package that is capable of computing literally everything. The reverse mode was used in the experiments.

There were four test cases of differing sizes with 5/10/20/30 input and 1/5/10/15 intermediate variables. The results were normalized, in such way that the slowest execution was scored with 100 time units, therefore less is better.

Looking at the results one can see that we did not create an ultimate AD tool - not yet. In smaller tests both ADEL modes greatly overrun their rivals. When the number of variables increases, the performance of ADEL [Forward] degrades because of the large memory overhead per variable. ADEL [Reverse] shows the best average performance; FADBAD comes second; ADOL-C received no prize at all.

Automatic differentiation is a very promising but extremely underused tool. The majority of the community still does not recognize it. This article is just a small but solid step into the bright future of AD tools. As a final word, we encourage you to download the source code and to experiment with ADEL on your own. If there is something we can improve, feel free to contact us.

Introducing Intel Xeon Phi - hardware and installation

Few days ago we have got Intel Xeon Phi 5110P Coprocessor - 60-core x86 accelerator that is direct competitor of GPGPU on HPC market and is based on Intel Many-Integrated Core (MIC) architecture.
As Xeon Phi is a quite new piece of hardware there are very few information in the Internet regarding its usage, installation and compatible hardware (especially motherboards). We would like to fill that gap and write about our experience of assembling custom server with Xeon Phi installed.
Technically, Xeon Phi looks very much like high-end 2-slot PCIe graphics card, except of the fact that it doesn't have DVI output on its back. We have 5110P model with passive cooling, hence there are models with active cooling. Without additional cooling there might be some temperature issues, so make sure that you have a case with additional cooler. Our server case has additional cooling installed and is placed in server room with powerful air conditioning, so no problems for us there.
Not every single motherboard can support this kind of accelerator. The trick is that it has to support 64-bit PCIe addressing and you can find this feature only on modern high-end motherboards. As Intel doesn't have official list of supported motherboards it could be quite tricky to choose motherboard for your Xeon Phi server, in case you buy it separately and assemble by yourself. We managed to run the accelerator with ASUS P9X79 WS motherboard (thanks ASUS for clear information that it supports Xeon Phi). Our motherboard had outdated version of BIOS, so we had to flush the latest version and change settings there in order to enable Xeon Phi support. If you visited a link above you might have noticed a statement that it supports only Xeon Phi 3100 series with active cooling, but in fact it works perfectly with 5110P, just don't forget about additional cooling in your case.
At peak loads Xeon Phi 5110P consumes 225W of power. 3100 series devices use even more, up to 300W, so make sure that your power-supply unit is able to feed it with that amount of energy.
Officially Xeon Phi works only with RHEL and SLES, but in our case it works perfectly with CentOS (open source clone of RHEL).
We are going to write more posts about Xeon Phi very soon. Stay tuned!

Speeding up MR Image Reconstruction with GeneRalized Autocalibrating Partially Parallel Acquisitions (GRAPPA)

This time we would like to share some details about realization of one of our projects in medical imaging.
Currently one of the bottlenecks in MR image reconstruction is speed improvement. Improving the speed of image reconstruction is difficult from algorithmic point of view. But it’s becoming more popular to improve algorithm performance using GPU.

Introduction

In magnetic resonance (MR) image reconstruction raw data measured from the scanner correspond to the Fourier coefficients of the target image, so the fast Fourier transform (FFT) is used for reconstruction. An important procedure called the density compensation is necessary in the very beginning to account for the non-uniform sampling. GeneRalized Autocalibrating Partially Parallel Acquisitions (GRAPPA) is a partially parallel acquisition (PPA) method which is used to reduce scan times. In this case only partial data is acquired and missing data is mathematically calculated from available data.

Subject for optimization

Original implementation was based on FFTW library for Fast Fourier transforms and adapted Singular Value Decomposition (SVD) algorithm from ACM Collected Algorithms for GRAPPA preprocessing. These two algorithms are the most computationally intensive parts of the whole image reconstruction. FFTW library is claimed to be the fastest CPU implementation using all possible optimizations like SSE2/3 and hyper threading, however it does not leverages the power of modern GPU cards. SVD algorithm was done on CPU as well. It is known to be badly parallelizable for small matrices, but in case of GRAPPA algorithm we have many image frames with same size which can be processed in parallel. Besides there are many intermediate steps which consume a lot of CPU and they can be easily parallelized on GPU.

Technical analysis

FFTW library performance is comparable with Intel MKL implementation. NVidia provides comparison for their CUDA based cuFFT library with MKL implementation (Figure 1):

Figure 1 Comparison of CUDA based cuFFT library with MKL implementation

According to this we should achieve up to 8-10x faster FFT processing  when using GPU accelerated cuFFT library from Nvidia. GPU accelerated SVD algorithm is also available, for example CULA library by EM Photonics. However, current CULA library implementation does not support batch mode, so we will need to process all image frames as a sequence. Brief testing showed that 64 image frames (256*256) are processed even slower than CPU based version. Since we haven’t found any good alternative to CULA library we decided to implement our own GPU accelerated version of SVD algorithm.

Implementation

FFT part of image reconstruction when using cuFFT library was straightforward, however we had to deal with image frames which does not fit into available GPU memory. We had to write algorithm to run FFT over portions of the large data frame with subsequent aggregation. Figure 2.1 below shows case when all data fits into GPU memory.

Figure 2.1 

Figure 2.2 illustrates the case when huge data is processed. Solid lines in figure below show measured performance, dashed lines show estimated time in case data fits into GPU memory.

Figure 2.2

Much more interesting was to implement GPU accelerated SVD algorithm with batch mode. All implementations we had found are focusing on maximum parallelization of a single SVD run, hence we had to change approach. Basically SVD algorithm consists of HouseHolder Reduction, QR Diagonalization and Back Transformation steps. All are iterative processes when next step depends on results from previous step. In case of small matrices each CUDA kernel can’t effectively utilize all parallel processing units of modern GPU. So we had to write kernels in a way when every iteration for all matrices is processed by a single kernel run. This way in case of 64 matrices with 128x128 size each we can process 64*128 elements at a time instead of 128. Figures 3.1 and 3.2 show performance comparison for CULA Library and our implementation. 

Figure 3.1

Figure 3.2

With more than 8 frames per batch our implementation shows much better performance comparing to sequential CULA calls, although it is not so efficient for a single frame.

Results

As a result we have developed a pure C++ library with a set of specialized functions which perform various stages of image reconstruction. It requires only CUDA runtime libraries and free cuFFT library provided by NVidia. In addition we have implemented lightweight C# wrapper for convenient usage. Also we have run a lot of benchmarks with various GPU cards and on different platforms. On test cases provided by customer we received up to 150x speedup comparing to original single-threaded CPU implementation. However significant part of received speedup was due to poorly optimized original code which was completely rewritten and ported to CUDA whenever possible. 

While it is usually understood what FFT stage does in image reconstruction, GRAPPA stage is not so obvious. Due to parallel acquisition of different frames arises distortion of acquired data which is effectively eliminated. Figure 4 shows visual representation of images before and after reconstruction. 

Figure 4 The image before the reconstruction (left), image after reconstruction (right)

Additionally, you can find a case-study on ELEKS website or download it in PDF. Stay tuned!

/by Volodymyr Kozak, Principal Developer, ELEKS R&D

20x performance for the same HW price: GPGPU-based Machine Learning at Yandex

Russian search-engine Yandex has disclosed some details about their machine-learning framework, FML. The most interesting detail is that it runs on 80 TFLOPS cluster powered by GPGPU. This is quite unusual application for GPU, as ML algorithms are usually hard to be paralleled. However they have managed to adapt their decision tree algorithm for high-level of parallelism. As a result Yandex has achieved more than 20x speed-up for the same hardware price.
They are going to upgrade their cluster to 300 TFLOPS. Yandex expects its cluster to be in the list of top 100 most powerful supercomputers in the world after that upgrade.

Disruptor.NET

Disruptor is a framework for building high-performance applications that allows you to process large amount of messages concurrently without locks. It was initially invented by LMAX software development team members - Martin Thompson, Mike Barker and Dave Farley. It is said that with Disruptor you can achieve 100K transactions per second performance with less than 1 ms latency. I won't describe all Disruptor concepts as there are plenty of great articles about Disruptor in the Internet. My favorite ones are Martin Fowler's article, LMAX presentation at InfoQ and absolutely awesome series of blog-posts by Ruslan Cheremin (in Russian). Initially Disruptor is written with Java and released under Apache license. You can find source code along with some useful info and links on Github. There also is pretty good port to .NET CLR: disruptor-net (and today I'm going to tell more about it).

At the first glance Disruptor is yet another name for old good ring buffer. But if you look deeper there are several important features that make the difference. Well-written Disruptor-based applications can be incredibly fast. This speed is achieved by a great mix of high-level architecture choices and low-level optimizations based on deep knowledge of modern CPU and memory architecture. Almost all decisions by Disruptor authors were made in order to make applications cache-efficient, lock-free and GC friendly.

Today I'm going to show you how to write very simple Disruptor-based application. And just to make some difference with other introduction articles I'll use disruptor-net, so it will be Disruptor application written with C#. First of all, let's define the scope. Imagine there are lots of moving objects (e.g. cars). Each object has GPS installed on-board and sends its location via mobile network once per some period of time. We want to track the path of every single object. It is very important to log every single message to disk in order to have ability to restore all the movement history. We also want to aggregate the distance traveled by each object as well as the total distance traveled by all objects. Sounds like a very simple application, right?

Our moving object will have three properties: latitude, longitude and some unique ID. Something like that:

As I don't have a fleet of thousands moving objects with GPS sensors we will simulate their movement with very simple random simulator (see code on Github).

The simplest way to use disruptor in your .NET application is to install it via nuget package manager:

Let's create the Disruptor instance:

Here we pass factory method that preallocate the ring buffer with messages (do you remember that one of the main ideas of Disruptor is to reduce memory allocation?), set its size (1024, it must be a power of 2) and tell it to use default task scheduler. Now we need some handlers. Our handler graph will look like that:

Lets start with some simple persist handler that writes objects into text file, line by line:

As you can see it implements IEventHandler<T> interface. This interface has only one method void OnNext(T data, long sequence, bool endOfBatch). Disruptor call this method for each item in the ring buffer and pass the item as the first parameter. Second parameter is the sequence number for the item. Third parameter is the most interesting: as Disruptor allows batch processing it  indicates if current item is the last one in the batch. We use it in order reduce the number of flush calls - we call flush only when it is end of the batch. Leaping ahead, as ObjectPersistHandler is the slowest one in our pipeline, batch processing allowed us to achieve 2x performance boost in our benchmark. It looks like there is some cheating here, but in fact everything is ok - other handlers won't start processing before persist handler finish the batch, so there will be no situation when we process some message before it is stored to disk. Please note, we use simple string serialization here, in real-world application some more sophisticated serialization method will be used.
We have few other handlers: DistanceHandler, ConsoleLogHandler and AknowledgementHandler. First one is kind of business logic - it performs some distance calculation for objects:

I think this code is pretty straightforward, so no additional explanations are needed. The only thing you might notice here is that it doesn't use any locks or other synchronization primitives, as each handler is single threaded and doesn't share any mutable data with other handlers.

Aknowledgement handler actually tells simulator that we finished the processing and it is time to make next object move. In real world code there will be some kind of network message send where we will tell the client that its message was handled so it can send next one. After we move the object it is published into Disruptor using following code:

Please note, we do not store new object into the ring buffer, but copy its properties into the existing entry. It makes us sure that no data is shared between different threads and guarantee lower GC load.
The final thing we need to do is to configure disruptor to use our handlers. It could be done with fluent syntax like this:

Please note that 'Then()' method creates synchronization barrier, so we can guarantee that handlers will be executed strictly in required order. Console log handler and distance handler are included into one handler group, so they will be called in parallel threads.
The only thing we need is to start execution:

Let's measure the performance. On my work machine (Core i5@3.2 GHz, 8 GB RAM and regular spin HDD) I've got following results:

Of course in real-world app there will be some additional latency introduced by network communication, but still, 340K events per second is a pretty impressive throughput. Moreover, I think it could be much better in case we set the goal to achieve maximum possible performance here.

Disruptor gives you flexible way to create high-throughput and low-latency event processing applications. It is very simple but powerful tool with lots of great ideas inside and really good implementation. You can use .NET CLR port in case there are some reasons why you don't want to use original Java implementation.
As usual, you can find demo code on Github. Stay tuned!

NVIDIA Tesla K20 benchmark: facts, figures and some conclusions

Newest GPGPU flagman, Tesla K20 was announced by NVIDIA at Supercomputing conference in Salt Lake City yesterday (BTW, you can meet Roman Pavlyuk, ELEKS' CTO and Oleh Khoma, Head of HPC Unit there). Due to partnership with NVIDIA we got access to K20 couple of months ago and did lots of performance tests. Today we're going to tell you more about it's performance in comparison with several other NVIDIA accelerators that we have here at ELEKS. 

Test environment

We implemented set of synthetic micro-benchmarks that measure performance of following basic GPGPU operations:

• Host/Device kernel operations latency
• Reduction time (SUM)
• Dependent/Independent FLOPs
• Memory management
• Memory transfer speed
• Device memory access speed
• Pinned memory access speed

You can find more information and benchmark results below. Our set of tests is available on GitHub, so that you can run them on your hardware if you want. We ran these tests on seven different test configurations:

• GeForce GTX 580 (PCIe-2, OS Windows, physical box)
• GeForce GTX 680 (PCIe-2, OS Windows, physical box)
• GeForce GTX 680 (PCIe-3, OS Windows, physical box)
• Tesla K20Xm (PCIe-3, ECC ON, OS Linux, NVIDIA EAP server)
• Tesla K20Xm (PCIe-3, ECC OFF, OS Linux, NVIDIA EAP server)
• Tesla M2050 (PCIe-2, ECC ON, OS Linux, Amazon EC2)
• Tesla M2050 (PCIe-2, ECC ON, OS Linux, PEER1 HPC Cloud)

One of the goals was to determine the difference between K20 and older hardware configurations in terms of overall system performance. Another goal: to understand the difference between virtualized and non-virtualized environments. Here is what we got:

Host/Device kernel operations latency

One of the new features of K20 is Dynamic Parallelism that allows you to execute kernels from each other. We did a benchmark that measure latency of kernel schedule and execution with and without DP. Results without DP look like that:

Surprisingly, new Tesla is slower than old one and GTX 680, probably because of the driver which was in beta version at the moment we measured performance.  It is also obvious that AWS GPU instances are much slower than closer-to-hardware PEER1 ones, because of virtualization.

Then we tried to run similar benchmark with DP on:

Obviously we couldn't run these tests on older hardware because it doesn't support DP. Surprisingly, DP scheduling is slower than traditional one, but DP execution time is pretty much the same with ECC ON and traditional is faster with ECC OFF. We expected that DP latency would be less than traditional. It is hard to say what is the reason of such slowness. We suppose that it could be a driver, but it is just our assumption.

Reduction time (SUM)

Next thing we tried to measure was reduce execution time. Basically we calculated array sum. We did it with different arrays and grid sizes (Blocks x Threads x Array size):

Here we got expected results. New Tesla K20 is slower on small data sets, probably because of less clock frequency and not fully-fledged drivers. It becomes faster when we work with big arrays and use as many cores as possible. 

Regarding virtualization, we found that virtualized M2050 is comparable with non-virtualized one on small data sets, but much slower on large data sets. 

Dependent/Independent FLOPs

Peak theoretical performance is one of the most misunderstood properties of computing hardware. Some people says it means nothing, some says it is critical. The truth is always somewhere between these points. We tried to measure performance in FLOPs using several basic operations. We measured two types of operations, dependent and independent in order to determine if GPU does automatic parallelization of independent operations. Here's what we got:

Surprisingly, but we haven't got better results with independent operations. Probably we have some issue with our tests or misunderstood how does automatic parallelization work in GPU, but we couldn't implement the test where independent operations are automatically paralleled.
Regarding overall results, Teslas are much faster than GeForces when you work with double precision floating point numbers, which is expected: consumer accelerators are optimized for single precision because double is not required in computer games, primary software they were designed for. FLOPs are also highly dependent on clock speed and number of cores, so newer cards with more cores are usually faster, except of one case with GTX 580/680 and double precision: 580 is faster because of higher clock frequency.
Virtualization doesn't affect FLOPs performance at all.

Memory management

Another critical thing for HPC is basic memory management speed. As there are several memory models available in CUDA it is also critical to understand all the implications of using each of them. We wrote a test that allocate and release 16 b, 10 MB and 100 MB blocks of memory in different models. Please note: we got quite a different results in this benchmark, so it makes sense to show them on charts with logarithmic scale. Here they go:

Device memory is obviously the fastest option in case you allocate big chunk of memory. And GTX 680 with PCIe-3 is our champion in device memory management. Teslas are slower than GeForces in all the tests. Virtualization seriosly affects Host Write Combined memory management. PCIe-3 is better than PCIe-2 which is also obvious.

Memory transfer speed

Another important characteristics of an accelerator is speed of data transfer from one memory model to other. We measured it by copying 100 MB blocks of data between Host and GPU memory in both directions using regular, page locked and write combined memory access models. Here's what we got:

Obviously, PCIe-3 configurations are much faster than PCIe2. Kepler devices (GTX 680 and K20) are faster than other. If you use Page Locked and Write Combined models it makes your transfer speed faster. Virtualization slightly affects regular memory transfer speed, and doesn't affect others at all. We also tested internal memory transfer speed (please note, we haven't multiplied it by 2 as NVIDIA does usually in their tests):

Tesla K20s are faster than GeForce, but difference is not so big. M2050 are almost two times slower then their succesors.

Device memory access speed

We also measured device memory access speed for each configuration we have. Here they go:

Alligned memory access is way faster than non-aligned (almost 10 times difference). Newer accelerators are better than older. Double precicion read/write is faster than single for all the configurations. Virtualization doesn't affect memory access speed at all.

Pinned memory access speed

Last metric we measured was pinned memory access speed when device interacts with host memory. Unfortunately we weren't able to run these tests on GTX 680 with PCIe-3 due to issue with big memory blocks allocation in Windows. 

New Tesla is faster then old one. PCIe-3 is obviously faster. Aligned access is almost ten times faster and if you read double precision floats your memory access speed is two times bigger than if you work with single precision floats. Virtualized environment is slower than non-virtualized.

Conclusions

All-in-all new Tesla K20 performs slightly faster than their predecessors. There is no revolution. There is evolution - we got better performance, new tools that make programmer's life easier. There also are several things that are not mentioned in this benchmark, like better support of virtualization and as a result cloud-readyness of K20. Some results were surprising. We expect better results of K20 in several months when new, optimized version of drivers will be available (NVIDIA always has some issues with new drivers just after release, but usually fix them after several updates).

You can find spreadsheet with complete results at Google Docs. Benchmark sources are available at our GitHub.