Thanks for providing code. I’ve written code like that before and the first thing that jumps out at me is that you have for loops in Python code which isn’t slow but will never be as fast as loops in a static compiled language. Sometimes for loops are the most readable way of expressing a calculation but as a numerical computation person my instinct is to look for opportunities to vectorize by rewriting loops into matrix notation on paper and then expressing them as array calculations.

Otherwise you’re right — whenever you have “for” loops you’re always going to end up ahead if you rewrite those in C, Rust, Julia or Fortran. And that’s exactly what authors of high performance numerical codes do.

Some of that Python is strange for numerical Python practices, yes, but I was aiming for close floating-point output equivalence. Rust ndarray provides the standard set of matrix ops, but I keep getting different results (\propto 1e-5 on a stable basis, which isn't quite a proof of incorrectness, but....) when I use too much numpy stuff like linalg.norm, etc. (ndarray slices in Rust are annoying to use because of lifetimes and such.)

I did my dissertation on symplectic exponential Runge-Kutta schemes with stable numerics for the big matrix exponentials needed (which have this specific structure that allow for theorems to be proven). But I didn't have time to write any code at all. I wonder if open source ODE solvers are getting good high-order symplectic methods by now...

>I wonder if open source ODE solvers are getting good high-order symplectic methods by now...

Up to 10th order https://diffeq.sciml.ai/stable/solvers/dynamical_solve/#Symp... . Also Magnus methods https://diffeq.sciml.ai/stable/solvers/nonautonomous_linear_... and exponential integrators https://diffeq.sciml.ai/stable/solvers/split_ode_solve/#Ordi... . And the expmv implementations are highly optimized as well: https://github.com/SciML/ExponentialUtilities.jl . A lot of this specialized matrix exponential stuff is rather fun, for example see https://github.com/SciML/ExponentialUtilities.jl/pull/64

I don't know about other languages, but in Julia I've heard people often say that loops end up faster than the equivalent vectorized code. So while this is true for Python/Matlab I don't think it is good universal advice.

That said, matrix notation can sometimes be the more readable way of expressing a calculation.

I might have used an early beta of Julia circa 2018 or something, but the chorus that it performs like $static_fastlang doesn't match the experience I had.

Do you have code to share and look at? All we can do is point to real-world code and benchmarks. For this specific case, see LoopVectorization.jl results (https://juliasimd.github.io/LoopVectorization.jl/latest/exam...), and the corresponding effects on stiff ODE solver benchmarks against C and Fortran packages (https://benchmarks.sciml.ai/html/StiffODE/Hires.html).

I'd don't know what you mean by "performs like" but it's definitely significantly closer in runtime speed after the first run to compiled languages than interpreted (in particular Python/Matlab).

Anyways, my point was more in response to this from the person I responded to: "my instinct is to look for opportunities to vectorize by rewriting loops into matrix notation on paper and then expressing them as array calculations". In some languages (Julia in particular) that is slower than the equivalent loop based code.