Getting a 10-1000x or more improvement on existing code is very common without putting in a ton of effort if the code was not already heavily optimized. These are listed roughly in order of importance, but performance is often such a non-consideration from most developers that a little effort goes a long way.

1. Most importantly, is the algorithm a good choice? Can we eliminate some work entirely? (this is what algo interviews are testing for)

2. Can we eliminate round trips to the kernel and similar heavy operations? The most common huge gain here is replacing tons of malloc calls with a custom allocator.

3. Can we vectorize? Explicit vector intrinsics like in the blog post are great, but you can often get the same machine code by reorganizing your data into arrays / struct of arrays rather than arrays of structs.

4. Can we optimize for cache efficiency? If you already reorganized for vectors this might already be handled, but this can get more complicated with parallel code if you can't isolate data to one thread (false sharing, etc.)

5. Can we do anything else that's hardware specific? This can be anything from using intrinsics to hand-coding assembly.

There's always going to be some fixed overhead and latency (and there's a great article about the impact of latency on performance called "It's the latency, stupid" that's worth a read: http://www.stuartcheshire.org/rants/latency.html) but sending far more data than is needed over a network connection will sooner or later kill performance.

Overall though, I agree with your considerations, and in roughish terms the order of them.

From your description though, it doesn't sound like something I'd classify as a network issue. That's just classic orm nonsense. I guess I don't know what you mean by "distributed query", but it sounds terrible.

The most classic network performance issue is forgetting to disable nagle's algorithm.

The most classic sql performance issue is not using an index.

When considering performance issues, in my experience it's a mistake not to explicitly consider the network.

I think OP's #1 agrees that there's nothing "technically" wrong with such a query (or an algo). It just generated work you didn't have to do. Work takes time. So you used time you didn't have to use. I also think this is the number 1 way of improving performance in general (computer, life).

A perfectly valid query of fetching 1M rows turned into 99.xxx% unnecessary work when you only needed a handful of rows. The query wasn't slow, it was just generating more work than you actually needed. The network also wasn't slow, it simply had to transfer (even at peak theoretical efficiency) a lot of data you never used.

You then used an equally valid query that wasn't even necessarily fast, it just generated much less work. This query (quote from #1) "eliminate[d] some work entirely", the work of carrying over unnecessary data.

SIMD intrinsics aren't needed for micro-kernel vectorization, as a decent C compiler will fully vectorize and unroll it. BLIS' pure C micro-kernel gets >80% of the performance of the hand-optimized implementation on Haswell with appropriate block sizes. The difference is likely to be due to prefecth, but I don't properly understand it.

- The matrices weren't that big. Numpy has cffi overhead.

- The perf difference was much more noticeable with _peak_ throughput rather than _mean_ throughput, which matters for almost no applications (a few, admittedly, but even where "peak" is close to the right measure you usually want something like the mean of the top-k results or the proportion with under some latency, ...).

- The benchmarking code they displayed runs through Python's allocator for numpy and is suggestive of not going through any allocator for the C implementation. Everything might be fine, but that'd be the first place I checked for microbenchmarking errors or discrepancies (most numpy routines allow in-place operations; given that that's known to be a bottleneck in some applications of numpy, I'd be tempted to explicitly examine benchmarks for in-place versions of both).

- Numpy has some bounds checking and error handling code which runs regardless of the underlying implementation. That's part of why it's so bleedingly slow for small matrices compared to even vanilla Python lists (they tested bigger matrices too, so this isn't the only effect, but I'll mention it anyway). It's hard to make something faster when you add a few thousand cycles of pure overhead.

- This was a very principled approach to saturating the relevant caches. It's "obvious" in some sense, but clear engineering improvements are worth highlighting in discussions like this, in the sense that OpenBLAS, even with many man-hours, likely hasn't thought of everything.

And so on. Anyone can rattle off differences. A proper explanation requires an actual deep-dive into both chunks of code.

