# For which statistical methods are GPUs faster than CPUs?

I have just installed a Nvidia GT660 graphic card on my desktop and, after some struggle, I manage to interface it with R.

I have been playing with several R packages that use GPUs, especially gputools, and I was comparing the time taken by my GPU and CPU to perform some basic operations:

• inverting matrices (CPU faster)
• qr decomposition (CPU faster)
• big correlation matrices (CPU faster)
• matrix multiplication (GPU much faster!)

Notice that I have experimented mainly with gputools so maybe other packages perform better.

In broad terms my question is: what some routine statistical operations that might be worth executing on a GPU rather than a CPU?

• Anything involving lots of matrix multiplication? :) GPUs are quite popular in the neural nets community. – larsmans Feb 24 '13 at 23:26
• you need to provide the size of the matrices involved. For example, last i checked (admittedly 2 years ago) inversion and decomposition were only faster on GPU starting from large matrices (2^9 times 2^9 and upwards) – user189035 Feb 24 '13 at 23:29
• I used matrices of around $10^3 \times 10^3$ for inversion, qr and matrix multiplication, while for correlations I have used around 10^4 observations of vectors of size 100. For matrix inversion the GPU was much slower, while for qr decomposition it was slower but comparable to the CPU. – Jugurtha Feb 24 '13 at 23:38
• this is a very good question but i think you'll get better answers by having it migrated to stackoverflow (through i think similar questions have been asked there before) – user189035 Feb 25 '13 at 0:06
• GPU's advantage of regular CPU's is the fact that they can be "massively" parallel, not that they are faster per core. As such, for jobs that require a lot of "housekeeping" like Cholesky factorization etc. you need to use block algorithms and so forth to get significant speed-up; this is not trivial and I assume it will take a while before GPU's take over such operations. What is definitely going the GPU way is MCMC-ing (and Random Number generation). Sampling from a posterior has "parallelization" written all over it... And sparse matrices computations; they are already "blocked" anyway... – usεr11852 says Reinstate Monic Feb 25 '13 at 7:13

In broad terms, algorithms that run faster on the GPU are ones where you are doing the same type of instruction on many different data points.

An easy example to illustrate this is with matrix multiplication.

Suppose we are doing the matrix computation

$A \times B = C$

A simple CPU algorithm might look something like

//starting with C = 0

for (int i = 0; i < C_Width; i++)
{
for (int j = 0; j < C_Height; j++)
{
for (int k = 0; k < A_Width; k++)
{
for (int l = 0; l < B_Height; l++)
{
C[j, i] += A[j, k] * B[l, i];
}
}
}
}


The key thing to see here is that there are a lot of nested for loops and each step must be executed one after the other.

See a diagram of this

Notice that the calculation of each element of C does not depend on any of the other elements. So it does not matter what order the calculations are done in.

So on the GPU, these operations can be done concurrently.

A GPU kernel for calulating a matrix multiplication would look something like

__kernel void Multiply
(
__global float * A,
__global float * B,
__global float * C
)
{
const int x = get_global_id(0);
const int y = get_global_id(1);
for (int k = 0; k < A_Width; k++)
{
for (int l = 0; l < B_Height; l++)
{
C[x, y] += A[x, k] * B[l, y];
}
}
}


This kernel only has the two inner for loops. A program sending this job to the GPU will tell the GPU to execute this kernel for each data point in C. The GPU will do each of these instructions concurrently on many threads. Just like the old saying "Cheaper by the dozen" GPUs are designed to be faster doing the same thing lots of times.

There are however some algorithms which will slow the GPU down. Some are not well suited for the GPU.

If for example, there were data dependencies, ie: imagine the computation of each element of C depended on the previous elements. The programmer would have to put a barrier in the kernel to wait for each previous computation to finish. This would be a major slow down.

Also, algorithms which have a lot of branching logic ie:

__kernel Foo()
{
if (somecondition)
{
do something
}
else
{
do something completely different
}
}


tend to run slower on the GPU because the GPU is no longer doing the same thing in each thread.

This is a simplified explanation because there are many other factors to consider. For example, sending data between the CPU and GPU is also time consuming. Sometimes it is worth doing a computation on the GPU even when its faster on the CPU, just to avoid the extra send time (And vice versa).

Also many modern CPUs support concurrency now as well with hyperthreaded multicore processors.

GPU's also seem to be not so good for recursion, see here which probably explains some of the problems with the QR algorithm. I believe that one has some recursive data dependencies.

• It's officially SX-naughty to comment on an answer just to say that it's a terrific answer, but I don't give a rat's perinæum about negs: this is a delightful and informative answer. One of SX's great injustices is the lack of kudos to people who give exquisitely-informative answers on 'old' (in internet time) questions. (Plus, I'm giving a thumbs-up to an 'old' (in internet time) answer: I know, right? META). – GT. May 16 '15 at 7:26
• An important consideration is whether there is actually a library to do the computation: e.g. to my knowledge, there are no sparse x dense GPU implementations of matrix multiplication, certainly not through R packages. If you are prepared to work with writing GPU C code, then good luck. – Jack Wasey Jun 4 '18 at 9:00

GPUs are sensitive beasts. Although Nvidia's beefiest card can theoretically execute any of the operations you listed 100x faster than the fastest CPU, about a million things can get in the way of that speedup. Every part of the relevant algorithm, and of the program which runs it, has to be extensively tweaked and optimized in order to get anywhere near that theoretical maximum speedup. R is generally not known to be a particularly fast language, and so it doesn't surprise me that its default GPU implementation is not that great, at least in terms of raw performance. However, the R GPU functions may have optimization settings that you can tweak in order to regain some of that missing performance.

If you're looking into GPUs because you've found that some calculation that you need to run is going to take weeks/months to finish, it may be worth your while to migrate from R to a more performance-friendly language. Python isn't too much harder to work with than R. The NumPy and SciPy packages have most of the same stat functions as R, and PyCuda can be used to implement your own GPU based functions in a fairly straightforward way.

If you really want to increase the speed at which your functions run on GPUs, I would consider implementing your own functions in a combination of C++ and CUDA. The CUBLAS library can be used to handle all of the linear algebra-related heavy lifting. However, keep in mind that it can take quite a while to write such code (especially if it's your first time doing so), and so this approach should be reserved only for those computations that take an extremely long time to run (months) and/or that you're going to be repeating hundreds of times.

For all of the applications you mentioned, GPUs should be more capable (from a hardware perspective) than CPUs for sufficiently large matrices. I don't know anything about R's implementation, but I've used cuBLAS and Magma with great success for inversions around $n = 2^{10}$ and multiplication/correlation for rectangular matrices with $n,m \approx 2^{10}, k \approx 2^{14}$. It is an especially big surprise to me that large correlation matrices would be faster on the CPU using R.

More broadly, I suspect most statistical operations that spend most of their time in dense linear algebra (BLAS, Lapack functionality) can be efficiently implemented on the GPU.

Multiple Imputation methods for Missing Data? Like those in Alice-II (R).

I think those tend to be often embarrassingly parallel and hence suitable to a GPU architecture. Never tried it myself though.