Derivatives of the incomplete gamma functions with respect to parameter a are missing

[devmotion]

Copied from #305:

I added ChainRules definitions for gamma(a, x), loggamma(a, x), and gamma_inc according to https://functions.wolfram.com/GammaBetaErf/GammaRegularized/introductions/Gammas/ShowAll.html . Similar to the Bessel functions, derivatives with respect to the first argument a are not implemented since they are given in terms of a (regularized) hypergeometric function.

This issue is a reminder that these derivatives are missing. Moreover, the issue can be linked in the error message displayed to users when they try to access and perform computations with these missing derivatives.

[lrnv]

Looks like the error message does not link this issue yet, I had to find it by hand after hitting the "no derivative w.r.t a" problem in my code.

From https://en.wikipedia.org/wiki/Incomplete_gamma_function#Lower_incomplete_gamma_function it looks like there exists a way to compute this derivative through the Meijer G function. Is this funciton implmeented somewhere ? Is there a way to patch this even a little dirty ?

Edit : Might have to do with #265. Also, the package HypergeometricFunctions.jl might help if the regularized hypergeometric function proposed by wolframalpha for the derivative is indeed in this package.

Edit: Link to discussion in discourse : https://discourse.julialang.org/t/gamma-inc-derivatives/93148

Edit: I got it :

using SpecialFunctions, HypergeometricFunctions
function Γ(a,x)
    return gamma_inc(a,x)[2]
end
function ∂Γ_∂a(a,x)
    g = gamma(a)
    dg = digamma(a)
    G = Γ(a,x)
    lx = log(x)
    r = pFq([a,a],[a+1,a+1],-x)
    return a^(-2) * x^a * r/g + (1 - G)*(dg - lx)
end
function emp(a,x)
    ϵ = 1e-5
    return (Γ(a+ϵ,x) - Γ(a-ϵ,x))/2ϵ
end


a = 1.9
x = 1.8
Γ(a,x)
(∂Γ_∂a(a,x), emp(a,x))

