Special half-integer branch for besselk

[cgeoga]

Here is a PR that closes the discussion of #25. I've tried to format things in a way that is consistent with the rest of the code in terms of four space tabs, not lining up the equals signs, putting the two consts into the appropriate file, and adding tests. But I'm sure there's something I missed so please let me know what should be fixed.

I recognize that having another check of isinteger(nu-1/2) in every Bessels.besselk call is probably not super exciting to you guys, but at least for stat people it is a branch that will get hit often and the speedup is pretty nice.

After getting the hang of making these PRs with this easy one I'll follow up and make another one for the asymptotic expansion to use in the case of small nu and large x, which I think will have some convincing speed gains.

[cgeoga]

I see that one test is failing, but it doesn't have to do with the code in this PR and all tests relevant to this new code pass. I'm not really sure how to handle that.

[heltonmc]

Thanks for contributing this!!! A few general ideas on how we want to handle this. Because this is so linked to the modified spherical bessel functions it would also be nice to have those defined. I've already given routines for the regular spherical bessel functions (https://github.com/JuliaMath/Bessels.jl/blob/master/src/sphericalbessel.jl) so having these linked would be nice (but probably in a separate PR).

I guess having to call something like sphericalbesselk to get specialized half order routines is not ideal? For your applications it is typically standard to use the half integer of besselk? One thing is I currently don't have a branch for half integer order for besselj and bessely I just assume that if you want those you can call sphericalbesselj instead. But either way it is probably good to be consistent (I'm open to either).

One thing is for very large orders it's probably more efficient to use the asymptotic expansion? I have defined something similar at (

Bessels.jl/src/sphericalbessel.jl

Lines 136 to 145 in 0302204

	function sphericalbessely_positive_args(nu::Integer, x::T) where T
	if besseljy_debye_cutoff(nu, x)
	# for very large orders use expansion nu >> x to avoid overflow in recurrence
	return SQPIO2(T) * besseljy_debye(nu + one(T)/2, x)[2] / sqrt(x)
	elseif nu < 250
	return sphericalbessely_forward_recurrence(nu, x)[1]
	else
	return sphericalbessely_generic(nu, x)
	end
	end

). The asymptotic expansion is so fast (hopefully we can figure out the other PR on vectorization) that we might want to consider that for orders less than 100 or so...

[cgeoga]

Good points. This function I'm submitting literally is sphericalbesselk with an offset of 1/2 on the order. Would you like me to just rename it (and change the shift as appropriate) and move it to the spherical file? It looks like you don't have a sphericalbesselk yet, so that could do two birds with one stone, and I don't think anybody will care what the function being called is name so long as they're getting the fast branch. I certainly don't have strong feelings about it.

Also, a good call about checking for v sufficiently large. That does make sense to do and I can do a few back of the envelope timings and pick a cutoff for that.

[heltonmc]

I think this definitely gives a solution for the modified spherical bessel function (only for integer orders) so we should probably define that as well. But I think I'm generally in favor of having that as a separate file (modifiedsphericalbessel.jl) and we can have the specialized half integer routines there.

I think the general question is should we dispatch for half integer orders within the besselk routine.... I think I'm open to whatever users would reasonably expect to happen. We already have a branch for integer orders

Bessels.jl/src/besselk.jl

Line 233 in 0302204

isinteger(nu) && return besselk_up_recurrence(x, besselk1(x), besselk0(x), 1, nu)[1]

so should that branch just be integer or half integer orders and then the forward recurrence can use the starting values either with besselk0 or besselkhalf.

Or would users look to the sphericalbesselk routines if they are using half integers of besselk? 🤷‍♂️

Also it does look like that test failure is related. It is failing at (v, x) = (51.5, 97.85) which should be dispatching to this new routine? Could you add some more higher orders into your test cases?

[heltonmc]

I think testing for larger orders and for more x values will be beneficial. These tests are so cheap and if we add other branches it will be nice to test around and after the cutoff.

[cgeoga]

Oof, you're right about the tests. Apologies for that.

With regard to whether the branch should be in besselk or not, I think to some degree that's up to you. But if you don't put it in there, I assume downstream users (like me) will just create their own wrappers and add it in ourselves. The statisticians who need this will hopefully use BesselK.jl to get the auto-derivatives w.r.t. v, and so I would be perfectly happy to just add the branch in my code and keep your implementation that little bit leaner. I do intend, by the way, to move to your besselk from SpecialFunctions.besselk when I get a spare minute.

With regard to the comparison with the whole integer order, you're right that it's the same. Maybe better to re-use the code a bit more efficiently and use the existing recurrence function. I can look into that.

Also, about modified spherical Bessels, I have never seen a non-integer order. NIST 10.47, the section on them, explicitly comments on the order only being an integer. Are there fractional-order spherical ones to worry about? If not, I would still be happy to move this function to a file and name it the spherical name (as well as the other checks you mention for v, adding more tests, etc).

I'm covered in meetings until the afternoon, but I won't let this PR languish, so will be on it soon.

[oscardssmith]

this will automatically inline.

[cgeoga]

Good call. Will address in next commit.

[codecov]

Codecov Report

Merging #46 (a3bfb4c) into master (0302204) will decrease coverage by 0.03%.
The diff coverage is 92.00%.

@@            Coverage Diff             @@
##           master      #46      +/-   ##
==========================================
- Coverage   93.79%   93.76%   -0.04%     
==========================================
  Files          17       18       +1     
  Lines        1467     1491      +24     
==========================================
+ Hits         1376     1398      +22     
- Misses         91       93       +2     

