RFC: Fix to how eps() works on complex values and types

[marcusps]

The current implementation of eps on complex values relies on abs. This leads to a problem for complex integer types, as abs returns a float, and thus eps(1im) does not fail like eps(1) (and it should). Similarly, eps(Complex{Float64}) and related expressions simply fail (not being defined).

The changes in this request fix both of these issues, which become important when trying to compare complex floating point values. The docs (doc/stdlib/base.rst specifically) mentions only eps of real Float values & types is sensible, but it strikes me as natural to support complex Floats.

[StefanKarpinski]

I'm pretty dubious about eps for complex numbers in general – it's kind of unclear what it should even mean. For floats is very clear – it's the difference between x and the next representable float. What's the use case for complex eps? I wonder if there isn't a saner thing to do, or at least having a use case might motivate what the "correct" behavior ought to be.

[marcusps]

I agree, it wasn't clear to me what was the proper definition in the complex case. Since it was already implemented, however, I just wanted to make it consistent with the definition for real Floats and Integers.

The only use case I have in mind is checking for approximate equality. Thinking of complex values are pairs of real values leads to a natural approach to comparison: if both the real and imag parts are approximately equal, then so are the complex values (this is somewhat analogous to the infinity norm for vectors).

[JeffBezanson]

It seems to me the distance between complex numbers is abs(z1-z2), which is probably why abs was used.

[ViralBShah]

Cc: @alanedelman

[marcusps]

I think abs was being used before just to get a real floating point number before. If it was relying on the distance between complex numbers, I would have expected it would have used abs(eps(real(z))+1.im*eps(imag(z))), but this also strikes me as somewhat odd as it would imply that the precision of complex numbers is worse than the precision of real numbers.

It all goes back to @StefanKarpinski's point, I suppose, and what kind of distance is more reasonable to use (if any, in the complex case).

[ViralBShah]

Also, there's an inconsistency in the current implementation, where eps works for complex integers and not for real integers. It should not work on complex integers.

julia> eps(1+im)
2.220446049250313e-16

julia> eps(1)
ERROR: no method eps(Int64,)

I am not even sure what eps means for complex numbers, and doing eps(abs(z)) just feels weird. I'd much rather not have it.

[jiahao]

As mentioned already, it really depends on how you want to define distance between complex numbers.

1. IEEE 754 does not define complex floating point, nor does it define arithmetic on it, nor does it define machine epsilon.
2. If we take as an operational definition that used by LAPACK, which appears to be the most common one, then the use of abs() for real numbers extends naturally to the Euclidean L²-norm of complex numbers, and is thus always real-valued. The latter is also conventionally written |...|. If we take this definition, then eps(abs(z)) would not compute epsilon for complex numbers correctly, and eps(Complex{Float64}) should not fail, but instead return a real-valued float that is ~ sqrt(2) eps(Float64).

The Euclidean norm is the most useful for complex operations, but there are other norms (Wikipedia) one could take such as the L^p norms. These range in absolute value between the L^∞-norm where |z1-z2|=max(abs(real(z1-z2)), abs(imag(z1-z2))) and L¹-norm where |z1-z2|=abs(real(z1-z2)) + abs(imag(z1-z2)). However to my knowledge these are not often used for complex numbers.

So perhaps the correct thing to do is to have eps() extendable to various norms that define different distance metrics between complex numbers, but default it to something sensible, say the Euclidean norm, in which case the existing abs() implementation would be correct.

[jiahao]

@marcusps "the precision of complex numbers is worse than the precision of real numbers."

It is not necessarily that surprising considering the number of operations on real-valued quantities one has to perform to do complex arithmetic, since there are twice as many real-valued numbers to work with in the latter.

[marcusps]

This seems to be getting into murky territory. Here are what appear to be the clear issues, suggested fixed

1. eps handles integer types inconsistently, so disallow eps on any integer type, whether it is real or complex
2. Define eps only for real values. Approximate equality of complex values and arrays could be handled by extending extras/nearequal.jl (or its successor), which internally should rely on eps of Floating point representations

I think this meshes well with the reservations that @StefanKarpinski haves, addresses some of the issues @ViralBShah and I have pointed out, and avoids the murkier questions brought up by myself, @jiahao and @JeffBezanson about how eps relates to norms (or whether it should relate to norms at all) by defining it as a property of the finite precision representation.

Thoughts?

[StefanKarpinski]

Getting into norm stuff feels like the wrong direction to me – eps isn't about math, it's about the inherently granular nature of floating-point arithmetic. In general, I find nextfloat and prevfloat to be much more useful, but people seem to like eps. We need to figure what the use cases are and see if they make sense when we extend them to complex values.

I feel like one of the classic use cases is to make an epsilon-ball around one value to compare it with another, to check for containment in the epsilon ball. I guess that's where the norm stuff comes in, since you typically use a norm like abs to check for containment.

[jiahao]

@marcusps (2) makes sense. I did a cursory search of the numerical analysis literature and didn't turn up anything particularly definitive on complex machine epsilon. It's not too surprising, since any sensible definition of eps for complex floating-point must necessarily be representation specific. A polar (angle, magnitude) representation would have very different machine-epsilon behavior to the usual (real, imaginary) representation.

@StefanKarpinski I think the discussion about norm is entirely relevant, since the whole point of eps is to measure distance between machine-representable numbers, and norms are about measuring distances, and it is really a question of how to go sanely beyond the real number line. I think there are three reasonable ways to extend eps:

1. Don't extend it beyond the real line,
2. Return a vector of finite-difference derivatives with the smallest measurable step in every possible direction,
3. Return some reasonable norm of the vector above.

Personally I only care about machine epsilon to the extent that I need to test for approximate equalities, and secondarily in error propagation to make sure that roundoff error isn't washing out all my calculations. Both of these cases are handled by epsilon-ball.

[jiahao]

A quick inspection of LAPACK also shows that eps is calculated only in the single- and double-precision real floating-point routines {S,D}LAMCH, which are just wrappers around the EPSILON intrinsic in Fortran, which is in turn only defined for REAL(*).

[ViralBShah]

Removing eps(Complex) is the right thing to do here.