0

We can not compare doubles directly using == so one of the recommended methods is to use a threshold epsilon in order to determine equality based on our precision standard as defined by the threshold.
What I have been noticing is that when having 2 doubles the check usually recommended in relevant posts is if Math.abs(a - b) < epsilon which is the part that is confusing to me.
In my understanding the comparison of the diff should be taking into account the magnitude of the numbers and not using the direct diff.
Example:
Assuming the threshold is 0.000001 in the following cases where we try to establish equality with a

double a = 57.33;
double b = 57.32973229;
double c = 57.33000002;

b would be rejected since 57.33 - 57.32973229 = 0.00026771 > 0.000001 but to me it sounds quite unreasonable given the fact that this is less than 0.0004 error (0.00026771/57.33). This is even more obvious with larger and larger (or smaller and smaller) numbers.
c would be accepted since 57.33000002 - 57.33 = 0.00000002 < 0.000001

It seems quite impractical unless in very specific situations to accept as equal only c.
What am I missing/misunderstanding here?

Update:
So why isn't the recommended approach (a - b)/Max(a,b) < epsilon instead?

Jim
  • 3,845
  • 3
  • 22
  • 47
  • Could you clearly state the question please? – Dropout Aug 04 '20 at 15:39
  • 1
    "0.00026771 < 0.000001" should be "0.00026771 > 0.000001". I tried to edit it, but edits must be at least 6 characters. – Alexander Guyer Aug 04 '20 at 15:40
  • I think you reason right. For large numbers small epsilon will result in inequality always. Can you tell more about your calculations? – hal Aug 04 '20 at 15:41
  • If you are suggesting that two floating point numbers should be considered equal if they are within some relative difference (rather than absolute difference) of each other, then you can of course do this. For instance, you could accept a 0.001% margin of error, if desired. However, floating point precision issues generally occur on an absolute scale (rather than relative) due to only having so many bits and power-of-two denominators. From my understanding, if you multiple two semi large numbers, the resulting information loss will not be proportionally higher. Someone correct me if I'm wrong. – Alexander Guyer Aug 04 '20 at 15:44
  • 1
    @Dropout: I updated the post. Does that help? – Jim Aug 04 '20 at 16:25
  • @hal: I updated the post. – Jim Aug 04 '20 at 16:26
  • @Nerdizzle: Corrected the post. Thanks for that – Jim Aug 04 '20 at 16:28
  • It depends on your problem. You can not tell that one is better from the other. – hal Aug 04 '20 at 16:30
  • @Jim absolutely – Dropout Aug 04 '20 at 16:32
  • @hal: Comparing the diff against the epsilon is giving as close as possible equality. So in that sense one could say it is better regardless of context. But on the other hand i.e. comparing the ratio instead seems more practical. I am confused on all the parameters one should take to reach a decision – Jim Aug 04 '20 at 16:42
  • If you do not comare large numbers then choose: [Math.abs(a - b) < epsilon]. – hal Aug 04 '20 at 16:50
  • “We can not compare doubles directly using `==`” is incorrect. The `==` operation works perfectly with no errors. It evaluates to true if and only if its two operands are equal. The actual problem is producing desired operands to test and/or knowing the characteristics of the values you have produced, not in testing them. [There is no general solution for comparing floating-point numbers that contain errors from previous operations.](https://stackoverflow.com/a/21261885/298225). All “solutions” are situation-dependent, so a specific answer cannot be provided without specific information. – Eric Postpischil Aug 04 '20 at 17:56
  • @EricPostpischil:In the post you linked the term "tolerance" means the epsilon I use in my post? – Jim Aug 04 '20 at 21:01
  • @Jim: Yes. “Epsilon” should be avoided because it has multiple meanings. It is used in general for any small number, in floating-point to refer to the step size from 1 to the next greater representable value (also called the ULP, unit of least precision, of 1), and in comparisons like this, as well as its meaning with delta-epsilon proofs in limits (calculus). – Eric Postpischil Aug 04 '20 at 21:08
  • `Math.abs(a - b) < epsilon` should be fine unless `epsilon` is too small. That is, too close to 2^-53. If it is too close, it may fail to allow two "close" values to count as "=". – Rick James Aug 04 '20 at 22:14
  • @EricPostpischil:Is the red flag you raise on this the term "epsilon" or the fact that a threshold number (termed epsilon) is usually recommended for check of equality? – Jim Aug 04 '20 at 22:34
  • @RickJames: `2^-53` is `0.00000000000000011102230246251565404236316680908203125`. If that is the threshold then it accentuates what I mention in the post since only something like `57.330000000000000001` would be considered as equal – Jim Aug 04 '20 at 22:38
  • @RickJames: Unless I misunderstood you point completely – Jim Aug 05 '20 at 08:31

1 Answers1

3

In my understanding the comparison of the diff should be taking into account the magnitude of the numbers and not using the direct diff.

In one floating-point operation, the computed result is the real-number result rounded to the nearest representable value. The distances between nearby values are scaled according to the magnitude of the result (with some quantization: the steps stay a constant size throughout one binade [all the numbers represented with the same exponent] and then jump as you transition to a neighboring binade). Therefore, for one floating-point operation, there is a known bound on the error that is proportional to the magnitude of the result.

So why isn't the recommended approach (a - b)/Max(a,b) < epsilon instead?

When a computation involves multiple operations, there may be a rounding error in each operation. That rounding error has a known bound that is proportional to the magnitude of its results. But as a computation proceeds through multiple steps, errors from previous steps may be compounded (or canceled) by various operations. So the final error is not proportional to Max(a, b) but is affected by all the magnitudes of the intermediate results as well as interactions between the errors and operations.

As an example, consider some very large number x that resulted from previous operations. It may have some error e that is large, but still proportional to x. If we subtract from x a number y that is near x in magnitude, we get a small result, a, includes that error e. And that error e is proportional to x, so it may be huge compared to a, very disproportionate compared to earlier errors.

This is part of the reason there is no general solution for accepting as equal floating-point numbers that are not equal: Just from the final results a and b alone, we cannot know how much error could have accumulated in them. Max(a, b) is insufficient. The error could be anything between zero and infinity, inclusive, and it could also be NaN. Designing a solution for a particular situation requires knowledge of the computations used to obtain the results and the input data.

Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
  • The part `with some quantization: the steps stay a constant size throughout one binade ...` would it be possible to explain it a bit more? – Jim Aug 04 '20 at 20:57
  • @Jim: In the `double` format, numbers are represented as ±x•2^e, where x is a 53-bit number less than 2 (such as 1.1010000…000, which would be 1⅝). So the steps between representable numbers near such a number are 0.0000…001•2^e. Since that x part, called the significand, has 53 bits, the last bit is at position 2^52, so 0.0000…001•2^e = 2^−52•2^e = 2^(e−52). All numbers with exponent e are spaced 2^(e-52) apart. That step size stays the same until we get to a point where the exponent changes. – Eric Postpischil Aug 04 '20 at 21:05
  • @Jim: For example, the steps between all numbers from 16 to 32 are 2^(4-52)=2^-48. The steps between all numbers from 32 to 64 are 2^(5-52) = 2^-47. When a calculation is done, the rounding error is at most ½ the step size, because we round to the nearest result, so it is at most ½ step away in one direction or the other. – Eric Postpischil Aug 04 '20 at 21:06
  • `the steps stay a constant size throughout one binade... the steps between all numbers` what is the definition of steps in this context? Is there some already specified algorithm that define what a step is? – Jim Aug 05 '20 at 08:32
  • If I understand your answer to its entirety (very insightful by the way, thank you very much) it is an error to define one generic function taking 2 double parameters and a threshold and return true if they are within the threshold? I.e. something like this https://commons.apache.org/proper/commons-math/javadocs/api-2.2/org/apache/commons/math/util/MathUtils.html#equals(double,%20double,%20double) since it doesnt include domain specific knowledge of the accumulated error? – Jim Aug 05 '20 at 08:35
  • @Jim: By “step” I simply refer to the distance between representable numbers. Suppose we had a two-digit decimal floating-point format, so it could represent numbers like 1.0, 1.1, 1.2, or 3.7•10^−1, 3.8•10^−1, 3.9•10^−1, or 7.4•10^5, 7.5•10^5, 7.6•10^5. Then the step between the adjacent representable numbers 1.1 and 1.2 is .1, the step between adjacent representable numbers 3.7•10^−1 and 3.8•10^−1 is .1•10^−1, and hte step between adjacent representable numbers 7.5•10^5 and 7.6•10^5 is .1•10^5. – Eric Postpischil Aug 05 '20 at 10:08
  • 1
    @Jim: I discourage pushing routines like that and argue that it is bad design and promotes incorrect use of floating-point arithmetic. – Eric Postpischil Aug 05 '20 at 10:10