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

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…

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