JuliaMath:master ← ChrisRackauckas:dual
opened 03:02PM - 20 Oct 22 UTC
The previous behavior led to https://discourse.julialang.org/t/error-with-forwar…ddiff-in-turing/88633/16 . Dual numbers would just hit an infinite recursion because `float(::Dual)::Dual !<: AbstractFloat`. This should handle more cases where float(T) !<: AbstractFloat too, like TrackedReal.
Though I know that this probably isn't the best way to handle it. In a higher level interface like Integrals.jl, I think it would be better to detect such differentiation by boundary terms and specialize it. It's actually kind of funny and would be a good homework problem. If you put a dual number into the boundary with dual part one, i.e. , so then you integrate (f + df/dt eps) dt and by linearity of the integral it's just fundamental theorem of calculus approximated by GK! However, the argument above assumes that the dual part generated in the rules is exactly 1 (or rather, matches the partials of the boundary term so that it weighs it correctly). The algorithm doesn't exactly work like that, but it ends up being mathematically the same. What happens naturally in the code is that the rule values are multiplied by segment sizes as weights`f(a + (1-x[2i-1])*s)` and it is there that the `x` is "imparted with" the correct dualness to then give the integral of the derivative.
Of course, algorithmically it should probably just see this is the case, integrate f, and then slap f into the dual part. I'm not sure if QuadGK.jl should be mucking with dual numbers and autodiff overloads though, so I'll just do that at a higher level. But it's something to think about here.
The reason to not do a complete overload is because for differentiating with respect to other quantities like parameters, you can fuse the primal and derivative calls rather than having them as two separate integrals, so it's definitely faster to differentiate the algorithm here. Except for handling the boundary.
I don't know, it's pretty cool haha. Such a beautifully convoluted way to approximate f(x).