ForwardDiff.Dual recurses with logabsgamma

[matbesancon]

(TestIndicator) pkg> add ForwardDiff SpecialFunctions
  Updating registry at `~/.julia/registries/General`
  Updating git-repo `https://github.com/JuliaRegistries/General.git`
 Resolving package versions...
  Updating `~/Documents/CTests/TestIndicator/Project.toml`
  [f6369f11] + ForwardDiff v0.10.3
  [276daf66] + SpecialFunctions v0.8.0
  Updating `~/Documents/CTests/TestIndicator/Manifest.toml`
 [no changes]

julia> import ForwardDiff, SpecialFunctions

julia> ForwardDiff.Dual(3.5)
Dual{Nothing}(3.5)

julia> d = ForwardDiff.Dual(3.5)
Dual{Nothing}(3.5)

julia> SpecialFunctions.logabsgamma(d)
ERROR: StackOverflowError:
Stacktrace:
 [1] logabsgamma(::ForwardDiff.Dual{Nothing,Float64,0}) at /home/mbesancon/.julia/packages/SpecialFunctions/ne2iw/src/gamma.jl:614 (repeats 80000 times)

[Sumegh-git]

@simonbyrne any suggestions?

[simonbyrne]

I think the solution is to add the rules to https://github.com/JuliaDiff/DiffRules.jl/blob/a9b00c81ac0d103888b124b2b9f58a7c815cebf9/src/rules.jl#L67-L70.

[torfjelde]

May I ask, how do you go about adding a diffrule for logabsgamma, which returns a tuple rather than a single number?

[torfjelde]

Seems like we need the following:

using ForwardDiff: Dual, value, partials
import SpecialFunctions: logabsgamma
function logabsgamma(x::Dual{T}) where T
    val, sgn = logabsgamma(value(x))

    return (Dual{T}(val, partials(x) * DiffRules._abs_deriv(value(x)) * digamma(abs(value(x)))), )
end

Is this something that should go into SpecialFunctions.jl?

EDIT: The above together with adding in

    @define_diffrule SpecialFunctions.loggamma(x) =
        :(  SpecialFunctions.digamma($x)  )

where @simonbyrne suggested, results in

julia> ForwardDiff.derivative(z -> first(logabsgamma(z)), -1.0)
0.5772156649015315

julia> ForwardDiff.derivative(z -> first(logabsgamma(z)), 1.0)
-0.5772156649015315

julia> ForwardDiff.derivative(z -> loggamma(z), 1.0)
-0.5772156649015315

which seems good 👍

[ararslan]

I believe @andreasnoack has already fixed this in DiffRules.

[torfjelde]

Hah, didn't notice; great! But I don't believe that change fixes the logabsgamma issue?

EDIT: yeah I just tried with the updated master of DiffRules; logabsgamma of a Dual still causes StackOverflowError

[torfjelde]

From what I understand, it seems like either the above needs to be overloaded in ForwardDiff.jl, or

SpecialFunctions.jl/src/gamma.jl

Line 614 in 61f9058

logabsgamma(x::Real) = logabsgamma(float(x))

needs to be removed?

[tpapp]

I think the underlying issue is this conversion combined with the fallback of float (cf JuliaLang/julia#26552).

[tpapp]

Cf JuliaDiff/ForwardDiff.jl#419

[andreasnoack]

I've only fixed loggamma and not logabsgamma but it's the former that is needed in Distributions. The derivative definition of the latter is a bit more involved but I guess the one you proposed works. Please open a PR against DiffRules.

While stack overflows are annoying I'm not sure that it's necessarily a bug. It's quite convenient to have a fallback definition that converts Reals with float. The issue here is the missing method in DiffRules.

I think the underlying issue is this conversion combined with the fallback of float (cf JuliaLang/julia#26552)

I don't think that method is causing the recursion. I think it's the method right above.

[tpapp]

Technically it is of course always a specific method (eg logabsgamma) that recurses, because it is expecting something from float that it doesn't provide. Cf my comment here JuliaLang/julia#26552 (comment)

[torfjelde]

I've only fixed loggamma and not logabsgamma but it's the former that is needed in Distributions.

Ah, sweet! Very much looking forward to getting these versions all bumped up.

Please open a PR against DiffRules.

I'm not too familiar with how DiffRules works, but seems to me that you can only define a rule for functions with a single return value which should be <: Real? Since this has multiple return values, we only really want to define how things propagate through logabsgamma. That's why I thought I had to reach for ForwardDiff (which means it can't be put into DiffRules, since ForwardDiff itself depends on DiffRules). So should the PR be to ForwardDiff maybe?

While stack overflows are annoying I'm not sure that it's necessarily a bug. It's quite convenient to have a fallback definition that converts Reals with float.

Sorry, I meant as in the cause of the issue in this case. I definitively see the convience of this method. But the PR @tpapp has submitted to ForwardDiff seems like the way to go for fixing the issue of logabsgamma resulting in a StackOverflowError, though we'd still have to define logabsgamma(x::Dual) to be able to propagate derivative-information through this method, right?

[LoganAMorrison]

Was there any resolution to this issue? I'm encountering the same issue with besselkx.

[andreasnoack]

Generally, the solution is to add derivative definitions to DiffRules if this happens.

[LoganAMorrison]

@andreasnoack Okay, so something like this needs to be added?

@define_diffrule SpecialFunctions.besselkx(n, x) =
    :(SpecialFunctions.besselkx($n, $x) - 1 / 2 * SpecialFunctions.besselkx($n - 1, $x) - 1 / 2 * SpecialFunctions.besselkx($n + 1, $x))

[devmotion]

With SpecialFunctions 1.7.0 (or more concretely #347) the original example throws a MethodError instead of a StackOverflowError:

julia> using ForwardDiff, SpecialFunctions

julia> d = ForwardDiff.Dual(3.5)
Dual{Nothing}(3.5)

julia> SpecialFunctions.logabsgamma(d)
ERROR: MethodError: no method matching _logabsgamma(::ForwardDiff.Dual{Nothing, Float64, 0})
Closest candidates are:
  _logabsgamma(::Float64) at /home/david/.julia/packages/SpecialFunctions/5CocL/src/gamma.jl:601
  _logabsgamma(::Float32) at /home/david/.julia/packages/SpecialFunctions/5CocL/src/gamma.jl:606
  _logabsgamma(::Float16) at /home/david/.julia/packages/SpecialFunctions/5CocL/src/gamma.jl:611
  ...
Stacktrace:
 [1] logabsgamma(x::ForwardDiff.Dual{Nothing, Float64, 0})
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/5CocL/src/gamma.jl:599
 [2] top-level scope
   @ REPL[9]:1

[andreasnoack]

I don't think there is more to do here. Any further action should happen in DiffRules.