It's NOT OK to compare floating-points using epsilons - Comments

It's NOT OK to compare floating-points using epsilons

natano

Thanks for posting this here. It's a good read!

doug-moen

Although I am not a fan of comparing floats using epsilons myself, I think it is good to read a defence of the idea written by people who have thought about this idea deeply, and who have done far more analysis than the OP.

The APL language has a concept of "tolerated comparison", which is like using an epsilon, but which uses a different formula than the OP, because it is quite tricky to implement this idea "correctly". In short, the epsilon used is not the same value for every comparison: it depends on the values being compared.

The problem is that these epsilons are pretty much always simply guessed, and there is no correct way to choose one epsilon out of many. If you have several such comparisons, no combination of epsilons will ensure that everything works correctly.

APL addresses this, by automatically computing an epsilon based on the magnitude of the values being compared. But:

Another problem with epsilons is that such comparison isn't transitive.

APL tolerated comparison isn't transitive, and I have a beef with this.

dpercy

APL tolerated comparison isn't transitive, and I have a beef with this.

I have thought about this a little bit! If it’s transitive then it can’t be tolerant everywhere.

If it’s transitive, then there must be equivalences classes, like islands, which must have boundaries. Points on a boundary can’t tolerate any error, at least in the direction that nudges it into the next island.

Conversely if every point has some neighborhood of values near it that it compares equal to, then comparison can’t be transitive—because the overlapping neighborhoods would chain together and make everything equal.

ptrettner

it depends on the values being compared.

This has also bitten me a lot in the past. Because the correct tolerance is almost never a local property. The most simple case is "a ~ b" vs "a - b ~ 0", but this can nest arbitrarily deep. When working on meshes, we often have a global notion of "this mesh was acquired in around this precision and changes below X are probably irrelevant to our users". But if you work locally on a face, suddenly vertex positions can be quite large compared to edge vectors, even though the formulas are often mathematically equivalent. So suddenly it matters a lot if you compare vertex positions or edge vectors or other derived properties. Any automatic local solution will fail and is a nightmare to debug.

