add derivatives for exponential integrals

[stevengj]

Just added to SpecialFunctions.jl: JuliaMath/SpecialFunctions.jl#236

The derivative with respect to x is a simple recurrence: https://en.wikipedia.org/wiki/Exponential_integral#Derivatives

[oxinabox]

If I understand @gxyd this is blocked by SpecialFunctions not having defined a meijerg function

[stevengj]

@oxinabox, that is only if you want derivatives with respect to order. For now, I would treat it like the Bessel functions, where you can only differentiate with respect to the second argument and not with respect to the order at the moment.

[oxinabox]

Oh, we should not be doing that, that is not good.
I was not aware we did that.
I thought we only were bailing out on things that were not differentiable due to requiring integers

[stevengj]

Oh, we should not be doing that, that is not good.

The perfect seems like the enemy of the good here? There are lots of applications for differentiating with respect to the argument alone, so it is useful to be able to do that automatically even if you can't handle differentiating the order.

For both Bessel functions and exponential integrals, in most applications the order will be a constant (usually an integer), not something that you would consider continuous perturbations of.

(Though it's theoretically possible to differentiate both Bessel functions and exponential integrals with respect to the order, and useful in rare applications, these both require special functions that are not currently implemented efficiently in any finite-precision open-source code that I'm aware of.)

[gxyd]

For reference: Here's the besselj function derivative (for both nu and z): https://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/20/ShowAll.html

For exponential integral: https://www.wolframalpha.com/input/?i=diff%28expint%28nu%2C+z%29%2C+nu%29 (here G is MeijerG function)

Though it's theoretically possible to differentiate both Bessel functions and exponential integrals with respect to the order, and useful in rare applications

Agreed.

[gxyd]

We might want to document the meaning of NaN to refer to whether it means "derivative doesn't exist" or "derivative can't be represented in the classical sense".

[sethaxen]

Oh, we should not be doing that, that is not good.
I was not aware we did that.
I thought we only were bailing out on things that were not differentiable due to requiring integers

We also do this for getproperty(::SVD, ::Symbol):

ChainRules.jl/src/rulesets/LinearAlgebra/factorization.jl

Lines 28 to 29 in 2eb203f

	# TODO: https://github.com/JuliaDiff/ChainRules.jl/issues/106
	throw(ArgumentError("Vt is unsupported; use V and transpose the result"))

.
I in general also don't think we should do this except in cases where no current Julia AD can AD through the function. Then there's no harm in writing a rule that ignores a differential that will be almost never used, as long as we can throw a sensible error should someone try. And of course we can always remove the rule if some future AD can handle the primal function. My two cents.

[oxinabox]

We also shouldn't be doning that for SVD, there is a (long standing) issue open about that. #106

We might want to document the meaning of NaN to refer to whether it means "derivative doesn't exist" or "derivative can't be represented in the classical sense".

NaN doesn't mean that.
I think the reason SpecialFunction's rules are like that is because they were ported from DiffRules very early in ChainRule's history, and in DiffRules NaN meant something like that.
ChainRuleCore.DoesNotExist() means that there is no meaningful (co-)tangent space associated with the primal space (e.g. for integer or Boolean valued arguments, that are not embedded in a floating point space).

I in general also don't think we should do this except in cases where no current Julia AD can AD through the function

I agree with this.
Is that the case here?

If so we we should do (kind of like SVD does) @thunk error("Not Implemented")
and also should do so for the besselj etc.

[stevengj]

If so we we should do (kind of like SVD does) @thunk error("Not Implemented")
and also should do so for the besselj etc.

Or you could even have a predefined ChainRulesCore.unimplemented = @thunk error("unimplemented") return value so that you are not creating a different Thunk on each call? (Or better yet, define a new type ChainRulesCore.Unimplemented <: AbstractThunk, which will provide more information to the caller.)

[oxinabox]

Nice think about @thunk in the right place is that it will give the the actual location in code where it was used.
Which is often a problem with AD, tracking down where some invalid value came from.
On that basis have actually considered replacing more things with macros that do the same trick, e.g. replaceing DoesNotExist.
Hopefully we don't have so many things that are not implemented that this is a problem.

Thunks are cheap to create I hope, they live in the stack.

[gxyd]

Or you could even have a predefined ChainRulesCore.unimplemented = @thunk error("unimplemented")

There was something similar #28 sometime back, though probably not exactly same.

[devmotion]

For reference: Here's the besselj function derivative (for both nu and z): functions.wolfram.com/Bessel-TypeFunctions/BesselJ/20/ShowAll.html

Related issue: #208

[stevengj]

Should be transferred to SpecialFunctions (#319).

[stevengj]

Note that Zygote seems to be capable of accurately differentiating our expint implementation, both with respect x and with respect to the order! At least, it worked for the few points I tried. (ForwardDiff, in contrast, fails.)

Differentiation with respect to x, however, seems to be about 7000x slower than using the analytical derivative formula.

julia> using Zygote, SpecialFunctions, BenchmarkTools

julia> f(x) = expint(3,x)

julia> fp(x) = -expint(2,x) # analytical derivative

julia> f'(3.7)
-0.004566575240288366

julia> fp(3.7)
-0.00456657524028863

julia> @btime f'(3.7);
  825.321 μs (5983 allocations: 173.81 KiB)

julia> @btime fp(3.7);
  118.380 ns (0 allocations: 0 bytes)

julia> g(ν) = expint(ν,4.2)

julia> g'(4)    # matches high-order finite-difference result to 12 digits
-0.0002319482036089421

It would be nice if there were some way to declare a chain rule for only one of the arguments, and let Zygote handle the other one if needed.

[oxinabox]

It would be nice if there were some way to declare a chain rule for only one of the arguments, and let Zygote handle the other one if needed

I wonder if we could solve this with #68 if we can do that without triggering an infinite loop.

[nickrobinson251]

Xref JuliaDiff/ChainRulesCore.jl#68

[mzgubic]

Can anyone transfer this to SpecialFunctions.jl? We are deleting the special functions rules in the next breaking release

[oxinabox]

It is not possible to transfer issues between orgs.
I will close this and open a new one