While discussing the disappointing performance of my Futhark DCT on my “retro GPU”(Nvidia NVS 4200M) with Troels Henriksen, it came up that the Python backend has quite some calling overhead.
Futhark can compile high-level functional code to very fast OpenCL, but Futhark is meant to be embedded in larger programs. So it provides a host library in C and Python that set up the GPU, transfer the memory, and run the code. It turns out the the Python backend based on PyOpenCL is quite a bit slower at this than the C backend.
I wondered why the Python backend did not use the C one via FFI, and Troels mentioned that someone had done this for a specific program and saw modest performance gains. However, this does require a working compiler and OpenCL installation, rather than just a pip install PyOpenCL, so he argued that PyOpenCL is the easiest solution for the average data scientist.
I figured I might be able to write a generic wrapper for the generated C code by feeding the generated header directly to CFFI. That worked on the first try, so that was nice. The hard part was writing a generic, yet efficient and Pythonic wrapper around the CFFI module.
The first proof of concept required quite a few fragile hacks (pattern matching on function names and relying on the type and number of arguments to infer other things) But it worked! My DCT ran over twice as fast. Then, Troels, helpful as always, modified the generated code to reduce the number of required hacks. He then proceeded to port some of the demos and benchmarks, request some features, and contribute Python 2 support.
futhark-ffi now supports all Futhark types on both Python 2 and 3, resulting in speedups of anywhere between 20% and 100% compared to the PyOpenCL backend. Programs that make many short calls benefit a lot, while programs that call large, long-running code benefit very little. The OpenCL code that runs is the same, only the calling overhead is reduced.
One interesting change suggested by Troels is to not automatically convert Futhark to Python types. For my use case I just wanted to take a Numpy array, pass it to Futhark, and get a Numpy array back. But for a lot of other programs, the Futhark types are passed between functions unchanged, so not copying them between the GPU and CPU saves a lot of time. There is even a compatibility shim that lets you use futhark-ffi with existing PyOpenCL code by merely changing the imports. An example of this can be seen here
After installing Futhark, you can simply get my library with pip. (working OpenCL required)
pip install futhark-ffi
Usage is as follows. First generate a C library, and build a Python binding for it
futhark-opencl --library test.fut
build_futhark_ffi test
From there you can import both the CFFI-generated module and the library to run your Futhark code even faster!
If you search for Loefflers algorithm you get a few academic papers, and for Futhark you get the Germanic runes. This post is a SEO masterpiece.
Discrete Cosine Transform is a variation on the Discrete Fourier Transform. It is used in basically every lossy compression format ever. The reason DCT is preferred is that discrete transforms are cyclic. So the DFT has a jump at the edges of the data, where it wraps around. (this is why windowing is frequently used in DFT) This jump at the edges leads to a fat tail in the frequency spectrum, which does not compress well.
The DCT constructs an “even” signal (mirrored around the 0 axis), so the signal is continuous at the edges. This leads to much lower high frequency coefficients. Lossy compression basically works by quantizing/masking/thresholding those coefficients, which produces many zeros at high frequencies. Long runs of zeros compress really well, so that’s what happens in most compression algorithms.
I was playing a bit with compression, but found that scipy.fftpack.dct was not fast enough to my liking. Since I had recently discovered Futhark, which is an amazing ML-like functional programming language for GPU programming that compiles to OpenCL, I thought it’d be fun to implement the DCT in Futhark. Little did I know what I was getting myself into.
After some searching, I found that Loefflers algorithm is the way to go. It’s what everyone seems to be using, because for an 8-point DCT it obtains the theoretical lower bound of 11 multiplications. After chasing some references in more recent papers, I found the original: Practical fast 1-D DCT algorithms with 11 multiplications, and after days of struggling, I almost understood it.
I knew that a Fast Fourier Transform is based on taking the DFT equation, and splitting it up in odd and even parts. If you keep doing this recursively (called decimation in time/decimation in frequency), you end up with this “butterfly” structure, which are additions of two “branches” scaled by some factor. For the DCT there are also butterflies, but also rotation blocks.
It took a few mental leaps to understand that you can write the DCT or DFT in matrix form, express elementary row operations in matrix form, use those to factorize the DCT matrix, and derive an optimal implementation from this matrix factorization.
The Futhark side of things was really fun. If you know a bit of functional programming, it’s really not hard, and you don’t need to know anything about OpenCL or GPU’s. I hopped on Gitter, and Troels Henriksen was super helpful. I’d come up with a problem, and a few hours later I’d git pull and the compiler got better.
There were a few surprises though. Many array operations are basically free, by returning a view of the same array. But there is an underlying assumption that arrays are big, heap allocated, and parallelized relentlessly. Tuples, on the other hand, are assumed to be small, and register allocated. By rewriting my inner DCT structure from (tiny) arrays to using tuples, performance more than doubled.
At first I tied to optimize my code to use in-place updates, but this was actually significantly slower than out-of-place. By doing in-place updates, I force the compiler to do the operations completely sequentially, while normally it could do a lot in parallel. It turns out that moving data around is by far the slowest thing, and arithmetic is basically free. So the best way to write fast code is to move less data, not to worry about every addition.
Actually implementing the DCT was really hard though. As I mentioned, searching for it brought up only academic papers, barely any working code. I managed to find two eventually: dct_simd and mozjpeg. I actually ported the first one to Python to compare intermediate results with my own implementation.
From this code and staring at the paper, I learned a few things. First of all figure 1 is wrong. The rotate block should be sqrt(2)c6 instead of sqrt(2)c1. Another small detail is the dashed lines, meaning that some butterflies are upside down. Another one is the rotate block symbol. It says kcn, which are the k and n in the block equation, not the one in the DCT equation, which confused me a lot. So for sqrt(2)c6 you just substitute sqrt(2) and 6 in the rotation block. I noted down some more insights in response to a two year old question about the paper on the DSP StackExchange
Having implemented the forward DCT from the paper, I moved on to the inverse. All information the paper has about this is “just do everything backwards”. Thanks, paper. It turns out you use the same blocks, but in the reverse order, except… the rotate block n becomes -n. The inverse cosine transform has a negative angle, and this translates to cos(-x)=cos(x), sin(-x)=-sin(x).
type octet 't = [8]t
type f32octet = octet f32
let butterfly (a: f32) (b: f32) : (f32, f32) =
(a+b, a-b)
let mk_coef (k:f32) (n:i32) : f32 =
k*f32.cos ((r32 n)*f32.pi/16)
let coefr = map (mk_coef (f32.sqrt 2)) <| iota 8
let coef1 = map (mk_coef 1) <| iota 8
let rotate (sin_coef: f32) (cos_coef: f32) (a: f32) (b: f32): (f32, f32) =
(a*cos_coef + b*sin_coef,
b*cos_coef - a*sin_coef)
entry fdct8 (a: f32octet) : f32octet =
-- stage 1
let (st1_0, st1_7) = butterfly a[0] a[7]
let (st1_1, st1_6) = butterfly a[1] a[6]
let (st1_2, st1_5) = butterfly a[2] a[5]
let (st1_3, st1_4) = butterfly a[3] a[4]
-- even part, stage 2
let (st2_0, st2_3) = butterfly st1_0 st1_3
let (st2_1, st2_2) = butterfly st1_1 st1_2
-- stage 3
let (y0, y4) = butterfly st2_0 st2_1
let (y2, y6) = rotate coefr[2] coefr[6] st2_2 st2_3
-- odd part, stage 2
let (st2_4, st2_7) = rotate coef1[5] coef1[3] st1_4 st1_7
let (st2_5, st2_6) = rotate coef1[7] coef1[1] st1_5 st1_6
-- stage 3
let (st3_4, st3_6) = butterfly st2_4 st2_6
let (st3_7, st3_5) = butterfly st2_7 st2_5
-- stage 4
let (y1, y7) = butterfly st3_7 st3_4
let y3 = f32.sqrt(2)*st3_5
let y5 = f32.sqrt(2)*st3_6
in [y0, y4, y2, y6, y7, y3, y5, y1]
entry idct8 (a: f32octet) : f32octet =
-- odd part, stage 4
let (st4_7, st4_4) = butterfly a[7] a[4]
let st4_5 = f32.sqrt(2)*a[5]
let st4_6 = f32.sqrt(2)*a[6]
-- stage 3
let (st3_4, st3_6) = butterfly st4_4 st4_6
let (st3_7, st3_5) = butterfly st4_7 st4_5
-- stage 2
let (st2_4, st2_7) = rotate (-coef1[5]) coef1[3] st3_4 st3_7
let (st2_5, st2_6) = rotate (-coef1[7]) coef1[1] st3_5 st3_6
-- even part, stage 3
let (st3_0, st3_1) = butterfly a[0] a[1]
let (st3_2, st3_3) = rotate (-coefr[2]) coefr[6] a[2] a[3]
-- stage 2
let (st2_0, st2_3) = butterfly st3_0 st3_3
let (st2_1, st2_2) = butterfly st3_1 st3_2
-- stage 1
let (st1_0, st1_7) = butterfly st2_0 st2_7
let (st1_1, st1_6) = butterfly st2_1 st2_6
let (st1_2, st1_5) = butterfly st2_2 st2_5
let (st1_3, st1_4) = butterfly st2_3 st2_4
in [st1_0/8, st1_1/8, st1_2/8, st1_3/8, st1_4/8, st1_5/8, st1_6/8, st1_7/8]
I try to keep my Linux system as clean as possible, but for almost every project I inadvertently have to install some packages with apt-get, build some from source with make install, or even install some proprietary program that in turn requires more packages. So over time my system acquires more and more junk.
To solve this problem, I made “virualenv for system packages”, a little script that makes a chroot with an unholy combination of OverlayFS and bind mounts. From inside the chroot, it looks and behaves exactly like your /, with no overhead or isolation whatsoever. The only difference is that writes to all system directories go to an overlay directory.
$ sudo bash activate.sh
(env)$ sudo apt-get install everything
(env)$ sudo make install
(env)$ echo "Hello world" > hello.txt
(env)$ ^D
$ everything
Command not found
$ cat hello.txt
Hello world
In the Python world this is pretty much the standard. You make a virtualenv, pip install all the things you need, and delete then environment after your project is done. I want this kind of behavior for all of my software.
Of course you can run everything in a VM, a Docker image, or a chroot. But these typically provide isolation that I don’t want or need. They also have a lot of overhead in RAM, disk space, and most importantly, effort on my part. However, the git-like overlay filesystems used by Docker gave me an idea.
What if I made a chroot, but instead of putting an entire Debian installation inside it, make an overlay on my own system. The only trouble is that special directories like /proc, /sys, and /dev should work as usual, and preferably my /home folder should also persist outside the chroot.
The solution turns out to be relatively simple: My /home and all the special directories are mounted with mount --bind. All system directories like /lib and /bin are mounted with mount -t overlay.
It works great for from-source installs and proprietary software. It also works for apt-get, but I’m sure weird things will happen once you upgrade your system. Once it gets wonky, just nuke the env and start over.
