#481 iszero workaround

[cgeoga]

Like a couple other people here, I am experiencing some issues with #481. My issue is in BesselK.jl, where I have some very specific methods that branch on an iszero(x) that were broken (for example).

I understand the logic of the change and clearly an argument can be made for it, but it is disappointing that the only solution I can come up with now is to bring ForwardDiff into the dependency tree and define ugly methods

  _iszero(x) = iszero(x)
  _iszero(x::ForwardDiff.Dual) = _iszero(ForwardDiff.value(x))

I am sure that it's too late to argue for just reverting #481, but it would be nice if there were some way to check that ForwardDiff.value(x) is zero without having to bring ForwardDiff into the tree. Is there some other generic comparison method in Base or something that ForwardDiff can add methods to that checks for zero values? Or is there some getter function that one can use to get ForwardDiff.value(x)? I'm very happy to write a different check than iszero---but I think it would be nice for the check to be possible without bringing the internal getters of ForwardDiff.jl into scope.

Excuse me if I'm using any incorrect terminology and thanks in advance for any thoughts you can share on this.

[mcabbott]

Sorry about breaking things.

What's the smallest example you can write of a sort-of user-facing function which will exercise this line? Something which ought to have a derivative.

[cgeoga]

No apologies necessary at all---I worship this package and really appreciate how much time people like you put into making it great.

Here's an example call where changing my linked _iszero to iszero creates an error:

using BesselK, ForwardDiff
ForwardDiff.derivative(_v->adbesselkxv(_v, 0.25), 2.0)

EDIT: I tweaked the example to use the function that end-users would use, not the internal method that it hits. Sorry to misread your question.

[mcabbott]

OK so the "before" picture is

julia> ForwardDiff.derivative(_v->adbesselkxv(_v, 0.25), 2.0)
2.225724617743265

julia> FiniteDifferences.grad(central_fdm(7, 1), _v->adbesselkxv(_v, 0.25), 2.0)[1]
2.225724617743085

(jl_Pl1N2w) pkg> st
Status `/private/var/folders/yq/4p2zwd614y59gszh7y9ypyhh0000gn/T/jl_Pl1N2w/Project.toml`
⌃ [432ab697] BesselK v0.5.2
  [26cc04aa] FiniteDifferences v0.12.25
⌃ [f6369f11] ForwardDiff v0.10.32

and the "after" one is

julia> ForwardDiff.derivative(_v->adbesselkxv(_v, 0.25), 2.0)
ERROR: Term tolerance 1.0e-12 not reached in 100 iters.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] temme_pair(v::ForwardDiff.Dual{ForwardDiff.Tag{var"#7#8", Float64}, Float64, 1}, z::Float64, maxit::Int64, tol::Float64, modify::Bool)
   @ BesselK ~/.julia/packages/BesselK/GkDiL/src/besk_temme.jl:176
 [3] _besselk_temme(v::ForwardDiff.Dual{ForwardDiff.Tag{var"#7#8", Float64}, Float64, 1}, z::Float64, maxit::Int64, tol::Float64, mod::Bool)
   @ BesselK ~/.julia/packages/BesselK/GkDiL/src/besk_temme.jl:192
 [4] _besselkxv
   @ ~/.julia/packages/BesselK/GkDiL/src/besk.jl:29 [inlined]
 [5] adbesselkxv
   @ ~/.julia/packages/BesselK/GkDiL/src/besk.jl:50 [inlined]

(jl_uQwbSf) pkg> st
Status `/private/var/folders/yq/4p2zwd614y59gszh7y9ypyhh0000gn/T/jl_uQwbSf/Project.toml`
⌃ [432ab697] BesselK v0.5.2
  [26cc04aa] FiniteDifferences v0.12.25
  [f6369f11] ForwardDiff v0.10.33

I would have to look much more closely to see what's going wrong! I am curious to know though why this isn't a check of the type 481 targets. A tolerance check sounds like it would use <, but I presume the branch is before that. Is there valid input to adbesselkxv which would take the other branch?

One way you can hack around this without a dep is something like this:

julia> Dual(0, 1) == 0
false

julia> 0 ≤ Dual(0, 1) ≤ 0
true

(That has a branch in it from &&; could be re-written with & if you care.)

But [edit] this is a pretty weird hack, and I'm not sure it's guaranteed to keep working! To anyone else reading this -- please don't just copy-paste this and move on. Please make an issue explaining what exactly your code was trying to do.

[cgeoga]

Oh, hey, that's a slick workaround! I did look through the PR and saw you explain something similar but it didn't click to me. I'm very happy to just do that.

Since you asked, the actual source of the issue is the branch on that _iszero---when iszero now gives false for a Dual(0, [not zero]), the f0_expansion function was using a function representation that gives a NaN when given zero. So it isn't actually that something is slower to converge---it's that a NaN is slowly propagating through the callstack and automatically failing all the checks.

I should say, though, that I did write this package with a very close eye on ForwardDiff and tuned several things to make it work and dispatch on whether functions are given a Dual or not, but I always used the type bound of v::AbstractFloat to put bounds on when non-AD methods should be called. But as of now all other AD tools give incorrect answers in at least some of the branched methods, so maybe even if there weren't another way to check the old iszero thing I shouldn't be complaining because everybody using that package already has ForwardDiff.jl in their tree anyway.

In any case, I'll go ahead with the double inequality check and I'm a happy camper. Thanks for taking a look and providing guidance!

[tpapp]

It would be great to have a section in the manual that explains the change. While I understand that it was extensively discussed in #481, not all users are prepared to wade through that long discussion.

Practical tips for diagnosis and workarounds would help a lot.

[mcabbott]

If you have a case which needs a workaround, please please make an issue about it. This iszero hack is awful and I don't think we should jump to encouraging such things. It's also not guaranteed not to change.

I do think https://juliadiff.org/ForwardDiff.jl/stable/user/limitations/ should be modernised. Ax_mul_Bx vanished some time ago. "must be written generically enough" should probably be expanded to mean that you must not be relying on tricks aimed at FloatNN behaviour, e.g. this log1p example.