switch to Bessels.jl for real bessel functions

[oscardssmith]

https://github.com/JuliaMath/Bessels.jl has pure Julia bessel functions for real arguments and orders that are frequently significantly faster than the ones currently used (provided by AMOS). The package was just registered, but the implementations are solid and well tested. Also, on the roadmap for the future are the complex orders and arguments.

[heltonmc]

Thanks for starting this discussion! I am very interested to see everyone's thoughts on this and open to whatever the consensus of the community is on this and would volunteer time to make this happen.

Though, it does appear that discussions like this have frequently occurred here and on discourse about what functions to include in SpecialFunctions.jl, the future of the package, finding maintainers, and reviewing PRs. It does appear like there is a lot of uncertainty around PRs about if they should be included or just their general quality.

Of course, I think Bessel functions are special enough to include here, but there would be a couple of concerns I would have.

1. Testing. Currently Bessels.jl tests many values against the AMOS routines provided by SpecialFunctions.jl. We would have to figure out a different way to test (ideally still having the AMOS routines available in a different package - openspeclib?).
2. Many packages depend on SpecialFunctions.jl which could create some instability for a large part of the julia ecosystem. The code in Bessels.jl to my knowledge is much faster (3x-10x) but it does take a lot longer to compile, is still new and being tested, and being refactored. Inclusion of Bessels.jl would probably over double the source code and increase maintainer burden.
3. Intermixing native julia code and non julia code. One of the advantages I think of Bessels.jl is it has no dependencies and is written entirely in Julia.

All of these can of course be alleviated with some effort, but it may also be worth considering a different approach. I think (speaking very idealistic) it would be great if SpecialFunctions.jl was able to provide only native julia implementations whereas we could put the fortran libraries into some other package (openspeclib.jl?) to test against.

