# Unreasonable slow speed in numerical integration

**URL:** https://discourse.julialang.org/t/unreasonable-slow-speed-in-numerical-integration/89910
**Category:** Performance
**Tags:** question, package, quadgk, integral
**Created:** [November 7, 2022, 10:17pm UTC](https://discourse.julialang.org/t/unreasonable-slow-speed-in-numerical-integration/89910 "2022-11-07T22:17:38Z")
**Posts on this page:** 1
**Showing post:** 2

<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 7, 2022, 11:35pm UTC](https://discourse.julialang.org/t/unreasonable-slow-speed-in-numerical-integration/89910/2 "2022-11-07T23:35:40Z")

</div>

> [@swish47](#):
>
> not comparable with Matlab and other programming language (one can check that Matlab only spend 1e-2 seconds on intergrating).

First, you should realize that the [Matlab `integrate` function](https://www.mathworks.com/help/matlab/ref/integral.html) defaults to `1e-6` relative tolerance (and `1e-10` absolute tolerance), whereas `quadgk` and `hcubature` default to `√ε ≈ 1.5e-8` relative tolerance, so you are asking for 100x times more accurate an answer from Julia than from Matlab, or even more if the integral is small (so that the absolute tolerance comes into play).

Second you are doing nested numerical integration, which is generally a bad idea or at least something one should be extremely cautious about:

1. Because your integrand is itself an adaptive numerical integral to `1e-8` relative tolerance, that effectively can mean it is “noisy” at the `1e-8` level. This can cause problems if you _also_ try to converge the “outer” integral to `1e-8` tolerance, because then it wastes a lot of function evaluations trying to converge the integral of the “noisy” _error_ in the interior integral. If you must use nested numerical integration, you should really be more careful of the tolerances (the outer integral should have a larger tolerance than the inner integral), and you may also need to set an absolute tolerance in order to avoid spending lots of time in regions where the integrand is small.

2. Instead of nesting numerical integrations, it is sometimes much better to do a single multidimensional integral. (Here, a single 3d integral.)

Third, you realize that the _first_ time you run Julia code, it spends a bunch of time compiling it, so that the _second_ (and subsequent) calls are fast, right? See the [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/#Measure-performance-with-%5B@time%5D(@ref)-and-pay-attention-to-memory-allocation) on measuring time. You also need to be careful about benchmarking in global scope, since Julia is oriented towards [performance-critical code in functions](https://docs.julialang.org/en/v1/manual/performance-tips/#Performance-critical-code-should-be-inside-a-function).

Note also that your [argument-type and return-type declarations don’t impact performance](https://docs.julialang.org/en/v1/manual/functions/#Argument-type-declarations) and it is probably clearer to omit them. Similarly with the explicit calls to `Float64(...)`.

(Also, I would recommend passing vectors to `hcubature` as [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) so that the compiler knows the dimensionality of the integral.)

I looked into optimizing your code, but I find that I cannot replicate your results. I get:

```julia
julia> @time @fastmath spinonsd(1.0)
 21.642360 seconds (453.33 M allocations: 18.803 GiB, 6.77% gc time)
0.08557197546327772

```

which is considerably different from the numbers you reported above. Did you miscopy something in your code?

---

_[View the full topic](https://discourse.julialang.org/t/unreasonable-slow-speed-in-numerical-integration/89910)._