Unreasonable slow speed in numerical integration

General Usage Performance
, , ,
15

So in this problem, because of the way the integrals separate, it turns out to be much much more efficient to integrate first over omega and then integrate over kx and ky. I’m assuming you want your frequency integral to go over the range -Inf to Inf, but if you actually only wanted to integrate from -3 to 3 the timings are the same.

#+begin_src julia
using HCubature, QuadGK, BenchmarkTools
function spinonenergy(kx, ky)
    2 * 0.05 * (cos(kx)+cos(-kx/2+sqrt(3)/2*ky)+cos(kx/2+sqrt(3)/2*ky)) + 0.04
end

function Sigma(η)
    hcubature((0.0,0.0),(2*pi,4*pi/sqrt(3))) do (kx, ky)
        ϵ = spinonenergy(kx, ky)
        quadgk(-Inf, Inf) do ω
            -1/pi*sqrt(3)/2*1/(2*pi)^2 * imag(inv(ω - ϵ + im*η))
        end[1]
    end[1]
end

@benchmark Sigma(0.001) # the complex offset you originally used
#+end_src

#+RESULTS:
: BenchmarkTools.Trial: 10000 samples with 1 evaluation.
:  Range (min … max):  193.428 μs …  2.202 ms  ┊ GC (min … max): 0.00% … 89.97%
:  Time  (median):     194.439 μs              ┊ GC (median):    0.00%
:  Time  (mean ± σ):   196.813 μs ± 44.045 μs  ┊ GC (mean ± σ):  0.49% ±  2.01%
: 
:   ▃▇█▇▆▅▃▃▂▁▁  ▂▃▃▄▄▄▃▃▂▂▁▁                                    ▂
:   █████████████████████████▇█▇▆▆▆▅▇▇██▇▇▇█▇▆▅▆▅▄▅▆▄▅▅▅▆▄▅▄▄▄▂▅ █
:   193 μs        Histogram: log(frequency) by time       210 μs <
: 
:  Memory estimate: 29.78 KiB, allocs estimate: 71.

So here, I’m able to do the whole integral over kx, ky, and omega in a tenth of a milisecond, and do it with the incredibly high default accuracies that QuadGK and HCubature default to.

Furthermore, we can use Richardson.jl to extrapolate this result to η = 0:

#+begin_src julia
using Richardson
@btime extrapolate(Sigma, 1.0)
#+end_src

#+RESULTS:
:   105.489 μs (114 allocations: 18.56 KiB)
: (0.9999999999998697, 1.3034018309099338e-13)

which is actually faster than plugging in η = 0.001 explicitly, and more accurate.

4 Likes