ChrisRackauckas

SciPy+Numba odeint vs Julia DifferentialEquations.jl vs NumbaLSODA Summary

All are solved at reltol=1e-3, abstol=1e-6 using the fastest ODE solver of the respective package for the given problem.

Absolute Performance Numbers:

• SciPy LSODA through odeint takes ~489μs
• SciPy LSODA through odeint with Numba takes ~257μs
• NumbaLSODA takes ~25μs
• DifferentialEquations.jl Rosenbrock23 takes ~9.2μs

ChrisRackauckas

using OrdinaryDiffEq, RecursiveArrayTools, LinearAlgebra, Test, SparseArrays, SparseDiffTools, Sundials

# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0

ChrisRackauckas

# Atom:
Version: 1.44.0
Dev Mode: false
Official Release: true
{
  "http_parser": "2.8.0",
  "node": "10.11.0",
  "v8": "6.9.427.31-electron.0",
  "uv": "1.23.0",
  "zlib": "1.2.11",

ChrisRackauckas

multithreadedfjac = (var"##MTIIPVar#581", var"##MTKArg#579")->begin
        @inbounds begin
                let (u₁ˏ₁ˏ₁, u₂ˏ₁ˏ₁, u₃ˏ₁ˏ₁, u₄ˏ₁ˏ₁, u₅ˏ₁ˏ₁, u₆ˏ₁ˏ₁, u₇ˏ₁ˏ₁, u₈ˏ₁ˏ₁, u₉ˏ₁ˏ₁, u₁₀ˏ₁ˏ₁, u₁₁ˏ₁ˏ₁, u₁₂ˏ₁ˏ₁, u₁₃ˏ₁ˏ₁, u₁₄ˏ₁ˏ₁, u₁₅ˏ₁ˏ₁, u₁₆ˏ₁ˏ₁, u₁ˏ₂ˏ₁, u₂ˏ₂ˏ₁, u₃ˏ₂ˏ₁, u₄ˏ₂ˏ₁, u₅ˏ₂ˏ₁, u₆ˏ₂ˏ₁, u₇ˏ₂ˏ₁, u₈ˏ₂ˏ₁, u₉ˏ₂ˏ₁, u₁₀ˏ₂ˏ₁, u₁₁ˏ₂ˏ₁, u₁₂ˏ₂ˏ₁, u₁₃ˏ₂ˏ₁, u₁₄ˏ₂ˏ₁, u₁₅ˏ₂ˏ₁, u₁₆ˏ₂ˏ₁, u₁ˏ₃ˏ₁, u₂ˏ₃ˏ₁, u₃ˏ₃ˏ₁, u₄ˏ₃ˏ₁, u₅ˏ₃ˏ₁, u₆ˏ₃ˏ₁, u₇ˏ₃ˏ₁, u₈ˏ₃ˏ₁, u₉ˏ₃ˏ₁, u₁₀ˏ₃ˏ₁, u₁₁ˏ₃ˏ₁, u₁₂ˏ₃ˏ₁, u₁₃ˏ₃ˏ₁, u₁₄ˏ₃ˏ₁, u₁₅ˏ₃ˏ₁, u₁₆ˏ₃ˏ₁, u₁ˏ₄ˏ₁, u₂ˏ₄ˏ₁, u₃ˏ₄ˏ₁, u₄ˏ₄ˏ₁, u₅ˏ₄ˏ₁, u₆ˏ₄ˏ₁, u₇ˏ₄ˏ₁, u₈ˏ₄ˏ₁, u₉ˏ₄ˏ₁, u₁₀ˏ₄ˏ₁, u₁₁ˏ₄ˏ₁, u₁₂ˏ₄ˏ₁, u₁₃ˏ₄ˏ₁, u₁₄ˏ₄ˏ₁, u₁₅ˏ₄ˏ₁, u₁₆ˏ₄ˏ₁, u₁ˏ₅ˏ₁, u₂ˏ₅ˏ₁, u₃ˏ₅ˏ₁, u₄ˏ₅ˏ₁, u₅ˏ₅ˏ₁, u₆ˏ₅ˏ₁, u₇ˏ₅ˏ₁, u₈ˏ₅ˏ₁, u₉ˏ₅ˏ₁, u₁₀ˏ₅ˏ₁, u₁₁ˏ₅ˏ₁, u₁₂ˏ₅ˏ₁, u₁₃ˏ₅ˏ₁, u₁₄ˏ₅ˏ₁, u₁₅ˏ₅ˏ₁, u₁₆ˏ₅ˏ₁, u₁ˏ₆ˏ₁, u₂ˏ₆ˏ₁, u₃ˏ₆ˏ₁, u₄ˏ₆ˏ₁, u₅ˏ₆ˏ₁, u₆ˏ₆ˏ₁, u₇ˏ₆ˏ₁, u₈ˏ₆ˏ₁, u₉ˏ₆ˏ₁, u₁₀ˏ₆ˏ₁, u₁₁ˏ₆ˏ₁, u₁₂ˏ₆ˏ₁, u₁₃ˏ₆ˏ₁, u₁₄ˏ₆ˏ₁, u₁₅ˏ₆ˏ₁, u₁₆ˏ₆ˏ₁, u₁ˏ₇ˏ₁, u₂ˏ₇ˏ₁, u₃ˏ₇ˏ₁, u₄ˏ₇ˏ₁, u₅ˏ₇ˏ₁, u₆ˏ₇ˏ₁, u₇ˏ₇ˏ₁, u₈ˏ₇ˏ₁, u₉ˏ₇ˏ₁, u₁₀ˏ₇ˏ₁, u₁₁ˏ₇ˏ₁, u₁₂ˏ₇ˏ₁, u₁

ChrisRackauckas

fjac = (var"##MTIIPVar#585", var"##MTKArg#583")->begin
        @inbounds begin
                let (u₁ˏ₁ˏ₁, u₂ˏ₁ˏ₁, u₃ˏ₁ˏ₁, u₄ˏ₁ˏ₁, u₅ˏ₁ˏ₁, u₆ˏ₁ˏ₁, u₇ˏ₁ˏ₁, u₈ˏ₁ˏ₁, u₉ˏ₁ˏ₁, u₁₀ˏ₁ˏ₁, u₁₁ˏ₁ˏ₁, u₁₂ˏ₁ˏ₁, u₁₃ˏ₁ˏ₁, u₁₄ˏ₁ˏ₁, u₁₅ˏ₁ˏ₁, u₁₆ˏ₁ˏ₁, u₁ˏ₂ˏ₁, u₂ˏ₂ˏ₁, u₃ˏ₂ˏ₁, u₄ˏ₂ˏ₁, u₅ˏ₂ˏ₁, u₆ˏ₂ˏ₁, u₇ˏ₂ˏ₁, u₈ˏ₂ˏ₁, u₉ˏ₂ˏ₁, u₁₀ˏ₂ˏ₁, u₁₁ˏ₂ˏ₁, u₁₂ˏ₂ˏ₁, u₁₃ˏ₂ˏ₁, u₁₄ˏ₂ˏ₁, u₁₅ˏ₂ˏ₁, u₁₆ˏ₂ˏ₁, u₁ˏ₃ˏ₁, u₂ˏ₃ˏ₁, u₃ˏ₃ˏ₁, u₄ˏ₃ˏ₁, u₅ˏ₃ˏ₁, u₆ˏ₃ˏ₁, u₇ˏ₃ˏ₁, u₈ˏ₃ˏ₁, u₉ˏ₃ˏ₁, u₁₀ˏ₃ˏ₁, u₁₁ˏ₃ˏ₁, u₁₂ˏ₃ˏ₁, u₁₃ˏ₃ˏ₁, u₁₄ˏ₃ˏ₁, u₁₅ˏ₃ˏ₁, u₁₆ˏ₃ˏ₁, u₁ˏ₄ˏ₁, u₂ˏ₄ˏ₁, u₃ˏ₄ˏ₁, u₄ˏ₄ˏ₁, u₅ˏ₄ˏ₁, u₆ˏ₄ˏ₁, u₇ˏ₄ˏ₁, u₈ˏ₄ˏ₁, u₉ˏ₄ˏ₁, u₁₀ˏ₄ˏ₁, u₁₁ˏ₄ˏ₁, u₁₂ˏ₄ˏ₁, u₁₃ˏ₄ˏ₁, u₁₄ˏ₄ˏ₁, u₁₅ˏ₄ˏ₁, u₁₆ˏ₄ˏ₁, u₁ˏ₅ˏ₁, u₂ˏ₅ˏ₁, u₃ˏ₅ˏ₁, u₄ˏ₅ˏ₁, u₅ˏ₅ˏ₁, u₆ˏ₅ˏ₁, u₇ˏ₅ˏ₁, u₈ˏ₅ˏ₁, u₉ˏ₅ˏ₁, u₁₀ˏ₅ˏ₁, u₁₁ˏ₅ˏ₁, u₁₂ˏ₅ˏ₁, u₁₃ˏ₅ˏ₁, u₁₄ˏ₅ˏ₁, u₁₅ˏ₅ˏ₁, u₁₆ˏ₅ˏ₁, u₁ˏ₆ˏ₁, u₂ˏ₆ˏ₁, u₃ˏ₆ˏ₁, u₄ˏ₆ˏ₁, u₅ˏ₆ˏ₁, u₆ˏ₆ˏ₁, u₇ˏ₆ˏ₁, u₈ˏ₆ˏ₁, u₉ˏ₆ˏ₁, u₁₀ˏ₆ˏ₁, u₁₁ˏ₆ˏ₁, u₁₂ˏ₆ˏ₁, u₁₃ˏ₆ˏ₁, u₁₄ˏ₆ˏ₁, u₁₅ˏ₆ˏ₁, u₁₆ˏ₆ˏ₁, u₁ˏ₇ˏ₁, u₂ˏ₇ˏ₁, u₃ˏ₇ˏ₁, u₄ˏ₇ˏ₁, u₅ˏ₇ˏ₁, u₆ˏ₇ˏ₁, u₇ˏ₇ˏ₁, u₈ˏ₇ˏ₁, u₉ˏ₇ˏ₁, u₁₀ˏ₇ˏ₁, u₁₁ˏ₇ˏ₁, u₁₂ˏ₇ˏ₁, u₁₃ˏ₇ˏ₁, u₁₄ˏ₇ˏ

ChrisRackauckas

:((var"##MTIIPVar#414", var"##MTKArg#412")->begin
          @inbounds begin
                  @sync begin
                          let (u₁, u₂, u₃, u₄, u₅, u₆, u₇, u₈, u₉, u₁₀) = (var"##MTKArg#412"[1], var"##MTKArg#412"[2], var"##MTKArg#412"[3], var"##MTKArg#412"[4], var"##MTKArg#412"[5], var"##MTKArg#412"[6], var"##MTKArg#412"[7], var"##MTKArg#412"[8], var"##MTKArg#412"[9], var"##MTKArg#412"[10])
                              begin
                                  Threads.@spawn begin
                                          (var"##MTIIPVar#414").nzval[1] = -2
                                          (var"##MTIIPVar#414").nzval[2] = 1
                                          (var"##MTIIPVar#414").nzval[3] = 1
                                          (var"##MTIIPVar#414").nzval[4] = -2

ChrisRackauckas

torchsde vs DifferentialEquations.jl / DiffEqFlux.jl (Julia)

This example is a 4-dimensional geometric brownian motion. The code for the torchsde version is pulled directly from the torchsde README so that it would be a fair comparison against the author's own code. The only change to that example is the addition of a dt choice so that the simulation method and time step matches between the two different programs.

The SDE is solved 100 times. The summary of the results is as follows:

ChrisRackauckas

torchdiffeq vs Julia DiffEqFlux Neural ODE Training Benchmark

The spiral neural ODE was used as the training benchmark for both torchdiffeq (Python) and DiffEqFlux (Julia) which utilized the same architecture and 500 steps of ADAM. Both achived similar objective values at the end. Results:

• DiffEqFlux defaults: 7.4 seconds
• DiffEqFlux optimized: 2.7 seconds
• torchdiffeq: 288.965871299999 seconds

ChrisRackauckas

DifferentialEquations.jl Runge-Kutta Integrators vs Taylor Integration

This is against a Taylor integration method that utilizes taylor-mode automatic differentiation in a specifically optimized way for ODE integration via its @taylorize macro (for more details, see the documentation). This demonstration is done on the classic Pleiades n-body problem, a non-stiff ODE since the Taylor integration methods are only for non-stiff equations.

ChrisRackauckas

Conservative Package Compiler Results

1. Using a ton of standard pervasive packages only helps 10%
2. Getting DiffEqBase helps 50%
3. Getting OrdinaryDiffEq without extra precompiles gets to 1.3 seconds, vs 0.6 seconds on the second run. That's pretty good!