Nested QuadGK: Long compilation time gets short when plotting the argument beforehand, but not when calling it

General Usage
1

Hey everyone, I stumbled upon a weird behaviour of compilation time.
I have the following function definitions:

using QuadGK                                                                                                                                                                                                                                                        
integrate(f,t) = quadgk(f,0,t)[1]                                                                                                                                                                                                                                                                                                                                       
g(t) = t^2                                                                                
f(t) = integrate(g,t)

Now if I call quadgk(f,0,1), it takes over two minutes to compile:

julia> @time quadgk(f,0,1);
142.877898 seconds (305.06 M allocations: 23.536 GiB, 4.36% gc time)

julia> @time quadgk(f,0,1);
  0.000242 seconds (552 allocations: 13.219 KiB)

If instead I plot f before, the compilation time is only half a second:

julia> @time using Plots; @time plot(f,0,1);
  4.834115 seconds (9.95 M allocations: 550.155 MiB, 4.37% gc time)
 22.190020 seconds (56.53 M allocations: 2.786 GiB, 7.81% gc time)

julia> @time quadgk(f,0,1);
  0.598319 seconds (1.36 M allocations: 83.927 MiB, 6.51% gc time)

julia> @time quadgk(f,0,1);
  0.000156 seconds (537 allocations: 12.984 KiB)

Just calling e.g f(0.1) does not have an impact on the compilation time.
Can anyone explain to me what’s going on here? What does plot(f,0,1) do that f(0.1) doesn’t?

I am using Julia 1.1.0

2 Likes
2

I am also surprised the first compilation @time quadgk(f,0,1) takes 140 seconds (150 on my machine).

3

I don’t have an explanation for the compilation time — you might file an issue at JuliaLang/julia — but I would strongly recommend against using nested 1d adaptive quadrature, which can easily waste a lot of time trying to refine 1d integrals in regions that don’t contribute much to the overall integral.

Instead, you should use a “cubature” routine (like in HCubature.jl, Cubature.jl, or Cuba.jl) that knows how to do multiple integrals in a holistic way. Even though these routines typically integrate over box domains, you can integrate over other domains with a change of variables. For example, your problem is a triangular domain, which can be mapped to a box with:

\int_0^1 dx \int_0^x dy \, F(x,y) = \int_0^1 dx \int_0^1 du \, xF(x,ux)

via the change of variables y = ux.

5 Likes