# Modified Bessel Function of the first kind of order zero

To calculate the Kaiser window, I would like to implement the Modified Bessel Function of the first kind of order zero (I0). To make this work both with `f32` and `f64`, I would like to use the `num-traits` crate.

This is my code so far:

``````use num_traits::{Float, NumAssignOps};

#[allow(non_snake_case)]
pub fn bessel_I0<Flt>(x: Flt) -> Flt
where
Flt: Float + NumAssignOps,
{
if !(x >= Flt::zero()) {
return Flt::nan();
}
let base = x * x / Flt::from(4).unwrap();
let mut power = Flt::one();
let mut coefficient = Flt::one();
let mut sum = Flt::one();
for j in 1.. {
power *= base;
if !power.is_finite() {
return Flt::infinity();
}
coefficient /= Flt::from(j * j).unwrap();
let addend = coefficient * power;
break;
}
}
sum
}
``````

(Playground is currently buggy, so I'll provide a Playground link later)

I would like to know if my implementation is idiomatic or if it is bad style to implement it as I did. I particularly struggle a bit with the constants, such as `4`, which I have to write like `Flt::from(4).unwrap()`, which uses `num_traits::cast::NumCast::from`.

Another thing that bugs me a bit is the `num_traits::NumAssignOps` trait. If I don't add it, I can't use the `*=` and `+=` operators. But adding this bound makes all generic functions which use this function require that bound as well (which bloats up function signatures a bit). Why is `NumAssignOps` a separate trait? Does it have to do with certain implementations of floats which are immutable by design? I think they could somehow be wrapped to support `+=`. Afterall, `Copy` is a supertrait of `num_traits::float::Float`.

Also happy on any comments on numerical stability (if there are any).

The code seems fine except for the backwards check for `x < 0` at the very beginning, as well as the the zero check of `addend` which could be moved up one line, directly checking whether `coefficient` is zero.

That's probably the only sensible way to create numeric constants of a generic type.

There's no such thing as an immutable type in Rust because mutability is a property of a binding, not a property of a type. I'm not an authority on `num-traits`, but it looks like the crate tries to follow the principle of composition/single responsibility by breaking up requirements into separate traits, instead of shoving everything in one giant trait.

1 Like

Oh, that's indeed a mistake (don't even remember why I added it; I think because I looked at the graph and thought I0(x) is undefined for negative x, but that's wrong). Thanks for catching it.

Oh, right; at that point it's known that `power` is finite.

Here's the fixed function:

``````use num_traits::{Float, NumAssignOps};

#[allow(non_snake_case)]
pub fn bessel_I0<Flt>(x: Flt) -> Flt
where
Flt: Float + NumAssignOps,
{
let base = x * x / Flt::from(4).unwrap();
let mut power = Flt::one();
let mut coefficient = Flt::one();
let mut sum = Flt::one();
for j in 1.. {
power *= base;
if !power.is_finite() {
return Flt::infinity();
}
coefficient /= Flt::from(j * j).unwrap();
if coefficient.is_zero() {
break;
}
let addend = coefficient * power;
}
sum
}
``````

I found `numeric_literals`, but not sure if it's worth adding a dependency for a single constant (other than zero and one).

Hmmm, right. `f32` and `f64` don't come with "inner mutability" either, but you can reassign a variable of type `f32` or `f64`. Sorry for my confusion.

But `Float` has a bunch of other supertraits as well:

Most of these are inherent to floating point numbers, but `Copy` isn't, I think? Also, `NumOps` is a supertrait, but `NumAssignOps` isn't. I wonder why.

My opinion is that generic functions working with primitive integer and floating-point types are more trouble than they're worth. It is much easier to split your function into two separate concrete functions `bessel_I0_f32` and `bessel_I0_f64`. You don't need to bother with the absurd generic bounds that Rust requires, with by-reference and by-value implementations of operations, with possibly missing at the generic level functionality of concrete types (new methods on primitives are regularly added, are they also added to `num-traits`?). The precision of floating-point types may also affect your algorithms, and you can use the best possible algorithm in each case without trying to shoehorn it into a generic mold.

The implementation within your function really looks simple enough that it can be duplicated without any trouble. If you use later a more complex evaluation method, you can try extracting separate parts of it into private generic functions with more specific and simple trait bounds, or you can turn the body of the function into a macro (which is, admittedly, also a bit painful to work with).

By the way, I don't see a reason to use `coefficient` and `power` separately. You can directly compute the `addend` on the next step via the recurrence, and you will save a multiplication. I also don't see a good reason to check `addend != 0`. This is really only possible if `x == 0`, which is easier to handle separately at the start of the function. Either that, or you expect the loop run long enough to underflow the floating point type, which doesn't make any sense. It is way longer than required to compute the series up to any reasonable precision, and many of the terms in the series won't even affect the sum, since they will be smaller than the minimal representable difference from `1`.

Also, if this is not just a learning problem and you actually need to compute `bessel_I0` at many points, consider using a faster converging series than the naive Taylor expansion

2 Likes

What do you mean exactly?

You mean that certain numeric types should be used by-reference? (`RefNum`, `NumRef`) Does that also apply to working with `num_traits::float::Float`? I see that `Copy` is a supertrait of `Float`, though I would guess you always work by value? Or is that a false conclusion?

Sounds like you'd worry that 3rd party crates could easily be outdated in that matter?

I wonder how other crates that perform numeric operations deal with this? I looked at `rustfft`, and it provides an own trait `FftNum` for the API (that doesn't say anything about the internal implementation though, which I didn't look at).

I think providing two different functions for the API isn't so nice. But apparently `rustfft` didn't use `num_traits::float::Float` for some reason, but their own trait (that's only implemented for `f32` and `f64`).

For such a small function, yes, but I will add more functions that will use this function. Should I duplicate them too? Eventually I'd end up with all code duplicated, I guess. That doesn't seem nice.

I thought this might help with numerical stability. Keeing the `coefficient` in a separate variable will ensure that `coefficient` will reach `Flt::zero()` quickly (in which case the iteration is stopped). I could consider putting `coefficient` and `power` into one variable (`addend`), which would then look like this:

``````use num_traits::Float;

#[allow(non_snake_case)]
pub fn bessel_I0<Flt: Float>(x: Flt) -> Flt {
let base = x * x / Flt::from(4).unwrap();
let mut sum = Flt::one();
for j in 1.. {
break;
}
return Flt::infinity();
}
}
sum
}
``````

