Need help understanding the allocations

General Usage Performance
,
1

Consider the following code

using BenchmarkTools
using Random
using QuadGK

f(p, xm, x) = -x^2 / ((p[1] - xm[1]*x)^2 + (p[2] - xm[2]*x)^2)
g(p, xm, x) = -x^2 / ((p[1] - xm[1])^2 + (p[2] - xm[2])^2)
integralf(p, xm) = quadgk(x -> f(p, xm, x), 0.0, 1.0)
integralg(p, xm) = quadgk(x -> g(p, xm, x), 0.0, 1.0)

function test(p, xm, integral)
    s = 0.0
    for i in 1:1000
        rand!(p)
        rand!(xm)
        res, err = integral(p, xm)
        s += res
    end
    return s
end
p  = zeros(Float64, 2)
xm = zeros(Float64, 2)
test(p, xm, integralf)
test(p, xm, integralg)

When I run the benchmarks, integralf is allocating, whereas integralg ( which only has x^2 in it) is not, which I do not understand.

@benchmark test($p, $xm, $integralf)

image
image704×185 18.9 KB

@benchmark test($p, $xm, $integralg)
image

Can someone help me understanding what is going on?

EDIT: I think it is related to the fact that there are 1/(a-x) kind of terms in the integrand, as shown by allocation profiler. With the division replaced by multiplication in f(p, xm, x), it vanishes.

2

One thing is that you better define your test function like this:

julia> function test(p, xm, integral::F) where {F<:Function}
           s = 0.0
           for i in 1:1000
               rand!(p)
               rand!(xm)
               res, err = integral(p, xm)
               s += res
           end
           return s
       end

because Julia by default may not specialize to function arguments, and this can cause allocations (although here you are using the function, so it should).

That being said, here this does not solve the problem. And since the direct call to each integral is not allocating, I think you reached some weird bug.

3

I believe the allocations come from the algorithm performing further subdivisions in the case of f, while it can early return for g, because the error estimate is sufficiently small already.

Step through it with Debugger and you’ll see what I mean.

2 Likes
4

Or use staticarrays

5

Running the integral functions outside the loop doesn’t allocate anything (can that be dependent on the random point given?)

6

It can be different. Quadgk is looking to see if the error bound is being met and that will be parameter dependent.

1 Like
7

Indeed, sometimes it allocates, sometimes it doesn’t:

julia> @btime integralf($(rand(2)), $(rand(2)))
  55.014 ns (0 allocations: 0 bytes)
(-0.3320182264630196, 0.0)

julia> @btime integralf($(rand(2)), $(rand(2)))
  307.376 ns (2 allocations: 368 bytes)
(-2.1198229226184617, 9.467969463994308e-10)

julia> @btime integralf($(rand(2)), $(rand(2)))
  517.234 ns (2 allocations: 368 bytes)
(-4.848361174937609, 1.5858105731347827e-8)

julia> @btime integralf($(rand(2)), $(rand(2)))
  627.907 ns (2 allocations: 368 bytes)
(-18.88726883755866, 1.4078685466045737e-7)

julia> @btime integralf($(rand(2)), $(rand(2)))
  57.118 ns (0 allocations: 0 bytes)
(-0.4933392928396032, 7.264673307361136e-10)

Does not work here:

julia> @btime integralf($(rand(SVector{2,Float64})), $(rand(SVector{2,Float64})))
  307.844 ns (2 allocations: 368 bytes)
(-2.7657212478611743, 4.964859356970663e-10)
8

Indeed as @skleinbo and @Oscar_Smith pointed out, it is because quadgk is using adaptive methods to meet the error bounds on the integral. The allocation profiler output also says so.

image
image1190×234 43.9 KB

Changing the rtol parameter for quadgk changes the total allocations. Higher the rtol, lower the allocations for integralf.
integralf(p, xm) = quadgk(x -> f(p, xm, x), 0.0, 1.0, rtol=1.e-4)
rtol=1.e-4 : Memory estimate: 112.02 KiB, allocs estimate: 609.
rtol=1.e-2 : Memory estimate: 42.19 KiB, allocs estimate: 222.

Is there a way to preallocate memory that quadgk can reuse? The integrand and the limits do not change, only the parameters p and xm do. I tried segbuf option but it increases the allocations.

9

Yes, you can pass that via the segbuf keyword. See the docstrings of quadgk and alloc_segbuf for how to do it.

However, I presume the algorithm is subdividing a lot because the integrand has a singularity, which may not be what you intend.