# The test isequal(+0.0,-0.0) returns FALSE

so as `isequal(+0.0im,-0.0im)`.

I discovered this during a painful debug process.

I really want to understand the philosophy behind this designing, and I wish to get some advice on numerical stability. Thank you @everyone !

If you want `+0.0` and `-0.0` to be equal, use `==`, which implements IEEE floating-point equality rules. `isequal` is used more for hash tables and other purposes that require finer-grained distinctions.

(In Julia, there are three different equality predicates: `==`, `===`, and `isequal`, to express slightly different notions of equality.)

7 Likes

As an aside, I learned recently from Valentin that `-0.0` is the additive identity, while `0.0` is not:

``````julia> mymulv1(a,b) = muladd(a,b,0.0)
mymulv1 (generic function with 1 method)

mymulv2 (generic function with 1 method)

julia> @code_native debuginfo=:none mymulv1(1.2, 3.4)
.text
vxorpd  %xmm2, %xmm2, %xmm2
vfmadd213sd     %xmm2, %xmm1, %xmm0 # xmm0 = (xmm1 * xmm0) + xmm2
retq
nopw    (%rax,%rax)

julia> @code_native debuginfo=:none mymulv2(1.2, 3.4)
.text
vmulsd  %xmm1, %xmm0, %xmm0
retq
nopw    %cs:(%rax,%rax)
nop
``````

So the compiler eliminates `-0.0`, but not `0.0`.
It would be nicer if it were the other way around, because `0.0` is faster to create (via `xor` instructions), making it the better choice if you can’t trust it to be eliminated. It’s also what you get from `zero(::Float64)` or `zero(::Type{Float64})`.

Although the people came up with the IEEE rules thought about this a lot more than I have, so there’s likely a good reason I don’t know about.

9 Likes

Thanks Steven!
In fact I was struggling with a version of `unique()` with certain tolerance like `1e-16`. I tracked the Julia implementation of `unique()` in `base.jl` line 1403:

`unique(A::AbstractArray; dims::Union{Colon,Integer} = :) = _unique_dims(A, dims)`

In the lines below that, I saw the actual code of `_unique_dims()`. It seems to me that the implementation of `_unique_dims()` is more “hash”.

My version of `unique()` with tolerance is
`@inline uniquen(v,n) = unique( v .|> (x->round(x,digits=n)) )`

I simply round the digits beyond certain precision.

Do you think this would be stable enough?

PS:
I just realized that
`unique([-0.0,+0.0])`
is a `Vector{Float64}` with TWO elements `[0.0,0.0]`

1 Like

are you from Mars? 1 Like

See the discussion here: https://github.com/JuliaLang/julia/issues/19147

A fundamental issue is that approximate equality is not an equivalence relation, so a `unique` function with a tolerance would have to give implementation-dependent results, or rather its exact behavior would be so complicated to describe that it would effectively be defined by its implementation.

1 Like

thank you, I’ve learnt this way of thinking now!

What you might do is something like:

``````unique(x -> round(x, sigdigits=13), myarray)
``````

which would perform the `unique` operation based on the first 13 significant digits. (You could use `digits` instead of `sigdigits` if you want an absolute tolerance).

But this might not necessarily be what you want, e.g. `unique(x -> round(x, sigdigits=13), [1.000000000000501, 1.0000000000005])` gives `[1.000000000000501, 1.0000000000005]` because those two values round differently in the 13th decimal digit. (Again, the fact that approximate equality is not an equivalence relation — it’s not transitive — means that any implementation of `unique` with a tolerance will have some odd behaviors.)

4 Likes

I see … I will investigate further with actual examples. I am dealing with degenerate eigenvalues, whose behaviour depends largely on the input matrix.

To get rid of the `unique([-0.0,+0.0])` trouble, I am currently using the following implementation:

`@inline correct_minus_00(x) = (real(x)==0.0 ? +0.0+imag(x)*im : x)`