The problem here is that I didn't overlook how big `j` can practically grow. `j * j` needs to fit into an `i32` to not cause overflow. Is it possible to prove that this algorithm is halting eventually (or before `j * j` overflows) for every input? Maybe this is more a theoretical problem and wouldn't cause any problem in practice.

But I now see you propose to provide a different termination condition anyway:

So I guess it's better to check if the result doesn't change anymore? Like this:

``````use num_traits::Float;

#[allow(non_snake_case)]
pub fn bessel_I0<Flt: Float>(x: Flt) -> Flt {
let base = x * x / Flt::from(4).unwrap();
let mut sum = Flt::one();
for j in 1.. {
let old = sum;
if sum == old || !sum.is_finite() {
break;
}
}
sum
}
``````

Regarding `+=` and `*=`, I wonder if there is any downside in writing `sum = sum + addend` instead of `sum += addend`, beside a more verbose syntax. It would let me get rid of the `NumAssignOps` bound. Are there any floating point number implementations that wouldn't want to implement `+=` or `*=` ?

I'm not an expert in numerical computation. Do you have any helpful links in that matter?

I mostly did this to learn how to use `num-traits` (and to see how well it works or if there is any trouble in writing generic numerical functions). I could also use an existing trait for my problem, but not sure if it's worth adding a big dependency.

One more question: What's your opinion on whether it's reasonable to violate the snake_case convention here by making the `I` in `bessel_I0` capital?

