Automatic differentiation of function with special points

I am using ForwardDiff.jl to compute derivatives of functions which often contain functions of the form:

f(x,y) = x/[exp(x/y)-1]

which one can show using L’Hôpital’s rule has the value f(0,y) = y. One way to construct such a function is:

f(x, y) = (x == 0.0 ? y : x/(exp(x/y)-1) )

This is continuous with respect to x, but ForwardDiff understandably gives an x-derivative of zero at x=0. The actual derivative at x=0 should be f’(0,y) = -1/2. If I only required the first derivative, then I could put in a workaround to take this into account (e.g. x == 0.0 ? y-x/2 : ... ). However, I also intend to use the second and third derivatives wrt x and I’d like to know how to generalise to similar functions.

Is there a way to do this by defining custom derivatives? Or put another way, why does ForwardDiff work on the sinc function when it also contains a special point similar to this at x=0?

ForwardDiff should ignore branches like x == 0.0 as these often lead to wrong answers (like here). For various reasons this is only on master branch right now. It leads to NaN instead of a wrong answer.

One way around this is to splice in a Taylor expansion on a finite region. A very quick attempt (without thinking carefully about how it will behave at very small nonzero x) looks like this:

julia> f1(x, y) = x/(exp(x/y)-1);
julia> f2(x, y) = x == 0.0 ? y : x/(exp(x/y)-1);
julia> f3(x, y) = abs2(x)>eps() ? x/(exp(x/y)-1) : y - x/2 + x^2/(12*y) - x^4/(720*y^3);

julia> using ForwardDiff  # master, i.e. v0.11

julia> ForwardDiff.derivative(x -> f1(x, 1), 0)

julia> ForwardDiff.derivative(x -> f2(x, 1), 0)  # was 0 on v0.10

julia> ForwardDiff.derivative(x -> f3(x, 1), 0)

This has an explicit rule in DiffRules.jl. Note that ForwardDiff reads these rules once while loading, so you can’t easily add rules for your functions this way. You can (as mentioned) add methods like sinc(x::Dual) :

julia> ForwardDiff.derivative(sinc, 0.0)  # this calls cosc

julia> @which sinc(ForwardDiff.Dual(0, 1))
sinc(d::Dual{T}) where T
     @ ForwardDiff ~/.julia/dev/ForwardDiff/src/dual.jl:238

julia> sinc1(x) = sin(x)/x;
julia> sinc2(x) = x==0 ? one(x) : sin(x)/x;
julia> sinc3(x) = abs2(x) > eps() ? sin(x)/x : 1 - x^2/6 + x^4/120;

julia> ForwardDiff.derivative(sinc1, 0.0)

julia> ForwardDiff.derivative(sinc2, 0.0)  # was 0.0 on ForwardDiff v0.10

julia> ForwardDiff.derivative(sinc3, 0.0)  # thanks to taylor series

julia> using ForwardDiff: Dual, value, partials

julia> sinc4(x::Dual{Z}) where Z = Dual{Z}(value(x), cosc(value(x)) * partials(x));

julia> ForwardDiff.derivative(sinc4, 0.0)

This is mostly unrelated to the ForwardDiff question, but for things like this expm1 is very much your friend.

1 Like

This is interesting, thank you I have used a Taylor series approach to resolve this issue for now. Is there any particular advantage in defining function(x::Dual{Z}) where ... as opposed to/in addition to function(x) = abs(x) > eps() ? ...? Ultimately the definition of cosc(x) still requires a value to be defined at x=0, right?

Indeed it is, thanks for reminding me of this.