DLRGrid throws ArgumentError in the high temperature limit

[HugoStrand]

Dear Kun,

Thank you for the nice Julia implementation of DLR. I am trying to use Lehmann.jl in a project but have stumbled on an issue when calling the DLRGrid builder at high temperatures.

Here is an example code that breaks throwing an ArgumentError

using Lehmann: DLRGrid

for n in range(1, 10)
    β = 1000.0/2^n
    dlr = DLRGrid(Euv=1., β=β, isFermi=true, rtol=1e-9, rebuild=true, verbose=false)
    @show β, n, length(dlr.τ)
end

For the output please see below.

Do you think it would be possible to "harden" the behaviour of the grid builder so that it performs ok also for high temperatures?

Cheers, Hugo

(β, n, length(dlr.τ)) = (500.0, 1, 41)
(β, n, length(dlr.τ)) = (250.0, 2, 37)
(β, n, length(dlr.τ)) = (125.0, 3, 31)
(β, n, length(dlr.τ)) = (62.5, 4, 25)
(β, n, length(dlr.τ)) = (31.25, 5, 20)
(β, n, length(dlr.τ)) = (15.625, 6, 16)
(β, n, length(dlr.τ)) = (7.8125, 7, 13)
ArgumentError: reducing over an empty collection is not allowed

Stacktrace:
  [1] _empty_reduce_error()
    @ Base ./reduce.jl:301
  [2] reduce_empty(op::Function, #unused#::Type{Float64})
    @ Base ./reduce.jl:311
  [3] mapreduce_empty(#unused#::typeof(identity), op::Function, T::Type)
    @ Base ./reduce.jl:345
  [4] reduce_empty(op::Base.MappingRF{typeof(identity), typeof(max)}, #unused#::Type{Float64})
    @ Base ./reduce.jl:331
  [5] reduce_empty_iter
    @ ./reduce.jl:357 [inlined]
  [6] mapreduce_empty_iter(f::Function, op::Function, itr::Vector{Float64}, ItrEltype::Base.HasEltype)
    @ Base ./reduce.jl:353
  [7] _mapreduce(f::typeof(identity), op::typeof(max), #unused#::IndexLinear, A::Vector{Float64})
    @ Base ./reduce.jl:402
  [8] _mapreduce_dim
    @ ./reducedim.jl:330 [inlined]
  [9] #mapreduce#725
    @ ./reducedim.jl:322 [inlined]
 [10] mapreduce
    @ ./reducedim.jl:322 [inlined]
 [11] #_maximum#743
    @ ./reducedim.jl:894 [inlined]
 [12] _maximum
    @ ./reducedim.jl:894 [inlined]
 [13] #_maximum#742
    @ ./reducedim.jl:893 [inlined]
 [14] _maximum
    @ ./reducedim.jl:893 [inlined]
 [15] #maximum#740
    @ ./reducedim.jl:889 [inlined]
 [16] maximum
    @ ./reducedim.jl:889 [inlined]
 [17] testInterpolation(dlrGrid::DLRGrid{Float64, :none}, τ::Lehmann.Discrete.CompositeChebyshevGrid, ω::Lehmann.Discrete.CompositeChebyshevGrid, kernel::Matrix{Float64}, print::Bool)
    @ Lehmann.Discrete ~/.julia/packages/Lehmann/v7X4o/src/discrete/kernel.jl:137
 [18] build(dlrGrid::DLRGrid{Float64, :none}, print::Bool)
    @ Lehmann.Discrete ~/.julia/packages/Lehmann/v7X4o/src/discrete/builder.jl:119
 [19] _build!(dlrGrid::DLRGrid{Float64, :none}, folder::Nothing, filename::String, algorithm::Symbol, verbose::Bool)
    @ Lehmann ~/.julia/packages/Lehmann/v7X4o/src/dlr.jl:338
 [20] DLRGrid(Euv::Float64, β::Float64, rtol::Float64, isFermi::Bool, symmetry::Symbol; rebuild::Bool, folder::Nothing, algorithm::Symbol, verbose::Bool, dtype::Type)
    @ Lehmann ~/.julia/packages/Lehmann/v7X4o/src/dlr.jl:172
 [21] #DLRGrid#7
    @ ~/.julia/packages/Lehmann/v7X4o/src/dlr.jl:189 [inlined]
 [22] top-level scope
    @ ./In[20]:5

[kunyuan]

Dear Hugo,

I apologize for the late reply. I have been traveling recently.

Thank you very much for creating the issue. It is a bug that needs to be fixed. We will release a bug fix in a couple of days.

[HugoStrand]

Dear Kun,

Thank you for the reply. Is there something I can do to fix this?

I am guessing that one need to have some minimum number of dyadic refinements in the initial discretization of the analytic continuation kernel. Where should one start to dig, to tweak that?

Cheers, Hugo

[fsxbhyy]

Dear Hugo,

Sorry for the late reply. Thank you for your comment on this issue! We have made a pull request to fix this issue, which will be soon updated once it is finalized. Meanwhile, I also what to comment that at high temperature, if you do not rebuild the DLR grid, it will automatically identify the existing DLR files and use the one that has the closest Lambda.

For example, if the smallest DLR file you can get is Lambda = 1000, rtol = 1e-9, this line below

dlr = DLRGrid(Euv=1., β=10, isFermi=true, rtol=1e-9, verbose=false)

will generate the grid with Lambda = 1000 accordingly. The number of grid points will be larger then rebuilding a new grid at Lambda = Euv * β = 10. But normally at high temperature Lambda is pretty small already, and given the logarithmic scaling of DLR rank with Lambda, the efficiency drop should be pretty small.

For this reason, in many situation it is not necessary to rebuild DLR grid each time. Also, if you encounter more problems in the future, you can always use this option as a backup solution when we try to fix the problem.

Please let me know if you have any other thoughts!

Cheers,

Tao Wang

[fsxbhyy]

Dear Hugo,

The new changes have been officially released. Please let us know if you have any questions!

Cheers,

Tao Wang