# Double integration with hcubature slower than nested quadgk

It is probably related to this topic this topic.

I find the performances of HCubature.jl not intuitive.

Lets take the following function as a test case for our MWE:

``````internal(k,x) = exp(-x^2)*sin(k*x)^2
``````

I tried performing the integration in 2 ways.

``````using QuadGK

function integrand(k)
f1(x) = internal(k,x)
(sol, er)  = quadgk(f1, -Inf, Inf)
return sol
end

function final(rtol)
(sol, err) = quadgk(integrand, 0.0, 100, rtol = rtol)
return sol, err
end

``````

And the “right way” using the multidimensional integration from HCubature.jl

``````function final2(rtol)
f2(k,t) = internal(k,t/(1-t^2))*(1+t^2)/(1-t^2)^2
f3(p) = f2(p[1],p[2])
(sol, err) = hcubature(f3, (0.0, -1), (100,1), rtol=rtol)
return sol,err
end

``````
``````julia> using Benchmarktools
julia> @time final(1e-9)
0.014318 seconds (789 allocations: 1.185 MiB)
(87.83729389116553, 8.133626228090662e-8)

julia> @time final2(1e-9)
12.298572 seconds (166.92 M allocations: 5.849 GiB, 9.48% gc time)
(87.83729438189894, 8.78372122849345e-8)
``````

I guess the result from final2() is more accurate than from final() but still I find the difference in time and allocations huge between both approaches.
It’s not clear to me how we can compare the accuracy for both methods.

Olivier

1 Like

Yes, I see that as well. If you increment a counter variable in `internal` (the innermost integrand), you’ll see it is evaluated 383745 for `quadgk` and `77742309` for `hcubature`.

You have to be careful about accuracy here because in your nested `quadgk` calls you are not setting a tolerance in the inner integrand, so it is actually lower than the tolerance in your outer integrand. I would tend to use a 10x smaller tolerance in the inner integrand. If I do this, I get `441285` evaluations for `final(1e-9)`, i.e. nested `quadgk`.

To assess the accuracy, I simply ran `final(1e-12)` (with a tolerance of `1e-13` in the inner integral). Compared to that, `final(1e-9)` has a relative error of `3.1e-12`, and `final2(1e-9)` has a relative error of `5e-9`.

I suspect there are three things going on here:

1. `quadgk` is defaulting to a higher-order quadrature rule than `hcubature`, so it can converge faster for smooth integrands. (Here, it isn’t smooth at the `x` endpoints because of the singularity in the coordinate transform, but it is smooth elsewhere.)

2. The main advantage of cubature over nested 1d quadrature arises in cases where your function is badly behaved in a small portion of the domain, so that choosing the grids one dimension at a time is suboptimal. That isn’t really the case for your test integrand here.

3. IIRC `quadgk` is currently more optimized than `hcubature` for minimizing allocations, which matters a lot for cheap integrands; I think we can probably cut down on the latter.

4 Likes

Hi Oliver. It seems that `QuadGK` is well optimized an even supports intergration over infinity. Maybe in most of problems, using nested `QuadGK` is better choice.
And I tried to use `QuadGG` (adaptive Simpson method) to repeat your result however it doesn’t work. Maybe `QuadGG` doesn’t support infinity intergration. Here are my attempts:

``````using BenchmarkTools,HCubature,QuadGG,QuadGK
internal(k,x) = exp(-x^2)*sin(k*x)^2
function integrand3(t)
f(k,t) = internal(k,t/(1-t^2))*(1+t^2)/(1-t^2)^2
ff(k)=f(k,t)
return sol
end

function final3(rtol)
sol = adaptsim(integrand3, -1.0, 1.0; tol = rtol)
return sol
end

function integrand4(k)
f1(x) = internal(k,x)
(sol, er)  = adaptsim(f1, -Inf, Inf)
return sol
end

function final4(rtol)
sol = adaptsim(integrand4, -1.0, 1.0; tol = rtol)
return sol
end

@btime final3(1e-9)
DomainError with -Inf:
sin(x) is only defined for finite x.

@btime final4(1e-9)
DomainError with Inf:
sin(x) is only defined for finite x.
``````

The variable changing trick fails here.

The basic problem is that you have to be more careful with the coordinate-transform trick (which is what `quadgk` uses internally for infinite integration domains) with the methods from QuadGG, because QuadGG’s algorithms evaluate the integrand at the endpoints (where your coordinate transform gives 0/0 integrands). You would need to take the limit to get the correct endpoint values.

In contrast, `quadgk` and `hcubature` only evaluate the integrand in the interior of the domain, so they can handle (integrable) singularities at the endpoints without any special handling.

