relative_eq() has incorrect behavior for small values

[dsroche]

Consider this small example:

let diff: f32 = 1.0e-7;
let x: f32 = 1.0e-5;
let y = x - diff;
assert_ne!(x, y);
relative_eq!(x, y) // returns true, but should be false

Note that x and y are normal, full-precision floating-point values (not "subnormals"), above the epsilon threshold, and their relative error is (x - y)/x, which equals 0.01, far above the default threshold. However, they are marked as "relatively equal" in the current implementation.

The problem is in the logic of the provided implementation of RelativeEq::relative_eq() for f32 and f64, namely these lines:

approx/src/relative_eq.rs

Lines 65 to 70 in 308176d

	let abs_diff = $T::abs(self - other);
	// For when the numbers are really close together
	if abs_diff <= epsilon {
	return true;
	}

While the precise definition of relative equality and the parameters epsilon and max_relative are never clearly specified in the documentation, based on the referenced blog posts, it seems that the intention is to use epsilon as a zero-closeness threshold, and max_relative as the relative error threshold.

The issue stems from not checking the two input numbers (self and other in this case) individually against epsilon for a zero-threshold, but instead testing their absolute difference against epsilon. As a result, relative equality works as expected for large numbers, but in practice defaults to absolute-difference equality for values less than 1.

[Stargateur]

Rust follow IEEE 754 all is specified here. (but I have no idea if this will answer you)

[imjasonmiller]

@dsroche does the code below illustrate your point? I'm not that familiar with relative floating point comparisons:

fn relative_eq(lhs: f32, rhs: f32, epsilon: f32, max_relative: f32) -> bool {
    // Handle same infinities
    if lhs == rhs {
        return true;
    }

    // Handle remaining infinities
    if f32::is_infinite(lhs) || f32::is_infinite(rhs) {
        return false;
    }

    let abs_diff = f32::abs(lhs - rhs);

    // For when the numbers are really close together
    if abs_diff <= epsilon {
        return true;
    }

    let abs_lhs = f32::abs(lhs);
    let abs_rhs = f32::abs(rhs);

    let largest = if abs_rhs > abs_lhs { abs_rhs } else { abs_lhs };

    println!("This part is never reached");
    // Use a relative difference comparison
    abs_diff <= largest * max_relative
}

fn main() {
    let diff = 1.0e-7;
    let x: f32 = 1.0e-5;
    let y: f32 = x - diff;

    println!(
        "relative_eq!(x, y) = {}",
        relative_eq(x, y, f32::EPSILON, f32::EPSILON) // true
    );

    println!(
        "abs_diff <= largest * max_relative = {}",
        (x - y).abs() <= 1.0e-5 * f32::EPSILON, // false
    );
}

So one has to pass their own epsilon in order for it to work for small values? E.g:

use approx::relative_eq;

relative_eq!(x, y, epsilon = 1.0e-8f32)

[brendanzab]

Oooh thanks for the report! I may have got something wrong in the implementation - it's been a long time since I looked inside. And to be honest floating point numerics are not something I'm an expert on. I'd be interested if we can find a fix for this though!

[kitt159]

From my experiments, the code of relative_eq looks correct (it is also the same as here), but the default value for epsilon is incorrect. It should probably be f64::MIN_POSITIVE instead of f64::EPSILON. f64::MIN_POSITIVE is exactly the border between normal and subnormal floats - this is the place where the relative test can fail because of insufficient precision.

[jonaspleyer]

I have not come around to fix this bug. Is there any progress/new insights to be gained here?

[kitt159]

I will try to fix it as I described above. I will see if it works, and eventually create PR.

[kitt159]

I investigated this more deeply and I believe that my fix replacing EPSILON with MIN_POSITIVE is correct.

(For simplicity, I will describe the behavior only for positive numbers represented as f64; however, the same applies in general to negative numbers and to f32.)

The current implementation performs relative comparisons only in the range from EPSILON to MAX (approximately from 2.22e-16 to 1.79e308). When using MIN_POSITIVE, relative comparison is applied in the range from MIN_POSITIVE to MAX (approximately from 1.79e-308 to 1.79e308). In my opinion, absolute comparison is still necessary for numbers smaller than MIN_POSITIVE, because precision is lost in this range and relative comparison would not be reliable.

While experimenting, I also noticed the same issue in the ulps_eq implementation.
For large numbers, the code works correctly:

let num = 1e100f64;
let mut shifted = num;
for _ in 0..4 {
    shifted = shifted.next_up();
}
assert_ulps_eq!(num, shifted, max_ulps = 4);    // shifted number is exactly 4 ulps from num -> equal
assert_ulps_ne!(num, shifted.next_up(), max_ulps = 4);  // shifted number is 5 ulps from num -> not equal

However, the same code fails for small numbers:

let num = 1e-100f64;
let mut shifted = num;
for _ in 0..4 {
    shifted = shifted.next_up();
}
assert_ulps_eq!(num, shifted, max_ulps = 4);    // equal
assert_ulps_ne!(num, shifted.next_up(), max_ulps = 4); // fails - also equal?!

This bug is again caused by absolute comparison using the EPSILON constant. In my opinion, code performing ULP-based comparisons should not contain absolute comparison at all. It should rely purely on the bit representation of floating-point numbers and their distance in the float space.

However, there is a complication that makes fixing this harder. This would be a significant change in behavior (especially for ulps_eq), and changing it could break dependent code. In the case of relative_eq, this could be handled by informing users that to preserve the original behavior, they can explicitly specify epsilon = f64::EPSILON, which may or may not be correct depending on the use case.

For ulps_eq, the situation is more complex, because at the moment it is not even clear to me how the behavior should be fixed correctly. We could again use comparison against MIN_POSITIVE, but in my opinion a more correct solution would be to remove absolute comparison entirely. This could be achieved by changing the default value of the epsilon parameter to 0, effectively disabling absolute comparison. Another option would be to remove the epsilon parameter altogether, which would be cleaner but would introduce a backward-incompatible change to the public API.

[jonaspleyer]

As far as I can tell, the Ulps trait is the least used one of all so breaking something there would probably have minor consequences in practice. The upcoming release 0.6.0 introduces quite a lot of new features and thus it would make sense to also include this change. So far, I do not see another major release coming very soon.

I will try to read about this issue myself but I am very restrained with my time at the moment. So this may take a while. But no need to rush. I believe it is better to sort this out properly instead of just pursuing a quick hack.

[jonaspleyer]

@kitt159 Would you be willing to repurpose some of the explanations used here to include them in the documentation of the crate. It was criticized in the past that the documentation with respect to the functionality of the crate was lacking.

see #54

[kitt159]

I didn't have time to look at it this weekend, maybe the next one... I would wait with the documentation until the fix is implemented.