However, it might be best to have SpecialFunctions.jl serve as an umbrella package that could reexport smaller lighter-weight dependencies. For example, I'm sure lots of people might just need a lightweight gamma functions but don't need Bessel or elliptic functions. This also helps with CI and streamlining testing. I think SpecialFunctions.jl could serve as a home and major resource to point users to these packages that were under some sort of purveyorship of the JuliaMath group. A great example is (https://juliadiff.org/) that provides a great overview of the automatic differentiation tools in the ecosystem.

I would then pose that classes of special functions be their own separate packages that consists of groups of similar functions: 1. gamma, inverse gamme, incomplete gamma, polygamma, zeta, 2. bessel, airy, hankel, spherical bessel functions, 3. kelvin, lommel, Anger, Weber, Struve, Scorer, 4. hypergeometric functions, 5. elliptic function, 6. orthogonal polynomials.

Now several of these packages already exist (e.g. HypergeomtricFunctions.jl, Elliptic.jl, ClassicalOrthogonalPolynomials.jl, Bessels.jl) but it would be nice if SpecialFunctions.jl served as a more centralized point and curator of the julia mathematical ecosystem. Of course this takes a lot of work, but I would be happy to volunteer time right now to make that happen.

Edited for conciseness.

[stevengj]

I think it would make a lot of sense to have some subsidiary packages for particular special functions, while keeping this package as an umbrella package (but changing it to simply import from the specialized packages), and encourage downstream dependencies to require the subsidiary package in order to lighten their dependencies.

[devmotion]

I agree, I think moving some classes of special functions to (or keeping them in) separate packages and only reexporting them in SpecialFunctions makes sense. We moved the log/exp functions from StatsFuns to a separate package LogExpFunctions a while ago and it seems development became more active, but the biggest benefit was probably that it made it possible for many packages to depend on and use LogExpFunctions that would never have taken a dependency on StatsFuns (and Rmath).

Generally, I guess it would be simpler though to extract things to separate packages before refactoring or rewriting them (e.g., in native Julia code). That would allow to preserve git history (we did that in LogExpFunctions) and would ensure that no breaking changes are accidentally introduced in SpecialFunctions during the split.

The code in Bessels.jl to my knowledge is much faster (3x-10x)

Are there any benchmarks one can check?

but it does take a lot longer to compile

Can we fix the compile time issues before using Bessels.jl in SpecialFunctions? SpecialFunctions is a central package in the ecosystem (currently 2181 direct and indirect dependents), so it would be very unfortunate to introduce any major compile time or runtime regressions.

[heltonmc]

We moved the log/exp functions from StatsFuns to a separate package LogExpFunctions a while ago and it seems development became more active

This seems like a really great model to follow. I would be very interested in something like that to happen here. Particular for gamma and related functions as it looks like a lot of those already have julia implementations. We could then port the Bessels.gamma(x::Float64) implementation over to there and have Bessels.jl depend on that instead.

Are there any benchmarks one can check?

Of course! I used that qualifier because performance really depends on the specific argument and which branch you fall in. Take for example the Hankel function (which has a lot of different branches in both AMOS and Bessels.jl)

julia> @benchmark Bessels.hankelh1(2.3, x) setup=(x=rand()*100)
BenchmarkTools.Trial: 10000 samples with 940 evaluations.
 Range (min … max):  102.350 ns … 805.496 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     103.856 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   178.015 ns ± 182.483 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █    ▂▃▂               ▁                         ▁▂  ▁ ▁ ▁▁▁  ▁
  █▇▅▅████▄▁▁▃▁▁▁▁▁▁▁▁▁▁▁█▅▄▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁███████████████ █
  102 ns        Histogram: log(frequency) by time        731 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark SpecialFunctions.hankelh1(2.3, x) setup=(x=rand()*100)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.175 μs …   4.062 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.404 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.581 μs ± 366.964 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

        ▁██▆▅                                                  
  ▂▂▂▂▃▆██████▇▅▄▄▄▄▄▃▂▂▂▁▁▁▂▂▂▁▁▁▂▁▁▂▁▂▂▂▂▂▂▂▂▂▃▆▇▆▅▄▃▂▂▂▂▂▂ ▃
  1.18 μs         Histogram: frequency by time         2.5 μs <

 Memory estimate: 32 bytes, allocs estimate: 2.

Of course the longest branch is still shorter than the shortest branch in the amos routine, but it just depends. Looking at an example that doesn't allocate (MPFR routine) for very large orders..

julia> @benchmark SpecialFunctions.bessely(1e6, x) setup=(x=rand()*1e4 + 1e6 )
BenchmarkTools.Trial: 2254 samples with 1 evaluation.
 Range (min … max):  2.189 ms …  2.452 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.222 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.216 ms ± 20.992 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▂▇    ▂                   ▅ █                              
  ▆██▃▃▅▅█▅▃▄▃▃▃▃▃▃▃▃▃▃▃▅▅▆▅▄█▇██▄▃▄▄▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂ ▃
  2.19 ms        Histogram: frequency by time        2.27 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Bessels.bessely(1e6, x) setup=(x=rand()*1e4 + 1e6 )
BenchmarkTools.Trial: 10000 samples with 909 evaluations.
 Range (min … max):  119.637 ns …   2.546 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     120.783 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   238.958 ns ± 394.606 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▆▃                                                           ▁
  ███▄▁▁▁▁▆▇▆▆▆▆▆▆▇▅▆▅▆▆▆▇▅▆▆▇▆▆▆▅▆▆▆▆▆▆▇▆▆▆▆▆▆▆▇▆▆▆▆▅▆▅▆▇▆▆▆▆▆ █
  120 ns        Histogram: log(frequency) by time       2.28 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

So again there is potential for a large amount of performance gains depending on the use case. Maybe a more widely used example with only 2 branches would be besselk of order 0.

julia> @benchmark Bessels.besselk(0, x) setup=(x=rand()*100)
BenchmarkTools.Trial: 10000 samples with 997 evaluations.
 Range (min … max):  19.642 ns … 43.214 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     26.119 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   26.119 ns ±  0.932 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

      ▁                                  ▁▆█▄                 ▁
  ▅▄▃▁█▆▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄████▅▄▅▆▅▅▄▃▄▃▅▃▁▃▄▇ █
  19.6 ns      Histogram: log(frequency) by time      28.9 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark SpecialFunctions.besselk(0, x) setup=(x=rand()*100)
BenchmarkTools.Trial: 10000 samples with 723 evaluations.
 Range (min … max):  173.351 ns …  2.379 μs  ┊ GC (min … max): 0.00% … 89.15%
 Time  (median):     228.043 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   250.504 ns ± 63.841 ns  ┊ GC (mean ± σ):  0.17% ±  1.27%

        ▄█▄                                                     
  ▁▁▁▁▁▁████▇▆▄▃▂▂▂▄▄▄▃▃▂▂▂▂▂▂▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  173 ns          Histogram: frequency by time          502 ns <

 Memory estimate: 16 bytes, allocs estimate: 1.

So there are definitely a lot of potential advantages in runtime for these as well as none of the function allocate.

Can we fix the compile time issues before using Bessels.jl in SpecialFunctions?

I agree with this! I have opened an issue at JuliaMath/Bessels.jl#41. This really only affects two functions (besselj and bessely) and is limited to one function that completely unrolls an asymptotic expansion. The first run for these functions is definitely noticeable (2s) but I will work on this today to improve. Edit: the culprit functions accounts for over half the first run time so hopefully we can find a way to improve this so first run is under 1 second.

[t-bltg]

Cc @juliohm for JuliaPlots/UnicodePlots.jl#291.

[juliohm]

An alternative path would be to develop a separate package named SpecialFuns.jl (or any similar name), which supersedes Bessels.jl and related functionality without dependencies or any promise of backwards compatibility. That way maintenance would be easier for the people working on it, and the switch from SepcialFunctions.jl to SpecialFuns.jl would be trivial downstream. No playing with version numbers, but a simple string replacement s/SpecialFunctions/SpecialFuns/g

Unless SpecialFunctions.jl aims to be a lightweight (zero-dependency) package for special functions, the effort of plugging Bessels..jl into it doesn't seem to pay off in the long-term. Of course I may be wrong.

[ParadaCarleton]

An alternative path would be to develop a separate package named SpecialFuns.jl (or any similar name), which supersedes Bessels.jl and related functionality without dependencies or any promise of backwards compatibility. That way maintenance would be easier for the people working on it, and the switch from SepcialFunctions.jl to SpecialFuns.jl would be trivial downstream. No playing with version numbers, but a simple string replacement s/SpecialFunctions/SpecialFuns/g

Unless SpecialFunctions.jl aims to be a lightweight (zero-dependency) package for special functions, the effort of plugging Bessels..jl into it doesn't seem to pay off in the long-term. Of course I may be wrong.

There would be a huge payoff. Right now, automatic differentiation doesn't work with Bessel functions, which is a huge problem. There's about half a dozen distributions that HMC can't work with in Turing alone. Switching to Bessels.jl would be a huge improvement.

That being said, it looks like compilation time has already been fixed, amazingly enough:

julia> @time using Bessels
  0.021716 seconds (37.68 k allocations: 3.757 MiB)

julia> @time besselix(1, 1)
  0.000001 seconds

This is some amazing work @heltonmc !

[devmotion]

automatic differentiation doesn't work with Bessel functions

SpecialFunctions has contained ChainRules definitions for Bessel functions for quite some time, and some bug fixes in ForwardDiff and ReverseDiff made it possible to finally merge JuliaDiff/DiffRules.jl#79 yesterday which added similar functionality to DiffRules-based AD backends such as ForwardDiff and ReverseDiff. Of course, that won't help if you want to differentiate Bessel functions wrt to the order (see eg #160 (comment)) but basically everything else should work now.

Generally, I also assume the best approach for AD in this case are custom rules, regardless of whether or not the primal function is implemented in native Julia or not.

[heltonmc]

Generally, I also assume the best approach for AD in this case are custom rules, regardless of whether or not the primal function is implemented in native Julia or not.

I agree. That is the best approach when nu=0 as they can be expressed with a single call for nu=1 which is also optimized. I originally thought the oscillating functions besselj and bessely might be less accurate for arbitrary orders due to the subtraction. But after some quick tests it looks like they can be calculated around ~5 ULP error when the Bessel functions are correctly rounded. There was a few cases where I saw Ulp > 40..... One thing to consider is the derivative then requires 2 Bessel function calls. If ForwardDiff could potentially be just as accurate and give the derivative in one call that could be an advantage.

Of course, that won't help if you want to differentiate Bessel functions wrt to the order (see eg #160 (comment)) but basically everything else should work now.

Also, just to mention that Bessels.jl doesn't even allow for dual numbers right now so doesn't work with ForwardDiff. This should be a simple fix to just allow both Float64 and Dual types to propagate. I don't have that much experience with AD so that would be a very nice addition to Bessels.jl if anyone would like to contribute. I should say that we shouldn't also just expect to drop in AD and get accurate derivatives wrt to order. I think besselj and besseli would be the easiest as their power series works for integer and non-integer orders. And I just checked that in the case of besselj the derivative also has stable recurrence... these two functions would be interesting cases to test. I'm less confident for the modified forms (https://github.com/cgeoga/BesselK.jl) as we depend on different algorithms for integer nu and non-integer nu. I'm unsure if derivative information around those branches would be conserved. Though, this would be a very interesting investigation and a great addition!

[heltonmc]

An alternative path would be to develop a separate package named SpecialFuns.jl (or any similar name), which supersedes Bessels.jl and related functionality without dependencies or any promise of backwards compatibility.

I do particularly like this idea.... I think the challenge would be finding maintainers. It would make sense if @devmotion and some of the contributors here could help maintain. But they are also some of the most prolific contributors in the entire Julia ecosystem. I'm not really sure how much time they would have to take on additional things considering it would be a big project. I would like to see the functions here in gamma.jl be pruned into a new repository with no external dependencies. I briefly played around with that trying to maintain git history but couldn't quite get it right. If someone wants to do that I would be happy to maintain it. Bessels.jl has it's own gamma implementation for Float64 that is accurate to <1.5 ULP and I'm working on loggamma for real arguments with that same accuracy. That package would have a lot of nice functionality already that the proposed new special functions could be built on.

[devmotion]

I don't think a separate SpecialFunctions package is the best way to go - I assume since SpecialFunctions is such a core part of the Julia ecosystem in most cases downstream projects and packages will only end up depending on both SpecialFunctions and SpecialFunctions2, im- or explicitly. And if there's a clearly better implementation, why should we keep the other one around instead of letting the whole Julia ecosystem benefit from the improved version? I think the approach outlined in #409 (comment) and #409 (comment) is better.

[devmotion]

And yeah, it's already hard finding maintainers and reviewers for SpecialFunctions, so I assume it would become even more difficult if some focus on one, and some on the other alternative of SpecialFunctions.

[heltonmc]

I think the approach outline in #409 (comment) and #409 (comment) is better.

Yes, sorry I wasn't clear I would be very supportive of that outline and would volunteer to help maintain here and other packages. One thing is figuring out tests as we test against SpecialFunctions (I enjoy having the Amos routines so easily accessible but I'm sure I can set up my tests directly to openspecfun). I would be in favor though of separating out Julia code and having SpecialFunctions.jl bridge the native Julia and external libraries until eventually those are replaced. Ideally, Bessels.jl could depend on a lightweight hypothetical Gamma.jl package that SpecialFunctions could also reexport.

[stevengj]

Onet path forward would be to split more functionality out of SpecialFunctions.jl into subpackages, e.g. GammaFunctions.jl, ErrorFunctions.jl, ExponentialIntegrals.jl, which are loaded and re-exported by SpecialFunctions.jl for backward compatibility. (i.e. the total amount of code wouldn't change).

Then we could encourage downstream packages to rely on only the subpackage(s) that they need rather than on all of SpecialFunctions, and add new special functions in the future to their own subpackage without necessarily including them in SpecialFunctions.jl (so that this package stops growing larger and larger).