Comment by rpep | Hacker News Reader

Comment by rpep | Hacker News Reader

I've been through the journey with compiling extensions to Python in lower languages (primarily C and C++) and it's really worth evaluating the benefit. In general, the control that it gives you over memory layout and chaining of operations make it very difficult to say from toy examples what benefit you'll get.

If you test a single array operation with NumPy then it will compare favourably to a low level implementation, but where it generally will lose out is when you’re chaining together multiple operations which a compiler will vectorise more efficiently. For e.g. doing things like add two arrays then multiply them by another array and subtract a constant. If I put that in a C or Rust function I’d expect that to be auto-vectorized using fewer operations than what NumPy will do, and you of course also remove the overhead of the foreign function interface and control returning to the interpreter. On top of that, it’s usually trivial to drop in thread level parallelism in a lower level language, and NumPy doesn’t operate with core level parallelism out of the box for most operations, you have to use something like numexpr to achieve it.

I've personally used OpenMP + SIMD intrinsics or auto-vectorisation to great effect in scientific software and got performance >100x that of NumPy, but it's certainly not a free lunch. The question you have to ask yourself is normally “is it worth it”, and that’s something that can only be answered by profiling, looking at the overhead of maintaining the architecture, and understanding whether the increased friction with debugging and building is worth the hassle. Additionally, if you're doing linear algebra then generally if you're dropping to a lower level language you still will want to be using the same libraries that NumPy etc. uses unless you're exploiting some property of the matrix structure that can make operations more efficient in memory or compute. I worked on one codebase for e.g. where there was a block matrix structure and a hand written implementation of a solver was written because much more efficient operations could be performed with knowledge of that structure.

I've written a couple of C-extensions for a Python service that needed to run in real-time. We started out 100x slower than real-time and managed to make it slightly faster than real-time in a hectic couple of weeks.

The main bottleneck that C-extensions solved for us then was marshaling. Instead of doing something like Numpy --> custom algorithm --> Numpy --> custom algorithm in Python, we moved it all to C. This way, we decreased the number of times we had to marshal data & we reduced the amount of data that needed to be marshaled.

IIRC we gained 10-100x improved performance in all scenarios where we used C-extensions. Before doing so we tried to use Numba, as it seemed to offer a "free lunch", but we never got it to work.

The final cherry on top that let us surpass real-time speed was replacing a graph algorithm with another that had a better time complexity for larger graphs. Before doing that switch, we'd still be faster than real-time for _most_ inputs but would lag behind severely whenever we had a dense input graph come in.

It's not really meaningful to say "100x slower than real-time" or "slightly faster than real-time", "real-time" isn't some threshold of sub-second precision. It usually refers to relative (though maybe absolute wall-clock) time of execution for instructions that are typically controlling hardware that needs precise timings (ie. the time between turning this thing on and turning it back off needs to be exactly 150us).

It is meaningful, but maybe I didn't give enough context.

The service ingests data from a timespan T spanning M minutes. Before the next timespan T+1 is ingested, the service has to process the data from timespan T in <M minutes otherwise we fall behind. What I mean with "faster than real-time" is that it has to process the data faster than it is ingested.

NumPy is insanely fast for most use cases, I've (grudgingly) come to the conclusion that if NumPy doesn't hack it then I should re-think the problem or my way of solving it rather than to try to optimize that particular bit of code if it isn't meant for something that is going to be run in production on a large number of machines. Likely there are better uses of my time. It's interesting how what is nominally a scripting language can perform so well for compute intensive tasks. Of course, technically you are just stringing together highly optimized low level functions together using the Python language but the advantages of doing it that way have more often than not surprised me.

This

Op is better off recompiling numpy to target his specific hardware thank trying to speed it up using rust. This would be to try different underlying math implementations that are linked with numpy ( blas, mkl, eigen). Each of these haha several internal simulation frameworks to optimize not just the actual instructions but memory layout for various math kernels.

How do you handle it in numpy if you want something like ((x - y)/2 + z) * w all elementwise over the arrays? Naively that's 3 intermediate unused arrays being created.

There's plenty of CPU and GPU numpy accelerators available.

* Numba: https://numba.pydata.org/

* JAX: https://jax.readthedocs.io/en/latest/notebooks/quickstart.ht...

numba - instead of writing in Rust, you write it in numba, which is also almost like another language. Not bad, but needs to be taken into account and is not pure numpy.

Inplace operators?

You use Julia and let loop fusion handle it. You get 1 allocation if you want to store the result in a new array, or zero allocations if you already have the array allocated.

https://julialang.org/blog/2017/01/moredots/

Isn't that what Theano/Aesara is for? Still based on a NumPy interface AIUI, but automagically compiling from NumPy-based to efficient code on the CPU or GPU.

Is there an optimal solution to that for all sizes of the arrays?

E.g. I’d expect for small/medium then re-writing as a fold over the input arrays would be fastest because you’d only traverse once & everything would fit in cache.

However if the 3 arrays combined size is larger than L1 cache, I’d be strongly tempted to bet on the naïve 3 operation approach to be faster. You’d save so much time on cache line flush / reload activities.

But then if 2 arrays are larger than L3 i’d expect going back to the fold / single traversal to win again because the iterating 2 arrays at a time behaviour is no longer any different to iterating 3 at a time.

Untested hypothesis.

You could use NumExpr: https://github.com/pydata/numexpr#what-is-numexpr

Implementing the hot path in e.g. C, compiling it into a lib an calling it using e.g. ctypes is very little effort and can yield dramatic runtime improvements / reduced bill. It's been an excellent use of my time, at least.

I'd recommend using Numba first, though.

e.g. comes from "exempli gratia" - roughly "for the sake of an example", or shortly "for example". In other words the "gratia" part already includes the "for", you don't have to write it again. E.g. this is how I would use it in a sentence.

Including another for would spell out to "for for example"

> but where it generally will lose out is when you’re chaining together multiple operations which a compiler will vectorise more efficiently

I think the bigger issue is remaining cache-friendly in the face of multiple operations. If you do the operations one at a time, you’re doing multiple passes over the data.

In any case, the approach that Polars[0] took seems like a pretty good solution, if difficult to implement well.

Polars is a dataframe library, and has a “lazy” API. This builds up a graph of operations, runs it through a query optimizer, and then executes them all at once. This allows for parallelization, optimizing expressions, and fewer passes over the data. The downside is that it’s nontrivial to implement a query optimizer.

[0]: https://pola.rs/

Reply on news.ycombinator.com