Impacted Files	Coverage Δ
src/modifiedsphericalbessel.jl	90.00% <90.00%> (ø)
src/besselk.jl	98.80% <100.00%> (+0.01%)	⬆️
src/constants.jl	79.31% <100.00%> (+0.73%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[heltonmc]

I'm wondering if there's a better way to combine the large order checks? As in we check for large orders here besselk_vhalfint_check and then after we check for large orders in besselik_debye_cutoff so we are kind of doubling the checks.

It might make sense to have the debye_cutoff check come before, however, we don't want to dispatch for large arguments (x) though so that would require some care.

[cgeoga]

Not crazy. I'll tinker with this and try to get rid of the extra check.

[heltonmc]

One way is to just have the large order check within the besselk_vhalfint which would avoid a check in the top level routine of besselk and also make besselk_vhalfint more generic (right now it is specialized to always do forward recurrence). So potentially if you want to call Bessels.besselk_vhalfint you'd have the full routine?

Just some thoughts on this if we want to keep the half order check in the top level routine.

[heltonmc]

Do you know why that test seemed to fail? I notice it is passing now but it's because it's not dispatching to that method anymore.

Funny enough - mathematica doesn't even give functions for SphericalBesselK or SphericalBesselI just the j and y variants. But yes there are non-integer orders for these but they lose their meaning and don't have any special expressions. But we still should just fall back to computing besselk for any order at that point, but I doubt anyone would call these really for noninteger orders.

I think the path forward (the code in general looks good) to merge this. Is to

1. Move the half integer to a modifiedsphericalbessel.jl file
2. The syntax should probably follow

   Bessels.jl/src/sphericalbessel.jl

   Line 113 in 0302204

   sphericalbessely(nu::Real, x::Real) = _sphericalbessely(nu, float(x))

   that from bessely (as forward recurrence is always stable for both of these.
3. Decide if we want to have the half order check in main routine. Though users could always call Bessels.besselk_vhalfint if they wanted to and just need to document that this function should be stable.

[cgeoga]

Makes sense. I've just pushed another commit that moves things over. In the end if it's going to be its own function I re-implemented the spherical modified k function with the simpler recursion. It does have an extra branch in the beginning which I know you guys won't like, but since the symmetry is broken I don't see a way around it. I'd love to be wrong though.

I'm confused about the CI failing, because all tests pass on my local branch.

Also, with regard to whether the branch goes into the main besselk or not, I feel like I've kind of made my argument above. I also think it is a funny line to draw to say that it's worth it to have a special branch for integer orders and not half-integer orders. But again, if you decide not to put it in then downstream folks like me will just add the control flow in our own code which is easy, so it doesn't seem like a big deal either way to me. Maybe you should just decide and let me know what you want?

With regard to why the tests were failing for a few Float32s, I actually don't have a great explanation. I looked into it using ArbNumerics and the current values are more accurate. But only in some reasonably specific parts of the domain. So kind of a stumper to me.

[heltonmc]

Really thanks for making these changes (I know it's a pain). I have something that will occupy me the rest of the day but one thought looking at the code now is it looks like it will give wrong results for negative orders if you call the routine directly.

In the tests you put in you are calling the general besselk routine which filters the arguments first to only call the routine for positive orders and arguments (so this will give correct results as it makes appropriate transformation). But what if you call sphericalbesselk directly for negative orders? On a quick inspection this looks like it would give you wrong results as we are only iterating in the forward direction towards positive numbers? Could you check on this or add tests? Just a general intuition without testing them carefully.

Edit: Of course now there's a bunch of method errors 🤷‍♂️

[cgeoga]

It's my pleasure to work this out. I spend all day fitting Matern functions and so I'll probably benefit the most of all the users of this package.

Good catch on the negative order. I added a special check for that in the manual integer order routines where the negative sign wouldn't otherwise be handled correctly and also added some tests.

Also, apologies for the method errors---my fault from the last overly quick commit.

[heltonmc]

What do we want to do for negative arguments? Do you need these routines for complex arguments? We should probably throw domain error for now instead of giving a method error.

[cgeoga]

Argh, you're right. Apologies, I'm just not in the headspace of implementing this as a full function that will get called on its own and stuff. I don't need this for complex args and I'd be happy to put in the same check/error you use on besselk.

[heltonmc]

No problem! Does NaN propagate through this as well? But I think the general structure of besselk for this is good. And we can add complex methods as we get to them.

[cgeoga]

Another good catch. Just pushed something up that I think addresses both, although yet another branch due to an isnan to avoid just getting the domain error.

[heltonmc]

Ok I've had time to look through this. I've made some changes but had some git issues which created a new branch. I'll merge this then fix those changes after.

There are a couple of things I tried to fix.

1. Returned NaN for infinity arguments (when order was integer)
2. With that cutoff there seemed like a small subset of numbers that would dispatch from besselk to sphericalbesselk then back to besselk. I think that should be avoided.
3. If you fix 2 by using same cutoffs you can just call the sphericalbesselk_int routine directly to avoid additional checks etc.
4. The constants are already defined in math_constants.jl so those are duplicates.
5. There is a premature overflow possibility in the forward recurrence because of doing the scaling by exp(-x) at the very end instead of scaling these first.
6. I've also added some docs and spacing issues.
7. Added some more tests for these cases
8. Handle Float16

I'll try to see if I can merge this into here real quick (I thought I was directly editing this.... but I guess not)

[heltonmc]

Thanks for contributing this @cgeoga !!! Looking forward to the large argument expansions 😉

[cgeoga]

Amazing, thanks for all the improvements! Preparing the large argument expansion one now.