[BoundaryValueDiffEq.jl] Reducing allocations

New to Julia
1

Hi,

I started to implement my Python pipeline for solving boundary value problems in Julia using BoundaryValueDiffEq.jl. I tried to optimize the functions but I am wondering if it is possible to increase performance further by reducing the number of allocations. In the minimal reproducible example below I can’t get below 440 k allocations. Is there any additional optimization I could try or are the numbers of allocations I am getting as good as it gets?

import BoundaryValueDiffEq as BVP


function distribution_fun_test!(dydt::AbstractVector{T}, y::AbstractVector{T},  p::Any, #
                                t::Float64) where {T<:Real}

    y1, y2, y3, y4 = y

    dydt[1] = 0.001 * y2
    dydt[2] = 10 * y1
    dydt[3] = -0.03 * y2 + 0.0001 * y4
    dydt[4] = 10 * y3
                    
    uptake = 0.1 * y3 / (1000 * y3 + 1)
    dydt[3] += uptake
end


function distribution_bc_test!(residual::AbstractVector{}, u::Any, p::Any, t::Vector{Float64})

    residual[1] = (u[4, end] - 0.001)
    residual[2] = u[1, end]
    residual[3] =   (- u[3, end])
    residual[4] =   (u[1, 1])

end

function test_solve()
    L = 100.0
    mesh_size = 500
    tspan::Tuple{Float64, Float64} = (0.0, Float64(L))

    y_initial::Vector{Float64} = [1, 0.001/0.03, 1, 0.001]

    bvp_problem = BVP.BVProblem(distribution_fun_test!, 
                                distribution_bc_test!, 
                                y_initial, tspan)
        
    @time BVP.solve(bvp_problem, BVP.MIRK5(); adaptive = false, dt=Float64(L/mesh_size))
    @time BVP.solve(bvp_problem, BVP.MIRK5(); adaptive = false, dt=Float64(L/mesh_size))
end

test_solve()

The output of the timing I am getting is:
3.606669 seconds (4.35 M allocations: 216.794 MiB, 1.00% gc time, 96.31% compilation time: 100% of which was recompilation) 0.132779 seconds (440.05 k allocations: 26.994 MiB)

I am using BoundaryValueDiffEq v5.18.0 and Julia v1.12.

Thanks a lot for your help in advance!
Max

2

Welcome to the Julia community.

It looks like… um… I think the performance bottleneck is in the BoundaryValueDiffEq package (rather than the code you wrote). The type instability is also caused by this package. :thinking:

3

Thanks! And thanks for your reply!

It’s a pitty that there seems to be nothing I can change in my functions.

Maybe there is a way to optimize how I set up BVProblem instead?
I tried setting the the parameter nlls to false or true for type stability as suggested in the docs - but this seemed to increase allocations and increased execution time.
I also tried returning a staticarray for both functions instead of in place changes (I supplied y_array as a staticarray, too). But this increased allocations massively.

As these now are quite specific points for the boundaryvaluediffeq package - should I post the question somewhere else?

4

You can profile exactly where the allocations are coming from by using @profview_allocs (if you’re using VSCode)

5

I am using VSCode, however, for some reason @profview_allocs did not work for me. I used the alternative PProf instead. But to be honest I couldn’t really find out which other options to test based on the output, as I am not familiar with the package structure of BoundaryValueDiffEq.jl.

It’s a pitty - this Julia implementation is not faster than my previous Python implementation. I expected the performance to be noticeably better in Julia after some optimization.

6

BoundaryValueDiffEq.jl is still getting a ton of work, with its core adaptivity algorithms and such just completed last summer. I’ll take a look at some of the performance stuff here. Note that generally the BVP solvers will need to allocate in order to adapt, and in almost all cases you want to have that adaptivity on. But as a smoke test we should get the non-adapting steps to be non-allocating, so I’ll spin up some tests on that to see what’s going on here.

This is better as an issuei n the package.

What type instability?

7

I found that the return type of both bvp_problem and BVP.solve is Any using Cthulhu’s @descend.

> @descend test_solve()
test_solve() @ Main REPL[5]:1
 1 function test_solve()::Any
 2     L::Core.Const(100.0) = 100.0::Core.Const(100.0)
 3     mesh_size::Core.Const(500) = 500::Core.Const(500)
 4     (tspan::Core.Const((0.0, 100.0))::Tuple{Float64, Float64})::Core.Const(true) = (0.0, Float64(L::Core.Const(100.0))::Core.Const(100.0))::Core.Const((0.0, 100.0))
 6 
 6     (y_initial::Vector{Float64}::Vector{Float64})::Core.Const(true) = [1, (0.001/0.03)::Core.Const(0.03333333333333333), 1, 0.001]::Vector{Float64}
 8 
 8     bvp_problem::Any = BVP.BVProblem::Type{SciMLBase.BVProblem}(distribution_fun_test!::Core.Const(Main.distribution_fun_test!), 
 9                                 distribution_bc_test!::Core.Const(Main.distribution_bc_test!), 
10                                 y_initial::Vector{Float64}, tspan::Core.Const((0.0, 100.0)))::Any
11         
12     @time BVP.solve::Core.Const(CommonSolve.solve)(bvp_problem::Any, BVP.MIRK5::Type{BoundaryValueDiffEqMIRK.MIRK5}()::Core.Const(MIRK5(jac_alg = BVPJacobianAlgorithm(), defect_threshold = 0.1, max_num_subintervals = 3000)); adaptive = false, dt=Float64((L::Core.Const(100.0)/mesh_size::Core.Const(500))::Float64)::Tuple{Any, Float64, Int64, Float64, Base.GC_Diff, Int64, Float64, Float64}::NamedTuple{(:value, :time, :bytes, :gctime, :gcstats, :lock_conflicts, :compile_time, :recompile_time), <:Tuple{Any, Float64, Int64, Float64, Base.GC_Diff, Int64, Float64, Float64}})::Any
13     @time BVP.solve::Core.Const(CommonSolve.solve)(bvp_problem::Any, BVP.MIRK5::Type{BoundaryValueDiffEqMIRK.MIRK5}()::Core.Const(MIRK5(jac_alg = BVPJacobianAlgorithm(), defect_threshold = 0.1, max_num_subintervals = 3000)); adaptive = false, dt=Float64((L::Core.Const(100.0)/mesh_size::Core.Const(500))::Float64)::Tuple{Any, Float64, Int64, Float64, Base.GC_Diff, Int64, Float64, Float64}::NamedTuple{(:value, :time, :bytes, :gctime, :gcstats, :lock_conflicts, :compile_time, :recompile_time), <:Tuple{Any, Float64, Int64, Float64, Base.GC_Diff, Int64, Float64, Float64}})::Any
14 end
8

Can you open an issue for that? That should get addressed.

1 Like
9

Thanks a lot for your reply, Chris. And thanks a lot for developing and maintaining such awesome Julia packages!

Very cool to know that BoundaryValueDiffEq.jl is still being worked on so much.

I just opened an issue.

For my BVP I did find that turning adaptivity on slowed down the solve a bit but did not seem to change the result. I find that I seem to have to set the grid size manually since very small grid sizes lead to inaccurate results - both with “adaptive” turned on or off.

10

Yes it’s still a highly active project. It’s the solver library that doesn’t have a paper out yet even, so it’ll be the next “big” SciML solver paper over the next year. The JuliaCon video describes a lot of what was done, so now lots of the algorithm is done but it needs lots of the little bits to complete it. So it needs issues, and then just lots of the little optimizations and bug fixes and interface things.