`@inline uniquen(myarray,n) = unique(x -> correct_minus_00.(round(x, digits=n)), myarray)`

which adds a sign correction to your suggestion Hi Elrod, do you have more elegant solution of the `-0.0` problem? I mean, `-0.0` is practically equal to `0.0` in my code, so I use this
`@inline correct_minus_00(x) = (real(x)==0.0 ? +0.0+imag(x)*im : x)`
to correct the sign of `0.0`. But it looks silly …

If the sign of zero doesn’t matter in your code, why do you need to “correct” it?

I am dealing with degenerate eigenvalues, whose behaviour depends largely on the input matrix.

I would very strongly recommend that you re-think your algorithms to not rely on detecting exactly degenerate eigenvalues … sounds numerically unstable. What are you actually trying to compute?

1 Like

I am trying to identify degenerate eigenvectors within a certain tolerance. These degenerate eigenvectors means quantum states of same energy, they are physically very interesting. I just need to identify and collect them, but my current algorithm is not very stable and generates wrong collections. I will flatten all our discussions in my mind …

I have the following list of vectors

``````testv = [   Complex{Float64}[-1.0+1.11022e-15im, -1.0-1.04777e-15im],
Complex{Float64}[-1.0+1.22125e-15im, -1.0-1.35308e-15im],
Complex{Float64}[1.0+1.55431e-15im, 1.0-1.30451e-15im],
Complex{Float64}[0.0+1.0im, -0.0+1.0im, 0.0-1.0im, -0.0-1.0im],
Complex{Float64}[-1.0+1.11022e-15im, -1.0-1.04777e-15im, -1.0-4.79e-16im],
Complex{Float64}[-1.0+0.0im, -1.0-0.0im],
Complex{Float64}[-1.0+1.11022e-15im, -1.0-1.04777e-15im],
Complex{Float64}[-1.0+0.0im, -1.0-0.0im],
Complex{Float64}[-1.0+1.22125e-15im, -1.0-1.35308e-15im],
Complex{Float64}[-1.0+0.0im, -1.0-0.0im],
Complex{Float64}[-1.0+1.22125e-15im, -1.0-1.35308e-15im],
Complex{Float64}[-1.0+0.0im, -1.0-0.0im],   ]
``````

within a tolerance of `1e-14`, I hope to reduce them to distinct values. Because `unique([-0.0,0.0])` gives two elements, I couldn’t use it in an intuitive manner.

But I’ve got my solution now:

``````@inline uniquen(myarray,n) = unique(myarray .|> (x->(round(x,digits=n)+(0.0+0.0im))))

uniquen.(testv,12)
``````

yields desired results.

BTW,

``````f1(myarray) = unique(x->round(x,digits=12), myarray)
``````

and

``````f2(myarray) = unique(myarray .|> (x->round(x,digits=12)))
``````

behaves differently. For example:

``````julia> tv = Complex{Float64}[-1.0+1.11022e-15im, -1.0-1.04777e-15im, -1.0-4.79e-16im] ;

julia> f1(tv)
2-element Array{Complex{Float64},1}:
-1.0 + 1.11022e-15im
-1.0 - 1.04777e-15im

julia> f2(tv)
2-element Array{Complex{Float64},1}:
-1.0 + 0.0im
-1.0 - 0.0im
``````

I wished that `f1` could give me the result of `f2` (shown above), but it didn’t.

If you are looking for non-accidental degeneracies (for a Hamiltonian that commutes with some symmetry group), then a more reliable way is to project the eigenvector onto irreps of the symmetry group.

3 Likes

yes I am doing that in the afternoon Here’s my obligatory mention of one of my favorite-titled papers of all time, Branch cuts for complex elementary functions, or much ado about nothing’s sign bit, by the estimable William Kahan.

6 Likes

Thanks for the paper, I realized that handling zeros is a highly non-trivial task indeed.

1 Like