Double integration with hcubature slower than nested quadgk

General Usage
1

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.
First the “bad way” by using nested quadgk calls

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.

Thanks in advance
Olivier

1 Like
2

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
3

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)
    sol  = adaptsim(ff, 0.0, 100.0)
    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.

4

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
5

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)
    sol= adaptsim(f1, 0.0, 100.0)
    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
    adaptsim(f1, -ran, 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)
    sol  = adaptsim(ff, 0.0, 100.0)
    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?

6

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
7

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?

Matlab document makes argument for advantages and disadvantages for different quad algorithm (maybe not true):
The list below contains information to help you determine which quadrature function in MATLAB to use:

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.

The quadv function vectorizes quad for an array-valued fun.

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.
quad

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.

8

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

1 Like
9

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