# Performance of large sparse discrete map problems in DifferentialEquations

**URL:** <https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382>\
**Category:** Performance\
**Tags:** diffeq, ode, sparse\
**Created:** [February 14, 2022, 2:42am UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382 "2022-02-14T02:42:12Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![roryh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roryh/32/23094_2.png) [@roryh](https://discourse.julialang.org/u/roryh)\
**Post date:** [February 14, 2022, 2:42am UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/1 "2022-02-14T02:42:12Z")

</div>

Hello,

I normally deal with dynamics on networks, so I have giant sparse systems of coupled ODE’s. The problem I have at the moment is that for a simple discrete map, a naive implementation is significantly faster than solving the DiscreteProblem provided by DIfferentialEquations.jl. For example, the code below is a MWE of what my problems look like,

```julia
using DifferentialEquations
using Random
using SparseArrays

# parametes
rng = MersenneTwister(0)
N = 80000
M = 150000
T = 10
β = 0.4
A = sprand(rng, N, M, 0.0001, (r, x) -> rand(r, [1], x))
p = (β, A)
###

function map!(unp1, un, p, t)
    β, A = p
    unp1 .= β .* un .* A
end

function naive_test(u0, tspan, p)
    unp1 = similar(u0)
    for i in tspan[1]:tspan[2]
        map!(unp1, u0, p, i)
    end
    return unp1
end

# test
u0 = sprand(N, M, 0.0001)
prob = DiscreteProblem(map!, u0, (0, T), p)

@time solve(prob; save_on=false);

@time naive_test(u0, (0, T), p);

```

For the `solve` function I turn off saving of intermediate results to better match the naive test and I get the following,

```julia
julia> @time solve(prob; save_on=false);
205.226954 seconds (16.25 M allocations: 1.241 GiB, 0.10% gc time, 1.79% compilation time)

```

and for the `naive_test` implementation I get

```julia
julia> @time naive_test(u0, (0, T), p);
  0.142919 seconds (16.55 k allocations: 20.360 MiB, 7.48% compilation time)

```

I was wondering if I’m missing something and this sort of performance is expected, I know there is a significant amount of machinery behind diffeq and it might not suit sparse arrays. Any help or explanation would be appreciated.

Using DifferentialEquations v7.1.0 and julia 1.7.1.

Rory.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [February 14, 2022, 7:07am UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/2 "2022-02-14T07:07:21Z")

</div>

Reducing the problem size and testing it like so

```julia
using LinearAlgebra
using SimpleDiffEq
using DifferentialEquations
using OrdinaryDiffEq
using Random
using SparseArrays
using BenchmarkTools

# parametes
const rng = MersenneTwister(0)
const N = 8000
const M = 15000
const T = 10
const β = 0.4
const A = sprand(rng, N, M, 0.0001, (r, x) -> rand(r, [1], x))
const p = (β, A)
###

function map!(unp1, un, p, t)
    β, A = p
    unp1 .= β .* un .* A
end

function naive_test(u0, tspan, p)
    unp1 = similar(u0)
    for i in tspan[1]:tspan[2]
        map!(unp1, u0, p, i)
    end
    return unp1
end

# test
const u0 = sprand(N, M, 0.0001)
const prob = DiscreteProblem(map!, u0, (0, T), p)

y1 = solve(prob, FunctionMap(); save_on=false);
y2 = solve(prob, SimpleFunctionMap());
println(norm(y2[end] - y1[end]))
y3 = naive_test(u0, (0, T), p);
println(norm(y3 - y2[end]))

@btime solve(prob, FunctionMap(); save_on=false);
@btime solve(prob, SimpleFunctionMap());
@btime naive_test(u0, (0, T), p);
println()

```

yields something like this for me

```julia
0.0
0.40618462634054964
  2.144 s (58 allocations: 1.72 MiB)
  2.158 ms (67 allocations: 3.28 MiB)
  3.318 ms (6 allocations: 305.64 KiB)

```

I’m wondering about the deviation in results and if the performance difference of `FunctionMap` and `SimpleFunctionMap` is justified in this case?

Edit:

- The deviation in results is probably due to a different number of rounds
- Runtime in `FunctionMap` is dominated by `_any(::typeof(DiffEqBase.NAN_CHECK), ::SparseMatrixCSC{Float64, Int64}, ::Colon)`, which seems like a bad idea…

Edit: Here is a flame graph

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/5/2556f7884a1b4f3a94202621e8b99f390e769acf.png)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 14, 2022, 11:23am UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/3 "2022-02-14T11:23:35Z")

</div>

> [@goerch](#):
>
> - Runtime in `FunctionMap` is dominated by `_any(::typeof(DiffEqBase.NAN_CHECK), ::SparseMatrixCSC{Float64, Int64}, ::Colon)` , which seems like a bad idea…

Well NaN check in and of itself is a good idea: you want to look for NaNs and if you see any then you want to stop calculating and throw the user a warning. What seems to be the issue is that `any(isnan,A::SparseMatrixCSC)` seems like it’s looping over all values, not all non-zero values, which for a very sparse matrix would take considerably longer than any other specialized calculation. We should specialize this in SciMLBase. Could you open an issue there?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 14, 2022, 11:25am UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/4 "2022-02-14T11:25:06Z")

</div>

Try this PR branch: [https://github.com/SciML/DiffEqBase.jl/pull/732](https://github.com/SciML/DiffEqBase.jl/pull/732)

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [February 14, 2022, 12:34pm UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/5 "2022-02-14T12:34:11Z")

</div>

> [@ChrisRackauckas](#):
>
> Try this PR branch: [Specialize NAN\_CHECK for nonzeros of states by ChrisRackauckas · Pull Request #732 · SciML/DiffEqBase.jl · GitHub](https://github.com/SciML/DiffEqBase.jl/pull/732)

Nice, this does the trick for me:

```julia
  2.081 ms (58 allocations: 1.72 MiB)
  2.172 ms (67 allocations: 3.27 MiB)
  3.266 ms (6 allocations: 304.77 KiB)

```

Just a remark:

```julia
@inline NAN_CHECK(x::SparseArrays.AbstractSparseMatrixCSC) = any(NAN_CHECK, nonzeros(x))

```

seems to work, too.

---

<div class="post-metadata">

**Author:** ![roryh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roryh/32/23094_2.png) [@roryh](https://discourse.julialang.org/u/roryh)\
**Post date:** [February 14, 2022, 1:01pm UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/6 "2022-02-14T13:01:02Z")

</div>

Thanks for fixing the error, I edited my original post so it should at least run. I should have run it in a clean repl before posting.

---

<div class="post-metadata">

**Author:** ![roryh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roryh/32/23094_2.png) [@roryh](https://discourse.julialang.org/u/roryh)\
**Post date:** [February 14, 2022, 1:03pm UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/7 "2022-02-14T13:03:28Z")

</div>

Just tested this branch and it works perfectly, I’m getting the same results as @goerch. Thanks a million, I have to say Julia’s biggest strength is it’s community 😁

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 14, 2022, 1:18pm UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/8 "2022-02-14T13:18:57Z")

</div>

> [@goerch](#):
>
> `@inline NAN_CHECK(x::SparseArrays.AbstractSparseMatrixCSC) = any(NAN_CHECK, nonzeros(x))`

Yeah, that’s probably better.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 14, 2022, 2:12pm UTC](https://discourse.julialang.org/t/performance-of-large-sparse-discrete-map-problems-in-differentialequations/76382/9 "2022-02-14T14:12:55Z")

</div>

Merged.