1: https://github.com/adewes/python-qubit-setup/blob/master/lib...
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.
* Numba: https://numba.pydata.org/
* JAX: https://jax.readthedocs.io/en/latest/notebooks/quickstart.ht...
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.
I'd recommend using Numba first, though.
Including another for would spell out to "for for example"
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.
Would you mind sharing what the stumbling blocks were for your usage of numba? We've been using it extensively for numerical computations with great success, provided any branching logic can be factored out.
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.
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/
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.