Over the past week, I’ve been working on a small image convolution kernel in C++. It applies a box blur to a grayscale image. I wanted a small problem I could use to learn performance engineering: read the assembly, predict what a change would do, and figure out why my prediction was wrong. I also wanted to get more familiar with ARM SIMD so I could apply some of this to mohg, my synthesizer.

I kept the scope to direct convolution on a single grayscale image, running on one M4 performance core. The code is on GitHub.

Estimating the theoretical limit

Before optimizing anything, I wanted a rough number to compare against. I used one output pixel as the unit of work and assumed I’d vectorize across 4 output pixels at a time.

For a \(k \times k\) kernel, each output pixel requires \(k^2\) multiply-adds. Assuming 4 floats per NEON instruction and 4 such instructions per cycle, I get an arithmetic throughput bound of:

\[\frac{k^2}{16} \text{ cycles/pixel}\]

For \(k=3\), that’s about 0.56 cycles/pixel. I also need to load the inputs. In this model, each 4-wide float load handles one tap for 4 output pixels, so I need \(\frac{k^2}{4}\) loads per pixel. At an assumed 3 loads/cycle, the load throughput bound is:

\[\frac{k^2}{12} \text{ cycles/pixel}\]

For \(k=3\), that’s 0.75 cycles/pixel, or about 0.17 ns/pixel at 4.4 GHz. The load bound is higher than the arithmetic bound for all the kernel sizes I’m testing.

I’m leaving out loop overhead, address calculations, dependencies, and other instructions, so this is an optimistic estimate. I wanted a starting point, but I don’t know how close I can get to it.

Results

I benchmarked on a \(741 \times 488\) grayscale image. These are kernel-only timings in ns/pixel.

variant change \(k=3\) \(k=9\)
v0 baseline, runtime \(k\), nested vectors 2.153 17.452
v1 templated \(k\) 1.939 6.342
v2 flattened image 1.654 6.015
v3 converted image to floats beforehand 1.578 5.495
v4 inverted loops 0.823 7.523
vImage Apple’s vImageConvolve_PlanarF 0.477 3.236
theoretical limit   0.17 1.53

For \(k=3\), my fastest version is about 5x slower than my estimated limit and 1.7x slower than Apple’s vImage. Even vImage is about 2–3x slower than the estimate. I don’t know how much of that gap comes from costs I’ve left out of the model.

Most of my predictions were wrong

I tried to predict the effect of each change before implementing it. The numbers below include earlier measurements from those experiments.

Vectorization

The first thing I tried was making \(k\) a template parameter so the compiler would know the kernel size at compile time. This helped, but the compiler still wasn’t vectorizing the loop. I checked the vectorization remarks and found:

cannot prove it is safe to reorder floating-point operations

I tried -ffast-math, expecting a 3–4x speedup because one 4-wide vector instruction could do the work of four scalar instructions. I got 1.13x for \(k=3\) and 2.8x for \(k=9\).

My prediction only accounted for the arithmetic. I hadn’t accounted for the loads and other instructions that would remain. Looking at the assembly, I also found that the compiler was vectorizing across kernel taps for one output pixel. I had expected more benefit from the four lanes without checking how the compiler would use them.

Flattening the image

I started with a vector<vector<uint8_t>>. I thought replacing it with a flat array wouldn’t make much difference because the extra pointer loads seemed small compared to the convolution itself. But flattening the image gave me a 1.29x speedup for \(k=3\).

I went back to the assembly and counted 22 loads per output pixel in v1 and 17 in v2. Assuming 3 loads/cycle, I get load throughput bounds of about 7.3 and 5.6 cycles/pixel. That’s a difference of about 1.7 cycles/pixel, close to the measured difference of 1.86. I checked \(k=5\) and \(k=7\), and those estimates were close too.

That made me suspect load throughput was part of the problem. I hadn’t predicted the speedup beforehand, but I could now account for much of the difference with the load counts. That doesn’t prove loads were the bottleneck; I’d still need to check what else changed.

Too many instructions

For v3, I converted the input to floats before running the convolution. I then checked the hardware counters and measured an IPC of about 5.4, with around 36 instructions per output pixel. My model accounted for only 5–6 instructions per pixel, so I wanted to understand where the rest came from.

I took the IPC as a reason to investigate instruction count, though it wasn’t enough to tell me which execution units were limiting performance. Looking at the loop, I was still calculating one output pixel at a time. I wanted to try vectorizing across outputs.

Inverting the loops

I moved the output-pixel loop to the inside. For each kernel tap, I would now loop across the output image and add that tap’s contribution. The compiler could then vectorize across output pixels. For \(k=3\), this gave me almost a 2x speedup over v3, but performance got worse as \(k\) increased. At \(k=9\), v4 was slower than v3.

My guess is that I’m moving too much data. Each kernel tap makes another pass over the input image and reads and writes the output buffer. Each float image is around 1.4 MB, much larger than L1 cache, and \(k=9\) means 81 passes. I haven’t checked the counters yet, so I don’t know whether that explains the slowdown.

What I’m timing

For \(k=3\), v4 was 2.62x faster than the baseline when I timed only the convolution. Including allocations and the conversion from floats back to bytes brought the speedup down to 1.72x. I’d been concentrating on the hot loop and hadn’t given much thought to the work around it.

Applying this to mohg

Individual voices in my synthesizer are independent, much like output pixels. I want to try vectorizing across voices and compare that with vectorizing the processing within one voice. I’ll also need to look at how I lay out and access the voice data, especially if the extra memory traffic explains what happened with v4.

What’s next?

I want to check the counters for v4 at \(k=3\) and \(k=9\) to see whether memory traffic explains the regression. I’m also curious what vImage does differently. It’s about 1.7x faster than my best version for \(k=3\), and I have a larger gap to close at \(k=9\).

I’m still getting comfortable reading assembly. For now, I want to keep making small changes, write down what I expect, and work through the cases where the results don’t match.