ITP implementation beaten by Roots.jl

[tecosaur]

Hello!

I wrote the original code that was developed into the current ITP implementation (#184) and recently for fun thought I'd compare it to what's in Roots.jl using the example problem in https://docs.sciml.ai/SciMLBenchmarksOutput/stable/IntervalNonlinearProblem/simpleintervalrootfind/.

Curiously, it seems like the Roots.jl implementation is ~15% faster. Adding an extra Roots.ITP() argument to find_zero, I find that (using Chairmarks and a Zen2 chip):

julia> @b froots(out, levels, (0.0, 2.0))
81.623 ms (100000 allocs: 15.259 MiB, 10.37% gc time)

julia> @b f(out, levels, (0.0, 2.0))
94.051 ms

From this, I suspect there's some performance on the table with the current implementation. I don't have time to dig into this myself, but thought I might as well let you folk know 🙂.

[ChrisRackauckas]

Huh interesting... and there's an alloc there? We should definitely dive in.

[oscardssmith]

The allocation is in the Roots.jl version. The surprising thing is that despite that, Roots.jl is faster than NonlinearSolve,

[oscardssmith]

Interestingly, this doesn't reproduce for me. I'm seeing better results on NonlinearSolve... @tecosaur what versions are you running?

julia> @benchmark f(out, levels, (0.0, 2.0))
BenchmarkTools.Trial: 44 samples with 1 evaluation per sample.
 Range (min … max):  113.582 ms … 121.557 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     114.558 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   114.902 ms ±   1.335 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

        █▆▂                                                      
  ▆▄▄▆█▆███▆▄█▆▆▆▁▁▄▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▁
  114 ms           Histogram: frequency by time          122 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark froots(out, levels, (0, 2))
BenchmarkTools.Trial: 32 samples with 1 evaluation per sample.
 Range (min … max):  153.048 ms … 162.735 ms  ┊ GC (min … max): 0.00% … 4.06%
 Time  (median):     156.048 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   156.316 ms ±   2.438 ms  ┊ GC (mean ± σ):  0.34% ± 0.84%

         █ ▃         ▃   █                                    ▃  
  ▇▁▇▁▇▇▇█▇█▁▁▁▇▇▇▇▁▇█▁▇▇█▇▁▇▁▁▁▇▁▇▁▁▇▁▁▁▁▁▇▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  153 ms           Histogram: frequency by time          163 ms <

 Memory estimate: 7.63 MiB, allocs estimate: 100000.

[oscardssmith]

That said, #547 will save a couple ns

[tecosaur]

The allocation is in the Roots.jl version. The surprising thing is that despite that, Roots.jl is faster than NonlinearSolve,

Indeed, that was exactly what caught my eye.

Interestingly, this doesn't reproduce for me. I'm seeing better results on NonlinearSolve... @tecosaur what versions are you running?

Let's see, that comment was from late last year so I suspect 1.10 and it would have been on an x86_64 Linux system. I'm guessing from the difference you see now that compiler changes may have helped the ITP implementation here.

[oscardssmith]

Yeah, that's plausible. The code is fairly textbook, and is relying on the compiler to do some fixes (e.g. replace division by 2 with multiplication by .5. In any case, since this aparently no longer reproduces, I think it's safe to close.

[ChrisRackauckas]

I also cannot reproduce:

using Roots, BracketingNonlinearSolve, BenchmarkTools

const N = 100_000;
levels = 1.5 .* rand(N);
out = zeros(N);
myfun(x, lv) = x * sin(x) - lv
function froots(out, levels, u0)
    for i in 1:N
        out[i] = find_zero(x->myfun(x, levels[i]), u0)
    end
end

@btime froots(out, levels, (0, 2)) # 182.550 ms (100000 allocations: 7.63 MiB)

function f(out, levels, u0)
    for i in 1:N
        out[i] = solve(IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun),
                u0, levels[i]), ITP()).u
    end
end

@btime f(out, levels, (0.0, 2.0)) # 110.998 ms (0 allocations: 0 bytes)

So I'll call this case closed.

[tecosaur]

Brilliant, glad to see everything's working as expected now 🙂