Optimizing DiffEq Code

Chris Rackauckas

In this notebook we will walk through some of the main tools for optimizing your code in order to efficiently solve DifferentialEquations.jl. User-side optimizations are important because, for sufficiently difficult problems, most of the time will be spent inside of your f function, the function you are trying to solve. "Efficient" integrators are those that reduce the required number of f calls to hit the error tolerance. The main ideas for optimizing your DiffEq code, or any Julia function, are the following:

  • Make it non-allocating

  • Use StaticArrays for small arrays

  • Use broadcast fusion

  • Make it type-stable

  • Reduce redundant calculations

  • Make use of BLAS calls

  • Optimize algorithm choice

We'll discuss these strategies in the context of small and large systems. Let's start with small systems.

Optimizing Small Systems (<100 DEs)

Let's take the classic Lorenz system from before. Let's start by naively writing the system in its out-of-place form:

function lorenz(u,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 [dx,dy,dz]
end
lorenz (generic function with 1 method)

Here, lorenz returns an object, [dx,dy,dz], which is created within the body of lorenz.

This is a common code pattern from high-level languages like MATLAB, SciPy, or R's deSolve. However, the issue with this form is that it allocates a vector, [dx,dy,dz], at each step. Let's benchmark the solution process with this choice of function:

using DifferentialEquations, BenchmarkTools
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  10.83 MiB
  allocs estimate:  101434
  --------------
  minimum time:     3.308 ms (0.00% GC)
  median time:      3.375 ms (0.00% GC)
  mean time:        3.887 ms (13.02% GC)
  maximum time:     9.092 ms (33.26% GC)
  --------------
  samples:          1286
  evals/sample:     1

The BenchmarkTools package's @benchmark runs the code multiple times to get an accurate measurement. The minimum time is the time it takes when your OS and other background processes aren't getting in the way. Notice that in this case it takes about 5ms to solve and allocates around 11.11 MiB. However, if we were to use this inside of a real user code we'd see a lot of time spent doing garbage collection (GC) to clean up all of the arrays we made. Even if we turn off saving we have these allocations.

@benchmark solve(prob,Tsit5(),save_everystep=false)
BenchmarkTools.Trial: 
  memory estimate:  9.47 MiB
  allocs estimate:  88652
  --------------
  minimum time:     2.909 ms (0.00% GC)
  median time:      2.917 ms (0.00% GC)
  mean time:        3.259 ms (9.71% GC)
  maximum time:     6.673 ms (30.76% GC)
  --------------
  samples:          1533
  evals/sample:     1

The problem of course is that arrays are created every time our derivative function is called. This function is called multiple times per step and is thus the main source of memory usage. To fix this, we can use the in-place form to ***make our code non-allocating***:

function lorenz!(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
lorenz! (generic function with 1 method)

Here, instead of creating an array each time, we utilized the cache array du. When the inplace form is used, DifferentialEquations.jl takes a different internal route that minimizes the internal allocations as well. When we benchmark this function, we will see quite a difference.

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  1.38 MiB
  allocs estimate:  12996
  --------------
  minimum time:     778.454 μs (0.00% GC)
  median time:      785.712 μs (0.00% GC)
  mean time:        851.452 μs (6.89% GC)
  maximum time:     3.895 ms (72.64% GC)
  --------------
  samples:          5859
  evals/sample:     1
@benchmark solve(prob,Tsit5(),save_everystep=false)
BenchmarkTools.Trial: 
  memory estimate:  5.22 KiB
  allocs estimate:  54
  --------------
  minimum time:     346.300 μs (0.00% GC)
  median time:      349.081 μs (0.00% GC)
  mean time:        349.778 μs (0.11% GC)
  maximum time:     4.335 ms (91.49% GC)
  --------------
  samples:          10000
  evals/sample:     1

There is a 4x time difference just from that change! Notice there are still some allocations and this is due to the construction of the integration cache. But this doesn't scale with the problem size:

tspan = (0.0,500.0) # 5x longer than before
prob = ODEProblem(lorenz!,u0,tspan)
@benchmark solve(prob,Tsit5(),save_everystep=false)
BenchmarkTools.Trial: 
  memory estimate:  5.22 KiB
  allocs estimate:  54
  --------------
  minimum time:     1.738 ms (0.00% GC)
  median time:      1.745 ms (0.00% GC)
  mean time:        1.746 ms (0.00% GC)
  maximum time:     1.854 ms (0.00% GC)
  --------------
  samples:          2862
  evals/sample:     1

since that's all just setup allocations.

But if the system is small we can optimize even more.

Allocations are only expensive if they are "heap allocations". For a more in-depth definition of heap allocations, there are a lot of sources online. But a good working definition is that heap allocations are variable-sized slabs of memory which have to be pointed to, and this pointer indirection costs time. Additionally, the heap has to be managed and the garbage controllers has to actively keep track of what's on the heap.

However, there's an alternative to heap allocations, known as stack allocations. The stack is statically-sized (known at compile time) and thus its accesses are quick. Additionally, the exact block of memory is known in advance by the compiler, and thus re-using the memory is cheap. This means that allocating on the stack has essentially no cost!

Arrays have to be heap allocated because their size (and thus the amount of memory they take up) is determined at runtime. But there are structures in Julia which are stack-allocated. structs for example are stack-allocated "value-type"s. Tuples are a stack-allocated collection. The most useful data structure for DiffEq though is the StaticArray from the package StaticArrays.jl. These arrays have their length determined at compile-time. They are created using macros attached to normal array expressions, for example:

using StaticArrays
A = @SVector [2.0,3.0,5.0]
3-element StaticArrays.SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
 2.0
 3.0
 5.0

Notice that the 3 after SVector gives the size of the SVector. It cannot be changed. Additionally, SVectors are immutable, so we have to create a new SVector to change values. But remember, we don't have to worry about allocations because this data structure is stack-allocated. SArrays have a lot of extra optimizations as well: they have fast matrix multiplication, fast QR factorizations, etc. which directly make use of the information about the size of the array. Thus, when possible they should be used.

Unfortunately static arrays can only be used for sufficiently small arrays. After a certain size, they are forced to heap allocate after some instructions and their compile time balloons. Thus static arrays shouldn't be used if your system has more than 100 variables. Additionally, only the native Julia algorithms can fully utilize static arrays.

Let's ***optimize lorenz using static arrays***. Note that in this case, we want to use the out-of-place allocating form, but this time we want to output a static array:

function lorenz_static(u,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 @SVector [dx,dy,dz]
end
lorenz_static (generic function with 1 method)

To make the solver internally use static arrays, we simply give it a static array as the initial condition:

u0 = @SVector [1.0,0.0,0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz_static,u0,tspan)
@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  466.28 KiB
  allocs estimate:  2595
  --------------
  minimum time:     336.309 μs (0.00% GC)
  median time:      340.596 μs (0.00% GC)
  mean time:        353.354 μs (3.13% GC)
  maximum time:     2.111 ms (78.97% GC)
  --------------
  samples:          10000
  evals/sample:     1
@benchmark solve(prob,Tsit5(),save_everystep=false)
BenchmarkTools.Trial: 
  memory estimate:  3.33 KiB
  allocs estimate:  28
  --------------
  minimum time:     237.341 μs (0.00% GC)
  median time:      240.594 μs (0.00% GC)
  mean time:        241.256 μs (0.00% GC)
  maximum time:     277.143 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

And that's pretty much all there is to it. With static arrays you don't have to worry about allocating, so use operations like * and don't worry about fusing operations (discussed in the next section). Do "the vectorized code" of R/MATLAB/Python and your code in this case will be fast, or directly use the numbers/values.

Exercise 1

Implement the out-of-place array, in-place array, and out-of-place static array forms for the Henon-Heiles System and time the results.

Optimizing Large Systems

Interlude: Managing Allocations with Broadcast Fusion

When your system is sufficiently large, or you have to make use of a non-native Julia algorithm, you have to make use of Arrays. In order to use arrays in the most efficient manner, you need to be careful about temporary allocations. Vectorized calculations naturally have plenty of temporary array allocations. This is because a vectorized calculation outputs a vector. Thus:

A = rand(1000,1000); B = rand(1000,1000); C = rand(1000,1000)
test(A,B,C) = A + B + C
@benchmark test(A,B,C)
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  2
  --------------
  minimum time:     1.102 ms (0.00% GC)
  median time:      1.124 ms (0.00% GC)
  mean time:        1.214 ms (7.41% GC)
  maximum time:     2.059 ms (37.65% GC)
  --------------
  samples:          4100
  evals/sample:     1

That expression A + B + C creates 2 arrays. It first creates one for the output of A + B, then uses that result array to + C to get the final result. 2 arrays! We don't want that! The first thing to do to fix this is to use broadcast fusion. Broadcast fusion puts expressions together. For example, instead of doing the + operations separately, if we were to add them all at the same time, then we would only have a single array that's created. For example:

test2(A,B,C) = map((a,b,c)->a+b+c,A,B,C)
@benchmark test2(A,B,C)
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  8
  --------------
  minimum time:     2.652 ms (0.00% GC)
  median time:      2.665 ms (0.00% GC)
  mean time:        2.755 ms (3.20% GC)
  maximum time:     3.586 ms (25.01% GC)
  --------------
  samples:          1813
  evals/sample:     1

Puts the whole expression into a single function call, and thus only one array is required to store output. This is the same as writing the loop:

function test3(A,B,C)
    D = similar(A)
    @inbounds for i in eachindex(A)
        D[i] = A[i] + B[i] + C[i]
    end
    D
end
@benchmark test3(A,B,C)
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  2
  --------------
  minimum time:     1.135 ms (0.00% GC)
  median time:      1.190 ms (0.00% GC)
  mean time:        1.470 ms (8.36% GC)
  maximum time:     3.864 ms (23.62% GC)
  --------------
  samples:          3386
  evals/sample:     1

However, Julia's broadcast is syntactic sugar for this. If multiple expressions have a ., then it will put those vectorized operations together. Thus:

test4(A,B,C) = A .+ B .+ C
@benchmark test4(A,B,C)
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  2
  --------------
  minimum time:     1.098 ms (0.00% GC)
  median time:      1.124 ms (0.00% GC)
  mean time:        1.215 ms (7.55% GC)
  maximum time:     2.096 ms (43.88% GC)
  --------------
  samples:          4097
  evals/sample:     1

is a version with only 1 array created (the output). Note that .s can be used with function calls as well:

sin.(A) .+ sin.(B)
1000×1000 Array{Float64,2}:
 0.840378  0.485274  0.924572  0.757179  …  1.39955   0.773387  0.657392
 1.41323   0.793729  0.831187  0.950668     0.654319  1.46885   1.52804
 1.20182   0.653114  0.482528  0.95313      0.913951  0.222537  0.610174
 1.49242   0.868178  1.01561   1.13976      1.24413   0.399607  0.527826
 1.07206   0.813371  0.816984  0.811787     0.258578  1.17443   1.42978
 0.842218  1.32533   1.08801   0.582951  …  0.926491  1.38897   0.849803
 0.989953  1.06931   1.33685   0.549195     1.49118   0.709916  0.906874
 0.469665  0.601914  0.377434  0.467431     0.862483  1.00223   0.840039
 0.883291  0.758317  0.660366  1.20652      0.890935  0.584826  0.691214
 0.298193  1.22886   0.60293   1.57381      1.08244   0.462226  1.58976
 ⋮                                       ⋱                      
 0.840968  0.84119   1.10298   0.778372     0.654766  0.461317  0.76752
 0.931224  0.68087   0.549049  0.484713     1.43144   0.627723  0.267985
 0.351612  0.310498  0.975875  0.843578     1.16872   1.29762   0.729312
 1.37635   0.510036  1.50359   1.48308      1.2985    0.79823   0.774335
 1.41766   0.566584  0.609064  0.492235  …  0.735346  0.145164  0.604026
 0.999136  1.06297   0.936562  0.660093     1.0534    1.39612   0.868007
 0.615487  0.504845  1.18105   0.849568     0.160437  1.15834   1.34111
 0.5038    1.14377   1.06801   0.857203     1.06988   0.765729  1.17879
 0.578855  1.22678   0.919697  1.1925       0.538303  0.809692  0.89885

Also, the @. macro applys a dot to every operator:

test5(A,B,C) = @. A + B + C #only one array allocated
@benchmark test5(A,B,C)
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  2
  --------------
  minimum time:     1.100 ms (0.00% GC)
  median time:      1.123 ms (0.00% GC)
  mean time:        1.217 ms (7.67% GC)
  maximum time:     2.102 ms (44.50% GC)
  --------------
  samples:          4089
  evals/sample:     1

Using these tools we can get rid of our intermediate array allocations for many vectorized function calls. But we are still allocating the output array. To get rid of that allocation, we can instead use mutation. Mutating broadcast is done via .=. For example, if we pre-allocate the output:

D = zeros(1000,1000);
1000×1000 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 ⋮                        ⋮              ⋱            ⋮                   
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0
.0

Then we can keep re-using this cache for subsequent calculations. The mutating broadcasting form is:

test6!(D,A,B,C) = D .= A .+ B .+ C #only one array allocated
@benchmark test6!(D,A,B,C)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.044 ms (0.00% GC)
  median time:      1.051 ms (0.00% GC)
  mean time:        1.053 ms (0.00% GC)
  maximum time:     1.328 ms (0.00% GC)
  --------------
  samples:          4725
  evals/sample:     1

If we use @. before the =, then it will turn it into .=:

test7!(D,A,B,C) = @. D = A + B + C #only one array allocated
@benchmark test7!(D,A,B,C)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.042 ms (0.00% GC)
  median time:      1.048 ms (0.00% GC)
  mean time:        1.050 ms (0.00% GC)
  maximum time:     1.334 ms (0.00% GC)
  --------------
  samples:          4738
  evals/sample:     1

Notice that in this case, there is no "output", and instead the values inside of D are what are changed (like with the DiffEq inplace function). Many Julia functions have a mutating form which is denoted with a !. For example, the mutating form of the map is map!:

test8!(D,A,B,C) = map!((a,b,c)->a+b+c,D,A,B,C)
@benchmark test8!(D,A,B,C)
BenchmarkTools.Trial: 
  memory estimate:  32 bytes
  allocs estimate:  1
  --------------
  minimum time:     2.283 ms (0.00% GC)
  median time:      2.293 ms (0.00% GC)
  mean time:        2.296 ms (0.00% GC)
  maximum time:     2.554 ms (0.00% GC)
  --------------
  samples:          2174
  evals/sample:     1

Some operations require using an alternate mutating form in order to be fast. For example, matrix multiplication via * allocates a temporary:

@benchmark A*B
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  2
  --------------
  minimum time:     6.649 ms (0.00% GC)
  median time:      6.747 ms (0.00% GC)
  mean time:        6.827 ms (1.30% GC)
  maximum time:     7.731 ms (12.68% GC)
  --------------
  samples:          732
  evals/sample:     1

Instead, we can use the mutating form mul! into a cache array to avoid allocating the output:

using LinearAlgebra
@benchmark mul!(D,A,B) # same as D = A * B
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     6.171 ms (0.00% GC)
  median time:      6.207 ms (0.00% GC)
  mean time:        6.210 ms (0.00% GC)
  maximum time:     6.837 ms (0.00% GC)
  --------------
  samples:          805
  evals/sample:     1

For repeated calculations this reduced allocation can stop GC cycles and thus lead to more efficient code. Additionally, ***we can fuse together higher level linear algebra operations using BLAS***. The package SugarBLAS.jl makes it easy to write higher level operations like alpha*B*A + beta*C as mutating BLAS calls.

Example Optimization: Gierer-Meinhardt Reaction-Diffusion PDE Discretization

Let's optimize the solution of a Reaction-Diffusion PDE's discretization. In its discretized form, this is the ODE:

\[ \begin{align} du &= D_1 (A_y u + u A_x) + \frac{au^2}{v} + \bar{u} - \alpha u\\ dv &= D_2 (A_y v + v A_x) + a u^2 + \beta v \end{align} \]

where $u$, $v$, and $A$ are matrices. Here, we will use the simplified version where $A$ is the tridiagonal stencil $[1,-2,1]$, i.e. it's the 2D discretization of the LaPlacian. The native code would be something along the lines of:

# Generate the constants
p = (1.0,1.0,1.0,10.0,0.001,100.0) # a,α,ubar,β,D1,D2
N = 100
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
Ax[2,1] = 2.0
Ax[end-1,end] = 2.0
Ay[1,2] = 2.0
Ay[end,end-1] = 2.0

function basic_version!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = r[:,:,1]
  v = r[:,:,2]
  Du = D1*(Ay*u + u*Ax)
  Dv = D2*(Ay*v + v*Ax)
  dr[:,:,1] = Du .+ a.*u.*u./v .+ ubar .- α*u
  dr[:,:,2] = Dv .+ a.*u.*u .- β*v
end

a,α,ubar,β,D1,D2 = p
uss = (ubar+β)/α
vss = (a/β)*uss^2
r0 = zeros(100,100,2)
r0[:,:,1] .= uss.+0.1.*rand.()
r0[:,:,2] .= vss

prob = ODEProblem(basic_version!,r0,(0.0,0.1),p)
ODEProblem with uType Array{Float64,3} and tType Float64. In-place: true
timespan: (0.0, 0.1)
u0: [11.028183604723862 11.034140483848649 … 11.04431642440116 11.004246829
655633; 11.078060648317335 11.038529617964667 … 11.02954445637886 11.010586
806215018; … ; 11.022685024396019 11.066408552141311 … 11.064943225759754 1
1.020455881763404; 11.072762015164164 11.074104020463531 … 11.0357794065046
16 11.096153551663113]

[12.100000000000001 12.100000000000001 … 12.100000000000001 12.100000000000
001; 12.100000000000001 12.100000000000001 … 12.100000000000001 12.10000000
0000001; … ; 12.100000000000001 12.100000000000001 … 12.100000000000001 12.
100000000000001; 12.100000000000001 12.100000000000001 … 12.100000000000001
 12.100000000000001]

In this version we have encoded our initial condition to be a 3-dimensional array, with u[:,:,1] being the A part and u[:,:,2] being the B part.

@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  186.88 MiB
  allocs estimate:  8551
  --------------
  minimum time:     47.767 ms (2.97% GC)
  median time:      49.226 ms (5.63% GC)
  mean time:        49.188 ms (5.56% GC)
  maximum time:     49.334 ms (5.61% GC)
  --------------
  samples:          102
  evals/sample:     1

While this version isn't very efficient,

We recommend writing the "high-level" code first, and iteratively optimizing it!

The first thing that we can do is get rid of the slicing allocations. The operation r[:,:,1] creates a temporary array instead of a "view", i.e. a pointer to the already existing memory. To make it a view, add @view. Note that we have to be careful with views because they point to the same memory, and thus changing a view changes the original values:

A = rand(4)
@show A
B = @view A[1:3]
B[2] = 2
@show A
A = [0.2409747399305775, 0.7345305736521792, 0.5717443440936183, 0.20656756
482662297]
A = [0.2409747399305775, 2.0, 0.5717443440936183, 0.20656756482662297]
4-element Array{Float64,1}:
 0.2409747399305775
 2.0
 0.5717443440936183
 0.20656756482662297

Notice that changing B changed A. This is something to be careful of, but at the same time we want to use this since we want to modify the output dr. Additionally, the last statement is a purely element-wise operation, and thus we can make use of broadcast fusion there. Let's rewrite basic_version! to ***avoid slicing allocations*** and to ***use broadcast fusion***:

function gm2!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  Du = D1*(Ay*u + u*Ax)
  Dv = D2*(Ay*v + v*Ax)
  @. du = Du + a.*u.*u./v + ubar - α*u
  @. dv = Dv + a.*u.*u - β*v
end
prob = ODEProblem(gm2!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  119.55 MiB
  allocs estimate:  7081
  --------------
  minimum time:     38.166 ms (3.92% GC)
  median time:      38.400 ms (3.95% GC)
  mean time:        38.856 ms (5.13% GC)
  maximum time:     40.026 ms (7.45% GC)
  --------------
  samples:          129
  evals/sample:     1

Now, most of the allocations are taking place in Du = D1*(Ay*u + u*Ax) since those operations are vectorized and not mutating. We should instead replace the matrix multiplications with mul!. When doing so, we will need to have cache variables to write into. This looks like:

Ayu = zeros(N,N)
uAx = zeros(N,N)
Du = zeros(N,N)
Ayv = zeros(N,N)
vAx = zeros(N,N)
Dv = zeros(N,N)
function gm3!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  mul!(Ayu,Ay,u)
  mul!(uAx,u,Ax)
  mul!(Ayv,Ay,v)
  mul!(vAx,v,Ax)
  @. Du = D1*(Ayu + uAx)
  @. Dv = D2*(Ayv + vAx)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm3!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  29.75 MiB
  allocs estimate:  5317
  --------------
  minimum time:     33.330 ms (0.00% GC)
  median time:      33.477 ms (0.00% GC)
  mean time:        33.944 ms (1.39% GC)
  maximum time:     35.020 ms (4.19% GC)
  --------------
  samples:          148
  evals/sample:     1

But our temporary variables are global variables. We need to either declare the caches as const or localize them. We can localize them by adding them to the parameters, p. It's easier for the compiler to reason about local variables than global variables. ***Localizing variables helps to ensure type stability***.

p = (1.0,1.0,1.0,10.0,0.001,100.0,Ayu,uAx,Du,Ayv,vAx,Dv) # a,α,ubar,β,D1,D2
function gm4!(dr,r,p,t)
  a,α,ubar,β,D1,D2,Ayu,uAx,Du,Ayv,vAx,Dv = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  mul!(Ayu,Ay,u)
  mul!(uAx,u,Ax)
  mul!(Ayv,Ay,v)
  mul!(vAx,v,Ax)
  @. Du = D1*(Ayu + uAx)
  @. Dv = D2*(Ayv + vAx)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm4!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  29.66 MiB
  allocs estimate:  1053
  --------------
  minimum time:     26.884 ms (0.00% GC)
  median time:      27.187 ms (0.00% GC)
  mean time:        27.634 ms (1.70% GC)
  maximum time:     28.720 ms (4.98% GC)
  --------------
  samples:          181
  evals/sample:     1

We could then use the BLAS gemmv to optimize the matrix multiplications some more, but instead let's devectorize the stencil.

p = (1.0,1.0,1.0,10.0,0.001,100.0,N)
function fast_gm!(du,u,p,t)
  a,α,ubar,β,D1,D2,N = p

  @inbounds for j in 2:N-1, i in 2:N-1
    du[i,j,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end

  @inbounds for j in 2:N-1, i in 2:N-1
    du[i,j,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds for j in 2:N-1
    i = 1
    du[1,j,1] = D1*(2u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
            a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for j in 2:N-1
    i = 1
    du[1,j,2] = D2*(2u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            a*u[i,j,1]^2 - β*u[i,j,2]
  end
  @inbounds for j in 2:N-1
    i = N
    du[end,j,1] = D1*(2u[i-1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
           a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for j in 2:N-1
    i = N
    du[end,j,2] = D2*(2u[i-1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
           a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds for i in 2:N-1
    j = 1
    du[i,1,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for i in 2:N-1
    j = 1
    du[i,1,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
              a*u[i,j,1]^2 - β*u[i,j,2]
  end
  @inbounds for i in 2:N-1
    j = N
    du[i,end,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for i in 2:N-1
    j = N
    du[i,end,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds begin
    i = 1; j = 1
    du[1,1,1] = D1*(2u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[1,1,2] = D2*(2u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
              a*u[i,j,1]^2 - β*u[i,j,2]

    i = 1; j = N
    du[1,N,1] = D1*(2u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[1,N,2] = D2*(2u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]

    i = N; j = 1
    du[N,1,1] = D1*(2u[i-1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[N,1,2] = D2*(2u[i-1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]

    i = N; j = N
    du[end,end,1] = D1*(2u[i-1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[end,end,2] = D2*(2u[i-1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]
   end
end
prob = ODEProblem(fast_gm!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  29.62 MiB
  allocs estimate:  466
  --------------
  minimum time:     8.027 ms (0.00% GC)
  median time:      8.179 ms (0.00% GC)
  mean time:        8.618 ms (5.18% GC)
  maximum time:     9.634 ms (14.07% GC)
  --------------
  samples:          580
  evals/sample:     1

Lastly, we can do other things like multithread the main loops, but these optimizations get the last 2x-3x out. The main optimizations which apply everywhere are the ones we just performed (though the last one only works if your matrix is a stencil. This is known as a matrix-free implementation of the PDE discretization).

This gets us to about 8x faster than our original MATLAB/SciPy/R vectorized style code!

The last thing to do is then ***optimize our algorithm choice***. We have been using Tsit5() as our test algorithm, but in reality this problem is a stiff PDE discretization and thus one recommendation is to use CVODE_BDF(). However, instead of using the default dense Jacobian, we should make use of the sparse Jacobian afforded by the problem. The Jacobian is the matrix $\frac{df_i}{dr_j}$, where $r$ is read by the linear index (i.e. down columns). But since the $u$ variables depend on the $v$, the band size here is large, and thus this will not do well with a Banded Jacobian solver. Instead, we utilize sparse Jacobian algorithms. CVODE_BDF allows us to use a sparse Newton-Krylov solver by setting linear_solver = :GMRES (see the solver documentation, and thus we can solve this problem efficiently. Let's see how this scales as we increase the integration time.

prob = ODEProblem(fast_gm!,r0,(0.0,10.0),p)
@benchmark solve(prob,Tsit5())
BenchmarkTools.Trial: 
  memory estimate:  2.76 GiB
  allocs estimate:  41632
  --------------
  minimum time:     1.491 s (47.78% GC)
  median time:      1.748 s (27.82% GC)
  mean time:        1.714 s (32.45% GC)
  maximum time:     1.904 s (24.70% GC)
  --------------
  samples:          3
  evals/sample:     1
using Sundials
@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES))
BenchmarkTools.Trial: 
  memory estimate:  44.93 MiB
  allocs estimate:  7051
  --------------
  minimum time:     277.657 ms (0.00% GC)
  median time:      279.860 ms (0.29% GC)
  mean time:        279.719 ms (0.31% GC)
  maximum time:     281.571 ms (0.61% GC)
  --------------
  samples:          18
  evals/sample:     1
prob = ODEProblem(fast_gm!,r0,(0.0,100.0),p)
# Will go out of memory if we don't turn off `save_everystep`!
@benchmark solve(prob,Tsit5(),save_everystep=false)
BenchmarkTools.Trial: 
  memory estimate:  2.90 MiB
  allocs estimate:  74
  --------------
  minimum time:     5.685 s (0.00% GC)
  median time:      5.685 s (0.00% GC)
  mean time:        5.685 s (0.00% GC)
  maximum time:     5.685 s (0.00% GC)
  --------------
  samples:          1
  evals/sample:     1
@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES))
BenchmarkTools.Trial: 
  memory estimate:  267.86 MiB
  allocs estimate:  48458
  --------------
  minimum time:     1.992 s (0.00% GC)
  median time:      2.025 s (1.39% GC)
  mean time:        2.042 s (2.43% GC)
  maximum time:     2.109 s (5.72% GC)
  --------------
  samples:          3
  evals/sample:     1

Now let's check the allocation growth.

@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)
BenchmarkTools.Trial: 
  memory estimate:  3.23 MiB
  allocs estimate:  40351
  --------------
  minimum time:     1.962 s (0.00% GC)
  median time:      1.962 s (0.00% GC)
  mean time:        1.962 s (0.00% GC)
  maximum time:     1.964 s (0.00% GC)
  --------------
  samples:          3
  evals/sample:     1
prob = ODEProblem(fast_gm!,r0,(0.0,500.0),p)
@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)
BenchmarkTools.Trial: 
  memory estimate:  4.19 MiB
  allocs estimate:  57133
  --------------
  minimum time:     2.783 s (0.00% GC)
  median time:      2.785 s (0.00% GC)
  mean time:        2.785 s (0.00% GC)
  maximum time:     2.787 s (0.00% GC)
  --------------
  samples:          2
  evals/sample:     1

Notice that we've elimated almost all allocations, allowing the code to grow without hitting garbage collection and slowing down.

Why is CVODE_BDF doing well? What's happening is that, because the problem is stiff, the number of steps required by the explicit Runge-Kutta method grows rapidly, whereas CVODE_BDF is taking large steps. Additionally, the GMRES linear solver form is quite an efficient way to solve the implicit system in this case. This is problem-dependent, and in many cases using a Krylov method effectively requires a preconditioner, so you need to play around with testing other algorithms and linear solvers to find out what works best with your problem.

Conclusion

Julia gives you the tools to optimize the solver "all the way", but you need to make use of it. The main thing to avoid is temporary allocations. For small systems, this is effectively done via static arrays. For large systems, this is done via in-place operations and cache arrays. Either way, the resulting solution can be immensely sped up over vectorized formulations by using these principles.

Appendix

This tutorial is part of the SciMLTutorials.jl repository, found at: https://github.com/SciML/SciMLTutorials.jl. For more information on doing scientific machine learning (SciML) with open source software, check out https://sciml.ai/.

To locally run this tutorial, do the following commands:

using SciMLTutorials
SciMLTutorials.weave_file("introduction","03-optimizing_diffeq_code.jmd")

Computer Information:

Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
  JULIA_LOAD_PATH = /builds/JuliaGPU/DiffEqTutorials.jl:
  JULIA_DEPOT_PATH = /builds/JuliaGPU/DiffEqTutorials.jl/.julia
  JULIA_CUDA_MEMORY_LIMIT = 2147483648
  JULIA_NUM_THREADS = 8

Package Information:

Status `/builds/JuliaGPU/DiffEqTutorials.jl/tutorials/introduction/Project.toml`
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.5.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.15.0
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 5.5.0
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 1.6.3
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.12.4
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 4.3.0
[37e2e46d-f89d-539d-b4ee-838fcccc9c8e] LinearAlgebra