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

General Usage Performance
,
6

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

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

2 Likes