Fast implementation of incomplete extended gamma function, QuadGK?


I am interested in developing a fast routine to compute the following extended incomplete gamma function for various pairs (x, q).

\gamma(\alpha, \kappa, p, q)(x) = \int_{0}^x t^{\alpha-1} \exp( -\frac{t^\kappa}{p} - \frac{q}{t}) \mathrm{d}t with \alpha \in \mathbb{R}, p>0, \kappa \in \mathbb{R}_{\neq 0}, q \geq 0

I would also be interested in computing its inverse \gamma^{-1}(\alpha, \kappa, p, q)(z) for z \geq 0.

I attach my current implementation based on QuadGK. I would appreciate any feedback to accelerate computations or reduce allocations.

Can it be beneficial to define the function \gamma as the solution of the ODE \frac{\mathrm{d} \gamma(\alpha, \kappa, p, q)(x) }{\mathrm{d} x} = x^{\alpha-1} \exp( -\frac{x^\kappa}{p} - \frac{q}{x}) , then pass this ODE to DifferentialEquations.jl.
Then we can compute its inverse by using a root-finding algorithm combined with the interpolation for free tools from DifferentialEquations

unnorm_pdf!(out, α::Float64, κ::Float64, p::Float64, q::Float64, t) = out .= t^(α-1)*exp(-t^κ/p - q/t)

# Routine to compute the un-normalized probability density function
function unnorm_pdf(α::Float64, κ::Float64, p::Float64, q::Float64, t)
    out = zeros(eltype(t),1)
    unnorm_pdf!(out, α, κ, p, q, t)
    return out

# Routine to compute the normalizing constant of the probability density function
function Z!(out, α::Float64, κ::Float64, p::Float64, q::Float64) 
    quadgk!((b,a)->unnorm_pdf!(b,α, κ, p, q, a),out, 0.0, Inf)

# Routine to compute the normalizing constant of the probability density function
function Z(α::Float64, κ::Float64, p::Float64, q::Float64)
    out = zeros(eltype(t),1)
    Z!(out, α, κ, p, q)
    return out

struct ExtendedGamma
    q::MVector{1, Float64}
    # Normalizing constant
    Z::MVector{1, Float64}
    cache::MVector{1, Float64}

function ExtendedGamma(α::Float64, κ::Float64, p::Float64, q::Float64)
    cache = MVector{1, Float64}(0.0)
    # Compute normalizing constant
    Z!(cache,α, κ, p, q)
    ZΓ = copy(cache)
    return ExtendedGamma(α, κ, p, MVector{1, Float64}(q), ZΓ, cache)

unnorm_pdf!(out, Γ::ExtendedGamma, t::Real) = unnorm_pdf!(out, Γ.α, Γ.κ, Γ.p, Γ.q[1], t)

function unnorm_pdf(Γ::ExtendedGamma, t::Real)
    unnorm_pdf!(Γ.cache, Γ, t)
    return Γ.cache
function pdf!(out, Γ::ExtendedGamma, t::Real) 
    unnorm_pdf!(out, Γ, t)
    out ./= Γ.Z[1]

function pdf(Γ::ExtendedGamma, t::Real)
    pdf!(Γ.cache, Γ, t)
    return Γ.cache

Z!(out, Γ::ExtendedGamma) = Z!(out, Γ.α, Γ.κ, Γ.p, Γ.q[1]) 

function Z(Γ::ExtendedGamma)
    Z!(Γ.cache, Γ)
    return Γ.cache

update_Z!(Γ::ExtendedGamma) = Z!(Γ.Z, Γ)

function cdf!(out, Γ::ExtendedGamma, x) 
    quadgk!((b,a)-> unnorm_pdf!(b, Γ, a),out, 0.0, x)[1]
    out ./= Γ.Z[1]

function cdf(Γ::ExtendedGamma, x) 
    cdf!(Γ.cache, Γ, x)
    return Γ.cache

# Complementary of the CDF to 1
function ccdf!(out, Γ::ExtendedGamma, x) 
    quadgk!((b,a)-> unnorm_pdf!(b, Γ, a),out, x, Inf)[1]
    out ./= Γ.Z[1]

# Complementary of the CDF to 1
function ccdf(Γ::ExtendedGamma, x) 
    ccdf!(Γ.cache, Γ, x)
    return Γ.cache


using LinearAlgebra
using QuadGK
using StaticArrays
using BenchmarkTools

α = 2.0
κ = 2.0
p = 1.0
q = 0.8
Γstruct = ExtendedGamma(α, κ, p, q)

# Time evaluation of the incomplete expanded Gamma function
@btime begin 
    cdf!(Γstruct.cache, Γstruct, 1.0)


 2.185 μs (18 allocations: 672 bytes)
1-element MVector{1, Float64} with indices SOneTo(1):

Basically quadrature is never used in optimized special-function implementations, though it’s great to have for comparison and bootstrapping purposes.

Instead, optimized special functions usually employ a combination of polynomial and rational expansions, e.g. continued-fraction expansions for large argument and Taylor series near roots. The tricky part is figuring out what approximation to use where, and how to stitch them together.

Unfortunately, this becomes more and more difficult to implement well as the special functions have more parameters.

Generally, I would recommend against ODE solvers to perform integrals, since that discards a lot of structure from the problem.

If you fix the parameters \alpha, \kappa, p, q and want to compute \gamma(x) and its inverse for some range of x values, then you can do much better. That is, suppose you are computing \gamma^{-1}(y) many times for given \alpha, \kappa, p, q, within a finite range of x and y values, and you can afford some startup cost.

Then, instead of using the DifferentialEquations interpolant, I would just construct a smooth polynomial approximation from a set of quadrature evaluations, e.g. a Chebyshev interpolant using a package like FastChebInterp.jl or ApproxFun.jl. This converges exponentially fast for smooth functions, much more accurate than an ODE interpolant. Given this interpolant for \gamma(x), you can then construct an interpolant for \gamma^{-1}(y) by applying Newton’s method to y values at Chebyshev points in your interval of interest. Thereafter, computing \gamma^{-1}(z) for any z in your interval is just a fast polynomial evaluation.

See also Approximate inverse of a definite integral - #3 by stevengj — as noted there, you can just do Newton’s method directly on the QuadGK results and then compute a Chebyshev interpolant of \gamma^{-1}(y) directly, skipping the intermediate step of constructing an interpolant for interpolant for \gamma(x).


Is there some kind of standard référence on this subject that we could follow for some new special function ? Only out of personal curiosity, I would love to read such a book :slight_smile: Something that would guide me into crafting, let say, a good Gamma function implémentation (and that would be applicable to over functions too ofc)

Thank you @stevengj for your detailed answer.

Using ApproxFun, we build interpolant for $d

using ApproxFun
chebΓ_integrand = Fun(x->pdf(Γstruct, x)[1], Interval(0,+Inf))
chebΓ = cumsum(chebΓ_integrand )

The function chebΓ_integrand decays to 0 at +\infty and is always positive. Is there a preferred option in ApproxFun.jl for this class of functions, or Interval(0, +Inf) sufficient?

Is it true that the operator chebΓ = cumsum(chebΓ_integrand ) defines the integrated function \int_0^x \text{chebΓ_integrand}(t) dt. Is the offset constant set such that the integral operator equals to 0 at the lower bound of the domain?

Can you elaborate on this point? I am not sure how to build the interpolant for \gamma^{-1}(y) from the interpolants chebΓ and chebΓ_integrand.

Are you suggesting to apply Newton’s method with x - chebΓ(x)/chebΓ_integrand(x)? P.S. the domain is [0, + \infty)

I have seen this book mentioned in another post: