# Double integration with hcubature slower than nested quadgk

**URL:** <https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943>\
**Category:** General Usage\
**Created:** [November 8, 2022, 4:38pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943 "2022-11-08T16:38:19Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Olivier\_Merchiers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivier_merchiers/32/4073_2.png) [@Olivier\_Merchiers](https://discourse.julialang.org/u/Olivier_Merchiers)\
**Post date:** [November 8, 2022, 4:38pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/1 "2022-11-08T16:38:19Z")

</div>

It is probably related to this topic [this topic](https://discourse.julialang.org/t/unreasonable-slow-speed-in-numerical-integration/89910).

I find the performances of HCubature.jl not intuitive.

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

```julia
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

```julia
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

```julia
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
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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 8, 2022, 5:35pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/2 "2022-11-08T17:35:03Z")

</div>

> [@Olivier\_Merchiers](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![swish47](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/swish47/32/43921_2.png) [@swish47](https://discourse.julialang.org/u/swish47)\
**Post date:** [November 9, 2022, 7:06pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/3 "2022-11-09T19:06:51Z")

</div>

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 9, 2022, 7:20pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/4 "2022-11-09T19:20:05Z")

</div>

> [@swish47](#):
>
> Maybe `QuadGG` doesn’t support infinity intergration. Here are my attempts:

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.

---

<div class="post-metadata">

**Author:** ![swish47](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/swish47/32/43921_2.png) [@swish47](https://discourse.julialang.org/u/swish47)\
**Post date:** [November 9, 2022, 9:51pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/5 "2022-11-09T21:51:40Z")

</div>

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 9, 2022, 10:37pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/6 "2022-11-09T22:37:12Z")

</div>

> [@swish47](#):
>
> I try to combine package `Richardson` with `QuadGG` to compute abnormal intergration, however there are still some unexpected errors.

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](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) provides this, but it’s not adaptive so you will have to pick a particular order yourself to achieve a desired tolerance.

---

<div class="post-metadata">

**Author:** ![swish47](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/swish47/32/43921_2.png) [@swish47](https://discourse.julialang.org/u/swish47)\
**Post date:** [November 10, 2022, 4:50pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/7 "2022-11-10T16:50:08Z")

</div>

> [@stevengj](#):
>
> I would just analytically plug in zero for the coordinate-transformed integrand at the endpoints, since that is the limit of the _integrand_.

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?

> [@stevengj](#):
>
> However, I don’t quite see the point. I don’t see what advantage there is to Gauss–Lobatto or Simpson’s rule (quadgk) over Gauss-Kronrod (quadgk).

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](https://www.mathworks.com/help/releases/R2021a/matlab/ref/quad.html?s_tid=doc_ta)

> [@stevengj](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 10, 2022, 7:07pm UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/8 "2022-11-10T19:07:35Z")

</div>

> [@swish47](#):
>
> But in this example for `QuadGK`, there isn’t any coordinate transform. And it can directly integrate over an infinite domain

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

---

<div class="post-metadata">

**Author:** ![Olivier\_Merchiers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivier_merchiers/32/4073_2.png) [@Olivier\_Merchiers](https://discourse.julialang.org/u/Olivier_Merchiers)\
**Post date:** [January 19, 2023, 9:27am UTC](https://discourse.julialang.org/t/double-integration-with-hcubature-slower-than-nested-quadgk/89943/9 "2023-01-19T09:27:34Z")

</div>

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