3 Likes

Thank for your great suggestion. I try to combine package `Richardson` with `QuadGG` to compute abnormal intergration, however there are still some unexpected errors.

``````using BenchmarkTools,QuadGG,Richardson

function integrand5(x)
f1(k) = internal(k,x)
return sol
end

function intergration(rtol,ran)
val = adaptsim(integrand5, -ran,ran;tol = rtol)
return val
end

function final5(rtol)
sol=extrapolate(100000.0, x0=Inf) do x
@show x
intergration(rtol,x)
end
return sol
end

function integrand6(k)
f1(x) = internal(k,x)
(sol, er)=extrapolate(100000.0, x0=Inf) do ran
@show ran
end
return sol
end

function final6(rtol)
sol = adaptsim(integrand4, -1.0, 1.0; tol = rtol)
return sol
end

function integrand7(t)
f(k,t) = internal(k,t/(1-t^2))*(1+t^2)/(1-t^2)^2
ff(k)=f(k,t)
return sol
end

function intergration7(rtol,ran)
val = adaptsim(integrand7, -ran,ran;tol = rtol)
return val
end

function final7(rtol)
(sol,err)=extrapolate(0.5, x0=1.0) do x
@show x
intergration7(rtol,x)
end
return (sol,err)
end

@btime final5(1e-9)
x = 99999.99999999999
x = 799999.9999999999
x = 99999.99999999999
x = 799999.9999999999
x = 99999.99999999999
x = 799999.9999999999
(0.0, 0.0)

@btime final6(1e-9)
DomainError with Inf:
sin(x) is only defined for finite x.

@btime final7(1e-9)
x = 1.5
x = 1.0625
x = 1.0078125
x = 1.0009765625
x = 1.0001220703125
x = 1.0000152587890625
x = 1.0000019073486328
x = 1.000000238418579
x = 1.0000000298023224
x = 1.5
x = 1.0625
``````

If I intergrate finite part first then fun just always oscillate between 99999 and 79999 and if I reverse order there is still Inf error. Moreover if I change coordinate, the intergration never converge but `x` will be in a loop from 1.5 to 1.0000000298023224, let alone intergrating fastly. Do you have some idea to make `QuadGG` fit to intergration with infinity?

I wouldn’t do this, I would just analytically plug in zero for the coordinate-transformed integrand at the endpoints, since that is the limit of the integrand.

However, I don’t quite see the point. I don’t see what advantage there is to Gauss–Lobatto or Simpson’s rule (quadgg) over Gauss-Kronrod (quadgk).

If you really want to do better, you would use a quadrature rule specialized for the particular asymptotic form of your integrand, e.g. Gauss–Hermite quadrature for e^{-x^2} decay as in this example. FastGaussQuadrature.jl provides this, but it’s not adaptive so you will have to pick a particular order yourself to achieve a desired tolerance.

2 Likes

But in this example for `QuadGK`, there isn’t any coordinate transform. And it can directly integrate over an infinite domain. Did you set auto transform in source codes?

The quad function may be most efficient for low accuracies with nonsmooth integrands.

The quadl function may be more efficient than quad at higher accuracies with smooth integrands.

The quadgk function may be most efficient for high accuracies and oscillatory integrands. It supports infinite intervals and can handle moderate singularities at the endpoints. It also supports contour integration along piecewise linear paths.

If the interval is infinite, [a,∞), then for the integral of fun(x) to exist, fun(x) must decay as x approaches infinity, and quadgk requires it to decay rapidly. Special methods should be used for oscillatory functions on infinite intervals, but quadgk can be used if fun(x) decays fast enough.

The quadgk function will integrate functions that are singular at finite endpoints if the singularities are not too strong. For example, it will integrate functions that behave at an endpoint c like log|x-c| or |x-c|p for p >= -1/2. If the function is singular at points inside (a,b), write the integral as a sum of integrals over subintervals with the singular points as endpoints, compute them with quadgk, and add the results.

I just feel confused that the support of infinite integration by `QuadGK` is because of the Gauss-Kronrod algorithm (unnecessary to evaluate on endpoints) or you make some special treatment?
Hope my naive questions are not disturbing you.

Internally it does the same coordinate transform to map to a finite domain, and then applies Gauss-Kronrod quadrature.

1 Like

I’m so sorry for not reacting earlier to your very detailed response.

Thank you very much again. Very instructive as always.

I wanted to come up with a better function to showcase the difference in behavior, but I haven’t found the time yet.

I hope I’ll be able to give some follow up on this soon.

Thanks again,

Olivier