Symbolic differentiation is failed at specific value

[eozw]

Hi,

Symbolic differentiation of power operator (**) is failed at 0.0 .

#I "../packages"
#r "DiffSharp/lib/net46/DiffSharp.dll"
#r "FSharp.Quotations.Evaluator/lib/net40/FSharp.Quotations.Evaluator.dll"

open DiffSharp.Symbolic.Float64

[<ReflectedDefinition>]
module m =
  let y x = x**2.0

let dy = diff <@ m.y @>

// following statement will return nan (0.0 is expected)
dy 0.0
// val it : float = nan

// but this is OK
dy 1.0
// val it : float = 2.0

[dsyme]

@barak @gbaydin The generalized power rule only applies when the x is non-zero. What's the right way to treat this in the variant of symbolic differentiation you're using? thx

[dsyme]

@gbaydin @barak I took a look at fixing this more out of curiosity than anything else. Is adding the conditional check the right fix? https://github.com/DiffSharp/DiffSharp/pull/31/files#diff-419be3dee4710aa72c56268abcf97b05L63

[barak]

The issue is the x**y is continuous and differentiable except at (x,y)=(0,0). But the equation needs to be special-cased for some of the cases x=0 or y=0. Depending on which is "active", to save unnecessary computation.

This is done reasonably carefully in the "FAD" package for Haskell, as I recall.

[barak]

See lines 382 and following in https://github.com/NUIM-BCL/fad/blob/master/Numeric/FAD.hs

[dsyme]

OK, thx

[barak]

Actually looking at that Haskell FAD code again, the division by x on line 391 is wrong when x==0. Needs special cases for that as well, basically x^(n-1) cannot be calculated as (x^n)/x when x==0, so even when x^n is in hand one needs to calculate x^(n-1) from scratch. In the FAD context, this is a particular issue when x is not just a dual number, but an infinite power series, as is often the case there. I suppose the right thing in that context is to first calculate x^(n-1) and then multiply it by x to get x^n.

This is the kind of fiddly numerics that makes me avoid numeric analysis...

[dsyme]

Closing because this is out of date for DiffSHarp 1.x, which is currently planned to be AD only