I was just discussing with colleagues some of the classical theoretical oddities of interval arithmetic, in particular the fact that, for an interval X = [a,b], we have X-X \neq 0 (we have instead only the inclusion 0 \in X-X, since X-X = [a-b, b-a]). References: Wikipedia’s Interval arithmetic and the Interval arithmetic tutorial from IntervalArithmetic.jl documentation).
This brought me to recheck what were the available option for interval computations in Julia. I was surprised to find two packages which seems to come from two points of views which are far enough so that they are, to the best of my knowledge, not cross-referencing each others:
So I try to assemble here my notes for other people that may come across this potential choice. I welcome comments about what I may have missed or what I got wrong.
IntervalArithmetic.jl
IntervalArithmetic.jl, which seems to come from a Validated Numerics and Global optimization background is perhaps the closest package to what I understand to be the classical definition of Interval arithmethic. Coherent with this, we have the non-cancellation when substracting intervals:
julia> using IntervalArithmetic
julia> x = 5±1
[4, 6]
julia> x-x
[-2, 2]
Measurements.jl
On the other hands Measurements.jl seems to come from a Physics background (e.g. it belongs to the JuliaPhysics GitHub organization). The focus is not on Validated Numerics but rather on linear error propagation theory (which is by principle an approximation from the true error). Although the syntax can look exactly the same as with IntervalArithmetic (see code below), the manipulated objects are variables with some error (or combinations of such variables) rather than interval sets. This semantic changes allows Measurement to track what its author calls functionnally correlated quantities (see Giordano’s 2016 manuscript). With this in mind, we get:
julia> using Measurements
julia> x = 2±1
2.0 ± 1.0
julia> x-x
0.0 ± 0.0
Of course with both packages, we have 0*x somehow equal to zero (either represented 0.0 ± 0.0 or [0, 0] which are subtly different).
5 Likes
One other major difference I missed in my first post is the interpretation of the sum of intervals/variables:
IntervalArithmetic.jl
The sum of intervals is understood in the worst case sense (i.e. a Minkowski sum of the sets):
julia> using IntervalArithmetic
julia> x = 0±1
[-1, 1]
julia> y = 0±1
[-1, 1]
julia> x+y
[-2, 2]
Measurements.jl
The error term of each variable is understood is a standard deviation of a random variable (while the nominal value is the expectation). Then, the rule for defining the error for the sum of two variables with some error is the rule for summing random variables**.
If the random variables are independant, the square of the standard deviations add up (i.e. sum of variances):
julia> using Measurements
julia> x = 0±1
0.0 ± 1.0
julia> y = 0±1
0.0 ± 1.0
julia> x+y
0.0 ± 1.4
However, if the two variables getting summed are 100% correlated, the error gets doubled:
julia> x+x
0.0 ± 2.0
7 Likes
You might also be interested in GitHub - baggepinnen/MonteCarloMeasurements.jl: Propagation of distributions by Monte-Carlo sampling: Real number types with uncertainty represented by samples. which keeps track of a collection of samples rather than a mean and standard deviation.
Also GitHub - kalmarek/Arblib.jl: Thin, efficient wrapper around Arb library (http://arblib.org/) which is a wrapper around the C library Arb and is more similar to interval arithmetic (but using a midpoint-radius representation).
1 Like
It’s worth noting that midpoint-radius is pretty much strictly superior to interval representation since it allows storing much tighter intervals.
2 Likes
why do you say that? given you can losslessly convert between the two, I don’t really see how that’s possible, unless you mean the way they are usually implemented or something like that
You can’t lossless convert in finite precision arithmatic. The interval 1.0 +- eps(1.0)/4 isn’t representable using an interval.
3 Likes
Oh! I see, good point. So since there’s more floats near 0, you can represent very small radii, and the density of floats might be much lower out where the endpoints of the interval are?
3 Likes
Here’s an introduction to Affine arithmetic (AA), which captures perfectly the linear relationships between variables (thus provides the optimal range).
An over-approximation can be computed for nonlinear primitive functions. Overall, affine arithmetic is a powerful enclosure tool (although more expensive than naive interval arithmetic).
IIRC, Jordan Ninin combined (= intersected) both the IA and the AA enclosures in his thesis.
5 Likes
Since this is somewhat a misunderstanding sometimes, I want to stress there is no correct way to handle these things: as you mentioned it already, IntervalArithmetic.jl and Measurements.jl (+ MonteCarloMeasurements.jl) are generally used in different fields, for different purposes. And even within Physics, there are cases where linear error propagation is sufficient (in which cases MonteCarlo propagation would be an overkill, but still consistent) and others where it’d return pretty garbage results (for example whenever derivatives are very small relative to some scale, think of stationary points of functions, in these cases something more advanced as second order derivatives or MonteCarlo propagation is necessary).
I’d also like to mention that most of the complexity in Measurements.jl comes from tracking the “functional correlation” between quantities (I couldn’t find a better term to distinguish it from statistical correlation): it’s pretty easy to implement a linear error propagation system without tracking the correlations (this is what the very first release of Measurements.jl did), but then you start getting inconsistent results like x + x != 2 * x, x * x != x ^ 2, and so on, which are also related to the already mentioned x - x != 0. Not having basic identities to hold, or changing the error estimation depending on whether you apply addition or multiplication, when you mean the same thing, would be pretty annoying.
4 Likes
To (sometimes) address these problems, here’s a powerful but tedious way: if an expression is monotonic wrt a variable, the result is formed by taking the endpoint evaluations only.
Example:
Let f(x) = x - x and X = [-5, 5]. The natural (= naive) evaluation is F_N(X) = [-5, 5] - [-5, 5] = [-10, 10]. That’s bad.
However, f is monotonic (increasing or decreasing) wrt x because f' = 0. So the result is (for example) [f(\underline{X}), f(\underline{X})], which is [0, 0]. Obviously, this also works for more complex expressions.
2 Likes
Thanks for pointing me to this 2003 manuscript about Affine arithmetic by Stolfi which I didn’t know. Very interesting indeed. Yesterday I was wondering what would @giordano’s package Measurements.jl look like with a new option to treat error sums in a worst case sense rather than in its current root mean square sense.
After reading Stolfi’s paper, I think that this is exactly what AA is doing (with an extra bonus for managing both the roundoff and the truncation of Taylor series at the first order). I guess that some of the references at the end, like the “first-order Taylor arithmetic” are also implementing this idea.
I wonder if there are AA implementations in Julia.
Perhaps a new option within Measurements.jl may be able to do it (without the extra bonus of managing truncation errors which would be pretty outside the original philosophy of this package, I think).
As an alternative, since the definition of Affine Arithmetic (eq. 3.1 in Stolfi 2003) looks very similar to dual numbers (extended to manage several parameters), it may be possible to perform the AA computation using ForwardDiff to compute the input-output first-order sensitivities (again the truncation errors would not be handled). Then, the section 3.2 “Joint range of affine forms” seems to be readily applicable.
1 Like
Thanks for pointing me to the “ball arithmetic” keyword. I don’t have floating-point representation issues, but otherwise @Oscar_Smith’s explanation indeed makes the point about the advantage the midpoint-radius representation.
For the reference, I just quickly read Arb lib features overview where an additional point is made: the radius can be stored using a lower precision floating point.
3 Likes