4 ms·

I recently learned just enough Rust to use its ndarray crate and rewrite bottleneck bits of code. I used to think the whole "write inner loops in $fastlang, use

I recently learned just enough Rust to use its ndarray crate and rewrite bottleneck bits of code. I used to think the whole "write inner loops in $fastlang, use $scriptlang as glue" was something of a lame "cope" (and I say this as someone who doesn't know $fastlang), but I'm seeing upwards of 100X time performance improvements for code that can't be hammered into numpy idioms. I can provide code examples. Edit: code example https://gist.github.com/asemic-horizon/2830ed3637cfd278e7937f56450fd271 https://gist.github.com/asemic-horizon/2830ed3637cfd278e7937...

wenc 5y ago
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.
- fasterbynight 5y ago
  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.
  - prionassembly 5y ago
    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.
    - ChrisRackauckas 5y ago
      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/examples/matrix_multiplication/ 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 https://benchmarks.sciml.ai/html/StiffODE/Hires.html).
    - fasterbynight 5y ago
      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.
  - deleted 5y ago
    [deleted]
- prionassembly 5y ago
  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...
  - ChrisRackauckas 5y ago
    >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/#Symplectic-Integrators https://diffeq.sciml.ai/stable/solvers/dynamical_solve/#Symp... . Also Magnus methods https://diffeq.sciml.ai/stable/solvers/nonautonomous_linear_ode/#State-Independent-Solvers https://diffeq.sciml.ai/stable/solvers/nonautonomous_linear_... and exponential integrators https://diffeq.sciml.ai/stable/solvers/split_ode_solve/#OrdinaryDiffEq.jl-2 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 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 https://github.com/SciML/ExponentialUtilities.jl/pull/64

wenc 5y ago

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.

fasterbynight 5y ago
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.
- prionassembly 5y ago
  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.
  - ChrisRackauckas 5y ago
    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/examples/matrix_multiplication/ 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 https://benchmarks.sciml.ai/html/StiffODE/Hires.html).
  - fasterbynight 5y ago
    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.
- deleted 5y ago
  [deleted]
prionassembly 5y ago
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...
- ChrisRackauckas 5y ago
  >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/#Symplectic-Integrators https://diffeq.sciml.ai/stable/solvers/dynamical_solve/#Symp... . Also Magnus methods https://diffeq.sciml.ai/stable/solvers/nonautonomous_linear_ode/#State-Independent-Solvers https://diffeq.sciml.ai/stable/solvers/nonautonomous_linear_... and exponential integrators https://diffeq.sciml.ai/stable/solvers/split_ode_solve/#OrdinaryDiffEq.jl-2 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 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 https://github.com/SciML/ExponentialUtilities.jl/pull/64