Could be MKL (i believe the conda version comes with it) but it could also be an ancient version of OpenBLAS you already had installed. So yeah, being faster than np.matmul probably just means your NumPy is not installed optimally.

By that reasoning, wouldn't it make more sense to wrap their C code and maybe even make it operate on numpy's array representation, so it can be called from Python?

Suppose that this blog post were part of a series that questions the axiom (largely bolstered by academic marketing) that one needs Python to do array computing. Then it is valid to compare C directly to NumPy.

It isn't even far fetched. The quality of understanding something after having implemented it in C is far greater than the understanding gained by rearranging PyTorch or NumPY snippets.

That said, the Python overhead should not be very high if M=1000, N=1000, K=1000 was used. The article is a bit silent on the array sizes, this is somewhere from the middle of the article.

It would be a huge detriment and time sink for these data scientist to take the time to learn to write an equivalent C program if their ultimate goal is to do data science.

The author also says "(...) implementation follows the BLIS design", but then proceeds to compare *only* with OpenBLAS. I'd love to see a more thorough analysis, and using C directly would make it easier to compare multiple BLAS libs.

Look, I'm trying to be charitable to the authors, hard as that might be.

I actually do have that question niggling in the back of my mind when I use something like NumPy. I don't necessarily care exactly _where_ the overhead comes from, I might just be interested whether it's close to ideal or not.

What they are doing here is patting themselves on the back after handicapping the competition. Not to mention that they have given themselves the chance to cherry pick the very best hyperparameters for this particular comparison while BLAS is limited to using heuristics to guess which of their kernels will suit this particular combination of hardware and parameters.

The authors need to be called out for this contrived comparison.

Very moot nitpick though, given that this is for only one column of the matrix, the following loops of maskload/maskstore will take significantly more time (esp. store, which is still slow on Zen 4[1] despite the AVX-512 instruction (whose only difference is taking the mask in a mask register) being 6x faster), and clang autovectorizes the shifting anyways (maybe like 2-3x slower than my suggestions).

[1]: https://uops.info/table.html?search=vmaskmovps&cb_lat=on&cb_...

Regarding "creating a constant global array and loading from it" - if I recall correctly, I've tested this approach and it was a bit slower than bit mask shifting. But let me re-test this to be 100% sure.

"Comparing a constant vector {0, 1, 2, 3, 4, ...} with broadcasted m and m-8" - good idea, I will try it!

Note you can keep int8_t elements in that array, and sign extend bytes into int32_t while loading. The _mm_loadu_si64 / _mm256_cvtepi8_epi32 combo should compile into a single vpmovsxbd instruction with a memory operand. This way the entire constant array fits in a single cache line, as long as it’s aligned properly with alignas(32)

This is good fit for the OP’s use case because they need two masks, the second vpmovsxbd instruction will be a guaranteed L1D cache hit.

> Important! Please don’t expect peak performance without fine-tuning the hyperparameters, such as the number of threads, kernel and block sizes, unless you are running it on a Ryzen 7700(X). More on this in the tutorial.

I think I'll need a TL;DR on what to change all these values to.

I have a Ryzen 7950X and as a first test I tried to only change NTHREADS to 32 in benchmark.c, but matmul.c performed worse than NumPy on my machine.

So I took a look at the other values present in the benchmark.c, but MC and NC are already calculated via the amount of threads (so these are probably already 'fine-tuned'?), and I couldn't really understand how KC = 1000 fits for the 7700(X) (the author's CPU) and how I'd need to adjust it for the 7950X (with the informations from the article).

Also, if you were designing for smaller cases, say MNK=16 or 32, how would you approach it differently? I'm implementing neural ODEs and this is one point I've been considering.

If you want to convert to Python lists its is going to take time. Not sure about Python arrays.

cursed code. I checked and every single invocation is just `int`, so why do this? you can just write a function:

    func min(x, y int) int {
       if x < y { return x }
       return y
    }

Great job! Self publishing things like this were a hallmark of the early internet I for one sorely miss.