# Problem with small integrals

I am computing 2 integrals over the real line using quadgk from the package of the same name. The goal is to estimate the ratio of the integrals. Often both integrals are quite small, especially in some extreme cases I’m using to test. Getting reasonable values for the numerator has been tricky, and I’m hoping for some advice on how to handle it.

The denominator is the integral of a positive function, and seems to work well with only a relative tolerance. But the numerator might be 0, and so using an absolute tolerance seems advisable. I have tried scaling it by the size of the denominator, but sometimes the result is much smaller than it should be. Removing the absolute tolerance fixes things, but that re-opens the possible bad behavior for true 0.

Sample code:

# illustrates problems with premature termination of integration

"conditional probability of n successes out of n tries"
function conddens(n::Int, z::Float64)
η = -2.0+0.25*z  # transform z to mean -2, sd 0.25
return (1.0+exp(-η))^(-n)
end

"combined normal density and weight function, ignoring constants"
function dens2(z::Float64)
λ = 0.48
return exp((2*λ-1)*z^2/2.0)
end

function trouble(n::Int)
f0(z) = conddens(n, z)*dens2(z)
r0, r0err = quadgk(f0, -Inf, Inf, atol=0.0)

f1(z) = z*conddens(n, z)*dens2(z) #same as z*f0(z)
tol1 = r0*sqrt(eps())
r1, r1err = quadgk(f1, -Inf, Inf, atol=tol1)
r1b, r1berr = quadgk(f1, -Inf, Inf, atol=0.0)

