# Fast evaluation of multivariate integrals in Julia for integrands depending on many parameters

**URL:** https://discourse.julialang.org/t/fast-evaluation-of-multivariate-integrals-in-julia-for-integrands-depending-on-many-parameters/71815
**Category:** Performance
**Tags:** package, integral
**Created:** [November 20, 2021, 8:17am UTC](https://discourse.julialang.org/t/fast-evaluation-of-multivariate-integrals-in-julia-for-integrands-depending-on-many-parameters/71815 "2021-11-20T08:17:59Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![vicdolean](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vicdolean/32/30134_2.png) [@vicdolean](https://discourse.julialang.org/u/vicdolean)
#### Post date: [November 20, 2021, 8:17am UTC](https://discourse.julialang.org/t/fast-evaluation-of-multivariate-integrals-in-julia-for-integrands-depending-on-many-parameters/71815/1 "2021-11-20T08:17:59Z")

</div>

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!

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 20, 2021, 11:16am UTC](https://discourse.julialang.org/t/fast-evaluation-of-multivariate-integrals-in-julia-for-integrands-depending-on-many-parameters/71815/2 "2021-11-20T11:16:49Z")

</div>

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.

---

<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 20, 2021, 1:14pm UTC](https://discourse.julialang.org/t/fast-evaluation-of-multivariate-integrals-in-julia-for-integrands-depending-on-many-parameters/71815/3 "2021-11-20T13:14:31Z")

</div>

> [@vicdolean](#):
>
> two dimensional domain D that could be R^2, given than f(x) is a exponentially decreasing function

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](https://github.com/stevengj/cubature/blob/master/README.md#infinite-intervals).

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](https://github.com/JuliaMath/HCubature.jl) (which is also callable via Quadrature.jl).

---

<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 20, 2021, 1:28pm UTC](https://discourse.julialang.org/t/fast-evaluation-of-multivariate-integrals-in-julia-for-integrands-depending-on-many-parameters/71815/4 "2021-11-20T13:28:58Z")

</div>

(And, of course, you should read the [Julia performance tips](https://docs.julialang.org/en/v1/manual/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.)

---

<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 20, 2021, 3:02pm UTC](https://discourse.julialang.org/t/fast-evaluation-of-multivariate-integrals-in-julia-for-integrands-depending-on-many-parameters/71815/5 "2021-11-20T15:02:35Z")

</div>

> [@vicdolean](#):
>
> The integrand depends on 8 parameters, so this kind of integral needs to be computed O(10^6) times.

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.