(I'm in the "exact comparisons by default and be explicit about every use of epsilon" camp)

mzl

I would agree that the APL version of fuzzy equality is a useful version of it, and the properties of being scale-invariant are good (and that it isn't transitive is as you say a big reason to new treat it as an actual comparison).

However, in the few cases where fuzzy comparisons might be useful, it is not (in my experience) uncommon for there to be real-world reasons for the specific fuzzyness that is needed in each case, so these automatic choices would be wrong then.

doug-moen

The 2nd & 3rd links are written by Marshall Lochbaum, the author of BQN, who explains why BQN doesn't have tolerated comparison:

APL and J use approximate comparison in primitives, controlled by a value called the comparison tolerance (⎕CT in APL). The choice of which primitives use it and how is kind of arbitrary, and a nonzero comparison tolerance can lead to confusion, bugs, and unavoidable performance problems in some primitives. Nonetheless, I planned to add tolerant comparison to BQN—until I realized that after a year spent programming in BQN I'd hardly noticed its absence, and no one had asked for it either.

cespare

I'm also a certified approxEqual hater.

A few years ago I ran into a case where I was comparing the results of operations in two different languages that were intended to be equivalent. The math seemed the same, but the results were a tiny bit different.

Ignoring temptation to write an approxEqual function, I hunted down the source of the difference. Turns out Java and Go sometimes differed in the last bit of the exponential function, which is normal for transcendentals. Then subsequent operations compounded the error slightly, but not beyond 4 ULP. So instead of an approxEqual function, I wrote down a ulpDiff function which takes a number of ULP difference to tolerate. That felt like a more defensible approach to me.

Another opinion I've formed is that if you are going to write an approxEqual, it's usually much better to test the relative error, not the absolute error, as demonstrated in the post:

bool approxEqual(float x, float y)
{
    return abs(x - y) < 1e-4f;
}

If you are dealing with tiny numbers, 1e-4 might be an enormous difference!

Corbin

A well-done article. The error analysis is nice to see. I like the variety of examples.

You've probably seen [the square-root-of-dot-product algorithm for the Euclidean norm] a dozen times, so what's wrong with it? Nothing really, except that for really small vectors (those having length around 1e-19) it returns zero, and for larger vectors we get a significant precision loss. Even if the result is a perfectly valid floating-point value, its square (i.e. the expression under the square root) can be too small to be adequately representable in floating-point. … However, an invariant saying that only the zero vector has zero length would be extremely convenient, and would simplify writing code in a lot of cases, especially when we want to treat floating-points carefully.

There's actually a clever fix for this, origin unknown. For small or large vectors, the norm is fairly well-approximated by the largest (absolute value of each) coordinate in the vector. We can do a quick range check to see whether the coordinates are likely to under- or overflow and switch to the alternative algorithm. Using symbols instead of words but not literally using an APL relative:

len(v) = if all(UNDER <= v < OVER) then √sum(v²) else max(abs(v))

A complete example is in Monte's reference implementation, where I left an ASCII diagram and some comments in Double.euclidean/1 which correspond to the 2D norm. I don't know how to compute the bounds, but Herbie gave me magic constants for the 2D double-precision case.

jcelerier

The article mentions one case of uncertainty in the beginning but does not really discusses it afterwards: CPU (and I'd say more generally) platform differences. One of the case that bit me most is assuming that given the same input and the same source code, a FP computation will have the same result on two computers, which is definitely not true if for instance one does everything in x86 FPU 80-bit math and the other in SSE 64-bit math. If you have long-running state, error can really compound, and then if you're trying to sync your machines over the network (for instance for game state) you're in for a surprise

mzl

Nice to see a good explanation of lots of problems with using epsilons with floating poins.

I used to work with floating points in 2D geometry, and one of the issues we had was with sorting points: a common case was that we wanted to sort points to be left most, bottom-most. Using epsilons in this intuitively looks like a good idea, since small variations in the x-coordinate should not trump large variations in the y-coordinate. However, adding an epsilon to the comparison function (by using x-equality with epsilon) when sorting does not work, since it no longer fulfills the contract of a comparison function for sorting algorithms. This is not a theoretical issue either, it was easy to see happen in practice (many modern sorting algorithms will helpfully crash for you).

In addition to sorting algorithms not working with a 2D comparison function that uses epsilons, you could also get fully non-sensical results. It is easy to construct cases where you have a set of coordinates that intuitively looks ordered from right-most to left-most, but where each pair is very close in the x-direction, and clearly ordered correctly in the y-direction. This means that visually it looks like the order is reversed from the intended result, but the comparison function says it is already sorted.

wink

The only time I had to compare a lot of floats was in a C++ backend that got a lot of input from a web frontend, and thus JS, with its own fp problems, and even then you're talking language boundaries and somewhat arbitrary user input (for geo coordinates).

I still don't see a way to compare this without epsilon, but I'm not working in that codebase anymore.

mxey

How much human time have we collectively wasted because our computer's default representation of numbers with decimal points in them is inaccurate? :/

mqudsi

How much time have we saved by being able to approximate decimals with such a huge, dynamic range in a standardized fashion with strict rules for handling the most pernicious edge cases, while being able to fully accelerate the mathematics thereof with hardware execution units built into the CPU?

mqudsi

Definitely with reading, but not necessarily with sticking to; certainly evaluate your constraints and requirements before following the approaches offered!

The author makes an (accidental?) omission in the “computing vector length” section: they compare the “naive” calculation approach with the “patched-but-has-another-bug” solution to argue the recommended approach is “only” twice as many assembly instructions (already ignoring that they’re slow, floating point division operations being added!) but, crucially, then points out that to use the patched code you need to actually also introduce a branch into the logic!

natano

Thanks for posting this here. It's a good read!

doug-moen

Although I am not a fan of comparing floats using epsilons myself, I think it is good to read a defence of the idea written by people who have thought about this idea deeply, and who have done far more analysis than the OP.

The APL language has a concept of "tolerated comparison", which is like using an epsilon, but which uses a different formula than the OP, because it is quite tricky to implement this idea "correctly". In short, the epsilon used is not the same value for every comparison: it depends on the values being compared.

The problem is that these epsilons are pretty much always simply guessed, and there is no correct way to choose one epsilon out of many. If you have several such comparisons, no combination of epsilons will ensure that everything works correctly.

APL addresses this, by automatically computing an epsilon based on the magnitude of the values being compared. But:

Another problem with epsilons is that such comparison isn't transitive.

APL tolerated comparison isn't transitive, and I have a beef with this.

dpercy

APL tolerated comparison isn't transitive, and I have a beef with this.

I have thought about this a little bit! If it’s transitive then it can’t be tolerant everywhere.

If it’s transitive, then there must be equivalences classes, like islands, which must have boundaries. Points on a boundary can’t tolerate any error, at least in the direction that nudges it into the next island.

Conversely if every point has some neighborhood of values near it that it compares equal to, then comparison can’t be transitive—because the overlapping neighborhoods would chain together and make everything equal.

ptrettner

it depends on the values being compared.

This has also bitten me a lot in the past. Because the correct tolerance is almost never a local property. The most simple case is "a ~ b" vs "a - b ~ 0", but this can nest arbitrarily deep. When working on meshes, we often have a global notion of "this mesh was acquired in around this precision and changes below X are probably irrelevant to our users". But if you work locally on a face, suddenly vertex positions can be quite large compared to edge vectors, even though the formulas are often mathematically equivalent. So suddenly it matters a lot if you compare vertex positions or edge vectors or other derived properties. Any automatic local solution will fail and is a nightmare to debug.

(I'm in the "exact comparisons by default and be explicit about every use of epsilon" camp)

mzl

I would agree that the APL version of fuzzy equality is a useful version of it, and the properties of being scale-invariant are good (and that it isn't transitive is as you say a big reason to new treat it as an actual comparison).

However, in the few cases where fuzzy comparisons might be useful, it is not (in my experience) uncommon for there to be real-world reasons for the specific fuzzyness that is needed in each case, so these automatic choices would be wrong then.

doug-moen

The 2nd & 3rd links are written by Marshall Lochbaum, the author of BQN, who explains why BQN doesn't have tolerated comparison:

APL and J use approximate comparison in primitives, controlled by a value called the comparison tolerance (⎕CT in APL). The choice of which primitives use it and how is kind of arbitrary, and a nonzero comparison tolerance can lead to confusion, bugs, and unavoidable performance problems in some primitives. Nonetheless, I planned to add tolerant comparison to BQN—until I realized that after a year spent programming in BQN I'd hardly noticed its absence, and no one had asked for it either.

cespare

I'm also a certified approxEqual hater.

A few years ago I ran into a case where I was comparing the results of operations in two different languages that were intended to be equivalent. The math seemed the same, but the results were a tiny bit different.

Ignoring temptation to write an approxEqual function, I hunted down the source of the difference. Turns out Java and Go sometimes differed in the last bit of the exponential function, which is normal for transcendentals. Then subsequent operations compounded the error slightly, but not beyond 4 ULP. So instead of an approxEqual function, I wrote down a ulpDiff function which takes a number of ULP difference to tolerate. That felt like a more defensible approach to me.

Another opinion I've formed is that if you are going to write an approxEqual, it's usually much better to test the relative error, not the absolute error, as demonstrated in the post:

bool approxEqual(float x, float y)
{
    return abs(x - y) < 1e-4f;
}

If you are dealing with tiny numbers, 1e-4 might be an enormous difference!

Corbin

A well-done article. The error analysis is nice to see. I like the variety of examples.

You've probably seen [the square-root-of-dot-product algorithm for the Euclidean norm] a dozen times, so what's wrong with it? Nothing really, except that for really small vectors (those having length around 1e-19) it returns zero, and for larger vectors we get a significant precision loss. Even if the result is a perfectly valid floating-point value, its square (i.e. the expression under the square root) can be too small to be adequately representable in floating-point. … However, an invariant saying that only the zero vector has zero length would be extremely convenient, and would simplify writing code in a lot of cases, especially when we want to treat floating-points carefully.

There's actually a clever fix for this, origin unknown. For small or large vectors, the norm is fairly well-approximated by the largest (absolute value of each) coordinate in the vector. We can do a quick range check to see whether the coordinates are likely to under- or overflow and switch to the alternative algorithm. Using symbols instead of words but not literally using an APL relative:

len(v) = if all(UNDER <= v < OVER) then √sum(v²) else max(abs(v))

A complete example is in Monte's reference implementation, where I left an ASCII diagram and some comments in Double.euclidean/1 which correspond to the 2D norm. I don't know how to compute the bounds, but Herbie gave me magic constants for the 2D double-precision case.

jcelerier

The article mentions one case of uncertainty in the beginning but does not really discusses it afterwards: CPU (and I'd say more generally) platform differences. One of the case that bit me most is assuming that given the same input and the same source code, a FP computation will have the same result on two computers, which is definitely not true if for instance one does everything in x86 FPU 80-bit math and the other in SSE 64-bit math. If you have long-running state, error can really compound, and then if you're trying to sync your machines over the network (for instance for game state) you're in for a surprise

mzl

Nice to see a good explanation of lots of problems with using epsilons with floating poins.

I used to work with floating points in 2D geometry, and one of the issues we had was with sorting points: a common case was that we wanted to sort points to be left most, bottom-most. Using epsilons in this intuitively looks like a good idea, since small variations in the x-coordinate should not trump large variations in the y-coordinate. However, adding an epsilon to the comparison function (by using x-equality with epsilon) when sorting does not work, since it no longer fulfills the contract of a comparison function for sorting algorithms. This is not a theoretical issue either, it was easy to see happen in practice (many modern sorting algorithms will helpfully crash for you).

In addition to sorting algorithms not working with a 2D comparison function that uses epsilons, you could also get fully non-sensical results. It is easy to construct cases where you have a set of coordinates that intuitively looks ordered from right-most to left-most, but where each pair is very close in the x-direction, and clearly ordered correctly in the y-direction. This means that visually it looks like the order is reversed from the intended result, but the comparison function says it is already sorted.

wink

The only time I had to compare a lot of floats was in a C++ backend that got a lot of input from a web frontend, and thus JS, with its own fp problems, and even then you're talking language boundaries and somewhat arbitrary user input (for geo coordinates).

I still don't see a way to compare this without epsilon, but I'm not working in that codebase anymore.

mxey

How much human time have we collectively wasted because our computer's default representation of numbers with decimal points in them is inaccurate? :/

mqudsi

How much time have we saved by being able to approximate decimals with such a huge, dynamic range in a standardized fashion with strict rules for handling the most pernicious edge cases, while being able to fully accelerate the mathematics thereof with hardware execution units built into the CPU?

mqudsi

Definitely with reading, but not necessarily with sticking to; certainly evaluate your constraints and requirements before following the approaches offered!

The author makes an (accidental?) omission in the “computing vector length” section: they compare the “naive” calculation approach with the “patched-but-has-another-bug” solution to argue the recommended approach is “only” twice as many assembly instructions (already ignoring that they’re slow, floating point division operations being added!) but, crucially, then points out that to use the patched code you need to actually also introduce a branch into the logic!