println("Denominator = $(r0) ±$(r0err).")
println("Numerator = $(r1) ±$(r1err) with atol = $(tol1).") println("Numerator =$(r1b) ± $(r1berr) with no atol.") end trouble(100)  Produces Denominator = 1.7081181337842018e-5 ± 9.259586775755586e-14. Numerator = 7.122927197595122e-22 ± 7.494964918478972e-22 with atol = 2.5452943649652627e-13. Numerator = 0.0003717149762667195 ± 3.760194450016882e-12 with no atol.  In these problem cases the mode of the function being integrated can be relatively far out, like at 10 to 15. I presume in the failure the region with the mode is not discovered. Another way to fix the problem, at least in this case, is to specify order=30 to the integral with the absolute tolerance. The documentation says “It is advisable to increase the integration order in rough proportion to the precision, for smooth integrands.” I can’t turn that into an operational rule, beyond “increase the order if you increase precision.” I’d guess “precision” refers to the absolute value of the power of 10 for the desired accuracy, and that the proportion is relative to the default of 7 for the order and a tolerance of … sqrt(machine epsilon)? All of these approaches seem very brute force compared to starting with a line search for the mode. I understand quadgk can’t assume anything about the function it is integrating, but in my case they seem to be fairly smooth, unimodal, functions that would lend themselves well to a line search. It is conceivable there might be 2 humps, however. The presence of most of the mass far out makes Gauss-Hermite quadrature unappealing despite the presence of a normal density in the function being evaluated. Why are you setting your absolute error tolerance so high? If I change the tol1 definition to, e.g., tol1 = floatmin(Float64) * 1e50, I get: julia> trouble(100) Denominator = 1.7081181337842018e-5 ± 9.259586775755586e-14. Numerator = 0.00037171497626672603 ± 2.3738119338954454e-19 with atol = 2.2250738585072016e-258. Numerator = 0.0003717149762667195 ± 3.760194450016882e-12 with no atol.  No, precision refers to the arithmetic precision, e.g. if you are using BigFloat numbers. (If you want e.g. 100 digits of accuracy, you’re not going to achieve it in a reasonable amount of time without exploiting exponential convergence via high-order quadrature. See QuadGK examples: Arbitrary-precision integrals.) 1 Like Maybe there is a sharp peak somewhere, and sometimes it “misses” it? Some things to try: • If the problem is resolved by setting atol=0, so that it is just using rtol, then maybe you just need to use a smaller atol. e.g. atol = eps(r0) * 10 (assuming your denominator r0 sets the scale of the accuracy you are interested in). Or you could get a better estimate of the achievable absolute tolerance by computing e.g. the L1 or L2 norm of the numerator integrand (or by doing some analytical estimate). • If the problem is that there are sharp localized features and that the initial quadrature points just “miss” the localized peak, then one solution is indeed to increase the order parameter to sample the domain more finely. (Or, if you know where the sharp peak is, you can put an additional endpoint there. Or do something even more clever if you know something analytically about your function.) • If your exponentially decay a Gaussian exp(-(z/L)^2), then I would advise changing variables to u=z/L so that the decay lengthscale is of order 1. See improper integrals in the QuadGK manual. 3 Likes That absolute tolerance will not actually be achievable — it will essentially be using rtol instead (which defaults to sqrt(eps())). Cancellation error means that you can’t usually calculate \int f(x) to an absolute tolerance better than \epsilon \int |f(x)|. In the OP’s case, it seems like the denominator sets the scale of the errors that they are interested in, in which case an atol proportional to the denominator makes sense? 4 Likes Specifically the difference here is that you want eps(denom) rather than sqrt(eps(denom)) because you only want the atol to kick in for points small enough that it won’t affect the result. 2 Likes Even with tol1=eps(r0) things don’t work: Numerator = 7.122927197595122e-22 ± 7.494964918478972e-22 with atol = 3.3881317890172014e-21. So precision is same as significant digits? Isn’t that more what the relative tolerance is about? That’s a very interesting point. In this case the weighting adjustment turns the term in dens2 to e^{-0.02z^2}\approx e^{-(z/7.07)^2}, which is one reason the modes are so far out. On top of that the random variable is multiplied by 0.25 in conddens(), which may or may not be important for integral–it’s certainly another reason the mode is so far out. The extraordinarily improbable outcome requires a success probability of almost exactly 1, which requires a very large \eta, which requires a very large z since z is divided by 4 before contributing to \eta. With double precision (Float64), you can’t get more than about 15 digits, which the default order=7 can usually attain in a reasonable number of evaluations for smooth functions. It’s when you try to go to 100 or 1000 digits of accuracy with BigFloat that you may have no choice but to go to high-order schemes. PS. Your very specific argument type declarations like function dens2(z::Float64) and function trouble(n::Int) don’t make Julia any faster. See Argument-Type Declarations in the manual. @Ross_Boylan If you are willing to sacrifice the performance of QuadGK for the numerical guarantees Arblib might help you: using QuadGK using Arblib "conditional probability of n successes out of n tries" function conddens(n::Int, z) # changed to integers not to pollute the results with "beyond 16 digits garbage" η = -2 + z / 4 return (1 + exp(-η))^(-n) end "combined normal density and weight function, ignoring constants" function dens2(z, λ = oftype(z, 48) / 100) # likewise return exp((2 * λ - 1) * z^2 / 2) end function trouble(n::Int, radius) f0(z) = conddens(n, z) * dens2(z) r0, r0err = quadgk(f0, -radius, radius; atol = 0.0) f1(z) = z * conddens(n, z) * dens2(z) #same as z*f0(z) tol1 = r0 * sqrt(eps()) r1, r1err = quadgk(f1, -radius, radius; atol = tol1) r1b, r1berr = quadgk(f1, -radius, radius; atol = 0.0) println("Denominator =$(r0) ± $(r0err).") println("Numerator =$(r1) ± $(r1err) with atol =$(tol1).")
println("Numerator = $(r1b) ±$(r1berr) with no atol.")
return (r0, r0err), (r1, r1err), (r1b, r1err)
end

λ = Arb(48; prec = p) / 100
f0(z) = conddens(n, z) * dens2(z, λ)
f1(z) = z * conddens(n, z) * dens2(z, λ)

@info "" Denominator = r0 Numerator = r1
return (real(r0), real(r1))
end

r0, r1, r2 = trouble(100, Inf)
a0, a1 = trouble(100, Arb(1000, prec=64))


Caveats:

• Arblib does not support improper integrals (hence: radius). The total integral beyond (-1000, 1000) should fit in the eps i.e. don’t change the result though.

I guess that Arblib also supports interval arithmetic, which might be helpful sometimes.

1 Like

yes it does, but with BigFloats it doesn’t seem to be competitive with Arblib.jl:

julia> using QuadGK, Arblib

julia> f(x) = sin(x + exp(x))
f (generic function with 1 method)