# works correctly for x = 1 or a=1 or a=2. Does not work correclty for other valyes i tried. 
# which might mean that my formula is somewhat wrong. 
a = 0.1:0.1:5
x = 0.1:0.1:5
diff(a,x) = abs(∂Γ_∂a(a,x) - emp(a,x))
maximum(diff.(a,x')) # 2e-10 on my machine. 

Is there a chance that something like this could be added to the package ?

[sethaxen]

With package extensions, perhaps it would be possible to add a SpecialFunctionsChainRulesCoreHypergeometricFunctionsExt that defines the rules with respect to the first argument. Then either a fallback rule could be added in the SpecialFunctionsChainRulesCoreExt extension to inform the user they need to manually load HypergeometricFunctions to use that functionality, or perhaps Base.Experimental.register_error_hint could be used to warn the user.

[devmotion]

Is this really possible? HypergeometricFunctions depends on SpecialFunctions, so I wonder if it can be a weak dependency at all (I am sure it can't be a proper dependency at least since it would create circular dependencies?). Maybe it could be an extension of HypergeometricFunctions but this seems a bit more far fetched.

[sethaxen]

Is this really possible? HypergeometricFunctions depends on SpecialFunctions, so I wonder if it can be a weak dependency at all (I am sure it can't be a proper dependency at least since it would create circular dependencies?).

Ah, no, it seems circularity is not okay even for weak deps:

┌ Warning: Circular dependency detected. Precompilation will be skipped for:
│   DualNumbers [fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74]
│   SpecialFunctionsChainRulesCoreHypergeometricFunctionsExt [58512b41-50c9-55ec-9df1-0539e334355e]
│   HypergeometricFunctions [34004b35-14d8-5ef3-9330-4cdb6864b03a]
└ @ Pkg.API ~/.julia/juliaup/julia-1.9.1+0.x64.linux.gnu/share/julia/stdlib/v1.9/Pkg/src/API.jl:1234
Precompiling project...
  1 dependency successfully precompiled in 1 seconds. 19 already precompiled. 3 skipped due to circular dependency.

Maybe it could be an extension of HypergeometricFunctions but this seems a bit more far fetched.

I agree, this is strange.

All we really need is $_2F_2$, which in HypergeometricFunctions uses the default pFq, which requires only 3 methods:

• https://github.com/JuliaMath/HypergeometricFunctions.jl/blob/089e208d5869dd09b0f12010dd28d1026a18a0ca/src/generalized.jl#L15-L21
• https://github.com/JuliaMath/HypergeometricFunctions.jl/blob/089e208d5869dd09b0f12010dd28d1026a18a0ca/src/specialfunctions.jl#L529-L541
• https://github.com/JuliaMath/HypergeometricFunctions.jl/blob/089e208d5869dd09b0f12010dd28d1026a18a0ca/src/weniger.jl#L199-L346

I wonder if these 3 methods could be duplicated here as internal functions to break the circularity.

[lrnv]

I remember that at the time, in this issue JuliaMath/HypergeometricFunctions.jl#50 , @MikaelSlevinsky gave a faster version for this derivative (much faster than the upper code). Still depending on HypergeometricFunctions.jl though so I do not really know where it should be implemented.

[sethaxen]

Thanks for the pointer, @lrnv . Agreed it is likely faster to differentiate through a hypergeometric function if supported than to implement the derivative separately.

@devmotion any objections to the plan proposed in #317 (comment)?

[devmotion]

Generally, I'm not a fan of duplicating code, in particular not in dependent packages. But I don't have a better idea right now (and, even though a bit weird, the duplication issue could be resolved if HypergeometricFunctions starts re-exporting the functionality in SpecialFunctions). I'm not sure if e.g. breaking up SpecialFunctions is smaller packages (which are bundled by SpecialFunctions) would help in this case - if HypergeometricFunctions continues depending on SpecialFunctions I assume we would still run into these circular dependency issues.

Based on the previous comment, I think we should benchmark your suggestion against using HypergeometricFunctions.M (performance-wise and numerically) as suggested in JuliaMath/HypergeometricFunctions.jl#50. This alternative approach would also only require to copy a few lines.

[oschulz]

Just curious, any news on this? (I recently ran into a case that needs both partial derivatives of incomplete gamma.)

[lrnv]

@oschulz you should use the proposal from JuliaMath/HypergeometricFunctions.jl#50 (comment) that's the best options right now.

[oschulz]

Thanks @lrnv !

[mleprovost]

Hello,

I have a related question: I would like to perform AD for the function gamma_inc_inv at https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/src/gamma_inc.jl:
ForwardDiff.derivative(x -> gamma_inc_inv(x, 0.1, 0.9), 0.1)
Are there routines in HypergeometricFunctions.jl to make this calculation that support AD?

Thank you for your help,

[MikaelSlevinsky]

Hi @mleprovost, it looks like the inverse incomplete gamma function uses series reversion-like techniques (e.g.

SpecialFunctions.jl/src/gamma_inc.jl

Lines 744 to 770 in ff76c3a

	"""
	gamma_inc_inv_psmall(a, logr)
	Compute `x0` - initial approximation when `p` is small.
	Here we invert the series in Eqn (2.20) in the paper and write the inversion problem as:
	```math
	x = r\\left[1 + a\\sum_{k=1}^{\\infty}\\frac{(-x)^{n}}{(a+n)n!}\\right]^{-1/a},
	```
	where ``r = (p\\Gamma(1+a))^{1/a}``
	Inverting this relation we obtain ``x = r + \\sum_{k=2}^{\\infty}c_{k}r^{k}``.
	"""
	function gamma_inc_inv_psmall(a::Float64, logr::Float64)
	r = exp(logr)
	ap1 = a + 1.0
	ap1² = ap1*ap1
	ap1³ = ap1*ap1²
	ap1⁴ = ap1²*ap1²
	ap2 = a + 2.0
	ap2² = ap2*ap2
	ck1 = 1.0
	ck2 = 1.0/(1.0 + a)
	ck3 = 0.5*(3*a + 5)/(ap1²*(a + 2))
	ck4 = (1.0/3.0)*(@horner(a, 31, 33, 8))/(ap1³*ap2*(a + 3))
	ck5 = (1.0/24.0)*(@horner(a, 2888, 5661, 3971, 1179, 125))/(ap1⁴*ap2²*(a + 3)*(a + 4))
	x0 = @horner(r, 0.0, ck1, ck2, ck3, ck4, ck5)
	return x0
	end

) for small or large arguments. Unless there's a hypergeometric structure to those reverted series, then I'm not sure how the HypergeometricFunctions package could help.

[mleprovost]

Thank you for your help. Should I open a separate issue?