Best way to use/implement Debye function in Julia?

General Usage
21

Right, improving this one cut computation time by a factor 2 for the x < pi branch of the if-clause. The other branch was also sped up a bit by substituting the powers, but that improvement was rather minor. I did have to introduce a check whether to use big or not, because for certain x it was needed to achieve the desired accuracy.
Still, performance is a far cry from TOMS or quadrature. One thing I don’t understand is why the benchmark tool reports several “allocations” (29 in the first, 237 in the second branch of the if block), but none in the TOMS and quadrature functions. What is it that is being allocated here, and may this be the cause for the poor performance?

22

I haven’t looked at your code, but are you aware that most BigFloat operations allocate a new BigFloat object? If you want to avoid this, use MutableArithmetics.jl.

23

The quadrature-based solution can be arbitrarily accurate, but gets slower the further you go from zero. quadgk from QuadGK.jl takes tolerance parameters which are used to set the desired accuracy.

An important thing that was missed in the discussion above was that one should use expm1(x) instead of exp(x) - 1 in the integrand, as that’s more accurate for tiny x and allows quadgk to be faster.

Polynomial approximations may be useful as an even faster alternative to QuadGK.jl, depending on what interval you’re targeting and what’s your desired accuracy.

Could you give an estimate for an upper bound for your interval? The interval doesn’t necessarily need to include all your inputs. If, say, 99% of your inputs were solved using a polynomial approximation, and the other 1% using quadrature, that could still provide a speed-up compared to using quadrature for all inputs.

Also, what is your desired accuracy? Usually that’s expressed in terms of maximum allowed relative error, or, alternatively, as the required number of accurate bits in the significand of the floating-point number.

24

Here’s a quadrature-based solution for D_1, D_2, D_3 and D_4:

module DebyeQuadrature

import QuadGK, MutableArithmetics

pow(x, ::Val{1}) = x
pow(x, ::Val{2}) = x*x
pow(x, ::Val{3}) = x^3
pow(x, ::Val{4}) = x^4

function pow(x::BigFloat, ::Val{3})
  o! = MutableArithmetics.operate!
  y = x*x
  o!(*, y, x)::BigFloat
end

function pow(x::BigFloat, ::Val{4})
  o! = MutableArithmetics.operate!
  y = x*x
  o!(*, y, y)::BigFloat
end

function evaluate(
  x, dord::Val{DOrd}, relative_tolerance, quadrature_order::Int,
) where {DOrd}
  DOrd::Int
  f = t -> pow(t, dord)/expm1(t)
  if iszero(x)
    one(x)
  else
    let z = false
      i = QuadGK.quadgk(
        f, z, x, maxevals = Inf, atol = z, rtol = relative_tolerance,
        order = quadrature_order,
      )
      DOrd/pow(x, dord)*first(i)
    end
  end
end

end

The quadrature order affects performance, a small value like 7 should be enough for Float64, while BigFloat performance benefits from larger values like 35, or even in the hundreds. Larger values increase compilation time.

25

For my purposes, I rarely expect values for x above 7; in most cases, they should barely exceed 2 and indeed often be smaller than 1. This is because my applications do not generally target cryogenic temperatures but rather the intermediate-to-high-T realm.

1 Like
26

What about pre-computing sufficiently much intermeadiate points to maximal desired accuracy, and integrate from the nearest one?

27

generally you can get better results just by finding the series expansion around infinity and using that for large values.

1 Like
28

Yes, for a given fixed accuracy (e.g. Float64), typically high-performance implementations of special functions stitch together a set of different approximations — Taylor series near zeros, minimax polynomials or rational approximants in finite real intervals, continued-fractions and/or asymptotic series for large arguments. It takes a lot of tuning to find the right set of approximations and to find good crossover points to switch between approximations.

29

A Julia package for evaluating Debye functions is being registered, it will live at DebyeFunctions.jl once it’s registered:

Git repo at Gitlab: Neven Sajko / DebyeFunctions.jl · GitLab

Please test it out and maybe suggest some improvements.

EDIT: well, an obvious and easy improvement would be to add a few more polynomial approximations (for each Debye function order). The trouble is choosing when to stop adding more polynomials, I guess.

1 Like
30

for orders 1,2,3, there are specific formulations that are defined in terms of the polylogarithm and the riemann’s zeta function:

import PolyLog,SpecialFunctions,LogExpFunctions
function debye3(x)
  N = 3
  ζ = SpecialFunctions.zeta(N+1)
  c1 = LogExpFunctions.log1mexp(x) #log(1-exp(x))
  c2 = 3*PolyLog.reli2(ex)
  c3 = -6*PolyLog.reli3(ex)
  c4 = 6*PolyLog.reli4(ex)
  ex = exp(x)
  return x*x*x*(-0.25*x + c1) + c2 + c3 + c4 - 6*ζ
end

(from https://mathworld.wolfram.com/DebyeFunctions.html)