julia> let y = 10.0, prec=128
setprecision(prec)
@time a = quadgk(f, big(0.0), big(y))
@time b = real(Arblib.integrate(f, Arb(0, prec=prec), Arb(y, prec=prec)))
@info "" a b
m,r = a
@assert all((m-r ≤ b, b ≤ m+r))
end
4.784665 seconds (32.15 M allocations: 979.900 MiB, 3.84% gc time)
0.120403 seconds (185.12 k allocations: 11.659 MiB, 74.81% compilation time)
┌ Info:
│   a = (0.3473008379090081427883270393701828570867, 2.661171977296800224531813650930867080159e-20)
└   b = [0.347300837909008142788327039370183 +/- 8.85e-34]

julia> let y = 10.0, prec=128
setprecision(prec)
@time a = quadgk(f, big(0.0), big(y))
@time b = real(Arblib.integrate(f, Arb(0, prec=prec), Arb(y, prec=prec)))
@info "" a b
m,r = a
@assert all((m-r ≤ b, b ≤ m+r))
end
3.932064 seconds (28.29 M allocations: 817.716 MiB, 4.95% gc time)
0.026477 seconds (131.86 k allocations: 8.338 MiB)
┌ Info:
│   a = (0.3473008379090081427883270393701828570867, 2.661171977296800224531813650930867080159e-20)
└   b = [0.347300837909008142788327039370183 +/- 8.85e-34]


Also Arblib strives to provide mathematically correct answers by tracking numerical errors (i.e. uses ball arithmethic) in every function, and using only provably correct algorithms (which doesn’t mean that the software is proved though…). On the other hand quadgk uses IEEE 754 arithmetic and therefore can’t provide such guarantees. This means that even if (v, e) = quadgk(f, a, b) it is not true that v-e \leqslant \int_{a}^{b}f(x)dx \leqslant v+e must hold mathematically.

@stevengj please correct me if what I claim about quadgk is wrong.

1 Like

quadgk uses whatever numeric type you give it, not necessarily IEEE 754. (It can compute quadrature nodes and weights in arbitrary precision.) It’s true that its error estimate is only a rough estimate (usually an upper bound for analytic functions), but that’s due mainly to truncation error, not to floating-point error. That being said, Arblib is definitely (a) more optimized for arbitrary-precision integration and (b) designed for rigorous error bounds. On the other hand, this kind of arithmetic is often an unaffordable luxury in high-performance applications.

Thanks for the clarification. I’m not sure what do you mean by truncation error, do you mean underflow/overflow with fixed precision input? E.g. in the other thread you mentioned this example:

julia> f(y) = (abs(y) < 10 ? expm1(y) * exp(-2.5y) : (exp(-1.5y) - exp(-2.5y))) / y^1.2
f (generic function with 1 method)

(0.6790524782732376, 9.565804913317641e-9)


but I noticed that depending on how I break the (0, Inf) interval I get vastly different results:

