4 ms·

This looks like a nice write-up and implementation. I'm left wondering what is the "trick"? How does it manage to beat OpenBLAS, which is assembly+C optimized o

by ks2048 2y ago

This looks like a nice write-up and implementation. I'm left wondering what is the "trick"? How does it manage to beat OpenBLAS, which is assembly+C optimized over decades for this exact problem? It goes into detail about caching, etc - is BLAS is not taking advantage of these things, or is this more tuned to this specific processor, etc?

sbstp 2y ago
Maybe -march=native gives it an edge as it compiles for this exact CPU model whereas numpy is compiled for a more generic (older) x86-64. -march=native would probably get v4 on a Ryzen CPU where numpy is probably targeting v1 or v2. https://en.wikipedia.org/wiki/X86-64#Microarchitecture_levels https://en.wikipedia.org/wiki/X86-64#Microarchitecture_level...
- stingraycharles 2y ago
  Doesn’t numpy have runtime SIMD dispatching and whatnot based on CPU flags? E.g. https://github.com/numpy/numpy/blob/main/numpy/_core/src/common/npy_cpu_dispatch.h https://github.com/numpy/numpy/blob/main/numpy/_core/src/com...
  - KeplerBoy 2y ago
    np.matmul just uses whatever blas library your NumPy distribution was configured for/shipped with. 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.
hansvm 2y ago
- OpenBLAS isn't _that_ optimized for any specific modern architecture. - 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.
- robxorb 2y ago
  To your third point - it looks as if the lines of mean values were averaged, this posts code would still be a clear winner.
ipsum2 2y ago
Comparison with numpy 2.0 should be better for numpy because it integrates Google highway for better simd across different microarchitectures.
- janwas 2y ago
  Highway TL here. The integration is a work in progress; what has been integrated so far is VQSort :)
  - ipsum2 2y ago
    Reading the numpy 2.0 change log I assumed Highway was deeply integrated, thanks for the correction.
teleforce 2y ago
To beat OpenBLAS is neither suprising nor it's not unheard of, for example D language library Mir for linear algebra did just that many years ago [1]. For C++ and C implementation please check these metaprogramming approaches [2], [3]. What's really surprising is that many modern languages are still depending on OpenBLAS for examples Matlab, Julia, Mojo, etc. But to be fair they probably have their own reasons, right? [1] Numeric age for D: Mir GLAS is faster than OpenBLAS and Eigen (2016): http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/glas-gemm-benchmark.html http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/... [2] Vastly outperforming LAPACK with C++ metaprogramming (2018): https://wordsandbuttons.online/vastly_outperforming_lapack_with_cpp_metaprogramming.html https://wordsandbuttons.online/vastly_outperforming_lapack_w... [3] Outperforming LAPACK with C metaprogramming (2018): https://wordsandbuttons.online/outperforming_lapack_with_c_metaprogramming.html https://wordsandbuttons.online/outperforming_lapack_with_c_m...

by ks2048 2y ago

sbstp 2y ago

Maybe -march=native gives it an edge as it compiles for this exact CPU model whereas numpy is compiled for a more generic (older) x86-64. -march=native would probably get v4 on a Ryzen CPU where numpy is probably targeting v1 or v2. https://en.wikipedia.org/wiki/X86-64#Microarchitecture_levels https://en.wikipedia.org/wiki/X86-64#Microarchitecture_level...

stingraycharles 2y ago
Doesn’t numpy have runtime SIMD dispatching and whatnot based on CPU flags? E.g. https://github.com/numpy/numpy/blob/main/numpy/_core/src/common/npy_cpu_dispatch.h https://github.com/numpy/numpy/blob/main/numpy/_core/src/com...
- KeplerBoy 2y ago
  np.matmul just uses whatever blas library your NumPy distribution was configured for/shipped with. 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.

hansvm 2y ago

- OpenBLAS isn't _that_ optimized for any specific modern architecture. - 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.

robxorb 2y ago
To your third point - it looks as if the lines of mean values were averaged, this posts code would still be a clear winner.

ipsum2 2y ago

Comparison with numpy 2.0 should be better for numpy because it integrates Google highway for better simd across different microarchitectures.

janwas 2y ago
Highway TL here. The integration is a work in progress; what has been integrated so far is VQSort :)
- ipsum2 2y ago
  Reading the numpy 2.0 change log I assumed Highway was deeply integrated, thanks for the correction.

teleforce 2y ago

To beat OpenBLAS is neither suprising nor it's not unheard of, for example D language library Mir for linear algebra did just that many years ago [1]. For C++ and C implementation please check these metaprogramming approaches [2], [3]. What's really surprising is that many modern languages are still depending on OpenBLAS for examples Matlab, Julia, Mojo, etc. But to be fair they probably have their own reasons, right? [1] Numeric age for D: Mir GLAS is faster than OpenBLAS and Eigen (2016): http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/glas-gemm-benchmark.html http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/... [2] Vastly outperforming LAPACK with C++ metaprogramming (2018): https://wordsandbuttons.online/vastly_outperforming_lapack_with_cpp_metaprogramming.html https://wordsandbuttons.online/vastly_outperforming_lapack_w... [3] Outperforming LAPACK with C metaprogramming (2018): https://wordsandbuttons.online/outperforming_lapack_with_c_metaprogramming.html https://wordsandbuttons.online/outperforming_lapack_with_c_m...