[ANN] PrecisionCarriers.jl: Easy Detection of Floating Point Precision Loss

Package Announcements
, ,
1

I recently registered a new package I’ve been working on: PrecisionCarriers.jl. It’s designed to let you check existing code very quickly and easily for loss of precision in its calculations.

I added a macro inspired by BenchmarkTools.jl 's @benchmark macro, called @bench_epsilons. This lets you write a function call, interpolate values from your current scope into it, and define search ranges for variables that are then searched. The output gives a histogram of the observed precision loss, a printout of the number of “infinite” precision losses, if any, and a list of the worst observed precision losses, together with the function calls to reproduce them. Here’s an example of what it currently looks like in a particle physics application:

image
image1277×433 42 KB

Internally, this works by wrapping the inputs in a PrecisionCarrier object, which does every computation with both the given value and type, and a BigFloat version of it. At the end of the calculation, the values can be compared, and an estimated amount of machine epsilons of error can be estimated (or the number of remaining significant digits in the result).

Precision is a difficult thing to get right or have a feeling for, and I think it’s helpful to have a good guess about how precise the results of some calculations actually are, especially in scientific code. Contrary to popular belief, there are many more possible ways to introduce precision loss than just catastrophic cancellation or very long chains of arithmetic.

There are a lot of things that can still be improved here, and I have a few ideas in mind already, but I’d also be interested in what others think. Maybe some people want to give it a try and have feedback.

A few packages that I found are similar, but don’t seem intended to actually detect precision lossy code:

  • FloatTracker.jl specifically finds NaN and Inf values
  • DoubleFloats.jl, MultiFloats.jl, ArbFloats.jl, etc. provide floating point types with more precision
  • ValidatedNumerics.jl provides a way to have known error bounds on calculations, which is probably the closest to what I wanted to achieve with my package
11 Likes
2

I think you might want to use MultiFloats.jl rather than BigFloats for this. The result will be a lot faster.

3 Likes
3

That might be a good idea, I can look into that. I was also considering using DoubleFloats.jl, because those allow specific backend floating point types to be used, which should allow usage on GPUs, when they don’t have fp64 support (like one api GPUs).

4

I’ve looked at MultiFloats.jl a bit more, and I think I will not be using it as a backend for now. I don’t think speed is too much of an issue right now, and especially for a package concerned with finding precision errors, I dislike that MultiFloats don’t follow IEEE754 for the sake of speed. They say in their Readme that precision is lost in some cases, even for addition or subtraction (which, unless I’m misunderstanding, to me means that it’s not guaranteed that they round correctly), and they do not correctly propagate Inf and NaN values. This would lead to my package reporting infinite deviation between Inf and NaN, for example, when there might have otherwise been no deviation (Inf vs Inf).

However, I did find something quite fascinating (and initially very confusing) with my go-to unstable function:

function unstable(x, N)
    y = abs(x)
    for i in 1:N y = sqrt(y) end
    w = y
    for i in 1:N w = w^2 end
    return w
end

In theoretical perfect arithmetic, this always yields abs(x) regardless of N, but for any approximation, it will tend towards 1 (usually).
For example, with standard Float64 or Float32:

julia> for i in 1:5000 if unstable(1.5, i) == 1.0 print(i); break; end end
51
julia> for i in 1:5000 if unstable(1.5f0, i) == 1.0 print(i); break; end end
22

For the standard BigFloat:

julia> for i in 1:5000 if unstable(big(1.5), i) == 1.0 print(i); break; end end
254

So for MultiFloats, I was expecting something like ~100 for a Float64x2. Surprisingly:

julia> for i in 1:5000 if unstable(Float64x2(1.5), i) == 1.0 print(i); break; end end
1074

That’s a lot more than I expected. And also, it doesn’t change for the other FloatxN types:

julia> for i in 1:5000 if unstable(Float64x8(1.5), i) == 1.0 print(i); break; end end
1074

I’m assuming this happens because the hi part just approaches 1, while the lo part(s) become smaller and smaller but are not constrained to be the “extended mantissa”, like they are in BigFloat. So this approach only fails when the exponent of the lo part becomes too small.

Kinda cool and maybe interesting for use cases where you need to represent numbers very close to 1.

1 Like