julia> let R = 2^16
end
[ Info: (2.4577749995638187e-8, 3.2622514628063114e-16)
[ Info: (0.6790524563328935, 4.826450691855635e-9)
[ Info: (0.0, 0.0)

julia> let R = 2^17
end
[ Info: (8.107633888319233e-9, 1.076141657259337e-16)
[ Info: (0.0, 0.0)
[ Info: (0.0, 0.0)


Is this a manifestation of truncation error that you mentioned?

that’s why I said that Arblib can be an option if you’re willing to sacrifice speed for correctness If you can’t and the numerics align, you’re absolutely better of with quadgk!

No, I mean the error that comes from sampling the integrand at only a finite number of points.

2 Likes

The function being integrated has a sharp peak that is “far” out, e.g., 20. So I looked into using LineSearches or Optim to find that point, figuring I would use it as an interval boundary for quadgk to get it to pay attention to that region. Many of those methods were also misled by the function (the derivative at 0 is ~10^-92), but some worked, expensively.

However, when I investigated adding a point to the center of the integration interval I discovered it didn’t seem to matter where the point was; I always got the right answer. That was great until I looked at the number of evaluations, which always seems maximal when I specify a central point. Does anyone know why that’s happening?

function trouble(n::Int)
f0(z) = conddens(n, z)*dens2(z)
r0, r0err, r0n = quadgk_count(f0, -Inf, Inf, atol=0.0)

f1(z) = z*conddens(n, z)*dens2(z)
tol1 = eps(r0)
r1, r1err, r1n = quadgk_count(f1, -Inf, Inf, atol=tol1)
r1b, r1berr, r1bn = quadgk_count(f1, -Inf, Inf, atol=0.0)

println("Denominator = $(r0) ±$(r0err) in $(r0n) iterations.") println("Numerator =$(r1) ± $(r1err) in$(r1n) iterations with atol = $(tol1).") println("Numerator =$(r1b) ± $(r1berr) in$(r1bn) iterations with no atol.")

println("\nIntegrals with additional central point x, n iterations and atol = $(tol1)") println(" x n estimate error") for x in -10.0:2.0:28.0 r1, r1err, r1n = quadgk_count(f1, -Inf, x, Inf; atol=tol1) println("$(x): $(r1n) :$(r1) ± \$(r1err)")
end
end


Yields

julia> trouble(100)
Denominator = 1.7081181337842018e-5 ± 9.259586775755586e-14 in 405 iterations.
Numerator = 7.122927197595122e-22 ± 7.494964918478972e-22 in 15 iterations with atol = 3.3881317890172014e-21.
Numerator = 0.0003717149762667195 ± 3.760194450016882e-12 in 375 iterations with no atol.

Integrals with additional central point x, n iterations and atol = 3.3881317890172014e-21
x   n   estimate   error
-10.0: 10000020 : 0.00037171497626672635 ± 2.9595962555413576e-19
-8.0: 10000020 : 0.00037171497626672684 ± 2.1600719988808766e-19
-6.0: 10000020 : 0.0003717149762667157 ± 2.9726608171516286e-19
-4.0: 10000020 : 0.0003717149762667145 ± 2.4210675385061672e-19
-2.0: 10000020 : 0.0003717149762667205 ± 3.083896044957896e-19
0.0: 10000020 : 0.0003717149762667261 ± 2.3738077979923826e-19
2.0: 10000020 : 0.0003717149762667186 ± 2.9872191578574917e-19
4.0: 10000020 : 0.00037171497626671866 ± 2.7164869127828906e-19
6.0: 10000020 : 0.00037171497626672245 ± 2.550440724768667e-19
8.0: 10000020 : 0.0003717149762667208 ± 2.3541032093695264e-19
10.0: 10000020 : 0.00037171497626672 ± 2.895442867024353e-19
12.0: 10000020 : 0.0003717149762667288 ± 3.0422940735211036e-19
14.0: 10000020 : 0.0003717149762667211 ± 3.081018957170925e-19
16.0: 10000020 : 0.0003717149762667235 ± 2.866930863194379e-19
18.0: 10000020 : 0.00037171497626672457 ± 2.991667677761489e-19
20.0: 10000020 : 0.00037171497626671893 ± 2.653119664175651e-19
22.0: 10000020 : 0.00037171497626671524 ± 2.792125425445524e-19
24.0: 10000020 : 0.00037171497626671416 ± 2.761971806481448e-19
26.0: 10000020 : 0.0003717149762667148 ± 3.023613886361402e-19
28.0: 10000020 : 0.0003717149762667188 ± 2.501072137376253e-19


Evaluating each subinterval separately showed that the problem is not from the specification of 3 points; one side or the other, and sometimes both, hit the maximum number of evaluations. It seemed whichever side included the central hump got lots of evaluations.

julia> trouble(100)
Denominator = 1.7081181337842018e-5 ± 9.259586775755586e-14 in 405 iterations.
Numerator = 7.122927197595122e-22 ± 7.494964918478972e-22 in 15 iterations with atol = 3.3881317890172014e-21.
Numerator = 0.0003717149762667195 ± 3.760194450016882e-12 in 375 iterations with no atol.

Integrals with additional central point x, n iterations and atol = 3.3881317890172014e-21
x   nLeft/nRight   left estimate/right   left error/ right error
-10.0: 15/10000005 : -6.616613168484678e-198/0.00037171497626671394 ± 1.9310739523417615e-201/1.4156288842082823e-19
-8.0: 15/10000005 : -2.8025193660849574e-176/0.0003717149762667112 ± 8.188815068545353e-180/1.5280888278869534e-19
-6.0: 15/10000005 : -6.079081756737245e-155/0.00037171497626671725 ± 1.7283075216732401e-158/1.537757561270762e-19
-4.0: 15/10000005 : -4.888003583802134e-134/0.000371714976266727 ± 1.2536466214988912e-137/1.5186588665236096e-19
-2.0: 15/10000005 : -8.132367412857363e-114/0.00037171497626671985 ± 1.4713447141469133e-117/1.4370218323664223e-19
0.0: 15/10000005 : -8.734072958846782e-96/0.0003717149762667247 ± 1.6222416597199295e-98/1.3196442938555292e-19
2.0: 15/10000005 : 1.133257247403499e-75/0.0003717149762667249 ± 2.1873453751257476e-80/1.662689468159731e-19
4.0: 15/10000005 : 1.4559276623308361e-58/0.00037171497626671513 ± 1.2396026409848553e-62/1.4361761334472752e-19
6.0: 15/10000005 : 9.313314424970754e-44/0.0003717149762667102 ± 1.2061183561148093e-49/1.2458181602759214e-19
8.0: 15/10000005 : 1.4111997813440777e-31/0.00037171497626672007 ± 5.115226963974558e-36/1.3921970551775055e-19
10.0: 15/10000005 : 3.749925720016767e-22/0.00037171497626671925 ± 1.9768566294976234e-26/1.510881367833836e-19
12.0: 45/10000005 : 2.5666846591536296e-15/0.0003717149762641551 ± 1.3787349706002912e-23/1.2141714046415605e-19
14.0: 135/10000005 : 1.1570851658940074e-10/0.00037171486055820693 ± 8.814646622734698e-22/1.0297856421686485e-19
16.0: 255/10000005 : 1.1050356724945566e-7/0.0003716044726994736 ± 4.313835242770197e-22/1.233645584387334e-19
18.0: 16515/10000005 : 6.762913356334931e-6/0.00036495206291038744 ± 3.387255740603334e-21/1.3574159688758056e-19
20.0: 10000005/10000005 : 6.542834531431433e-5/0.0003062866309523992 ± 2.597688132474256e-20/1.1245063331650594e-19
22.0: 10000005/10000005 : 0.00019820756381742585/0.00017350741244929097 ± 7.283721011603446e-20/5.959921023105253e-20
24.0: 10000005/10000005 : 0.00031082759737572133/6.088737889099611e-5 ± 9.131485170758274e-20/2.1357474022209636e-20
26.0: 10000005/10000005 : 0.0003577959418693771/1.3919034397343316e-5 ± 1.1395477966649482e-19/4.4500672742009865e-21
28.0: 10000005/495 : 0.00036947231384273834/2.242662423977671e-6 ± 1.1831071251315108e-19/3.0311160292058338e-21


Yes, because that is the side that will need to be refined the most.

To reduce the number of evaluations further, you could play with various subdivisions and/or increase the order of the quadrature rule. To do much better, you’d ideally need to know something analytically about the peak, as in this example.

The huge number of evaluations arose from setting the tolerance essentially to 0: tol1 = eps(r0), with the value used for an absolute tolerance. That was a holdover from some of the earlier experiments. tol1 = abs(r0)*√eps() worked much better, resulting in ~450 function evaluations. The number of evaluations didn’t seem terribly sensitive to the choice of cutpoint, as long as there was a cutpoint, and it didn’t vary in any consistent way, e.g., more or less repetitions as the cutpoint moved toward the actual mode.

This still seems like a lot of evaluations, albeit ~10^5 times better than before, but even integrating e^{-x^2} took around 220 evaluations. The latter expression is what you get automatically with Gauss-Hermite quadrature, which ordinarily does pretty well with 11 evaluations.

Yes, because integrating with infinite limits in QuadGK involves a singular coordinate transformation, which slows down convergence. In contrast, with Gauss–Hermite quadrature the singularity is built-in to the quadrature scheme, and hence the infinite bounds are “free” (and it integrates e^{-x^2} times low-degree polynomials exactly with only a few points) … at the price of not being adaptive.