add exponential scaling support to Bessel function family

[jdrugo]

Would the Julia development team consider supporting the option to return exponentially scaled versions of the Bessel function? This option already exists in the used Amos library through setting the KODE argument, but is currently not used in Julia.

Supporting this option would allow evaluating the Bessel function family for larger arguments, for which the standard return values would be out of bounds. This is, for example, useful for computing various properties of the von Mises distribution (see JuliaStats/Distributions.jl#223 ).

I could imagine it being either implemented as a flag to the Bessel functions, for example as

function besseli(nu::Float64, z::Complex128, scaled::Bool = false)

as in R, or by a set of separate methods

besseliscaled
besselhscaled
...

I could provide the required code if you are happy with this request.

[simonbyrne]

I think this should be in Base. Matlab also makes it available as an option to the bessel* function, but I think a more julian approach would be to create separate functions.

There is already erfcx which scales erfc by exp(x^2): we could call them bessel*x?

[stevengj]

bessel[jyikh]x sounds good to me. These kinds of functions are quite useful in practice, by allowing you to deal separately (analytically) with the exponential factors in algorithms where they arise.

Note that this is needed for besselk too, to avoid underflow, and for bessel{jy} in the complex plane.

[jiahao]

The core Julia team have tried very hard to not introduce more function names, and these functions, while arguably important, would introduce a large number of one-off functions.

I think we need a better nomenclature for these things.

[stevengj]

@jiahao, we could always use besselk(ν,x, scaled=true).

[JeffBezanson]

If these names are reasonably standard it could be ok to add them.

[jdrugo]

@JeffBezanson, both R and MATLAB support scaling through optional arguments rather than separate function names.

[jiahao]

I would much prefer the scaled=true form. The ___x names are unfamiliar to me.

[JeffBezanson]

Sounds good.

[kmsquire]

Is it worth calling the argument "logscaled" yo make it perfectly clear?
To the uninformed (me), scaled implies linear scaling. (Of course, I don't
use these functions.)

On Tuesday, April 29, 2014, Jeff Bezanson [email protected] wrote:

Sounds good.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/6665#issuecomment-41724269
.P

[stevengj]

I don't think scaled implies any particular type of scaling; you will have to read the documentation to find out precisely what it means (but it is easy to guess if you use these kinds of functions).

[StefanKarpinski]

How about log=true as the keyword name?

[jiahao]

logscaled or log suggests that the log of the function (or something) is being computed, which isn't the case.

[StefanKarpinski]

I'm unclear on how this scaling is done, which makes it hard to come up with an apt name.

[jiahao]

The scale factor is there to damp the exponentially growing asymptotics of the special function. See, for example, the comment in amos:besseli and compare that with the asymptotic formula from DLMF. Note that it isn't fully accounting for the zeroth order asymptotics, just the exponential part.

[StefanKarpinski]

Ok, well, then. I will leave it to those better informed on the subject than I am to pick a good keyword.

[jdrugo]

scaled still seems to be the best name, as it does not induce the illusion that all these functions are scaled similarly. If no one disagrees, I will proceed with the implementation using function signatures such as

besselj(nu, z, scaled=true)

Or are there good reasons to make scaled a keyword argument?

[jiahao]

Or are there good reasons to make scaled a keyword argument?

Yes.

[stevengj]

As noted in #4967, there is a substantial performance penalty to using keyword arguments.

[jdrugo]

Speed is also my concern with keyword arguments. In the light of #4967, I would prefer scaled to be an optional positional argument.

Specifically, I suggest the following:

• Leave the current bessel[jyikh] functions (almost) as they are.

• Introduce another set of bessel[jyikh]x (or _bessel[jyikh]x) functions that are not exported, and implement the scaled versions of bessel[jyikh]. As for bessel[jyikh], these functions would be vectorizable with @vectorize_[12]arg.

• Introduce a separate dispatching function, for example,

  besselj{T1 <: Real, T2 <: Number}(nu::Union(T1,AbstractArray{T1}), 
  z::Union(T2,AbstractArray{T2}), scaled::Bool) = 
  scaled ? besseljx(nu, z) : besselj(nu, z)

Up to the last point, this would roughly correspond to the current state of #6688, except that bessel[jyikh]x would not be exported.

This has the following advantages:

• Calls to bessel[jyikh] will be as efficient as they are now.
• Implementation of bessel[jyikh] and bessel[jyikh]x for different argument types is kept separate. This is useful, as there, for example, exist more efficient implementations for besselj(nu, x) when nu is an integer.
• vectorize_[12]arg can still be used for vectorization despite the additional argument.

Possible disadvantages:

• Calling bessel[jyikh] with the scaled argument causes an additional function call.
• Another set of bessel[jyikh]x functions are introduced in the local namespace.

Any thoughts or suggestions on the above?

[jdrugo]

I have now updated #6688 according to the above suggestion. The scaled function versions are _bessel[hijky]x and _airyx, and are vectorized, but not exported. Instead, a f(..., scaled::Bool) dispatch function is introduced for all Bessel functions, that chooses the scaled version if scaled is true.

The scaled argument is introduced for airy(k, z), besselj(nu, z), bessely(nu, z), besseli(nu, z), and besselk(nu, z), all of which can handle vector arguments. As besselh(nu, k, z) is not vectorized, I added the scaled argument to the vectorized hankelh1(nu, z) and hankelh2(nu, z) instead.

The scaled argument has not been added to more specialized variants, such as airyai(z), airyaiprime(z), airybi(z), airybiprime(z), airyprime(z), besselj0(z), besselj1(z), bessely0(z), and bessely1(z). This can be easily changed, if deemed necessarily.

[stevengj]

I'd rather have besseljx(n, x) than besselj(n, x, true).

[jdrugo]

This is unfortunately not for me to decide. I will postpone any further changes of the code until the function signature has been decided upon.

[StefanKarpinski]

Yes, that seems better. If we eventually get faster keywords, we can deprecate there for the keyword version.

On May 5, 2014, at 8:11 AM, "Steven G. Johnson" [email protected] wrote:

I'd rather have besseljx(n, x) than besselj(n, x, true).

—
Reply to this email directly or view it on GitHub.

[jdrugo]

This is obviously not a high priority issue, but I would still like to make some progress towards finalizing the code. However, at this point I cannot make any further progress until the Julia developers have agreed on how the functions are supposed to be called.

To summarize, the three main options seem to be

besselj(nu, z, scaled::Bool = false)
besselj(nu, z; scaled::Bool = false)
besseljx(nu, z)

where the first (positional argument) option seems to preferred due to

• not introducing any new function names (as raised by @jiahao), and
• not being slower than the current implementation.

Option two (keyword argument)

• might be preferable from the syntactic point-of-view (besselj(nu, z, scaled=true) being more explicit than besselj(nu, z, true)), but
• comes at the cost of being slower (as raised by @stevengj).

The third option is preferred by @stevengj, and might be admissible if the function names were sufficiently standard (@JeffBezanson), which they unfortunately aren't.

Which one shall it be?

[staticfloat]

@jdrugo thank you for your patience working through this issue. It looks like @stevengj and @StefanKarpinski are both in agreement for option 3 as you've listed it up there. Until we get faster keyword arguments, I think we're just going to have to live with the suffixes on function names.

[jdrugo]

@staticfloat OK, I have updated #6688 to use suffixes on function names.

Again, I did not export besselhx, but rather hankelh1x and hankelh2x to facilitate vectorization by @vectorize_2arg.