Unreasonable slow speed in numerical integration

General Usage Performance
, , ,
2

First, you should realize that the Matlab integrate function 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 on measuring time. You also need to be careful about benchmarking in global scope, since Julia is oriented towards performance-critical code in functions.

Note also that your argument-type and return-type declarations don’t impact performance 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 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> @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?

38 Likes