# How to improve the implementation of a function involving a numerical integration?

**URL:** <https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444>\
**Category:** Performance\
**Tags:** quadgk, integral\
**Created:** [October 26, 2023, 3:37pm UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444 "2023-10-26T15:37:24Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![lsablon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lsablon/32/47083_2.png) [@lsablon](https://discourse.julialang.org/u/lsablon)\
**Post date:** [October 26, 2023, 3:37pm UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/1 "2023-10-26T15:37:24Z")

</div>

Hello,

I am trying to optimize my code, but I don’t really know in which direction I should go (but I know I don’t want to use multithreading).

Here is my goal: I want to implement the following function

k\_\perp(x,y; l) = k(x,y;l) - \frac{\int\_{-1}^{1} p(w)k(x,z;l)dw \times \int\_{-1}^{1} p(z)k(z,y;l)dz}{ \int\_{-1}^{1} \int\_{-1}^{1} p(w)p(z)k(w,z;l)dwdz}

where p(x) = \frac{1}{\pi} \sqrt{1-x^2} and k(x,y;l) = \exp{-\frac{(x-y)^2}{2l^2}} .  
(this function will later be used to build a kernel and some optimization will be performed on l)

Unfortunately, this integral has no analytical solution (at least according to Mathematica).

Therefore, I wrote this little code (see below) using QuadGK to compute the function, but I am not really happy regarding the performance and I am pretty sure it could be faster. I also saw that ApproxFun could probably do something interesting using some orthogonal polynomials but did not manage to write a clean code…

```julia
using QuadGK, BenchmarkTools
k(x,y,l) = exp(-(x-y)^2/(2*l^2))
p(x) = sqrt(1-x^2)/π

ES(x,l) = quadgk(y -> k(x,y,l) * p(y), -1, 1)[1]
@benchmark ES(0.0, 1e1)

ES2(l) = quadgk(x -> p(x) * ES(x,l), -1, 1)[1]
@benchmark ES2(1e1)

k_orth(x,y,l) = @. k(x,y,l) - ES(x,l) * ES(y,l) / ES2(l)

x = range(-1.0, 1.0, length=100);
y = range(-1.0, 1.0, length=100);

@benchmark k_orth(x, y, 1e1)

contourf(x, y, (x,y) -> k_orth(x, y, 1e1))

```

Thanks,  
L.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [October 26, 2023, 3:48pm UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/2 "2023-10-26T15:48:33Z")

</div>

> [@lsablon](#):
>
> I am pretty sure it could be faster

why?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 26, 2023, 4:10pm UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/3 "2023-10-26T16:10:45Z")

</div>

One thing worth trying is using a package for multi-dimensional integration (e.g. [GitHub - JuliaMath/HCubature.jl: pure-Julia multidimensional h-adaptive integration](https://github.com/JuliaMath/HCubature.jl)) rather than nesting 1d integration. Doing so can be a lot more efficient since you get better error control.

---

<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:** [October 26, 2023, 4:54pm UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/4 "2023-10-26T16:54:50Z")

</div>

> [@lsablon](#):
>
> `ES(x,l) = quadgk(y -> k(x,y,l) * p(y), -1, 1)[1]`

Note that you are using the default tolerances here — by default, `quadgk` tries to get at least 8 significant digits. Do you need that kind of accuracy?

What is the typical value of the Gaussian width, l? If it is small (i.e. the function is sharply peaked) then you might want to do something more specialized.

Another basic issue here is that you have a square-root singularity at the edges of your domain, via the p(y) factors. This will slow down convergence substantially unless you employ a specialized quadrature scheme. (I’d eventually like you to be able to specify such endpoint singularities analytically to `quadgk`.) @lxvm has an experimental package with some support for this called [QuadJGK.jl](https://github.com/lxvm/QuadJGK.jl).

Even simpler, do a change of variables y=\cos\theta that eliminates the singularity:

\int\_{-1}^{+1} f(y) p(y) dy = \frac{1}{\pi} \int\_0^\pi f(\cos\theta) \sin(\theta)^2 d\theta \, ,

since p(\cos\theta) = \sin(\theta)/\pi for \theta \in [0,\pi].

For example, with x=0.1 and l = 10, this change of variables gives the same answer (actually, an even more accurate answer) with 15x fewer function evaluations:

```julia
julia> quadgk_count(y -> k(0.1,y,1e1) * p(y), -1, 1)
(0.49935087498090297, 6.368550962729288e-9, 645)

julia> quadgk_count(θ -> k(0.1,cos(θ),1e1) * sin(θ)^2/π, 0,pi)
(0.49935087461297634, 5.176886697100258e-12, 45)

```

If you lower the tolerance to `1e-4` then you can get another factor of 3 on this example (and it actually is still accurate to machine precision because the integrand is so smooth!):

```julia
julia> quadgk_count(θ -> k(0.1,cos(θ),1e1) * sin(θ)^2/π, 0,pi, rtol=1e-4)
(0.4993508746129765, 6.152897258360213e-8, 15)

```

You can also try [pre-allocating buffers](https://juliamath.github.io/QuadGK.jl/stable/api/#QuadGK.alloc_segbuf) to speed things up, but I would focus first on giving QuadGK (or any integration routine!) “nicer” functions to integrate.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 26, 2023, 11:23pm UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/5 "2023-10-26T23:23:10Z")

</div>

Wouldn’t a change of variables to use [erf](https://en.wikipedia.org/wiki/Error_function) be the way to go for `ES`?

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [October 27, 2023, 12:25am UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/6 "2023-10-27T00:25:44Z")

</div>

> [@stevengj](#):
>
> ```julia
> julia> quadgk_count(θ -> k(0.1,cos(θ),1e1) * sin(θ)^2/π, 0,pi, rtol=1e-4)
> (0.4993508746129765, 6.152897258360213e-8, 15)
> 
> ```

The change of variables trick works about as well as the specialized quadrature for Jacobi weight functions with QuadJGK.jl. I get:

```julia
julia> using QuadJGK
[Info: Precompiling QuadJGK [f27e79a0-bbaa-47b5-b2d3-8db3cef8b089]

julia> k(x,y,l) = exp(-(x-y)^2/(2*l^2))
k (generic function with 1 method)

julia> p(x) = sqrt(1-x^2)/π
p (generic function with 1 method)

julia> quadjgk_count(y -> k(0.1,y,1e1) * p(y), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)))
(0.49935087461297667, 1.1102230246251565e-16, 15)

```

For this example, machine precision is still possible with 7 integrand evaluations

```julia
julia> quadjgk_count(y -> k(0.1,y,1e1) * p(y), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)), order=3)
(0.4993508746129767, 1.6239870559431324e-10, 7)

```

And the same is true for the second integral

```julia
julia> ES(x,l) = quadjgk(y -> k(x,y,l) * p(y), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)))[1]
ES (generic function with 1 method)

julia> quadjgk_count(x -> p(x) * ES(x,1e1), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)), order=3)
(0.24937694744330427, 8.052006283953972e-11, 7)

```

So `7^2 = 49` function evaluations should be sufficient to compute the normalization for this large value of `l=1e1`. For \ell \ll 1, HCubature.jl is probably the better option for the 2d integral since the integrand evaluations scale as \mathcal{O}(\log(\ell)) as \ell \to 0 for cubatures of a Gaussian localized at the origin, whereas for iterated integration it will be \mathcal{O}(\log^2(\ell)) (i.e. one factor for each integral).

Update: by loosening the error tolerance, it seems like only 5 integrand evaluations are needed for each 1d integral, and the result is still accurate to machine precision

```julia
julia> quadjgk_count(y -> k(0.1,y,1e1) * p(y), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)), order=2, rtol=1e-4)
(0.4993508746129768, 3.8987719608307714e-7, 5)

julia> quadjgk_count(x -> p(x) * ES(x,1e1), JacobiSpace(0.5, 0.5, OpenInterval(-1, 1)), order=2, rtol=1e-4)
(0.24937694744330427, 1.9377299406708737e-7, 5)

```

Edit: I wasn’t thinking: k(x,y;l) is localized on the line x=y, not at a point, so the scaling of cubature is like \mathcal{O}(\log(l)/l), which makes iterated integration more efficient as l \to 0. It would be hard to get rid of the logarithmic scaling without doing something special

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [October 27, 2023, 2:52am UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/7 "2023-10-27T02:52:51Z")

</div>

You could precompute ES and ES2 and then simply do interpolation. ApproxFun is pretty good for that, but in reality any interpolation package can do the Job.

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [October 27, 2023, 4:10am UTC](https://discourse.julialang.org/t/how-to-improve-the-implementation-of-a-function-involving-a-numerical-integration/105444/8 "2023-10-27T04:10:03Z")

</div>

You can also do a lot of sneaky things given that you have a gaussian.  
You can identify 2 domains

First, you can chop a lot of the domain when l is small enough. In ES(x,l), when x \in (-1,1), solve for w in exp( -(x-w)^2/(2l^2) ) \< \epsilon\_{tol}, then change your integration limits to be (\max(-1, w\_min), \min(1, w\_max)).

When l is large you can use the FFT, since ES is actually a convolution between a gaussian (that we will denote by K(x-y, l)) and p (where you have to extend p by 0 outside of (-1,1)). Note that \int\_{-1}^1 exp(-i2 \pi \xi x) p(x) = \frac{J\_1(2\pi \xi )}{2\pi \xi }. So ES(x,l) = \int\_{-\infty}^\infty \frac{J\_1(2\pi \xi )}{2\pi \xi } \hat{K}(\xi, l) e^{2\pi i \xi x} d\xi. The bandwidth of \hat{K} is small, so you can use the same trick of choping of the integral to the desired tolerance.

For the second integral, if l is small, since everything is concentrated in the diagonal, probably you are best with higher dimensional techniques, since you don not want to spend that much time off the diagonal.

For the the case of large l, you can do the fourier transform 2 times, so you are only left with 1 integral instead of 2. If I got everything right, ES2(l) = \int\_{-\infty}^\infty \left(\frac{J\_1(2 \pi \xi)}{2 \pi \xi} \right)^2 \hat{K}(\xi, l) d\xi

Note: When doing the integration in fourier space, after doing the chop of, just use the trapezoidal rule, it will give you exponential convergence.