# Multidimentional integrals in julia: speed and automatic differentiation

**URL:** <https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375>\
**Category:** Numerics\
**Created:** [June 14, 2020, 2:18pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375 "2020-06-14T14:18:44Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [June 14, 2020, 2:18pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/1 "2020-06-14T14:18:44Z")

</div>

I would like to differentiate through a two-dimensional integral.

- seems to be not possible for c++ wrappers (e.g. with the great [Cuba.jl](https://github.com/giordano/Cuba.jl))
- Julia native integrators are way slower ([HCubature.jl](https://github.com/JuliaMath/HCubature.jl), [NIntegration.jl](https://github.com/pabloferz/NIntegration.jl))
- nested [QuadGK.jl](https://github.com/JuliaMath/QuadGK.jl) is also not fast enough.

Therefore, two questions:

1. Are there tricks to differentiate `Cuba`?
2. Is there a way to speed up Julia native implementations?

Here is MWE and tests

```julia
using BenchmarkTools

using HCubature
using NIntegration
using QuadGK
using Cuba

function model(cosθ,ϕ; pars)
    α1 = α2 = pars[1]-1
    v = (exp(-α1*(cosθ+1)) + exp(-α2*(1-cosθ)))*(1+sin(ϕ))^(1.2)
    return abs2(v)
end
#
# using Plots
# plot(cosθ -> model(cosθ,π/4; pars = [2.3]), -1, 1)
#
I_quadgk(pars) = quadgk(cosθ->quadgk(ϕ->model(cosθ,ϕ;pars=pars),-π,π)[1],-1.0, 1.0)[1]
I_nintegration(pars) = nintegrate((cosθ,ϕ,z)->model(cosθ,ϕ;pars=pars), (-1.0, -π, 0.0), (1.0, π, 1.0))[1]
I_hcubature(pars) = hcubature(x->model(x[1],x[2]; pars= pars), [-1.0, -π], [1.0, π])[1]
I_cuba(pars) = cuhre((x,f)->f[1]=model(2x[1]-1, π*(2x[2]-1);pars=pars),2,1).integral[1]*(4π)
#
@btime I_cuba([2.3])
@btime I_quadgk([2.3])
@btime I_hcubature([2.3])
@btime I_nintegration([2.3])

```

Results are:

```julia
1) Cuba: 72.500 μs (1305 allocations: 61.28 KiB)
2) 2x QuadGK: 160.800 μs (827 allocations: 19.44 KiB)
3) NIntegration: 303.800 μs (209 allocations: 9.91 KiB)
4) HCubature: 2.215 ms (23787 allocations: 850.86 KiB)

```

I found it quite surprising with respect to the number of github stars on the packages 🙂

EDIT: replaced parameters to a vector.

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [June 14, 2020, 2:20pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/2 "2020-06-14T14:20:57Z")

</div>

Pin the experts: @stevengj, @giordano, @PabloZubieta.  
Any thoughts are welcome!

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [June 14, 2020, 2:28pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/3 "2020-06-14T14:28:06Z")

</div>

Shouldn’t it be fairly straightforward to just use an analytical derivative here? If you’re using this inside other code using e.g. Zygote, you could also write a custom adjoint for Cuba.jl.

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [June 14, 2020, 3:00pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/4 "2020-06-14T15:00:40Z")

</div>

My working problem is more complicated than that.  
However, trying `Zygote` is an interesting idea. Could maybe sketch an example?  
Let’s say, the parameter is a vector of dim2. Thanks in advance!

---

<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:** [June 14, 2020, 3:00pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/5 "2020-06-14T15:00:54Z")

</div>

> [@simeonschaub](#):
>
> Shouldn’t it be fairly straightforward to just use an analytical derivative here? If you’re using this inside other code using e.g. Zygote, you could also write a custom adjoint for Cuba.jl.

Yes. FWIW, we made Quadrature.jl parameter explicit to make it easy to write the adjoint, I just haven’t gotten to it yet. It would be fairly straightforward to do though.

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [June 14, 2020, 3:05pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/6 "2020-06-14T15:05:22Z")

</div>

With Zygote, you could solve it like this, by differentiating through the integrand:

```julia
using Cuba, Zygote

function model(cosθ,ϕ; αp1 = 2.3)
    α1 = α2 = αp1-1.0
    v = (exp(-α1*(cosθ+1)) + exp(-α2*(1-cosθ)))*(1+sin(ϕ))^(1.2)
    return abs2(v)
end

integrand(x, αp1) = model(2x[1]-1, π*(2x[2]-1);αp1=αp1)
I_cuba(αp1) = cuhre((x, f) -> f[1]=integrand(x, αp1),2,1).integral[1]*(4π)
Zygote.@adjoint I_cuba(αp1) = I_cuba(αp1), Δ -> (Δ * cuhre((x, f) -> f[1]=Zygote.pullback(a -> integrand(x, a), αp1)[2](1)[1], 2, 1).integral[1] * (4π),)

```

And then get the gradient:

```julia
julia> Zygote.gradient(x -> I_cuba(x[1]), [2.3])
([-13.368478696191737],)

```

It should be fairly easy to replace `αp1` with a vector if you have multiple parameters.

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [June 14, 2020, 3:08pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/7 "2020-06-14T15:08:38Z")

</div>

If you only a few parameters, it might be also worth considering using `ForwardDiff` for the integrand instead of `Zygote.pullback`, since that will then probably be faster.

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [June 14, 2020, 3:10pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/8 "2020-06-14T15:10:23Z")

</div>

Nice, I will try that. Thanks a lot

---

<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:** [June 14, 2020, 7:17pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/9 "2020-06-14T19:17:29Z")

</div>

Regarding the benchmark — it isn’t an entirely fair comparison, since you probably aren’t computing the integrals to the same accuracy in the different routines. For example, in NIntegration I think the default relative error is `1e-6`, whereas in HCubature the default is `sqrt(ε) ≈ 1.5e-8`

(Even if you set the same error tolerances in all the approaches, because they estimate errors in different ways you might get substantially different accuracy. One way to compare would be to compute the “exact” integral using a very low tolerance in one of the routines, and then check out what tolerance is needed in each routine to roughly obtain a given relative error, e.g. 1e-6.)

That being said, I suspect that there may have been a performance regression in HCubature — it didn’t use to require so many allocations if I recall correctly (it uses StaticArrays internally so its temporary arrays are not supposed to be heap-allocated.) I’ll have to look into it.

Since you have a smooth integrand in 2d, you might also want to try the `pcubature` function from Cubature.jl.

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [June 14, 2020, 8:16pm UTC](https://discourse.julialang.org/t/multidimentional-integrals-in-julia-speed-and-automatic-differentiation/41375/10 "2020-06-14T20:16:14Z")

</div>

A good point, I have not thought about accuracy, indeed.  
I tend to assume pragmatically that the returned value is exact, but of course it is not the case.

Thank you for suggestion trying `Cubature`, I was puzzled by the difference between `Cubature` and `HCubature`.