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)
NaN
julia> ForwardDiff.derivative(x -> f2(x, 1), 0) # was 0 on v0.10
NaN
julia> ForwardDiff.derivative(x -> f3(x, 1), 0)
-0.5
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
-0.0
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)
NaN
julia> ForwardDiff.derivative(sinc2, 0.0) # was 0.0 on ForwardDiff v0.10
NaN
julia> ForwardDiff.derivative(sinc3, 0.0) # thanks to taylor series
0.0
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)
-0.0
2 Likes
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.