I recently had to compute many inner products with a given matrix $\Ab$ for many different vectors $\xb_i$, or $\xb_i^T \Ab \xb_i$. Each vector $\xb_i$ represents a shoe from Zappos and there are 50k vectors $\xb_i \in \R^{1000}$. This is computation took place behind a user-facing web interface and during testing had a delay of 5 minutes. This is clearly unacceptable; how can we make it faster?1

I spent a couple hours trying to get the best possible performance from my functions… and through this, I found a speed optimization2 that put most of the computation on NumPy’s shoulders. After I made this change, the naïve for-loop and NumPy were about a factor of 2 apart, not enough to write a blog post about.

Use of a NVIDIA GPU significantly outperformed NumPy. Given that most of the optimization seemed to be focused on a single matrix multiplication, let’s focus on speed in matrix multiplication.

We know that matrix multiplication has computational complexity of something like $\bigO{n^{2.8074}}$3, but very likely greater than $\bigO{n^{2.375477}}$4 when multiplying two $n\times n$ matrices. We can’t get around this without diving into theory, but we can change the constant that dictates exactly how fast these algorithms run.

The tools I’ll test are

• the default NumPy install, with no MKL (even though it’s now provided by default with Anaconda)
• Intel MKL, a tool that provides acceleration for BLAS/LAPACK
• the GPU. To do this, I’ll need an Amazon AWS machine and the NVIDIA CUDA Toolkit. An easy interface is available through cudamat but scikit-cuda and Accelerate also have nice interfaces and provide more access.

I had planned to test other tools but these tests didn’t pan out for reasons in the footnotes. My test script can be summarized as follows:

import numpy as np
import cudamat as cm

n, p = int(2e3), int(40e3)
A = np.random.randn(n, p)
B = np.random.randn(p, n)
%timeit A @ B

cm.cublas_init()
cm.CUDAMatrix.init_random()
A_cm = cm.empty((n, p)).fill_with_randn()
B_cm = cm.empty((p, n)).fill_with_randn()
%timeit A_cm.dot(B_cm)
cm.cublas_shutdown()


When doing this, I generate the following graph:

Environment NumPy + no MKL NumPy + MKL cudamat
Time (seconds) 7.18 4.057 0.2898

Under the default Anaconda environment (i.e., with MKL), we see that our script runs 80% slower without MKL and has a 14x speedup under cudamat!

This simple test shows that using the GPU is powerful, but cudamat is limited in that it only provides basic mathematical capability for the GPU (the dot product is as far as it goes).

However, the libraries Accelerate and scikit-cuda use the GPU and both provide more complex mathematical functions, including fft, svd and eig.

Accelerate and scikit-learn are both fairly similar. In choosing whether to use Accelerate or scikit-learn, there are two obvious tradeoffs:

• scikit-cuda has access to linear algebra functions (e.g., eig) and Anaconda does not. However, access to these higher level mathematical functions comes through CULA, another framework that requires a license (free academic licenses are available).
• Anaconda can accept raw ndarrays and scikit-cuda needs to have gpuarrays passed in (meaning more setup/cleanup).

edit: Other GPU libraries numba.cuda.jit, numba.hsa.jit also exist.

Whichever is chosen, large speed enhancements exist. I have timed a common function (fft) over different values of n; there is some overhead to moving to the GPU and I wanted to see where that is. I provide a summary of my testing script in the appendix.

CULA has benchmarks for a few higher-level mathematical functions (source: the CULA Dense homepage):

## Appendix

### Other GPU libraries

Anaconda has published a good overview titled “Getting started with GPU computing”. I think I would start with Numba: it has debugging and supports some notion of kernels. [updated 2017-11]

…and of course I didn’t optimize any loop-based functions. To do optimize loop speed, I would look at numba first and then possibly Cython.

### FFT timing script summary

In this script, I show preparing for the FFT and preparing for linear algebra functions (e.g., cilinalg.init()). I found that it’s useful to look at the scikit-cuda demos.

import numpy as np
from accelerate.cuda.blas import Blas
import accelerate.cuda.fft as acc_fft
import pycuda.autoinit
import pycuda.gpuarray as gpuarray
import skcuda.fft as cu_fft
import skcuda.linalg as culinalg
import skcuda.misc as cumisc

# for scikit-learn
culinalg.init()

# for accelerate when calling wrapped BLAS functions (e.g., blas.dot)
blas = Blas()

def fft_accelerate(x, y):
f = acc_fft.FFTPlan(shape=x.shape, itype=x.dtype, otype=y.dtype)
f.forward(x, out=y)  # note: we're passing np.ndarrays
return y

def fft_scikit(x, y):
plan_forward = cu_fft.Plan(x.shape, np.float32, np.complex64)
cu_fft.fft(x, y, plan_forward)
return y.get()

n = int(40e4)
x = np.random.randn(n).astype('float32')
y = np.zeros(n).astype('complex64')  # needed because fft has complex output
%timeit fft_accelerate(x, y)

x = gpuarray.to_gpu(x)
y = gpuarray.empty(n//2 + 1, np.complex64)
%timeit fft_scikit(x, y)

1. Note: I link to the libraries I discovered in “Other GPU libraries”. Update 2017-11: Look here for libraries.

2. which was calculating $\Ab \Xb^T$ outside the loop

3. using the Strassen algorithm