In `argmin` I created a new trait `ArgminFloat` with all the trait bounds that I found reasonable for my use case and implemented this trait for all types which fulfill these trait bounds. I think this is similar to the approach in rustfft. You can find the implementation here. A simple `float!(..)` macro avoids having to write `F::from(...).unwrap()` everywhere (which still isn't great, but at least it leads to fewer distracting unwraps).

In `linfa` they did it in a similar way, but they implemented their `Float` trait on `f32` and `f64` explicitly and they added a convenient `cast` method.

I'm certainly not an authority on this topic, but I have found this to work pretty well.

1 Like

So practically a "trait alias" (see also nightly feature) for all the traits you need. It seems to be common practice then, that each provider of floating point functions will provide their own trait, which specifies which functionality for floating point types is needed for a particular library.

I could create one for `Float + NumAssignOps`, I guess. But maybe desugaring `x += y` to `x = x + y` is not a big problem in my case. However, `+=` (`AddAssign`) works by reference on the left hand side and might thus be faster for some implementations of `Float`s, I guess? That would be a good reason to provide my own trait. I think there might be an advantage to implement the trait on `f32` and `f64` explicitly (and maybe even make it a sealed trait?). If you later need to add more supertraits to your own trait, this would break semver compatibility. See also what @afetisov said above:

I didn't consider semver compatibility in that matter.

Edit: Maybe a reasonable solution would be to simply document that new bounds might be added to whichever `Float` trait (alias) is being used, and that this doesn't break semver compatibility (i.e. you can only rely on `f32` and `f64` implementing it).

I see three and a half ways to go (which may be different for the API and internal/private functions):

• Not use generics and only work on `f32` and `f64` directly
• either by duplicating code manually
• or through macros;
• Using a trait alias (or trait with blanket implementation on stable Rust) that combines something like `num_traits::float::Float` with any other traits that are needed for a particular use case;
• Using an own (possibly sealed?) trait that's directly implemented for `f32` and `f64`.

Hmmmm… lots of possible trouble. Not sure what I should do.

Side question:

Does it also contain a (multivariate) Estimation of Distribution Algorithm as one possible strategy, and what do you think about those algorithms?

If I understand it right, the macro requires the type to be named `F` (source). I think that should be documented (currently doesn't seem to be documented). Maybe I'll create a private macro for my purpose. Something like:

``````use num_traits::Float;

macro_rules! flt {
(\$x:expr) => {
Flt::from(\$x).unwrap()
};
}

#[allow(non_snake_case)]
pub fn bessel_I0<Flt: Float>(x: Flt) -> Flt {
let base = x * x / flt!(4);
let mut sum = Flt::one();
for j in 1.. {
let old = sum;
if sum == old || !sum.is_finite() {
break;
}
}
sum
}

fn main() {
for x in [-1.0, 0.0, 0.5, 1.0, 2.0, 13.4, 500.0, 1000.0] {
println!("I_0({}) = {}", x, bessel_I0(x));
}
}
``````

Taking the next step, and calculating the actual Kaiser window, I wrote this:

``````pub fn relative_kaiser_beta<Flt: Float>(beta: Flt, x: Flt) -> Flt {
bessel_I0(beta * Flt::sqrt(Flt::one() - x * x))
}

pub fn kaiser_fn_with_beta<Flt: Float>(beta: Flt) -> impl Fn(Flt) -> Flt {
let scale = bessel_I0(beta).recip();
move |x| relative_kaiser_beta(beta, x) * scale
}

fn main() {
let kaiser = kaiser_fn_with_beta(2.0);
for i in -8..=8 {
let x = 0.0 + i as f64 / 8 as f64;
println!("kaiser({}) = {}", x, kaiser(x));
}
}
``````

`relative_kaiser_beta` calculates the Kaiser window with an unknown scaling factor (which in some cases is okay because you perform some sort of normalization anyway).

`kaiser_fn_with_beta` returns a closure that returns values of a normalized (peak at `1.0`) Kaiser window with a given `beta`. I used Currying here, so that I can write

``````kaiser_fn_with_beta(beta)(x)
``````

but also store any intermediate values (that are the same for different `x`) as part of the closure:

``````let kaiser = kaiser_fn_with_beta(beta);
let y1 = kaiser(x1);
let y2 = kaiser(x2);
/* … */
``````

I.e. this doesn't calculate `scale` more than once.

Is Currying a reasonable technique to use here? Or should I rather have a `struct` with a method?

Another problem: If I increase `beta` too much, I run into numerical problems: Playground. Again, I'm not experienced well-enough with numerics (or math in general) to come up with a way to work around it. Maybe I need to find a simplified representation of the scaled result (perhaps using a computer algebra software). In practice this doesn't matter, as `beta` should be reasonably small in most cases.

P.S.: Naïvely feeding this into WolframAlpha doesn't work.

P.S: Naïvely feeding this into WolframAlpha doesn't seem to work. I still get a result where a factor in the numerator grows quickly: β0, β2, β4, β6, …

That's true, it assumes `F` for the type.

I actually didn't consider sermver compatibility either. Implementing the trait explicitly of course solves that problem.

argmin currently does not contain such algorithms but they do sound very interesting and are certainly in the scope of the library. I'm happy to accept a PR if you (or anyone else) is interested in implementing Not sure if implementing the trait explicitly solves the issue. Looking at your earlier example of `linfa::Float`, this is a public trait, and I believe it's non-sealed (though not sure if all its supertraits are really in public modules). I could use `linfa` with my floating point type for which I provide a `linfa::Float` implementation. Now let's assume `std` extends `f32` or `f64` with a new method. That method is added to `linfa::Float` as well. This would mean `linfa` would (theoretically) do a major version bump (see Rust reference). Adding an inherent item in contrast would be considered only potentially breaking (see Rust reference).

So the only solution would be to have a sealed trait or to specify that users of the crate may not rely on the trait being stable between new minor or patch versions (but only rely on being implemented by `f32` and `f64`).

I mentioned these because they seem to be rarely used (though I'm no expert on that domain) but likely have some interesting properties, especially when functions behave a bit "weird". If I get again into programming something where I have an optimization problem, I'll take a look at `argmin` and consider adding these. Not sure if/when they are the best choice though.

This topic was automatically closed 90 days after the last reply. We invite you to open a new topic if you have further questions or comments.