Approximate equality

#1

I read somewhere that the approximation binary operator (\approx) checks if the number is within machine precision. Whatever it does exactly, it does not appear to behave consistently:

julia> ℯ^(-π*im)
-1.0 - 1.2246467991473532e-16im

julia> ℯ^(-π*im) == -1
false


This is probably normal behaviour, and I guess that \sim 10^{-16} is machine or type precision. So I will use \approx to get a trueth, but reordening of symbols, or considering only the imaginary part, results in a falsehood:

julia> ℯ^(-π*im) ≈ -1
true

julia> ℯ^(-π*im) + 1 ≈ 0
false

julia> imag(ℯ^(-π*im))
-1.2246467991473532e-16

julia> imag(ℯ^(-π*im)) ≈ 0
false


Tip: typing \app[Tab][Tab] in the interpreter results in TeX’s \approx and then Unicode.

Most of the above numerical expressions are of the same type Complex{Float64}, so precision should be the same for everything. The last two imag() expressions are Float64 and shows the problem as well. What happens here?

I am using Julia 0.7.0-DEV.3686, same on 0.6.2.

isposdef inconsistent results, numerically unstable?
#2

What is happening is shown in the docs for isapprox.

isapprox(x, y; rtol::Real=sqrt(eps), atol::Real=0, nans::Bool=false, norm::Function)

Inexact equality comparison: true if norm(x-y) <= atol + rtol*max(norm(x), norm(y))


Why the rules are set like that, I’m not sure.

#3

There was much discussion about this at one time. For isapprox to work reasonably over pairs that may be very close, close, far or very far from each other , it is necessary to combine both the absolute difference and the relative inexactness and that is usually done additively. Norms are taken to normalize the metric over various kinds of input.

#4

I have been breaking my head over Julia’s definition of isapprox(). I was missing an important ingredient, namely the fact that computers approximate real numbers by floats! For floats the precision depends on the magnitude of the number (and it depends on how many bits the float has). For instance,

julia> isapprox(1 + 1.2e-16, 1)
true

julia> isapprox(1.2e-16, 0)
false


To some extend, this is consistent how for instance physicists think about approximation. Something like a \approx b \Leftrightarrow a/b \approx 1. Approximation is relative. This implies that “approximately zero” doesn’t make sense (ratio is ill-defined), also not in Julia:

julia> isapprox(1.2e-300, 0)
false


#5

I didn’t know that Julia has such a function, and I propose that it should be dropped from the language. For users with experience and understanding of floating-point, this function is superfluous since they can use norm and abs appropriately. For users lacking experience, this function gives them a false sense of security as they enter a mine-field.

#6

I’m not sure if I agree. If I see an “\approx” somewhere of which I don’t know its definition yet, I won’t assume any type of security at all. If I don’t know what it means I either don’t use it, or I try to understand what it means first.

If Julia’s isapprox(a, b) is indeed consistent with what is usually understood by a/b\approx 1, I propose to leave it as-is. Otherwise, I don’t know (it may be there for a different reason).

#7

I understand the mathematical justification for isapprox(a,0.0) being ill-defined, but it is pretty confusing for users.

Is there a meaningful way to define
isapproxzero(a)? and if so could the a \approx 0 be written to support that particular form for the integer 0 literal?. i.e. something like

function isapprox(x, zero_literal???; rtol::Real=sqrt(eps), atol::Real=0, nans::Bool=false, norm::Function)
return norm(x) < rtol
end


At this point, since x \approx 0 is pretty much useless, I don’t see what the harm is?

#8

This was discussed in https://github.com/JuliaLang/julia/issues/23376 and the consensus was not to do it.

#9

I think I understand. So the issue is that the units of the default tolerance only make sense for the particular norm you are interested in?

#10

The issue is that when you say that a is approximately equal to b, it means that |b-a| is small compared to the characteristic magnitude of a and b, which is taken to be max(a,b) for lack of a better way. If b is zero, then that procedure fails. The alternative would be to define the smallness absolutely (which you can do by yourself with isapprox(rtol=Inf, atol=your_tol)), but this really assumes that the quantities you’re interested in have a natural unit of about one, which might or might not be true. 1e-10 meters is small when you’re comparing the length of an atom to the length of your arm, it’s not so small when you’re comparing the length of an atom to the Planck length or something like that. Really in any well-behaved computation you should normalize everything to a relevant set of units, and everything will be about unity, so this would make sense, but it’s probably not a decision that the base language should make (see stevengj’s comment in the linked issue)

#11

Makes complete sense. Is there ever a time when a \approx 0 or a \approx 0.0 is meaningful? If not, is there a way to give a warning to users who may not have thought this through, so it silently fail? I have personally made the mistake before

#12

My understanding is that isapprox(a, 0.0) == false for an a != 0.0. I prefer |y - x| \le \epsilon\cdot(1+|x|+|y|) in practice. But no single predicate will satisfy all expectations, and it is useful to have something in Base.

#13

If it is not something then at least it should warn the user that it is not doing what they think it is doing. i.e. sounds like if min(abs(a),abs(b)) == 0 but a != b EXACTLY (which often fails for floating point operations) then the function should do a warning or error

#14

Approximate comparison to an explicit zero is equivalent to exact comparison but you can have x ≈ y where both end up being exactly zero and it will succeed.

Idle thought: would it make any sense to consider subnormal floats to be approximately zero? Or some subset of subnormal floats?

#15

We could tweak it to be something like

abs(x-y) <= atol + rtol*max(abs(x), abs(y), realmin(typeof(x)), realmin(typeof(y)))


That way by default any Float64 of magnitude less than realmin()*sqrt(eps()) == 3.3156184e-316 would be approximately equal to zero.

#16

That seems like it would fit the “half the significant digits” criterion.

#17

half the sigdigs when you are focusing on values which by definition are shy some bits of significance (the subnormals) seems overly fine a comb.

any float that is < realmin(typeof(float)) is already zeroish given it preceeds the realmin and is >= 0

#18

I do agree that the current def is less elegant than it might be.

#19

One option would be to make the default value of atol to be the max of the realmins?

#20

that sounds decent – if the max(realmin) is conditioned on the type of the values being approxequaled rather than
maximum(realmin(Float64), realmin(Float32), realmin(Float16))