Hi,
I need to evaluate an oscillatory integral of the form
\int_D f(x) exp(\omega g(x)) dx
on a two dimensional domain D that could be R^2, given than f(x) is a exponentially decreasing function equal to an exponential time some smooth function and g(x) is simply an affine function.
The integrand depends on 8 parameters, so this kind of integral needs to be computed O(10^6) times. Matlab black box integral2 takes 0.1 s per integral → more than 24h in total !
Some of the integrals are very small and can be neglected for some parameter but I still need to evaluate an important number leading to prohibiting solution times.
Can someone recommend the fastest possible Julia routine that could provide a reasonable solution time?
Thanks a lot!
Use Quadrature.jl for the in-place form or static array form and try a few methods like CubatureHJL or cuhre and see what works out.
Note that you’ll not want to treat this by simply using a large domain, because then your integrand will be sharply peaked in a tiny portion of the domain.
One approach is to use a change of variables to transform the domain to a finite one.
Even better, you should analytically estimate the asymptotic decay rate of your integrand and exploit it if possible. At the simplest level, if you know that your integrand decays over a lengthscale L_i in variable x_i, you can rescale variables to x_i' = x_i / L_i (before the infinite → finite change above), so that the decay length in the new variables is always \sim 1. A more sophisticated approach would be to use a tensor product of something like Gauss–Legendre or Gauss–Hermite quadrature rules, but this requires you to know more about your integrand and about quadrature.
For example, if your integral \int F(x_1,x_2) is over \mathbb{R}^2 and you have characteristic decay lengths L_i in the two directions (determined analytically!), you could transform it to:
\int_{-\infty}^\infty dx_1 \int_{-\infty}^\infty F(x_1, x_2) dx_2 = L_1 L_2 \int_{-\infty}^\infty dx_1' \int_{-\infty}^\infty F(x_1' L_1, x_2' L_2) dx_2' \\
= L_1 L_2 \int_{-1}^1 dt_1 \int_{-1}^1F\left(\frac{t_1}{1-t_1^2} L_1, \frac{t_2}{1-t_2^2} L_2\right) \frac{(1+t_1^2) (1+t_2^2) }{ (1-t_1^2)^2 (1-t_2^2)^2 } dt_2
and then apply some adaptive 2d quadrature rule like hcubature from GitHub - JuliaMath/HCubature.jl: pure-Julia multidimensional h-adaptive integration (which is also callable via Quadrature.jl).
5 Likes
(And, of course, you should read the Julia performance tips to make sure your integrand is fast. Make sure your integrand doesn’t depend on global variables—use a closure to pass in parameters—and make sure it is type-stable and non-allocating.)
1 Like
The other thing, of course, is to re-examine why you need to evaluate this integral so many times. Are you optimizing over those 8 parameters, or trying to construct an interpolant? If so, there might be ways to sample it fewer times.
1 Like